% ********************************************************************************
% ** FICHE F.3.  LOW-STORAGE MD PROGRAMS USING THE LEAPFROG VERLET ALGORITHM    **
% ** This FORTRAN code is intended to illustrate points made in the text.       **
% ** To our knowledge it works correctly.  However it is the responsibility of  **
% ** the user to test it, if it is to be used in a research application.        **
% ********************************************************************************

%     *******************************************************************
%     ** TWO SEPARATE PARTS: FORTRAN AND BASIC VERSIONS.               **
%     *******************************************************************

%     *******************************************************************
%     ** FICHE F.3 - PART A                                            **
%     ** FORTRAN PROGRAM USING THE LEAPFROG ALGORITHM.                 **
%     *******************************************************************

%     *******************************************************************
%     ** FORTRAN PROGRAM TO CONDUCT MOLECULAR DYNAMICS OF ATOMS.       **
%     **                                                               **
%     ** A SPECIAL LOW-STORAGE FORM OF THE LEAPFROG VERLET ALGORITHM   **
%     ** IS USED AND WE TAKE THE LENNARD-JONES POTENTIAL AS AN EXAMPLE.**
%     **                                                               **
%     ** REFERENCE:                                                    **
%     **                                                               **
%     ** FINCHAM AND HEYES, CCP5 QUARTERLY, 6, 4, 1982.                **
%     **                                                               **
%     ** ROUTINES REFERENCED:                                          **
%     **                                                               **
%     ** SUBROUTINE READCN ( CNFILE )                                  **
%     **    READS IN CONFIGURATION                                     **
%     ** SUBROUTINE FORCE ( DT, SIGMA, RCUT, NEWV, NEWVC, NEWW )       **
%     **    CALCULATES THE ACCELERATIONS AND ADVANCES THE VELOCITIES   **
%     **    FROM T - DT/2 TO T + DT/2. ALSO CALCULATES POTENTIAL       **
%     **    ENERGY AND VIRIAL AT TIMESTEP T.                           **
%     ** SUBROUTINE MOVE ( DT )                                        **
%     **    ADVANCES POSITIONS FROM T TO T + DT.                       **
%     ** SUBROUTINE KINET ( NEWK )                                     **
%     **    CALCULATES KINETIC ENERGY.                                 **
%     ** SUBROUTINE WRITCN ( CNFILE )                                  **
%     **    WRITES OUT CONFIGURATION                                   **
%     **                                                               **
%     ** PRINCIPAL VARIABLES:                                          **
%     **                                                               **
%     ** INTEGER N                 NUMBER OF ATOMS                     **
%     ** REAL    DT                TIMESTEP                            **
%     ** REAL    RX(N),RY(N),RZ(N) ATOMIC POSITIONS                    **
%     ** REAL    VX(N),VY(N),VZ(N) ATOMIC VELOCITIES                   **
%     ** REAL    ACV,ACK ETC.      AVERAGE VALUE ACCUMULATORS          **
%     ** REAL    AVV,AVK ETC.      AVERAGE VALUES                      **
%     ** REAL    ACVSQ, ACKSQ ETC. AVERAGE SQUARED VALUE ACCUMULATORS  **
%     ** REAL    FLV,FLK ETC.      FLUCTUATION AVERAGES                **
%     **                                                               **
%     ** USAGE:                                                        **
%     **                                                               **
%     ** THE LEAPFROG ALGORITHM IS OF THE FORM                         **
%     ** VX(T + DT/2) = VX(T - DT/2) + DT * AX(T)  (SIMILARLY FOR Y,Z) **
%     ** RX(T + DT)   = RX(T) + DT * VX(T + DT/2)  (SIMILARLY FOR Y,Z) **
%     ** TO SAVE STORAGE IN THIS PROGRAM WE ACCUMULATE VALUES AX(T)    **
%     ** DIRECTLY ONTO THE VELOCITIES VX(T - DT/2) IN THE FORCE LOOP.  **
%     ** THIS MEANS THAT AN APPROXIMATION MUST BE USED FOR THE KINETIC **
%     ** ENERGY AT EACH TIME STEP:                                     **
%     ** K = ( OLDK + NEWK ) / 2 + ( OLDV - 2 * V + NEWV ) / 8         **
%     ** WHERE K    = K( T        ), V    = V( T )                     **
%     **       OLDK = K( T - DT/2 ), OLDV = V( T - DT )                **
%     **       NEWK = K( T + DT/2 ), NEWV = V( T + DT )                **
%     ** AT THE START OF A STEP THE FOLLOWING VARIABLES ARE STORED:    **
%     ** R    : R(STEP)         V    : V(STEP+1/2)                     **
%     ** OLDV : V(STEP-2)       V    : V(STEP-1)        NEWV : V(STEP) **
%     ** OLDK : K(STEP-1/2)     NEWK : K(STEP+1/2)                     **
%     ** THIS PROGRAM USES UNIT 10 FOR CONFIGURATION INPUT AND OUTPUT  **
%     **                                                               **
%     ** UNITS:                                                        **
%     **                                                               **
%     ** THE PROGRAM USES LENNARD-JONES REDUCED UNITS FOR USER INPUT   **
%     ** AND OUTPUT BUT CONDUCTS SIMULATION IN A BOX OF UNIT LENGTH.   **
%     ** SUMMARY FOR BOX LENGTH L, ATOMIC MASS M, AND LENNARD-JONES    **
%     ** POTENTIAL PARAMETERS SIGMA AND EPSILON:                       **
%     **                OUR PROGRAM           LENNARD-JONES SYSTEM     **
%     ** LENGTH         L                     SIGMA                    **
%     ** MASS           M                     M                        **
%     ** ENERGY         EPSILON               EPSILON                  **
%     ** TIME           SQRT(M*L**2/EPSILON)  SQRT(M*SIGMA**2/EPSILON) **
%     ** VELOCITY       SQRT(EPSILON/M)       SQRT(EPSILON/M)          **
%     ** PRESSURE       EPSILON/L**3          EPSILON/SIGMA**3         **
%     *******************************************************************
% ***********************************************************************
fprintf('Running the program. Please be patient.\n');
% ***********************************************************************
% Declaring set input variables for use. (These are the variables that you can change)
NSTEP = 200;       % Number of steps. (At least 200 or the system won't restore to equilibrium)
IPRINT = 5;        % Frequency of outputting data. (Eg: 5 means output data every 5 steps). 
                   % This will decide the frequency of the data shown on the screen and saved)
