Born-Oppenheimer Molecular Dynamics

This tutorial aims at presenting how to perform a Born-Oppenheimer Molecular Dynamics (BOMD) simulation with Octopus, as well as the related tools that can be used to analyze the data generated by such a run. A typical BOMD simulation goes as follows. One first computes the ground state of the system under interest. Then, one performs the actual time evolution, in which the ions are evolved, using a velocity Verlet algorithm, and the electrons remain in their ground state at each step of the time evolution. A SCF cycle is then done for each step of the molecular dynamics, after which the electronic forces are computed.

Ground-state input file

As always, we will start with a simple input file. As an example, in this tutorial we will look at the case of a Water molecule:

CalculationMode = gs
FromScratch = yes

BoxShape = parallelepiped
Lsize = 7

%Coordinates
 'O'  |  0.000000  | -0.553586  |  0.000000
 'H'  |  1.429937  |  0.553586  |  0.000000
 'H'  | -1.429937  |  0.553586  |  0.000000
%

Spacing = 0.5

Eigensolver = chebyshev_filter
ExtraStates = 4

Now run Octopus using the above input file. This produces the usual output, as described in previous tutorials.

From the static/info file, one can find the usual information of the ground-state of the water molecule. As H2O is a polar molecule, we can in particular find the dipole moment of the molecule, which is here of about 1.33 Debye, quite far from the experimental value of 1.855 Debye, see Ref.¹, as the calculation is not converged here.

Born-Oppenheimer molecular dynamics

Let us now perform the BOMD calculation. For this, we will use the following input file:

CalculationMode = td
FromScratch = yes

BoxShape = parallelepiped
Lsize = 7

%Coordinates
 'O'  |  0.000000  | -0.553586  |  0.000000
 'H'  |  1.429937  |  0.553586  |  0.000000
 'H'  | -1.429937  |  0.553586  |  0.000000
%

Spacing = 0.5

Eigensolver = chebyshev_filter
ChebyshevFilterNIter = 0
OptimizeChebyshevFilterDegree = no
ExtraStates = 4

TDDynamics = bo
TDTimeStep = 20.0
TDPropagationTime = 50*fs
MoveIons = yes

ParStates = no

In order to perform the BOMD, we changed few thing in the above input file. First, we change the CalculationMode to td, as we are now perforing a time evolution. In order to tell Octopus that we want a BOMD calculation, we set the variable TDDynamics to bo. We further activated the motion of the ions using MoveIons = yes.

In order to make the run fast, we also here deactivated the parallelization over states, using ParStates = no, and tuned some parameters of the Chebyshev filtering to avoid recomputing initial bounds of the filters, as each BOMD step provides a good initial starting point for the next one.

It should be noted that the parameters used in this tutorial (box size, grid spacing, tolerance on the relative density…) are not converged, and should be adjusted for any practical calculation.

Analysis of the results

In order to analyze the motion of the molecule, we now compute the velocity correlation function and compute its Fourier transform. This is done using the tool oct-vibrational_spectrum. As we performed here a relatively short calculation, we want to use all the time steps for the analysis, so we use

VibrationalSpectrumTimeStepFactor = 1

After running the utility, we obtain the files td.general/velocity_autocorrelation and td.general/vibrational_spectrum. Plotting the file td.general/vibrational_spectrum shows the spectrum of the velocity correlation function, displaying peaks at the vibrational modes of H2O, see Fig. 1. We observe a main peak around 3000 cm^-1, as well as other side structures.

Generate a movie from the actual data

Let us finally produce a movie from the coordinates generated during the BOMD run. For this, one can use the utility oct-xyz-anim. This generates a file called td.general/movie.xyz, that can be later converted to an animated GIF using for instance XCrysden.

xcrysden --xyz td.general/movie.xyz

The resulting movie is shown in Fig. 2 below, for a longer simulation of 100fs, using a tighter relative density convergence criterim of 1e-7, and otherwise the same input parameters as the above input. It is possible to use the variable AnimationSampling to select how many frame are used for the movie. In order to produce the movie below, we used AnimationSampling=1.

Adding a Nose-Hoover thermostat

The previous section performs a Born-Oppenheimer molecular dynamics run without explicitly controlling the ionic temperature during the propagation. This is a useful first calculation, since it allows one to check that the time step and the electronic convergence are stable. In many applications, however, one wants to equilibrate the ions around a target temperature. This can be done by adding an ionic thermostat to the BOMD run.

As a simple example, we now repeat the BOMD calculation using a Nose-Hoover thermostat at 300 K. Starting from the previous BOMD input file, we can use the same ground-state, grid, box, and time-propagation settings, and add the following lines:

ExperimentalFeatures = yes

RandomVelocityTemp = 300
Thermostat = nose_hoover
ThermostatMass = 1.0

TemperatureFunction = "temperature"

%TDFunctions
  "temperature" | tdf_cw | 300
%

Here, we start the trajectory with random ionic velocities corresponding to 300 K, given an Eistein distribution. For a fully reproducible tutorial run, one should replace this by an explicit Velocities block.

Here Thermostat = nose_hoover activates the Nose-Hoover thermostat, while TemperatureFunction points to a time-dependent function giving the target temperature. In this example the function is constant, using tdf_cw, and the target temperature is 300 K. The value of ThermostatMass controls how quickly the thermostat variable reacts to temperature fluctuations. A very small value can lead to aggressive temperature oscillations, while a very large value can make equilibration slow. The value used here is only a reasonable tutorial starting point, not a generally converged parameter.

The ionic temperature in a thermostatted molecular dynamics simulation should fluctuate around the target value; it should not be expected to remain exactly equal to 300 K at every time step. It is therefore useful to compare this run with the previous BOMD trajectory without a thermostat. For a more systematic exercise, one can for instance repeat the calculation with ThermostatMass = 0.1, 1.0, and 10.0, and compare the temperature and energy files written in td.general.

A temperature ramp can also be used instead of a constant temperature. For example, the following function ramps the target temperature from 50 K to 300 K during the first 10 fs, and then keeps it at 300 K:

TemperatureFunction = "temperature"

%TDFunctions
  "temperature" | tdf_from_expr | "min(300, 50 + 250*t/(10*fs))"
%

The same analysis described above can also be used for the thermostatted trajectory. When using a thermostat, the initial part of the trajectory should normally be regarded as equilibration and omitted from quantitative analysis. The short trajectories in this tutorial are intended to demonstrate how to use the functionality and should not be considered for producing converged vibrational spectra or thermodynamic averages.

References