# Gaussian DFT (BLYP) single point energy

from chemsh       import *
from chemsh.utils import testutils

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

# default functional is 'blyp'
nwchem = NWChem(frag=water, method='tddft', basis='3-21g', nstates=4, rpa=True, mult=3)

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

sp.run()

#ecalc_abs = sp.result.energy
## NB: Python starts with 0!
#ecalc_exc = sp.result.excitations['energies'][1]
#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)
#testutils.validate(fcalc_exc, fref_exc)
#

# ground state
eref_0 = -75.668746375925
testutils.validate(sp.result.energy, eref_0)
# excited state 1
eref_1 = -75.575852258471
testutils.validate(sp.result.excitations['energies'][0], eref_1)
# excited state 2
eref_2 = -75.558669752449
testutils.validate(sp.result.excitations['energies'][1], eref_2)
# excited state 3
eref_3 = -75.438350640740
testutils.validate(sp.result.excitations['energies'][2], eref_3)
# excited state 4
eref_4 = -75.212445289312
testutils.validate(sp.result.excitations['energies'][3], eref_4)

