## Code

```
= 1.661e-27
kg_per_amu = 25000; # Semicolon to suppress code output from final document num_steps
```

julia

physics

chemistry

Author

Alicia

Published

March 23, 2024

Molecular dynamics models the motion of atoms within molecules using classical mechanics calculations. Molecular dynamics envisions atoms as balls on springs and calculates the accelerations of atoms within the system in a potential energy force field that has terms for several types of motions in the molecule. In this post, I will use the letter *U* to denote potential energy. Typical force fields model the following potentials: \[ U = U_{stretch} + U_{angle} + U_{dihedral} + U_{improper} \] For a complete description of the force field terms and motions, I refer the reader to Chapter 3 of the Cramer text referenced at the bottom of this document. In this project, I only model bond stretch energy so I could keep the code simple. The notation I will use in this post and code mostly comes from Chapter 3 of the Cramer text.

The stretch potential energy of a bond between atoms A and i is denoted U_{stretch} and has the following form: \[ U_{stretch} = {1 \over 2}k_{Ai}(r_{Ai}-r_{Ai,eq})^2 \] Where k_{Ai} is the spring constant of the bond, r_{Ai,eq} is the equilibrium distance between atoms A and B and r_{AB} is the instantaneous distance between the atoms. Since I am only modeling the stretch energy in this project, I will shorten U_{stretch} to U from here onward.

The first derivative of the stretch energy which is used to calculate the stretch forces acting upon an atom by the other atoms it is bonded to is: \[ {\partial U \over \partial x_A}=\sum_{i} {\partial U \over \partial r_{Ai}} {\partial r_{Ai} \over \partial x_A} \] Where the summation is over the atoms bonded to A and i refers to the other atoms. Similarly to x axis potential, the y and z potentials are: \[ {\partial U \over \partial y_A}=\sum_{i} {\partial U \over \partial r_{Ai}} {\partial r_{Ai} \over \partial y_A} \] \[ {\partial U \over \partial z_A}=\sum_{i} {\partial U \over \partial r_{Ai}} {\partial r_{Ai} \over \partial z_A} \] The partial derivatives of r_{Ai} for each x, y, z coordinate are: \[ {\partial r_{Ai} \over \partial x_A} = {x_A-x_i \over \sqrt{(x_A-x_i)^2 + (y_A-y_i)^2 + (z_A-z_i)^2}} \] \[ {\partial r_{Ai} \over \partial y_A} = {y_A-y_i \over \sqrt{(x_A-x_i)^2 + (y_A-y_i)^2 + (z_A-z_i)^2}} \] \[ {\partial r_{Ai} \over \partial z_A} = {z_A-z_i \over \sqrt{(x_A-x_i)^2 + (y_A-y_i)^2 + (z_A-z_i)^2}} \] Finally, the partial derivative of stretch energy with respect to bond distance r_{Ai} is: \[ {\partial U \over \partial r_{Ai}} = k_{Ai}(r_{Ai}-r_{Ai,eq}) \]

The analytical period of a harmonic oscillator that models the interatomic bonds in this simulation is: \[ {1 \over \mu} = {1 \over m_1} + {1 \over m_2} \]

\[ T = \sqrt{\mu \over k_b} \] where mu is the reduced mass, m_{1} and m_{2} are the masses of the atoms (in kg) and k_{b} is the spring constant of the bond. T is the period of the oscillation in seconds. While I do not use this equation in the simulation, I do use it to verify the simulation is producing physically valid output.

To propagate the trajectory over time, I use the velocity Verlet algorithm. At each time step, I perform the algorithm for each atom. The steps of the algorithm are as follows: \[ \mathbf x_{t+1} = \mathbf x_t + \mathbf {v}_t \Delta t + \mathbf a_t (\Delta t)^2 \] \[ \mathbf v_{t+{1 \over 2}} = \mathbf v_t + {1 \over 2} \mathbf a_i \Delta t \] \[ \mathbf a_{i+1} = -{\mathbf \nabla_i \over m} \] \[ \mathbf v_{i+1} = \mathbf v_{t+{1 \over 2}} + {1 \over 2}\mathbf a_{i+1} \Delta t \]

As diagnostic metrics to assess the proper behavior of the simulation, I calculate the kinetic, potential, and total energy of the system at each time step. For the kinetic energy, I just sum over the kinetic energies of all the atoms to arrive at the total kinetic energy E_{k}: \[ E_k = \sum_i^N {1 \over 2} m_i v_i^2 \] Since the only energy is stretch energy in this simulation, I sum over the stretch energy over all atoms to find the total potential energy E_{p}: \[ E_p = \sum_i^N U_i \]

