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]
- class ase.calculators.harmonic.HarmonicCalculator(harmonicforcefield)[source]
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)[source]
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
thisHarmonicForceField
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 viaget_q_from_x
andget_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 forcartesian = True
. Forcartesian = 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 tohessian_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 toref_atoms
might be different (typically for molecules). Set to False to speed up the calculation whencartesian = 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
where the second term corresponds to the integral over the \(\lambda\)-path
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]\)
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
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
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 (H2O) 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()