Manual 05 for MD++
Molecular Dynamics Simulations

Keonwook Kang and Wei Cai

[First written, [ Jan 19 ]], 2007

[Last Modified, [ Mar 03 ]], 2008

The case studies described in Manuals 02-04 only deal with atomic positions as degrees of freedom. For example, in Manual 02, we learned how to create a perfect crystal by MD++ command makecrystal. In Manual 03, we varied the size of the perfect crystal to find the equilibrium lattice constants. In Manual 04, we introduced a vacancy to the perfect crystal and let the atomic positions relax to a local energy minimum. In this manual, we will learn how to perform finite temperature, Molecular Dynamics (MD) simulations. For this purpose, we will need to deal with atomic velocities as degrees of freedom, as well as a number of control variables, such as integrator type, time step, and output controls. We will learn to specify the simulation ensembles from available choices: microcanonical (NVE), canonical (NVT), isoenthalpic-isobaric (NPH) and constant pressure/temperature (NPT) ensembles. Different quantities are conserved or controlled in different ensembles. For example, the total number of atoms (N), volume (V) and energy (E) are conserved in the (NVE) ensemble.

It would be a good start explaining MD simulation of NVE ensemble to understand what you do with MD simulations. Imagine a system of N-particle gas inside of a rigid container which is well insulated that the system is isolated from the environment. There is no heat and mass transfer through the container wall. Then you expect that the number of particles, the container volume and the total energy will not change. In MD simulations of NVE ensemble, you do not see a visible wall but you may think of a closed system with fictitious walls.

Now I give the example script mo_NVE.tcl below, written in Tcl version. You can run the script by typing

$ bin/fs_gpp scripts/mo_NVE.tcl

and it ﬁrst creates a perfect BCC crystal of molybdenum (Mo) with a supercell of 5[1 0 0] X 5[0 1 0] X 5[0 0 1] and then runs MD simulations using the Finnis-Sinclair (FS) potential.

# -*-shell-script -*-
# run NVE MD simulation of perfect crystal of Mo.
#
#*******************************************
# Definition of procedures
#*******************************************
proc initmd { } { MD++ {
#setnolog
setoverwrite
dirname = runs/mo-example
zipfiles = 1   # zip output files
#
#--------------------------------------------
# Create Perfect Crystal
element0 = Mo
crystalstructure = body-centered-cubic
latticeconst = 3.1472           # in Angstrom for Mo
latticesize = [ 1 0 0 5
                0 1 0 5
                0 0 1 5]
}}
#------------------------------------------------------------
proc readMoPot { } { MD++ {
# Read the potential file
potfile = ~/Codes/MD++/potentials/mo_pot readpot
}}
#-------------------------------------------------------------
proc openwindow { } { MD++ {
# Plot Configuration
atomradius = 1.0 bondradius = 0.3 bondlength = 2.725
atomcolor = orange highlightcolor = purple
bondcolor = red backgroundcolor = gray70
plotfreq = 10 rotateangles = [ 0 0 0 1.25 ]
openwin alloccolors rotate saverot eval plot
}}
#--------------------------------------------
proc setup_md { } { MD++ {
#MD settings
ensemble_type = "NVE" integrator_type = "VVerlet"
T_OBJ = 300        # (in Kelvin) Desired Temperature
atommass = 95.94   # (in g/mol)
timestep = 0.001   # (in ps)
totalsteps = 2000
# save property every 10 steps
saveprop = 1 savepropfreq = 10
savecn = 0   savecnfreq = 100
writeall = 1 DOUBLE_T = 1 
randseed = 12345 srand48 #randomize random number generator

#*******************************************
# Main program starts here
#*******************************************
initmd
readMoPot
MD++ makecrystal finalcnfile = perf.cn writecn
openwindow
#---------------------------------------------
# run MD
setup_md
MD++ initvelocity finalcnfile = init.cn writecn
#format of thermodynamic property file
MD++ {output_fmt = "curstep EPOT KATOM Tinst"}
MD++ outpropfile = thermo.out openpropfile
MD++ run closepropfile
MD++ finalcnfile = 300K-5X5X5.cn writecn eval
MD++ sleep quit
}}

In the fuction setup_md { } in the above script, important MD++ parameters for MD simulations are assigned. First, you choose which ensemble you impose by specifying ensemble_type to be NVE, NVT, NPH or NPT. You also decide the integrator between Gear 6-th order predictor-corrector and velocity Verlet by saying integrator_type = "Gear6" or integrator_type = "VVerlet". The predictor-corrector integrator has higher order of accuracy but the Verlet integrator is symplectic and is believed to be resilient to the larger time step.

