function xyz = LJ7findmin()
global rstar
rstar = 2^(1/6); % argument of the minimum of the Lennard-Jones pair potential V(r) = 6^(-12) - r^(-6)

% Four lical minima of LJ7:
% f1 = -16.50538417 Pentagonal bipyramid 
% f2 = -15.93504306 Capped octahedron 
% f3 = -15.59321094 Tricapped tetrahedron 
% f4 = -15.53306005 Bicapped trigonal bipyramid
Na = 7; % the number of atoms

for model = 1 : 4
    xyz = zeros(3,Na);

    %% Set the initial guess (start with polynedrons with edges equal to 1)
    if model == 1 % Pentagonal bipyramid
        p5 = 0.4*pi;
        pp = p5/2;
        p0 = pi*1.5 - pp;
        % set coordinates of the center so that position satisfy (*)
        x0 = sin(pp);
        y0 = cos(pp);
        z0 = 0;
        he = sqrt(1 - 0.25/(sin(pp))^2);
        for k = 1 : 5
            xyz(1,k) = x0 + cos(p0 + (k-1)*p5); 
            xyz(2,k) = y0 + sin(p0 + (k-1)*p5); 
            xyz(3,k) = z0;  
        end
        xyz(:,6) = [x0, y0, he]';
        xyz(:,7) = [x0, y0, -he]';
    end

    if model == 2 % Capped octahedron
        r = 1/sqrt(2);
        p4 = 0.5*pi;
        pp = p4/2;
        p0 = pi*1.5 - pp;
        % set coordinates of the center so that position satisfy (*)
        x0 = sin(pp);
        y0 = cos(pp);
        z0 = 0;
        for k = 1 : 4
            xyz(1,k) = x0 + r*cos(p0 + (k-1)*p4);
            xyz(2,k) = y0 + r*sin(p0 + (k-1)*p4);
            xyz(3,k) = z0;
        end
        xyz(:,5) = [x0, y0, z0 + r]';
        xyz(:,6) = [x0, y0, z0 - r]';
        xyz(:,7) = [3*x0, y0, z0 + r]';
    end

    if model == 3  % Tricapped tetrahedron
        p3 = 2*pi/3;
        pp = p3/2;
        r = 1/sqrt(3);
        beta = 0.5*pi -asin(1/3) - acos(1/sqrt(3));
        r1 = cos(beta);
        p0 = 1.5*pi - pp;
        % set coordinates of the center so that position satisfy (*)
        x0 = sin(pp);
        y0 = cos(pp);
        z0 = 0;
        for k = 1 : 3
            xyz(1,k) = x0 + r*cos(p0 + (k-1)*p3);
            xyz(2,k) = y0 + r*sin(p0 + (k-1)*p3);
            xyz(3,k) = z0;
            xyz(1,k + 3) = x0 + r1*cos(p0 + pp + (k-1)*p3);
            xyz(2,k + 3) = y0 + r1*sin(p0 + pp + (k-1)*p3);
            xyz(3,k + 3) = z0 + sqrt(2/3) - sin(beta);
        end
        xyz(:,7) = [x0, y0, z0 + sqrt(2/3)]';
    end

    if model == 4 % Bicapped trigonal bipyramid
        p3 = 2*pi/3;
        pp = p3/2;
        r = 1/sqrt(3);
        beta = 0.5*pi -asin(1/3) - acos(1/sqrt(3));
        r1 = cos(beta);
        p0 = 1.5*pi - pp;
        % set coordinates of the center so that position satisfy (*)
        x0 = sin(pp);
        y0 = cos(pp);
        z0 = 0;
        for k = 1 : 3
            xyz(1,k) = x0 + r*cos(p0 + (k-1)*p3);
            xyz(2,k) = y0 + r*sin(p0 + (k-1)*p3);
            xyz(3,k) = z0;
        end
        xyz(:,4) = [x0 + r1*cos(p0 + pp), y0 + r1*sin(p0 + pp), z0 + sqrt(2/3) - sin(beta)]';
        xyz(:,5) = [x0 + r1*cos(p0 + pp + p3), y0 + r1*sin(p0 + pp+p3), z0 - sqrt(2/3) + sin(beta)]';
        xyz(:,6) = [x0, y0, z0 - sqrt(2/3)]';
        xyz(:,7) = [x0, y0, z0 + sqrt(2/3)]';
    end
    xyz = xyz*rstar;

    %% Set the initial vector 

    x_aux = xyz(:); % make xyz into column vector
    % x = [x2 x3 y3 x4 y4 z4 ... xNa yNa ZNa]';
    x = x_aux([4,7,8,10 : Na*3]);
    N = length(x);
    fprintf('Ndim = %d\n',N);

    tic % start cpu timer
    
    f = LJpot(x);
    df = LJgrad(x);
    nor = norm(df);
    fprintf('Initially, V = %d\t, |grad V| = %d\n',f,nor); 
    
    %% Newton's method
    tol = 1e-10;
    c1 = 0.1; % parameter in Wolfe's condition for sufficient decrease
    fac = 0.9;
    iter = 0;
    feval = 0;
    flag = 0;
    pmax = 1e-1;
    while nor > tol
        % find Newton's direction
        H = LJhess(x);
        p_newton = -H\df; 
        ang = -p_newton'*df/(norm(df)*norm(p_newton));
