比较分子动力学模拟中Verlet算法、Speed Verlet算法及辛算法的精度
原理
分子动力学模拟是一种研究凝聚态系统的有力方法,对于经典体系,分子动力学假定所有粒子的运动服从牛顿运动方程,只要给出粒子的初始位置、速度(初始条件),边界条件,以及粒子之间相互作用的势函数,就能得到系统在未来任意时刻的状态信息。
由于多原子体系的相互作用非常复杂,难以用解析法求解体系的运动方程,因此用数值方法求解运动方程的近似解。这个过程就不可避免地引入了误差,包括截断误差和舍入误差。前者与不同算法有关,如有限差分法采用TaylorTaylorTaylor展开后截断某些高阶项会引入截断误差,后者来源于计算机本身的精度。两种误差都能通过减小时间步长δt\delta tδt来减小,当δt\delta tδt较大时,截断误差起主要作用,但随着δt\delta tδt的减小会迅速降低,舍入误差随δt\delta tδt变化不大,并在δt\delta tδt较小时起主要作用。
进行数值求解的算法有很多种,对于差分法,其共同点是将粒子的位置,速度,加速度展开为TaylorTaylorTaylor级数。下面我们来分析常用的VerletVerletVerlet算法和更精确一些的SpeedSpeedSpeed VerletVerletVerlet算法的具体步骤以及它们的截断误差。
VerletVerletVerlet AlgorithmAlgorithmAlgorithm
VerletVerletVerlet算法在分子动力学模拟中应用极为广泛,而且原理也很简明。这种方法利用粒子在ttt时刻的位置和加速度,以及t−δtt-\delta tt−δt时刻的位置,计算出t+δtt+\delta tt+δt时刻的位置。
将粒子的位置用TaylorTaylorTaylor公式展开,忽略四阶以上小量:
r(t+δt)=r(t)+ddtr(t)δt+12!d2dt2r(t)(δt)2+13!d3dt3r(t)(δt)3+o[(δt)4]r(t+\delta t)=r(t)+\frac{d}{dt} r(t) \delta t+\frac{1}{2!} \frac{d^2}{dt^2}r(t)(\delta t)^2+\frac{1}{3!} \frac{d^3}{dt^3}r(t)(\delta t)^3+o[(\delta t)^4]r(t+δt)=r(t)+dtdr(t)δt+2!1dt2d2r(t)(δt)2+3!1dt3d3r(t)(δt)3+o[(δt)4]
将式中的δt\delta tδt换为−δt-\delta t−δt得:
r(t−δt)=r(t)−ddtr(t)δt+12!d2dt2r(t)(δt)2−13!d3dt3r(t)(δt)3+o[(δt)4]r(t-\delta t)=r(t)-\frac{d}{dt} r(t) \delta t+\frac{1}{2!} \frac{d^2}{dt^2}r(t)(\delta t)^2-\frac{1}{3!} \frac{d^3}{dt^3}r(t)(\delta t)^3+o[(\delta t)^4]r(t−δt)=r(t)−dtdr(t)δt+2!1dt2d2r(t)(δt)2−3!1dt3d3r(t)(δt)3+o[(δt)4]
以上两式相加,并忽略o[(δt)4]o[(\delta t)^4]o[(δt)4]小量得:
r(t+δt)=2r(t)−r(t−δt)+d2dt2r(t)(δt)2r(t+\delta t)=2r(t)-r(t-\delta t)+\frac{d^2}{dt^2} r(t) (\delta t)^2r(t+δt)=2r(t)−r(t−δt)+dt2d2r(t)(δt)2
两式相减,并忽略o[(δt)2]o[(\delta t)^2]o[(δt)2]小量得:
v(t)=drdt=12δt[r(t+δt)−r(t−δt)]v(t)=\frac{dr}{dt}=\frac{1}{2\delta t}[r(t+\delta t)-r(t-\delta t)]v(t)=dtdr=2δt1[r(t+δt)−r(t−δt)]
至此,我们得到了t+δtt+\delta tt+δt时刻粒子位置和ttt时刻速度的表达式,并且知道了位置的误差正比于o[(δt)4]o[(\delta t)^4]o[(δt)4],速度的误差正比于o[(δt)2]o[(\delta t)^2]o[(δt)2]。VerletVerletVerlet算法的数值稳定性比简单的欧拉方法高很多,并保持了物理系统中的时间可逆性与相空间体积元体积守恒的性质。然而,这种算法得出的轨迹与速度无关(无法与热浴耦联),并且不能给出同一时刻的位置和速度。
SpeedSpeedSpeed VerletVerletVerlet AlgorithmAlgorithmAlgorithm
在VerletVerletVerlet AlgorithmAlgorithmAlgorithm中,我们发现其无法给出同一时刻的速度和位置,SpeedSpeedSpeed VerletVerletVerlet(有些地方也叫VelocityVelocityVelocity VerletVerletVerlet)正是用来克服这个缺陷的。
将位置展开到二阶小量:
r(t+δt)=r(t)+ddtr(t)δt+12!d2dt2r(t)(δt)2r(t+\delta t)=r(t)+\frac{d}{dt} r(t) \delta t+\frac{1}{2!} \frac{d^2}{dt^2}r(t)(\delta t)^2r(t+δt)=r(t)+dtdr(t)δt+2!1dt2d2r(t)(δt)2
引入t+δt2t+\frac{\delta t}{2}t+2δt时刻的速度:
v(t+δt2)=v(t)+12ddtv(t)δtv(t+\frac{\delta t}{2})=v(t)+\frac{1}{2}\frac{d}{dt} v(t) \delta tv(t+2δt)=v(t)+21dtdv(t)δt
则t+δtt+\delta tt+δt时刻的速度可以表示为
v(t+δt)=v(t+δt2)+12ddtv(t+δt2)δtv(t+\delta t)=v(t+\frac{\delta t}{2})+\frac{1}{2}\frac{d}{dt} v(t+\frac{\delta t}{2}) \delta tv(t+δt)=v(t+2δt)+21dtdv(t+2δt)δt
因此,t+δtt+\delta tt+δt时刻的速度又可表示为:
v(t+δt)=v(t)+12ddtv(t)δt+12ddtv(t+δt2)δtv(t+\delta t)=v(t)+\frac{1}{2}\frac{d}{dt} v(t) \delta t+\frac{1}{2}\frac{d}{dt} v(t+\frac{\delta t}{2}) \delta tv(t+δt)=v(t)+21dtdv(t)δt+21dtdv(t+2δt)δt
其中加速度项通过a(t)=ddtv(t)=−1m∇V(r(t))a(t)=\frac{d}{dt} v(t)=-\frac{1}{m}\nabla V(r(t))a(t)=dtdv(t)=−m1∇V(r(t))得出。因此,SpeedSpeedSpeed VerletVerletVerlet能够给出同一时刻的速度和位置。
SymplecticSymplecticSymplectic AlgorithmAlgorithmAlgorithm
经典动力学方程的求解算法主要包括各级各阶的RK算法,中心差分法,WilsonWilsonWilson θ\thetaθ方法和NewmarkNewmarkNewmark方法等。然而,这些算法本身耗散系统的能量,并使得动态响应的相位滞后,因此其长期跟踪能力有限。
冯康等[1][1][1]建立了哈密顿系统的辛几何算法,从理论上阐明了传统数值方法的耗散来自于其截断项,而辛算法对应的差分格式严格保持哈密顿系统的辛结构,有限阶辛算法的截断部分不会导致系统能量发生线性变化。
辛格式有许多种。Kane等证明了隐式中点格式和NewmarkNewmarkNewmark格式都是变分积分格式,对某些参数的NewmarkNewmarkNewmark格式为辛形式。邢誉峰[2][2][2]等证明了δ=0.5\delta=0.5δ=0.5和α=0.25\alpha=0.25α=0.25的NewmarkNewmarkNewmark算法就是EulerEulerEuler中点隐式差分格式。
由于哈密顿方程可写为q.=Tp\stackrel{.}{q}=Tpq.=Tpp.=−Vq\stackrel{.}{p}=-Vqp.=−Vq对其进行差分可得到1h(qm+3/2−qm+1/2)=TPm+1\frac{1}{h}\Big(q_{m+3/2}-q_{m+1/2}\Big)=TP_{m+1}h1(qm+3/2−qm+1/2)=TPm+11h(pm+1−qm)=−Vqm+1/2\frac{1}{h}\Big(p_{m+1}-q_{m}\Big)=-Vq_{m+1/2}h1(pm+1−qm)=−Vqm+1/2可以证明,上述表示是一种辛格式,即可分线性哈密顿体系的交叉显式辛格式。选用上式进行模拟。
设计思路
一维双谐振子体系
首先通过对一维双谐振子体系的模拟,比较VerletVerletVerlet AlgorithmAlgorithmAlgorithm及SpeedSpeedSpeed VerletVerletVerlet AlgorithmAlgorithmAlgorithm及辛算法的速度及位置精度。
对单个一维线性谐振子,有V(x)=12kx2V(x)=\frac{1}{2}kx^2V(x)=21kx2
md2xdt2=−kxm\frac{d^2x}{dt^2}=-kxmdt2d2x=−kx
对于两个谐振子的情形,如下图所示,给予它们共线的初速度,考虑它们运动的情形
假设两球质量mA=mB=1kgm_A=m_B=1kgmA=mB=1kg,弹簧劲度系数k=0.5N/mk=0.5N/mk=0.5N/m,初始位置xA=−5m,xB=5mx_A=-5m,x_B=5mxA=−5m,xB=5m,初始速度xA=1m/s,xB=−1m/sx_A=1m/s,x_B=-1m/sxA=1m/s,xB=−1m/s,初始加速度aA=aB=0m/s2a_A=a_B=0 m/s^2aA=aB=0m/s2,时间步长δt=0.01\delta t=0.01δt=0.01,步数N=1000−10000N=1000-10000N=1000−10000(谐振子周期在此条件下为T=2πT=2\piT=2π,时间步长应取特征时间的10−310^{-3}10−3~10−210^{-2}10−2,因此取δt=0.01\delta t=0.01δt=0.01是在范围内的)。
运算结果
对于双振子中的一个,比较VerletVerletVerlet AlgorithmAlgorithmAlgorithm,SpeedSpeedSpeed VerletVerletVerlet AlgorithmAlgorithmAlgorithm,辛算法与精确解的轨迹图有(pythonpythonpython实现):
将它们放在一起比较:
放大局部区域可见
可见VerletVerletVerlet AlgorithmAlgorithmAlgorithm与SpeedSpeedSpeed VerletVerletVerlet AlgorithmAlgorithmAlgorithm以及SymplecticSymplecticSymplectic AlgorithmAlgorithmAlgorithm与精确解的轨迹都重合得很好,这是可以预期的,因为前两种算法位置的误差都是正比于o[(δt)4]o[(\delta t)^4]o[(δt)4]的,而辛算法的耗散在理论上应该非常小。
为了比较三种算法速度的精度,作x−px-px−p相空间轨迹图如下:
放大局部区域可见
可见SpeedSpeedSpeed VerletVerletVerlet AlgorithmAlgorithmAlgorithm以及辛算法与精确解的相空间轨迹图重合得很好。
为了更细致地比较,列出N取不同值时各个算法得到的最后时刻的位置表(δt=0.01\delta t =0.01δt=0.01,保留8位小数):
最后时刻的速度表:
二维三谐振子体系
除一维双谐振子体系之外,还进行了二维三谐振子体系的模拟。但由于这种体系的精确解难以求解,就没有进行几种算法与精确解的对照,只采用VerletVerletVerlet AlgorithmAlgorithmAlgorithm以及SpeedSpeedSpeed VerletVerletVerlet AlgorithmAlgorithmAlgorithm进行了模拟。
假设三球质量mA=mB=mC=1kgm_A=m_B=m_C=1kgmA=mB=mC=1kg,弹簧劲度系数k=1N/mk=1N/mk=1N/m,
初始位置(xA,yA)=(5,0.5)m,(xB,yB)=(−5,0)m,(xC,yC)=(0,1)m(x_A,y_A)=(5,0.5)m,(x_B,y_B)=(-5,0)m,(x_C,y_C)=(0,1)m(xA,yA)=(5,0.5)m,(xB,yB)=(−5,0)m,(xC,yC)=(0,1)m,
初始速度(xA,yA)=(−1,0)m/s,(xB,yB)=(0,0)m/s,(xC,yC)=(0.5,−0.5)m/s(x_A,y_A)=(-1,0)m/s,(x_B,y_B)=(0,0)m/s,(x_C,y_C)=(0.5,-0.5)m/s(xA,yA)=(−1,0)m/s,(xB,yB)=(0,0)m/s,(xC,yC)=(0.5,−0.5)m/s,
初始加速度全为零,时间步长δt=0.01\delta t=0.01δt=0.01,步数N=100N=100N=100。这里取较小的步数是为了得到清晰易于观察的图像。
下面三幅图从上到下依次是VerletVerletVerlet AlgorithmAlgorithmAlgorithm与SpeedSpeedSpeed VerletVerletVerlet AlgorithmAlgorithmAlgorithm模拟三谐振子体系的x−t,y−t,x−yx-t,y-t,x-yx−t,y−t,x−y图。
下图是两种算法模拟三谐振子体系的相空间轨迹图,选取了其中一个谐振子的x−px-px−p轨迹。
可见在上述模拟参数下,两种算法得到的位置和速度没有明显区别。
为了进一步验证VerletVerletVerlet AlgorithmAlgorithmAlgorithm与SpeedSpeedSpeed VerletVerletVerlet AlgorithmAlgorithmAlgorithm得到的速度和位置的一致性,取N=1000N=1000N=1000时最后一帧的速度和位置如下表(δt=0.01,\delta t =0.01,δt=0.01,保留8位小数):
小结
通过对一维双谐振子体系的模拟,比较了VerletVerletVerlet AlgorithmAlgorithmAlgorithm,SpeedSpeedSpeed VerletVerletVerlet AlgorithmAlgorithmAlgorithm以及辛算法的精度,三者和精确解重合度都很高,而且在速度和位置上的精度在测试的模拟参数下没有明显的区别。辛算法由于具有低耗散的优势,其优势在长时段的模拟中应该能够体现出来。
通过对二维三谐振子体系的模拟,验证了VerletVerletVerlet AlgorithmAlgorithmAlgorithm与SpeedSpeedSpeed VerletVerletVerlet AlgorithmAlgorithmAlgorithm得到的速度及位置是相当一致的。
参考文献
[1]冯康,秦孟兆.哈密顿系统的辛几何算法.[M].杭州:浙江科学技术出版社,
[2]邢誉峰,杨蓉.动力学平衡方程的Euler中点辛差分求解格式[J].力学学报,(01):100-105.
python 代码
一维双谐振子体系
from pylab import *import osdt = 0.01N=1000x0=5v0=1a0=0m1=1m2=1k=0.5# 用verlet算法模拟谐振子的运动轨迹def trj_v():t = [0,dt]a1 = [a0]a2 = [a0]vx1 = [-v0]vx2 = [v0]x1=[x0,x0+vx1[0]*dt+ k/m1*a1[0]*1/2*dt**2]x2=[-x0,-x0+vx2[0]*dt+ k/m2*a2[0]*1/2*dt** 2]while(t[-1]<N):a1.append(-k / m1 * (x1[-1] - x2[-1] - (x1[0] - x2[0])))a2.append(-k / m2 * (x2[-1] - x1[-1] - (x2[0] - x1[0])))x1.append(2*x1[-1] -x1[-2] + ( a1[-1]*dt**2))x2.append(2*x2[-1] -x2[-2] + (a2[-1]*dt**2))vx1.append((x1[-1]-x1[-3])/(2*dt))vx2.append((x2[-1]-x2[-3])/(2*dt))t.append(t[-1]+dt)return x1[:-1],x2[:-1],t[:-1],vx1,vx2#删去最后一个元素,便于和speed verlet的结果相比较#输出坐标文件x1_v=trj_v()[0]vx1=trj_v()[3]t_v=trj_v()[2]coordinate = open("verlet", 'w')for i in range(len(t_v)):coordinate.write(str(t_v[i]) + "\t" + str(x1_v[i]) + "\t" + str(vx1[i]) + "\n")coordinate.close()#画verlet算法相应的轨迹图def plot_v_trj():plot(t_v,x1_v,'yellowgreen',label='verlet algorithm')grid(True)#画verlet算法的相空间轨迹图def plot_v_xp():x1 = trj_v()[0]p1 = m1*trj_v()[3]plot(x1, p1, 'yellowgreen',label='verlet algorithm',marker="<")grid(True)######################精确解的轨迹def plot_exact():t=trj_v()[2]x=-sin(array(t))+5 #适用于以上参数情形v = -cos(array(t))plot(t,x,'orange',label='exact solution')return x,t,vdef plot_exact_xp():t = trj_v()[2]x = -sin(array(t))+ 5p = -cos(array(t))plot(x, p, 'orange', label='exact solution',marker="x")t_ex=plot_exact()[1]x_ex=plot_exact()[0]v_ex=plot_exact()[2]#输出坐标文件coordinate = open('Exact Solution', 'w')for i in range(len(t_ex)):coordinate.write(str(t_ex[i]) + "\t" + str(x_ex[i]) + "\t" + str(v_ex[i]) + "\n")coordinate.close()###################### 用 speed verlet算法模拟谐振子的运动轨迹def trj_sv():t = [0]x1=[x0]x2=[-x0]vx1 = [-v0]vx2=[v0]a1=[a0]a2=[a0]k=0.5while(t[-1]<N):x1.append(x1[-1] + vx1[-1]*dt + a1[-1]*1/2*dt**2)x2.append(x2[-1] + vx2[-1]*dt + a2[-1]*1/2*dt**2)a1.append(-k / m1 * (x1[-1] - x2[-1]-(x1[0]-x2[0])))a2.append(-k / m2 * (x2[-1] - x1[-1]-(x2[0]-x1[0])))# 定义加速度时要小心vx1.append(vx1[-1] + 1 / 2 * dt * (a1[-1] + a1[-2]))vx2.append(vx2[-1] + 1 / 2 * dt * (a2[-1] + a2[-2]))t.append(t[-1]+dt)return x1,x2,t,vx1,vx2#画speed verlet算法相应的轨迹图def plot_sv_trj():plot(t,x1,'steelblue',label='speed verlet algorithm')grid(True)#画speed_verlet算法的相空间轨迹图def plot_sv_xp():x1 = trj_sv()[0]p1 = trj_sv()[3]plot(x1, p1, 'r', label='speed verlet algorithm',marker=">")grid(True)#输出坐标文件x1=trj_sv()[0]x2=trj_sv()[1]t=trj_sv()[2]vx1=trj_sv()[3]coordinate = open("speed_verlet", 'w')for i in range(len(t)):coordinate.write(str(t[i]) + "\t" + str(x1[i]) + "\t" + str(vx1[i]) + "\n")coordinate.close()###################辛算法,可分线性哈密顿体系的交叉显式辛格式def trj_e2():t = [0]vx1 = [-v0]vx2 = [v0]x1 = [x0]x2 = [-x0]#w = 2 ** (1 / 2)k=1while (t[-1] < N):x1.append(x1[-1] +dt*vx1[-1])x2.append(x2[-1] +dt*vx2[-1])# a1.append(-k / m1 *(x1[-1] - x2[-1] - (x1[0] - x2[0])))# a2.append(-k / m2 * (x2[-1] - x1[-1] - (x2[0] - x1[0])))vx1.append(vx1[-1]-dt*k*(x1[-1]-x1[0])/m1)vx2.append(vx1[-1]-dt*k*(x2[-1]-x2[0])/m2)t.append(t[-1] + dt)#print(x1)return x1, x2, t, vx1, vx2def plot_e2_trj():plot(t_e2,x1_e2,'r',label='Sympletic Algorithm')grid(True)def plot_e2_xp():x1 = trj_e2()[0]p1 = m1*trj_e2()[3]plot(x1, p1, 'g', label='Sympletic Algorithm',marker="*")grid(True)#输出坐标文件x1_e2=trj_e2()[0]x2_e2=trj_e2()[1]t_e2=trj_e2()[2]vx1_e2=trj_e2()[3]coordinate = open('Sympletic Algorithm', 'w')for i in range(len(t_e2)):coordinate.write(str(t_e2[i]) + "\t" + str(x1_e2[i]) + "\t" + str(vx1_e2[i]) + "\n")coordinate.close()################################ 使得图像可以随鼠标滚轮而缩放#思路:#使用fig.canvas.mpl_connect()函数来绑定相关fig的滚轮事件#利用事件event的inaxes属性获取当前鼠标所在坐标系ax#使用get_xlim()函数获取坐标系ax的x/y轴坐标刻度范围#使用set()函数对坐标系ax进行放大/缩小import numpy as npfig = figure()def call_back(event):axtemp = event.inaxesx_min, x_max = axtemp.get_xlim()increment = (x_max - x_min) / 10if event.button == 'up':axtemp.set(xlim=(x_min + increment, x_max - increment))print('up')elif event.button == 'down':axtemp.set(xlim=(x_min - increment, x_max + increment))print('down')fig.canvas.draw_idle() # 绘图动作实时反映在图像上fig.canvas.mpl_connect('scroll_event', call_back)fig.canvas.mpl_connect('button_press_event', call_back)# plot_v_trj()# plot_sv_trj()# plot_e2_trj()# plot_exact()# legend(loc="right")# xlabel('t',fontsize="16")# ylabel('x',fontsize="16")# title('Twins Harmonic Oscillators',fontsize="16")# show()subplot(4,1,1)plot_v_trj()title('Twins Harmonic Oscillators',fontsize="16")legend(loc="right")ylabel('x',fontsize="16")subplot(4,1,2)plot_sv_trj()ylabel('x',fontsize="16")legend(loc="right")subplot(4,1,3)plot_e2_trj()legend(loc="right")ylabel('x',fontsize="16")subplot(4,1,4)plot_exact()legend(loc="right")xlabel('t',fontsize="16")ylabel('x',fontsize="16")show()fig = figure()def call_back(event):axtemp = event.inaxesx_min, x_max = axtemp.get_xlim()increment = (x_max - x_min) / 10if event.button == 'up':axtemp.set(xlim=(x_min + increment, x_max - increment))print('up')elif event.button == 'down':axtemp.set(xlim=(x_min - increment, x_max + increment))print('down')fig.canvas.draw_idle() # 绘图动作实时反映在图像上fig.canvas.mpl_connect('scroll_event', call_back)fig.canvas.mpl_connect('button_press_event', call_back)plot_v_xp()plot_sv_xp()plot_e2_xp()plot_exact_xp()legend(loc="right")xlabel('x',fontsize="16")ylabel('p',fontsize="16")title('Twins Harmonic Oscillators',fontsize="16")show()
二维三谐振子体系
from matplotlib.pyplot import *from vpython import *class Harmonics:def __init__(self,_m1=1,_x1=5,_y1=0.5,_vx1=-1,_vy1=0,_ax1=0,_ay1=0,_m2=1,_x2=-5,_y2=0,_vx2=1,_vy2=0,_ax2=0,_ay2=0,_m3=1,_x3=0,_y3=1,_vx3=0.5,_vy3=-0.5,_ax3=0,_ay3=0,_t=0,_dt=0.01,_k=1,_n=1000):self.x1=[]self.x1.append(_x1)self.vx1=[]self.vx1.append(_vx1)self.ax1=[]self.ax1.append(_ax1)self.y1=[]self.y1.append(_y1)self.vy1=[]self.vy1.append(_vy1)self.ay1 = []self.ay1.append(_ay1)self.x2=[]self.x2.append(_x2)self.vx2=[]self.vx2.append(_vx2)self.ax2 = []self.ax2.append(_ax2)self.y2=[]self.y2.append(_y2)self.vy2=[]self.vy2.append(_vy2)self.ay2 = []self.ay2.append(_ay2)self.x3=[]self.x3.append(_x3)self.vx3=[]self.vx3.append(_vx3)self.ax3 = []self.ax3.append(_ax3)self.y3=[]self.y3.append(_y3)self.vy3=[]self.vy3.append(_vy3)self.ay3 = []self.ay3.append(_ay3)self.m1=_m1self.m2=_m2self.m3=_m3self.dt=_dtself.n=_nself.t=[]self.t.append(_t)self.k=_kreturndef calculate_speedverlet(self):while self.t[-1]<self.n:self.x1.append(self.x1[-1] + self.vx1[-1] * self.dt + self.ax1[-1] * 1 / 2 * self.dt ** 2)self.x2.append(self.x2[-1] + self.vx2[-1] * self.dt + self.ax2[-1] * 1 / 2 * self.dt ** 2)self.x3.append(self.x3[-1] + self.vx3[-1] * self.dt + self.ax3[-1] * 1 / 2 * self.dt ** 2)self.y1.append(self.y1[-1] + self.vy1[-1] * self.dt + self.ay1[-1] * 1 / 2 * self.dt ** 2)self.y2.append(self.y2[-1] + self.vy2[-1] * self.dt + self.ay2[-1] * 1 / 2 * self.dt ** 2)self.y3.append(self.y3[-1] + self.vy3[-1] * self.dt + self.ay3[-1] * 1 / 2 * self.dt ** 2)self.ax1.append(-self.k / self.m1 * (self.x1[-1] - self.x2[-1] - (self.x1[0] - self.x2[0]))-self.k / self.m1 * (self.x1[-1] - self.x3[-1] - (self.x1[0] - self.x3[0])))self.ax2.append(-self.k / self.m2 * (self.x2[-1] - self.x1[-1] - (self.x2[0] - self.x1[0]))-self.k / self.m2 * (self.x2[-1] - self.x3[-1] - (self.x2[0] - self.x3[0])))self.ax3.append(-self.k / self.m3 * (self.x3[-1] - self.x1[-1] - (self.x3[0] - self.x1[0]))-self.k / self.m3 * (self.x3[-1] - self.x2[-1] - (self.x3[0] - self.x2[0])))self.ay1.append(-self.k / self.m1 * (self.y1[-1] - self.y2[-1] - (self.y1[0] - self.y2[0]))- self.k / self.m1 * (self.y1[-1] - self.y3[-1] - (self.y1[0] - self.y3[0])))self.ay2.append(-self.k / self.m2 * (self.y2[-1] - self.y1[-1] - (self.y2[0] - self.y1[0]))- self.k / self.m2 * (self.y2[-1] - self.y3[-1] - (self.y2[0] - self.y3[0])))self.ay3.append(-self.k / self.m3 * (self.y3[-1] - self.y1[-1] - (self.y3[0] - self.y1[0]))- self.k / self.m3 * (self.y3[-1] - self.y2[-1] - (self.y3[0] - self.y2[0])))self.vx1.append(self.vx1[-1] + 1 / 2 * self.dt * (self.ax1[-1] + self.ax1[-2]))self.vx2.append(self.vx2[-1] + 1 / 2 * self.dt * (self.ax2[-1] + self.ax2[-2]))self.vx3.append(self.vx3[-1] + 1 / 2 * self.dt * (self.ax3[-1] + self.ax3[-2]))self.vy1.append(self.vy1[-1] + 1 / 2 * self.dt * (self.ay1[-1] + self.ay1[-2]))self.vy2.append(self.vy2[-1] + 1 / 2 * self.dt * (self.ay2[-1] + self.ay2[-2]))self.vy3.append(self.vy3[-1] + 1 / 2 * self.dt * (self.ay3[-1] + self.ay3[-2]))self.t.append(self.t[-1] + self.dt)return self.t , self.x1,self.vx1def plot_tx(self):plot(self.t,self.x1,'r',lw=1,label='Speed Verlet Algorithm')plot(self.t,self.x2,'yellowgreen',lw=1.5,label='Speed Verlet Algorithm')plot(self.t,self.x3,'orange',lw=1,label='Speed Verlet Algorithm')xlabel('t')ylabel('x')grid(True)legend(loc='right')returndef plot_ty(self):plot(self.t,self.y1,'r',lw=1,label='Speed Verlet Algorithm')plot(self.t,self.y2,'yellowgreen',lw=1.5,label='Speed Verlet Algorithm')plot(self.t,self.y3,'orange',lw=1,label='Speed Verlet Algorithm')xlabel('t')ylabel('y')grid(True)legend(loc='right')returndef plot_2d(self):plot(self.x1,self.y1,'r',lw=1,label='Speed Verlet Algorithm')plot(self.x2,self.y2,'yellowgreen',lw=1.5,label='Speed Verlet Algorithm')plot(self.x3,self.y3,'orange',lw=1,label='Speed Verlet Algorithm')xlabel('x')ylabel('y')grid(True)legend(loc='right')returndef plot_xp(self):plot(self.x1, self.m1*self.vx1, 'steelblue', label='speed verlet algorithm',lw="2")grid(True)class Harmonics_2(Harmonics):def calculate_verlet(self):self.t.append(self.dt)self.x1.append(5 + self.vx1[0]*self.dt+ self.k/self.m1*self.ax1[0]*1/2*self.dt**2)self.y1.append(0.5 + self.vy1[0] * self.dt + self.k / self.m1 * self.ay1[0] * 1 / 2 * self.dt ** 2)self.x2.append(-5 + self.vx2[0] * self.dt + self.k / self.m2 * self.ax2[0] * 1 / 2 * self.dt ** 2)self.y2.append(0 + self.vy2[0] * self.dt + self.k / self.m2 * self.ay2[0] * 1 / 2 * self.dt ** 2)self.x3.append(0 + self.vx3[0] * self.dt + self.k / self.m3 * self.ax3[0] * 1 / 2 * self.dt ** 2)self.y3.append(1 + self.vy3[0] * self.dt + self.k / self.m3 * self.ay3[0] * 1 / 2 * self.dt ** 2)while self.t[-1]<self.n:self.ax1.append(-self.k / self.m1 * (self.x1[-1] - self.x2[-1] - (self.x1[0] - self.x2[0]))- self.k / self.m1 * (self.x1[-1] - self.x3[-1] - (self.x1[0] - self.x3[0])))self.ax2.append(-self.k / self.m2 * (self.x2[-1] - self.x1[-1] - (self.x2[0] - self.x1[0]))- self.k / self.m2 * (self.x2[-1] - self.x3[-1] - (self.x2[0] - self.x3[0])))self.ax3.append(-self.k / self.m3 * (self.x3[-1] - self.x1[-1] - (self.x3[0] - self.x1[0]))- self.k / self.m3 * (self.x3[-1] - self.x2[-1] - (self.x3[0] - self.x2[0])))self.ay1.append(-self.k / self.m1 * (self.y1[-1] - self.y2[-1] - (self.y1[0] - self.y2[0]))- self.k / self.m1 * (self.y1[-1] - self.y3[-1] - (self.y1[0] - self.y3[0])))self.ay2.append(-self.k / self.m2 * (self.y2[-1] - self.y1[-1] - (self.y2[0] - self.y1[0]))- self.k / self.m2 * (self.y2[-1] - self.y3[-1] - (self.y2[0] - self.y3[0])))self.ay3.append(-self.k / self.m3 * (self.y3[-1] - self.y1[-1] - (self.y3[0] - self.y1[0]))- self.k / self.m3 * (self.y3[-1] - self.y2[-1] - (self.y3[0] - self.y2[0])))self.x1.append(2 * self.x1[-1] - self.x1[-2] + (self.ax1[-1] * self.dt ** 2))self.x2.append(2 * self.x2[-1] - self.x2[-2] + (self.ax2[-1] * self.dt ** 2))self.x3.append(2 * self.x3[-1] - self.x3[-2] + (self.ax3[-1] * self.dt ** 2))self.y1.append(2 * self.y1[-1] - self.y1[-2] + (self.ay1[-1] * self.dt ** 2))self.y2.append(2 * self.y2[-1] - self.y2[-2] + (self.ay2[-1] * self.dt ** 2))self.y3.append(2 * self.y3[-1] - self.y3[-2] + (self.ay3[-1] * self.dt ** 2))self.vx1.append((self.x1[-1] - self.x1[-3]) / (2 * self.dt))self.vx2.append((self.x2[-1] - self.x2[-3]) / (2 * self.dt))self.vx3.append((self.x3[-1] - self.x3[-3]) / (2 * self.dt))self.vy1.append((self.y1[-1] - self.y1[-3]) / (2 * self.dt))self.vy2.append((self.y2[-1] - self.y2[-3]) / (2 * self.dt))self.vy3.append((self.y3[-1] - self.y3[-3]) / (2 * self.dt))self.t.append(self.t[-1] + self.dt)return self.t ,self.x1,self.vx1def plot_tx(self):plot(self.t,self.x1,'r',label='Verlet Algorithm',linestyle=':',lw=2.5)plot(self.t,self.x2,'yellowgreen',label='Verlet Algorithm',linestyle=':',lw=2.5)plot(self.t,self.x3,'orange',label='Verlet Algorithm',linestyle=':',lw=2.5)xlabel('t')ylabel('x')grid(True)legend(loc='right')returndef plot_ty(self):plot(self.t,self.y1,'r',label='Verlet Algorithm',ls=':',lw=2.5)plot(self.t,self.y2,'yellowgreen',label='Verlet Algorithm',linestyle=':',lw=2.5)plot(self.t,self.y3,'orange',label='Verlet Algorithm',linestyle=':',lw=2.5)xlabel('t')ylabel('y')grid(True)legend(loc='right')returndef plot_2d(self):plot(self.x1,self.y1,'r',label='Verlet Algorithm',ls=':',lw=2.5)plot(self.x2,self.y2,'yellowgreen',label='Verlet Algorithm',ls=':',lw=2.5)plot(self.x3,self.y3,'orange',label='Verlet Algorithm',ls=':',lw=2.5)xlabel('x')ylabel('y')grid(True)legend(loc='right')returndef plot_xp(self):x1_cut = self.x1[:-1]plot(x1_cut, self.m1*self.vx1, 'r', label='verlet algorithm',ls=':',lw=2.5)grid(True)a=Harmonics()a.calculate_speedverlet()fig = figure()def call_back(event):axtemp = event.inaxesx_min, x_max = axtemp.get_xlim()increment = (x_max - x_min) / 10if event.button == 'up':axtemp.set(xlim=(x_min + increment, x_max - increment))print('up')elif event.button == 'down':axtemp.set(xlim=(x_min - increment, x_max + increment))print('down')fig.canvas.draw_idle() # 绘图动作实时反映在图像上fig.canvas.mpl_connect('scroll_event', call_back)fig.canvas.mpl_connect('button_press_event', call_back)subplot(3,1,1)a.plot_tx()title("Three Harmonic Oscillators")#title("speed verlet algorithm")subplot(3,1,2)a.plot_ty()subplot(3,1,3)a.plot_2d()##############################b=Harmonics_2()b.calculate_verlet()subplot(3,1,1)b.plot_tx()#title("verlet algorithm")subplot(3,1,2)b.plot_ty()subplot(3,1,3)b.plot_2d()show()###############################fig = figure()def call_back(event):axtemp = event.inaxesx_min, x_max = axtemp.get_xlim()increment = (x_max - x_min) / 10if event.button == 'up':axtemp.set(xlim=(x_min + increment, x_max - increment))print('up')elif event.button == 'down':axtemp.set(xlim=(x_min - increment, x_max + increment))print('down')fig.canvas.draw_idle() # 绘图动作实时反映在图像上fig.canvas.mpl_connect('scroll_event', call_back)fig.canvas.mpl_connect('button_press_event', call_back)a.plot_xp()b.plot_xp()xlabel('x',fontsize="15")ylabel('p',fontsize="15")title("Three Harmonic Oscillators")legend()show()#输出坐标文件t=a.calculate_speedverlet()[0]x1=a.calculate_speedverlet()[1]vx1=a.calculate_speedverlet()[2]coordinate = open("speed_verlet_three", 'w')for i in range(len(t)):coordinate.write(str(t[i]) + "\t" + str(x1[i]) + "\t" + str(vx1[i]) + "\n")coordinate.close()#输出坐标文件t_v=b.calculate_verlet()[0]x1_v=b.calculate_verlet()[1]vx1_v=b.calculate_verlet()[2]coordinate = open("verlet_three", 'w')for i in range(len(t)):coordinate.write(str(t_v[i]) + "\t" + str(x1_v[i]) + "\n")coordinate.close()coordinate = open("verlet_three_v", 'w')for i in range(len(vx1)):coordinate.write( str(vx1_v[i])+ "\n")coordinate.close()
如果觉得《比较分子动力学模拟中Verlet算法 Speed Verlet算法及辛算法的精度》对你有帮助,请点赞、收藏,并留下你的观点哦!