1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
|
.. module:: ase.calculators.harmonic
.. _harmonic:
===================
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 :mod:`~ase.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
:mod:`~ase.vibrations`) normal modes and harmonic vibrational frequencies can
be obtained.
The following :class:`HarmonicCalculator` can be used to compute energy and forces
with a Hessian-based harmonic force field (:class:`HarmonicForceField`).
Moreover, it can be used to compute Anharmonic Corrections to the
Harmonic Approximation. [1]_
.. [1] 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`.
.. autoclass:: ase.calculators.harmonic.HarmonicCalculator
.. autoclass:: ase.calculators.harmonic.HarmonicForceField
.. 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
.. math::
A_1 = A_0 + \Delta A_{0 \rightarrow 1}
where the second term corresponds to the integral over the `\lambda`-path
.. math::
\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]`
.. math::
\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
.. math::
\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
.. math::
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: :class:`~ase.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.
.. code-block:: python
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 :meth:`~ase.Atoms.get_forces` and
:meth:`~ase.Atoms.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 (:mol:`H_2O`) stored in
``ref_atoms``.
.. literalinclude:: ../../../ase/test/calculator/test_harmonic.py
:language: python
:start-after: start doc example 3
:end-before: end doc example 3
.. literalinclude:: ../../../ase/test/calculator/test_harmonic.py
:language: python
:start-after: test_internals():
:end-before: atoms = setup_water(calc) # distorted copy of ref_atoms
:dedent: 4
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]_
.. literalinclude:: ../../../ase/test/calculator/test_harmonic.py
:language: python
:start-after: test_thermodynamic_integration():
:end-before: assert -0.005 < dA < 0.005
:dedent: 4
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 :class:`HarmonicCalculator` is coupled with
a calculator that can compute interactions beyond the Harmonic Approximation,
e.g. :mod:`~ase.calculators.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()
|
