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 abovementioned classes:
from collections import OrderedDict
from itertools import product
import numpy as np
from ..core.results import Results
from ..core.basejob import MultiJob
from ..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):
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)
An example usage:
mol = Molecule('methanol.xyz')
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.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(r.get_hessian(mass_weighted=True))