DENS = 0.9;        % Density (0.7-0.8 for liquid)
TEMPINPUT = 2.8;   % Temperature Input. (The density and temperature should match, 
                   % I would like to recommend Dens=0.9 and T=2.8)
RCUT = 3;          % Potential cut off distance (2.5 - 5 as you like)
DT = 0.01;         % Time step.

% Declaring and initializing some variables for use by the program. (You should not touch)  
N = 500;           % Number of atoms.
FREE = 3.0;
RX = zeros(N,1);   % Atomic position for x-coord.
RY = zeros(N,1);   % Atomic position for y-coord.
RZ = zeros(N,1);   % Atomic position for z-coord.
VX = zeros(N,1);   % Atomic velocity for x-dir.
VY = zeros(N,1);   % Atomic velocity for y-dir.
VZ = zeros(N,1);   % Atomic velocity for z-dir.

% Setting up tracking arrays.
PX = zeros(NSTEP);
PY = zeros(NSTEP);
PZ = zeros(NSTEP);

% Setting up the velocity distribution array.
VDist = zeros(20,1);

% Initializing the Velocity vectors.
for i = 1:N
  VX(i,1)=2*(rand-0.5);
  VY(i,1)=2*(rand-0.5);
  VZ(i,1)=2*(rand-0.5);
end;

% Read configuration into common / Block1 / Variables.
NC = (N/4)^(1/3);
XL = 1.0/NC;
Y = .5*XL;

