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 the Velocity-Verlet integration scheme in the NVE ensemble. The positions are evolved according to:
and the velocities are evolved as
Note that to compute the acceleration, the most common isotope masses are used (rather than the average isotope mass). A simple example of a 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=10,
energy_output="BOMD.md.energies",
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.energy_output.close()
myintegrator.trajectory_output.close()
This will create 2 files, BOMD.md.energies and BOMD.md.xyz that will contain the energy and structures for each frame along the trajectory. The BOMD.md.energies will look like
time Epot Ekin Etot
0.00 -1.497083671814E+02 0.000000000000E+00 -1.497083671814E+02
5.00 -1.497083678220E+02 6.403560409065E-07 -1.497083671817E+02
10.00 -1.497083697402E+02 2.557810719094E-06 -1.497083671824E+02
15.00 -1.497083729250E+02 5.741509150521E-06 -1.497083671835E+02
20.00 -1.497083773584E+02 1.017340316838E-05 -1.497083671850E+02
25.00 -1.497083830154E+02 1.582841563588E-05 -1.497083671870E+02
30.00 -1.497083898640E+02 2.267458853778E-05 -1.497083671894E+02
35.00 -1.497083978654E+02 3.067323581032E-05 -1.497083671922E+02
40.00 -1.497084069752E+02 3.977943396938E-05 -1.497083671958E+02
45.00 -1.497084171413E+02 4.994201038141E-05 -1.497083671993E+02
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