function Atom()
    fplot(@(x)(1/x^12-2/x^6),[0.8 2]);
    waitforbuttonpress
    
    X=fminsearch(@F,[0; 0; 0; 1; 0; 0; 0 ; 1 ;0 ; 0 ; 0; 1]);
    X=reshape(X,3,4)
    plot3(X(1,:),X(2,:),X(3,:),'X');
    waitforbuttonpress


    X=rand(30,1);
    X=fmincon(@F,X,[],[],[],[],[],[],@mycon);
    F(X)
    mycon(X)
    X=reshape(X,3,10)
    
    t=-1:0.005:1;
    z=t;
    r=(1-t.^2).^(1/2);
    x=r.*sin(20*t); 
    y=r.*cos(20*t);
    plot3(x,y,z);
    hold on;
    plot3(X(1,:),X(2,:),X(3,:),'X');
    hold off;
    waitforbuttonpress

    options = optimset('LargeScale','off','MaxFunEvals',100000);
    X=rand(30,1);
    X=fmincon(@F,X,[],[],[],[],[],[],@mycon,options);
    F(X)
    mycon(X)
    X=reshape(X,3,10)

    t=-1:0.005:1;
    z=t;
    r=(1-t.^2).^(1/2);
    x=r.*sin(20*t); 
    y=r.*cos(20*t);
    
    plot3(X(1,:),X(2,:),X(3,:),'X');
    hold on;
    plot3(x,y,z,'g');
    hold off;
    waitforbuttonpress
    

function y=F(x)
    [sx,sy] = size(x);
    x=reshape(x,3,sx/3)';
    y=0;
    for i=1:sx/3
        for j=i+1:sx/3
            y=y+LJP(x(i,:),x(j,:));
        end
    end

function p=LJP(x,y)
    d=sum((x-y).^2);
    if d==0 
        p=inf;
    else
        p=(1/d^6-2/d^3);
    end

function [c,ceq] = mycon(x)
    [sx,sy] = size(x);
    x=reshape(x,3,sx/3)';
    y=0;
    c=zeros(sx/3,1);
    for i=1:sx/3
        c(i) = sum(x(i,:).^2)-1;  
    end

    ceq = [];