Week 07 Lecture: Molecular Dynamics Concepts 2
Contents
Week 07 Lecture: Molecular Dynamics Concepts 2¶
temperature
forcefields
non-bonded interactions
periodic boundaries
Links:
Kahn-Academy Maxwell-Boltzmann
MolSSI MMTools
A molecular definition of temperature¶
What is temperature? What really makes a cup of cold water different from a cup of hot water?
Temperature is a speedometer¶
The hotter water molecules are moving faster. The relationship between the temperature and the mean squared velocity of the particles (\(\langle{v^2}\rangle\)) is:
where \(n\) is the number of spatial dimensions, \(m\) is the particle mass, and \(k_B\) is Boltzmann’s constant.
\(T\) here is in Kelvin, showing that at absolute zero (-273.15 C), the average squared velocity is zero and all particles are perfectly stationary.
Maxwell-Boltzmann distribution¶
The distribution of molecular speeds is given by:
where \(P(|v|)\) is the probability of a particular speed (or velocity magnitude, \(|v|\)).
import matplotlib.pyplot as plt
import numpy as np
def mb_prob(v,kt_over_m):
return (1/(2*np.pi*kt_over_m))**(1.5)*v**2*np.exp(-np.square(v)/(2*kt_over_m))
v = np.linspace(0,20,100)
for kt_over_m in [1,2,5]:
plt.plot(v,mb_prob(v,kt_over_m),label=f'kT/m = {kt_over_m}')
plt.xlabel("v")
plt.ylabel("Probability")
plt.legend()
<matplotlib.legend.Legend at 0x7fe09d01bd60>
A particular component of the velocity (say, \(v_x\)) is distributed according to:
Notice the difference in the exponent (1/2 vs 3/2) and that the \(v^2\) term is missing.
def mb_prob_x(vx,kt_over_m):
return (1/(2*np.pi*kt_over_m))**(0.5)*np.exp(-np.square(vx)/(2*kt_over_m))
vx = np.linspace(-10,10,100)
for kt_over_m in [1,2,5]:
plt.plot(vx,mb_prob_x(vx,kt_over_m),label=f'kT/m = {kt_over_m}')
plt.xlabel("v_x")
plt.ylabel("Probability")
plt.legend()
<matplotlib.legend.Legend at 0x7fe09d23d670>
Initializing particle velocities¶
How do we initialize particle velocities in a way that is random, but still consistent with the MB distribution?
Use the cumulative distribution function for \(P(v_x)\):
This is equal to:
where \(\text{erf}\) is the “error function”.
from scipy.special import erf
def cdf(vx,kt_over_m):
return 0.5*(erf(np.sqrt(1/(2*kt_over_m))*vx) + 1)
vx = np.linspace(-5,5,100)
for kt_over_m in [1,2,5]:
plt.plot(vx,cdf(vx,kt_over_m),label=f'kT/m = {kt_over_m}')
plt.xlabel("V")
plt.ylabel("Prob (vx < V)")
plt.legend()
<matplotlib.legend.Legend at 0x7fe09d103a90>
How does this help us initialize our velocities?
If we generate a series of random numbers in the range \([0,1]\), we can map these to a series of velocity values using the CDF!
but first we need to solve this equation for \(V_{gen}\):
where \(\text{erfinv}\) is – you guessed it – the “inverse error function”.
from scipy.special import erfinv
import numpy as np
def init_vel(kt_over_m, shape):
R = np.random.random(shape) # generates a set of numbers uniformly in [0,1]
return np.sqrt(2*kt_over_m)*erfinv(2*R-1)
Let’s check a distribution of \(v_x\) values:
counts, bins = np.histogram(init_vel(1,(10000)), bins=50)
cbins = 0.5*(bins[1:] + bins[:-1]) # get bin midpoints
probs = counts / counts.sum() # normalize counts -> probabilities
mb_probs = mb_prob_x(cbins,1)
mb_probs /= mb_probs.sum() # normalize by sum for proper comparison
plt.plot(cbins, probs, label='Generated vx values')
plt.plot(cbins, mb_probs, label='Maxwell-Boltzmann probs')
plt.xlabel('vx')
plt.ylabel('Probability')
plt.legend()
<matplotlib.legend.Legend at 0x7fe09d627310>
Thermostats¶
In this course we will only be using Langevin dynamics which acts as its own thermostat. But here are some others:
Algorithm |
Description |
|---|---|
Nosé-Hoover |
adds extra degree of freedom (\(s\)) that scales all velocities; \(s\) has a mass and moves to keep simulation at desired \(T\) |
Andersen |
simulates “collisions” of randomly chosen particles and resamples their velocity according to Maxwell-Boltzmann |
Berendsen |
rescales all velocities at each step towards the target temperature; can result in the Flying Ice Cube Effect |
“Friends don’t let friends use Berendsen.” - @jchodera tweet
Force fields¶
Last week we used a simple potential energy function for our one-particle simulations. For biomolecular systems these need to contain all of the information about each atom (its charge, radius, etc) and how it interacts with other atoms in the system.
The interactions can be neatly separated into bonded and non-bonded interactions:
Reminder: the variable \({\bf x}\) holds all of the positions of all the atoms in the system. Today let’s consider \({\bf x}\) to be a 2D tensor (or, a 2D array) of size \((N,3)\).
The bonded energy function describes interactions that are due to covalent bonds:
The non-bonded energy function is composed of electrostatic and Lennard-Jones interactions:
\(U_{\text{elec}}\) is the interaction between charged groups from Coulomb’s Law:
where \(k_e\) is Coulomb’s constant, \(q_a\) is the “partial charge” on atom \(a\), and \(r_{ij}\) is the distance between atoms \(i\) and \(j\).
\(U_{\text{LJ}}\) is the interaction between particles that doesn’t account for their net charge. It does a couple of things:
prevents atomic electron clouds from overlapping in space (which is not allowed by the Pauli Exclusion principle)
allows for favorable “induced dipole” interactions (also referred to as “van der Waals” interactions, “London” forces, or “dispersion” forces)
The commonly used equation to describe these forces is the “Lennard-Jones” potential:
where \(\epsilon_{ij} = \sqrt{\epsilon_i \epsilon_j}\), and \(\sigma_{ij} = \frac{1}{2}(\sigma_i + \sigma_j)\).
\(\sigma_a\) is the van der Waals radius of atom \(a\) and \(\epsilon\) is the depth of its van der Waals energy well.
def ULJ(r_over_sigma, epsilon):
sigma_over_r = 1/r_over_sigma
sor6 = sigma_over_r**6
sor12 = sor6**2
return 4*epsilon*(sor12-sor6)
r_o_s = np.linspace(0.98,2.5,100)
plt.plot(r_o_s, ULJ(r_o_s, 0.5), label='epsilon = 0.5')
plt.plot(r_o_s, ULJ(r_o_s, 0.8), label='epsilon = 0.8')
plt.plot(r_o_s, ULJ(r_o_s, 1.0), label='epsilon = 1.0')
plt.xlabel("r/sigma")
plt.ylabel("ULJ(r)")
plt.legend()
<matplotlib.legend.Legend at 0x7fe09d57ae50>
Periodic boundaries¶
When conducting molecular simulations we can only afford to include so many particles. A common way to minimize the effects of a finite system size is to include periodic boundaries.
This essentially assumes that the simulation cell is surrounded in all directions by other identical simulation cells.
The implementation of periodic boundaries seems simple enough:
if you exit through the left side, you enter through the right
if you exit through the top, you enter through the bottom
etc.
But periodic boundaries can be tricky and are probably the cause of most errors and bugs in molecular dynamics simulations.
“Wrapping” coordinates¶
To apply those rules above and make sure particles do not leave the simulation box, you need to periodically “wrap” the coordinates back into the box.
Below is trajectory that leaves a simulation box, defined in two dimensions as \(x \in [-0.5,0.5]\), \(y \in [-0.5,0.5]\):
# an imaginary trajectory
traj = np.zeros((20,2))
traj[:,0] = np.linspace(0,0.7,20)
traj[:,1] = 0.2*np.sin(8*traj[:,0])
# used for visualizing the box
box_x = [-0.5,0.5,0.5,-0.5,-0.5]
box_y = [-0.5,-0.5,0.5,0.5,-0.5]
plt.figure(figsize=[5,5])
plt.plot(traj[:,0],traj[:,1],'bo')
plt.plot(box_x,box_y)
plt.xlim([-1,1])
plt.ylim([-1,1])
(-1.0, 1.0)
Here is a function that applies boundary conditions to a given trajectory frame and returns the “wrapped” positions back inside the box:
def apply_pbc(positions, box_size):
# positions and box size are both arrays with
# a dimensionality equal to the number of spatial dimensions
# i.e. if we are in 2D space, each array has two elements
half_box_size = box_size * 0.5
pbc_wrap_positions = np.zeros_like(positions)
old_positions = positions
# this while loop allows adjustments to be made until
# the positions of the particles stops changing
while np.any(pbc_wrap_positions != old_positions):
#
# Note on np.where syntax:
# np.where(CONDITION, value_if_true, value_if_false)
#
pbc_wrap_positions = np.where(positions > half_box_size,
positions - box_size,
positions)
pbc_wrap_positions = np.where(positions <= -half_box_size,
positions + box_size,
pbc_wrap_positions)
old_positions = pbc_wrap_positions
return pbc_wrap_positions
Let’s try this out by applying it to all of the frames in our trajectory:
new_frames = []
box_size = np.array([1,1])
for frame in traj:
new_frames.append(apply_pbc(frame,box_size))
new_traj = np.array(new_frames)
plt.figure(figsize=[5,5])
plt.plot(new_traj[:,0],new_traj[:,1],'bo')
plt.plot(box_x,box_y)
plt.xlim([-1,1])
plt.ylim([-1,1])
(-1.0, 1.0)
Note that we can use the same function with 3D data as well:
unwrapped_trajectory = np.random.random(size=(100,10,3))*10 # nframes, nparticles, ndim
print("Max x coordinate:",unwrapped_trajectory[:,:,0].max())
print("Max y coordinate:",unwrapped_trajectory[:,:,1].max())
print("Max z coordinate:",unwrapped_trajectory[:,:,2].max())
Max x coordinate: 9.997272451544207
Max y coordinate: 9.995211595558018
Max z coordinate: 9.997304830041116
box_size = np.array([5,5,5])
wrapped_trajectory = apply_pbc(unwrapped_trajectory,box_size)
print("Max x coordinate (after wrapping):",wrapped_trajectory[:,:,0].max())
print("Max y coordinate (after wrapping):",wrapped_trajectory[:,:,1].max())
print("Max z coordinate (after wrapping):",wrapped_trajectory[:,:,2].max())
Max x coordinate (after wrapping): 4.997272451544207
Max y coordinate (after wrapping): 4.995211595558018
Max z coordinate (after wrapping): 4.9973048300411165
Displacement vectors through periodic boundaries¶
Wrapping is important, but it is also important that inter-particle forces are computed properly through periodic boundaries. In other words, even though the positions of two particles are far apart, their “images” in adjacent cells could be close together. If these interactions aren’t taken into account, particles could experience discontinuities in their forces as they are wrapped through the periodic boundary. For instance, a particle would know that it is sitting almost on top of another particle until it crosses the boundary, at which point it suddenly experiences an ultra-high repulsive force, and then your system explodes!
To avoid this catastrophic scenario, we will make sure to compute forces using the closest possible displacement vector.
Consider two particles \(i\) and \(j\):
The difference in their \(x\) coordinates can be calculated three ways:
directly (\(x_j - x_i\))
using the “right hand image” (\(x_j - (x_i + L_x)\))
using the “left hand image” (\(x_j - (x_i - L_x)\))
where \(L_x\) is the length of the box in the \(x\) direction.
Clearly if we want \(i\) and \(j\) to “feel” each other across the periodic boundary, we need to use the difference with the smallest magnitude when we are calculating the forces.
The following function determines the shortest displacement vector between two atoms, even considering images that differ by more than one box length.
def shortest_disp_vector(posA, posB, box_size):
r = posA - posB # has one element for each spatial dimension
# loop over spatial dimensions
for i in range(len(posA)):
# if difference is greater than half the box length, adjust
if r[i] < -0.5*box_size[i]:
# determine number of adjustments needed
n = int((r[i] + 0.5*box_size[i])/box_size[i]) + 1
# make them
r[i] += n*box_size[i]
elif r[i] > 0.5*box_size[i]: # same thing in the negative direction
n = int((r[i] - 0.5*box_size[i])/box_size[i]) + 1
r[i] -= n*box_size[i]
return r
Now let’s test this out:
positions = np.array([[0.01,0.52,0.55],
[0.90,0.51,0.2]]) # two particles in 3D space
box_size = [1,1,1]
# compute direct displacement vector
direct_disp_vec = positions[0] - positions[1]
print("Direct displacement vector:",direct_disp_vec,"has magnitude:",np.sqrt(np.sum(np.square(direct_disp_vec))))
# compute displacement vector through the boundaries
closest_disp_vec = shortest_disp_vector(positions[0],positions[1],box_size)
print("Shortest displacement vector:",closest_disp_vec,"has magnitude:",np.sqrt(np.sum(np.square(closest_disp_vec))))
Direct displacement vector: [-0.89 0.01 0.35] has magnitude: 0.9563994981178107
Shortest displacement vector: [0.11 0.01 0.35] has magnitude: 0.3670149860700514
In this week’s lab we will be using a periodic box to simulate a set of Liquid Argon atoms.