Tutorial for CPU & GPU- The unified “middle” thermostat scheme in AMBER for efficient configurational sampling

Introduction

The tutorial offers an introduction to the unified “middle” scheme of thermostatting algorithms in AMBER.

The prerequisites are listed as below:

• AMBER software package (2018 or 2019 version) - sander and its parallel version, sander.MPI, are used to run simulations. In addition, LEaP (tleap and its graphical user interface, xleap) is utilized to generate the topology file of the simulated system.

The “middle” scheme offers a unified framework to develop efficient thermostatting algorithms for configurational sampling for the canonical ensemble, as described in Refs. [1–5] . It can be implemented for performing molecular dynamics (MD) or path integral molecular dynamics (PIMD), either with or without holonomic constraints. The “middle” scheme allows the use of much larger time intervals (i.e., time stepsizes) \(\Delta t\) to maintain the same accuracy, which significantly improves the configurational sampling efficiency. That is, it is efficient for calculating structural properties and thermodynamic observables that depend on coordinate variables. Most thermostats control the temperature by updating momenta of the system. Some prevailing thermostats include stochastic ones (such as the Andersen thermostat and Langevin dynamics) and deterministic ones (such as the Nosé-Hoover thermostat and Nosé-Hoover chain). In the “middle” scheme, immediately after the coordinate-updating step for half a time interval, the thermostat process for a full time interval takes place, which then followed by the coordinate-updating step for another half time interval [1,5] .

Here we present a brief introduction to the “middle” scheme. For many thermostats, the integration in one time step \(\Delta t\) can be splitted into three parts, the steps for updating coordinates, momenta and thermostat, denoted as “x”, “p” and “T”, respectively. In this case the “equations of motion” may be expressed as

\[ \begin{bmatrix} \mathrm{d}\mathbf{x}_t\\ \mathrm{d}\mathbf{p}_t \end{bmatrix} = \underbrace{ \begin{bmatrix} \mathbf{M}^{-1}\mathbf{p}_t\\ 0 \end{bmatrix}\mathrm{d}t }_{\text{x}} + \underbrace{ \begin{bmatrix} 0 \\ -\nabla_{\mathbf{x}_t}U(\mathbf{x}_t) \end{bmatrix}\mathrm{d}t }_{\text{p}} + \underbrace{ \begin{bmatrix} \text{thermostat} \\ \end{bmatrix} }_{\text{T}} \] (1) Here, \(U\) is the potential energy, \(\mathbf{M}\) is the diagonal mass matrix, \(\mathbf{x}\) and \(\mathbf{p}\) are the vectors of coordinate and momentum, respectively. Equation (1) is, however, not convenient to do the analysis.

A more useful approach is to employ the forward Kolmogorov equation to express the evolution of the density distribution in the phase space \(\rho(\mathbf{x},\mathbf{p})\) . \[ \begin{aligned} \frac{\partial}{\partial t}\rho =& \mathcal{L}\rho\\ =& (\mathcal{L}_{\text{x}} + \mathcal{L}_{\text{p}} + \mathcal{L}_{\text{T}})\rho \end{aligned} \] (2) The relevant Kolmogorov operators for the 1st and 2nd terms of the right-hand side (RHS) are \[ \mathcal{L}_{\text{x}}\rho = -\mathbf{p}^T\mathbf{M}^{-1}\nabla_{\mathbf{x}}\rho \] (3) \[ \mathcal{L}_{\text{p}}\rho = \nabla{}_{\mathbf{x}}U(\mathbf{x})\cdot\nabla{}_\mathbf{p}\rho \] (4) respectively. The definition of \(\mathcal{L}_{\text{T}}\) depends on the specific thermostat. The phase space propagators for a time interval \(\Delta t\) for the three parts are \(e^{\mathcal{L}_{\text{x}}\Delta t}\) , \(e^{\mathcal{L}_{\text{p}}\Delta t}\) , and \(e^{\mathcal{L}_{\text{T}}\Delta t}\) , respectively.

