function xyz = LJ7_NonlinCG_random_initial_conf()
global rstar feval

% You need to insert Polak - Ribiere and Fletcher - Reeves 
% nonlinear conjugate methods into designated places


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

Na = 7; % the number of atoms
hR = 0.01;

Nruns = 10; % the number of experiments with random initial configurations

% result: format
% run, PR: energy, PR: # iter, PR:# feval, FR: energy, FR: #iter, FR #feval
% # iter = the number of iterations
% # feval = the number of calls of LJ7 that evaluates f and grad f
result = zeros(Nruns,7);

for run = 1 : Nruns
    %% set up a random initial configuration
    xyz = zeros(3,Na);
    xyz(:,1) = [0;0;0];
    a = randn(3,Na - 1);
    rad = sqrt(a(1,:).^2 + a(2,:).^2 + a(3,:).^2);
    a = a*diag(1./rad);
    for i = 2 : Na
        clear rad
        clear x
        rad = sqrt(xyz(1,1 : i - 1).^2 + xyz(2,1 : i - 1).^2 + xyz(3,1 : i - 1).^2);
        R = max(rad) + rstar;
        xa = R*a(:,i - 1)
        x = [xyz(:,1 : i - 1), xa];
        [f, df] = LJ(x(:));
        fnew = f;
        while 1
            R = R - hR;
            xa = R*a(:,i - 1);
            x = [xyz(:,1 : i - 1), xa];
            f = fnew;
            [fnew, dfnew] = LJ(x(:));
            if fnew > f
                break;
            end
        end
        xyz(:,i) = xa;
    end
    cmass = mean(xyz,2);
    xyz = xyz - cmass*ones(1,Na);
    drawconf(xyz);

    %% Set the initial vector 

    x0 = xyz(:); % make it into a column vector
    N = length(x0);
    fprintf('Ndim = %d\n',N);

    tol = 1e-6; % tolerance for the stopping criterion of the while-cycle
    % parameter in strong Wolfe's conditions 
    c1 = 0.1;
    c2 = 0.3;
    fac = 0.9; % parameter for backtracking

    %% Insert Polak - Ribiere Nonlinear congugate gradient methods here
    tic % start cpu timer
    feval = 0;
    x = x0;
    [f, g] = LJ(x); % evaluate f and grad f
    nor = norm(g);
    fprintf('Initially, V = %d\t, |grad V| = %d\n',f,nor); 
    iter = 0;
    p = -g; % Set the initial direction to be the steepest descend
    while nor > tol

    
    % INSERT YOUR CODE HERE
    
    
    end
    
    toc % stop cpu timer
    
    fprintf('Result: f = %.8f\t |df| = %d\t niter = %d\t feval = %d\n',f,nor,iter,feval);
    xyz = reshape(x,3,Na);
    drawconf(xyz);
    
    % record the result into the table
    result(run,1 : 4) = [run, f, iter, feval];

%%   Insert Fletcher - Reeves Nonlinear congugate gradient methods here 
    tic % start cpu timer
    feval = 0;
    x = x0;
    [f, g] = LJ(x); % evaluate f and grad f
    nor = norm(g);
    fprintf('Initially, V = %d\t, |grad V| = %d\n',f,nor); 
    iter = 0;
    p = -g; % Set the initial direction to be the steepest descend
    while nor > tol

    
    % INSERT YOUR CODE HERE
    
    
    end
    
    toc % stop cpu timer
    
    fprintf('Result: f = %.8f\t |df| = %d\t niter = %d\t feval = %d\n',f,nor,iter,feval);
    xyz = reshape(x,3,Na);
    drawconf(xyz);
    
    % record the result into the table
    result(run,5 : 7) = [f, iter, feval];
    
end
for run = 1 : Nruns
    fprintf('run %d: V = %.8f\tNiter = %d\t feval = %d\n',...
        run,result(run,2),result(run,3),result(run,4),result(run,5),result(run,6),result(run,7));
end
end

 %% Evaluates LJ energy and its gradient
function [v, dv] = LJ(xyz) % 
global feval
feval = feval + 1;
m = length(xyz);
Na = m/3;
x = reshape(xyz,3,Na);
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 = 2*sum(sum(L));


g = zeros(size(xyz));
L = -6*(2./r.^6 - 1)./r.^8;
for k = 1 : Na
    Lk = L(:,k);
    g(1 + (k-1)*3 : k*3) = (x(:,k)*ones(1,Na) - x)*Lk;
end

dv = 4*g;
end

%% Draws configuration of atoms given by xyz
function drawconf(xyz)
global rstar
r = 0.5*rstar;
x = xyz(1,:); y = xyz(2,:); z = xyz(3,:);
na = length(x);
Mx=max(max(x))+1;
mx=min(min(x))-1;
My=max(max(y))+1;
my=min(min(y))-1;
Mz=max(max(z))+1;
mz=min(min(z))-1;
[x1,x2,x3]=sphere;
x1=x1*r; x2=x2*r;x3=x3*r;
figure;
hold on;
for j=1:na
    surf(x1+x(j),x2+y(j),x3+z(j)); 
end
set(gca,'DataAspectRatio',[1 1 1]);
axis([mx Mx my My mz Mz]);
drawnow
end