Molecular Dynamics

In this chapter performing molecular dynamics simulations using MAPP will be discussed.

Supported Ensembles,dt) NVT ensemble,T,dt) N\mathbf{S}T ensemble,T,dt,gas_element,seed) \mu VT ensemble

NHC Thermostat

The thermostat that is employed in MAPP is the well-known Nose-Hoover chains (NHC). Here we will touvh on the equations of motions of NHC breifly and its implementation in MAPP in order to clarify the meaning of some of the variables that users have control over, namely L, t_{\mathrm{relax}}, and n_{\mathrm{iters}}. Interested reader might refer to [1] and [2] for a full discription and mathematical formulation of of NHC governing equations. With that being said the main idea behind NHC is as follows: In order to simulate the coupling of a system with an external thermal bath, we couple its kinetic energy with kinetic enery of a series of fictitious object with fictitious masses, namely links in NHC chain. With index i denoting i th link in the chain.


Schematic representation of NHC of length L .

Suppose we have a general d dimeensional system with N_f deggrees of freedom (they can be positions of atoms or anything else), namely y_0, \cdots, y_{N_f -1}. The equations of motion are as follows:

\ddot{y}_i&=-\frac{1}{m_i}\frac{\partial U}{\partial y_i}-\dot{y}_i\dot{\xi}_0 \quad  0 \le i <N_f\\
\ddot{\xi}_0&=\frac{N_f k_B}{Q_0} \left(T_{\mathrm{int}}-T_{\mathrm{ext}}\right)-\dot{\xi}_0\dot{\xi}_1\\
\ddot{\xi}_i&=\frac{Q_{i-1}}{Q_i}\dot{\xi}_{i-1}^2-\frac{k_BT_{\mathrm{ext}}}{Q_i}-\dot{\xi}_i\dot{\xi}_{i+1} \quad 0 < i < L-1\\

where U is the potential energy, and m_i is the mass associated i th degree of freedom; \xi_i and Q_i are degree of freedom of and mass of i th link in Nose-Hoover chain, respectively. Here \dot{\square} denotes derivative with respect to time. T_{\mathrm{ext}} is the external (bath) temperature; T_{\mathrm{int}} is the internal (current) temperature of system, calculated using

T_{\mathrm{int}}=\frac{1}{k_B N_f}\sum_{i=0}^{N_f-1} m_i\dot{y}_i^2

As it is apparent from the equations above \xi_i is uinitless and Q_i is of unit of \mathrm{ML^2}. Let us intorduce the following change of variable

\theta_i&=\frac{Q_i}{Q_{i+1}}  \quad 0 < i < L-1\\
t_{\mathrm{relax}}^2 &= \frac{Q_0}{N_f k_B T_{\mathrm{ext}}}

where \theta_i is a unitless parameter and t_{\mathrm{relax}} is a positive time scale, which in turn simplify thermostat equations of motion:

\ddot{\xi}_0&=t_{\mathrm{relax}}^{-2} \left(\frac{T_{\mathrm{int}}}{T_{\mathrm{ext}}}-1\right)-\dot{\xi}_0\dot{\xi}_1\\
\ddot{\xi}_i&=\theta_{i-1}\dot{\xi}_{i-1}^2-\frac{t_{\mathrm{relax}}^{-2}}{N_f}\prod_{j=0}^{i-1}\theta_j-\dot{\xi}_i\dot{\xi}_{i+1} \quad 0 < i < L-1\\

in MAPP as a matter of convention \theta_i are chosen as follows:

\theta_0=N_f, \quad \theta_i=1,\quad 0 < i < L,

which gives us the final form of equations of motion implemented in MAPP:

\ddot{\xi}_0&=t_{\mathrm{relax}}^{-2} \frac{T_{\mathrm{int}}}{T_{\mathrm{ext}}}-t_{\mathrm{relax}}^{-2}-\dot{\xi}_0\dot{\xi}_1\\
\ddot{\xi}_1&=N_f\dot{\xi}_{0}^2-t_{\mathrm{relax}}^{-2}-\dot{\xi}_1\dot{\xi}_{2} \\
\ddot{\xi}_i&=\dot{\xi}_{i-1}^2-t_{\mathrm{relax}}^{-2}-\dot{\xi}_i\dot{\xi}_{i+1} \quad 1 < i < L-1\\

For completeness here are the equations of motions for L=2:

\ddot{\xi}_0&=t_{\mathrm{relax}}^{-2} \frac{T_{\mathrm{int}}}{T_{\mathrm{ext}}}-t_{\mathrm{relax}}^{-2}-\dot{\xi}_0\dot{\xi}_1\\

