Numerical gradients

This module implements a simple numerical differentiation scheme with respect to molecular coordinates. We define a new job type NumGradJob 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 accuracy of the numerical differentiation scheme (npoints) and size and unit of displacement step.

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 NumGradJob are directly passed to children jobs, so creating a NumGradJob strongly resembles creating a regular single point job.

The dedicated Results subclass for NumGradJob takes care of extracting the gradient from the results of all single point jobs. Any function that takes results of a single point job and returns a single number can be used with the get_gradient method defined in NumGradResults

The source code of the whole module with both abovementioned classes:

from collections import OrderedDict
from itertools import product

from scm.plams.core.basejob import MultiJob
from scm.plams.core.results import Results

__all__ = ["NumGradJob", "NumGradResults"]  # names exported to the main namespace


class NumGradResults(Results):

    def _get_gradient_component(self, atom, coord, extract):
        """Get gradient component for *atom* along *coord*.

        *atom* should be integer, *coord* a single character string with 'x', 'y' or 'z'.
        *extract* should be a function that takes results of a single point job and
        returns a single number.
        """
        energies = [extract(self.job.children[(atom, coord, i)].results) for i in self.job.steps]
        coeffs, denom = self.job._coeffs[self.job.npoints]
        return sum([c * e for c, e in zip(coeffs, energies)]) / (denom * self.job.step)

    def get_gradient(self, extract):
        """Get the full gradient vector. Returns a list of length 3N,
        where N is the number of atoms in the calculated molecule.

        *extract* should be a function that takes results of a single point job and
        returns a single number.
        """
        return [
            self._get_gradient_component(atom, coord, extract)
            for atom, coord in product(range(1, 1 + len(self.job.molecule)), "xyz")
        ]


class NumGradJob(MultiJob):
    """A class for calculating numerical gradients of energy
    (or some other single-valued property) with respect to
    cartesian displacements.

    The differentiation is done using the final difference method
    with 2, 4, 6 or 8 points. The length of the step can be adjusted.
    """

    _result_type = NumGradResults

    # finite difference coefficients for different number of points:
    _coeffs = {
        2: ([-1, 1], 2),
        4: ([1, -8, 8, -1], 12),
        6: ([-1, 9, -75, 75, -9, 1], 60),
        8: ([3, -32, 168, -672, 672, -168, 32, -3], 840),
    }

    def __init__(self, molecule, npoints=2, step=0.01, unit="angstrom", jobtype=None, **kwargs):
        # initialize parent class and store all additional constructor arguments
        MultiJob.__init__(self, children=OrderedDict(), **kwargs)
        self.molecule = molecule
        self.npoints = npoints
        self.step = step
        self.unit = unit
        self.jobtype = jobtype  # who is going to calculate single points

    def prerun(self):  # noqa F811
        self.steps = list(range(-(self.npoints // 2), 0)) + list(range(1, self.npoints // 2 + 1))

        for atom, axis, i in product(range(1, 1 + len(self.molecule)), "xyz", self.steps):
            v = (
                self.step * i if axis == "x" else 0,
                self.step * i if axis == "y" else 0,
                self.step * i if axis == "z" else 0,
            )
            newmol = self.molecule.copy()
            newmol[atom].translate(v, self.unit)
            newname = "{}_{}_{}".format(atom, axis, i)
            # settings of NumGradJob are directly passed to children single point jobs
            self.children[(atom, axis, i)] = self.jobtype(name=newname, molecule=newmol, settings=self.settings)

An example usage:

#!/usr/bin/env plams
import numpy as np
from scm.plams.recipes.numgrad import NumGradJob

config.default_jobrunner = JobRunner(parallel=True, maxjobs=8)

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[""]


s = Settings()
s.input.AMS.Task = "SinglePoint"
s.input.ForceField.Type = "UFF"
s.runscript.nproc = 1

j = NumGradJob(name="numgrad", molecule=mol, settings=s, jobtype=AMSJob)
r = j.run()
r.wait()

s_analytical = s.copy()
s_analytical.input.AMS.Properties.Gradients = "Yes"
analytical_job = AMSJob(name="analytical", molecule=mol, settings=s_analytical)
analytical_job.run()

numerical_gradients = np.array(r.get_gradient(AMSResults.get_energy)).reshape(-1, 3)
analytical_gradients = np.array(analytical_job.results.get_gradients()).reshape(
    -1, 3
) * Units.convert(1.0, "bohr^-1", "angstrom^-1")

print("Numerical Gradients (NumGradJob), hartree/angstrom:")
print(numerical_gradients)
print("Analytical Gradients, hartree/angstrom:")
print(analytical_gradients)