function LJ6in2Dsetup()
% Finds a Minimum Energy Path for the move of the central atom to the
% periphery using the string method. Then uses the dimer method to nail down the
% saddles along the path.
%% SET UP
global rstar N
rstar = 2^(1/6); % the minimizer of LJpair(r) = r^{-12} - r^{-6}
N = 6; % the number of atoms

%% Define the initial configuration
% the initial guess
y = zeros(2,N);
y(1,6) = 0; y(2,6) = 0; 
y(1,4) = 1; y(2,4) = 0; 
y(1,3) = 2; y(2,3) = 0; 
y(1,5) = 0.5; y(2,5) = sqrt(3)*0.5; 
y(1,2) = 1.5; y(2,2) = sqrt(3)*0.5; 
y(1,1) = 1; y(2,1) = sqrt(3); 
y = rstar*y;
x = remove_rotations_translations(y); % x is 1-by-(2N - 3)
% Find the potential minimum
[xi Vi] = fminsearch(@LJ,x);  % xi is 1-by-(2N - 3), the found minimizer, Vi is the found minimum
showpath(xi,1); % the first argument is the configuration or an array of configurations,
% the second argument is the figure number
fprintf('Initial Minimum: V = %d\n',Vi);
%% Define the final confiruration
y(1,6) = 0; y(2,6) = 0; 
y(1,4) = 1; y(2,4) = 0; 
y(1,2) = 2; y(2,2) = 0; 
y(1,5) = 0.5; y(2,5) = sqrt(3)*0.5; 
y(1,1) = 1.5; y(2,1) = sqrt(3)*0.5; 
y(1,3) = 1.5; y(2,3) = -sqrt(3)*0.5; 
y = rstar*y;
x = remove_rotations_translations(y);
% Find the potential minimum
[xf Vf] = fminsearch(@LJ,x);
showpath(xf,2);
fprintf('Final Minimum: V = %d\n',Vf);

%% Connect the initial and the final configurations 
% by a segment of straingt line: this segment will be the initial path
Npath = 51; % Set Npath images along the path
t = linspace(0,1,Npath)';  % the parameter along the path
path = interp1([0;1],[xi;xf],t);
% path: Npath-by-(2N-3) matrix
showpath(path,3);
end

%% AUXILIARY FUNCTIONS
function V = LJ(x1) 
[m n] = size(x1);
Na = (n + 3)/2;
V = zeros(m,1);
% x is a matrix Npath-by-(2N-3) 
% Computes the Lennard-Jones potential 
% V = 4*sum_(i < j) [rij^{-12} - r^{-6}]
% at image of the path of m images
V = zeros(m,1); 
for j = 1 : m
    x = restore(x1(j,:));
    g = zeros(size(x));
    r = zeros(Na);
    for k = 1 : Na
        r(k,:) = sqrt(sum((x - x(:,k)*ones(1,Na)).^2,1));
    end
    r = r + diag(ones(Na,1));
    aux = 1./r.^6;
    L = (aux - 1).*aux;
    L = L - diag(diag(L)); 
    V(j) = 2*sum(sum(L)); 
end
end
%%
function dV = LJgrad(x1)
[m n] = size(x1);
Na = (n + 3)/2;
% x is a matrix Npath-by-(2N-3) 
% Computes the gradient of the Lennard-Jones potential 
% dV/dx_i = 4*sum_{i\neq j}(-12r_{ij}^{-13} + 6r_{ij}^{-7})*(x_i/r_{ij})
dV = zeros(m,n); 
for j = 1 : m
    x = restore(x1(j,:));
    g = zeros(size(x));
    r = zeros(Na);
    for k = 1 : Na
        r(k,:) = sqrt(sum((x - x(:,k)*ones(1,Na)).^2,1));
    end
    r = r + diag(ones(Na,1));
    r = max(r,0.95);
    L = -6*(2./r.^6 - 1)./r.^8;
    for k = 1 : Na
        Lk = L(:,k);
        g(:,k) = (x(:,k)*ones(1,Na) - x)*Lk;
    end
    g = g(:)';
    g([1 2 4]) = [];
    dV(j,:) = 4*g;
end
end

%%
function showpath(path1,fig)
global rstar N
rst = rstar/2;
figure(fig);
clf;
hold on;
axis equal;
axis off;
[m n] =size(path1);
path = [zeros(m,2), path1(:,1), zeros(m,1), path1(:,2 : end)];
[m n] = size(path);
t = linspace(0,2*pi,100);
% Draw the initial configuration
x = path(1,:);
col = parula(N+3);
for i = 1 : N
     conf(i) = plot(x(2*i - 1)+rst*cos(t),x(2*i)+rst*sin(t),'Linewidth',3,'color',col(i,:));
end
drawnow
if m > 1
    for j = 2 : m
        x = path(j,:);
        for i = 1 : N
            set(conf(i),'xdata',x(2*i-1)+rst*cos(t),'ydata',x(2*i)+rst*sin(t));
        end
        pause(0.02)
        drawnow
    end
end
end
%%
function plot_potential(path,fig,col)
[m n] = size(path);
g = linspace(0,1,m);
V = LJ(path);
figure(fig);
plot(g,V,'Linewidth',2,'color',col);
drawnow;
end
%%
function x = remove_rotations_translations(xy)
% Places atom 1 at the origin, place atom 2 to the x-axis, removes
% coordinates that are always 0.
% The output is a row vector;
[m Na] = size(xy);
if m ~= 2 || Na < 2
    fprintf('Error in remove_rotations_translations: [m = %d, Na = %d]\n',m,Na);
    return
end
    % shift atom 1 to the origin;
    xy = xy - xy(:,1)*ones(1,Na);
    % Perform rotation around the origin to place atom 3 onto the xy-plane
    R = eye(2);
    a = xy(:,2);
    r = norm(a);
    if r > 1e-12
        R(1,1) = a(1)/r;
        R(2,2) = R(1,1);
        R(1,2) = a(2)/r;
        R(2,1) = -R(1,2);
        xy = R*xy;
    end
    %% prepare input vector
    x_aux = xy(:); % make xyz into column vector
    % x = [x2 x3 y3 x4 y4 ... xNa yNa]';
    x = x_aux([3,5 : Na*2])';
%     fprintf('Ndim = %d\n',N);

end

%%
function xy = restore(x)
% Restores x into 2-by-Na array
N = (length(x) + 3)/2;
xy = [0,0,x(1),0,x(2:end)]';
xy = reshape(xy,2,N);
end
%%

