Molecular Dynamics

SIESTA is capable of running molecular dynamics simulations under different ensembles (NVT, NVE, NPT). This tutorial will show you around the different possibilities for SIESTA MD.

Typically, most of the input options for the fdf are prepended with MD.*, for example MD.TypeOfRun. This very option has you choose what kind of dynamics you wish to run:

  • Coordinate relaxations include are CG, Broyden, and FIRE (which are covered in Structural optimization using forces and stress).

  • Options for regular molecular dynamics can be chosen through Verlet (NVE), Nose (NVT), ParrinelloRahman (NPE), and NoseParrinelloRahman (NPT). These are new options that will be explored in this tutorial.

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

Note

For molecular dynamics, the file *.ANI contains all of the structures generated during the relaxation in XYZ format. It can be processed by various graphical tools, such as Ovito, XCrysden or VMD.

If you want to plot energy, pressure and temperature across the molecular dynamics run, you can use the scripts available in plotscripts under this same tutorial folder.

Time Integrated Molecular Dynamics (MD)

We will cover here three types of molecular dynamics: Verlet (NVE), Nose (NVT) and Parrinello-Rahman (NPE). The last option, Nose-Parrinello-Rahman (NPT), is just the combination of the NPE and NVT options.

Verlet Dynamics (NVE)

Hint

Go to the directory 01-Verlet.

Take a close look at the .fdf file within, in particular this section:

MD.TypeOfRun          Verlet
MD.InitialTimeStep    1
MD.FinalTimeStep      200
MD.LengthTimeStep     3.0 fs
MD.InitialTemperature 300.0 K

WriteMDHistory T

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; in essence, this means that the atoms are assigned random velocities drawn from the Maxwell-Bolzmann distribution for that temperature. A constraint of zero center of mass velocity is imposed.

Also note the WriteMDHistory option, which is toggled to true. This means that 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.

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, such as VERLET or NOSE) 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, but it will be essentially re-setting the information related to thermostats, barostats and velocity propagators.

Running the simulations

Run the simulation with the provided pseudopotential and .fdf file. It should not take very long. Now look at the resulting files in the directory: the .VERLET_RESTART and .XV files have been produced, so if you re-run SIESTA, it will add more configurations and information to the trajectory and .MDE files, respectively.

Plot the total and Kohn-Sham energies from the si8.MDE file as a function of timestep. 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.

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 (or back up your previous results), so you can compare plots from separate files.

Nose Dynamics (NVT)

Hint

Go to the directory 02a-Nose-300-100.

Now we move on to constant volume, constant temperature dynamics. We will use the Nose thermostat, which provides a thermal bath coupled to the system. The magnitude of this coupling is controlled by the MD.NoseMass variable, as you can find in the fdf:

MD.TypeOfRun            Nose
MD.InitialTimeStep      1
MD.FinalTimeStep        200
MD.LengthTimeStep       3.0 fs
MD.InitialTemperature   300.0 K
MD.TargetTemperature    300.0 K
MD.NoseMass             100.0 Ry*fs**2

See that we also have a new option here, MD.TargetTemperature, which, as the name implies, sets the desired (average) temperature for the system. In a real production run, MD.NoseMass should be carefully chosen by testing how the system responds to it. This reference can be helpful in determining a reasonable value for it:

\(Q = g * K * T * t^{2} / (2 \pi^{2})\)

Here, g is 3 times the number of atoms (i.e. the degrees of freedom for the system), T is the target temperature, K is the Boltzmann constant, and t is the timescale of the dominant modes for the system itself (usually 1-10 fs for lighter atoms, 10-100 fs por larger atoms).

Running the simulations

Run the system, and see how the temperature in the .MDE file behaves.

Hint

Go to the directory 02b-Nose-300-10 and later to the directory 02c-Nose-300-1.

This is the same system as before, but changing the NoseMass to 10 Ry*fs**2 and 1 Ry*fs**2 respectively. Run these examples.

Once the simulations are completed, plot the temperature as a function of time-step 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?

Hint

Now do the same for 02d-Nose-600-10.