The propagation in each time step with the velocity Verlet (VV) algorithm is performed as \[ e^{\mathcal{L}\Delta t} \approx e^{\mathcal{L}^{\text{VV}}_{\text{middle}}\Delta t} = e^{\mathcal{L}_{\text{p}}\Delta t/2} e^{\mathcal{L}_{\text{x}}\Delta t/2} e^{\mathcal{L}_{\text{T}}\Delta t} e^{\mathcal{L}_{\text{x}}\Delta t/2} e^{\mathcal{L}_{\text{p}}\Delta t/2} \] (5) The phase space propagator for the thermostat part \(e^{\mathcal{L}_{\text{T}}\Delta t}\) is designed in the middle. Equation (5) is denoted as “VVMiddle”. The numerical algorithm reads \[ \begin{aligned} \text{Update Momenta for half a step:} \quad &\mathbf{p} \leftarrow \mathbf{p} - \frac{\partial U}{\partial \mathbf{x}} \frac{\Delta t}{2}\\ \text{Update Coordinates for half a step:} \quad &\mathbf{x} \leftarrow \mathbf{x} + \mathbf{M}^{-1}\mathbf{p} \frac{\Delta t}{2}\\ \text{Thermostat for a full time step:} \quad &\textit{thermostat_step} \\ \text{Update Coordinates for another half step:} \quad &\mathbf{x} \leftarrow \mathbf{x} + \mathbf{M}^{-1}\mathbf{p} \frac{\Delta t}{2}\\ \text{Update Momenta for another half step:} \quad &\mathbf{p} \leftarrow \mathbf{p} - \frac{\partial U}{\partial \mathbf{x}} \frac{\Delta t}{2}\\ \end{aligned} \] (6) where is the subroutine for the thermostat process, which is determined according to the thermostat method of choice.

