Chemistry Reference
In-Depth Information
% v trajectory of velocities
% eout energy along computed trajectory
%
%Eric Barth
%Kalamazoo College, 2007
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,ffast,fslow]
¼
feval(fun,x0,0);
%eout(1)
¼
e
þ
.5*sum(M.*v0.*v0);
eout(1)
¼
e;
for i
¼
1:nstep
xx
¼
x(:,i);
vv
¼
v(:,i);
for j
¼
1:tau
rf
¼
dynlang1(n,gamma,temp,h,M);
fext
¼
ffast
þ
fslow-rf;
vv
¼
(vv
þ
(h/2)*fext./M)/(1
þ
h/2*gamma);
xx
¼
xx
þ
h*vv;
if (j
¼¼
ceil(tau/2)),
[e,f,ffast,fslow]
¼
feval(fun,xx,0);
else
[e,f,ffast,fdum]
¼
feval(fun,xx,1);
end
fext
¼
ffast
þ
fslow-rf;
vv
¼
vv*(1
(h/2)*gamma)
þ
(h/2)*fext./M;
end
x(:,i
þ
1)
¼
xx;
v(:,i
þ
1)
¼
vv;
eout(i
þ
1)
¼
e; %report only the bond energy
%[e,f,ffast,fslow]
¼
feval(fun,xx,0);
end
%%%%%%%%%%end of function %%%%%%%%%%%%%%%%%%%%%%
function [e,f,ffast,fslow]
¼
twowaters(x,splitflag);
%[e,f,ffast,fslow]
¼
twowaters(x,splitflag);
%Energy and Force routine used by integration methods