atommass is the atomic mass in g/mol, timestep is $Δ t$ used in the integrator in the unit of picosecond, totalsteps is the total number of iterations in time integration. As long as it is guaranteed that the total energy of the system is conserved, the timestep (timestep) is usually taken as big as possible. The total time step (totalsteps) is typically chosen long enough for the system to reach its dynamical steady state.

The initial velocity is ﬁrst assigned as random numbers but then scaled to correspond to an instantaneous temperature of T_OBJ.<ref>T_OBJ is the target temperature in NVT ensemble. However, it is used here just to generate atoms' initial velocities and that does not necessarily mean the temperature will be maintained at the specified value during the simulation because there is no temperature we control in NVE ensemble.</ref> The MD++ command initvelocity initializes the velocity to each atom randomly according to the value of T_OBJ. For random number generation, srand48 is executed with an integer random seed randseed. Running the simulation again with the same randseed will give you the same result (which is good for debugging). Running the simulation again with a diﬀerent randseed will give you a diﬀerent initial velocity condition (which is good for collecting statistics). You may also use the MD++ command srand48bytime, which will use the current time as the random seed (in this way, your simulation will always have a diﬀerent initial velocity condition whenever you run it again).

If the parameter DOUBLE_T is set to be 1, the initial temperature becomes twice. You can see the effect of DOUBLE_T = 1 in Fig.1. When you see the temperature plot in Fig.1, you will realize that the temperature at equilibrated state seems to be half of the initial temperature. The reason why the temperature reduces by half is followed. Since the initial atomic conﬁguration at 0 K is the local-energy-minimum structure, part of the kinetic energy initially given is always converted to the potential energy during MD simulation and this portion amounts to be half of the kinetic energy.

$\Delta E_{\mathrm{pot}} = \Delta E_{\mathrm{kin}} = \frac{3}{2} N k_B T$ .

for 3D case due to the equipartition theorem.

Fig.1 The effect of DOUBLE_T
Temperature in NVE ensemble with different DOUBLE_T setting

The MD++ command run starts MD simulation with all these parameters. There are more parameters related with MD simulation. saveprop and savepropfreq will determine whether or not the thermodynamic properties will be written in file and how often they will be. If saveprop = 1 and savepropfreq = 10, the thermodynamic properties such as potential energy, kinetic energy and temperature will be dumped out at every 10 step only if the MD++ command openpropfile is executed.

The content of a property ﬁle is defined by output_fmt. Accordig to the example script, each line of the property file contains current step (curstep), potential energy (EPOT), kinetic energy (KATOM), and temperature (Tinst). If we do not specify output_fmt, the default array of thermodynamic data output is current step (curstep), kinetic energy (KATOM), potential energy (EPOT), pressure (PRESSURE), - $σ x y$ (TSTRESS_xy), - $σ y z$ (TSTRESS_yz), - $σ z x$ (TSTRESS_zx), HELM, the extended energy(HELMP), thermodynamic friction coeﬃcient (zeta), dEdlambda and the system volume (OMEGA). Units are given as eV for energy and eV/A˚^3 for stress. The sign of the stress TSTRESS is the opposite of the conventional definition in elasticity or continuum theory, because MD++ is implemented following Parrinello and Rahman's definition of stress.<ref>M. Parrinello and A. Rahman, "Polymorphic transitions in single crystals: A new molecular dynamics method", J. Appl. Phys 52 7182-7190 1981</ref>

The thermodynamic properties can be loaded and plotted by the external program such as Octave and Matlab.<ref>A simple way to plot the total energy using gnuplot (a free software on Unix/Linux) is to run the command.

plot "prop.out" u ($1):($2+$3) with line

</ref> In Fig.2, total energy per atom and temperature are plotted for the systems of two different size. The bigger system (10[100] X 10[010] X 10[001]) has eight times more atoms than the other (5[100] X 5[010] X 5[001]). Obviously, you see that the thermodynamic data in the bigger system show less fluctuation.

Fig.2 Energy and Temperature fluctuation of NVE ensemble
(a) Energy Fluctuation in NVE ensemble.
(b) Temperature in NVE ensemble