The stationary state distribution of “VVMiddle” for a harmonic system \(U(\mathbf{x})=\frac{1}{2}(\mathbf{x}-\mathbf{x}_{\text{eq}})^T \mathbf{A} (\mathbf{x}-\mathbf{x}_{\text{eq}})\) is \[ \begin{aligned} \rho^{\text{VV}}_{\text{middle}}(\mathbf{x},\mathbf{p})=&\frac{1}{Z_N}\exp \left\{-\beta\left[ \frac{1}{2}\mathbf{p}^T\left(\mathbf{M}-\mathbf{A}\frac{\Delta t^2}{4}\right)^{-1}\mathbf{p} \right. \right.\\ &\left. \left. +\frac{1}{2}(\mathbf{x}-\mathbf{x}_{\text{eq}})^T \mathbf{A} (\mathbf{x}-\mathbf{x}_{\text{eq}}) \right]\right\} \end{aligned} \] (7) as long as the thermostat process keeps the Maxwell (or Maxwell-Boltzmann) momentum distribution unchanged, i.e. \[ e^{\mathcal{L}_{\text{T}}\Delta t}\rho_{\text{MB}}(\mathbf{p}) = \rho_{\text{MB}}(\mathbf{p}) \] (8) where the Maxwell momentum distribution is \[ \rho_{\text{MB}}(\mathbf{p}) = \left(\frac{\beta}{2\pi}\right)^{3N/2}\lvert\mathbf{M}\rvert^{-1/2} \exp\left[-\beta\left(\frac{1}{2}\mathbf{p}^T\mathbf{M}^{-1}\mathbf{p}\right)\right] \] (9) Here \(\beta=\frac{1}{k_BT}\) with \(k_B\) as the Boltzmann constant, \(T\) is the temperature of the system. (\(N\) is the number of particles.) It is then easy to verify that the marginal distribution of coordinates for “VVMiddle” is exact in the harmonic limit. Many types of thermostats satisfy the criteria (thermostat process keeps the Maxwell momentum distribution unchanged, Equation (8), which include, but not limited to, the thermostats listed below.

• Andersen thermostat (real dynamics case)

In this thermostat, each particle of the system stochastically collides with a fictitious heat bath, and once the collision occurs, the momentum of this particle is chosen afresh from the Maxwell-Boltzmann momentum distribution. The explicit form for the thermostat process can be expressed as \[ \begin{aligned} \mathbf{p}^{(j)}&\leftarrow \sqrt{\frac{m_j}{\beta}}\boldsymbol{\theta}_j, \ (j=\overline{1,N})\\ &\text{if } \mu_j < \nu \Delta t \ ( \text{or more precisely } \mu_j < 1-e^{-\nu \Delta t}) \end{aligned} \] (10) Here \(\nu\) is the collision frequency, \(\boldsymbol{\theta}_j\) is a vector of independent Gaussian-distributed random numbers with zero mean and unit variance, \(m_j\) the mass for the \(j\) th atom, \(\mu_j\) is a uniformly distributed random number in the range (0,1). Here \(\mu_j\) may be different for each particle (\(j=\overline{1,N}\) ). In the current version of AMBER \(\mu_j\) is the same for all particles.

The phase space propagator for the thermostat process is \[ \begin{aligned} e^{\mathcal{L}_{\text{T}}\Delta t}\rho = &e^{-\nu \Delta t}\rho(\mathbf{x},\mathbf{p})\\ &+ (1-e^{-\nu \Delta t})\rho_{\text{MB}}(\mathbf{p}) \int_{-\infty}^{\infty}\rho(\mathbf{x},\mathbf{p})\mathrm{d}\mathbf{p} \end{aligned} \] (11)

• Andersen thermostat (virtual dynamics case)

The explicit form for the virtual dynamics case of the Andersen thermostat is expressed as \[ \left. \begin{aligned} \mathbf{p}^{(j)}\leftarrow & \sqrt{\frac{m_j}{\beta}}\boldsymbol{\theta}_j, \text{ if } \mu_j < 1-e^{-\nu \Delta t}\\ \mathbf{p}^{(j)}\leftarrow &- \mathbf{p}^{(j)}, \quad \text{otherwise} \end{aligned} \right\rbrace (j=\overline{1,N}) \] (12) The phase space propagator for the thermostat process is \[ \begin{aligned} e^{\mathcal{L}_{\text{T}}\Delta t}\rho = &e^{-\nu \Delta t}\rho(\mathbf{x},-\mathbf{p})\\ &+ (1-e^{-\nu \Delta t})\rho_{\text{MB}}(\mathbf{p}) \int_{-\infty}^{\infty}\rho(\mathbf{x},\mathbf{p})\mathrm{d}\mathbf{p} \end{aligned} \] (13)

• Langevin dynamics (real dynamics case)

The thermostat process is the solution to the Ornstein-Uhlenbeck (OU) process \[ \mathbf{p}\leftarrow e^{-\boldsymbol{\gamma}\Delta t}\mathbf{p} + \sqrt{\frac{1}{\beta}}\mathbf{M}^{1/2}(\mathbf{1}-e^{-2\boldsymbol{\gamma}\Delta t})^{1/2}\boldsymbol{\eta} \] (14) Here, \(\boldsymbol{\gamma}\) is the diagonal friction coefficient matrix. In the current version of AMBER all diagonal elements of \(\boldsymbol{\gamma}\) are set to be the same. (That is, the friction coefficient is the same for all particles.)

The phase space propagator for the thermostat process is \[ \begin{aligned} e^{\mathcal{L}_{\text{T}}\Delta t}\rho =&\left(\frac{\beta}{2\pi}\right)^{3N/2}\lvert \mathbf{M} (\mathbf{1}-e^{-2\boldsymbol{\gamma}\Delta t})\rvert^{-1/2} \\ &\cdot\int\mathrm{d}\mathbf{p}_0\,\rho(\mathbf{x},\mathbf{p}_0) \exp\left[-\frac{\beta}{2}(\mathbf{p}-e^{-\boldsymbol{\gamma}\Delta t}\mathbf{p}_0)^T\right. \\&\left. \cdot\mathbf{M}^{-1}(\mathbf{1}-e^{-2\boldsymbol{\gamma}\Delta t})^{-1}(\mathbf{p}-e^{-\boldsymbol{\gamma}\Delta t}\mathbf{p}_0)\right] \end{aligned} \] (15)

• Langevin dynamics (virtual dynamics case)

The virtual dynamics case represents another type of discrete evolution that may not correspond to a continuous, real dynamical counterpart of the Langevin equation. \[ \mathbf{p}\leftarrow -e^{-\boldsymbol{\gamma}\Delta t}\mathbf{p} + \sqrt{\frac{1}{\beta}}\mathbf{M}^{1/2}(\mathbf{1}-e^{-2\boldsymbol{\gamma}\Delta t})^{1/2}\boldsymbol{\eta} \] (16) The virtual dynamics case is also able to produce the desired stationary distribution.
The phase space propagator for the thermostat process is then \[ \begin{aligned} e^{\mathcal{L}_{\text{T}}\Delta t}\rho =&\left(\frac{\beta}{2\pi}\right)^{3N/2}\lvert \mathbf{M} (\mathbf{1}-e^{-2\boldsymbol{\gamma}\Delta t})\rvert^{-1/2} \\ &\cdot\int\mathrm{d}\mathbf{p}_0\,\rho(\mathbf{x},\mathbf{p}_0) \exp\left[-\frac{\beta}{2}(\mathbf{p}+e^{-\boldsymbol{\gamma}\Delta t}\mathbf{p}_0)^T\right. \\&\left. \cdot\mathbf{M}^{-1}(\mathbf{1}-e^{-2\boldsymbol{\gamma}\Delta t})^{-1}(\mathbf{p}+e^{-\boldsymbol{\gamma}\Delta t}\mathbf{p}_0)\right] \end{aligned} \] (17)

• Nosé-Hoover (NH) thermostat and Nosé-Hoover chain (NHC)

See Ref. [1] for more detailed discussions.

The “middle” scheme of a thermostat includes both real and virtual dynamics cases. (See Refs. [3,4] .) It is proved in Ref. [3] that, while the Langevin equation algorithm (BAOAB) proposed in Ref. [6] is simply only the real dynamics case of “VVMiddle”, another Langevin equation algorithm proposed (without employing the Lie-Trotter splitting) in Ref. [7] is equivalent to “VVMiddle” for Langevin dynamics. The real dynamics case for the Andersen thermostat and that for Langevin dynamics have been implemented in the current version of AMBER.

When the leapfrog algorithm, rather than the velocity-Verlet algorithm, is employed in the “middle” scheme, it is denoted as “LFMiddle” [5] . The propagation in each time step with the leapfrog (LF) algorithm is performed as \[ e^{\mathcal{L}\Delta t} \approx e^{\mathcal{L}^{\text{LF}}_{\text{middle}}\Delta t} = e^{\mathcal{L}_{\text{x}}\Delta t/2} e^{\mathcal{L}_{\text{T}}\Delta t} e^{\mathcal{L}_{\text{x}}\Delta t/2} e^{\mathcal{L}_{\text{p}}\Delta t} \] (18) The numerical algorithm of “LFMiddle” reads \[ \begin{aligned} \text{Update Momenta for a full time step:} \quad &\mathbf{p} \leftarrow \mathbf{p} - \frac{\partial U}{\partial \mathbf{x}} {\Delta t}\\ \text{Update Coordinates for half a step:} \quad &\mathbf{x} \leftarrow \mathbf{x} + \mathbf{M}^{-1}\mathbf{p} \frac{\Delta t}{2}\\ \text{Thermostat for a full time step:} \quad &\textit{thermostat_step} \\ \text{Update Coordinates for another half step:} \quad &\mathbf{x} \leftarrow \mathbf{x} + \mathbf{M}^{-1}\mathbf{p} \frac{\Delta t}{2} \end{aligned} \] (19) It has the same accuracy as “VVMiddle” for sampling the marginal distribution in the coordinate space. The marginal distribution of momentum is also accurate, which has an interval for half a time step with the coordinates. For simplicity and compatibility, only “LFMiddle” is integrated into AMBER.

The “middle” scheme with bond length constraints is also implemented. While MD with holonomic constraints are widely used in biological simulations, PIMD with holonomic constraints may help understand nuclear quantum effects of different motions in molecular systems. For instance, help assign spectral features as shown in Ref. [8] . In the “middle” scheme, when holonomic constraints are applied, the SHAKE [9] and RATTLE [10] algorithms are used for fixing coordinates and velocities, respectively. Particularly for the molecular system that contains water molecules, the analytical SETTLE algorithm [11] may be used to apply the constraint to the water molecule.

Besides Refs. [1–5] , one more paper is in preparation for the “middle” scheme [12] .

While the “middle” scheme for PIMD with the staging transformation (staging PIMD) was first demonstrated in Ref. [2] , that for PIMD with the normal-mode transformation (normal-mode PIMD) was first proposed in the supplemental material of Ref. [2] in 2016, which can be found via the URLs provided by the publisher:

• https://aip.scitation.org/doi/suppl/10.1063/1.4954990/suppl_file/supplemental+material-submitted.docx
• ftp://ftp.aip.org/epaps/journ_chem_phys/E-JCPSA6-145-007626

In addition, the arXiv preprint (that includes Ref. [2] and its supplemental material) is also available (https://arxiv.org/ftp/arxiv/papers/1611/1611.06331.pdf).
In the current version, the “middle” scheme is implemented for the primitive version of PIMD (PRIMPIMD) of AMBER. The staging PIMD or normal-mode PIMD algorithms in the “middle” scheme will also be implemented into AMBER soon.

Input parameters

In order to perform MD or PRIMPIMD simulations with the “middle” scheme, three additional flags should be added in the file, which are used to distinguish different methods.

ischeme	Flag for choosing an integration scheme for molecular dynamics.
	=0 (default) conventional scheme in AMBER.
	=1 “middle” scheme based on the leapfrog algorithm.
ithermostat	Flag for different thermostats when the “middle” scheme is employed. Two types of thermostats are currently available.
	=1 Langevin dynamics
	=2 Andersen thermostat
therm_par	The parameter used in a thermostatting method of the “middle” scheme, in the unit of ps\(^{-1}\) , which should always be a positive number. It refers to the friction coefficient for Langevin dynamics ( = 1) or the collision frequency for the Andersen thermostat ( = 2).

The recommended value for is related to the characteristic frequency (\(\tilde{\omega}\) ) of the specific system. The characteristic time of the potential energy autocorrelation function is \[ \tau_{UU} = \int_0^{\infty}\frac{\langle U(0)U(t)\rangle - \langle U \rangle^2}{\langle U^2\rangle - \langle U \rangle^2}\,\mathrm{d}t \] (20) The optimal value of the thermostat parameter that produces the minimum correlation time of the potential is \(\xi^{opt} \approx \tilde{\omega}\) for Langevin dynamics and \(\xi^{opt}\approx\sqrt{2}\tilde{\omega}\) for the Andersen thermostat, as the time interval \(\Delta t\) approaches zero. E.g. for a HO molecule, the frequency of the O-H stretch is around 3600 cm\(^{-1}\) (~680 ps\(^{-1}\) ), so one can choose 680 ps\(^{-1}\) as the value of when Langevin dynamics is used, or 960 ps\(^{-1}\) when the Andersen thermostat is employed. When the time interval \(\Delta t\) is finite in the two thermostatting methods, while the characteristic correlation time goes to infinity as the thermostat parameter approaches zero, the characteristic correlation time gradually reaches a plateau as the thermostat parameter increases. (Please see Refs. [3,4] for more discussions.) When condensed phase systems are simulated, it is not straightforward to estimate the optimal thermostat parameter(s) that could be related to the mixing of frequencies or time scales of the system [13] . Some numerical tests are necessary for obtaining the reasonable region for the thermostat parameter. For a liquid water system (216 water molecules in a cell with periodic boundary conditions) with no holonomic constraints, the thermostat parameter is usually chosen to be \(2-50\) ps\(^{-1}\) .

Test cases

Molecular dynamics (for classical statistics)

(a) MD input using the Langevin thermostat with the “LFMiddle” scheme for liquid water.

Test: $AMBERHOME/test/middle-scheme/MD_Unconstr_Langevin_water

MD: NVT simulation of liquid water
&cntrl    
ipimd = 0, nstlim = 10        ! MD for 10 steps    
ntx = 1, irest = 0            ! read coordinates    
temp0 = 300, tempi = 300      ! temperature: target and initial    
dt = 0.001                    ! time step in ps    
cut = 7.0                     ! non-bond cut off    
ischeme = 1                   !! leapfrog middle scheme
ithermostat = 1               !! Langevin thermostat   
therm_par = 5.0               !! thermostat parameter in 1/ps
ig = 1000                     ! random seed 
ntc = 1, ntf = 1              ! no constraints
ntpr = 1, ntwr = 5, ntwx = 5  ! output settings 
/ 

One can run either a serial job (using sander ):

$ sander -O -i md_LGV.in -p qspcfw216.top -c nvt.rst -o md_LGV.out \
  -r lgv.rst -info lgv.info

or a parallel job (using sander.MPI ):

$ mpirun -np 4 sander.MPI -O -i md_LGV.in -p qspcfw216.top -c nvt.rst \
  -o md_LGV.out -r lgv.rst -info lgv.info

(b) MD input using Langevin dynamics with the “LFMiddle” scheme for the liquid water. Lengths of the bonds having hydrogen atoms are constrained.

Test: $AMBERHOME/test/middle-scheme/MD_Constr_Langevin_water

MD: NVT simulation of liquid water
&cntrl
ipimd = 0, nstlim = 10        ! MD for 10 steps
ntx = 1, irest = 0            ! read coordinates
temp0 = 300, tempi = 300      ! temperature: target and initial
dt = 0.004                    ! time step in ps
cut = 7.0                     ! non-bond cut off
ischeme = 1,                  !! leapfrog middle scheme
ithermostat = 1,              !! Langevin thermostat, random seed is default value
therm_par = 5.0               !! thermostat parameter, in 1/ps
ntc = 2, ntf = 2              ! constrain lengths of the bonds having hydrogen atoms
ntpr = 1, ntwr = 5, ntwx = 5  ! output settings 
/

Run either a serial way (using sander ):

$ sander -O -i md_LGV.in -p qspcfw216.top -c nvt.rst -o md_LGV.out \
  -r lgv.rst -info lgv.info

or a parallel job (using sander.MPI ):

$ mpirun -np 4 sander.MPI -O -i md_LGV.in -p qspcfw216.top -c nvt.rst \
  -o md_LGV.out -r lgv.rst -info lgv.info

Path integral molecular dynamics (for quantum statistics)

(c) PRIMPIMD input using the Andersen thermostat with the “LFMiddle” scheme for the liquid water. Lengths of the bonds having hydrogen atoms are constrained.

Test: $AMBERHOME/test/middle-scheme/PIMD_Constr_Andersen_water

PRIMPIMD: NVT simulation of liquid water
&cntrl
ipimd = 1, nstlim = 10   ! PRIMPIMD for 10 steps
ntx = 5, irest = 0       ! read coordinates
temp0 = 300, tempi = 300 ! target temperature and initial temperature
dt = 0.002               ! time step in ps
cut = 7.0                ! non-bond cut-off
ischeme = 1,             !! leapfrog middle scheme
ithermostat = 2,         !! Andersen thermostat
therm_par = 8.0          !! thermostat parameter, in 1/ps
ig = 777                 ! random seed
ntc = 2,ntf = 2          ! constrain lengths of the bonds having hydrogen atoms
ntpr=1, ntwr=5, ntwx=5   ! output settings
/ 

(d) PRIMPIMD input using Langevin dynamics with the “LFMiddle” scheme for the liquid water. No constraints are applied.

Test: $AMBERHOME/test/middle-scheme/PIMD_Langevin_water

PRIMPIMD: NVT simulation of liquid water
&cntrl
ipimd = 1,nstlim = 10   ! PRIMPIMD for 10 steps
ntx = 5,irest = 0       ! read coordinates, 
                        ! and run as a new simulation.
temp0 = 300,tempi = 300 ! target and initial temperature
dt = 0.001              ! time step, in ps
cut = 7.0               ! non-bond cut off
ischeme = 1,            !! leapfrog middle scheme
ithermostat = 1,        !! Langevin thermostat
therm_par = 5.0         !! thermostat parameter, in 1/ps
ntc = 1                 ! no constraints, default
ntpr=1, ntwr=5, ntwx=5  ! output settings
/

When one runs PIMD in AMBER, a groupfile is needed. The groupfile gf_pimd may look like:

-O -i pimd.in -p qspcfw216.top -c nvt1.rst -o bead1.out -r bead1.rst 
 -x bead1.mdcrd -inf bead1.mdinfo 
-O -i pimd.in -p qspcfw216.top -c nvt2.rst -o bead2.out -r bead2.rst 
 -x bead2.mdcrd -inf bead2.mdinfo 
-O -i pimd.in -p qspcfw216.top -c nvt3.rst -o bead3.out -r bead3.rst 
 -x bead3.mdcrd -inf bead3.mdinfo 
-O -i pimd.in -p qspcfw216.top -c nvt4.rst -o bead4.out -r bead4.rst 
 -x bead4.mdcrd -inf bead4.mdinfo 

Note that each line starts with -O and ends with -inf <info> . The groupfile above contains 4 lines, which means 4 path integral beads are used.

sander.MPI is executed via the following command:

$ mpirun -np 8 sander.MPI -ng 4 -groupfile gf_pimd

The number of processes (8) that is specified by -np is a multiple of the number of groups (4). In this case 2 CPU processes are used on each path integral bead.

QM/MM molecular dynamics

(e) QM/MM MD input using the Langevin thermostat with the “LFMiddle” scheme for the alanine dipeptide solved in methol box. Lengths of the bonds having hydrogen atoms are constrained for the MM part, while no constraints are applied to the QM part.

Test: $AMBERHOME/test/middle-scheme/QMMM_Constr_ALA_Methol

constrained MD NVT: Alanine dipeptide in meoh (explicit solvent).
&cntrl
ipimd = 0, nstlim = 10,   ! MD for 10 steps
irest = 0, ntx = 1,       ! read coordinates
temp0 = 300               ! target temperature 
tempi = 300               ! initial temperature 
dt = 0.002,               ! time step, in ps 
cut = 8,                  ! non-bond cut off  
ig = 6666,                ! random seed for reproducing results 
ischeme = 1,              !! leapfrog middle scheme
ithermostat = 1           !! Langevin thermostat 
therm_par = 5.0,          !! thermostat parameter, in 1/ps 
ntc=2,ntf=2               ! constrain lengths of bonds having hydrogen atoms
 atoms
ntpr=1, ntwr=1, ntwx=1    ! output settings 
ifqnt=1                   ! switch on QM/MM coupled potential 
/ 
&qmmm qmmask=':ACE,ALA,NME', ! residues treated using QM 
qmcharge=0,            ! charge on QM region is 0 
qmshake=0,             ! no SHAKE for QM region 
qm_theory='PM3',       ! use the PM3 semi-empirical Hamiltonian 
qmcut=8.0              ! use 8 angstrom cut off for QM region 
/ 

One can run either a serial job (using sander ):

$ sander -O -i qmmm.in -p ala.top -c ala.crd -o qmmm.out -r qmmm.rst \
  -x qmmm.crd -info qmmm.info

or a parallel job (using sander.MPI ):

$ mpirun -np 4 sander.MPI -O -i qmmm.in -p ala.top -c ala.crd \
  -o qmmm.out -r qmmm.rst -x qmmm.crd -info qmmm.info

Replica exchange molecular dynamics

(f) REMD input using the Langevin thermostat with the “LFMiddle” scheme for the ACE-ALA-ALA-ALA-NME in vacuum.Lengths of the bonds having hydrogen atoms are constrained.

Test: $AMBERHOME/test/middle-scheme/REMD_Constr_ALA

Below is the input file for one of the replicas. The target temperatures are 300, 325, 350, and 400K for the 4 replicas, respectively.

REMD test with 4 replicas  
&cntrl         
imin = 0, nstlim = 100,     ! MD for 100 steps         
irest=1,ntx = 5,            ! read coordinates and velocities   
tempi = 0.0, temp0 = 300.0, ! initial and target temperature         
ischeme= 1,                 !! leapfrog middle scheme   
ithermostat = 1,            !! Langevin thermostat  
therm_par = 1.0,            !! thermostat parameter,in 1/ps   
dt = 0.002,                 ! time step, in ps  
ig=6666,                    ! random seed         
ntc = 2, ntf = 2,           ! constrain lengths of the bonds having hydrogen atoms         
ntwx = 50, ntwr =50, ntpr = 50,  ! output setting     
ntb=0,                      ! no periodicity         
cut = 99.0,                 ! non bond cut off         
numexchg=5,                 ! exchange frequency 
&end

When one runs REMD in AMBER, a groupfile is needed. The groupfile groupfile may look like:

-O -rem 1 -remlog rem.log -i rem.in.000 -p ala3.top -c mdrestrt -o rem.out.000
 -r rem.r.000 -inf reminfo.000
-O -rem 1 -remlog rem.log -i rem.in.001 -p ala3.top -c mdrestrt -o rem.out.001
 -r rem.r.001 -inf reminfo.001
-O -rem 1 -remlog rem.log -i rem.in.002 -p ala3.top -c mdrestrt -o rem.out.002
 -r rem.r.002 -inf reminfo.002
-O -rem 1 -remlog rem.log -i rem.in.003 -p ala3.top -c mdrestrt -o rem.out.003
 -r rem.r.003 -inf reminfo.003 

Note that each line starts with -O and ends with -inf <info> . The groupfile has 4 lines, which means 4 replicas are employed in REMD. sander.MPI is executed via the following command:

$ mpirun -np 4 sander.MPI -ng 4 -groupfile groupfile

The number of processes (4) that is specified by -np can be replaced by any multiple of the number of replicas used in REMD (4 in this case).

Pmemd examples for CPU and GPU

(g) Pmemd input for liquid water with 4096 water molecules in a cell

Test: $AMBRHOME/test/middle-scheme/4096wat

MD: NVT simulation of liquid water
&cntrl
ntx = 5, irest = 1, ! read coordinates
ntc = 2, ntf = 2, ! constrain lengths of bonds
tol = 0.0000001, ! having hydrogen atoms
nstlim = 10, ! MD for 10 steps
ntpr = 1, ntwr = 10000 ! output settings
dt = 0.001, ! timestep in ps
ig = 71277, ! random seed
cut = 7.0, ! non-bond cut off
ischeme = 1, !! Leapfrog middle scheme
ithermostat = 1, !! Langevin thermostat
therm_par = 5.0, !! thermostat parameter
midpoint = 1 ! use midpoint method, only for pmemd.MPI;
! remove this flag otherwise
/
&end

(h) Pmemd input for DNA A7 −T7 duplex solvated in a cell with 3351 water molecules and 12 sodium ions as counter-ions

Test: $AMBERHOME/test/middle-scheme/DNA7

MD: NVT simulation of DNA duplex
&cntrl
ntx = 5, irest = 1, ! read coordinates
ntc = 2, ntf = 2, ! constrain lengths of bonds
tol = 0.0000001, ! having hydrogen atoms
nstlim = 10, ! MD for 10 steps
ntpr = 1, ntwr = 10000 ! output settings
dt = 0.001, ! timestep in ps
ig = 71277, ! random seed
cut = 9.0, ! non-bond cut off
ischeme = 1, !! Leapfrog middle scheme
ithermostat = 1, !! Langevin thermostat
therm_par = 5.0, !! thermostat parameter
midpoint = 1 ! use midpoint method, only for pmemd.MPI;
! remove this flag otherwise
/

(i) Pmemd input for Ethaline Deep Eutectic Solvent(512 choline chloride and 1024 ethylene glycol molecules)

Test: $AMBERHOME/test/middle-scheme/ETH

MD: NVT simulation of Ethaline Deep Eutectic Solvent
&cntrl
ntx = 5, irest = 1, ! read coordinates
ntc = 2, ntf = 2, ! constrain lengths of bonds
tol = 0.0000001, ! having hydrogen atoms
nstlim = 10, ! MD for 10 steps
ntpr = 1, ntwr = 10000 ! output settings
dt = 0.001, ! timestep in ps
ig = 71277, ! random seed
cut = 10.0, ! non-bond cut off
ischeme = 1, !! Leapfrog middle scheme
ithermostat = 1, !! Langevin thermostat
therm_par = 5.0, !! thermostat parameter
midpoint = 1 ! use midpoint method, only for pmemd.MPI;
! remove this flag otherwise
/

&end

All those three pmemd examples can be executed in a serial job(using sander or pmemd):

$ sander -O -i mdin -o mdout -p prmtop -c inpcrd -r restrt

or

$ pmemd -O -i mdin -o mdout -p prmtop -c inpcrd -r restrt

or a parallel job(using sander.MPI or pmemd.MPI):

$ sander.MPI -O -i mdin -o mdout -p prmtop -c inpcrd -r restrt

or

$ pmemd.MPI -O -i mdin -o mdout -p prmtop -c inpcrd -r restrt

or a GPU-accelerated job(using pmemd.cuda):

$ pmemd.cuda_DPFP -O -i mdin -o mdout -p prmtop -c inpcrd -r restrt
$ pmemd.cuda_SPFP -O -i mdin -o mdout -p prmtop -c inpcrd -r restrt

References

Related links:

Jian Liu research group (Peking University)

AMBER (2018/2019 versions).