Molecular Dynamics

ASAP provides a flexible GUI to facilitate all the steps involved in the workflow of Molecular Dynamics simulations (MD).

ASAP provides two molecular dynamics drivers. They differ in where the equations of motion are integrated, not in the accuracy of the forces:

• MD-Multicode: The dynamics are integrated by ASE, while the forces at each step are supplied by an external calculator. Because the driver does not depend on the code providing the forces, it can be combined with the SIESTA, Quantum Espresso and MACE calculators (hence multicode).

• MD-SIESTA: The dynamics are integrated by SIESTA itself and are therefore available only with the SIESTA calculator.

Each driver offers a different set of ensembles and algorithms; see MD Algorithms (SIESTA only) and MD Algorithms (Multicode) below. With the Quantum Espresso or MACE calculators, MD-Multicode is the only option.

With the SIESTA calculator both are available, and the choice depends on what your simulation requires:

• MD-SIESTA reuses the wavefunctions and density matrix from the previous MD step as the starting point for the next one, which reduces the cost per step. Additionally, it avoids reinitializing the system infrastructure at each step, further speeding up simulation runs. It is also the only driver that provides simulated annealing.

• MD-Multicode offers finer control over the cell dynamics: you can choose the barostat mode, constrain individual cell degrees of freedom and apply pressure selectively along chosen axes.

See MD Algorithms (SIESTA only) and MD Algorithms (Multicode) for the full list of algorithms available in each.

The project starts with the definition of the atomic structure of the system. We refer the user to section Structure Modeling in ASAP to find more information on the ASAP structure builder/viewer interface.

Once the atomic structure is defined, select Molecular Dynamics from the list of project types.

After selecting the project type, click on the Parameters icon to open the Parameter widget.

You can use the Parameter widget to set up the input parameters related to the selected type of project (Settings), edit the masses, magnetic moments, atomic positions and velocities of the system (Variables).

We describe below the settings and variables relevant for MD calculations, e.g. Velocity Distribution, as it affects the systems dynamics. We refer the user to chapter Parameters for more information on all parameters that can be tuned in the Parameter widget.

Settings

The settings for each chosen ensemble are different, however, all MD ensemble share two options:

• Time step: Controls the length of the time step in the MD simulation. If the time step is too large, the dynamics of the system will show unphysical behavior: energy is not conserved, high-frequency modes become populated, and the system may even melt. It is always safer to choose a small time step, at the cost of computation time. The default value is 1 fs, a standard choice for most systems. Systems with light atoms (e.g. hydrogen) and/or strong bonds (e.g. C-C, C-H) require a smaller value, while systems composed of heavy elements can often tolerate a larger one. We advise examining the system under study and verifying energy conservation in a short NVE run before committing to a production value.

  The time units are selected at the right of this parameter. The options are femto-seconds (fs), i.e. 10\(^{-15}\) seconds, and pico-seconds (ps), i.e. 10\(^{-12}\) s.

• number of steps: Number of steps of the MD simulation. The number of steps multiplied by the time step will give a full-time MD calculation.

ASAP currently has different statistical ensembles integrated by the algorithms described in the following two sections.

MD Algorithms (SIESTA only)

Standard Verlet algorithm (NVE)

Microcanonical ensemble, where the number of atoms (N), volume (V) and total energy (E) are conserved.

• Initial temperature: Initial temperature for the MD run. Represents the initial kinetic energy distribution of the particles. The atoms are assigned random velocities drawn from the Maxwell-Boltzmann distribution with the corresponding temperature. The constraint of zero center of mass velocity is imposed. In units of K, eV or meV.

This algorithm assumes the system is isolated, this is especially useful when studying systems where energy conservation is crucial, such as understanding the long-term behavior of isolated systems or exploring energy landscapes.

Nosé thermostat (NVT)

MD with temperature controlled by means of a Nosé thermostat. Only available for MD-SIESTA

• Initial temperature: Initial temperature for the MD run. The atoms are assigned random velocities drawn from the Maxwell-Boltzmann distribution with the corresponding temperature. The constraint of zero center of mass velocity is imposed. In units of K, eV or meV.

• Target temperature: Target temperature for Nosé thermostat. In units of K, eV or meV.

