# Theoretical formulations

## Contents

### General form of empirical potential functions

- In classical molecular dynamics, the
**total potential energy**[math]U[/math] of a system can be written as the sum of**site potentials**[math]U_i[/math]:

$$U=\sum_{i=1}^N U_i(\{\vec{r}_{ij}\}_{j\neq i}).$$

- Here $$\vec{r}_{ij} \equiv \vec{r}_j - \vec{r}_i$$

is the **position difference vector** used throughout this online manual.

### General form of interatomic forces

- Recently, a well-defined force expression for
**general many-body potentials**that explicitly respects (the weak version of)**Newton's third law**has been derived as [1]:

$$ \vec{F}_{i} = \sum_{j \neq i} \vec{F}_{ij}; $$ $$ \vec{F}_{ij} = - \vec{F}_{ji} = \frac{\partial U_{i}}{\partial \vec{r}_{ij}} - \frac{\partial U_{j}}{\partial \vec{r}_{ji}} = \frac{\partial \left(U_{i} + U_{j}\right) }{\partial \vec{r}_{ij}}. $$

- Here, [math]\partial U_{i}/\partial \vec{r}_{ij}[/math] is a shorthand notation for a vector with cartesian components [math]\partial U_{i}/\partial x_{ij}[/math], [math]\partial U_{i}/\partial y_{ij}[/math], and [math]\partial U_{i}/\partial z_{ij}[/math].

Starting from the above force expression, one can derive expressions for the stress tensor and the heat current.

### Stress tensor

**Stress tensor**is an important quantity in MD simulations. It consists of two parts: a virial part which is related to the force and an ideal-gas part which is related to the temperature. The virial part must be calculated along with force evaluation.

- The validity of Newton's third law is crucial to derive the following expression of the
**virial tensor**[1]:

$$ \textbf{W} = \sum_{i} \textbf{W}_i; $$ $$ \textbf{W}_i = -\frac{1}{2} \sum_{j \neq i} \vec{r}_{ij} \otimes \vec{F}_{ij}, $$ where only relative positions [math]\vec{r}_{ij}[/math] are involved.

- After a little algebra, we can also express the virial as [2]

$$ \textbf{W} = \sum_{i} \textbf{W}_i; $$ $$ \textbf{W}_i = \sum_{j \neq i} \vec{r}_{ij} \otimes \frac{\partial U_j}{\partial \vec{r}_{ji}}. $$

- The ideal-gas part of the stress is isotropic, which is given by the ideal-gas pressure:

$$ p_{\text{ideal}}=\frac{Nk_{\rm B}T}{V}, $$ where [math]N[/math] is the number of particles, [math]k_{\rm B}[/math] is Boltzmann's constant, [math]T[/math] is the absolute temperature, and [math]V[/math] is the volume of the system.

- The total stress tensor is thus

$$ \sigma^{\alpha\beta} = -\frac{1}{2V} \sum_{i} \sum_{j \neq i} r_{ij}^{\alpha} F_{ij}^{\beta} + \frac{Nk_{\rm B}T}{V} \delta^{\alpha\beta}. $$

### Heat current

- Using the force expression, one can derive the following expression for the
**heat current**for the whole system ([math]E_i[/math] is the total energy of atom [math]i[/math]) [1]:

$$ \vec{J} = \vec{J}^{\text{pot}} + \vec{J}^{\text{kin}} = \sum_{i} \vec{J}^{\text{pot}}_i + \sum_{i} \vec{J}^{\text{kin}}_i; $$ $$ \vec{J}^{\text{kin}}_i = \vec{v}_i E_i; $$ $$ \vec{J}^{\text{pot}}_i = -\sum_{j \neq i} \vec{r}_{ij} \left(\frac{\partial U_i}{\partial \vec{r}_{ij}} \cdot \vec{v}_j -\frac{\partial U_j}{\partial \vec{r}_{ji}} \cdot \vec{v}_i\right). $$

