Harmonic calculator

Introduction

The local Harmonic Approximation of the potential energy surface (PES) is commonly applied in atomistic simulations to estimate entropy, i.e. free energy, at elevated temperatures (e.g. in ASE via thermochemistry). The term ‘harmonic’ refers to a second order Taylor series of the PES for a local reference configuration in Cartesian coordinates expressed in a Hessian matrix. With the Hessian matrix (e.g. computed numerically in ASE via vibrations) normal modes and harmonic vibrational frequencies can be obtained.

The following HarmonicCalculator can be used to compute energy and forces with a Hessian-based harmonic force field (HarmonicForceField). Moreover, it can be used to compute Anharmonic Corrections to the Harmonic Approximation. [1]

[1] (1,2,3,4,5,6)

Amsler, J. et al., Anharmonic Correction to Adsorption Free Energy from DFT-Based MD Using Thermodynamic Integration, J. Chem. Theory Comput. 2021, 17 (2), 1155-1169. doi: 10.1021/acs.jctc.0c01022.

class ase.calculators.harmonic.HarmonicCalculator(harmonicforcefield)

Class for calculations with a Hessian-based harmonic force field.

See HarmonicForceField and the literature. [1]

Parameters:

harmonicforcefield (HarmonicForceField) – Class for calculations with a Hessian-based harmonic force field.

class ase.calculators.harmonic.HarmonicForceField(ref_atoms, hessian_x, ref_energy=0.0, get_q_from_x=None, get_jacobian=None, cartesian=True, variable_orientation=False, hessian_limit=0.0, constrained_q=None, rcond=1e-07, zero_thresh=0.0)

Class that represents a Hessian-based harmonic force field.

Energy and forces of this force field are based on the Cartesian Hessian for a local reference configuration, i.e. if desired, on the Hessian matrix transformed to a user-defined coordinate system. The required Hessian has to be passed as an argument, e.g. predetermined numerically via central finite differences in Cartesian coordinates. Note that a potential being harmonic in Cartesian coordinates x is not necessarily equivalently harmonic in another coordinate system q, e.g. when the transformation between the coordinate systems is non-linear. By default, the force field is evaluated in Cartesian coordinates in which energy and forces are not rotationally and translationally invariant. Systems with variable orientation, require rotationally and translationally invariant calculations for which a set of appropriate coordinates has to be defined. This can be a set of (redundant) internal coordinates (bonds, angles, dihedrals, coordination numbers, …) or any other user-defined coordinate system.

Together with the HarmonicCalculator this HarmonicForceField can be used to compute Anharmonic Corrections to the Harmonic Approximation. [1]

Parameters:

• ref_atoms (Atoms object) – Reference structure for which energy (ref_energy) and Hessian matrix in Cartesian coordinates (hessian_x) are provided.

• hessian_x (numpy array) – Cartesian Hessian matrix for the reference structure ref_atoms. If a user-defined coordinate system is provided via get_q_from_x and get_jacobian, the Cartesian Hessian matrix is transformed to the user-defined coordinate system and back to Cartesian coordinates, thereby eliminating rotational and translational traits from the Hessian. The Hessian matrix obtained after this double-transformation is then used as the reference Hessian matrix to evaluate energy and forces for cartesian = True. For cartesian = False the reference Hessian matrix transformed to the user-defined coordinates is used to compute energy and forces.

• ref_energy (float) – Energy of the reference structure ref_atoms, typically in \(eV\).

• get_q_from_x (python function, default: None (Cartesian coordinates)) – Function that returns a vector of user-defined coordinates q for a given Atoms object ‘atoms’. The signature should be: get_q_from_x(atoms).

• get_jacobian (python function, default: None (Cartesian coordinates)) – Function that returns the geometric Jacobian matrix of the user-defined coordinates q w.r.t. Cartesian coordinates x defined as \(dq/dx\) (Wilson B-matrix) for a given Atoms object ‘atoms’. The signature should be: get_jacobian(atoms).

• cartesian (bool) – Set to True to evaluate energy and forces based on the reference Hessian (system harmonic in Cartesian coordinates). Set to False to evaluate energy and forces based on the reference Hessian transformed to user-defined coordinates (system harmonic in user-defined coordinates).

• hessian_limit (float) – Reconstruct the reference Hessian matrix with a lower limit for the eigenvalues, typically in \(eV/A^2\). Eigenvalues in the interval [zero_thresh, hessian_limit] are set to hessian_limit while the eigenvectors are left untouched.