• Nosé mass: Generalised mass of Nosé variable. This determines the time scale of the Nosé variable dynamics, and the coupling of the thermal bath to the physical system.

Parrinello-Rahman (NP)

MD with pressure controlled by the Parrinello-Rahman method. Only available for MD-SIESTA

• Initial temperature. Initial temperature for the MD run. The atoms are assigned random velocities drawn from the Maxwell-Boltzmann distribution with the corresponding temperature. The constraint of zero center of mass velocity is imposed (in units of K, eV or meV).

• Target pressure. Target pressure for Parrinello-Rahman method. In units of GPa, eV/Ångström\(^{3}\) or Ry/Bohr\(^{3}\).

• Parrinello-Rahman mass. Generalised mass of Parrinello-Rahman variable. This determines the time scale of the Parrinello-Rahman variable dynamics, and its coupling to the physical system.

Nosé Parrinello-Rahman (NPT)

Isothermal–isobaric ensemble, where the number of atoms (N), pressure (P) and temperature (T) are conserved. MD with temperature is controlled by means of a Nosé thermostat and pressure controlled by the Parrinello-Rahman method.

If your system does not have a cell when you try to select this algorithm the following warning will pop up:

Click Yes to continue with the default unit cell or define one in the 3D structure editor.

• Initial temperature: Initial temperature for the MD run. The atoms are assigned random velocities drawn from the Maxwell-Boltzmann distribution with the corresponding temperature. The constraint of zero center of mass velocity is imposed (in units of K, eV or meV).

• Target temperature: Target temperature for the Nosé thermostat. In units of K, eV or meV.

• Nosé mass: Generalised mass of Nosé variable. This determines the time scale of the Nosé variable dynamics, and the coupling of the thermal bath to the physical system.

• Target pressure: Target pressure for Parrinello-Rahman method. In units of GPa, eV/Ångström or Ry/Bohr.

• Parrinello-Rahman mass: Generalised mass of Parrinello-Rahman variable. This determines the time scale of the Parrinello-Rahman variable dynamics, and its coupling to the physical system.

Annealing

MD with annealing to a desired temperature and/or pressure. Only available for MD-SIESTA

Anneal. option: The annealing options are desired temperature and/or pressure.

Temperature

If the Temperature is selected, the additional settings are:

• Initial temperature: Initial temperature for the MD run. The atoms are assigned random velocities drawn from the Maxwell-Boltzmann distribution with the corresponding temperature. The constraint of zero center of mass velocity is imposed (in units of K, eV or meV).

• Target temperature: Target temperature for annealing MD. In units of K, eV or meV.

• Tau relax: Relaxation time to reach target temperature and/or pressure in annealing MD. Note that this is a relaxation time, and as such it gives a rough estimate of the time needed to achieve the given targets. As a normal simulation also exhibits oscillations, the actual time needed to reach the averaged targets will be significantly longer.

Pressure

If the Pressure is selected, the additional settings are:

• Initial temperature: Initial temperature for the MD run. The atoms are assigned random velocities drawn from the Maxwell-Boltzmann distribution with the corresponding temperature. The constraint of zero center of mass velocity is imposed (in units of K, eV or meV).

• Tau relax: Relaxation time to reach target temperature and/or pressure in annealing MD. Note that this is a relaxation time, and as such it gives a rough estimate of the time needed to achieve the given targets. As a normal simulation also exhibits oscillations, the actual time needed to reach the averaged targets will be significantly longer. In units of fs or ps.

• Target pressure: The target pressure is achieved by velocity and unit cell rescaling, in a given time determined by the variable Tau relax. In units of GPa, eV/Ångström\(^{3}\) or Ry/Bohr\(^{3}\).

• Bulk modulus estimation: Estimation (may be rough) of the bulk modulus of the system. This is needed to set the rate of change of cell shape to reach target pressure. In units of GPa, eV/Ångström\(^{3}\) or Ry/Bohr\(^{3}\).

Temperature and pressure

If both Temperature and Pressure are selected, the additional settings are:

• Initial temperature: Initial temperature for the MD run. The atoms are assigned random velocities drawn from the Maxwell-Boltzmann distribution with the corresponding temperature. The constraint of zero center of mass velocity is imposed (in units of K, eV or meV).

