Thermodynamic integration with harmonic reference: Difference between revisions

From VASP Wiki
No edit summary
No edit summary
 
(36 intermediate revisions by 2 users not shown)
Line 1: Line 1:
The Helmholtz free energy (<math>A</math>) of a fully interacting system (1) can be expressed in terms of that of harmonic system (0) as follows
The Helmholtz free energy (<math>A</math>) of a fully interacting system (1) can be expressed in terms of that of system harmonic in Cartesian coordinates (0,<math>\mathbf{x}</math>) as follows
:<math>
:<math>
A_{1} = A_{0} + \Delta A_{0\rightarrow 1}
A_{1} = A_{0,\mathbf{x}} + \Delta A_{0,\mathbf{x}\rightarrow 1}
</math>  
</math>  
where <math>\Delta A_{0\rightarrow 1}</math> is anharmonic free energy. The latter term can be determined by means of thermodynamic integration (TI)
where <math>\Delta A_{0,\mathbf{x}\rightarrow 1}</math> is anharmonic free energy. The latter term can be determined by means of thermodynamic integration<ref>[https://pubs.aip.org/aip/jcp/article-abstract/3/5/300/203696/Statistical-Mechanics-of-Fluid-Mixtures?redirectedFrom=fulltext J. G. Kirkwood, ''J. Chem. Phys''. '''3,''' 300–313 (1935)]</ref> (TI)
:<math>
:<math>
\Delta A_{0\rightarrow 1} = \int_0^1 d\lambda \langle V_1 -V_0 \rangle_\lambda
\Delta A_{0,\mathbf{x}\rightarrow 1} = \int_0^1 d\lambda \langle V_1 -V_{0,\mathbf{x}} \rangle_\lambda
</math>
</math>
with <math>V_i</math> being the potential energy of system <math>i</math>, <math>\lambda</math> is a coupling constant and <math>\langle\cdots\rangle_\lambda</math> is the NVT ensemble average of the system driven by the Hamiltonian
with <math>V_i</math> being the potential energy of system <math>i</math>, <math>\lambda</math> is a coupling constant and <math>\langle\cdots\rangle_\lambda</math> is the NVT ensemble average of the system driven by the Hamiltonian
:<math>
:<math>
\mathcal{H}_\lambda = \lambda \mathcal{H}_1 + (1-\lambda)\mathcal{H}_0
\mathcal{H}_\lambda = \lambda \mathcal{H}_1 + (1-\lambda)\mathcal{H}_{0,\mathbf{x}}
</math>
</math>


Line 19: Line 19:
configuration corresponding to the potential energy minimum with the  
configuration corresponding to the potential energy minimum with the  
atomic position vector <math>\mathbf{x}_0</math>,
atomic position vector <math>\mathbf{x}_0</math>,
the number of vibrational degrees of freedom <math>N_\mathrm{vib}</math>, and the angular frequency <math>\omega_i</math> of vibrational mode <math>i</math>.
the number of vibrational degrees of freedom <math>N_\mathrm{vib}</math>, and the angular frequency <math>\omega_i</math> of vibrational mode <math>i</math> obtained using the Hesse matrix <math>\underline{\mathbf{H}}^\mathbf{x}</math>.
The
Finally, the harmonic potential energy is expressed as
:<math>
:<math>
V_{0,\mathbf{x}}(\mathbf{x}) = V_{0,\mathbf{x}}(\mathbf{x}_0) + \frac{1}{2} (\mathbf{x} - \mathbf{x}_0)^T \underline{\mathbf{H}}^\mathbf{x} (\mathbf{x} - \mathbf{x}_0)
V_{0,\mathbf{x}}(\mathbf{x}) = V_{0,\mathbf{x}}(\mathbf{x}_0) + \frac{1}{2} (\mathbf{x} - \mathbf{x}_0)^T \underline{\mathbf{H}}^\mathbf{x} (\mathbf{x} - \mathbf{x}_0)
</math>
</math>
Thus, a conventional TI calculation consists of  the following steps:
# determine <math>\mathbf{x}_0</math> and <math>V_{0,\mathbf{x}}(\mathbf{x}_0)</math> in structural relaxation
# compute <math>\omega_i</math> in vibrational analysis
# use the data obtained in the point 2 to determine <math>\underline{\mathbf{H}}^\mathbf{x}</math> that defines the harmonic forcefield
# perform NVT MD simulations for several values of  <math>\lambda \in \langle0,1\rangle</math> and determine <math>\langle V_1 -V_{0,\mathbf{x}}  \rangle</math>
# integrate <math>\langle V_1 -V_{0,\mathbf{x}}  \rangle</math> over the <math>\lambda </math> grid and compute <math>\Delta A_{0,\mathbf{x}\rightarrow 1}</math>
Unfortunately, there are several problems linked with such a straightforward approach. First, the systems with rotational and/or translational degrees of freedom cannot be treated in a straightforward manner because <math>V_{0,\mathbf{x}}(\mathbf{x})</math> is not invariant under rotations and translations. Conventional TI is thus unsuitable for simulations of gas phase molecules or adsorbate-substrate systems. and this problem also imposes restrictions on the choice of thermostat used in NVT simulation (Langevin thermostat, for instance, does not conserve position of the center of mass and is therefore unsuitable for the use in conventional TI). Furthermore, if the Hesse matrix of the harmonic system has one or more eigenvalues that nearly vanish, the simulations with <math>\lambda \rightarrow </math> 0 is likely to generate unphysical configurations causing serious convergence issues. These problems have been addressed in series of works by Amsler et al.<ref>[https://pubs.acs.org/doi/abs/10.1021/acs.jctc.0c01022 J. Amsler, P. N. Plessow, F. Studt, T. Bucko, J. Chem. Theory  Comput. 17, 1155-1169 (2021)]</ref><ref>[https://pubs.acs.org/doi/abs/10.1021/acs.jctc.3c00169 J. Amsler, P. N. Plessow, F. Studt, T. Bučko, J. Chem. Theory Comput. 19, 2455-2468 (2023)]</ref>
First, the method was formulated in terms of rotationally and translationally invariant internal coordinates <math>\mathbf{q}=\mathbf{q}(\mathbf{x})</math>, whereby the free energy of interacting system is repartitioned as follows:
:<math>
A_1 = A_{0,\mathbf{x}} + \Delta A_{0,\mathbf{x} \rightarrow 0,\mathbf{q}} + \Delta A_{0,\mathbf{q} \rightarrow 1}
</math>
where <math>\Delta A_{0,\mathbf{x} \rightarrow 0,\mathbf{q}}</math> is the free energy change due to transformation from the system harmonic in <math>\mathbf{x}</math> into the system harmonic in <math>\mathbf{q}</math> and <math>\Delta A_{0,\mathbf{q} \rightarrow 1} </math> is that for the transformation of the latter into a fully interacting system. The force field for the system harmonic in <math>\mathbf{q}</math> is defined as:
:<math>
V_{0,\mathbf{q}}(\mathbf{q}) = V_{0,\mathbf{q}}(\mathbf{q}_0) + \frac{1}{2} (\mathbf{q} - \mathbf{q}_0)^T \mathbf{\underline{H}^q} (\mathbf{q} - \mathbf{q}_0)
</math>
where <math>\mathbf{q}_0=\mathbf{q}(\mathbf{x}_0)</math> and the Hesse matrix <math>\mathbf{\underline{H}}^\mathbf{q}</math> defined for a potential energy minimum is related to <math>\mathbf{\underline{H}}^\mathbf{x}</math> via
<math>
\mathbf{\underline{H}}^\mathbf{x} = \mathbf{\underline{B}}^T \mathbf{\underline{H}}^\mathbf{q} \mathbf{\underline{B}}
</math>
with
<math>
\mathbf{\underline{B}}_{i,j} = \frac{\partial q_i}{\partial x_j}
</math>
being the Wilson matrix. Note that the calculation of the term <math>\Delta A_{0,\mathbf{x} \rightarrow 0,\mathbf{q}}</math> is inexpensive as it corresponds to a force field to force field transformation. Furthermore, this term vanishes in the case of phase volume conserving coordinates, such as interatomic distances.
The TI calculations in internal coordinates are performed in NVT ensemble using any thermostat available in VASP. The coupling parameter <math>\lambda</math> is defined by setting the parameter {{TAG|TI_LAMBDA}} in the [[INCAR|INCAR]] file.
The set of internal coordinates used in the TI calculation are defined via the {{FILE|ICONST}} file by setting the status to 3. The Hesse matrix <math>\mathbf{\underline{H}}^\mathbf{x}</math> is provided in the file  {{FILE|HESSEMAT}} and its transformation into <math>\mathbf{\underline{H}}^\mathbf{q}</math> is performed by VASP. The potential energies of the system 1 and 0,<math>\mathbf{q}</math>, needed to compute  <math>\langle V_1 -V_{0,\mathbf{q}}  \rangle</math> used as integrant in the TI expression for <math>\Delta A_{0,\mathbf{q} \rightarrow 1} </math>, are written in the file {{FILE|REPORT}} in lines introduce by a string "e_ti>"
== References ==
[[Category:Molecular dynamics]]

Latest revision as of 06:30, 8 May 2024

The Helmholtz free energy () of a fully interacting system (1) can be expressed in terms of that of system harmonic in Cartesian coordinates (0,) as follows

where is anharmonic free energy. The latter term can be determined by means of thermodynamic integration[1] (TI)

with being the potential energy of system , is a coupling constant and is the NVT ensemble average of the system driven by the Hamiltonian

Free energy of harmonic reference system within the quasi-classical theory writes

with the electronic free energy for the configuration corresponding to the potential energy minimum with the atomic position vector , the number of vibrational degrees of freedom , and the angular frequency of vibrational mode obtained using the Hesse matrix . Finally, the harmonic potential energy is expressed as

Thus, a conventional TI calculation consists of the following steps:

  1. determine and in structural relaxation
  2. compute in vibrational analysis
  3. use the data obtained in the point 2 to determine that defines the harmonic forcefield
  4. perform NVT MD simulations for several values of and determine
  5. integrate over the grid and compute

Unfortunately, there are several problems linked with such a straightforward approach. First, the systems with rotational and/or translational degrees of freedom cannot be treated in a straightforward manner because is not invariant under rotations and translations. Conventional TI is thus unsuitable for simulations of gas phase molecules or adsorbate-substrate systems. and this problem also imposes restrictions on the choice of thermostat used in NVT simulation (Langevin thermostat, for instance, does not conserve position of the center of mass and is therefore unsuitable for the use in conventional TI). Furthermore, if the Hesse matrix of the harmonic system has one or more eigenvalues that nearly vanish, the simulations with 0 is likely to generate unphysical configurations causing serious convergence issues. These problems have been addressed in series of works by Amsler et al.[2][3]

First, the method was formulated in terms of rotationally and translationally invariant internal coordinates , whereby the free energy of interacting system is repartitioned as follows:

where is the free energy change due to transformation from the system harmonic in into the system harmonic in and is that for the transformation of the latter into a fully interacting system. The force field for the system harmonic in is defined as:

where and the Hesse matrix defined for a potential energy minimum is related to via with being the Wilson matrix. Note that the calculation of the term is inexpensive as it corresponds to a force field to force field transformation. Furthermore, this term vanishes in the case of phase volume conserving coordinates, such as interatomic distances.

The TI calculations in internal coordinates are performed in NVT ensemble using any thermostat available in VASP. The coupling parameter is defined by setting the parameter TI_LAMBDA in the INCAR file. The set of internal coordinates used in the TI calculation are defined via the ICONST file by setting the status to 3. The Hesse matrix is provided in the file HESSEMAT and its transformation into is performed by VASP. The potential energies of the system 1 and 0,, needed to compute used as integrant in the TI expression for , are written in the file REPORT in lines introduce by a string "e_ti>"

References