In this section I describe the code that I wrote to run the simulation.

In this simulation, I model ^{1}H^{35}Cl, as it is a classic molecule used to learn physical chemistry. The values I use to model the bond motion are:

Quantity | Variable | Value | Units |
---|---|---|---|

spring constant | k_{Ai} |
516.0 | N/m |

equilibrium bond distance | r_{Ai,eq} |
1.57 x 10^{-10} |
m |

mass of hydrogen 1 atom | m_{1} |
1 | amu |

mass of chlorine 35 atom | m_{2} |
35 | amu |

k_{Ai} and r_{Ai,eq} are from Atkins and de Paula (below), page 454. Masses for ^{1}H^{35}Cl are converted from atomic mass units to kilograms in the simulation code. The analytical period of oscillation of the ^{1}H^{35}Cl bond is **11 fs**, as calcualted below. This value will be used to analyze the period of oscillation in the results and discussion sections.

`kg_per_amu`

converts amu to kg. `num_steps`

is the total number of time steps to run the simulation (delta t is 1 x 10^{-18} s, as will be explained below.)

```
# Positions (xs), velocities (vs), and accelerations (accels) arrays:
# First axis is timestep
# Second axis are atoms
# Third axis is x,y,z (meters for position, m/s for velocities)
vs = zeros(Float64, num_steps, 2, 3)
xs = zeros(Float64, num_steps, 2, 3)
accels = zeros(Float64, num_steps, 2, 3)
# Masses: The masses of each atom (kg)
ms = zeros(Float64, 2)
# Total kinetic energies: The total kinetic energy of the system
# at each timestep.
tkes = zeros(Float64, num_steps)
# Total potential energies: Total potential energy of the system
# at each timestep
tpes = zeros(Float64, num_steps); # Semicolon to suppress code output from final document
```

Set the initial conditions of the simulation, including start positions and velocities of the atoms and their masses. Also set the topology to the ^{1}H^{35}Cl bond here.

```
# Equilibrium bond length for HCl
r_ab_eq_hcl = 1.57e-10
# Assume Cl is at 0,0,0 and H lies along the x-axis
# HCl equilibrium bond length
xs[1, 2, :] = [r_ab_eq_hcl*0.999, 0.0, 0.0]
# Masses, Cl first then H
ms[1] = 35 * kg_per_amu
ms[2] = 1 * kg_per_amu
# 1-2 Bonds
# Rows are bonds, columns are atoms participating in bond
# Note: This is specifying edges on a graph, so 1-2 also has 2-1
one_two_bonds = [1 2; 2 1]
# 1-2 Bonds, stretch constants
# Note: There is the same constant for each direction of the bond
# HCl bond constant 516 N/m according to Atkins and de Paula, pg. 454
one_two_bonds_kab = [516.0 516.0]
# 1-2 Bonds, equilibrium distances
# Note: There is a distance for each direction of the bond
one_two_bonds_req = [r_ab_eq_hcl r_ab_eq_hcl]; # Semicolon to suppress code output from final document
```

While the code above defines the data for the simulation to operate on, these functions are what perform the calculations on those data.

