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

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

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))