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)