Molecular Dynamics


Alec Wills (Stony Brook University)

This exercise will take you through various scenarios for the different choices of dynamics one can run in SIESTA.

At the time of writing, the Read the Docs manual links to a PDF verson of the user’s guide for version MaX-1.3.0 of SIESTA. Section 7 of the manual contains the options which this tutorial will give an overview of.

Most of the input options are prepended with MD.X, for some option X.

The MD.TypeOfRun option has you choose what kind of dynamics you wish to run. The coordinate relaxation options are CG, Broyden, and FIRE. Options for standard time-integrated molecular dynamics can be chosen through Verlet (NVE), Nose (NVT), ParrinelloRahman (NPE), and NoseParrinelloRahman (NPT).

Examples for running dynamics simulations are described below. The exercises utilize a Silicon supercell in the diamond structure. For efficiency on the virtual machine, the basis size is only a single-\(\zeta\) basis with short orbitals.

Verlet Dynamics

Look at the tutorial folder verlet, and examine the .fdf file within.

Note the MD options selected – the Verlet run will run for 200 steps with a timestep of 3.0 fs. The initial velocities are selected to correspond to a temperature of 300 K. Per the documentation:

The atoms are assigned random velocities drawn from the Maxwell-Bolzmann distribution with the corresponding temperature. The constraint of zero center of mass velocity is imposed.

Also note the WriteMDHistory option, which is toggled to true:

If true, SIESTA accumulates the molecular dynamics trajectory in the following files:

  • SystemLabel.MD : atomic coordinates and velocities (and lattice vectors and their time derivatives, if the dynamics implies variable cell). The information is stored unformatted for postprocessing with utility programs to analyze the MD trajectory.

  • SystemLabel.MDE: shorter description of the run, with energy, temperature, etc., per timestep.

As noted just beneath the quoted block in the manual, these files are cumulative, so you can always restart the simulation for more sampling if the appropriate files are in a given directory. The files needed for a seamless restart are the SystemLabel.X_RESTART file (where X is some label corresponding to chosen dynamics option) and the SystemLabel.XV file, which stores current velocities and predicted positions for the next step. One can restart with just the .XV file as well, it is not seamless (but can be used for restarting a simulation from a different ensemble).

Run the simulation with the provied pseudopotential and .fdf file. It should hopefully not take very long (on my laptop, it took 54 seconds). Now look at the resulting files in the directory – the .VERLET_RESTART and .XV files have been produced, so rerunning the siesta command will add more configurations and thermodynamic values to the trajectory and .MDE files, respectively.

Plot the total and Kohn-Sham energies from the si8.MDE file as a function of timestep. xmgrace is pre-installed on the virtual machine. One can use xmgrace -block DATA -bxy x:y to plot, where DATA is the columnar data file (i.e. the .MDE file), x is the index of the x-data column, and similar for y. Is energy conserved? How large are the fluctuations?

Now look at the temperature as a function of time. How does it change as the system evolves? Look also at the pressure, to see what happens to it.

Finally, examine the trajectory itself – the virtual machine comes with jmol pre-installed. Run this on the .ANI file (jmol si8.ANI – the program should identify the format as an .xyz format) and animate the atoms’ trajectory. (Tools–> Animate –> Loop).

Now try re-running the system with a larger timestep (say, 10 fs). How does the system respond? Make sure to change the SystemLabel before re-running this time, so you can compare plots from separate files.

Nose Dynamics

Now change to any of the nose-X-Y directories (where X is target temperature, Y is Nose mass), and examine again the .fdf file included. Note now that there MD.TypeOfRun is Nose, and the parameters MD.TargetTemperature and MD.NoseMass are specified. The TargetTemperature is self-explanatory, but the NoseMass affects the time scale of the temperature variation in the system now.

Run each of the different nose-X-Y simulations. Once complete, plot the temperature as a function of timestep to compare the effect the Nose mass has on a given target temperature. Try to plot the running average of the temperature over a given number of timesteps (i.e. 25). How does the average fluctuate across Nose masses?

