速度verlet算法简介及其程序实现

前面介绍的verlet算法存在的一个问题是不能同时获得当前时刻的位置和速度，那么也就不能获得当前状态下动能，温度等等。这在实际的积分过程中是有点纠结的。实际应用中，速度verlet（velocity-verlet）是更为实用的方法，大部分的分子动力学模拟软件实际采用的就是这种方法（lammps就是）。这种方法可以同时给出位置、速度和加速度，并且不牺牲精度。

速度verlet算法的介绍可以参考维基百科：velocity-verlet 。

速度verlet算法简要介绍

位置表达式



可以看出， 下一时刻的位置只依赖于当前时刻的位置、速度和加速度。

速度表达式



可以看出，下一时刻的速度依赖于当前时刻的速度、加速度和下一时刻的加速度。

速度verlet算法积分过程

启动速度verlet算法需要知道最初t 时刻的位置、速度和加速度。具体积分过程如下：

1. 根据 t 时刻的位置，速度和加速度获得 t+Δt 时刻的位置。

2. 根据 t+Δt 时刻的位置，获得 t+Δt 时刻的加速度。

3. 根据 t 时刻的速度和加速度，t+Δt 时刻的加速度，获得 t+Δt 时刻的速度。

这样，就获得了t+Δt 时刻的位置、速度和加速度（其中，加速度先获得）。这样，通过一个积分步，就将状态向前推进了一个时间步（可以对比红色部分）。

下图反映了这样一个过程。

速度verlet算法程序实现

仍然使用 verlet积分算法简介及程序实现 中的例子进行介绍。

重复部分，不再赘述，直接给出MATLAB代码：velocity_verlet.m

function [r,v] = velocity_verlet(dt,totstep)
% This script calculate the displacement, velocity, force
% of harmonic ossication using velocity-verlet methods.
% Input:
%	dt: timestep  (recommend: 0.01 ~ 2.5)
% 	totstep: total steps.
% Powered By Xianbao Duan
% Email: xianbao.d@gmail.com
% Website: www.52souji.net

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% log:
% 2012-06-05: Complete
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% basic parameters
k = 1;
m = 1;
r0 = 2;
v0 = 2*sqrt(3);

% velocity verlet algorithm
r(1) = r0;	v(1) = v0;	F(1) = -k*r(1);

for i = 1:(totstep-1)
	r(i+1) = r(i)+dt*v(i)+0.5*F(i)*dt*dt/m;
	F(i+1) = -k*r(i+1);
	v(i+1) = v(i)+0.5*(F(i+1)+F(i))*dt/m;
end

可以使用以下命令调用以上函数:

>> [r,v]=velocity_verlet(0.01,3000);
>> plot(1:3000,r)
>> grid on
>> xlabel('steps')
>> ylabel('displacement')
>> title('velocity-verlet algorithm, dt=0.01, totstep=3000')
>> figure
>> plot(r,v)
>> grid on
>> xlabel('displacement');
>> ylabel('velocity')
>> title('velocity-verlet algorithm, dt=0.01, totstep=3000')

下图是绘制得到的位置随着时间步的变化



下图是绘制得到的速度随着时间的变化



对比 verlet积分算法简介及程序实现 中给出的图片，两者几乎没有差别。

实际上，两个程序在给出结果上是完全等价的，这个是可以被严格证明的。