Molecular Dynamics

CONQUEST can perform molecular dynamics both when the density matrix is computed using diagonalisation and O(N), the latter allowing dynamical simulations of (but not limited to) tens of thousands of atoms. The equations of motion are integrated using the velocity Verlet method in the case of the microcanonical ensemble (NVE), and modifications thereof for the canonical (NVT) and isobaric-isothermal (NPT) ensembles, the details of which can be found in Molecular Dynamics: Theory. In addition to converging the parameters for the electronic structure calculations, the following points must also be considered.

Go to top.

Self-consistency tolerance and XL-BOMD

The convergence of the electronic structure is important in MD, as insufficient convergence can be responsible for “drift” in the conserved quantity of the dynamics. Although the molecular dynamics integrators used in CONQUEST are time reversible, the SCF procedure is not. Therefore tight convergence (minE.SCTolerance for diagonalisation, minE.LTolerance for linear scaling) is necessary. In the case of diagonalisation, an SCF tolerance of 1E-6 is typically enough to negate the drift. However, extended-Lagrangian Born-Oppenheimer MD (XL-BOMD) [MD1], currently only implemented for O(N), essentially makes the SCF component of the MD time-reversible by adding the electronic degrees of freedom to the Lagrangian, relaxing the constraint on minE.LTolerance — although it is still somewhat dependent on the ensemble. In the NVE and NVT ensembles, a L-tolerance of 1E-5 has been found to be sufficient to give good energy conservation, decreasing to 1E-6 in the NPT ensemble. The following flags are required for XL-BOMD:

DM.SolutionMethod ordern
AtomMove.ExtendedLagrangian T

Go to top.

Restarting

Assuming the calculation ended gracefully, it can easily be restarted by setting,

AtomMove.RestartRun T

This will do several things: it will read the atomic coordinates from md.position and read the md.checkpoint file, which contains the velocities and extended-system (Nose-Hoover chain and cell) variables. Depending on the value of DM.SolutionMethod, it will read the K-matrix files (diagon) or the L-matrix files (ordern). Finally, it will append new data to the md.stats and md.frames files, but it will overwrite all other files, including Conquest_out.

Unless explicitly overridden, AtomMove.RestartRun T changes the defaults of the following keywords to T:

General.LoadDM T
SC.MakeInitialChargeFromK T
AtomMove.ReadVelocity T

For XL-BOMD it also defaults XL.LoadX to T. For MSSF and blip-basis calculations it defaults Basis.LoadCoeffs to T, so that the saved support-function coefficients are loaded.

In addition to the files mentioned above, CONQUEST will try to read the K-matrix from Kmatrix2.i00.* when using diagonalisation or the L-matrix from Lmatrix2.i00.* when using O(N), and Xmatrix2.i0*.* if the extended-Lagrangian formalism is used. Note that metadata for these files is stored in InfoGlobal.i00 which is also required when restarting. If the calculation ended by hitting the walltime limit, the writing of these matrix files may have been interrupted, rendering them unusable. In this case, set General.LoadDM F and, for XL-BOMD, XL.LoadX F after setting AtomMove.RestartRun T. Setting General.MaxTime to some number of seconds less (say 30 minutes) than the calculation wall time limit will force the calculation to stop gracefully, preventing the aforementioned situation.

Go to top.

Visualising the trajectory

Setting the flag AtomMove.WriteXSF T dumps the coordinates to the file trajectory.xsf every AtomMove.XsfFreq steps (which defaults to AtomMove.OutputFreq). The .xsf file can be read using VMD. A small VMD script, view.vmd is included with the code, and can be invoked using,

vmd -e view.vmd

assuming the vmd executable is in your path.

Go to top.

TDEP output

CONQUEST molecular dynamics data can be used to perform lattice dynamical calculations using the Temperature Dependent Effective Potential (TDEP) code. Setting the flag MD.TDEP T will make conquest dump configurations, forces and metadata in a format readable by TDEP.

Go to top.

Non-Hamiltonian dynamics

Canonical (NVT) ensemble

The thermostat is set using the MD.Thermostat flag, and can take the values svr (stochastic velocity rescaling) and nhc (Nose-Hoover chain). These thermostats generate the correct canonical ensemble phase space distribution, and both give a conserved quantity that allows the quality of the dynamics to be monitored.

  1. Stochastic velocity rescaling

AtomMove.IonTemperature 300.0
MD.Ensemble nvt
MD.Thermostat svr
MD.tauT 10

While the NHC uses chaotic sensitivity to initial conditions to achieve better ergodicity, the SVR thermostat [MD2] uses a judiciously chosen stochastic force coupled to a weak scaling thermostat to correctly generate the canonical phase space distribution. The MD.tauT parameter gives the coupling timescale; the velocity scaling factor is modified by a factor \(\Delta t/\tau\), so a larger \(\tau\) results in a more slowly varying temperature. While some characterisation of the system is recommended, values of \(\tau\) around 20–200fs are reasonable. To reproduce a simulation, the random number generator seed can be set with the General.RNGSeed <integer> flag.

  1. Nose-Hoover chain

