Lennard Jones

%% Clear variables
clear all;

%% Parameters
N = 3000;    %Steps
Dt = 10^-12;   %Delta t
movie = 0; % create video of movement

%Potentiaal
eps = 1;
sigma = 1;
Rc = 2.5*sigma;
m = 1;
%Deeltjes en volume
Num = 20; % Number of particles
rho = 0.2;
L = (Num/rho)^(1/3);

Np = round(Num^(1/3)); %Number of planes in cubic lattice

%Beginsettings en herschalingsparam
T0 = 2;
tau = 10^2; %Characteristic time until Tkin ~ T0 !!
state = 1; %Rescaling on => state != 0, off => state == 0;

%arrays
x = zeros(N,Np^3);
y = x;
z = x;

Vx = zeros(N,Np^3);
Vy = Vx;
Vz = Vx;

Ax = zeros(1,Np^3);
Ay = Ax;
Az = Ax;

K = zeros(N,1);
P = K;
E = K;

T = zeros(N,1); %Kinetic temperature
%% Initialize positions
g = 1;
for ii = 1:Np
    for jj = 1:Np
        for kk = 1:Np
            x(1,g) = (ii-1)*L/(Np);
            y(1,g) = (jj-1)*L/(Np);
            z(1,g) = (kk-1)*L/(Np);

           g = g+1;
           if g == Np^3
               break;
           end
        end
        if g == Np^3
               break;
        end
    end
    if g == Np^3
               break;
    end

end

%% initialize velocity

Vx(1,:) = sqrt(T0)*randn(Np^3,1).';
Vy(1,:) = sqrt(T0)*randn(Np^3,1).';
Vz(1,:) = sqrt(T0)*randn(Np^3,1).';

som = sum(Vx(1,:))+sum(Vy(1,:))+sum(Vz(1,:));
shift = som/(3*Np^3)*ones(1,Np^3);
Vx(1,:) = Vx(1,:)-shift;
Vy(1,:) = Vy(1,:)-shift;
Vz(1,:) = Vz(1,:)-shift;

som = sum(Vx(1,:))+sum(Vy(1,:))+sum(Vz(1,:));

