Introduction to Molecular Dynamics¶
In this notebook you will use a Verlet scheme to simulate the dynamics of a 1-D Harmonic Oscillator and 1-D double well potential
#setup the notebook
%pylab inline
import numpy as np
Part 1, set up the potential and plot it¶
This function has to return the energy and force for a 1d harmonic potential. The potential is $U(x) = \frac{1}{2} k (x - x_0)^2$ and $F = -\frac{dU(x)}{dx}$
#this function returns the energy and force on a particle from a harmonic potential
def harmonic_oscillator_energy_force(x,k=1,x0=0):
#calculate the energy on force on the right hand side of the equal signs
energy = 0.5*k*(x-x0)**2
force = -k*(x-x0)
return energy, force
#this function will plot the energy and force
#it is very general since it uses a special python trick of taking arbitrary named arguments (**kwargs)
#and passes them on to a specified input function
def plot_energy_force(function, xmin=-3,xmax=3,spacing=0.1,**kwargs):
x_points = np.arange(xmin,xmax+spacing,spacing)
energies, forces = function(x_points,**kwargs)
label = 'U(x)'
for arg in kwargs:
label=label+', %s=%s'%(arg,str(kwargs[arg]))
p = plt.plot(x_points,energies,label=label)
plt.plot(x_points,forces,label='',color=p[0].get_color(),linestyle='--')
plt.legend(loc=0)
#we can plot the energy (solid) and forces (dashed) to see if you did it right.
#Is this what you expect???
#Play around with changing inputs to the energy, k and x0
plot_energy_force(harmonic_oscillator_energy_force,k=1)
plot_energy_force(harmonic_oscillator_energy_force,k=5)
plot_energy_force(harmonic_oscillator_energy_force,k=10)
plot_energy_force(harmonic_oscillator_energy_force,k=5,x0=1)
Part 2, code the exact solutions for the 1-d harmonic oscillator¶
When coding a new algorithm, like Molecular Dynamics, it is always good to start with an exactly solvable problem. We want to solve
$F = m a = m \frac{d^2 x}{dt^2}$, where a is the second time derivative of position, and $F = -\frac{dU(x)}{dx}$. For the harmonic oscillator potential given above, there is an exact solution.
It is:
$x(t) = A \cos(\omega t + \phi)$
$v(t) = \frac{d x(t)}{d t} = -A \omega \sin(\omega t + \phi)$
where $\omega=\sqrt{k/m}$
$\phi$ and $A$ are set by the initial conditions of your system. What we do is choose an initial energy $E$, which we expect to be conserved at all times. We need to know the maximum $x$, which comes when all energy is potential energy, and when $\cos=1$ (so $x_{max}=A$). Therefore $\frac{1}{2} k x_{max}^2 = \frac{1}{2} k A^2 = E$. Hence, $A=\sqrt{2E/k}$. $\phi$ depends on where the particle starts, since $x(0) = A \cos(\phi)$.
Check for yourself that $\frac{1}{2} m v_{max}^2 = E$ as well.
Below, you will write this into python code and then it will be plotted.
#fill in the formulas below. remember, to get trig functions you have to do things like np.cos()
def harmonic_oscillator_position_velocity(t, A=1, omega=1, phi=0):
position = A*np.cos(omega*t+phi)
velocity = -A*omega*np.sin(omega*t+phi)
return position, velocity
#this function will plot the energy and force for a harmonic oscilator
def plot_harmonic_oscillator(t_max=10,dt=0.1,**kwargs):
t_points = np.arange(0,t_max+dt,dt)
position, velocity = harmonic_oscillator_position_velocity(t_points,**kwargs)
label = 'x(t)'
for arg in kwargs:
label=label+', %s=%s'%(arg,str(kwargs[arg]))
p = plt.plot(t_points,position,label=label)
plt.plot(t_points,velocity,label='',color=p[0].get_color(),linestyle='--')
plt.legend(loc='upper right')
#if it is working, you will get a plot with position as a solid line and velocity as dashed
plot_harmonic_oscillator()
Look at the effect of changing harmonic oscillator parameters¶
Notice the effect of changing, $A$, $\omega$, and $\phi$.
Things to look at and check for yourself:
- Position and velocity are out of phase
- Changing the amplitude increases both the max speed and position
- Changing the frequency changes how often the peaks appear, and also the max velocity
- Changing the phase does not change much about the solution, just where the peaks appear
plot_harmonic_oscillator()
plot_harmonic_oscillator(A=2)
plt.title('Changing amplitude')
plt.figure()
plot_harmonic_oscillator()
plot_harmonic_oscillator(omega=2)
plt.title('Changing frequency')
plt.figure()
plot_harmonic_oscillator()
plot_harmonic_oscillator(phi=1)
plt.title('Changing phase')
Part 3, code velocity verlet algorithm¶
Now you will code the velocity verlet algorithm. You will know it is correct because you will get curves that overlap with the exact solutions, but at discrete times, with small intervals $dt$.
The following equations are alternated to move time forward
$x(t+dt) = x(t) + dt \cdot v(t) + \frac{dt^2}{2 m} F(t)$
$v(t+dt) = v(t) + \frac{dt}{2m} (F(t+dt) + F(t))$
Note that you have to recompute the forces after each position update.
In the following, you can set the mass $m=1$
def position_update(x,v,F,dt):
x_new = x + dt*v + 0.5*dt*dt*F
return x_new
def velocity_update(v,F_new,F_old,dt):
v_new = v + 0.5*dt*(F_old+F_new)
return v_new
#this function will take the initial energy as an input and run velocity verlet dynamics
def velocity_verlet_harmonic_oscillator(potential, max_time, dt, initial_position, initial_velocity,
save_frequency=3, **kwargs ):
x = initial_position
v = initial_velocity
t = 0
step_number = 0
positions = []
velocities = []
total_energies = []
save_times = []
while(t<max_time):
potential_energy, force = potential(x,**kwargs)
if step_number%save_frequency == 0:
e_total = .5*v*v + potential_energy
positions.append(x)
velocities.append(v)
total_energies.append(e_total)
save_times.append(t)
# update the positions
x = position_update(x,v,force,dt)
# recompute the energy and force
potential_energy2, force2 = potential(x,**kwargs)
# update the velocity
v = velocity_update(v,force2,force,dt)
t = t+dt
step_number = step_number + 1
return save_times, positions, velocities, total_energies
Now we will run our MD simulation of a harmonic oscillator¶
You can choose the initial energy and spring constant that you want and see how the results look. If it's working, then the energy should be conserved to much better than 1% and your simulation data (circles) should fit on the exact solution (lines)
initial_energy = 2
my_k = 3
my_max_time = 10
#set the initial conditions. this is easiest if you set x=x_max and v=0 based on the input energy E.
# Otherwise you will have to include a phase factor to show that you got the exact solution
initial_position = np.sqrt(2*initial_energy/my_k)
initial_velocity = 0
#let's set a good timestep based on dt = tau/100 where tau=2 pi/ omega
my_omega = np.sqrt(my_k)
tau = 2*np.pi/my_omega
my_dt = tau/100.
times, positions, velocities, total_energies = velocity_verlet_harmonic_oscillator(harmonic_oscillator_energy_force,
my_max_time, my_dt, \
initial_position, initial_velocity,\
k=my_k)
# What is the A value prefactor for the harmonic oscillator? See equations above
my_A = np.sqrt(initial_energy*2/my_k)
plt.plot(times,positions,marker='o',label='sim. position',linestyle='')
plt.plot(times,velocities,marker='s',label='sim. velocity',linestyle='')
plot_harmonic_oscillator(t_max=my_max_time,omega=my_omega, A=my_A)
xlabel('time')
legend(loc='upper center')
plt.figure()
plt.plot(times,total_energies,marker='o',linestyle='',label='Simulated E')
plt.axhline(initial_energy,color='black',label='Exact')
plt.ylim(0.95*initial_energy,1.05*initial_energy)
xlabel('time')
ylabel("Total Energy")
legend()
Histogram Position and Velocity¶
What is the probability of seeing a given position or velocity? Do the histograms look like what you expect?
def bin_centers(bin_edges):
return (bin_edges[1:]+bin_edges[:-1])/2.
#to get a good histogram, we need to run a lot longer than before
my_max_time = 1000
times, positions, velocities, total_energies = velocity_verlet_harmonic_oscillator(harmonic_oscillator_energy_force,
my_max_time, my_dt, \
initial_position, initial_velocity,\
k=my_k)
dist_hist, dist_bin_edges = np.histogram(positions,bins=20,density=True)
vel_hist, vel_bin_edges = np.histogram(velocities,bins=20,density=True)
plot(bin_centers(dist_bin_edges), dist_hist,marker='o',label='P(x)')
plot(bin_centers(vel_bin_edges), vel_hist,marker='s',label='P(v)')
legend(loc='upper center')
Part 4, simulate double well potential¶
For the harmonic oscillator, we knew the exact solution. Let's do another similar problem for practice.
One useful model potential is the double-well:
$U(x) = \frac{k}{4} (x-a)^2 (x+a)^2$
This potential has a minimum at $x=a$ and $x=-a$. It also has a barrier at $x=0$.
What is the height of this barrier?
Fill in the potential and force functions below and it will be plotted.
Is the height of the barrier what you predicted? Try changing $k$ and $a$ to see the effect.
#this function returns the energy and force on a particle from a double well
def double_well_energy_force(x,k,a):
#calculate the energy on force on the right hand side of the equal signs
energy = 0.25*k*((x-a)**2)*((x+a)**2)
force = -k*x*(x-a)*(x+a)
return energy, force
plot_energy_force(double_well_energy_force, xmin=-4,xmax=+4, k=1, a=2)
axhline(0,linestyle='--',color='black')
axvline(0,linestyle='--',color='black')
ylim(-10,10)
Run velocity verlet dynamics on the double well¶
Set your $k$ and $a$ values. Pick an initial position, velocity, and timestep. Then run the dynamics.
Try changing these parameters. What is your initial energy? What happens if the initial energy is less than the barrier height?
# Here the initial energy is 8 and the barrier height is 4, so it goes back and forth
my_k = 1
my_a = 2
my_initial_position = -2
my_initial_velocity = 4
my_dt = 0.001
my_max_time = 10
times, positions, velocities, total_energies = velocity_verlet_harmonic_oscillator(double_well_energy_force,
my_max_time, my_dt, \
my_initial_position, my_initial_velocity,\
k=my_k, a=my_a)
plt.plot(times,positions,marker='o',label='position',linestyle='')
plt.plot(times,velocities,marker='s',label='velocity',linestyle='')
xlabel('time')
legend(loc='upper center')
plt.figure()
initial_energy = total_energies[0]
plt.plot(times,total_energies,marker='o',linestyle='',label='Simulated E')
plt.axhline(initial_energy,color='black',label='Exact')
plt.ylim(0.95*initial_energy,1.05*initial_energy)
xlabel('time')
ylabel("Total Energy")
legend()
plt.figure()
#to get a good histogram, we need to run a lot longer than before
my_max_time = 1000
times, positions, velocities, total_energies = velocity_verlet_harmonic_oscillator(double_well_energy_force,
my_max_time, my_dt, \
my_initial_position, my_initial_velocity,\
k=my_k, a=my_a)
dist_hist, dist_bin_edges = np.histogram(positions,bins=25,density=True)
vel_hist, vel_bin_edges = np.histogram(velocities,bins=25,density=True)
plot(bin_centers(dist_bin_edges), dist_hist,marker='o',label='P(x)')
plot(bin_centers(vel_bin_edges), vel_hist,marker='s',label='P(v)')
legend(loc='upper center')
xlim(-4,4)
# here the initial energy is much less, hence the particle gets trapped (see the histogram).
#Then it's dynamics are like a simple harmonic oscillator centered at, in this case, the left basin center
my_k = 1
my_a = 2
my_initial_position = -2
my_initial_velocity = 1.5
my_dt = 0.001
my_max_time = 10
times, positions, velocities, total_energies = velocity_verlet_harmonic_oscillator(double_well_energy_force,
my_max_time, my_dt, \
my_initial_position, my_initial_velocity,\
k=my_k, a=my_a)
plt.plot(times,positions,marker='o',label='position',linestyle='')
plt.plot(times,velocities,marker='s',label='velocity',linestyle='')
xlabel('time')
legend(loc='upper center')
plt.figure()
initial_energy = total_energies[0]
plt.plot(times,total_energies,marker='o',linestyle='',label='Simulated E')
plt.axhline(initial_energy,color='black',label='Exact')
plt.ylim(0.95*initial_energy,1.05*initial_energy)
xlabel('time')
ylabel("Total Energy")
legend()
plt.figure()
#to get a good histogram, we need to run a lot longer than before
my_max_time = 1000
times, positions, velocities, total_energies = velocity_verlet_harmonic_oscillator(double_well_energy_force,
my_max_time, my_dt, \
my_initial_position, my_initial_velocity,\
k=my_k, a=my_a)
dist_hist, dist_bin_edges = np.histogram(positions,bins=25,density=True)
vel_hist, vel_bin_edges = np.histogram(velocities,bins=25,density=True)
plot(bin_centers(dist_bin_edges), dist_hist,marker='o',label='P(x)')
plot(bin_centers(vel_bin_edges), vel_hist,marker='s',label='P(v)')
legend(loc='upper center')
xlim(-4,4)
