%% Finite-temperature string method (FTS) code.
% applied to the rugged Mueller potential.

% Initialization
clear all

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

% temperature
mu=9;

% number of images
n1=40;


% frequency of plotting
nstepplot = 1e2;
% plot string every nstepplot if flag1 = 1 
flag1 = 1;


% number of images during prerelaxation
n2 = 1e1;
% number of prerelaxation steps
nstep1 = 2e3;
% weight for running averages
nstep2 = 1e2;
% total number of iterations
nstep3 = 2e3;
% number of iterations stored
nstepav = 5e2;
% time-step
h = 1e-4;

% string updating parameters:
% time-step
ha = 0.1;
% smoothin coefficient
kap = 0.1;



% end points of the initial string
xa = -1;
ya = 0.5;

xb = 0.7;
yb = 0.6;

% initialization
g1 = linspace(0,1,n1);
g2 = linspace(0,1,n2);
x = (xb-xa)*g2+xa;
y = (x-xa)*(yb-ya)/(xb-xa)-0.1*(x-xa).*(xb-x)+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(n2);

x = interp1(lxy,x,g2);
y = interp1(lxy,y,g2);
xi=x;
yi=y;

r1 = ones(1,n1);
r1(1) = 0;
r1(n1) = 0;

r2 = ones(1,n2);
r2(1) = 0;
r2(n2) = 0;

% parameters in rugged 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];

k = 5;
amp1 = 9;

xmin=-1.5;
xmax= 1.2;
ymin=-0.2;
ymax= 2.0;


[xx,yy] = meshgrid(xmin:0.01:xmax,ymin:0.01:ymax);


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;
V1 = V1 + amp1*sin(2*pi*k*xx).*sin(2*pi*k*yy);



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


[vx vy] = voronoi(x,y);


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


%% Pre-relaxation stage
xav = x;
yav = y;

xt = x;
yt = y;


tt = delaunay(xt,yt);


for nstep=1:nstep1
    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;

    dVx = dVx +  2*pi*k*amp1*cos(2*pi*k*x).*sin(2*pi*k*y);
    dVy = dVy +  2*pi*k*amp1*sin(2*pi*k*x).*cos(2*pi*k*y);

    x0 = x;
    y0 = y;
    
    
    x = x - h*dVx + sqrt(2*h*mu)*randn(1,n2);
    y = y - h*dVy + sqrt(2*h*mu)*randn(1,n2);     
   
    
    k1 = dsearch(xt,yt,tt,x,y);
   
    test1 = k1 == linspace(1,n2,n2);
    ntest1 = not(test1);
    x = x.*test1 + x0.*ntest1;
    y = y.*test1 + y0.*ntest1;   

    xav = (xav*(nstep+nstep2-1) + x)/(nstep+nstep2);
    yav = (yav*(nstep+nstep2-1) + y)/(nstep+nstep2);
    
    dxtmp = circshift(xt,[0 -1]) + circshift(xt,[0 1]) - 2*xt;
    dytmp = circshift(yt,[0 -1]) + circshift(yt,[0 1]) - 2*yt;
         
    xt = xt - (xt-xav)*ha + kap*ha*n2*r2.*dxtmp;
    yt = yt - (yt-yav)*ha + kap*ha*n2*r2.*dytmp;

    dxt = xt-circshift(xt,[0 1]);
    dyt = yt-circshift(yt,[0 1]);
    dxt(1) = 0;
    dyt(1) = 0;
    lxyt = cumsum(sqrt(dxt.^2+dyt.^2));
    lxyt = lxyt/lxyt(n2);

    xt = interp1(lxyt,xt,g2);
    yt = interp1(lxyt,yt,g2);
    
    tt = delaunay(xt,yt);
    k1 = dsearch(xt,yt,tt,x,y);
   
    test1 = k1 == linspace(1,n2,n2);
    ntest1 = not(test1);
    x = x.*test1 + xt.*ntest1;
    y = y.*test1 + yt.*ntest1;
    

    if and(flag1 == 1,mod(nstep,nstepplot) == 0)
        figure(2);clf;
        contourf(xx,yy,min(V1,200),40);
        hold on
        plot(xt,yt,'w.')
        [vx vy] = voronoi(xt,yt);
        plot(vx,vy,'w-')
        axis([xmin xmax ymin ymax]);
        set(gca,'XTick',-1.5:.5:1,'YTick',0:.5:2);
        xlabel('x','FontAngle','italic');
        ylabel('y','FontAngle','italic');
        drawnow
    end
    
    
end

fprintf('End Prerelaxation\n')

