Please upload one PDF file to Gradescope (Entry code: G64YK2).
For this homework you will implement code (on PrairieLearn) to perform molecular dynamics simulations using periodic boundary conditions. You specifically should implement functions that determine the coordinates of a particle, assuming periodic boundary conditions, that determine the vector that connects two particles (as well as its length), and functions for forces and energies (please see PrairieLearn for details and instructions).
Often, the goal of atomic simulations is to measure bulk properties of a system. “Bulk” means that one wishes to describe the many-body interactions of an infinite system (the so-called thermodynamic limit). A particle in a bulk system should interact with an infinite volume of particles around it; it cannot be at a surface or boundary, where it will experience different interactions.
Of course, we can simulate only a finite system on a computer, so we must make an approximation to the thermodynamic limit. Typically, the solution is to implement Periodic Boundary Conditions (PBC). This is similar to old video games where, when you come out of one side of the screen, you re-enter on the other side.
The effect of PBC is that a particle in the simulation box ``sees’’ all other particles replicated on all sides of the simulation box, thus creating the effective interactions of a bulk system. A technical limitation of PBC is that there will be artificial spatial correlations on the order of the box size. However, it is a standard simulation technique and this artifact can be analyzed and mitigated by studying the effect of systematically increasing the system size.
It is important to ensure that you are enforcing PBC correctly; it can be tricky to get this right!
The Lennard-Jones Fluid. Once you have implemented the corresponding functions correctly, you will implement a molecular dynamics code to model a Lennard-Jones gas. The Lennard-Jones potential is: \[ V\left(r\right)=4\epsilon\left[\left(\frac{\sigma}{r}\right)^{12}-\left(\frac{\sigma}{r}\right)^{6}\right]. \] Here, there is no external potential; the only force you need to calculate is from the pairwise interactions given by the Lennard-Jones potential. (Recall that the force \(\vec{F} = -\nabla V = -V'(r) \hat{\mathbf{r}}\).)
For simplicity, we will work in reduced units (see e.g. p.40 of Frenkel & Smit). This means lengths are expressed in terms of \(\sigma\) and the temperature is in units of the Lennard-Jones energy parameter, \(\epsilon\). Thus, the reduced temperature is \(T_* = \frac{k_BT}{\epsilon}\). Then the Lennard-Jones potential simplifies to \(V\left(r\right)=4[(\frac{1}{r})^{12}-(\frac{1}{r})^{6}]\) and we need only specify a density and temperature at which to perform our simulation.
We will simulate \(N=64\) particles at a number density of \(\rho = 0.8442\) and a reduced temperature of \(T_* = 0.728\). These values correspond to a box side length of approximately \(4.232317\). These initial conditions are already set in the code provided below. As described in Frenkel & Smit (p. 66), we initialize the particles with random velocities, adjust the system to have zero net momentum, then rescale all velocities so the system starts at the desired temperature of \(0.728\).
Note: You are welcome to spend some time experimenting with other ways to initialize your simulation, but it is non-trivial. For example, starting all your particles at the same spot is a bad idea. The Lennard-Jones potential has a repulsive core, so your system will start with effectively infinite potential energy. Because you are working in the micro-canonical ensemble, energy is conserved, so the kinetic energy will explode and the simulation will be unstable.
Potential truncation. We will implement the “cut and shifted” Lennard-Jones potential in our small simulation cell. That is, what we will actually use is \[ V_{cs}(r) = \begin{cases} V(r)-V(r_c), & r<r_c\\ 0 & r\geq r_c \end{cases} \] where \(v\) is the pair potential and \(r_c\) is the cutoff radius of interaction. The reason for the shift is to remove the discontinuity at \(r_c\). Technically, we should add a “tail correction” back to the total energy \[V_\text{tail} = \frac{1}{2} \rho N \int{ dr 4 \pi r^2 g(r) [ V(r) - V_{cs}(r)]}\] However, we will ignore this correction to simplify the code. Since this correction is a constant shift of the total energy, it does not affect the dynamics in the microcanonical ensemble.
Computing the Energy. You’ll also want to get some output from your simulation. We’re usually less interested in the precise positions of all the atoms as we are in some physical properties. The simplest one to examine is the energy. In case you’ve forgotten, let us remind you precisely what the energy means in this context: Firstly, the energy of the system is the sum of the kinetic and potential energies. The kinetic energy is simply the sum over all atoms of \(\frac{1}{2} m v_i^2\), where \(v_i\) is the velocity of the \(i\)-th atom and \(m\) its mass. The potential energy is the sum over all distinct pairs of atoms of the Lennard-Jones potential (described above) between the atoms. Just to be clear, let’s take an example: If we had 3 atoms in the box, we might write express the potential energy as \[ PE = V_{LJ}( \mathrm{Distance}[R1,R2] ) + V_{LJ}(\mathrm{Distance}[R1,R3] ) + V_{LJ}(\mathrm{Distance}[R2,R3] ), \] where \[ V_{LJ}(r) = 4.0 * [1.0/(r^{12}) - 1.0/(r^6)]. \] Since the total energy is a conserved quantity, we expect it will stay almost constant as our simulation proceeds.
To save you some effort, you will not have to implement the entire molecular dynamics code. You should peruse the following code to see the flow of logic in a simple molecular dynamics simulation:
import numpy as np
def minimum_image(r, L):
"""
required for: displacement_table(), advance()
args:
r: array of any shape
L: side length of cubic box
returns:
array of the same shape as r,
the minimum image of r
"""
pass
def cubic_lattice(tiling, L):
"""
required for: initialization
args:
tiling (int): determines number of coordinates,
by tiling^3
L (float): side length of simulation box
returns:
array of shape (tiling**3, 3): coordinates on a cubic lattice,
all between -0.5L and 0.5L
"""
pass
def get_temperature(mass, velocities):
"""
calculates the instantaneous temperature
required for: initial_velocities()
args:
mass (float): mass of particles;
it is assumed all particles have the same mass
velocities (array): velocities of particles,
assumed to have shape (N, 3)
returns:
float: temperature according to equipartition
"""
pass
def initial_velocities(N, m, T):
"""
initialize velocities at a desired temperature
required for: initialization
args:
N (int): number of particles
m (float): mass of particles
T (float): desired temperature
returns:
array: initial velocities, with shape (N, 3)
"""
pass
def displacement_table(coordinates, L):
"""
required for: force(), advance()
args:
coordinates (array): coordinates of particles,
assumed to have shape (N, 3)
e.g. coordinates[3,0] should give the x component
of particle 3
L (float): side length of cubic box,
must be known in order to compute minimum image
returns:
array: table of displacements r
such that r[i,j] is the minimum image of
coordinates[i] - coordinates[j]
"""
pass
def kinetic(m, v):
"""
required for measurement
args:
m (float): mass of particles
v (array): velocities of particles,
assumed to be a 2D array of shape (N, 3)
returns:
float: total kinetic energy
"""
pass
def potential(dist, rc):
"""
required for measurement
args:
dist (array): distance table with shape (N, N)
i.e. dist[i,j] is the distance
between particle i and particle j
in the minimum image convention
note that the diagonal of dist can be zero
rc (float): cutoff distance for interaction
i.e. if dist[i,j] > rc, the pair potential between
i and j will be 0
returns:
float: total potential energy
"""
pass
def force(disp, dist, rc):
"""
Warning: this computes the force on ALL particles
rather than just one, as in the PrairieLearn exercise
also: instead of taking the particle coordinates and length of box,
this takes the displacement and distance tables,
since these are computed externally
required for: advance()
args:
disp (array): displacement table,
with shape (N, N, 3)
dist (array): distance table, with shape (N, N)
can be calculated from displacement table,
but since there is a separate copy available
it is just passed in here
rc (float): cutoff distance for interaction
i.e. if dist[i,j] > rc, particle i will feel no force
from particle j
returns:
array: forces f on all particles, with shape (N, 3)
i.e. f[3,0] gives the force on particle i
in the x direction
"""
pass
def advance(pos, vel, mass, dt, disp, dist, rc, L):
"""
advance system according to velocity verlet
args:
pos (array): coordinates of particles
val (array): velocities of particles
mass (float): mass of particles
dt (float): timestep by which to advance
disp (array): displacement table
dist (array): distance table
rc (float): cutoff
L (float): length of cubic box
returns:
array, array, array, array:
new positions, new velocities, new displacement table,
and new distance table
"""
accel = force(disp, dist, rc) / mass
#move
vel_half = vel + 0.5*dt*accel
pos_new = pos + dt*vel_half
pos_new = minimum_image(pos_new, L) #not strictly necessary
disp_new = displacement_table(pos_new, L)
dist_new = np.linalg.norm(disp_new, axis=-1)
#repeat force calculation for new pos
accel = force(disp_new, dist_new, rc) / mass
#finish move
vel_new = vel_half + 0.5*dt*accel
return pos_new, vel_new, disp_new, dist_new
#for reproducibility
np.random.seed(466)
#parameters
num_particles = 64
temperature = 0.728
length = 4.232317
cutoff = 2.08
mass = 48
timestep = 0.01
#system
coordinates = cubic_lattice(4, length)
velocities = initial_velocities(num_particles, mass, temperature)
#tables required to compute quantities like forces, energies
displacements = displacement_table(coordinates, length)
distances = np.linalg.norm(displacements, axis=-1)
#advance and record energies
#it can also be useful to save coordinates and velocities
KE = []
PE = []
for _ in range(3000):
coordinates, velocities, displacements, distances = advance(coordinates,\
velocities, mass, timestep, displacements, distances, cutoff,\
length)
PE.append(potential(distances, cutoff))
KE.append(kinetic(mass, velocities))
#write to disk/analyze Note: The slowest part of a molecular dynamics simulation is computing the distances and the forces. If you are interested in speeding up your code, consider the following tips.
numba can help you implement functions that run as
fast as numpy code).You can vary the time step by changing the variable h just above the Verlet time-stepping functions. Using a time step of 0.01, run your simulation for at least 1000 steps (a run of 1000 steps should take a few seconds to minutes to complete). Plot an energy trace of the run (the total energy or energy per particle at each time step in the simulation) and include in your report. (As a check on your energy trace, look to see if your total energy remains approximately constant.)
Initially, we used a time step of \(0.01\), which is probably not optimal and may not even be stable! It is important to be in the habit of ensuring that your time step is neither too small (your system will never evolve from its initial configuration) nor too large (your numerical integration will be unstable and the velocities will explode).
Here, our goal is to find the largest time step that gives a stable mean and keeps total energy fluctuations within 5,% of the mean. This would be the largest acceptable time step, and you will generally want to work with a time step that is an order of magnitude smaller than this to mitigate the time step error.
Make a graph of the standard deviation of the total energy versus time step size. In other words, you are trying to see how the energy fluctuates with respect to your choice of time step. For this you need to make several simulations using as many as 10 different time steps.
What happens if you pick a very large time step? Why?
From your graph, determine the largest time step one can use and still maintain energy conservation (a stable mean and fluctuations within 5,% of the mean).
So far, we have primarily used the velocity Verlet integrator. This isn’t the most obvious choice. Instead, one might imagine that we could simply Taylor expand \[ r(t+\delta t)=r(t)+r'(t)\delta t + \frac{1}{2}r''(t)\delta t^2 + ...\rightarrow r(t+\delta t)=r(t)+v(t)\delta t + \frac{1}{2}a(t)\delta t^2 + ... \] Try using this. Discuss what happens (look at the energy)! Why don’t we use this integrator instead of Verlet?
Tip: Write separate functions for this integrator rather than changing the Verlet functions.
Prove that the Verlet algorithm is time-reversal invariant. That is, show that if we use \(r(t+h)\) and \(r(t)\) as inputs we get \(r(t-h)\) (up to round-off errors).
State two possible problems with finite-precision arithmetic about the Verlet time-stepping algorithm.