```
# Distance from atom a to b
function r_ab(a, b)
sqrt(sum((a-b).^2))
end
# Stretch energy for a 1-2 bond
function stretch_energy(a, b, k_ab, r_ab_eq)
0.5*k_ab*(r_ab(a, b)-r_ab_eq)^2
end
# Kinetic energy for an atom
function kinetic_energy(vs::Array{Float64}, ms::Array{Float64}, timestep, atom)
v = vs[timestep, atom, :]
m = ms[atom]
0.5 * m * sum(v.^2)
end
# Stretch energy gradient for a single 1-2 bond
function stretch_gradient(a, b, k_ab, r_ab_eq)
du_drab = 0.5 * k_ab * (2*r_ab(a, b)-2*r_ab_eq)
drab_dxa = (a[1]-b[1])/r_ab(a, b)
drab_dya = (a[2]-b[2])/r_ab(a, b)
drab_dza = (a[3]-b[3])/r_ab(a, b)
[drab_dxa, drab_dya, drab_dza] * du_drab
end
# Propagate 1-2 bond stretch trajectories
function stretch_velocity_verlet(xs, vs, accels, tkes, tpes, one_two_bonds, one_two_bonds_kab, one_two_bonds_req, ms, dt, num_steps)
for time_i in 1:num_steps-1
for bond_i in [1 2]
k_ab = one_two_bonds_kab[bond_i]
r_eq = one_two_bonds_req[bond_i]
atom_a = one_two_bonds[bond_i, 1]
atom_b = one_two_bonds[bond_i, 2]
xs[time_i+1, atom_a, :] = xs[time_i, atom_a, :] + vs[time_i, atom_a, :] * dt + accels[time_i, atom_a, :] * dt^2
v_mid = vs[time_i, atom_a, :] + 0.5 * accels[time_i, atom_a, :] * dt
accels[time_i+1, atom_a, :] = -stretch_gradient(xs[time_i, atom_a, :], xs[time_i, atom_b, :], k_ab, r_eq) / ms[atom_a]
vs[time_i+1, atom_a, :] = v_mid + 0.5 * accels[time_i+1, atom_a, :] * dt
end
tkes[time_i] = total_kinetic_energy(vs, ms, time_i)
tpes[time_i] = total_stretch_energy(xs, one_two_bonds, one_two_bonds_kab, one_two_bonds_req, time_i)
end
tkes[num_steps] = total_kinetic_energy(vs, ms, num_steps)
tpes[num_steps] = total_stretch_energy(xs, one_two_bonds, one_two_bonds_kab, one_two_bonds_req, num_steps)
return nothing
end
# Determine the total kinetic energy of the system
function total_kinetic_energy(vs::Array{Float64}, ms::Array{Float64}, timestep::Int)
sum([kinetic_energy(vs, ms, timestep, atom) for atom in eachindex(ms)])
end
# Determine total stretch energy of the system.
function total_stretch_energy(xs::Array{Float64}, one_two_bonds::Array{Int}, one_two_bonds_kab::Array{Float64}, one_two_bonds_req::Array{Float64}, timestep::Int)
stretch_energies = zeros(Float64, length(one_two_bonds_kab))
for i in eachindex(one_two_bonds_kab)
atom_a = one_two_bonds[i, 1]
atom_b = one_two_bonds[i, 2]
k_ab = one_two_bonds_kab[i]
r_eq = one_two_bonds_req[i]
pos_a = xs[timestep, atom_a, 1:3]
pos_b = xs[timestep, atom_b, 1:3]
stretch_energies[i] = stretch_energy(pos_a, pos_b, k_ab, r_eq)
end
# Divide by 2.0 to prevent double counting the 1-2 bonds
sum(stretch_energies) / 2.0
end; # Semicolon to suppress code output from final document
```

Now for the fun part: Run the simulation! The following code block calls the functions to operate on the data above and calculate the trajectory of ^{1}H^{35}Cl.

As referenced above and as a check of the simulation, calculate the analytical period of the oscillation of the ^{1}H^{35}Cl bond. Units are in seconds.

Since the ^{1}H^{35}Cl molecule lies on the x axis, to find the period of oscillation we can look at plot of the hydrogen atom’s x-coordinate versus time, as shown below:

The plot shows the period of oscillation is 11 fs, which agrees with the analytical calculation above.

Ideally, at each time step the total energy of the system should remain constant. To verify this, the following plot displays the kinetic, potential, and total energy in the simulation (please ignore the warning messages about plot ticks):

```
┌ Warning: No strict ticks found
└ @ PlotUtils C:\Users\amfke\.julia\packages\PlotUtils\8mrSm\src\ticks.jl:194
┌ Warning: No strict ticks found
└ @ PlotUtils C:\Users\amfke\.julia\packages\PlotUtils\8mrSm\src\ticks.jl:194
┌ Warning: No strict ticks found
└ @ PlotUtils C:\Users\amfke\.julia\packages\PlotUtils\8mrSm\src\ticks.jl:194
┌ Warning: No strict ticks found
└ @ PlotUtils C:\Users\amfke\.julia\packages\PlotUtils\8mrSm\src\ticks.jl:194
┌ Warning: No strict ticks found
└ @ PlotUtils C:\Users\amfke\.julia\packages\PlotUtils\8mrSm\src\ticks.jl:194
┌ Warning: No strict ticks found
└ @ PlotUtils C:\Users\amfke\.julia\packages\PlotUtils\8mrSm\src\ticks.jl:194
┌ Warning: No strict ticks found
└ @ PlotUtils C:\Users\amfke\.julia\packages\PlotUtils\8mrSm\src\ticks.jl:194
```