and L=1:

\ddot{\xi}_0&=t_{\mathrm{relax}}^{-2} \frac{T_{\mathrm{int}}}{T_{\mathrm{ext}}}-t_{\mathrm{relax}}^{-2}

The two parameters t_{\mathrm{relax}} and L can be adjusted by user through the following attributes

Another variable that can be adjusted is n_{\mathrm{iters}}, which is related to the numerical time integration. In order to increase the percision of numerical integration instead of eveolving the system through one step of \Delta t, it can be performed through n_{\mathrm{iters}} steps of \Delta t / n_{\mathrm{iters}}. n_{\mathrm{iters}} can be adjusted by user through the following attributes


Isothermal-isostress ensemble that is used in mapp is taken from [3]. Aside from external stress, temperature and thermostat related parameters (see NHC Thermostat), there is one other matrix parameter (\mathbf{\tau} ) that user has control over. Before going over the equations of motion it should be noted that our formultion is different from orginal presentation in three major ways:

  • stress is coupled to a separate thermostat rather than particle thermostat
  • unitcell matrix (\mathbf{H} ) is the transpose of the one suggested in the original formulation (\mathbf{h} )
  • unlike the original formulation unitcell mass is not scalar i.e. each of \mathbf{H} components have their own mass, and are defined using a time scale matrix \mathbf{\tau}

Cosider a d dimensional unitcell defined by \mathbf{H} tensor, containing N atoms, interacting via a given potential energy, where \mathbf{x}_i denotes the position of atom i. The particles in the system have N_f degrees of freedom (which is not necessarily equal to dN ) and are coupled with NHC thermostat. The system is subjected to an external temperature of T_{\mathrm{ext}} and an extrenal stress tensor \mathbf{S}_{\mathrm{ext}}. The equations of motions regarding the particles are as follows


here \dot{\square} represents derivative with respect to time. Please note that here \mathbf{v}_i is not the actual particle trajectory but rather \dot{\mathbf{x}}_i is. \mathbf{f}_i is the force exerted on atom i due to potential energy and m_i is its mass. \mathbf{V}^H is a tensor representing the equivalent of velocity for \mathbf{H}; \xi_0 is the particle thermostat related parameter, see NHC Thermostat. The equations of motion pertaining to \mathbf{H} are:

\dot{\mathbf{V}}^H_{ij}&=\frac{ \mathbf{A}_{ij}}{k_B T_{\mathrm{ext}} \mathbf{\tau}_{ij}^2}  -\dot{\xi}_0^H\mathbf{V}^H_{ij}, \quad 0 \le i ,j  <d\\

where \mathbf{\tau} is a given length scale matrix; \xi_0^H is the stress thermostat related parameter (see NHC Thermostat), and

\mathbf{A}&=v\left(s_{\mathrm{ext}}\mathbf{I}-\mathbf{S}_{\mathrm{int}}\right)+v_{\mathrm{ref}}\mathbf{H^T}\mathbf{H_{\mathrm{ref}}^{-T}}\mathbf{S}_{\mathrm{ext}}^{\mathrm{dev}}\mathbf{H_{\mathrm{ref}}^{-1}}\mathbf{H}+k_B T_{\mathrm{int}}\mathbf{I}

where v is the current volume (determinant of \mathbf{H}); \mathbf{H}_{\mathrm{ref}} and v_{\mathrm{ref}} describe the reference configuration and is its volume, respectively; and

s_{\mathrm{ext}}=\frac{\mathrm{Tr}[\mathbf{S}_{\mathrm{ext}}]}{d}, \quad \mathbf{S}_{\mathrm{ext}}^{\mathrm{dev}}=\mathbf{S}_{\mathrm{ext}}-s_{\mathrm{ext}}\mathbf{I}


[1]Mark E. Tuckerman. Statistical mechanics : theory and molecular simulation. Oxford graduate texts. Oxford ; New York : Oxford University Press, 2010., 2010. ISBN 978-0-19-852526-4.
[2]Daan Frenkel and Berend Smit. Understanding molecular simulation : from algorithms to applications. Computational science: v. 1. San Diego, Calif. ; London : Academic, c2002., 2002. ISBN 978-0-12-267351-1.
[3]W. Shinoda, M. Shiga, and M. Mikami. Rapid estimation of elastic constants by molecular dynamics simulation under constant stress. Physical Review B (Condensed Matter and Materials Physics), 69(13):134103–1–8, April 2004. doi:10.1103/PhysRevB.69.134103.