• Target temperature: The target temperature is achieved by velocity and unit cell rescaling, in a given time determined by the variable Tau relax. In units of K, eV or meV.

• Tau relax: Relaxation time to reach target temperature and/or pressure in annealing MD. Note that this is a relaxation time, and as such it gives a rough estimate of the time needed to achieve the given targets. As a normal simulation also exhibits oscillations, the actual time needed to reach the averaged targets will be significantly longer. In units of fs or ps.

• Target pressure: The target pressure is achieved by velocity and unit cell rescaling, in a given time determined by the variable Tau relax. In units of GPa, eV/Ångström\(^{3}\) or Ry/Bohr\(^{3}\).

• Bulk modulus estimation: Estimation (may be rough) of the bulk modulus of the system. This is needed to set the rate of change of cell shape to reach target pressure. In units of GPa, eV/Ångström\(^{3}\) or Ry/Bohr\(^{3}\).

MD Algorithms (Multicode)

Important: Initial velocities for MD (Multicode) are set through the Variables–>Velocity menu, as explained in section Variables (Velocities).

Verlet (NVE)

Microcanonical ensemble, where the number of atoms (N), volume (V) and total energy (E) are conserved.

• Logging frequency: defines how frequently the information of a running MD calculation is stored. During the MD simulation, atomic positions and velocities are stored. This is very useful to visualise and analyse the calculations.

This algorithm assumes the system is isolated, this is especially useful when studying systems where energy conservation is crucial, such as understanding the long-term behavior of isolated systems or exploring energy landscapes.

NVT Melchionna

Melchionna NVT is simply Melchionna NPT with the barostat disabled, leaving only the thermostat active.

If your system does not have a cell when you try to select this algorithm the following warning will pop up:

If your system is a bulk but the simulation box is not an upper triangular matrix please adjust the simulation box:

• Logging frequency: defines how frequently the information of a running MD calculation is stored. During the MD simulation, atomic positions and velocities are stored. This is very useful to visualise and analyse the calculations.

• Temperature: Target temperature for the thermostat. In units of K, eV or meV.

• T relaxation time: The time scale used to regulate temperature toward the target value in MD. Note that this is a relaxation time, and as such it gives a rough estimate of the time needed to achieve the given targets. As a normal simulation also exhibits oscillations, the actual time needed to reach the averaged targets will be significantly longer. In units of fs or ps.

NVT Berendsen

In NVT Berendsen simulations, the velocities are rescaled at each step so that the temperature relaxes exponentially toward its target value, at a rate set by T relaxation time.

• Fix center of mass: Check this box if you want the program to consider a fixed center of mass for the system.

• Logging frequency: defines how frequently the information of a running MD calculation is stored. During the MD simulation, atomic positions and velocities are stored. This is very useful to visualise and analyse the calculations.

• Temperature: Target temperature for the thermostat. In units of K, eV or meV.

• T relaxation time: The time scale used to regulate temperature toward the target value in MD. Note that this is a relaxation time, and as such it gives a rough estimate of the time needed to achieve the given targets. As a normal simulation also exhibits oscillations, the actual time needed to reach the averaged targets will be significantly longer. In units of fs or ps.

NVT Langevin

A stochastic thermostatting method. It controls system temperature by balancing a friction term (drag) with random fluctuating forces, rather than rescaling velocities directly.

• Fix center of mass: Check this box if you want the program to consider a fixed center of mass for the system.

• Logging frequency: defines how frequently the information of a running MD calculation is stored. During the MD simulation, atomic positions and velocities are stored. This is very useful to visualise and analyse the calculations.

• Temperature: Target temperature for the thermostat. In units of K, eV or meV.

• T relaxation time: The time scale used to regulate temperature toward the target value in MD. Note that this is a relaxation time, and as such it gives a rough estimate of the time needed to achieve the given targets. As a normal simulation also exhibits oscillations, the actual time needed to reach the averaged targets will be significantly longer. In units of fs or ps.

NVT Andersen

