''' CASTEP Geometry Optimisation of Silicon '''
from datetime import datetime 
from chemsh       import *
from chemsh.utils import testutils
from chemsh.data import unitconvert

begin=datetime.now()

si2 = Fragment(coords='Si2_opt.cjson')
castep_obj = CASTEP(frag=si2,
functional='pbe',nbands=10,maxiter=100,
kpoints=[4,4,4],mixing_scheme='pulay',rand_seed=214708793,
 )

opt = Opt(theory=castep_obj, algorithm='lbfgs',maxcycle=200,
trust_radius='energy')
          
opt.run()

ecalc = opt.result.energy[0]
 
end=datetime.now()

print(f'Total elapsed time={(end-begin)}') 
print(f'Optimised Energy:{ecalc} Hartree')
print(f'Optimised Energy:{ecalc*unitconvert.Ha2eV} eV')

eref = -12.527061705
# NB 01/12/2023
# Due to variations in the generated pseudopotential (.usp file) for Si
# from machine to machine (despite using same compilers and libaries),
# a more relaxed tolerance has to be used for now.
testutils.validate(ecalc, eref, tol=1.0e-5)