RX(1,1) = 0.0;
RY(1,1) = 0.0;
RZ(1,1) = 0.0;

RX(2,1) = 0.0;
RY(2,1) = Y;
RZ(2,1) = Y;

RX(3,1) = Y;
RY(3,1) = 0.0;
RZ(3,1) = Y;

RX(4,1) = Y;
RY(4,1) = Y;
RZ(4,1) = 0;

M = 0;
for i = 1:NC
  for j = 1:NC
    for t = 1:NC
      for ij = 1:4
        RX(ij+M,1) = RX(ij,1)+XL*(t-1);
	    RY(ij+M,1) = RY(ij,1)+XL*(j-1);
	    RZ(ij+M,1) = RZ(ij,1)+XL*(i-1);
      end;
      M = M + 4;
    end;
  end;
end;

% Convert input data into program units.
SIGMA = (DENS/N)^(1/3);
RCUT  = RCUT * SIGMA;
DT    = DT * SIGMA;
DENS  = DENS / (SIGMA^3);

% Calculate Long range corrections.
% NOTE: Specific to Lennard-Jones.
SR3    = (SIGMA/RCUT)^3;
SR9    = SR3^3;
SIGCUB = SIGMA^3;
VLRC = (8/9)*pi*DENS*SIGCUB*N*(SR9-3.0*SR3);
WLRC = (16/9)*pi*DENS*SIGCUB*N*(2*SR9-3.0*SR3);

% Zero Accumulators.
ACV  = 0.0;
ACK  = 0.0;
ACE  = 0.0;
ACEC = 0.0;
ACP  = 0.0;
ACT  = 0.0;

ACVSQ  = 0.0;
ACKSQ  = 0.0;
ACESQ  = 0.0;
ACECSQ = 0.0;
ACPSQ  = 0.0;
ACTSQ  = 0.0;

FLV  = 0.0;
FLK  = 0.0;
FLE  = 0.0;
FLEC = 0.0;
FLP  = 0.0;
FLT  = 0.0;

% Calculate Initial Values.
[VX,VY,VZ,NEWV,NEWVC,NEWW] = Force(-DT, SIGMA, RCUT, VX, VY, VZ, RX, RY, RZ, N);
[RX,RY,RZ] = Move(RX, RY, RZ, VX, VY, VZ, -DT, N);
[VX,VY,VZ,V,VC,W] = Force(-DT, SIGMA, RCUT, VX, VY, VZ, RX, RY, RZ, N);  
[VX,VY,VZ,V,VC,W] = Force(DT, SIGMA, RCUT, VX, VY, VZ, RX, RY, RZ, N);
OLDK = Kinet(VX, VY, VZ, N);
[RX,RY,RZ] = Move(RX, RY, RZ, VX, VY, VZ, DT, N);
[VX,VY,VZ,NEWV,NEWVC,NEWW] = Force(DT, SIGMA, RCUT, VX, VY, VZ, RX, RY, RZ, N);
NEWK = Kinet(VX, VY, VZ, N);

% Include long range corrections.
V = V + VLRC;
W = W + WLRC;
NEWV = NEWV + VLRC;
NEWW = NEWW + WLRC;
if (IPRINT <= 0) IPRINT = NSTEP + 1; end;

% Setting up output file procedures.
posFile = fopen('position.xyz', 'w');
kinFile = fopen('kinetic.txt', 'w');
velFile = fopen('velocity.txt', 'w');
preFile = fopen('pressure.txt', 'w');
temFile = fopen('temperature.txt', 'w');
eneFile = fopen('energy.txt', 'w');