• variable_orientation (bool) – Set to True if the investigated Atoms has got rotational degrees of freedom such that the orientation with respect to ref_atoms might be different (typically for molecules). Set to False to speed up the calculation when cartesian = True.

• constrained_q (list) – A list of indices ‘i’ of constrained coordinates \(q_i\) to be projected out from the Hessian matrix (e.g. remove forces along imaginary mode of a transition state).

• rcond (float) – Cutoff for singular value decomposition in the computation of the Moore-Penrose pseudo-inverse during transformation of the Hessian matrix. Equivalent to the rcond parameter in scipy.linalg.lstsq.

• zero_thresh (float) – Reconstruct the reference Hessian matrix with absolute eigenvalues below this threshold set to zero.

Note

The reference Hessians in x and q can be inspected via HarmonicForceField.hessian_x and HarmonicForceField.hessian_q.

Theory for Anharmonic Correction via Thermodynamic Integration (TI)

Thermodynamic integration (TI), i.e. \(\lambda\)-path integration, connects two thermodynamic states via a \(\lambda\)-path. Here, the TI begins from a reference system ‘0’ with known free energy (Harmonic Approximation) and the Anharmonic Correction is obtained via integration over the \(\lambda\)-path to the target system ‘1’ (the fully interacting anharmonic system). Hence, the free energy of the target system can be written as

\[A_1 = A_0 + \Delta A_{0 \rightarrow 1}\]

where the second term corresponds to the integral over the \(\lambda\)-path

\[\Delta A_{0 \rightarrow 1} = \int_0^1 d \lambda \langle H_1 - H_0 \rangle_\lambda\]

The term \(\langle ... \rangle_\lambda\) represents the NVT ensemble average of the system driven by the classical Hamiltonian \(\mathcal{H}_\lambda\) determined by the coupling parameter \(\lambda \in [0,1]\)

\[\mathcal{H}_\lambda = \lambda \mathcal{H}_1 + (1 - \lambda) \mathcal{H}_0\]

Since the Hamiltonians differ only in their potential energy contributions \(V_1\) and \(V_0\), the free energy change can be computed from the potentials

\[\Delta A_{0 \rightarrow 1} = \int_0^1 d \lambda \langle V_1 - V_0 \rangle_\lambda\]

The Cartesian coordinates x used in the common Harmonic Approximation are not insensitive to overall rotations and translations that must leave the total energy invariant. This limitation can be overcome by transformation of the Hessian in x to a suitable coordinate system q (e.g. internal coordinates). Since the force field of that Hessian which is harmonic in x is not necessarily equivalently harmonic in q, the free energy correction can be rewritten to

\[A_1 = A_{0,\mathbf{x}} + \Delta A_{0,\mathbf{x} \rightarrow 0,\mathbf{q}} + \Delta A_{0,\mathbf{q} \rightarrow 1}\]

The terms in this equation correspond to the free energy from the Harmonic Approximation with the reference Hessian (\(A_{0,\mathbf{x}}\)), the free energy change due to the coordinate transformation (\(\Delta A_{0,\mathbf{x} \rightarrow 0,\mathbf{q}}\)) obtained via TI (see Example 3) and the free energy change from the harmonic to the fully interacting system (\(\Delta A_{0,\mathbf{q} \rightarrow 1}\)) obtained via TI (see Example 4). Please see Amsler, J. et al. for details. [1]

Note

Anharmonicity is quantified by comparison of the total free energy \(A_1\) to the free energy contributions by the standard Harmonic Approximation with the unmodified Hessian. The reference Hessian and its free energy contribution \(A_{0,\mathbf{x}}\) have no meaning outside the TI procedure.

Examples

Prerequisites: Atoms object (ref_atoms), its energy (ref_energy) and Hessian (hessian_x).

Example 1: Cartesian coordinatates

In Cartesian coordinates, forces and energy are not invariant with respect to rotations and translations of the system.

import numpy as np
from ase.calculators.harmonic import HarmonicForceField, HarmonicCalculator
hff = HarmonicForceField(ref_atoms=ref_atoms, ref_energy=ref_energy,
                         hessian_x=hessian_x)
atoms = ref_atoms.copy()
atoms.calc = HarmonicCalculator(hff)

Note

Forces and energy can be computed via get_forces() and get_potential_energy() for any configuration that does not involve rotations with respect to the configuration of ref_atoms.

Warning

In case of system rotations, Cartesian coordinates return incorrect values and thus cannot be used without an appropriate coordinate system as demonstrated in the Supporting Information of Amsler, J. et al.. [1]

Example 2: Internal Coordinates

