# Gaussian DFT (BLYP) single point energy

from chemsh       import *
from chemsh.utils import testutils
from chemsh.data  import unitconvert

water = Fragment(coords="water.cjson")

# default functional is 'blyp'
# root state 2: the 2nd excited state
nwchem = NWChem(frag=water, method='tddft', basis="3-21g", nstates=4, target_state=2, rpa=True)

sp = SP(theory=nwchem, gradients=True)

sp.run()

ecalc_abs = sp.result.energy
# NB: Python starts with 0! (They are absolute energies in result.excitations)
ecalc_exc = (sp.result.excitations['energies'][1]-sp.result.E0)*unitconvert.Ha2eV
fcalc_exc = sp.result.excitations['f'][1]

# Test result
# Py-ChemShell/Gaussian16 as of 11/01/2021

# absolute energy (in a.u.) of the S2 (A1) state
eref_abs = -75.591158734700

# excitation energy (in eV) of the S2 (A1) state
eref_exc = 9.6799

# oscillator strength of the S2 (A1) state
fref_exc = 0.07990


testutils.validate(ecalc_abs, eref_abs)
testutils.validate(ecalc_exc, eref_exc, tol=0.0001)
testutils.validate(fcalc_exc, fref_exc, tol=0.0001)