Now pick a simulation with a certain Nose mass to increase the timestep itself. (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 (NPE)

Hint

Go to the directory 03-ParrinelloRahman.

Now, instead of maintaining the constant temperature, we will keep constant pressure. This means that the cell sizes will change during the MD. Examine the .fdf file here. Now the MD.TypeOfRun is ParrinelloRahman, with corresponding options MD.ParrinelloRahmanMass and MD.TargetPressure:

MD.TypeOfRun          ParrinelloRahman
MD.InitialTimeStep    1
MD.FinalTimeStep      200
MD.LengthTimeStep           3.0 fs
MD.ParrinelloRahmanMass  10.0 Ry*fs**2
MD.TargetPressure     0.0 GPa

Here the MD.ParrinelloRahmanMass plays a similar role to the mass of the Nose thermostat. The settings set the desired system pressure to 0.

Running the simulations

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.

Note

Having a “negative” instant pressure might be a bit confusing here. What this really means in MD is that we are taking into account the direction in which the force is exerted (i.e., towards inside or outside the simulation box).

Optional - visualization

Now we want to see the behavior of the cell, but for that we will need to change our files to aformat readable by other visualization programs. Run for example md2axsf in the PR directory. Specify your SystemLabel to si8 (or the one you are using if you changed its value). 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 your favorite visualization program (jmol, ovito, VMD, xcrysden, etc) 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?

Annealing - Si(001) + H2 and Constraints

Hint

Go to the directory 04-AnnealSiH2.

We will cover now one last option for MD, which works specially well for system equilibration to a desired temperature and/or pressure. Look at the .fdf file in this final example. 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:

MD.TypeofRun          Anneal
MD.InitialTimeStep    1
MD.FinalTimeStep      100
MD.LengthTimeStep     2.0 fs
MD.AnnealOption       Temperature
MD.TargetTemperature  300.0 K
MD.TauRelax           50.0 fs

This time, the dynamics will be simulated by annealing to the desired temperature, according to the specified relaxation time TauRelax. Also, note the MD.AnnealOption variable; it can take the values of Temperature, Pressure and TemperatureAndPressure. The annealing method slowly takes the system from the initial to the desired conditions, using TauRelax as the characteristic time for the process. Beware that annealing is only recommended to equilibrate a system before moving on to other methods; taking MD-related properties (such as radial distribution functions) is not recommended since annealing will yield unphysical results.

Independently from this, you will also find the following option in the fdf file:

MD.UseSaveXV   .true.

This means that SIESTA will attempt to find an .XV file to use as a restart, taking the initial coordinates and velocities from it. This will superceed the coordinates present in the fdf file and skip the initial random velocity generation (this means that MD.InitialTemperature will not have any effect if an .XV is found).

Hint

You can find a restart file called Si001+H2.XV under restart_files. Copy it to your working directory every time you run a new simulation, so that you are always guaranteed to start from the same conditions.

Take a look at the Si001+H2.XV file. Which atoms are moving, and where are they moving towards?

There is one last new thing in this fdf:

%block GeometryConstraints
      position  from   19  to  38
%endblock GeometryConstraints

This basically fixes the coordinates of the bottom half of the Si structure, so that only the surface interacting with H2 molecule is allowed to move, accelerating the calculations. This block is very flexible and other variants can be found in the manual.

Run the example and take a 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, using the Verlet algorithm as before (suggestion: change the name of SystemLabel to keep runs separate). Copy the Si001+H2.XV restart file and rename it to the appropriate filename; them, redo the simulation. Take another look at the .MDE file: How does it change? Is the energy conserved this time?

Relaxation with Quenched Molecular Dynamics advanced

There is an alternative relaxation method that uses a physically motivated scheme, rather than a purely mathematical search for a ‘zero forces and stress’ configuration. Imagine that we perform a molecular dynamics simulation in which, rather than relaxing, we move the atoms, and the cell vectors, according to the (classical) equations of motion, using the forces and stress (see this tutorial). The trick in this case is that, every time an atom senses a force ‘opposite’ its velocity (in the sense that their scalar product is negative), the velocity is set to zero. This roughly corresponds to the idea: “since the atom seems to be moving away from its equilibrium point, we rather stop it”. The same can be done with strains and stress in the case of variable cell.

Hint

Enter the directory 05-quenched_md

Look now at the “relaxation section” of the si8.fdf file:

MD.TypeOfRun             ParrinelloRahman
MD.InitialTimeStep       1
MD.FinalTimeStep         200
MD.LengthTimeStep        3.0 fs
MD.ParrinelloRahmanMass  10.0 Ry*fs**2
MD.Quench                T

The ParrinelloRahman scheme is a combined atoms+cell microcanonical scheme (refer to the MD tutorial). We allow it to run for up to 200 steps, with a time-step of three femtoseconds. Note also the appearance of a “mass” with the dimensions of energy*time^2: this is used to homogeneize the dynamics of the system, which has to deal, as we indicated earlier, with fundamentally different sets of variables. The final line requests the quenching of the molecular dynamics, setting the forces to zero in the direction opposite to that of the atom’s velocity.

If you now run the example, you will notice that it converges quite nicely, with a monotonic decrease in the energy and a cleaner evolution of the stress (even if not entirely monotonic). This method is always quite robust, and particularly efficient in the case of variable-cell. It can be generally counted on to bring systems closer to the optimal structure, and the final relaxation can be done afterwards with a faster method. We can even start with some extra kinetic energy, in the form of some ‘starting temperature’, to wiggle things around and free the system from any undesired local minimum.