- The potential part of the per-atom heat current can also be written in the following equivalent forms [1]:

$$ \vec{J}^{\text{pot}}_i = -\sum_{j \neq i} \vec{r}_{ij} \left(\frac{\partial U_i}{\partial \vec{r}_{ij}} \cdot \vec{v}_j\right); $$ $$ \vec{J}^{\text{pot}}_i = \sum_{j \neq i} \vec{r}_{ij} \left(\frac{ \partial U_j} {\partial \vec{r}_{ji}} \cdot \vec{v}_i\right). $$

- Therefore, the per-atom heat current can also be expressed in terms of the per-atom virial [2]:

$$
\vec{J}^{\text{pot}}_i = \mathbf{W}_i \cdot \vec{v}_i.
$$
where the per-atom virial tensor cannot be assumed to be symmetric and the fully tensor with 9 components should be used [2]. This result has actually been clear from the derivations in Ref. [1], but it was **wrongly stated** in Ref. [1] that the potential part of the heat current **cannot** be expressed in terms of the per-atom virial.

- One can also derive the following expression for the heat current from a subsystem [math]A[/math] to a subsystem [math]B[/math] [3]:

$$ Q_{A \rightarrow B} = -\sum_{i \in A} \sum_{j \in B} \left(\frac{\partial U_i}{\partial \vec{r}_{ij}} \cdot \vec{v}_j -\frac{\partial U_j}{\partial \vec{r}_{ji}} \cdot \vec{v}_i\right). $$

## Integration by one step

- The aim of time evolution is to find the
**phase trajectory**

$$
\{ \vec{r}_i(t_1), ~\vec{v}_{i}(t_1)\}_{i=1}^N,~
\{ \vec{r}_i(t_2), ~\vec{v}_{i}(t_2)\}_{i=1}^N,~
\cdots
$$
starting from the **initial phase point**
$$
\{ \vec{r}_i(t_0), ~\vec{v}_{i}(t_0)\}_{i=1}^N.
$$

- The time interval between two time points [math]\Delta t=t_1-t_0=t_2-t_1=\cdots[/math] is called the
**time step**.

- The algorithm for integrating by one step depends on the ensemble type and other external conditions. There are many ensembles used in MD simulations, but we only consider the following 3 in the current version:

### The [math]NVE[/math] ensemble

- This is also called the
**micro-canonical ensemble**. We use the**velocity-Verlet**integration method with the following equations:

$$ \vec{v}_i(t_{m+1}) \approx \vec{v}_i(t_{m}) + \frac{\vec{F}_i(t_m)+\vec{F}_i(t_{m+1})}{2m_i}\Delta t; $$ $$ \vec{r}_i(t_{m+1}) \approx \vec{r}_i(t_{m}) + \vec{v}_i(t_m) \Delta t + \frac{1}{2} \frac{\vec{F}_i(t_m)}{m_i} (\Delta t)^2. $$

Explanations:

- [math]\vec{v}_i(t_{m})[/math] is the velocity vector of particle [math]i[/math] at time [math]t_{m}[/math].
- [math]\vec{r}_i(t_{m})[/math] is the position vector of particle [math]i[/math] at time [math]t_{m}[/math].
- [math]\vec{F}_i(t_{m})[/math] is the force vector of particle [math]i[/math] at time [math]t_{m}[/math].
- [math]m_i[/math] is the mass of particle [math]i[/math].
- [math]\Delta t[/math] is the time step for integration.

### The [math]NVT[/math] ensemble

This is also called the **canonical ensemble**. We have implemented quite a few thermostats for the NVT ensemble.

#### The Berendsen thermostat

- The velocities are scaled in the
**Berendsen thermostat**[4] in the following way:

$$ \vec{v}_i^{\text{scaled}} = \vec{v}_i \sqrt{1 + \alpha_T \left(\frac{T_0}{T} - 1\right)}. $$