savecn and savecnfreq are used to determine how often the intermediate configuration will be saved as .cn file with the MD++ command openintercnfile just like saveprop and savepropfreq. You can specify names of the property file and the intermediate configuration files using MD++ parameter outpropfile and intercnfile, respectively. If we do not specify them, the default ﬁle name would be prop.out and inter.cn. The expression zipfiles = 1 in the 11th line will compress ﬁles such as .out and .cn. For example, the conﬁguration ﬁle init.cn and the thermodynamic data ﬁle thermo.out will be zipped as init.cn.gz and thermo.out.gz, respectively.

Calling writecn saves a conﬁguration (.cn) ﬁle, which contains the following informations; The ﬁrst line speciﬁes the total number of atoms (NP). After that there are NP lines of data and each line has three columns in default, which correspond to the scaled coordinate ( $s x$ , $s y$ , $s z$ ) of one atom. And then there are three lines of data that speciﬁes a 3 X 3 box matrix $\mathbf{H}$ . The real coordinates of an atom (x, y, z ) and the scaled coordinates are related by the following

$\left( \begin{array}{c} x \\ y \\ z \end{array} \right) = \left( \begin{array}{ccc} H_{11} & H_{12} & H_{13} \\ H_{21} & H_{22} & H_{23} \\ H_{31} & H_{32} & H_{33} \end{array} \right) \left( \begin{array}{c} s_x \\ s_y \\ s_z \end{array} \right)$ .

If writeall = 1, additional informations are added to NP lines of data such that each line has sacaled velocities (_VSR[i].x, _VSR[i].y and _VSR[i].z), individual potential energy (_EPOT_IND[i]), whether-or-not ﬁxed (fixed[i]), topology (_TOPOL[i]), atom species (species[i]), atom groups (group[i]) and images (image[i]) after the scaled coordinates(_SR[i].x, _SR[i].y and _SR[i].z for i=1 to NP) . You can check this diﬀerence by comparing two conﬁguration ﬁles, perf.cn and init.cn. Because perf.cn is written before setting writeall = 1, it has only number of atoms, scaled coordinates of each atom and the box matrix $\mathbf{H}$ while init.cn has all the informations. Alternatively, writevelocity = 1 can be used, which stores only the scaled coordinates and the scale velocities.

Commands and variables related with visualization in openwindow will be covered in a separate manual in detail.

After the simulation (run) has ﬁnished, the sleep command allows the graphics window to stay alive for a while so that we can look at the ﬁnal atomic structure. You can press ctrl-c in the terminal to exit. You can also comment out the sleep command so that MD++ quits right after ﬁnishing run.

Canonical(NVT) Ensemble

The following example script mo_NVT.tcl reads the structure file of BCC Mo and runs NVT MD simulation during which temperature is maintained at 300 K using Nose-Hoover thermostat<ref> S. Nose, “A Molecular Dynamics Method for Simulations in the Canonical Ensemble”, Molecular Physics, 52 255 (1984)</ref>. Run the script by typing

$ bin/fs_gpp scripts/mo_NVT.tcl

# -*-shell-script -*-
# run NVT MD simulation of perfect crystal of Mo.
#
#*******************************************
# Definition of procedures
#*******************************************
proc initmd { } { 
         ...
}
#------------------------------------------------------------
proc readMoPot { } { 
         ...
}
#-------------------------------------------------------------
proc openwindow { } { 
         ...
}
#--------------------------------------------
proc setup_md { } { MD++ {
#MD settings
ensemble_type = "NVT" integrator_type = "VVerlet"
implementation_type = 0 vt2 = 1e28
T_OBJ = 300       # (in Kelvin) Desired Temperature
atommass = 95.94  # (g/mol)
timestep = 0.001 # (ps)
totalsteps = 2000
saveprop = 1 savepropfreq = 10
savecn = 0 savecnfreq = 100
writeall = 1
}}
#*******************************************
# Main program starts here
#*******************************************
initmd
readMoPot
MD++ {
  incnfile = 300K-5X5X5.cn
  atommass = 95.94 timestep = 0.001
  readcn eval
}
openwindow
#---------------------------------------------
# run MD
setup_md
MD++ {output_fmt = "curstep EPOT KATOM Tinst HELMP"}
MD++ outpropfile = thermo-NVT.out openpropfile
MD++ run eval 
#reverse velocity (to test reversibility)
MD++ input = -1 multiplyvelocity
#MD simulation (reverse)
MD++ run eval
MD++ finalcnfile = 300K-5X5X5-NVT.cn writecn
MD++ sleep quit