The constant temperature is imposed by stochastic collisions with a heat bath. Due to the random decorrelation of velocities, the dynamics are unphysical and cannot represent dynamical properties like e.g. diffusion or viscosity. Moreover, the collisions are stochastic in nature, so repeating the simulation will not give exactly the same trajectory.

• Fix center of mass: Check this box if you want the program to consider a fixed center of mass for the system.

• Logging frequency: defines how frequently the information of a running MD calculation is stored. During the MD simulation, atomic positions and velocities are stored. This is very useful to visualise and analyse the calculations.

• Temperature: Target temperature for the thermostat. In units of K, eV or meV.

• Collision probability: The collision probability defines the average number of collisions per atom and time-step.

NVT Bussi

The velocities are rescaled at each step by a stochastic factor that drives the temperature toward its target value while preserving the correct canonical fluctuations, at a rate set by T relaxation time.

Unlike the closely related Berendsen thermostat, the random component ensures that the temperature fluctuates as expected in a canonical ensemble, which makes this thermostat suitable for production runs.

In this MD-algorithm, if the atom velocities are set to zero, a warning will pop up asking whether you want to generate initial velocities according to the Maxwell-Boltzmann distribution at the target temperature.

• Logging frequency: defines how frequently the information of a running MD calculation is stored. During the MD simulation, atomic positions and velocities are stored. This is very useful to visualise and analyse the calculations.

• Temperature: Target temperature for the thermostat. In units of K, eV or meV.

• T relaxation time: The time scale used to regulate temperature toward the target value in MD. Note that this is a relaxation time, and as such it gives a rough estimate of the time needed to achieve the given targets. As a normal simulation also exhibits oscillations, the actual time needed to reach the averaged targets will be significantly longer. In units of fs or ps.

NPT Melchionna

An implementation of NPT dynamics combining a Nosé-Hoover thermostat with a Parrinello-Rahman barostat following Melchionna. The cell degrees of freedom are integrated together with the atomic coordinates, at rates set by T relaxation time and P relaxation time.

Unlike the Berendsen coupling, this algorithm reproduces the correct fluctuations of the isothermal-isobaric ensemble, which makes it the appropriate choice for production runs.

If your system does not have a cell when you try to select this algorithm the following warning will pop up:

If your system is a bulk but the simulation box is not un upper triangular matrix please adjust the simulation box:

• Logging frequency: defines how frequently the information of a running MD calculation is stored. During the MD simulation, atomic positions and velocities are stored. This is very useful to visualise and analyse the calculations.

• Temperature: Target temperature for the thermostat. In units of K, eV or meV.

• T relaxation time: The time scale used to regulate temperature toward the target value in MD. Note that this is a relaxation time, and as such it gives a rough estimate of the time needed to achieve the given targets. As a normal simulation also exhibits oscillations, the actual time needed to reach the averaged targets will be significantly longer. In units of fs or ps.

• Barostat mode: determines how the target pressure is applied. All modes share the following parameters:

  • Pressure units: Allow the user to select the units for the pressure, eV/Ångström\(^{3}\), GPa, kBar, or MPa.

  • P relaxation time: The time scale used to regulate pressure toward the target value in MD. Note that this is a relaxation parameter, it provides a rough time scale rather than the exact equilibration time. Because pressure fluctuates and exhibits oscillations, actual convergence to the target average requires significantly longer. Expressed in fs.

  • Bulk modulus: The material’s resistance to uniform compression, used to estimate volume response under pressure. Expressed in eV/Ångström\(^{3}\).

  • P-factor: A read-only parameter, calculated automatically from the P relaxation time and the Bulk modulus (\(\text{P-factor} = \text{P relaxation time}^2 \times \text{Bulk modulus}\)). It dictates how rapidly the simulation box volume adjusts to maintain target pressure. Typical metallic bulk moduli are on the order of 0.6 eV/Ångström\(^{3}\).

They modes differ only in which components of the pressure tensor are active and whether they share a single target pressure:

Hydrostatic pressure: Uniform pressure is locked on. The active components are the diagonal stress tensor terms \(xx\), \(yy\), and \(zz\) constrained to be equal, so a single target pressure value is required. Hydrostatic pressure means applying an isotropic stress tensor, which requires the diagonal components \(xx\), \(yy\), and \(zz\) to be equal, causing the corresponding cell axes to scale together by the same factor.