%        fprintf('cos(angle(p_newton,df)) = %d\n',ang);
        % check if it is a descend direction
        if ang > 0 % descend direction
            p = p_newton;
%            alp = 1;
            fprintf('Direction: Newton; ');
        else
            p = -df; % steepest descend;
            fprintf('Direction: Steepest descend; ');
%            alp = 0.1;
        end
        nor = norm(p);
        if nor > pmax
            p = pmax*p/nor;
        end
        % line search 
        alp = 1;
        flag = 0;
        xnew = x + alp*p;
        fnew = LJpot(xnew);
        feval = feval + 1;
        % backtracking
        while fnew > f + c1*alp*df'*p 
            alp = alp*fac;
            xnew = x + alp*p;
            fnew = LJpot(xnew);
%            fprintf('fnew = %d\n',fnew);
            feval = feval + 1;
            if alp < 1e-14
                flag = 1;
                break;
            end
        end
        
        if flag == 1
            fprintf('\nLine search failed\n');
            break;
        else
            df = LJgrad(xnew);
            nor = norm(df);
            x = xnew;
            f = fnew;
            iter = iter + 1;
            fprintf('alpha = %d;',alp);
            fprintf('iter = %d; f = %d; |df| = %d\n',iter,f,nor);
        end
    end
    
    toc % stop cpu timer
    
    if flag == 0
        fprintf('Result: f = %.8f\t |df| = %d\t niter = %d\t feval = %d\n',f,nor,iter,feval);
        x_aux = [0;0;0;x(1);0;0;x(2:3);0;x(4:N)];
        xyz = reshape(x_aux,3,Na);
        draw_polyhedron(xyz);
    end
    
    H = LJhess(x);
    evals = sort(eig(H),'ascend');
    for j = 1 : N
        fprintf('%d\n',evals(j));
    end
    
end
end
  %% Lennard-Jones potential        

 function v = LJpot(xyz) % 
 [m n] = size(xyz);
 na = (m + 6)/3 ;
 x_aux = [0;0;0;xyz(1);0;0;xyz(2:3);0;xyz(4:m)];
 m = length(x_aux);
 % restore atomic coordinates
 x = x_aux(1 : 3 : m - 2); 
 y = x_aux(2 : 3 : m - 1); 
 z = x_aux(3 : 3 : m);
 r = zeros(na);
 v = 0;
 for k = 2 : na
     for j = 1 : k - 1
         r(j,k) = sqrt((x(k)-x(j)).^2+(y(k)-y(j)).^2+(z(k)-z(j)).^2);
         r(k,j) = r(j,k);
         v = v + r(j,k).^(-12) - r(j,k).^(-6);
     end
 end
 v = v*4;
 end
 %%
function dv = LJgrad(xyz) % 
[m n] = size(xyz);
na = (m + 6)/3 ;
x_aux = [0;0;0;xyz(1);0;0;xyz(2:3);0;xyz(4:m)];
m = length(x_aux);
% restore atomic coordinates
x = x_aux(1 : 3 : m - 2); 
y = x_aux(2 : 3 : m - 1); 
z = x_aux(3 : 3 : m);
r = zeros(na);
for k = 2 : na
 for j = 1 : k - 1
     r(j,k) = sqrt((x(k)-x(j)).^2+(y(k)-y(j)).^2+(z(k)-z(j)).^2);
     r(k,j) = r(j,k);
 end
end

dVx = zeros(1,na);
dVy = zeros(1,na);
dVz = zeros(1,na);
for k = 1 : na
    for j = 1 : na
        if j ~= k
            raux = r(k,j);
            lj = LJder(raux)./raux;
            dVx(k) = dVx(k) + lj.*(x(k) - x(j));
            dVy(k) = dVy(k) + lj.*(y(k) - y(j));
            dVz(k) = dVz(k) + lj.*(z(k) - z(j));
        end
    end
end
aux = 4*[dVx; dVy; dVz];
aux = aux(:);
dv = aux([4,7,8,10 : na*3]);
end

%%
 function dv = LJder(r)
 dv = -12*r.^(-13) + 6*r.^(-7);
 end
%%
function H = LJhess(x) % find the Hessian using finite differences
h = 1e-3;
n = length(x);
H = zeros(n);
for i = 1 : n
    di = zeros(size(x));
    di(i) = h;
    for j = 1 : i
        dj = zeros(size(x));
        dj(j) = h;
        if i == j
            f1 = LJpot(x + di);
            f2 = LJpot(x - di);
            f0 = LJpot(x);
            H(i,i) = (f1 - 2*f0 + f2)/h^2;
        else
            f1 = LJpot(x + di + dj);
            f2 = LJpot(x - di + dj);
            f3 = LJpot(x + di - dj);
            f4 = LJpot(x - di - dj);
            H(i,j) = 0.25*(f1 - f2 - f3 + f4)/h^2;
            H(j,i) = H(i,j);
        end
    end
end
end
%%