temp = 1/(3*Np^3)*sum(Vx(1,:).*Vx(1,:)+Vy(1,:).*Vy(1,:)+Vz(1,:).*Vz(1,:));
Vx(1,:) = Vx(1,:)/sqrt(temp/T0);
Vz(1,:) = Vy(1,:)/sqrt(temp/T0);
Vz(1,:) = Vz(1,:)/sqrt(temp/T0);
temp2 = 1/(3*Np^3)*sum(Vx(1,:).*Vx(1,:)+Vy(1,:).*Vy(1,:)+Vz(1,:).*Vz(1,:));
K(1) = 0.5*m*sum(Vx(1,:).*Vx(1,:)+Vy(1,:).*Vy(1,:)+Vz(1,:).*Vz(1,:));
E(1) = K(1);
delta = zeros(N,Np^3,Np^3);
%% Main loop
teller = 0;
h = waitbar(0,'Progress 0%%', 'Name', 'Time evolution');
for ii=2:N

      % Step 1 of velocity-verlet
      x(ii,:) = x(ii-1,:)+Vx(ii-1,:)*Dt+(Ax(1,:)/2)*Dt^2;
      y(ii,:) = y(ii-1,:)+Vy(ii-1,:)*Dt+(Ay(1,:)/2)*Dt^2;
      z(ii,:) = z(ii-1,:)+Vz(ii-1,:)*Dt+(Az(1,:)/2)*Dt^2;

      % Apply PBC

      x(ii,:) = x(ii,:)-L*floor(x(ii,:)/L);
      y(ii,:) = y(ii,:)-L*floor(y(ii,:)/L);
      z(ii,:) = z(ii,:)-L*floor(z(ii,:)/L);

      % Intermediate update of the velocities
      Vx(ii,:) = Vx(ii-1,:)+Dt*Ax(1,:)/2;
      Vy(ii,:) = Vy(ii-1,:)+Dt*Ay(1,:)/2;
      Vz(ii,:) = Vz(ii-1,:)+Dt*Az(1,:)/2;

      % Calculate new accelerations

      % Reset accelerations
      Ax = zeros(1,Np^3);
      Ay = Ax;
      Az = Ax;
      scale = zeros(N,1);
      % Interactions
      for jj = 1:Np^3
        for kk = jj+1:Np^3

                deltaX = x(ii,jj)-x(ii,kk);
                deltaY = y(ii,jj)-y(ii,kk);
                deltaZ = z(ii,jj)-z(ii,kk);

                % Second part of PBC
                deltaX = deltaX - round(deltaX/L)*L;
                deltaY = deltaY - round(deltaY/L)*L;
                deltaZ = deltaZ - round(deltaZ/L)*L;


                r = sqrt(deltaX^2+deltaY^2+deltaZ^2);
                delta(ii,jj,kk) = r;
                if  r < Rc
                    Ax(1,jj) = Ax(1,jj)+24*eps/m*(-sigma^6/r^8+2*sigma^12/r^14)*(deltaX);
                    Ay(1,jj) = Ay(1,jj)+24*eps/m*(-sigma^6/r^8+2*sigma^12/r^14)*(deltaY);
                    Az(1,jj) = Az(1,jj)+24*eps/m*(-sigma^6/r^8+2*sigma^12/r^14)*(deltaZ);

                    Ax(1,kk) = Ax(1,kk)-24*eps/m*(-sigma^6/r^8+2*sigma^12/r^14)*(deltaX);
                    Ay(1,kk) = Ay(1,kk)-24*eps/m*(-sigma^6/r^8+2*sigma^12/r^14)*(deltaY);
                    Az(1,kk) = Az(1,kk)-24*eps/m*(-sigma^6/r^8+2*sigma^12/r^14)*(deltaZ);

                    P(ii) = P(ii)+ (4*eps*((sigma/r)^12-(sigma/r)^6)-4*eps*((sigma/Rc)^12-(sigma/Rc)^6));
                end

        end
      end

      % Second update of velocity

      Vx(ii,:) = Vx(ii,:)+Ax(1,:)*Dt/2;
      Vy(ii,:) = Vy(ii,:)+Ay(1,:)*Dt/2;
      Vz(ii,:) = Vz(ii,:)+Az(1,:)*Dt/2;

      if Vx(ii,:)~=Vx(ii-1,:)
          teller = teller+1;
      end
      waitbar(ii/N,h,sprintf('Progress %d%%',round(ii/N*100)));
      K(ii) = 0.5*m*sum(Vx(ii,:).*Vx(ii,:)+Vy(ii,:).*Vy(ii,:)+Vz(ii,:).*Vz(ii,:),2);


      T(ii) = K(ii)/(3*Np^3); % Kinetic temperature before rescaling
      % rescaling
      if state == 0
          scale(ii) = 1;
      else
          scale(ii) = sqrt( 1 + 1/tau*(T0/T(ii)-1));
      end
      Vx(ii,:) = scale(ii)*Vx(ii,:);
      Vy(ii,:) = scale(ii)*Vy(ii,:);
      Vz(ii,:) = scale(ii)*Vz(ii,:);
      % recalculate Kinetic energy
      K(ii) = 0.5*m*sum(Vx(ii,:).*Vx(ii,:)+Vy(ii,:).*Vy(ii,:)+Vz(ii,:).*Vz(ii,:),2);

      E(ii) = K(ii)+P(ii);

end
close(h);
%% plotting
figure(1);clf;hold on;view(3)
scatter3(x(1,:),y(1,:),z(1,:)) % Show initial positions, the cubic grid at the center
scatter3(x(10,:),y(10,:),z(10,:),'fill') % Show positions at the end of the simulation
scatter3([0 0 0 0 L L 0 L L L/2],[0 L L 0 0 0 0 L L L/2], [0 0 L L 0 L L 0 L L/2],100, 'fill'); % Show corners of cubic volume
%axis([0 L 0 L 0 L])

%% Conservation of energy


figure(2);clf;hold on;
time = linspace(0,Dt*N,N);

subplot(4,1,1)
plot(time,E, '.r');
title('Total energy');
subplot(4,1,2)
plot(time,P, '.-b');
title('Scale factor');
subplot(4,1,3)
plot(time,K, 'm');
title('Kinetic energy');
subplot(4,1,4)
plot(time,T);
title('Temperature')

%% Movement

if movie > 0
figure(3)
hold off;

h = waitbar(0,'Progress 0%%', 'Name', 'Time evolution');
for ii = 1:N
    grid off;
    scatter3([0 0 0 0 L L 0 L L L/2],[0 L L 0 0 0 0 L L L/2], [0 0 L L 0 L L 0 L L/2],100, 'fill');
    scatter3(x(ii,:),y(ii,:),z(ii,:),'fill');
    axis([0 L 0 L 0 L])
    title('Lennard-Jones liquid');
    M(ii) = getframe;
    waitbar(ii/N,h,sprintf('Progress %d%%',round(ii/N*100)));
end
close(h);
end