Hydrostatic pressure for selected axes: Uniform pressure is locked on. You choose which of the stress tensor diagonal terms to relax (for instance \(xx\), \(yy\) while leaving \(zz\) fixed). All active components share the same target pressure, so a single value is still required.

Orthorhombic cell (fixed angles): Uniform pressure can be unchecked. The active components are the diagonal terms. With Uniform pressure checked, the selected axes share one target pressure; with it unchecked, each enabled component is set independently.

Triclinic cell (fully flexible): Uniform pressure can be unchecked. In addition to the diagonal terms, the off-diagonal shear components are available. With Uniform pressure checked, all active components share one target pressure; with it unchecked, each enabled component is set independently.

NPT Berendsen

The velocities and the simulation box are rescaled at each step so that temperature and pressure relax exponentially toward their target values, at rates set by T relaxation time and P relaxation time. The coupling is efficient for equilibration, but it does not reproduce the correct fluctuations of the isothermal-isobaric ensemble; for production runs where statistical properties matter, use NPT Melchionna.

If your system does not have a cell when you try to select this algorithm the following warning will pop up:

• Fix center of mass: Check this box if you want the program to consider a fixed center of mass for the system.

• Logging frequency: defines how frequently the information of a running MD calculation is stored. During the MD simulation, atomic positions and velocities are stored. This is very useful to visualise and analyse the calculations.

• Temperature: Target temperature for the thermostat. In units of K, eV or meV.

• T relaxation time: The time scale used to regulate temperature toward the target value in MD. Note that this is a relaxation time, and as such it gives a rough estimate of the time needed to achieve the given targets. As a normal simulation also exhibits oscillations, the actual time needed to reach the averaged targets will be significantly longer. In units of fs or ps.

• Barostat mode: determines how the target pressure is applied. Two modes are available, Hydrostatic pressure and Orthorhombic cell (fixed angles). All modes share the following parameters:

  • Pressure units: Allow the user to select the units for the pressure, eV/Ångström\(^{3}\), GPa, kBar, or MPa.

  • P relaxation time: The time scale used to regulate pressure toward the target value in MD. Note that this is a relaxation parameter, it provides a rough time scale rather than the exact equilibration time. Because pressure fluctuates and exhibits oscillations, actual convergence to the target average requires significantly longer. Expressed in fs.

  • Bulk modulus: The material’s resistance to uniform compression, used to estimate volume response under pressure. Expressed in eV/Ångström\(^{3}\).

  • Compressibility: Read-only, computed from the bulk modulus. Expressed in GPa\(^{-1}\).

The two modes differ in which components of the pressure tensor are active:

• Hydrostatic pressure: Uniform pressure locked on. Active components \(xx\), \(yy\), and \(zz\), constrained equal: one target value.

• Orthorhombic cell (fixed angles): Diagonal terms only. Uniform pressure can be unchecked: checked, the selected axes share one target pressure; unchecked, each is set independently.

Variables (Velocities)

The Velocity Distribution widget is only available when selecting Molecular dynamics (multicode).

The Velocity Distribution widget shows the list of elements and the components of their velocity vectors.

The initial velocities are initially set to zero, however, they can be changed manually. The constraint on zero center of mass velocity is imposed automatically. Random velocities drawn by a Maxwell-Boltzmann distribution can be assigned to atoms of the system. You can set a desired temperature in the Maxwell-Boltzmann box. Press the Initialize now button to set up the initial velocities.

To maintain the system at the target temperature select the Force tick-box.

Press the OK button to close the parameter widget once the parameters are properly set.

Click on the Calculator icon to select the computational engine to be used during the molecular dynamics run.

We refer the user to chapter Calculators for further information on ASAP available calculators.

Click on the Run icon to open the Run widget. Then click on the Run button to submit the molecular dynamics calculation. See chapter Advanced Configuration and Remote Execution for further information on the computational resources configuration in ASAP.

After submitting a job (run), the Calculator output tab in Run widget shows the complete calculation output in real time.

In addition, if you are running a MD Algorithms (Multicode) project, the Task output tab in the Run widget shows relevant information on the Molecular Dynamic output in real-time.