AtomMove.IonTemperature 300.0
MD.Ensemble nvt
MD.Thermostat nhc
MD.nNHC 5
MD.nYoshida 5
MD.tauT 30

When thermostatting using a Nose-Hoover chain [MD3, MD4, MD5], it may be necessary to set a couple more flags. MD.nNHC sets the number of thermostats in the chain (the default of 5 is generally sensible), and MD.nYoshida determines the order of Yoshida-Suzuki integration. This is essentially a higher level integration scheme that can improve energy conservation in cases when rapid changes in the Nose-Hoover thermostat velocity is causing integration errors. Note that MD.tauT means something different to the SVR case. A good guess is the time period of the highest frequency motion of the system in fs; however, in the NVT ensemble, the energy conservation is not very sensitive to this value. The NHC masses can also be set manually using the following block.

MD.CalculateXLMass F
MD.nNHC 5
%block MD.NHCmass
  5 1 1 1 1
%endblock

Go to top.

Isobaric-Isothermal (NPT) ensemble

Two extended-system barostats are implemented: Parrinello-Rahman [MD6] and Martyna-Tobias-Tuckerman-Klein (MTTK) [Ta1].

  1. Parrinello-Rahman

AtomMove.IonTemperature 300.0
AtomMove.TargetPressure 10.0
MD.Ensemble npt
MD.Thermostat nhc
MD.Barostat pr
MD.nNHC 5
MD.nYoshida 5
MD.tauT 100
MD.tauP 200
MD.PDrag 10.0

The Parrinello-Rahman barostat generates the correct ensemble, but can be subject to low frequency “ringing” fluctuations in the temperature and pressure that can destabilise the system or slow equilibration. Unlike in the NVT ensemble, this combination of barostat and thermostat is very sensitive to the choice of both MD.tauT and MD.tauP; note that their values are somewhat higher in this case, since integration errors in the NHC tend to be more severe due to coupling of the cell and atomic motions. They are dependent on the system, so it is advised that you find a combination of these parameters that gives the best energy conservation. The cell is thermostatted using a separate Nose-Hoover chain to the atoms by default, but they can be controlled with the same chain by setting MD.CellNHC F. The ad hoc damping controlled by MD.PDrag acts on the cell/barostat velocities and, when present, the separate cell-thermostat velocities. Each application uses the factor

\[1 - \frac{D_P\,\Delta t}{\tau_P n_{\mathrm{MTS}} n_{\mathrm{YS}}},\]

where \(D_P\) is MD.PDrag. This factor is applied within the integration substeps, potentially more than once per full MD timestep. A nonzero drag perturbs the formal NPT dynamics, but can help damp ringing fluctuations.

Note that the NPT ensemble can also be generated correctly by thermostatting using the SVR thermostat, although the meaning of the parameter MD.tauT is different in this case, as in NVT dynamics.

  1. Martyna-Tobias-Tuckerman-Klein

The MTTK barostat is implemented for isotropic volume fluctuations and is coupled to the Nose-Hoover-chain thermostat:

AtomMove.IonTemperature 300.0
AtomMove.TargetPressure 10.0
MD.Ensemble npt
MD.Thermostat nhc
MD.Barostat mttk
MD.CellConstraint volume
MD.tauT 100
MD.tauP 200

Use the Parrinello-Rahman barostat instead when the three orthorhombic cell lengths must vary independently with MD.CellConstraint xyz.

Postprocessing tools

The Python utility supplied for analysing CONQUEST molecular-dynamics output is described in Molecular-dynamics analysis.

Go to top.

[MD1]

A. M. N. Niklasson. Extended Born-Oppenheimer Molecular Dynamics. Phys. Rev. Lett., 100:123004, 2008. doi:10.1103/PhysRevLett.100.123004.

[MD2]

G. Bussi, D. Donadio, and M. Parrinello. Canonical sampling through velocity rescaling. J. Chem. Phys., 126:014101, 2007. doi:10.1063/1.2408420.

[MD3]

S. Nosé. A unified formulation of the constant temperature molecular dynamics methods. J. Chem. Phys., 81:511, 1984. doi:10.1063/1.447334.

[MD4]

W. G. Hoover. Canonical dynamics: Equilibrium phase-space distributions. Phys. Rev. A, 31:1695, 1985. doi:10.1103/PhysRevA.31.1695.

[MD5]

G. J. Martyna, M. L. Klein, and M. Tuckerman. Nosé–hoover chains: the canonical ensemble via continuous dynamics. J. Chem. Phys., 97:2635, 1992. doi:10.1063/1.463940.

[MD6]

M. Parrinello and A. Rahman. Polymorphic transitions in single crystals: A new molecular dynamics method. J. Appl. Phys., 52:7182–7190, December 1981. doi:10.1063/1.328693.

Go to top.