Molecular Dynamics

Introduction

QuantumATK makes Molecular Dynamics (MD) simulations very simple: Just add the desired calculator and a Molecular Dynamics block to your configuration and start the simulation! There are, however, some guidelines regarding which parameters and settings are the most suitable ones for different types of simulations.

MD simulation is a technique to simulate the motion of atoms and molecules under predefined conditions, such as temperature, pressure, stress, external forces, etc. MD simulations can therefore be used to study dynamical processes at the nanoscale and to calculate a broad range of properties, e.g. phase diagrams, diffusion coefficients, or various response functions, as well as static quantities such as radial distribution functions, coordination numbers, elastic moduli, etc.

Methodology

MD essentially utilizes the numerical solution of Newton’s equations of motion for a set of atoms from a given initial configuration. This is commonly achieved via numerical integration by discretizing time into small intervals called the time step using integrators such as the velocity verlet algorithm [1]. The interactions between the atoms, i.e. the interatomic forces, can be calculated based on various methods, ranging from Density Functional Theory (DFT) to classical potentials. These forces determine the acceleration of the atoms and allow the positions and velocities to propagate towards the next time step. Repeating this procedure many times yields a series of snapshots, describing the trajectory of the system in phase space, which can be analyzed to extract the desired properties.

Before setting up the simulation you need to decide what type of calculation you are interested in. Should the total energy be conserved, as in an isolated system? Should the temperature be kept constant to mimic the coupling of the system to a heat bath? Is the system exposed to any external pressure or stress? Based on these considerations, a suitable set of simulation parameters should be selected: Time step size, number of integration steps (duration of simulation), integration algorithm, initial temperature, constraints, etc.

Similar to a real experiment, some empirical knowledge and experience is required for performing MD simulations. In these notes, you will get to know the basic ingredients of MD simulations for running MD simulations with the QuantumATK MolecularDynamics module. A tutorial to try out MD can be found in this link: How to Setup Basic Molecular Dynamics Simulations.

NVE Simulations

As mentioned previously, MD simulations are based on the solution of Newton’s equations of motion. In the purest form, an MD simulation reproduces an NVE ensemble (also called the microcanonical ensemble) where the number of atoms N, the volume V, and the total energy E are conserved. These conditions correspond to a completely isolated system.

Setting Up the Geometry

The primitive unit cell of the material to be simulated could have non-orthogonal cell vectors. Although it is possible to perform MD simulation on such cells, it is often more convenient to use orthogonal cells (in particular when the cell size is supposed to change during the simulation).

Furthermore, it is useful to increase the size of the structure to include more atoms. This generally improves the average of the measured observables and reduces finite size effects due to interaction of atoms in a small simulation cell with their periodic images. (Ideally, the lengths of the cell vectors should be larger than twice the interaction range of the potential).

Setting up the Calculation Script

In the MolecularDynamics, the value of Log Interval determines how often the snapshots of the current configuration are written to disk. A too small value can decrease the speed of the simulation and, in the case of large systems, produce huge output files.

Tip

In most MD applications it is not necessary to record snapshots at a high frequency, unless you explicitly want to analyze high-frequency oscillations.

Setting Initial Velocity type to Maxwell-Boltzmann at a Temperature of T Kelvin will create the random initial velocities of the atoms drawn from a Maxwell-Boltzmann distribution corresponding to a temperature of T Kelvin. It is also useful to remove the centre of mass motion of the cell. The Time step is a crucial parameter in MD simulations as it determines the accuracy and efficiency of the numerical integration scheme.

The central point is that a larger time step, although it may at first glance improve the efficiency of the simulation, increases the error in the numerical integration scheme of the equations of motion. The underlying assumption that the atomic forces are approximately constant during one integration step is not valid any more. Essentially, the chosen time step should be small enough to resolve the highest vibrational frequencies of the atoms (i.e. it should be much smaller than the smallest vibrational period), so if you have light atoms (e.g. hydrogen), you will generally be required to use a smaller time step than if you have only heavy atoms (such as gold). A smaller time step size may also be necessary if you have different elements in your calculation, if the temperature is high, or if the atoms are far away from their equilibrium configuration, i.e if large forces act on the particles. For most systems a safe choice to start with, if you do not know what time step to use, is 1 fs. Larger time step values can then be assessed by monitoring the conservation of the total energy in an NVE simulation under the conditions of interest.

Another very important aspect of performing MD simulations is to be aware that it may take some time before a system is equilibrated to the chosen external temperature, pressure, etc. This is important because any measurement of observables should only be carried out after the system is equilibrated (unless you are specifically interested in non-equilibrium phenomena). In general, equilibration times of MD simulations can vary over a wide range up to several hundreds of nanoseconds (e.g. in biological molecules or polymer systems) and depend essentially on the longest relaxation time present in the system. You should therefore always check carefully that the observables of interest have reached a stationary state in your simulations.

NVT Simulations