% ********************************************************************************
% MAIN LOOP BEGINS.
% ********************************************************************************
for STEP = 1:NSTEP
  % Implementing algorithm.
  [RX,RY,RZ] = Move(RX, RY, RZ, VX, VY, VZ, DT, N);
  OLDV = V;
  V = NEWV;
  OLDVC = VC;
  VC = NEWVC;
  W = NEWW;
  [VX,VY,VZ,NEWV,NEWVC,NEWW] = Force(DT, SIGMA, RCUT, VX, VY, VZ, RX, RY, RZ, N);

  % Include Long-range corrections.
  NEWV = NEWV + VLRC;
  NEWW = NEWW + WLRC;

  % Calculate kinetic energy at current step.
  K = (NEWK+OLDK) / 2.0 + (NEWV-2.0*V+OLDV) / 8.0;
  KC = (NEWK+OLDK) / 2.0 + (NEWVC-2.0*VC+OLDVC) / 8.0;
  OLDK = NEWK;
  NEWK = Kinet(VX, VY, VZ, N);

  % Calculate instanttaneous values.
  E    = K+V;
  EC   = KC+VC;
  VN   = V/N;
  KN   = K/N;
  EN   = E/N;
  ECN  = EC/N;
  TEMP = 2.0 * KN/FREE;
  PRES = DENS*TEMP+W;

  % Convert to Lennard-Jones Units.
  PRES = PRES * SIGMA^3;

  % Increment Accumulators.
  ACE  = ACE  + EN;
  ACEC = ACEC + ECN;
  ACK  = ACK  + KN;
  ACV  = ACV  + VN;
  ACP  = ACP  + PRES;

  ACESQ  = ACESQ  + EN^2; 
  ACECSQ = ACECSQ + ECN^2;
  ACKSQ  = ACKSQ  + KN^2;
  ACVSQ  = ACVSQ  + VN^2;
  ACPSQ  = ACPSQ  + PRES^2;
  
  % Assigning the trajectory for one partcle.
  PX(STEP) = RX(244,1);
  PY(STEP) = RY(244,1);
  PZ(STEP) = RZ(244,1);

  % Setting up the output data.
  if (mod(STEP,IPRINT) == 0)
    % Output the position in a xyz format.
    fprintf(posFile, '%d\n%d\n', N, STEP);
    for i = 1:N
      fprintf(posFile, '%s\t %g\t %g\t %g\n', 'Ar', RX(i,1)*40, RY(i,1)*40, RZ(i,1)*40);
    end;
    fprintf(posFile, '\n');

    % Output KE.
    fprintf(kinFile, '%g\t %g\n', STEP, KN);

    % Output velocity.
    fprintf(velFile, '%g\t %g\n', STEP, VN);

    % Output pressure.
    fprintf(preFile, '%g\t %g\n', STEP, PRES);  

    % Output temperature.
    fprintf(temFile, '%g\t %g\n', STEP, TEMP);

    % Output energy.
    fprintf(eneFile, '%g\t %g\n', STEP, EN);
  end;
 
  BEITA = (TEMPINPUT/TEMP)^(1/2);
  for i = 1:N
    VX(i,1) = BEITA*VX(i,1);
    VY(i,1) = BEITA*VY(i,1);
    VZ(i,1) = BEITA*VZ(i,1);
  end;
end;   
% ********************************************************************************
% MAIN LOOP ENDS.
% ********************************************************************************
% Closing all output files.
fclose(posFile);
fclose(kinFile);
fclose(velFile);
fclose(preFile);
fclose(temFile);
fclose(eneFile);

% Output run averages.
NORM = NSTEP;
AVE  = ACE  / NORM;
AVEC = ACEC / NORM;
AVK  = ACK  / NORM;
AVV  = ACV  / NORM;
AVP  = ACP  / NORM;

ACESQ  = (ACESQ/NORM) - AVE^2;
ACECSQ = (ACECSQ/NORM) - AVEC^2;
ACKSQ  = (ACKSQ/NORM) - AVK^2;
ACVSQ  = (ACVSQ/NORM) - AVV^2;
ACPSQ  = (ACPSQ/NORM) - AVP^2;

