Numerical Hessian¶
This module implements a simple scheme for calculating a numerical Hessian matrix.
We define a new job type NumHessJob
by extending MultiJob
.
The constructor (__init__
) of this new job accepts several new arguments and simply stores them.
These new arguments define the initial Molecule
, the type of job used for single point calculations (jobtype
), the size and unit of displacement step and the way of extracting gradients from single point results.
Then the prerun()
method takes the given Molecule
and creates multiple copies of it, each one with one atom displaced along one axis.
For each of these molecules an instance of single point job is created and stored in the children
dictionary.
Settings of NumHessJob
are directly passed to children jobs, so creating a NumHessJob
strongly resembles creating a regular single point job.
The dedicated Results
subclass for NumHessJob
takes care of extracting the Hessian from results of all single point jobs.
The returned Hessian can optionally be mass weighted.
The source code of the whole module with both aforementioned classes:
from collections import OrderedDict
from itertools import product
from scm.plams.core.basejob import MultiJob
from scm.plams.core.results import Results
import numpy as np
from scm.plams.tools.units import Units
__all__ = ["NumHessJob", "NumHessResults"] # names exported to the main namespace
class NumHessResults(Results):
def get_hessian(self, mass_weighted=False):
j = self.job
n = len(j.molecule)
hessian = np.empty((3 * n, 3 * n))
for atom, axis in product(range(1, n + 1), range(3)):
v1 = np.array(j.gradient(j.children[(atom, axis, -1)].results))
v2 = np.array(j.gradient(j.children[(atom, axis, 1)].results))
hessian[3 * (atom - 1) + axis] = (v2 - v1) / (2 * Units.convert(j.step, "angstrom", "bohr"))
if mass_weighted:
masses = [at.mass for at in j.molecule]
weights = np.empty((n, n))
for i, j in product(range(n), repeat=2):
weights[i, j] = masses[i] * masses[j]
hessian /= np.kron(weights, np.ones((3, 3)))
return (hessian + hessian.T) / 2
class NumHessJob(MultiJob):
"""A class for calculating numerical hessian.
The length and unit of the geometry displacement step can be adjusted.
*gradient* should be a function that takes results of *jobtype* and
returns energy gradient in hartree/bohr.
"""
_result_type = NumHessResults
def __init__(self, molecule, step=0.01, unit="angstrom", jobtype=None, gradient=None, **kwargs):
MultiJob.__init__(self, children=OrderedDict(), **kwargs)
self.molecule = molecule
self.step = Units.convert(step, unit, "angstrom")
self.jobtype = jobtype # who is going to calculate single points
self.gradient = gradient # function extracting gradients from children
def prerun(self): # noqa F811
for atom, axis, step in product(range(1, 1 + len(self.molecule)), range(3), [-1, 1]):
vec = [0, 0, 0]
vec[axis] = self.step * step
newmol = self.molecule.copy()
newmol[atom].translate(vec)
newname = "{}_{}_{}".format(atom, axis, step)
self.children[(atom, axis, step)] = self.jobtype(name=newname, molecule=newmol, settings=self.settings)
To follow along, either
Download
NumHess.py
(run as$AMSBIN/amspython NumHess.py
).Download
NumHess.ipynb
(see also: how to install Jupyterlab in AMS)
Worked Example¶
Initialization¶
from scm.plams import *
# this line is not required in AMS2025+
init()
PLAMS working folder: /path/plams/examples/NumHess/plams_workdir
Define molecule¶
Normally you would read the molecule from an xyz file
# mol = Molecule('methanol.xyz')
mol = AMSJob.from_input(
"""
System
Atoms
C 0.000000000000 0.138569980000 0.355570700000
O 0.000000000000 0.187935770000 -1.074466460000
H 0.882876920000 -0.383123830000 0.697839450000
H -0.882876940000 -0.383123830000 0.697839450000
H 0.000000000000 1.145042790000 0.750208830000
End
End
"""
).molecule[""]
Setup and run job¶
Here we use the fast DFTB engine, for ADF it is recommended to disable symmetry
s = Settings()
s.input.ams.task = "SinglePoint"
s.input.ams.Properties.Gradients = "Yes"
# s.input.adf.basis.type = 'DZP'
# s.input.adf.symmetry = 'NOSYM'
# s.input.adf.xc.gga = 'PW91'
s.input.DFTB.Model = "GFN1-xTB"
s.runscript.nproc = 1
j = NumHessJob(name="test", molecule=mol, settings=s, jobtype=AMSJob, gradient=lambda x: x.get_gradients().reshape(-1))
r = j.run(jobrunner=JobRunner(parallel=True, maxjobs=8))
[06.03|13:55:22] JOB test STARTED
[06.03|13:55:22] JOB test RUNNING
[06.03|13:55:22] JOB test/1_0_-1 STARTED
[06.03|13:55:22] JOB test/1_0_1 STARTED
Print results¶
print(r.get_hessian(mass_weighted=True))
[06.03|13:55:22] JOB test/1_1_-1 STARTED
[06.03|13:55:22] Waiting for job test to finish
[06.03|13:55:22] JOB test/1_1_1 STARTED
[06.03|13:55:22] JOB test/1_2_-1 STARTED
[06.03|13:55:22] JOB test/1_2_1 STARTED
[06.03|13:55:22] JOB test/1_0_-1 RUNNING
[06.03|13:55:22] JOB test/2_0_-1 STARTED
[06.03|13:55:22] JOB test/2_0_1 STARTED
[06.03|13:55:22] JOB test/2_1_-1 STARTED
[06.03|13:55:22] JOB test/2_1_1 STARTED
[06.03|13:55:22] JOB test/2_2_-1 STARTED
... (PLAMS log lines truncated) ...
[[ 2.43775999e-03 -1.24945968e-10 8.35714013e-11 -2.16941328e-05
5.02435569e-12 -7.14677254e-12 -1.60687481e-02 1.15383294e-02
-7.21385604e-03 -1.60687464e-02 -1.15383280e-02 7.21385522e-03
3.47625119e-03 3.62021265e-11 -5.46656826e-11]
[-1.24945968e-10 2.43617119e-03 1.02498360e-05 4.74486821e-12
-2.19212819e-05 2.94502688e-05 1.14532498e-02 -3.28822540e-03
4.41476048e-03 -1.14532483e-02 -3.28822472e-03 4.41475997e-03
-8.03126599e-11 -2.20627127e-02 -9.41894828e-03]
[ 8.35714013e-11 1.02498360e-05 2.14233490e-03 7.19072217e-13
2.90465272e-05 -8.63468227e-04 -4.79807048e-03 3.01867574e-03
-3.68101848e-03 4.79806931e-03 3.01867515e-03 -3.68101833e-03
-4.67601951e-11 -6.61690618e-03 -4.46017458e-03]
[-2.16941328e-05 4.74486821e-12 7.19072217e-13 1.49619173e-04
-2.14387656e-12 8.46463904e-13 -3.85034451e-04 -5.25094480e-04
-1.04652467e-03 -3.85034464e-04 5.25094462e-04 1.04652463e-03
-1.35573946e-03 -4.86485419e-12 2.14361038e-11]
[ 5.02435569e-12 -2.19212819e-05 2.90465272e-05 -2.14387656e-12
1.50336369e-04 -3.08725780e-05 -5.36109667e-04 -1.06913031e-03
6.27011593e-04 5.36109618e-04 -1.06913034e-03 6.27011524e-04
1.40055968e-11 3.87414930e-06 -1.11023668e-03]
[-7.14677254e-12 2.94502688e-05 -8.63468227e-04 8.46463904e-13
-3.08725780e-05 1.04165349e-03 -7.28706894e-04 4.42296851e-04
-2.04725210e-03 7.28706950e-04 4.42296845e-04 -2.04725207e-03
1.96785022e-11 -7.46717101e-04 -2.14861643e-03]
[-1.60687481e-02 1.14532498e-02 -4.79807048e-03 -3.85034451e-04
-5.36109667e-04 -7.28706894e-04 2.30394724e-01 -1.34285111e-01
7.67026070e-02 -2.16065892e-02 5.58885425e-03 -1.11180652e-02
-1.12833086e-02 6.64824977e-04 3.15373514e-03]
[ 1.15383294e-02 -3.28822540e-03 3.01867574e-03 -5.25094480e-04
-1.06913031e-03 4.42296851e-04 -1.34285111e-01 8.15853116e-02
-4.64977897e-02 -5.58885505e-03 -7.11299026e-03 -1.05553100e-02
1.06516082e-02 -1.84750271e-02 1.40635538e-02]
[-7.21385604e-03 4.41476048e-03 -3.68101848e-03 -1.04652467e-03
6.27011593e-04 -2.04725210e-03 7.67026070e-02 -4.64977897e-02
6.64040738e-02 1.11180655e-02 -1.05553100e-02 4.44753112e-03
1.47911831e-02 -5.53928405e-03 5.49673857e-03]
[-1.60687464e-02 -1.14532483e-02 4.79806931e-03 -3.85034464e-04
5.36109618e-04 7.28706950e-04 -2.16065892e-02 -5.58885505e-03
1.11180655e-02 2.30394703e-01 1.34285094e-01 -7.67025954e-02
-1.12833089e-02 -6.64823699e-04 -3.15373388e-03]
[-1.15383280e-02 -3.28822472e-03 3.01867515e-03 5.25094462e-04
-1.06913034e-03 4.42296845e-04 5.58885425e-03 -7.11299026e-03
-1.05553100e-02 1.34285094e-01 8.15853035e-02 -4.64977824e-02
-1.06516072e-02 -1.84750267e-02 1.40635536e-02]
[ 7.21385522e-03 4.41475997e-03 -3.68101833e-03 1.04652463e-03
6.27011524e-04 -2.04725207e-03 -1.11180652e-02 -1.05553100e-02
4.44753112e-03 -7.67025954e-02 -4.64977824e-02 6.64040713e-02
-1.47911830e-02 -5.53928397e-03 5.49673857e-03]
[ 3.47625119e-03 -8.03126599e-11 -4.67601951e-11 -1.35573946e-03
1.40055968e-11 1.96785022e-11 -1.12833086e-02 1.06516082e-02
1.47911831e-02 -1.12833089e-02 -1.06516072e-02 -1.47911830e-02
2.46727294e-03 -2.63269039e-10 2.12920387e-11]
[ 3.62021265e-11 -2.20627127e-02 -6.61690618e-03 -4.86485419e-12
3.87414930e-06 -7.46717101e-04 6.64824977e-04 -1.84750271e-02
-5.53928405e-03 -6.64823699e-04 -1.84750267e-02 -5.53928397e-03
-2.63269039e-10 2.99742684e-01 1.01775207e-01]
[-5.46656826e-11 -9.41894828e-03 -4.46017458e-03 2.14361038e-11
-1.11023668e-03 -2.14861643e-03 3.15373514e-03 1.40635538e-02
5.49673857e-03 -3.15373388e-03 1.40635536e-02 5.49673857e-03
2.12920387e-11 1.01775207e-01 7.62490469e-02]]
Complete Python code¶
#!/usr/bin/env amspython
# coding: utf-8
# ## Initialization
from scm.plams import *
# this line is not required in AMS2025+
init()
# ## Define molecule
# Normally you would read the molecule from an xyz file
# mol = Molecule('methanol.xyz')
mol = AMSJob.from_input(
"""
System
Atoms
C 0.000000000000 0.138569980000 0.355570700000
O 0.000000000000 0.187935770000 -1.074466460000
H 0.882876920000 -0.383123830000 0.697839450000
H -0.882876940000 -0.383123830000 0.697839450000
H 0.000000000000 1.145042790000 0.750208830000
End
End
"""
).molecule[""]
# ## Setup and run job
# Here we use the fast DFTB engine, for ADF it is recommended to disable symmetry
s = Settings()
s.input.ams.task = "SinglePoint"
s.input.ams.Properties.Gradients = "Yes"
# s.input.adf.basis.type = 'DZP'
# s.input.adf.symmetry = 'NOSYM'
# s.input.adf.xc.gga = 'PW91'
s.input.DFTB.Model = "GFN1-xTB"
s.runscript.nproc = 1
j = NumHessJob(name="test", molecule=mol, settings=s, jobtype=AMSJob, gradient=lambda x: x.get_gradients().reshape(-1))
r = j.run(jobrunner=JobRunner(parallel=True, maxjobs=8))
# ## Print results
print(r.get_hessian(mass_weighted=True))