Julia: Molecular Dynamics of HCl

Background

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.

Stretch Potential Energy

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}) \]

Analytical Period

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₁ and m₂ 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.

Velocity Verlet

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 \]

Kinetic, Potential, and Total Energy

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 \]

Methods

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

¹H³⁵Cl Properties

In this simulation, I model ¹H³⁵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	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

k_Ai and r_Ai,eq are from Atkins and de Paula (below), page 454. Masses for ¹H³⁵Cl are converted from atomic mass units to kilograms in the simulation code. The analytical period of oscillation of the ¹H³⁵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.

Define Constants

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.)

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

Vectors and matrices holding simulation information.

# 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

Initialize Simulation

Set the initial conditions of the simulation, including start positions and velocities of the atoms and their masses. Also set the topology to the ¹H³⁵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

Functions to Run the Simulation

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

Velocity Verlet

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 ¹H³⁵Cl.

dt = 1e-18
stretch_velocity_verlet(xs, vs, accels, tkes, tpes, one_two_bonds, one_two_bonds_kab, one_two_bonds_req, ms, dt, num_steps);  # Semicolon to suppress code output from final document

Calculate Analytical Period

As referenced above and as a check of the simulation, calculate the analytical period of the oscillation of the ¹H³⁵Cl bond. Units are in seconds.

μ = 1/(1/ms[1]+1/ms[2])
2π * sqrt(μ/one_two_bonds_kab[1])

1.1115336262465321e-14

Results

X Coordinate of Hydrogen

Since the ¹H³⁵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:

x_axis = eachindex(xs[:, 2, 1]) / 1000
h_x_axis_trajectory = xs[:, 2, 1]
hydrogen_plot = plot(x_axis, h_x_axis_trajectory, xlabel="Time (fs)", ylabel="H x position (m)", legend=false)
display(hydrogen_plot)

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

Kinetic, Potential, and Total Energy

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):

plot(x_axis, tkes, xlabel="Time (fs)", ylabel="J", label="Ek")
plot!(x_axis, tpes, label="Ep")
energy_plot = plot!(x_axis, tpes + tkes, label="total")
display(energy_plot)

┌ 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

Discussion

This was my first project in the Julia programming language, and I am pleased with the results. The analytical period agrees with the simulation period of ¹H³⁵Cl. I find that very exciting! However, the total energy of the system increases slightly over the run of the trajectory, which means that energy is not being conserved by the algorithm. To combat this upward energy creep, I made the time step 1 x 10^-18 s. This seemed to help the energy creep, but did not ameliorate the situation entirely.

Future directions of this work could include:

1. Incorporating U_angle energy and simulating ¹H₂¹⁶O.

2. Simulating multiple molecules at once.

Thanks for reading!

References

Atkins, P. W. & De Paula, J. Atkins’ Physical Chemistry. (W.H. Freeman, New York, 2006).

Cramer, C. J. Essentials of Computational Chemistry: Theories and Models. (Wiley, Chichester, West Sussex, England ; Hoboken, NJ, 2004).