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 relaxation include are CG, Broyden, and FIRE (some of which are covered in Structural optimization using forces and stresses). We will explore a new variant to relax structures in this tutorial.

  • 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.

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. 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, 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, it is not seamless (but can be used for restarting a simulation from a different ensemble).

Run the simulation with the provided pseudopotential and .fdf file. It should hopefully not take very long. Now look at the resulting files in the directory: the .VERLET_RESTART and .XV files have been produced, so re-running 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. 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. A good estimation can be provided by

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

Where 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).

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 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?

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.

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?

Geometry Relaxations

Note

Reviewing Structural optimization using forces and stresses before running these might be worth the effort!

We will see now two examples of geometry optimization. While the first uses the traditional approach of relaxing the system, the second tries a different approach, relaxing via MD by cooling down the system.

Hint

Go to the directory 04a-Relax-CG.

Take a look at the .fdf file, and note the differences from earlier .fdf files:

MD.TypeOfRun          CG
MD.NumCGSteps         100
MD.VariableCell       T
MD.MaxForceTol        0.1 eV/Ang
MD.MaxStressTol       0.1 GPa
MD.TargetPressure     0.0 GPa

Notably, the MD.TypeOfRun is now CG, which does a coordinate optimization using the Conjugate Gradients method. By specifying MD.VariableCell T, we allow the lattice vectors to also vary and be relaxed through CG along with the atom coordinates. MaxForceTol sets the atomic force tolerance, MaxStressTol sets the cell stress tolerance (since this is variable cell CG). The optimization process will stop once the maximum force and the maximum stress fall below these tolerances.

Note

You can visualize the geometry optimization in the same way as the Parrinello-Rahman dynamics.

Hint

Go to the directory 04b-Relax-PR.

Now look take a look at this .fdf file:

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
MD.Quench                T

This file is largely similar to the Parrinello-Rahman dynamics from earlier, but with one key difference: MD.Quench is toggled on. This means that we will be doing a Quenched MD, in which we slowly cool down the system in order to relax it.

Run the example, and 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.

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

Annealing - Si(001) + H2

Hint

Go to the directory 05-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 may 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?