Molecular dynamics: Basics


Virtual NanoLab and ATK make molecular dynamics 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. This tutorial provides a basic overview of the MD functionalities in VNL-ATK and explains step-by-step how to set up a simulation correctly to obtain the desired results.

Molecular dynamics (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.


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. 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 to propagate the positions and velocities 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 this tutorial, you will get to know the basic ingredients of MD simulations and how you can apply these to run MD simulations with the ATK Molecular Dynamics module.

You will apply these techniques to bulk silicon as a simple and well-known test system. You will learn step-by-step how to select the time step size, how to control the temperature and the pressure, and how to fix the position individual atoms by setting constraints.

NVE Simulations

As mentioned in the introduction, 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.

In this section, you will gain experience with running NVE simulations, and how to assess their quality based on monitoring the total energy (which should be conserved). You will learn how to choose a suitable time step size for a given system.

Setting up the geometry

Start VNL, create a new project with a suitable name, and open it. Launch the Builder by clicking the builder_icon icon on the toolbar.

In the builder, click Add ‣ From Database. Search for “Silicon” and add it to the stash.

The primitive unit cell of silicon has 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). You can obtain an orthogonal silicon cell by using the Bulk Tools ‣ Supercell tool in the builder panel. Click Conventional and Transform to convert the primitive unit cell into a conventional, cubic cell.


Furthermore, it is useful to increase the size of the silicon 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). Increase the size of the structure by clicking Bulk Tools ‣ Repeat and choose 4 x 4 x 4 repetitions. This will result in a total number of 512 atoms.

Next send the structure to the Script Generator by using the sendto_icon icon in the lower right-hand corner of the builder window.

Setting up the calculation script

In the following, you will set up the calculation script that will be used to run the MD calculations. We will use a Tersoff-type potential for the Si-Si interactions [1].

In the Script Generator do the following

  • Add a New Calculator.
  • Add an Optimize > MolecularDynamics block.

Double-click the New Calculator and select the ATK-ForceField calculator with the Tersoff_Si_1988 potential [1]. Uncheck Save and Print, as we do not need to save or print additional informations about the original configuration.


For ATK-versions older than 2017, the ATK-ForceField calculator can be found under the name ATK-Classical.

Now, double-click the MolecularDynamics block. First of all, take a look at the Type combobox. This is the most important setting as it allows you to select the integration algorithm and the physical ensemble the system is simulated in. For the first simulation, set the type to NVE Velocity Verlet. This setting selects an ensemble with constant number of particles (N), constant volume (V), and constant total energy (E), and uses the velocity verlet integration algorithm [2].

Set the number of Steps to 1000 and the Log Interval to 5. The latter value 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.


In most MD applications it is not necessary to record snapshots at such a high frequency unless you explicitly want to analyze high-frequency oscillations. If you suspect that saving the snapshots to disc unduly slows down your simulation, you can uncheck the Save Trajectory checkbox to avoid writing to disc during the simulation. Instead you can use the Save option at the bottom, to write the trajectory only at the end of the simulation. For technical reasons, this is substantially faster but preserves the same amount of information. On the downside, you cannot use the intermediate snapshots to restart from the last snapshot, in case the simulation is interrupted.

Enter a suitable trajectory file name, e.g. Set Initial Velocity type to Maxwell-Boltzmann at a Temperature of 600 K. This will create the random initial velocities of the atoms drawn from a Maxwell-Boltzmann distribution corresponding to a temperature of 600 K. Additionally, as selected by the corresponding check-box, the momentum of the center of mass of the entire cell is removed from the distribution to avoid drifting of the system during the simulation.

The Time step is a crucial parameter in MD simulations as it determines the accuracy and efficiency of the numerical integration scheme. We will investigate its effect in detail in the following, but, to begin with, we keep the default value of 1 fs.

Finally, you can uncheck the Save and Print check boxes, as the trajectory is already being saved.

In the end, your MD settings should look like the following.


Click OK and send the script to the Job Manager (again using the sendto_icon button). Save the script as and start the simulation. It should take less than a minute to complete.

Analyzing the results