Now pick a simulation with a certain Nose mass to increase the timestep of (making sure to change the SystemLabel). Compare the results for different timesteps with the same Nose mass. How do the temperature values and the averages differ?

Parrinello-Rahman Dynamics

Now change to the PR directory, and examine the .fdf file. Now the MD.TypeOfRun is ParrinelloRahman, with corresponding options MD.ParrinelloRahmanMass and MD.TargetPressure. The settings set the desired system pressure to 0. Run the file and examine the resulting pressures and temperatures. Is zero pressure achieved? Check a running window average (with, say, 25 timesteps). How does the average look? Check the temperature and its average as well.

Now we want to see the behavior of the cell. Run md2axsf in the PR directory. Specify your SystemLabel if not the default. To examine the behavior of the unit cell, have it analyze the .MD file and output all the steps. Visualize the behavior of the atoms and the cell by running jmol on the resulting .AXSF file. Look at the steps where the cell seems particularly compressed, and compare it with the plot of pressure over time. Do the values make sense? What about steps where the cell is particularly expanded?


There are two example directories for relaxation: Relax-CG and Relax-PR.

First look at the Relax-CG directory’s .fdf file, and note the differences from earlier .fdf files. Notably, the MD.TypeOfRun is now CG, which does coordinate optimization. By specifying MD.VariableCell, we allow the lattice vectors to also vary and be relaxed through CG along with the atom coordinates. Run this simulation and visualize the relaxation process as we did in the Parrinello-Rahman dynamics. Note that the MaxForceTol sets the atomic force tolerance, MaxStressTol sets the cell stress tolerance (since this is variable cell CG).

Now look at the Relax-PR directory’s .fdf file. This file is largely similar to the PR dynamics directory from earlier, but with one key difference – MD.Quench is toggled. This option is not documented, but examine the .MDE file’s temperature as time goes on. Since we start the dynamics without specifying a temperature, the cell contraction will cause the stationary atoms to start moving, which is why the temperature spikes. The MD.Quench option essentially cools the atoms back down so that they remain within the relaxed box as the process continues.

Which method takes more iterations to achieve relaxation? How closely do the lattice parameters agree?

Si(001) + H2

Look at the .fdf in the final directory, H2_md. The system is now composed of Si atoms and H atoms, and the H atoms have been given the mass of deuterium to allow for a faster timestep. Note the MD options specified: dynamics will be simulated by Anneal-ing to the desired temperature with the specified relaxation time TauRelax.

Note that the .fdf file specifies to UseSaveXV to start a simulation. For this case, it functions as an overwrite for the velocity distribution – rather than using the InitialTemperature (which defaults to 0 K) to generate velocities, a file is manually provided to the program to zero all velocities except for an H2 dimer, which is made to approach the Si(001) surface. In this manner, we can control the behavior of atoms or molecules in the system, forcing certain reactions to occur. The corresponding files are provided in the FILES directory, so be sure to copy them into the working directory before running the simulation. (Note: you should copy these files before each subsequent run, as they’re overwritten by the output). Look at the .XV file. Which atoms are moving, and where are they moving towards?

Also note the %block of GeometryConstraints. This is Section 7.7 of the SIESTA manual, and you can do many things with this block. Specifically, in this simulation position from 19 to 38 fixes the coordinates of the back half of the Si structure, so that only the surface exposed to the H2 molecule can move and restructure itself after any reactions.

Visualize the atom trajectories using jmol on the .ANI file. What behavior occurs? Look at the .MDE file now – is the desired temperature reached? How does the energy behave as the simulation progresses?

Now change the .fdf file to do simulate the NVE ensemble, making sure to change the SystemLabel. Copy the relevant .XV file and rename it to the appropriate filename, and redo the simulation. Visualize the behavior – how does it change? Is the energy conserved this time? Why might the behavior you see be happening?