Instead of creating BCC Mo crystal, this script reads the conﬁguration ﬁle 300K_5X5X5.cn, which was generated by the previous script mo_NVE.tcl, using the MD++ command readcn. Note that when MD++ reads the configuration ﬁle, the timestep should be idnetical with that of the previous run because the scaled velocity is stored as dimensionless unit after being multiplied with the timestep. It is safe that you specify timestep and atommass when reading .cn ﬁle as shown in the script.

Setting ensemble_type = "NVT" activates Nose-Hoover thermostat, which basically adds an additional variable ( $ζ$ ) to the system such that the temperature is controlled by the feedback routine virtually mimicking heat transfer to the heat reservoir,

$\begin{array}{rcl} \ddot\mathbf{r}_i &=& -\frac{1}{m_i}\frac{\partial U(\{\mathbf{r}_i\})}{\partial \mathbf{r}_i} - \zeta \dot\mathbf{r}_i \\ \dot\zeta &=& \frac{3N k_B}{Q} ( T_{\mathrm{inst}} - T_{\mathrm{OBJ}} ) \end{array}$ .

when the degree of freedom is 3N. In MD++, Nose-Hoover thermostat are implemented in three diﬀerent way, each of which uses different integrator deﬁned through implementation_type. If it is zero, the implicit integrator will be chosen. If it is one, the explicit integrator based on Sto¨rmer-Verlet method will work. If it is two, another explicit integrator based on Liouville formulation will be activated.

Thermal coeﬃcient vt2 is related with thermal mass Q of Nose´-Hoover thermostat as

$Q = \frac{3 N k_B T_{\mathrm{OBJ}}}{\mathrm{vt2}}\times (1e+24)$ .

where N is number of atoms and k_B is the Boltzmann constant (8.617343e-5 eV/K). For convenience, we also defined thermal mass NHMass, which is exactly Nose'-Hoover thermal mass Q. When N = 250 and T = 300 K (or k_B T ~ 0.026 eV), vt2 = 1e28 corresponds to NHMass = 1.95e-03. The bigger vt2 is (or the smaller NHMass is), the faster the response of the system is to control the temperature. The eﬀect of vt2 on Helmholtz free energy and temperature ﬂuctuation is shown in Fig. 3.

Fig.3 The Effect of Nose-Hoover Mass in NVT MD
(a) Helmholtz energy fluctuation in NVT ensemble.
(b) Temperature in NVT ensemble

HELMP in the output_fmt is Helmholtz free energy which is the conserved quantity in NVT MD.

multiplyvelocity multiplies the number speciﬁed by the preceding input to all velocity components of all atoms. With input=-1, the MD++ command multiplyvelocity reverses the velocity direction of all the atoms and they are expected to be at the initial position after the reverse run. This is done to test the reversibility of the symplectic integrator. Comparing two configuration files 300K-5X5X5.cn and 300K-5X5X5-NVT.cn reveals that the numbers become different after 10^(-14) digit. If we consider the total step is 4000, the difference at each step would be O(10^(-17)), which corresponds to the machine precision. On the other hand, if you test integrator_type = Gear6, the numbers become different only after 10^(-6) digit. We can say the symplectic integrator is indeed reversible.

Canonical(NVT) Ensemble with Nose-Hoover Chain

Martyna et al<ref>G. J. Martyna, M. L. Klein and M. Tuckerman, "Nose-Hoover chains: The canonical ensemble via continuous dynamics", J. Chem. Phys. 97 2635--2643</ref> modified the original Nose-Hoover thermostat such that a chain of thermostats is added instead of one thermostat. The purpose of this chain method is to increase the size of the phase space and help the system to be ergodic even when the system is stiff or its dynamics is not ergodic. To implement the Nose-Hoover chain integrator, we refered another paper of Martyna's.<ref>G. J. Martyna, "Explicit reversible integrators for extended systems dynamics", Molecular Physics 87 1117-1157 (1996)</ref>

The MD++ script Mo_NVTC.tcl, given below, shows how to run Nose-Hoover chain method in MD++ to obtain a system of canonical ensemble.