Using the Movie Tool plugin, you can investigate the temperature as well as the kinetic, potential, and total energies of the system as a function of simulation time. Simultaneously, an animation of the motion of the atoms is displayed. Select the icon from the LabFloor and click the Movie Tool plugin.

If you take a look at the temperature, you will observe how the temperature drops immediately and then oscillates around 300 K. The reason is that the displayed temperature is the instantaneous temperature of the system, calculated from the kinetic energy of the atoms, via \(E_{kin}(t) = 3/2 N k_B T(t)\). Since the atoms are not free but interact with the neighboring atoms in the lattice, their initial kinetic energy will be partially converted into potential energy. As the initial configuration is close to the potential energy minimum approximately half of the kinetic energy is converted, resulting in a temperature drop of around 50 %. In fact, if you take a close look at the potential energy curve, you will see an opposite picture (potential energy quickly rising, then oscillating around a stationary value).

Since you simulated an NVE ensemble, a crucial parameter to monitor is the conservation of the total energy. With the zoom tool in the toolbar below the energy graph, zoom into the total energy curve. You should see a curve similar (but probably not exacly identical) to the one shown below.


Overall, the energy of the system is conserved up to very small fluctuations with an order of magnitude of less than \(10^{-5}\) relative to the average value. Importantly, the average value is stationary and no drift can be observed. These findings indicate properly chosen MD settings.

To see what happens if you choose a time step size that is too large, repeat the above simulation but with an initial temperature of 5000 K and a time step size of 5 fs. You can easily change these parameters by re-using your MolecularDynamics block or manually editing the script by dropping the script onto the Editor icon. Remember to also change the name of the trajectory files. Save this new, modified script to the file and run the job.

If you visualize the resulting simulation trajectory in the Movie Tool, you will find an increased motion of the atoms. The structure will even partially melt. If you again zoom into the total energy curve, you will observe a pronouced drift of this value. The energy is not conserved very well any more, as the timestep has been chosen too large for this temperature.


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, if you do not know what time step to use, is 1 fs. Larger timestep 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. In this example, it takes around 600 fs before the system temperature as well as the magnitude of the total energy total fluctuations have reached a stationary state. 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 the present case you should then wait at least 600 fs. 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

In the following section you will learn how to perform MD simulations in the NVT ensemble, i.e. simulations in which the temperature is controlled. 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.

You will use the same bulk silicon system with the ATK-ForceField calculator and the Tersoff-potential, as you did in the previous section.

Setting up the script

You will start with an equilibration of the structure at room temperature. Prepare the script in the same way and with the same blocks (NewCalculator and Optimize ‣ MolecularDynamics) as described above. You can, in fact, re-use your previous script if still stored in the ScriptGenerator.

