Code
= 1.661e-27
kg_per_amu = 25000; # Semicolon to suppress code output from final document num_steps
Alicia
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 Ustretch and has the following form: \[ U_{stretch} = {1 \over 2}k_{Ai}(r_{Ai}-r_{Ai,eq})^2 \] Where kAi is the spring constant of the bond, rAi,eq is the equilibrium distance between atoms A and B and rAB is the instantaneous distance between the atoms. Since I am only modeling the stretch energy in this project, I will shorten Ustretch 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 rAi 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 rAi 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, m1 and m2 are the masses of the atoms (in kg) and kb 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 Ek: \[ 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 Ep: \[ 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 1H35Cl, 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 | kAi | 516.0 | N/m |
equilibrium bond distance | rAi,eq | 1.57 x 10-10 | m |
mass of hydrogen 1 atom | m1 | 1 | amu |
mass of chlorine 35 atom | m2 | 35 | amu |
kAi and rAi,eq are from Atkins and de Paula (below), page 454. Masses for 1H35Cl are converted from atomic mass units to kilograms in the simulation code. The analytical period of oscillation of the 1H35Cl 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 1H35Cl 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 1H35Cl.
As referenced above and as a check of the simulation, calculate the analytical period of the oscillation of the 1H35Cl bond. Units are in seconds.
Since the 1H35Cl 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