The NVT ensemble is also called the canonical ensemble and the main difference to the microcanonical NVE ensemble is that the system is not isolated any more but can exchange energy with a surrounding virtual heat bath.

Setting Up the Script

QuantumATK offers four possibilities for NVT simulations: 1) NVT Berendsen [2] , 2) NVT Nose Hoover [3] , 3) Langevin [4] and 4) NVT Bussi Donadio Parinello [5]. Each of them implements the coupling to the heat bath in a different way. Although the best choice depends on the system and the particular application at hand, the Nose Hoover algorithm can generally be considered the most reliable method for this task.

NVT Nose Hoover Simulations

The Thermostat timescale value determines how quickly the system temperature approaches the reservoir temperature. The Heating rate parameter can be used to linearly increase or decrease the temperature during the simulation (cf. NVTNoseHoover in the reference manual). The Thermostat chain length parameter determines how many subsequent thermostats are invoked to control the temperature. The default value of 3 should work well in most cases - only if you experience a strong, persistent oscillation in the temperature should this value be increased to achieve a more stable simulation.

As mentioned, how fast the temperature of the system approaches the target temperature can be controlled via the thermostat timescale parameter. In general, a small thermostat timescale parameter causes a tighter coupling to the virtual heat bath and a system temperature that closely follows the reservoir temperature. Such a tight thermostat coupling is typically accompanied by a more pronounced interference with the natural dynamics of the particles, however. Therefore, when aiming at a precise measurement of dynamical properties, such as diffusion or vibrations, you should use a larger value of the thermal coupling constant or consider running the simulation in the NVE ensemble to avoid artifacts of the thermostat.

NVT Berendsen and Langevin Simulations

Aside from the Nose-Hoover chain thermostat you can also select two other NVT algorithms, NVT Berendsen (cf. NVTBerendsen in the reference manual) and Langevin (cf. Langevin in the reference manual).

The Berendsen thermostat implements an algorithm which is in some situations more stable, as it effectively suppresses temperature oscillations. The drawback is that in contrast to the Nose-Hoover thermostat, it does not exactly reproduce a canonical ensemble, although the deviations are typically small. Effectively this means that observables (such as the velocity distribution) do not have the correct distribution. For this reason it should, if at all, only be used for equilibration simulations, primarily in cases where a more robust temperature control is required.

The third thermostat type that you can choose for NVT simulations is the Langevin algorithm. This algorithm solves the Langevin equation [4], which explicitly includes friction as well as stochastic collision forces into Newton’s equations of motion to mimic the interaction with particles of the heat bath.

Unlike the other thermostat types, the Langevin thermostat individually couples each particle to the heat bath. This produces a very tight coupling, as if the system were immersed in a virtual, viscous medium, but it also suppresses the natural dynamics in a more pronounced way. For this reason, the Langevin thermostat should primarily be used to generate structures or to sample an ensemble rather than to calculate dynamical properties.

In a Langevin NVT simulation, the degree of coupling can be set by the Friction parameter. It is defined in inverse time units, meaning that it works in the opposite way compared to the coupling parameters of the Berendsen or Nose Hoover thermostat. Increasing the friction value therefore causes a tighter coupling to the heat bath and a more severe modification of the dynamics.

NVT Bussi Donadio Parrinello

As we mentioned earlier, NVT Berendsen thermostat doesn’t sample the correct canonical ensemble. However, a stochastic variant of NVT Berendsen is the NVT Bussi Donadio Parrinello thermostat (cf. NVTBussiDonadioParrinello in the reference manual), which correctly samples the canonical ensemble along with the other advantages of NVT Berendsen thermostat.

NPT Simulations

There are three possibilities to run isothermic-isobaric (NPT) simulations: The NPT Berendsen [2] (cf. NPTBerendsen in the reference manual), the NPT Martyna Tobias Klein [6] (cf. NPTMartynaTobiasKlein in the reference manual) and the NPT Bernetti Bussi [7] algorithm (cf. NPTBernettiBussi in the reference manual). The NPT Berendsen method works via the same principle as the corresponding NVT integrator and therefore does not exactly reproduce the correct physical ensemble, although, again, the deviations are typically small for large unit cells.

For production simulations, it is recommended not to use the NPT Berendsen barostat. The NPT Bernetti Bussi method is usually a better choice. It is a stochastic variant of the Berendsen barostat which rescales the unit cell stochasticaly so as to sample the NPT ensemble properly even for very small unit cells.

Similar to the thermostat time scale, the Barostat time scale parameter can be used to set how quickly the system pressure approaches and oscillates around the target pressure.