%% initial string for second stage: n1 images

xt = interp1(lxyt,xt,g1,'cubic');
yt = interp1(lxyt,yt,g1,'cubic');

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


% Main Loop

x = xt;
y = yt;

xav = xt;
yav = yt;





tt = delaunay(xt,yt);

cx = zeros(nstepav,n1);
cy = zeros(nstepav,n1);
flag1 = 1;

for nstep=1:nstep3
    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;

    dVx = dVx +  2*pi*k*amp1*cos(2*pi*k*x).*sin(2*pi*k*y);
    dVy = dVy +  2*pi*k*amp1*sin(2*pi*k*x).*cos(2*pi*k*y);

    x0 = x;
    y0 = y;
    
    x = x - h*dVx + sqrt(2*h*mu)*randn(1,n1);
    y = y - h*dVy + sqrt(2*h*mu)*randn(1,n1);     
    
    k1 = dsearch(xt,yt,tt,x,y);
   
    test1 = k1 == linspace(1,n1,n1);
    ntest1 = not(test1);
    x = x.*test1 + x0.*ntest1;
    y = y.*test1 + y0.*ntest1;   

    xav = (xav*(nstep+nstep2-1) + x)/(nstep+nstep2);
    yav = (yav*(nstep+nstep2-1) + y)/(nstep+nstep2);
    
    dxtmp = circshift(xt,[0 -1]) + circshift(xt,[0 1]) - 2*xt;
    dytmp = circshift(yt,[0 -1]) + circshift(yt,[0 1]) - 2*yt;
         
    xt = xt - (xt-xav)*ha + kap*ha*n1*r1.*dxtmp;
    yt = yt - (yt-yav)*ha + kap*ha*n1*r1.*dytmp;

    dxt = xt-circshift(xt,[0 1]);
    dyt = yt-circshift(yt,[0 1]);
    dxt(1) = 0;
    dyt(1) = 0;
    lxyt = cumsum(sqrt(dxt.^2+dyt.^2));
    lxyt = lxyt/lxyt(n1);

    xt = interp1(lxyt,xt,g1);
    yt = interp1(lxyt,yt,g1);
    
    tt = delaunay(xt,yt);
    k1 = dsearch(xt,yt,tt,x,y);
   
    test1 = k1 == linspace(1,n1,n1);
    ntest1 = not(test1);
    x = x.*test1 + xt.*ntest1;
    y = y.*test1 + yt.*ntest1;
    
    if nstep>nstep3-nstepav
        cx(nstep-nstep3+nstepav,:) = x;
        cy(nstep-nstep3+nstepav,:) = y;
    end;

    
    if and(flag1 == 1,mod(nstep,nstepplot) == 0)
        figure(1);clf;
        contourf(xx,yy,min(V1,200),40);
        hold on
        plot(xt,yt,'w.')
        [vx vy] = voronoi(xt,yt);
        plot(vx,vy,'w-')
        axis([xmin xmax ymin ymax]);
        set(gca,'XTick',-1.5:.5:1,'YTick',0:.5:2);
        drawnow
    end
    
end


%% Output
fprintf('\n')
fprintf('\n')
fprintf('FTS calculation with %d images\n',n1)


figure(1);clf;
contourf(xx,yy,min(V1,200),40);
hold on
plot(xt,yt,'wo','MarkerFaceColor','w','MarkerSize',4)
plot(xt,yt,'-w')
[vx vy] = voronoi(xt,yt);
plot(vx,vy,'w-','LineWidth',2)
axis([xmin xmax ymin ymax]);
set(gca,'XTick',-1.5:.5:1,'YTick',0:.5:2);
xlabel('x','FontAngle','italic');
ylabel('y','FontAngle','italic');
title('Final string and associated isocommittor surfaces');
drawnow

figure(2);clf;
contourf(xx,yy,min(V1,200),40);
hold on
plot(cx(1:nstepav,1:2:n1),cy(1:nstepav,1:2:n1),'w.')
plot(cx(1:nstepav,2:2:n1),cy(1:nstepav,2:2:n1),'.','Color',[0.75 0.75 0.75])
[vx vy] = voronoi(xt,yt);
plot(vx,vy,'-k')
contour(xx,yy,min(V1,200),40,'k');
plot(xt,yt,'ko','MarkerFaceColor','k','MarkerSize',4)
axis([xmin xmax ymin ymax]);
xlabel('x','FontAngle','italic');
ylabel('y','FontAngle','italic');

set(gca,'XTick',-1.5:.5:1,'YTick',0:.5:2);
title('Transition tube');
drawnow


