%% Zero-temperature string method (ZTS) code 
% applied to the Mueller potential.


clear all

set(0,'DefaultTextFontName','TimesRoman')
set(0,'DefaultAxesFontSize',16)


% max number of iterations
nstepmax = 2e3;
% frequency of plotting
nstepplot = 1e1;

% plot string every nstepplot if flag1 = 1 
flag1 = 1;

% parameter used as stopping criterion
tol1 = 1e-7;


% number of images along the string (try from  n1 = 3 up to n1 = 1e4)
n1 = 25;

% time-step (limited by the ODE step on line 83 & 84 but independent of n1)
h = 1e-4;

% end points of the initial string
% notice that they do NOT have to be at minima of V -- the method finds
% those automatically
xa = -1;
ya = 0.5;

xb = 0.7;
yb = 0.5;

% initialization
g1 = linspace(0,1,n1);
x = (xb-xa)*g1+xa;
y = (x-xa)*(yb-ya)/(xb-xa)+ya;

dx = x-circshift(x,[0 1]);
dy = y-circshift(y,[0 1]);
dx(1) = 0;
dy(1) = 0;
lxy = cumsum(sqrt(dx.^2+dy.^2));
lxy = lxy/lxy(n1);
x = interp1(lxy,x,g1);
y = interp1(lxy,y,g1);
xi = x;
yi = y;

% parameters in Mueller potential
aa = [-1 -1 -6.5 0.7];
bb = [0 0 11 0.6];
cc = [-10 -10 -6.5 0.7];
AA = [-200 -100 -170 15];

XX = [1 0 -0.5 -1];
YY = [0 0.5 1.5 1];

[xx,yy] = meshgrid(-1.5:0.01:1.2,-0.2:0.01:2);

V1 = AA(1)*exp(aa(1)*(xx-XX(1)).^2+bb(1)*(xx-XX(1)).*(yy-YY(1))+cc(1)*(yy-YY(1)).^2);
for j=2:4
    V1 =  V1 + AA(j)*exp(aa(j)*(xx-XX(j)).^2+bb(j)*(xx-XX(j)).*(yy-YY(j))+cc(j)*(yy-YY(j)).^2);
end;


figure(1);clf;
contourf(xx,yy,min(V1,200),40);
hold on
plot(xi,yi,'.-w','MarkerSize',14)
set(gca,'XTick',-1.5:.5:1,'YTick',0:.5:2);
xlabel('x','FontAngle','italic');
ylabel('y','FontAngle','italic');
title('Initial string');
drawnow





figure(2);clf;
contourf(xx,yy,min(V1,200),40);
whitedots = line(x,y,'linestyle','none','marker','.','color','w','MarkerSize',14);
set(gca,'XTick',-1.5:.5:1,'YTick',0:.5:2);
xlabel('x','FontAngle','italic');
ylabel('y','FontAngle','italic');
title('String evolution')
drawnow






%% Main loop

for nstep=1:nstepmax

    
% calculation of the x and y-components of the force, dVx and dVy respectively
    ee = AA(1)*exp(aa(1)*(x-XX(1)).^2+bb(1)*(x-XX(1)).*(y-YY(1))+cc(1)*(y-YY(1)).^2);
    dVx = (2*aa(1)*(x-XX(1))+bb(1)*(y-YY(1))).*ee;
    dVy = (bb(1)*(x-XX(1))+2*cc(1)*(y-YY(1))).*ee;
    for j=2:4
        ee = AA(j)*exp(aa(j)*(x-XX(j)).^2+bb(j)*(x-XX(j)).*(y-YY(j))+cc(j)*(y-YY(j)).^2);
        dVx = dVx + (2*aa(j)*(x-XX(j))+bb(j)*(y-YY(j))).*ee;
        dVy = dVy + (bb(j)*(x-XX(j))+2*cc(j)*(y-YY(j))).*ee;
    end;
   
    x0 = x;
    y0 = y;

    
% string steps:

% 1. evolve
    x = x - h*dVx;
    y = y - h*dVy;

% 2. reparametrize  
    dx = x-circshift(x,[0 1]);
    dy = y-circshift(y,[0 1]);
    dx(1) = 0;
    dy(1) = 0;
    lxy = cumsum(sqrt(dx.^2+dy.^2));
    lxy = lxy/lxy(n1);

    x = interp1(lxy,x,g1);
    y = interp1(lxy,y,g1);
 
    if and(flag1 == 1,mod(nstep,nstepplot) == 0)
        set(whitedots,'xdata',x,'ydata',y) 
        drawnow
        pause(0.025)
    end
    
    tol = (norm(x-x0)+norm(y-y0))/n1;
    
    if tol <= tol1; break; end;
    

end



%% Output
fprintf('\n')
    fprintf('\n')
    fprintf('ZTS calculation with %d images\n',n1)
if tol > tol1
    fprintf('The calculation failed to converge after %d iterations\n',nstep)
else
    fprintf('The calculation terminated after %d iterations\n',nstep)
end;



figure(2);clf;
contourf(xx,yy,min(V1,200),40);
hold on
plot(x,y,'.-w','MarkerSize',14)
set(gca,'XTick',-1.5:.5:1,'YTick',0:.5:2);
xlabel('x','FontAngle','italic');
ylabel('y','FontAngle','italic');
title('Final string');
drawnow


% Energy along MEP

tx = circshift(x,[0 -1])-circshift(x,[0 1]);
ty = circshift(y,[0 -1])-circshift(y,[0 1]);

% potential computed as integral of projection of gradV on string tangent
Vz=cumtrapz(tx.*dVx + ty.*dVy);
Vz=0.5*Vz;
Vz=Vz-min(Vz);


ntxy = sqrt(tx.*tx+ty.*ty);
tx = tx./ntxy;
ty = ty./ntxy;


% err is an estimate of the error between disrectized MEP and actual one 
% err scale as 1/n1
err = trapz(1-(tx.*dVx+ty.*dVy).^2./(dVx.*dVx+dVy.*dVy))/(n1-1);
fprintf('Estimate of difference between discretized MEP and actual one: %f\n',err)

%
% Reinterpolate string with lots of points to get accurate energy along it
g0 = linspace(0,1,1e3);
x0 = interp1(lxy,x,g0);
y0 = interp1(lxy,y,g0);
ee = AA(1)*exp(aa(1)*(x0-XX(1)).^2+bb(1)*(x0-XX(1)).*(y0-YY(1))+cc(1)*(y0-YY(1)).^2);
V0=ee;
dVx = (2*aa(1)*(x0-XX(1))+bb(1)*(y0-YY(1))).*ee;
dVy = (bb(1)*(x0-XX(1))+2*cc(1)*(y0-YY(1))).*ee;
for j=2:4
    ee = AA(j)*exp(aa(j)*(x0-XX(j)).^2+bb(j)*(x0-XX(j)).*(y0-YY(j))+cc(j)*(y0-YY(j)).^2);
    V0 = V0 + ee;
    dVx = dVx + (2*aa(j)*(x0-XX(j))+bb(j)*(y0-YY(j))).*ee;
    dVy = dVy + (bb(j)*(x0-XX(j))+2*cc(j)*(y0-YY(j))).*ee;
end;
V0=V0-min(V0);

figure(3);clf;
hold on
title('Energy along MEP');
plot(g1,Vz)
plot(g0,V0,'r')
box on
legend('Thermodynamic integration along string','Exact','Location','South')
legend('boxoff')


