: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 relaxation include are ``CG``, ``Broyden``, and ``FIRE`` (some of which are covered in :ref:`tutorial-basic-structure-optimization`). 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-:math:`\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 :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? 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?