The importance of step size
Timing: \(\sim\) 50 minutes per substance
A deuterium atom has one proton and one neutron in its nucleus, thus about twice as heavy as hydrogen.
Here, we will investigate the effects of the time step size for the integration of the equations of motion in a microcanonical molecular dynamics simulation. For a better illustration, we will not only consider the H5O2+ molecule, but also its heavier, deuterated counterpart D5O2+.
Checkpoint files for both H5O2+ and D5O2+ from previous simulations where the molecules were thermalized are provided in the
control.in with default self-consistency convergence criteria and the
geometry.in files from the previous exercise.
Both H5O2+ and D5O2+ have the exact same electronic structure.
This means that for converging the DFT self-consistent cycle, the answer will be the same no matter the name or mass of the nuclei 1.
Therefore, one can use the same FHI-aims input files for both molecules.
What will matter, though, is that in the evolution of the nuclei through Newton's equation of motion, the mass of the atom will enter in the propagation of the equations of motion.
Therefore we will use different masses (corresponding to hydrogen and deuterion) in the checkpoint file that i-PI is reading at its initialization.
Find the masses
If you wish to have a look at where this is specified, open the checkpoint files and look for the block named
<m shape=... >. (The units are atomic mass units.)
input.xml file is also provided in the directory. Please read and complete it in the appropriate way for each of the following tasks.
In order to speed up the exercise in the live tutorials, you will be asked to simulate either H5O2+ or D5O2+.
Instructions to set up the calculations
Copy the "default"
geometry.infiles from the last exercise.
Do not add any flags that we do not mention! They are not needed and might hinder the performance of the calculation.
This should be 800 steps
input.xmlfile to run a 0.40 ps MD run in the microcanonical ensemble, using a 0.0005 ps (\(\Delta t\) = 0.5 fs) time step.
Complete all the missing fields ('xxx') in the
Run the simulation.
When it is done, plot the total energy vs. the simulation time.
Increase the time step to 0.001 ps and change the number of steps accordingly to have again a total of 0.40 ps simulation time. Take the necessary precautions to not overwrite the files and run the simulation again.
Increase the time step to 0.002 ps and change the number of steps accordingly to have again a total of 0.40 ps simulation time. Take the necessary precautions to not overwrite the files and run the simulation again.
If FHI-aims aborts the run, you can create a file named EXIT in the folder where i-PI is running to stop it in a clean
way by typing
touch EXIT in your terminal window.
Analysis of the MD runs
Together with the plot for \(\Delta t\) = 0.0005 ps, plot the total energy vs. simulation time in
xmgracefor \(\Delta t\) = 0.001 ps and \(\Delta t\) = 0.002 ps. If you're simulating D5O2+ please try also \(\Delta t\) = 0.0033 ps. You can find a python plotting script in the solutions directory. Run it by the command
python plot_Energy.py fs_0.5/ex3.out fs_1.0/ex3.out fs_2.0/ex3.out
How do the energy fluctuations develop? What is happening for the \(\Delta t\) = 0.002 ps (or for D5O2+ \(\Delta t\) = 0.0033 ps) run?
From a practical point of view, a larger time step is desirable, since it allows to assess longer trajectories in shorter computational times. Notice, however, that the \(\Delta t\) = 0.002 ps simulation diverges for H5O2+.
What do you see upon inspecting the dynamics of the molecule by opening
ex3.pos.xyz in VMD (or Jmol, or Molden, or whichever program you prefer)?
In fact, the molecule dissociates.
The reason for the dissociation is that the integrator is unable to deal with these "big" time steps. This integrator uses a simple Verlet algorithm 2, where the error in the integrator goes with \(\Delta t^4\).
If you are simulating D5O2+, the \(\Delta t\) = 0.002 ps simulation does not diverge, although the energy fluctuations become very large. You can also look at the dynamics of this molecule in jmol or VMD.
Although such large energy fluctuations would already not be accurate enough for a production run with this molecule, the fact that it does not explode illustrates an important point: the largest \(\Delta t\) that can be used in a particular integration algorithm depends on the highest vibrational frequency of the system. Since the D atoms, being heavier, have a larger vibrational period (Do you understand why is this? Think about a harmonic oscillator), the \(∆t\) that can be used for the integrator can also be larger.