- Here, [math]\alpha_T[/math] is a dimensionless parameter, [math]T_0[/math]is the target temperature, and [math]T[/math] is the instant temperature calculated from the current velocities. The parameter [math]\alpha_T[/math] should be positive and not larger than 1.
- When [math]\alpha_T=1[/math], the above formula reduces to the simple velocity-scaling formula:

$$ \vec{v}_i^{\text{scaled}} = \vec{v}_i \sqrt{\frac{T_0}{T}}. $$

- A smaller [math]\alpha_T[/math] represents a weaker coupling between the system and the thermostat. Practically, any value of [math]\alpha_T[/math] in the range of [math][0.01, 1][/math] can be used.

- The above re-scaling is applied at each time step after the velocity-Verlet integration.

- This thermostat is usually used for reaching equilibrium and is
**not good**for sampling the canonical ensemble.

#### The Nose-Hoover chain thermostat

- In the
**Nose-Hoover chain**method, the equations of motion for the particles in the thermostatted region are (those for the thermostat variables are not presented):

$$ \frac{d \vec{r}_i}{dt} = \frac{\vec{p}_i}{m_i}; $$ $$ \frac{d \vec{p}_i}{dt} = \vec{F}_i - \frac{\pi_0}{Q_0} \vec{p}_i. $$

- Here,
- [math]\vec{r}_i[/math] is the position of atom [math]i[/math].
- [math]\vec{p}_i[/math] is the momentum of atom [math]i[/math].
- [math]m_i[/math] is the mass of atom [math]i[/math].
- [math]\vec{F}_i[/math] is the total force on atom [math]i[/math] resulted from the potential model used.
- [math]Q_0=N_{\rm f} k_{\rm B} T \tau^2[/math] is the
*mass*of the thermosat variable directly coupled to the system and [math]\pi_0[/math] is the corresponding*momentum*.- [math]N_{\rm f}[/math] is the degree of freedom in the thermostatted region.
- [math]k_{\rm B}[/math] is Boltzmann's constant and [math]T[/math] is the target temperature.

- [math]\tau[/math] is a time parameter, and we suggest a value from 0.01 ps to 1 ps.

- We implemented the integrator as described in Tuckerman's book [5].
- We used a fixed chain length of 4.

#### The Langevin thermostat

- In the
**Langevin**method, the equations of motion for the particles in the thermostatted region are

$$ \frac{d \vec{r}_i}{dt} = \frac{\vec{p}_i}{m_i}; $$ $$ \frac{d \vec{p}_i}{dt} = \vec{F}_i - \frac{\vec{p}_i}{\tau} + \vec{f}_i, $$

- Here,
- [math]\vec{r}_i[/math] is the position of atom [math]i[/math].
- [math]\vec{p}_i[/math] is the momentum of atom [math]i[/math].
- [math]m_i[/math] is the mass of atom [math]i[/math].
- [math]\vec{F}_i[/math] is the total force on atom [math]i[/math] resulted from the potential model used.
- [math]\vec{f}_i[/math] is a random force with a variation determined by the fluctuation-dissipation relation to recover the canonical ensemble distribution with the target temperature.
- [math]\tau[/math] is a time parameter. We suggest a value from 0.01 ps to 1 ps.

- We implemented the integrator proposed by Bussi and Parinello.

#### The Bussi-Donadio-Parrinello thermostat

- The Berendsen thermostat does not generate a true NVT ensemble. As an extension of the Berendsen thermostat, the
**Bussi-Donadio-Parrinello**(BDP) thermostat incorporates a proper randomness into the velocity re-scaling factor and generates a true NVT ensemble.

- In the BDP thermostat, the velocities are scaled in the following way:

$$ \vec{v}_i^{\text{scaled}} =\alpha \vec{v}_i; $$ $$ \alpha^2= e^{-\Delta t/\tau} + \frac{T_0}{TN_f} \left( 1-e^{-\Delta t/\tau} \right) \left( R_1^2 + \sum_{i=2}^{N_f}R_i^2 \right) + 2e^{-\Delta t/2\tau} R_1 \sqrt{\frac{T_0}{TN_f} \left( 1-e^{-\Delta t/\tau} \right) }. $$

- Here,
- [math]\vec{v}_i[/math] is the velocity of atom [math]i[/math] before the re-scaling.
- [math]N_{\rm f}[/math] is the degree of freedom in the thermostatted region.
- [math]T[/math] is instant temperature and [math]T_0[/math] is the target temperature.
- [math]\Delta t[/math] is the time step for integration.
- [math]\tau[/math] is a time parameter, and we suggest a value from 0.01 ps to 1 ps.
- [math]\{R_i\}_{i=1}^{N_f}[/math] are [math]N_{\rm f}[/math] Gaussian distributed random numbers with zero mean and unit variance.

- We used the functions provided by Bussi to calculate the stochastic velocity re-scaling factor [math]\alpha[/math].

- We have not thoroughly tested this thermostat yet.

### The [math]NPT[/math] ensemble

- There seems to be no simple name for this important ensemble, but it is usually called the isothermal-isobaric ensemble. We have only implemented the
**Berendsen barostat**[4].

- In the Berendsen barostat algorithm, the particle positions and box length in a given direction are scaled if periodic boundary conditions are applied to that direction.
- For example, the scaling of the positions in the [math]x[/math] direction reads

$$ x_i^{\text{scaled}} = x_i \left[ 1-\alpha_p (p_{x0} - p_x) \right]. $$

- Here, [math]p_{x0}[/math] ([math]p_{x}[/math]) is the target (instant) pressure in the [math]x[/math] direction.
- The parameter [math]\alpha_p[/math] is not dimensionless, and it requires some trial-and-error to find a good value of it for a given system. A harder/softer system requires a smaller/larger value of [math]\alpha_p[/math]. In the unit system adopted by
`GPUMD`

, it is recommended that [math]\alpha_p = 10^{-4} \sim 10^{-2}[/math]. - Only directions with periodic boundary conditions will be affected by the barostat. The target pressure in a non-periodic direction is thus irrelevant, although the code requires you to set one. You can just set it to zero.

## Heat transport

- GPUMD has many utilities to study heat transport, including

### The EMD method for heat transport

- In MD simulations, a popular approach of computing the lattice thermal conductivity is to use the Green-Kubo formula. In this method, the running lattice thermal conductivity along the [math]x[/math]-direction (similar expressions apply to other directions) can be expressed as an integral of the heat current autocorrelation (HAC):

$$ \kappa_{xx}(t) = \frac{1}{k_BT^2V} \int_0^{t} dt' \text{HAC}_{xx}(t'). $$

- Here, [math]k_{\rm B}[/math] is Boltzmann's constant, [math]V[/math] is the volume of the simulated system, [math]T[/math] is the absolute temperature, and [math]t[/math] is the correlation time. The HAC is

$$ \text{HAC}_{xx}(t)=\langle J_{x}(0)J_{x}(t)\rangle, $$ where [math]J_{x}(0)[/math] and[math]J_{x}(t)[/math] are the total heat current of the system at two time points separated by an interval of [math]t[/math]. The symbol [math]\langle \rangle[/math] means that the quantity inside will be averaged over different time origins.

- Related keyword in the run.in file: compute_hac

- Related output file: hac.out

- We only used the potential part of the heat current. If you are studying fluids, you need to output the heat currents (potential and kinetic part) using the compute keyword and calculated the HACs by yourself.

- We have decomposed the potential part of the heat current into in-plane and out-of-plane components [3]. If you do not need this decomposition, you can simply sum up some components in the hac.out file.

### The NEMD method for heat transport

