Extending the Time Scale in Atomically Detailed Simulations - Reviews in Computational Chemistry - page 402

M1 ¼ [15.994 1.008 1.008]';
M2 ¼ kron(M1,ones(2,1));
M ¼ [M2;M2];
gamma ¼ .1;
temp ¼ 300;
rand('state',0)
if method(1) ¼¼ 'v', %Velocity Verlet method
[t,x,v,eout] ¼ vv('twowaters',h,nstep,x0,v0,M);
elseif method(1) ¼¼ 'i', %Impulse MTS method
[t,x,v,eout] ¼ impulsemts('twowaters',h,tau,
nstep,x0,v0,M);
elseif method(1) ¼¼ 'b' %Langevin Dynamics
[t,x,v,eout] ¼ bbk1('twowaters',h,nstep,x0,v0,M,
gamma,temp);
elseif method(1) ¼¼ 'e' %Force Extrapolation MTS with
Langevin Dynamics
[t,x,v,eout] ¼ extrapmts('twowaters',h,nstep,x0,
v0,M,tau,gamma,temp);
else
disp('unrecognized method specification')
end
%%%%%%%%%%end of function %%%%%%%%%%%%%%%%%%%%%%
function [t,x,v,eout] ¼ vv(fun,h,nstep,x0,v0,M);
%[t,x,v,eout] ¼ vv(fun,h,nstep,x0,v0,M);
% velocity Verlet 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
%OUTPUT ARGUMENTS
% t vector of time values
% x trajectory of positions
% v trajectory of velocities
% eout energy along computed trajectory
%
%Eric Barth
%Kalamazoo College, 2007