Molecular Dynamics#

Modules: pyscf.md

Introduction#

Using any available method for computing gradients, molecular dynamics can be run to evolve the nuclei along their potential energy surfaces. The MD module in PySCF supports both NVE and NVT ensemble integrations. A simple example of an NVE Born-Oppenheimer MD simulation with CASSCF is given in examples/md/00-simple_nve.py

#!/usr/bin/env python

'''
A simple example to run an NVE BOMD simulation.
'''

import pyscf
import pyscf.md

mol = pyscf.M(
    atom = 'O 0 0 0; O 0 0 1.2',
    basis = 'ccpvdz',
    spin = 2)

myhf = mol.RHF().run()

# 6 orbitals, 8 electrons
mycas = myhf.CASSCF(6, 8)
myscanner = mycas.nuc_grad_method().as_scanner()

# Generate the integrator
# sets the time step to 5 a.u. and will run for 100 steps
# or for 50 a.u.
myintegrator = pyscf.md.NVE(myscanner,
                            dt=5,
                            steps=100,
                            data_output="BOMD.md.data",
                            trajectory_output="BOMD.md.xyz").run()

# Note that we can also just pass the CASSCF object directly to
# generate the integrator and it will automatically convert it to a scanner
# myintegrator = pyscf.md.NVE(mycas, dt=5, steps=100)

# Close the file streams for the energy and trajectory.
myintegrator.data_output.close()
myintegrator.trajectory_output.close()

This will create 2 files, BOMD.md.data and BOMD.md.xyz that will contain the energy and structures for each frame along the trajectory. The BOMD.md.data will look like

time          Epot                 Ekin                 Etot               T
    0.00  -1.497083671804E+02  0.000000000000E+00  -1.497083671804E+02   0.0000
    5.00  -1.497083678218E+02  6.402095343905E-07  -1.497083671816E+02   0.0674
   10.00  -1.497083697399E+02  2.557532001000E-06  -1.497083671823E+02   0.2692
   15.00  -1.497083729245E+02  5.741058994089E-06  -1.497083671835E+02   0.6043
   20.00  -1.497083773578E+02  1.017280141026E-05  -1.497083671850E+02   1.0708
   25.00  -1.497083830147E+02  1.582769986352E-05  -1.497083671870E+02   1.6660
   30.00  -1.497083898631E+02  2.267377586371E-05  -1.497083671893E+02   2.3866
   35.00  -1.497083978644E+02  3.067236427680E-05  -1.497083671921E+02   3.2285
   40.00  -1.497084069741E+02  3.977851890862E-05  -1.497083671956E+02   4.1870
   45.00  -1.497084171401E+02  4.994098786544E-05  -1.497083671991E+02   5.2567

and the BOMD.md.xyz will look like

2
MD Time 0
O           0.00000        0.00000        0.00000
O           0.00000        0.00000        1.20000
2
MD Time 5
O          -0.00000       -0.00000       -0.00001
O           0.00000        0.00000        1.20001
2
MD Time 10
O          -0.00000        0.00000       -0.00002
O           0.00000       -0.00000        1.20002
2
MD Time 15
O           0.00000        0.00000       -0.00006
O          -0.00000       -0.00000        1.20006
2
MD Time 20
O           0.00000        0.00000       -0.00010
O          -0.00000       -0.00000        1.20010
2
MD Time 25
O           0.00000        0.00000       -0.00015
O          -0.00000       -0.00000        1.20015
2
MD Time 30
O           0.00000        0.00000       -0.00022
O          -0.00000       -0.00000        1.20022
2
MD Time 35
O           0.00000        0.00000       -0.00030
O          -0.00000       -0.00000        1.20030
2
MD Time 40
O           0.00000        0.00000       -0.00039
O          -0.00000       -0.00000        1.20039
2
MD Time 45
O           0.00000        0.00000       -0.00050
O          -0.00000       -0.00000        1.20050

Integrators#

The NVE ensemble uses the Velocity-Verlet integration scheme. The positions are evolved according to:

\[\mathbf{r}(t_{i+1}) = \mathbf{r}(t_i) + (\delta t) \mathbf{v}(t_i) + \frac{1}{2}(\delta t)^2 \mathbf{a}(t_i)\]

and the velocities are evolved as

\[\mathbf{v}(t_{i+1}) = \mathbf{v}(t_i) + (\delta t) \frac{\mathbf{a}(t_{i+1}) + \mathbf{a}(t_i)}{2}\]

Note that to compute the acceleration, the most common isotope masses are used (rather than the average isotope mass).

The NVT ensemble only has the Berendsen thermostat implemented. The Berendsen thermostat updates the velocities by weakly coupling the system to a heat bath at some temperature (\(T_0\)). The temperature is therefore corrected such that the deviation exponentially decays with some time constant \(\tau\). \(\tau\) can also be understood as the strength of the coupling between the system and bath.

\[\frac{dT}{dt} = \frac{T_0 - T}{\tau}\]

At each time step, the velocities are rescaled according to

\[\mathbf{v}_\mathrm{scaled}(t) = \mathbf{v}(t)\lambda\]

where \(\lambda\) is given as

\[\lambda = \sqrt{1 + \frac{\delta t}{\tau}\left(\frac{T_0}{T} - 1\right)}\]

Note that we limit \(\lambda\) to reasonable values. That is, if \(\lambda\) is less 0.9, we set it to 0.9; and if \(\lambda\) is greater than 1.1, we set it to 1.1. See examples/md/04-mb_nvt_berendson.py for an example of how to use the Berendson theromstat.

Initial Velocity#

The initial velocity for the system can be set through the veloc kwarg in the construction of the integrator such as

integrator = pyscf.md.NVE(hf, dt=5, steps=10, veloc=init_veloc)

This can be useful for restarting a dynamics simulation with a prior velocity, or starting a simulation with a specific temperature. Additionally, the MD module also contains a function to sample the velocities from a Maxwell-Boltzmann distribution at a given temperature (in Kelvin).

init_veloc = md.distributions.MaxwellBotlzmannVelocity(mol, T=300)

The Maxwell-Botlzmann distribution at a specific temperature \(T\) for a particle of mass \(m\) is given as

\[f(v)dv = \sqrt{\frac{m}{2\pi k T}}\exp\left(-\frac{mv^2}{2kT}\right)\]

Since all of the distributions use the same pseudo-random number generator, the specific seed used can be set with the md.set_seed() function. The seed can also be set within the PYSCF_CONFIG_FILE as the SEED variable.