To compute forces and energy correctly even for rotated systems, a user-defined coordinate system must be provided. Within this coordinate system, energy and forces must be invariant with respect to rotations and translations of the system. For this purpose internal coordinates (distances, angles, dihedrals, coordination numbers and linear combinations thereof, etc.) are widely used. The following example works on a water molecule (H₂O) stored in ref_atoms.

dist_defs = [[0, 1], [1, 2], [2, 0]]  # define three distances by atom indices


def water_get_q_from_x(atoms):
    """Simple internal coordinates to describe water with three distances."""
    q_vec = [atoms.get_distance(i, j) for i, j in dist_defs]
    return np.asarray(q_vec)


def water_get_jacobian(atoms):
    """Function to return the Jacobian for the water molecule described by
    three distances."""
    pos = atoms.get_positions()
    dist_vecs = [pos[j] - pos[i] for i, j in dist_defs]
    derivs = get_distances_derivatives(dist_vecs)
    jac = []
    for i, defin in enumerate(dist_defs):
        dqi_dxj = np.zeros(ref_pos.shape)
        for j, deriv in enumerate(derivs[i]):
            dqi_dxj[defin[j]] = deriv
        jac.append(dqi_dxj.flatten())
    return np.asarray(jac)

parameters = {'ref_atoms': ref_atoms, 'ref_energy': ref_energy,
              'hessian_x': hessian_x, 'get_q_from_x': water_get_q_from_x,
              'get_jacobian': water_get_jacobian, 'cartesian': False}
hff = HarmonicForceField(**parameters)  # calculation in internals
calc = HarmonicCalculator(hff)

Example 3: Free Energy Change due to Coordinate Transformation

A transformation of the coordinate system may transform the force field. The change in free energy due to this transformation (\(\Delta A_{0,\mathbf{x} \rightarrow 0,\mathbf{q}}\)) can be computed via thermodynamic (\(\lambda\)-path) integration. [1]

parameters = {'ref_atoms': ref_atoms, 'ref_energy': ref_energy,
              'hessian_x': hessian_x, 'get_q_from_x': water_get_q_from_x,
              'get_jacobian': water_get_jacobian, 'cartesian': True,
              'variable_orientation': True}
hff_1 = HarmonicForceField(**parameters)
calc_harmonic_1 = HarmonicCalculator(hff_1)
parameters['cartesian'] = False
hff_0 = HarmonicForceField(**parameters)
calc_harmonic_0 = HarmonicCalculator(hff_0)
ediffs = {}  # collect energy difference for varying lambda coupling
lambs = [0.00, 0.25, 0.50, 0.75, 1.00]  # integration grid
for lamb in lambs:
    ediffs[lamb] = []
    calc_linearCombi = MixedCalculator(calc_harmonic_0, calc_harmonic_1,
                                       1 - lamb, lamb)
    atoms = ref_atoms.copy()
    atoms.calc = calc_linearCombi
    MaxwellBoltzmannDistribution(atoms, temperature_K=300, force_temp=True)
    Stationary(atoms)
    ZeroRotation(atoms)
    with Andersen(atoms, 0.5 * fs, temperature_K=300, andersen_prob=0.05,
                  fixcm=False) as dyn:
        for _ in dyn.irun(50):  # should be much longer for production runs
            e0, e1 = calc_linearCombi.get_energy_contributions(atoms)
            ediffs[lamb].append(float(e1) - float(e0))
        ediffs[lamb] = np.mean(ediffs[lamb])
dA = trapezoid([ediffs[lamb] for lamb in lambs], x=lambs)  # anharm. corr.

Integration of the mean energy differences (‘ediffs’) over the integration grid (\(\lambda\) path) leads to the change in free energy due to the coordinate transformation.

Example 4: Anharmonic Corrections

The strategy in Example 3 can be used to compute anharmonic corrections to the Harmonic Approximation when the HarmonicCalculator is coupled with a calculator that can compute interactions beyond the Harmonic Approximation, e.g. vasp.

Note

The obtained Anharmonic Correction applies to the Harmonic Approximation (\(A_{0,\mathbf{x}}\)) of the reference system with the reference Hessian which is generated during initialization of the Calculator and may differ from the standard Harmonic Approximation. The vibrations for the reference system can be computed numerically with high accuracy.

>>> from ase.vibrations import Vibrations
>>> atoms = ref_atoms.copy()
>>> atoms.calc = calc_harmonic_0  # with cartesian=True
>>> vib = Vibrations(atoms, nfree=4, delta=1e-5)
>>> vib.run()