Optional Lesson on Choosing intial bond-length and velocity

Overview

Teaching: 45 min
Exercises: 15 min

Questions

• How do we initialize the position and velocity of the atoms for our molecular dynamics simulation?

Objectives

• To demonstrate a reasonable strategy for initializing a molecular dynamics simulation of vibrational motion.

For further consideration: What makes a “sensible range of values” for position and velocity?

Since we have only computed the potential energy surface between around 1 and 4 atomic units, that provides a reasonable bound on the random initial positions we will choose. The question remains what represents a suitable range of values for the initial velocity? One strategy is to use the fact that we can estimate the expectation value of kinetic energy for a very similar system (the Harmonic oscillator) in the ground state as follows: \begin{equation} \langle T \rangle = \frac{1}{2} E_g, \end{equation} where \(E_g\) is the ground state of the Harmonic oscillator (this is making use of the Virial theorem). We can easily find the ground state energy in the Harmonic oscillator approximation of HF using our frequency calculation described above as \begin{equation} E_g = \frac{1}{2} h \nu, \end{equation} which implies the kinetic energy expectation value is \begin{equation} \langle T \rangle = \frac{h}{8 \pi} \sqrt{\frac{k}{\mu}}. \end{equation} Since we can say classically that the kinetic energy is given by \(T = \frac{1}{2}\mu v^2\), we can estimate the velocity of the bond stretch as follows: \begin{equation} v = \sqrt{\frac{2 \langle T \rangle}{\mu}} = \sqrt{ \frac{\hbar \sqrt{\frac{k}{\mu}}}{2\mu}} \end{equation} where we have simplified using the fact that \(\hbar = \frac{h}{2\pi}\) (\(\hbar\) has the value 1 in the atomic unit system we are using up to this point!). We will assume that a reasonable range of velocities spans plus or minus 3 times this “ground-state” velocity.

### define "ground-state" velocity for each level of theory
v_RHF = np.sqrt( np.sqrt(RHF_k/mu)/(2*mu))
v_MP2 = np.sqrt( np.sqrt(MP2_k/mu)/(2*mu))
v_CCSD = np.sqrt( np.sqrt(CCSD_k/mu)/(2*mu))


### get random position and velocity for RHF HF within a reasonable range
r_init = np.random.uniform(0.75*RHF_Req,2*RHF_Req)
v_init = np.random.uniform(-2*v_RHF,2*v_RHF)

### print initial position and velocity
print("Initial separation is ",r_init, "atomic units")
print("Initial velocity is   ",v_init, "atomic units")


### get initial force on the particle based on its separation
RHF_F_init = -1*RHF_Force(r_init)
print("Initial Force is ", RHF_F_init, "atomic units")

Key Points

• We can utilize the virial theorem from quantum mechanics to provide a reasonable value of the classical velocity associated with vibrational motion.