In the workflow block in NanoLab, the Reservoir Pressure group box allows you to set the pressure and choose whether you want to use isotropic or anisotropic pressure coupling. Isotropic pressure is a suitable choice if isotropic systems, such as liquids or crystals with cubic symmetry, are simulated. Here, the pressure is calculated as \(P = (P_{xx} + P_{yy} + P_{zz})/3\), with \(P_{ij} = - \sigma _{ij}\) being the components of the pressure tensor, obtained from the stress tensor \(\sigma\). All cell vectors will be changed by the same factor, which preserves the shape of the cell. For anisotropic systems, which may have different moduli in different directions, you should uncheck the Isotropic Pressure box. This means that all active components \(P_{ij}\) of the pressure tensor pressure are coupled independently, and the cell vectors can respond separately to the external pressure. Individual components can be switched on and off via the Pressure Coupling check boxes.

Moreover, instead of a single external pressure value, you can also specify independent reference stress values for each component to simulate for instance mechanical shear or creep experiments. The NPTMartynaTobiasKlein entry in the QuantumATK Reference Manual and the tutorial Simulating a creep experiment of polycrystalline copper provide a more detailed explanation of this functionality.

A typical application of NPT simulations is to simulate phase transitions, e.g. from \(\alpha\)- to \(\beta\)-quartz, by performing a series of simulations at different temperatures and measuring the average cell volume as a function of the temperature, as reported in Ref. [8].

Non-Equilibrium Simulations

In non-equilibrium simulations, a temperature gradient is enforced or strain is applied to study a material’s behaviour out of equilibrium.

The NonEquilibriumMomentumExchange algorithm simulates a heat flow from a hot region to a cold region and can be used to calculate the thermal conductivity (cf. NonEquilibriumMomentumExchange in the reference manual for a usage example).

In a stress-strain calculation the cell is extended or strained in one direction and the resulting stress in that direction is then recorded(cf. StrainConfigurationHook in the reference manual for a usage example).

MD Simulations with Constraints

It is possible to fix the positions of selected atoms while running MD simulations. The easiest way to set such constraints is to add a tag to the selected atoms to be constrained in the configuration (using the Builder tool or the Configuration block in the Workflow Builder).

Then, in the Molecular Dynamics widget, click on Add Constraints to select the tag and apply the needed constraints. You can select between Fixed atoms and Rigid body constraints. The former type fixes all atoms in the respective group to their initial positions or constrains their movement in the X, Y, or Z direction by selecting Fixed X, Fixed Y, or Fixed Z, respectively. The Rigid constraint type treats all atoms within this group as a rigid body. This means that these atoms move collectively according to the force on their center of mass. In the current implementation only translational motions are possible, whereas rigid body rotations are deactivated.

When using constraints, at least one atom must remain unconstrained. The thermostat will automatically adapt to the reduced number of degrees of freedom due to the constraints.

Note, however, that some features in the MD Analyzer are not designed to work with constraints.

Note

It is generally not recommended to perform constant pressure/stress (i.e. NPT) simulations in the presence of constrained atoms, as the stress contribution arising within the constrained parts of the system is not removed from the overall stress tensor.

Device Configurations

Performing an NVE/NVT simulation of a DeviceConfiguration essentially works in the same way as for a BulkConfiguration. The atoms inside the electrode extensions in the central region are automatically fixed by constraints to preserve compatibility with the electrodes. Although it is technically also possible to perform constant pressure/stress (i.e. NPT) simulations of DeviceConfigurations, such calculations are not recommended, due to the presence of constraints.

Molecules

Isolated molecules in gas phase can be simulated under NVE/NVT ensemble. Molecule simulations can be setup within MoleculeConfiguration or BulkConfiguration. When using BulkConfiguration, care must be taken to specify a box size large enough to account for the finite size effects when periodic boundary condition is applied.

References

[1]

W. C. Swope, H.C. Andersen, P.H. Berens, K.R. Wilson: A computer simulation method for the calculation of equilibrium constants for the formation of physical clusters of molecules: Application to small water clusters. J. Chem. Phys. 76, 637 (1982)

[2] (1,2)

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 (1982).

[3]

G. J. Martyna, M. L. Klein, M. Tuckerman: Nosé–Hoover chains: The canonical ensemble via continuous dynamics. J. Chem. Phys. 97, 2635 (1992)

[4] (1,2)

W. F. van Gunsteren, H. J. C. Berendsen: A leap-frog algorithm for stochastic dynamics. Mol. Sim. 1, 173 (1988)

[5]

Giovanni Bussi, Davide Donadio, and Michele Parrinello: Canonical sampling through velocity rescaling. J. Chem. Phys. 126(1):014101, (2007).

[6]

G. J. Martyna, D. J. Tobias, and M. L. Klein: Constant pressure molecular dynamics algorithms. J. Chem. Phys. 101, 4177 (1994)

[7]

M. Bernetti, and G. Bussi: Pressure control using stochastic cell rescaling. J. Chem. Phys. 153, 114107 (2020)

[8]

M. H. Müser, and K. Binder: Molecular dynamics study of the α-β transition in quartz: elastic properties, finite size effects, and hysteresis in the local structure. Phys. Chem. Minerals 28, 746 (2001)