# -*-shell-script -*-
# run NVE MD simulation of perfect crystal of Mo.
#
#*******************************************
# Definition of procedures
#*******************************************
proc initmd { } { MD++ {
    ....
}}
#------------------------------------------------------------
proc readMoPot { } { MD++ {
    ....
}}
#-------------------------------------------------------------
proc openwindow { } { MD++ {
    ....
}}
#--------------------------------------------
proc setup_md { } { MD++ {
#MD settings
ensemble_type = "NVTC" integrator_type = "VVerlet"
NHChainLen = 4                       # MAXNHCLEN = 20 in md.h
NHMass = [ 1.95e-4 2e-6 2e-6 2e-6 ]
T_OBJ = 300        # (in Kelvin) Desired Temperature
atommass = 95.94   # (in g/mol)
timestep = 0.001   # (in ps)
totalsteps = 5000
saveprop = 1 savepropfreq = 10
savecn = 0 savecnfreq = 100
writeall = 1
DOUBLE_T = 0 randseed = 12345 srand48
}}

#*******************************************
# Main program starts here
#*******************************************
initmd
readMoPot
MD++ {
  incnfile = ../mo-example/300K-5X5X5.cn
  atommass = 95.94 timestep = 0.001
  readcn eval
}
openwindow
#---------------------------------------------
# run MD
setup_md
MD++ {output_fmt = "curstep EPOT KATOM Tinst HELMP"}
MD++ outpropfile = thermo-NVTC.out openpropfile
MD++ run eval
MD++ sleep quit
}}

You choose the ensemble type to be NVTC. The number of chains and masses of each chain can be specified using NHChainLen and NHMass. According to my experience, NHMass needs to be chosen with caution, otherwise the system may lose its reversibility, even though Martyna et al stated that the choice of thermal mass NHMass is not critical. In fact, they also presented how to take reasonable choice of NHMass.

NHMass[0] ~  $N f k B T / ω 2$  where  is the number of degrees of freedom.
And NHMass[i] ~  $k B T / ω 2$  (for i = 2 to M, where M is the number of chains)

thermostat 1 to M -1 fluctuates at a frequency of omega.

Constant Stress (N $σ$ H) Ensemble

In constant stress ensemble, the stress is controlled to maintain the external stress (or the applied stress) using Parrinello-Rahman's method. The applied stress tensor $σ EXT$ can be decomposed as the sum of the hydrostatic term and the deviatoric term,

$\sigma_{\mathrm{EXT}} = \sigma_{\mathrm{EXT}}^{\mathrm{hydro}} + \sigma_{\mathrm{EXT}}^{\mathrm{devi}}$ .

where $\sigma_{\mathrm{EXT}}^{\mathrm{hydro}} = \frac{1}{3}(\sigma_{\mathrm{EXT}})_{kk}\cdot\mathrm{I}$ and $\sigma_{\mathrm{EXT}}^{\mathrm{devi}} = \sigma_{\mathrm{EXT}} - \sigma_{\mathrm{EXT}}^{\mathrm{hydro}}$ . Here we used the Einstein's convention, in which repeated indices are summed over from 1 to 3, to express the sum of the diagonal components of the stress tensor. Then the equations of motion of Parrinello-Rahman's method can be written as

$\begin{array}{rcl} \ddot\mathbf{r}_i &=& -\frac{1}{m_i}\frac{\partial U(\{\mathbf{r}_i\})}{\partial \mathbf{r}_i} - \mathbf{H}^{-T}\dot\mathbf{G}\mathbf{H}^{-1} \dot\mathbf{r}_i \\ \ddot\mathbf{H} &=& -\frac{1}{M} \left[ \left( \sigma_{\mathrm{Virial}} - \sigma_{\mathrm{EXT}}^{\mathrm{hydro}} \right)\Xi - \left( \mathbf{H}\mathbf{H}_0^{-1} \sigma_{\mathrm{EXT}}^{\mathrm{devi}} \right) \Xi_0 \right] \end{array}$ .

where $\mathbf{H}$ is the 3X3 box matrix at time t and $\mathbf{H}_0$ is the initial box matrix at t = 0. $\mathbf{G}$ is defined as $\mathbf{G} = \mathbf{H}^T\mathbf{H}$ and called the metric tensor. $\Xi = \Omega \mathbf{H}^{-T}$ and $\Xi_0 = \Omega_0 \mathbf{H}_0^{-T}$ are used to reduce the notation. ( $Ω$ is the volume of the box or $\Omega = \det{\mathbf{H}}$ .) M is the "box mass".

If we apply the constant pressure $P_{\mathrm{EXT}} = -\frac{1}{3}\left(\sigma_{\mathrm{EXT}}\right)_{kk}$ , the deviatoric part of the external stress is zero and the second equation becomes simply

