Example: Excited state gradients: plams¶
Download SP_LR-TDDFTB_gradients.run
#!/bin/sh
cp $AMSHOME/examples/dftb/SP_LR-TDDFTB_gradients/SP_LR-TDDFTB_gradients.plms .
cp -r $AMSHOME/examples/dftb/SP_LR-TDDFTB_gradients/molecules .
cp -r $AMSHOME/examples/dftb/SP_LR-TDDFTB_gradients/numgrad_precalc .
export NSCM=1
$AMSBIN/plams SP_LR-TDDFTB_gradients.plms
Download SP_LR-TDDFTB_gradients.plms
import numpy as np
import os.path
# our test molecules and their excitations for which we want to calculate the gradients
tests = [
('acetamide' , 'singlet', 2, False),
('acetone' , 'singlet', 5, False),
('adenine' , 'singlet', 2, False),
('benzene' , 'singlet', 1, False),
('benzoquinone' , 'singlet', 1, False),
('butadiene' , 'singlet', 1, False),
('carbonmonoxide' , 'singlet', 8, False),
('cyclopentadiene', 'singlet', 1, False),
('cyclopropene' , 'singlet', 1, False),
('cytosine' , 'singlet', 2, False),
('ethene' , 'singlet', 2, False),
('formaldehyde' , 'singlet', 4, False),
('formamide' , 'singlet', 2, False),
('furan' , 'singlet', 1, False),
('hexatriene' , 'singlet', 1, False),
('imidazole' , 'singlet', 2, False),
('naphthalene' , 'singlet', 1, False),
('nitrogen' , 'singlet', 8, False),
('norbornadiene' , 'singlet', 1, False),
('octatetraene' , 'singlet', 1, False),
('propanamide' , 'singlet', 2, False),
('pyrazine' , 'singlet', 3, False),
('pyridazine' , 'singlet', 4, False),
('pyridine' , 'singlet', 3, False),
('pyrimidine' , 'singlet', 6, False),
('pyrrole' , 'singlet', 3, False),
('tetrazine' , 'singlet', 6, False),
('thymine' , 'singlet', 3, False),
('triazine' , 'singlet', 6, False),
('uracil' , 'singlet', 3, False),
('betacarotene' , 'singlet', 1, False),
('acetamide' , 'triplet', 3, False),
('acetone' , 'triplet', 1, False),
('adenine' , 'triplet', 4, False),
('benzene' , 'triplet', 1, False),
('benzoquinone' , 'triplet', 4, False),
('butadiene' , 'triplet', 1, False),
('carbonmonoxide' , 'triplet', 3, False),
('cyclopentadiene', 'triplet', 2, False),
('cyclopropene' , 'triplet', 1, False),
('cytosine' , 'triplet', 3, False),
('ethene' , 'triplet', 1, False),
('formaldehyde' , 'triplet', 7, False),
('formamide' , 'triplet', 10, False),
('furan' , 'triplet', 2, False),
('hexatriene' , 'triplet', 2, False),
('imidazole' , 'triplet', 1, False),
('naphthalene' , 'triplet', 9, False),
('nitrogen' , 'triplet', 1, False),
('norbornadiene' , 'triplet', 2, False),
('octatetraene' , 'triplet', 1, False),
('propanamide' , 'triplet', 1, False),
('pyrazine' , 'triplet', 4, False),
('pyridazine' , 'triplet', 1, False),
('pyridine' , 'triplet', 5, False),
('pyrimidine' , 'triplet', 1, False),
('pyrrole' , 'triplet', 2, False),
('tetrazine' , 'triplet', 6, False),
('thymine' , 'triplet', 4, False),
('triazine' , 'triplet', 5, False),
('uracil' , 'triplet', 2, False),
('betacarotene' , 'triplet', 1, False),
('acetamide' , 'singlet', 2, True),
('acetone' , 'singlet', 5, True),
('adenine' , 'singlet', 2, True),
('benzene' , 'singlet', 1, True),
('benzoquinone' , 'singlet', 1, True),
('butadiene' , 'singlet', 1, True),
('carbonmonoxide' , 'singlet', 8, True),
('cyclopentadiene', 'singlet', 1, True),
('cyclopropene' , 'singlet', 1, True),
('cytosine' , 'singlet', 2, True),
('ethene' , 'singlet', 2, True),
('formaldehyde' , 'singlet', 4, True),
('formamide' , 'singlet', 2, True),
('furan' , 'singlet', 1, True),
('hexatriene' , 'singlet', 1, True),
('imidazole' , 'singlet', 2, True),
('naphthalene' , 'singlet', 1, True),
('nitrogen' , 'singlet', 8, True),
('norbornadiene' , 'singlet', 1, True),
('octatetraene' , 'singlet', 1, True),
('propanamide' , 'singlet', 2, True),
('pyrazine' , 'singlet', 3, True),
('pyridazine' , 'singlet', 4, True),
('pyridine' , 'singlet', 3, True),
('pyrimidine' , 'singlet', 6, True),
('pyrrole' , 'singlet', 3, True),
('tetrazine' , 'singlet', 6, True),
('thymine' , 'singlet', 3, True),
('triazine' , 'singlet', 6, True),
('uracil' , 'singlet', 3, True),
('betacarotene' , 'singlet', 1, True),
('acetamide' , 'triplet', 3, True),
('acetone' , 'triplet', 1, True),
('adenine' , 'triplet', 4, True),
('benzene' , 'triplet', 1, True),
('benzoquinone' , 'triplet', 4, True),
('butadiene' , 'triplet', 1, True),
('carbonmonoxide' , 'triplet', 3, True),
('cyclopentadiene', 'triplet', 2, True),
('cyclopropene' , 'triplet', 1, True),
('cytosine' , 'triplet', 3, True),
('ethene' , 'triplet', 1, True),
('formaldehyde' , 'triplet', 7, True),
('formamide' , 'triplet', 10, True),
('furan' , 'triplet', 2, True),
('hexatriene' , 'triplet', 2, True),
('imidazole' , 'triplet', 1, True),
('naphthalene' , 'triplet', 9, True),
('nitrogen' , 'triplet', 1, True),
('norbornadiene' , 'triplet', 2, True),
('octatetraene' , 'triplet', 1, True),
('propanamide' , 'triplet', 1, True),
('pyrazine' , 'triplet', 4, True),
('pyridazine' , 'triplet', 1, True),
('pyridine' , 'triplet', 5, True),
('pyrimidine' , 'triplet', 1, True),
('pyrrole' , 'triplet', 2, True),
('tetrazine' , 'triplet', 6, True),
('thymine' , 'triplet', 4, True),
('triazine' , 'triplet', 5, True),
('uracil' , 'triplet', 2, True),
('betacarotene' , 'triplet', 1, True),
]
# numpy setup
#np.set_printoptions(precision=8, suppress=True)
np.set_printoptions(formatter={'float': ' {: 0.8f} '.format})
# plams set up
config.log.stdout = 0
# common input for all tests
comin = Settings()
comin.input.ams.task = 'SinglePoint'
comin.input.ams.properties.gradients = True
comin.input.DFTB.model = 'SCC-DFTB'
comin.input.DFTB.resourcesdir = 'DFTB.org/mio-1-1'
comin.input.DFTB.properties.excitations.tddftb['print'] = 'evcontribs'
failedtests = []
for test in tests:
molname = test[0]
multi = test[1]
excit = test[2]
ldep = test[3]
if multi == 'singlet':
kfsec = 'SS'
else:
kfsec = 'ST'
if ldep:
ldpf = 'ldep'
else:
ldpf = 'noldep'
print('\nTESTING: '+molname+' '+multi+' '+str(excit)+' '+ldpf)
teststr = molname+'_'+multi+'_'+str(excit)+'_'+ldpf
comin_thistest = comin.copy()
comin_thistest.input.DFTB.properties.excitations.tddftb.calc = multi
comin_thistest.input.DFTB.properties.excitations.tddftb.lowest = excit
if ldep: comin_thistest.input.DFTB.scc.orbitaldependent = 'yes'
mol = Molecule(filename='./molecules/'+molname+'.xyz')
# numerical gradient
if os.path.isfile('./numgrad_precalc/'+teststr+'.npy'):
print('Precalculated numerical gradient found -> reading from file')
numgrad = np.load('./numgrad_precalc/'+teststr+'.npy')
else:
print('Precalculated numerical gradient not found -> calculating')
numgradjob = AMSNumGradJob(name=teststr+'_numgrad', molecule=mol, npoints=3, step=0.001)
numgradjob.settings.child = comin_thistest
numgradresults = numgradjob.run()
def exenergy(results):
if excit == 1:
return results.readrkf('Excitations '+kfsec+' A', 'excenergies', file='dftb')
else:
return results.readrkf('Excitations '+kfsec+' A', 'excenergies', file='dftb')[excit - 1]
numgrad = np.empty([len(mol), 3])
for n in range(1,len(mol)+1):
numgrad[n-1,0] = numgradresults.get_gradient(n,'x', func=exenergy)
numgrad[n-1,1] = numgradresults.get_gradient(n,'y', func=exenergy)
numgrad[n-1,2] = numgradresults.get_gradient(n,'z', func=exenergy)
numgrad = Units.conversion_ratio('bohr','angstrom') * numgrad
# write numerical gradient to file
print('Saving numerical gradient to file')
np.save('./numgrad_precalc/'+teststr+'.npy', numgrad)
print('Numerical gradient =')
print(numgrad)
# analytical gradient
job = AMSJob(name=teststr+'_anagrad', molecule=mol, settings=comin_thistest)
job.settings.input.DFTB.properties.excitations.tddftbgradients.excitation = excit
results = job.run()
anagrad = np.array(results.readrkf('Excitations '+kfsec+' A', 'gradient '+str(excit), file='dftb')).reshape((-1,3))
print('Analytical gradient =')
print(anagrad)
# print the difference between analytical and numerical gradients
diff = numgrad - anagrad
print('Deviation =')
print(diff)
# check if the difference is small enough
passed = np.allclose(numgrad, anagrad, atol=1e-4)
if passed:
print('TEST FINISHED: PASSED!')
else:
print('TEST FINISHED: FAILED!')
failedtests.append(test)
print('\nTESTS PASSED: '+str(len(tests)-len(failedtests))+'/'+str(len(tests)))
for test in failedtests:
molname = test[0]
multi = test[1]
excit = test[2]
ldep = test[3]
if ldep:
ldpf = 'ldep'
else:
ldpf = 'noldep'
print('FAILED: '+molname+' '+multi+' '+str(excit)+' '+ldpf)