if (ACESQ > 0.0 ) FLE  = sqrt(ACESQ); end;
if (ACECSQ > 0.0 ) FLEC = sqrt(ACECSQ); end;
if (ACKSQ > 0.0 ) FLK  = sqrt(ACKSQ); end;
if (ACVSQ > 0.0 ) FLV  = sqrt(ACVSQ); end;
if (ACPSQ > 0.0 ) FLP  = sqrt(ACPSQ); end;

AVT = AVK * 2.0 / FREE;
FLT = FLK * 2.0 / FREE;

% Outputting the averages.
aveFile = fopen('average.txt', 'w');
fprintf(aveFile, '%g\t %g\t %g\t %g\t %g\t %g\n', AVE, AVEC, AVK, AVV, AVP, AVT);
fclose(aveFile);        

% Outputting the fluctuations.
fluFile = fopen('fluctuations.txt', 'w');
fprintf(aveFile, '%g\t %g\t %g\t %g\t %g\t %g\n', FLE, FLEC, FLK, FLV, FLP, FLT);
fclose(fluFile);    

% Assigning the velocity distruibution array for X Direction.
for i = 1:N
  if (VX(i,1) <= -4.5) VDist(1) = VDist(1)+1;
  elseif (VX(i,1) > -4.5 && VX(i,1) <= -4) VDist(2) = VDist(2)+1;
  elseif (VX(i,1) > -4 && VX(i,1) <= -3.5) VDist(3) = VDist(3)+1;
  elseif (VX(i,1) > -3.5 && VX(i,1) <= -3) VDist(4) = VDist(4)+1;
  elseif (VX(i,1) > -3 && VX(i,1) <= -2.5) VDist(5) = VDist(5)+1;
  elseif (VX(i,1) > -2.5 && VX(i,1) <= -2) VDist(6) = VDist(6)+1;
  elseif (VX(i,1) > -2 && VX(i,1) <= -1.5) VDist(7) = VDist(7)+1;
  elseif (VX(i,1) > -1.5 && VX(i,1) <= -1) VDist(8) = VDist(8)+1;
  elseif (VX(i,1) > -1 && VX(i,1) <= -0.5) VDist(9) = VDist(9)+1;
  elseif (VX(i,1) > -0.5 && VX(i,1) <= 0) VDist(10) = VDist(10)+1;
  elseif (VX(i,1) > 0 && VX(i,1) <= 0.5) VDist(11) = VDist(11)+1;
  elseif (VX(i,1) > 0.5 && VX(i,1) <= 1) VDist(12) = VDist(12)+1;
  elseif (VX(i,1) > 1 && VX(i,1) <= 1.5) VDist(13) = VDist(13)+1;
  elseif (VX(i,1) > 1.5 && VX(i,1) <= 2) VDist(14) = VDist(14)+1;
  elseif (VX(i,1) > 2 && VX(i,1) <= 2.5) VDist(15) = VDist(15)+1;
  elseif (VX(i,1) > 2.5 && VX(i,1) <= 3) VDist(16) = VDist(16)+1;
  elseif (VX(i,1) > 3 && VX(i,1) <= 3.5) VDist(17) = VDist(17)+1;
  elseif (VX(i,1) > 3.5 && VX(i,1) <= 4) VDist(18) = VDist(18)+1;
  elseif (VX(i,1) > 4 && VX(i,1) <= 4.5) VDist(19) = VDist(19)+1;
  elseif (VX(i,1) > 4.5) VDist(20) = VDist(20)+1;
  end;
end;

% Plot the trajectory.
figure(1);
set(gca,'DataAspectRatio',[1 1 1]);
plot3(PX,PY,PZ, '.', PX,PY,PZ);
zoom(1.5);
grid on;
title('The trajectory of one particle');
xlabel('X');
ylabel('Y');
zlabel('Y');

% Plot the bar graph.
figure(2)
bar(VDist,'r');
title('Velocity Distribution for all particles at a time instance');
xlabel('Each column represents a bin size of 0.5 ranging from -4.5 to 4.5');
ylabel('The distribution of each bin');

% Signalling the end of the program.
fprintf('DONE\n');
% ********************************************************************************
% ********************************************************************************