Double-click on the MolecularDynamics block to open its settings. As we intend to run a simulation in the NVT ensemble, we need to select a different Type this time. ATK offers four possibilities for NVT simulations:

  1. NVT Berendsen [3]
  2. NVT Nose Hoover [4]
  3. Langevin [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

In the MolecularDynamics settings, select the type NVT Nose Hoover. Set the number of steps to 5000 and the log interval to 20. Save the trajectory to a file Select Maxwell-Boltzmann initial velocities at 300 K. Leave the time step at 1 fs, the Reservoir temperature (which specifies the temperature of the virtual heat bath) at 300 K, and the Thermostat timescale at 100 fs. The latter value determines how quickly the system temperature approaches the reservoir temperature. We will take a closer look at the effect of this parameter later on in this section. The Final temperature value will always be synchronized with the reservoir temperature, unless you explicitly specifiy a different value here. This parameter can be used to linearly increase or decrease the temperature during the temperature (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.


Click OK to close the molecular dynamics settings and run the script by sending it to the JobManager. If you used the Script Generator set up for the the previous calculation, remember to rename the script file and save it to

After the simulation has finished, the trajectory file will appear on the Lab floor. Open it with the Movie Tool.


Take a look at the temperature curve. You will find that after a short initial equilibration period the curve fluctuates slightly around the selected reservoir temperature of 300 K. In fact, it does not look very different from the results of the NVE simulation in the previous section. In contrast to the NVE simulation, however, the final equilibrium temperature of the system is independent of the initial value, set via the Maxwell-Boltzmann distribution, and is instead coupled to the reservoir temperature. Let us take a closer look at this now.

Open again the Script Generator of the previous simulation. Double-click the MolecularDynamics block and change the Reservoir temperature to 1000 K. Leave the initial velocities as a Maxwell-Boltzmann distribution at 300 K. Finally, change the trajectory filename to Close the settings and save and run the script.

The result should look similar to the following picture.


The temperature curve begins at an inital value of 300 K (as selected via the Maxwell-Boltzmann distribution). Subsequently, as the system is in contact with the heat bath, the temperature quickly increases towards the reservoir temperature of 1000 K.

You may also notice an increase in the total energy. In contrast to NVE simulations, this is not a sign of poorly chosen parameters, but reflects the fact that energy is exchanged with the virtual heat bath.

During the simulation the motion of the atoms increases compared to the previous simulation. If you increase the reservoir temperature to even higher values (e.g 3000 K), you will find that the crystal will begin to melt and become partially amorphous. Annealing crystals at high temperatures is in fact a common way to generate amorphous structures of a material. You can learn more about this method in the tutorial Generating Amorphous Structures.

As mentioned, how fast the temperature of the system approaches the target temperature can be controlled via the thermostat timescale parameter. Now you will repeat the same simulation with a different thermostat timescale value. To this aim, start from the same script as in the previous simulation and open the MolecularDynamics settings. Leave the settings at the same values but change the Thermostat timescale value to 500 fs. To fully capture the prolonged equilibration phase increase the Number of steps to 10 000 and change the filename to Close the settings and run the simulation.

The results of this simulation should look similar to the following picture.


Qualitatively, the results are very similar to the previous simulation. You will notice, however, that the approach of the temperature towards the reservoir temperature proceeds much slower now.

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 and Langevin.

The Berendsen thermostat implements an algorithm which is in some situations more stable, as it effectively supresses 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 [5], 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.

NPT Simulations

In this section you will learn how to perform constant pressure simulations.

Setting up the calculation

Set up a script containing the same blocks as for the NVT simulation (i.e. BulkConfiguration, NewCalculator, Optimize ‣ MolecularDynamics), or reuse the script from the previous simulation if you did not close it.

Open the MolecularDynamics settings. Under Type, you will find two possibilities to run NPT simulations: The NPT Berendsen and the NPT Martyna Tobias Klein algorithm.

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 production simulations, you are recommended to use the NPT Martyna Tobias Klein barostat [6] instead. Suitable settings to equilibrate the silicon supercell under constant temperature and constant pressure of 1 bar are shown in the picture below.


The Molecular Dynamics group box contains the general settings, which should be well-known by now. 10 000 integration steps will be sufficient to demonstrate the barostat. The log interval can be set to 100 steps which should be enough to resolve the fluctuations. Be aware that this produces rather large file sizes, and you may consider removing the trajectories after the evaluation.

The next group box contains the thermostat and barostat settings. 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. The values here can be left at default values for this simulation.

Finally, 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 ATK Reference Manual and the tutorial Simulating a creep experiment of polycrystalline copper provide a more detailed explanation of this functionality.

In this example, use the default option with Isotropic pressure checked and set to 1 bar.

Editing the script

Close the settings window and send the script to the Editor to add some extra functionality. To visualize the time dependent volume and pressure fluctuations, append the following lines of code at the end of the script:

# Get the number of snapshots in the trajectory
trajectory_length = md_trajectory.length()

# Extract the cell volume from the trajectory.
volume = [md_trajectory.image(i).bravaisLattice().unitCellVolume().inUnitsOf(Ang**3) for i in range(trajectory_length)]

# Extract the pressures
pressures = md_trajectory.pressures().inUnitsOf(bar)

# Extract the simulation time
times = md_trajectory.times().inUnitsOf(ps)

import pylab

# Plot the volume over time.
pylab.subplot(2, 1, 1)
pylab.plot(times, volume)
pylab.xlabel('Time (ps)')
pylab.ylabel('Volume (Ang**3)')

# Plot the pressure over time.
pylab.subplot(2, 1, 2)
pylab.plot(times, pressures)
pylab.xlabel('Time (ps)')
pylab.ylabel('Pressure (bar)')

These added lines will produce a plot of volume vs. time by extracting the cell volume as well as the pressure from each snapshot saved in the trajectory and plotting the values over the simulation times.

Simulation results

Save the script to and run the simulation by sending the script to the Job Manager. It will take a few minutes to complete. Opening the trajectory with the Movie Tool, you can monitor the changes in the size and shape of the cell during the simulation. The temperature fluctuates around the selected reservoir temperature of 300 K.

To inspect the volume fluctuations you can open the file Fluctuations_NPT.png. It should look similar to the following figure.


The cell volume in the upper plot shows pronounced oscillations. These oscillations are a typical outcome of the Nose-Hoover-like temperature and pressure control. As long as their magnitude stays within reasonable boundaries, and does not increase dramatically throughout the simulation, there is no need to worry. The period of these oscillations can be changed by the Barostat timescale setting. The smaller this value is chosen, the higher the frequency. The amplitude of the oscillations increases with decreasing barostat timescale.

The pressure during the simulation, as shown in the lower plot, shows a similar behaviour. It oscillates around the target pressure of 1 bar, but the amplitude is large. Such large fluctuations of the instantaneous pressure values are in fact not unusual in MD simulations, as the pressure is calculated based on the stress, which varies considerably with small changes of the volume. Moreover, pressure is a macrosopic property, and for such a small system it is more reasonable to consider the corresponding average rather than the instantaneous values.

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. [7].

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 open the Molecular Dynamics widget in the Script Generator and then click on Add constraints to open the Constraints editor widget.


In the widget you will find a table listing all atom groups that are marked by tags. You can add new groups by selecting the desired atoms in the presenter window (draw a rectangle or hold Ctrl to select more than one atom) and clicking Add tag from selection. Constraints can be assigned to each of these atom groups in the Constraint column. 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 constraints their movement in the X, Y, or Z direction by selecting Fixed X, Fixed Y, or Fixed Z, respectively. The Rigid constraint type 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, you have to take care that at least one atom remains 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 Movie Tool or the MD Analyzer are not designed to work with constraints. For instance, the temperature displayed in the Movie Tool will not reflect the actual temperature of the free atoms in the system, but rather indicate a lower temperature value when constraints are present during the MD simulation because it also includes the atoms with zero temperature (the fixed atoms). If you want to monitor the temperature, you can do this in a script via the kinetic energy, as demonstrated in the following example:

# Get the kinetic energies
kinetic_energies = md_trajectory.kineticEnergies()

# The number of constrained atoms
number_of_constraints = 10

# Calculate the degrees of freedom in the system
degrees_of_freedom = 3 * (bulk_configuration.numberOfAtoms() - number_of_constraints)

# Calculate the temperatures from the kinetic energies
temperatures = (kinetic_energies / (0.5 * degrees_of_freedom) / boltzmann_constant).convertTo(Kelvin)


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. Alternative approaches are available to optimize the cell of a DeviceConfiguration (see the tutorial Geometry optimization of interfaces)


[1](1, 2) J. Tersoff: New empirical approach for the structure and energy of covalent systems. Phys. Rev. B 37, 6991 (1988)
[2]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)
[3]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).
[4]G. J. Martyna, M. L. Klein, M. Tuckerman: Nosé–Hoover chains: The canonical ensemble via continuous dynamics. J. Chem. Phys. 97, 2635 (1992)
[5](1, 2) W. F. van Gunsteren, H. J. C. Berendsen: A leap-frog algorithm for stochastic dynamics. Mol. Sim. 1, 173 (1988)
[6]G. J. Martyna, D. J. Tobias, and M. L. Klein: Constant pressure molecular dynamics algorithms. J. Chem. Phys. 101, 4177 (1994)
[7]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)