Chemistry Reference
In-Depth Information
n
¼
length(x0);
t
¼
0:h:h*nstep;
x
¼
zeros(n,nstep
þ
1);
v
¼
zeros(n,nstep
þ
1);
eout
¼
zeros(1,nstep
þ
1);
x(:,1)
¼
x0;
v(:,1)
¼
v0;
[e,f]
¼
feval(fun,x0);
eout(1)
¼
e
þ
.5*sum(M.*v0.*v0);
for i
¼
1:nstep
rf
¼
dynlang1(n,gamma,temp,h,M);
f
¼
f-rf;
vhalf
¼
(v(:,i)
þ
(h/2)*f./M)/(1
þ
h/2*gamma);
x(:,i
þ
1)
¼
x(:,i)
þ
h*vhalf;
[e,f]
¼
feval(fun,x(:,i
þ
1));
f
¼
f-rf;
v(:,i
þ
1)
¼
vhalf*(1
h/2*gamma)
þ
(h/2)*f./M;
eout(i
þ
1)
¼
e
þ
.5*sum(M.*v(:,i
þ
1).*v(:,i
þ
1));
end
%%%%%%%%%%end of function %%%%%%%%%%%%%%%%%%%%%%
function [t,x,v,eout]
¼
extrapmts(fun,h,nstep,x0,v0,
M,tau,gamma,temp);
%[t,x,v,eout]
¼
extrapmts(fun,h,nstep,x0,v0,M,tau,
gamma,temp);
% Force Extrapolation Langevin dynamics MTS integrator
%
%INPUT ARGUMENTS
% fun string containing file name of energy and force
routine
% h timestep
% nstep number of steps
% x0 vector of initial atom positions
% v0 vector of initial atom velocities
% M vector of atomic masses
% tau number of steps between slow force updates
% gamma Langevin friction parameter
% temp target temperature for Langevin dynamics
%OUTPUT ARGUMENTS
% t vector of time values
% x trajectory of positions