Analysis

When the calculation is completed, select the Exit and analyse button to open the analysis widget.

The MD analysis widget is equivalent for Molecular dynamics (multicode) and Molecular dynamics (SIESTA only).

The MD series widget offers three plot types: Time series, Statistical properties and Radial distribution function.

At the top right corner of the widget, you can find the most important output parameters calculated during the MD run.

• \(\langle Ek \rangle\) and \(\langle T \rangle\) are the average kinetic energy and average temperature, respectively. These parameters are helpful to control the evolution of the system during dynamic run.

• \(\langle \Delta T^{2} \rangle\) represents the average temperature variation observed during an MD run. This parameter provides insight into the magnitude of temperature fluctuation that occur throughout the simulation. This measurement is particularly valuable for analyses conducted near phase transitions, as it effectively captures significant fluctuations that often manifest in such scenarios.

• 2\(\langle T \rangle^{2}\)/N is the averaged square temperature normalised by the degrees of freedom (N). It is an indirect measure of specific heat.

Check the Geometry tick-box, at the right side of the MD series widget, to visualise the geometry at any of the steps of the MD run.

The option View angles allows you to rotate the system around cartesian axes X, Y and Z by 90º.

Press the View in 3D… button to visualise the image with ASE GUI at the selected step.

Press the Play button to visualise the structure’s evolution during the MD run in three dimensions.

You can visualise the results of MD as time series (y-axis) of:

• kinetic energy

  

• kinetic and potential. For the visualisation of both kinetic and potential energies.

  

• potential energy

  

• pressure. The option is only available if provided by the calculator.

• temperature.

  

• total energy

  

Notice that you can modify the logging step range that will be shown in the visualisation,

and visualise the Time series as a function of steps or time (in unit of fs).

After editing a parameter, click the Update the plot button.

You have the option to visualize various statistical measures, explained below, over a range of steps or time (in fs units) by selecting ’Statistical Properties’ from the ’Plot Types’ menu.

Note that it is possible to select certain atoms from the system with the Range of atom indices feature. You can select the particles one by one, e.g. (0,1,2,3…) or in batches, e.g. (0-4).

Additionally, you can choose which chemical species information to display by clicking them in the Chemical species feature.

• Mean-square displacement: The average square distance taken over all particles. Gives the deviation of the position of a particle with respect to a reference position over a certain period of time (units of Ångström\(^{2}\)).

  

• Diffusion coefficient: obtained from the slope of the mean-square displacement versus time (Einstein relation) (units of Ångström\(^{2}\)/fs).

  

• Root-mean-square displacement: measures the variability of the positions of the particles with respect to a reference position over a certain period of time. It is particularly useful to compare two data sets of a system that has been transformed in some way (units of Ångström).

  

• Velocity correlation function (VCF): average of the product of the velocities of all the particles at two different times (t, t+\(\Delta\)t). This function quantifies the rate at which the particles return to equilibrium state after a small perturbation.

  

• Power spectrum: allows the user to visualize the frequency distribution of the velocity correlation function (units of fs). It helps to identify characteristic frequencies and modes of vibration of the system. Note that for this statistical property the options for the abscissa axis are not steps or (explicitly) time, but frequency as a function of time (units of fs) or inverse of time (Hz), inverse of wavelength (units of cm^-1) and energy (units of eV, meV, and Rydberg and Hartree units of energy).

  

Click the Subtract center of mass check-box to dismiss the contribution of the center of mass of the system (option available for Mean-square displacement, Diffusion coefficient and Root mean-square displacement).

The Radial distribution plot type allows you to visualize the Radial Distribution Function (RDF). RDF measures the probability of finding an atom at a distance r (units of Ångström) given that there is an atom at position r=0. Essentially, it is a histogram of interatomic distances.

You can tune the following parameters:

• Num. intervals: parameter defining the sampling of the x axis.

• Snapshot step: to select a specific moment in time of the MD run.

• Min radius: minimum radial distance for the RDF.

• Max radius: maximum radial distance for the RDF.

• Units: optionally units of Ångström, nanometers or Bohr radius.

• Plot style: optionally lines or bars.