Implementing the Euler integrator, and testing it against the Strmer/Verlet

In the lectures we also discussed, as an example of a poor algorithm for MD time integration, the Euler integrator. Let us now test its performance against the Strmer/Verlet scheme.

Task 2: Go through the below code cell, understand it, and then implement the Euler time integration algorithm by editing the (seven) lines marked by the comment #-- Edit this line.

Question 5: Note that here we take the initial velocities to be zero. Think: Does this correspond exactly to the same initial condition that we used in the Strmer/Verlet implemention above? Why / why not?

Answer 5: These two initial conditions are...??? Because...???

[ ]:

#-- initial positions:

x1 = 3.4

y1 = 2.1

z1 = 3.1

x2 = 2.7

y2 = 1.6

z2 = 3.7

#-- initial velocities:

vx1 = 0.0

vy1 = 0.0

vz1 = 0.0

vx2 = 0.0

vy2 = 0.0

vz2 = 0.0

#-- time integration details:

dt = 0.01 #-- time step

dt2 = dt**2 #-- time step squared

hdt2 = 0.5*dt2 #-- to save a multiplication in the innermost loop

startTime = 0 #-- from 

endTime = 1000 #-- until (in number of steps)

# INITIALIZE vectors for output:

tOutEU = []

rOutEU = []

#-- simulate:

for time in range(startTime, endTime, 1):

 #-- distance vector from the perspective of particle 1:

 x12 = x2 - x1

 y12 = y2 - y1

 z12 = z2 - z1

 #-- useful powers of the distance:

 r12squared = x12**2 + y12**2 + z12**2

 r12toSixth = r12squared**3

 r12toTwelfth = r12toSixth**2

 ljFactor = #-- Edit this line

 #-- force on particle 1:

 fx1 = x12*ljFactor

 fy1 = y12*ljFactor

 fz1 = z12*ljFactor

 #-- force on particle 2:

 fx2 = -fx1

 fy2 = -fy1

 fz2 = -fz1

 # MOVE PARTICLES using Euler time integration:

 #

 # calculate the New positions:

 x1 = #-- Edit this line

 y1 = #-- Edit this line

 z1 = #-- Edit this line

 x2 = #-- Edit this line

 y2 = #-- Edit this line

 z2 = #-- Edit this line

 # calculate the New velocities: 

 vx1 = vx1 + fx1*dt

 vy1 = vy1 + fy1*dt

 vz1 = vz1 + fz1*dt

 vx2 = vx2 + fx2*dt

 vy2 = vy2 + fy2*dt

 vz2 = vz2 + fz2*dt

 # save the time and the particle distance (so they can be plotted in the next cell):

 tOutEU.append(time*dt)

 r = np.sqrt(r12squared)

 rOutEU.append(r)

[ ]:

plt.plot(tOutEU,rOutEU, label='Euler')

#plt.plot(tOutSV,rOutSV, label='Strmer/Verlet')

plt.axis(xmin=0,xmax=endTime*dt)

plt.xlabel('time (LJ-units)')

plt.ylabel('distance between particles (LJ-units)')

plt.legend(loc='upper left')

plt.show()

Question 6: Run the simulation using the initial positions from Question 3: 1=3.4, 1=2.1, 1=3.1 for the first particle and 2=2.7, 2=1.6, 2=3.7 for the second particle. Interpret again the distance-vs.-time plot: How does the behavior of the particles differ from what you saw when you started from the same positions with Strmer/Verlet in Question 3? Why?

Answer 6: Now the particles...??? Because...???