:sequential_nav: next
.. _tutorial-molecular-dynamics:
Molecular Dynamics
******************
.. :Author: Alec Wills (Stony Brook University)
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 :ref:`tutorial-basic-structure-optimization`). See this same
tutorial for Quenched Molecular Dynamics.
* 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-:math:`\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:
:math:`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?
.. Geometry Relaxations
.. ====================
.. note:
.. Reviewing :ref:`tutorial-basic-structure-optimization` 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?