%   *******************************************************************
%   ** ROUTINE TO COMPUTE FORCES WITH CUTOFF & MINIMUM IMAGING.      **
%   **                                                               **
%   ** THE ROUTINE ACTUALLY RETURNS TWO DIFFERENT POTENTIAL ENERGIES.**
%   ** V IS CALCULATED USING THE UNSHIFTED LENNARD-JONES POTENTIAL.  **
%   ** WHEN LONG-RANGE TAIL CORRECTIONS ARE ADDED, THIS MAY BE USED  **
%   ** TO CALCULATE THERMODYNAMIC INTERNAL ENERGY ETC.               **
%   ** VC IS CALCULATED USING THE SHIFTED LENNARD-JONES POTENTIAL    **
%   ** WITH NO DISCONTINUITY AT THE CUTOFF.  THIS MAY BE USED TO     **
%   ** CHECK ENERGY CONSERVATION.                                    **
%   ** PROGRAM UNITS MAKE EPSILON = MASS = 1, BUT SIGMA IS NOT UNITY **
%   ** SINCE WE TAKE A BOX OF UNIT LENGTH.                           **
%   ** THE ROUTINE IS COMBINED WITH THE LEAPFROG ALGORITHM FOR       **
%   ** ADVANCING VELOCITIES, I.E. FORCES ARE ACCUMULATED DIRECTLY    **
%   ** ONTO VELOCITIES.                                              **
%   *******************************************************************
function [VX,VY,VZ,V,VC,W] = Force(DT, SIGMA, RCUT, VX, VY, VZ, RX, RY, RZ, N);

% Calculate square distance for inner loop.
SIGSQ  = SIGMA^2;
RCUTSQ = RCUT^2;

% Conduct outer loop over atoms.
V = 0.0;
W = 0.0;
NCUT = 0;

for i = 1:N-1        
  RXI = RX(i,1);
  RYI = RY(i,1);
  RZI = RZ(i,1);
  VXI = VX(i,1);
  VYI = VY(i,1);
  VZI = VZ(i,1);

  % Conduct inner loop over atoms.
  for j = i+1:N
    RXIJ = RXI - RX(j,1);
    RXIJ = RXIJ - round(RXIJ);

    if (abs(RXIJ) < RCUT)
      RYIJ = RYI - RY(j,1);
      RYIJ = RYIJ - round(RYIJ);
      RIJSQ = RXIJ^2 + RYIJ^2;
    
      if (RIJSQ < RCUTSQ)
        RZIJ = RZI - RZ(j,1);
        RZIJ = RZIJ - round(RZIJ);
        RIJSQ = RIJSQ + RZIJ^2;

        if (RIJSQ < RCUTSQ)
          SR2   = SIGSQ / RIJSQ;
          SR6   = SR2^3;
          SR12  = SR6^2;
          VIJ   = 4.0 * (SR12 - SR6);
          WIJ   = 24.0 * (2.0 * SR12 - SR6);
          VELIJ = WIJ * DT / RIJSQ;
          DVX   = VELIJ * RXIJ;
          DVY   = VELIJ * RYIJ;
          DVZ   = VELIJ * RZIJ;
          VXI   = VXI + DVX;
          VYI   = VYI + DVY;
          VZI   = VZI + DVZ;
          VX(j,1) = VX(j,1) - DVX;
          VY(j,1) = VY(j,1) - DVY;
          VZ(j,1) = VZ(j,1) - DVZ;
          V     = V + VIJ;
          W     = W + WIJ;
          NCUT  = NCUT + 1;
        end;
      end;
    end;
  end;   % End of inner Loop.

  VX(i,1) = VXI;
  VY(i,1) = VYI;
  VZ(i,1) = VZI;
end;     % End of outer Loop.

% Calculate Potential Shift.
SR2 = SIGSQ / RCUTSQ;
SR6 = SR2 * SR2 * SR2;
SR12 = SR6 * SR6;
VIJ = 4.0 * (SR12 - SR6);
VC = V - NCUT * VIJ;
W = W / 3.0;