$\ddot\mathbf{H} = -\frac{1}{M} \left[ \left( \sigma_{\mathrm{Virial}} + P_{\mathrm{EXT}}\mathbf{I} \right)\Xi \right]$ .

and this follows isoenthalpic-isobaric (NPH) Ensemble for large number of particles. The script Mo_NPH.tcl, given below, shows how we set parameters and commands to run NPH MD simulation using MD++.

# -*-shell-script -*-
# run NPH MD simulation of perfect crystal of Mo.
#
#*******************************************
# Definition of procedures
#*******************************************
proc initmd { } { MD++ {
        ...
}}
#------------------------------------------------------------
proc readMoPot { } { MD++ {
        ...
}}
#-------------------------------------------------------------
proc openwindow { } { MD++ {
        ...
}}
#--------------------------------------------
proc setup_md { } { MD++ {
#MD settings
ensemble_type = "NPH" integrator_type = "VVerlet"
wallmass = 1000 boxdamp = 0.1
saveH # Use current H as reference (H0), needed for specifying stress
stress = [ 1000    0    0
              0 1000    0
              0    0 1000 ] # compression in MPa
conj_fixboxvec = [ 0 1 1
                   1 0 1
                   1 1 0 ]
T_OBJ = 300                 # (in Kelvin) Desired Temperature
atommass = 95.94            # (in g/mol)
timestep = 0.001            # (in ps)
totalsteps = 10000
saveprop = 1 savepropfreq = 20
savecn = 0 savecnfreq = 100
writeall = 1
DOUBLE_T = 0 randseed = 12345 srand48
}}

#*******************************************
# Main program starts here
#*******************************************
initmd
readMoPot
MD++ makecrystal finalcnfile = perf.cn writecn
openwindow
#---------------------------------------------
# run MD
setup_md
MD++ initvelocity finalcnfile = init.cn writecn
MD++ {output_fmt = "curstep EPOT KATOM Tinst HELMP \
                    TSTRESSinMPa_xx TSTRESSinMPa_yy TSTRESSinMPa_zz \
                    TSTRESSinMPa_xy TSTRESSinMPa_xz TSTRESSinMPa_yz \
                    H_11 H_22 H_33 H_12 H_13 H_21 H_23 H_31 H_32" }
MD++ outpropfile = thermo-NPH.out openpropfile
MD++ run closepropfile
#MD++ finalcnfile = Mo-5X5X5-NPH.cn writecn eval
MD++ sleep quit
}}

First, you choose the ensemble type and the integrator type. wallmass corresponds to the box mass M of the Parrinello-Rahman (PR) method. boxdamp is used to make the system converge fast. You can see the effect of boxdamp in Fig.4. If boxdamp is nonzero, the box acceleration is redefined as

$\ddot\mathbf{H} := \ddot\mathbf{H}- (\mathrm{boxdamp}) \times \dot\mathbf{H}$

and the system becomes no longer symplectic. The external stress is assigned throuth MD++ variable stress in MPa. Note that the sign convention of stress is opposite of that from the elasticity theory. The positive stress in the MD++ script means that you apply compression to the system. conj_fixboxvec defines which component of the box matrix $H$ will be fixed. you need to specify $H 0$ by calling saveH before using the PR method. According to Ray and Rahman<ref>J. R. Ray and A. Rahman, "Statistical ensembles and molecular dynamics studies of anisotropic solids", J. Chem. Phys, 80 4423--4428 </ref>, " $H 0$ should be taken as the average value of $H$ when the stress is zero."

Fig.4 The Effect of Box Damp in NPH MD
(a) Extended energy fluctuation in NPH ensemble.
(b) Stress in NPH ensemble.
(c) Simulation box in NPH ensemble

Isothermal-isobaric (NPT) Ensemble

Manual 05 for MD++
Molecular Dynamics Simulations

Contents

Microcanonical (NVE) Ensemble

Canonical(NVT) Ensemble

Canonical(NVT) Ensemble with Nose-Hoover Chain

Constant Stress (N $σ$ H) Ensemble

Isothermal-isobaric (NPT) Ensemble

Manual 05 for MD++ Molecular Dynamics Simulations

Contents

Microcanonical (NVE) Ensemble

Canonical(NVT) Ensemble

Canonical(NVT) Ensemble with Nose-Hoover Chain

Constant Stress (NσH) Ensemble

Isothermal-isobaric (NPT) Ensemble

Manual 05 for MD++
Molecular Dynamics Simulations

Constant Stress (N $σ$ H) Ensemble