- Nonequilibrium molecular dynamics (NEMD) can be used to study thermal transport. In this method, two local thermostats at different temperatures are used to generate a nonequilibrium steady state with a constant heat flux.

- If the temperature difference between the two thermostats is [math]\Delta T[/math] and the heat flux is [math]Q/S[/math], the
**thermal conductance**[math]G[/math] between the two thermostats can be calculated as

$$ G=\frac{Q/S}{\Delta T}. $$ Here, [math]Q[/math] is the energy transfer rate between the thermostat and the thermostated region and [math]S[/math] is the cross-sectional area perpendicular to the transport direction.

- We can also calculate an
**effective thermal conductivity**[math]\kappa(L)[/math] for the finite system:

$$ \kappa(L) = GL = \frac{Q/S}{\Delta T/L}. $$ where [math]L[/math] is the length between the heat source and the heat sink. This is to say that the temperature gradient should be calculated as [math]\Delta T/L[/math], rather than that extracted from the linear part of the temperature profile away from the local thermostats. This is an important conclusion in Ref. [6].

- To generate the nonequilibrium steady state, one can use a pair of local thermostats. Base on Ref. [6], the Langevin thermostatting method is recommended. Therefore, the ensemble keyword with the first parameter of
`heat_lan`

should be used to generate the heat current.

- The compute keyword should be used to compute the temperature profile and the heat transfer rate [math]Q[/math].

- Related output file: compute.out

### The HNEMD method for heat transport

- The homogeneous nonequilibrium molecular dynamics (HNEMD) method for heat transport by Evans has been recently generalized to general many-body potentials [7]. This method is physically equivalent to the Green-Kubo (EMD) method but is computationally much faster.

- In this method, an external force of the form [7]

$$ \vec{F}_{i}^{\rm ext} = E_i \vec{F}_{\rm e} + \sum_{j \neq i} \left(\frac{\partial U_j}{\partial \vec{r}_{ji}} \otimes \vec{r}_{ij}\right) \cdot \vec{F}_{\rm e} $$ is added to each atom [math]i[/math], driving the system out of equilibrium. According to Ref. [2], it can also be written as $$ \vec{F}_{i}^{\rm ext} = E_i \vec{F}_{\rm e} + \vec{F}_{\rm e} \cdot \mathbf{W}_i $$

- Here,
- [math]E_i[/math] is the total energy of particle [math]i[/math].
- [math]U_i[/math] is the potential energy of particle [math]i[/math].
- [math]\vec{r}_{ij}\equiv\vec{r}_{j}-\vec{r}_{i}[/math], and [math]\vec{r}_i[/math] is the position of particle [math]i[/math].
- The parameter [math]\vec{F}_{\rm e}[/math] is of the dimension of inverse length and should be small enough to keep the system within the linear response regime.

- The driving force will induce a nonequilibrium heat current [math]\langle \vec{J} \rangle_{\rm ne}[/math] linearly related to [math]\vec{F}_{\rm e}[/math]:

$$ \frac{\langle J^{\mu}(t)\rangle_{\rm ne}}{TV} = \sum_{\nu} \kappa^{\mu\nu} F^{\nu}_{\rm e}, $$ where [math]\langle J^{\mu}(t)\rangle_{\rm ne}[/math] is the thermal conductivity tensor, [math]T[/math] is the system temperature, and [math]V[/math] is the system volume.

- A global thermostat should be applied to control the temperature of the system. For this, we recommend using the Nose-Hoover chain thermostat. So one should use the ensemble keyword with the first parameter of
`nvt_nhc`

. - The compute_hnemd keyword should be used to add the driving force and calculate the thermal conductivity.
- The computed results are saved to the kappa.out fiel.

### Spectral heat current

- In the framework of the NEMD and HNEMD methods, one can also calculate spectrally decomposed thermal conductivity (or conductance). In this method, one first calculates the following virial-velocity correlation function [2]:

