Molecular Dynamics¶
- Author
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?
Relaxation¶
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?