$$ \vec{K}(t) = \sum_{i} \sum_{j \neq i} \left\langle \mathbf{W}_i(0) \cdot \vec{v}_i (t) \right\rangle, $$ which reduces to the nonequilibrium heat current when [math]t=0[/math].

- Then one can define the following Fourier transform pairs [2]:

$$ \tilde{\vec{K}}(\omega) = \int_{-\infty}^{\infty} dt e^{i\omega t} K(t) $$ $$ \vec{K}(t) = \int_{-\infty}^{\infty} \frac{d\omega}{2\pi} e^{-i\omega t} \tilde{\vec{K}}(\omega) $$

- By setting [math]t=0[/math] in the equation above, we can get the following spectral decomposition of the nonequilibrium heat current:

$$ \vec{J} = \int_{0}^{\infty} \frac{d\omega}{2\pi} \left[2\tilde{\vec{K}}(\omega)\right]. $$

- From the spectral decomposition of the nonequilibrium heat current, one can deduce the spectrally decomposed thermal conductance in the NEMD method:

$$ G(\omega) = \frac{2\tilde{\vec{K}}(\omega)}{V\Delta T} \quad \text{with} \quad G = \int_{0}^{\infty} \frac{d\omega}{2\pi} G(\omega). $$ where [math]\Delta T[/math] is the temperature difference between the two thermostats and [math]V[/math] is the volume of the considered system or subsystem.

- One can also calculate the spectrally decomposed thermal conductivity in the HNEMD method:

$$ \kappa(\omega) = \frac{2\tilde{\vec{K}}(\omega)}{VTF_{\rm e}} \quad \text{with} \quad \kappa = \int_{0}^{\infty} \frac{d\omega}{2\pi} \kappa(\omega). $$ where [math]F_{\rm e}[/math] is the magnitude of the driving force parameter in the HNEMD method.

- This calculation is invoked by the compute_shc keyword and the results are saved to the shc.out file.

- This method seems to be only applicable to special systems. So far, we have only applied it to graphene and carbon nanotubes. We have obtained wired results in some other systems and are still not clear about the reason.

## References

- [1] Zheyong Fan, Luiz Felipe C. Pereira, Hui-Qiong Wang, Jin-Cheng Zheng, Davide Donadio, and Ari Harju,
*Force and heat current formulas for many-body potentials in molecular dynamics simulations with applications to thermal conductivity calculations*, Phys. Rev. B**92**, 094301 (2015). - [2] Alexander J. Gabourie and Zheyong Fan, In preparation.
- [3] Zheyong Fan, Luiz Felipe C. Pereira, Petri Hirvonen, Mikko M. Ervasti, Ken R. Elder, Davide Donadio, Tapio Ala-Nissila, and Ari Harju,
*Thermal conductivity decomposition in two-dimensional materials: Application to graphene*, Phys. Rev. B**95**, 144309 (2017). - [4] H. J. C. Berendsen, J. P. M. Postma, W. F. van Gunsteren, A. DiNola, and J. R. Haak,
*Molecular dynamics with coupling to an external bath*, J. Chem. Phys.**81**, 3684 (1984). - [5] Mark E. Tuckerman,
*Statistical Mechanics: Theory and Molecular Simulation (Oxford Graduate Texts)*, 1st Edition, Oxford University Press (2010). - [6] Zhen Li, Shiyun Xiong, Charles Sievers, Yue Hu, Zheyong Fan, Ning Wei, Hua Bao, Shunda Chen, Davide Donadio, and Tapio Ala-Nissila,
*Influence of Thermostatting on Nonequilibrium Molecular Dynamics Simulations of Heat Conduction in Solids*, J. Chem. Phys.**151**, 234105 (2019). - [7] Zheyong Fan, Haikuan Dong, Ari Harju, and Tapio Ala-Nissila,
*Homogeneous nonequilibrium molecular dynamics method for heat transport and spectral decomposition with many-body potentials*, Phys. Rev. B**99**, 064308 (2019).