* molec_dyn.f
*   Molecular dynamics for n particles, L-J potential, 3 dimensions.
*      Verlet's algorithm used to integrate equations of motion.
*      Reads data from input file md.in and writes results to md.out.
* Edited version omits date information
***********************************************************************
      program MD_LJ 
      IMPLICIT NONE
      INTEGER B1, B2, I, ITMP, K, NATOM, NATOMX, NSTEP, PRNFREQ
      PARAMETER (NATOMX = 50)

*****************************************************************************************
* B1,B2 - Buffer indices.
* ITMP - Temporary variable, for switching buffer index values.
* I - Loop index.
* K - Current time step.  (OUT)
* NATOM - number of atoms in chain. (IN)
* NATOMX - Maximum number of atoms permitted in chain.
* NSTEP - Number of time steps in computation.
* PRNFREQ - Number of time steps between printing.
******************************************************************************************

      DOUBLE PRECISION H, HO2, HP, PE, KE, TEND,
     $     X(NATOMX), Y(NATOMX), Z(NATOMX),
     $     U(NATOMX), V(NATOMX), W(NATOMX),
     $     FX(NATOMX,2), FY(NATOMX,2), FZ(NATOMX,2),
     $     POTENTIAL_LJ

           OPEN (UNIT = 2, FILE = 'md.in')
           OPEN (UNIT = 3, FILE = 'md.out')

********************************************************************************************
* H - Time step. (IN)
* HO2 - h/2.
* HP - Time interval between printing.   (IN)
* PE - Potential energy of the system.   (OUT)
* KE - Kinetic energy of the system.   (OUT)
* TEND - Time at end of computation: TEND = H*NSTEP.  (IN)
* FX(),FY(),FZ() - Force vectors in x,y,z directions.
* X(),Y(),Z() - Position vectors in x,y,z directions.  (IN/OUT)
* U(),V(),W() - Velocity vectors in x,y,z directions.  (IN/OUT)
* POTENTIAL_LJ - Function to compute LJ potential.
**********************************************************************************************

* While not done --
  10   CONTINUE
 
* Read computation parameters
      READ(2,*,END=999) NATOM, H, HP, TEND

* Write computation parameters
       WRITE(3, '(I8)') NATOM
       WRITE(3, '(I8)') INT(TEND/HP)

* Read initial positions, velocities
      DO  I = 1,NATOM
         READ(2,*) X(I), Y(I), Z(I)
         READ(2,*) U(I), V(I), W(I)
      END DO

* Set initial buffer indices
      B1 = 1
      B2 = 2

* Set parameter values
      NSTEP = NINT(TEND/H)
      PRNFREQ = NINT(HP/H)
      HO2 = H/2.0

* Initial forces
      CALL FORCE_LJ(NATOM,NATOMX,X,Y,Z,B1,FX,FY,FZ)

* Print initial coordinates
!      WRITE(*,*) 'Time   0   KineticEnergy  PotentlEnergy  TotalEnergy'
!      WRITE(*,*) 'Time  Atom  X-Coordinate  Y-Coordinate  Z-Coordinate'
!      WRITE(*,*) 'Time  Atom   X-Velocity    Y-Velocity    Z-Velocity'
!      WRITE(*,*) ' '
      K = 0
      PE = POTENTIAL_LJ(NATOM,X,Y,Z)
      KE = 0
      DO I = 1,NATOM
         KE = KE + (U(I)**2 + V(I)**2 + W(I)**2)/2
      END DO
      WRITE(3,99) K*0.0, 0, PE, KE, PE+KE
      DO I = 1,NATOM

         WRITE(3,99) K*H, I, X(I), Y(I), Z(I)
         WRITE(3,99) K*H, I, U(I), V(I), W(I)
      END DO

* Main loop
      DO K = 1, NSTEP
         DO I = 1,NATOM
            X(I) = X(I) + H*(U(I) + HO2*FX(I,B1))
            Y(I) = Y(I) + H*(V(I) + HO2*FY(I,B1))
            Z(I) = Z(I) + H*(W(I) + HO2*FZ(I,B1))
         END DO
         CALL FORCE_LJ(NATOM,NATOMX,X,Y,Z,B2,FX,FY,FZ)
         DO I = 1,NATOM
            U(I) = U(I) + HO2*(FX(I,B1) + FX(I,B2))
            V(I) = V(I) + HO2*(FY(I,B1) + FY(I,B2))
            W(I) = W(I) + HO2*(FZ(I,B1) + FZ(I,B2))
         END DO

         IF (MOD(K,PRNFREQ) .EQ. 0) THEN
            PE = POTENTIAL_LJ(NATOM,X,Y,Z)
            KE = 0
            DO I = 1,NATOM
               KE = KE + (U(I)**2 + V(I)**2 + W(I)**2)/2
            END DO
            WRITE(3,99)INT(1000.*K*H),0,PE,KE,PE+KE
            DO I = 1,NATOM
               WRITE(3,99)K*H,I,X(I),Y(I),Z(I)
               WRITE(3,99)K*H,I,U(I),V(I),W(I)
            END DO
         END IF

         ITMP = B1
         B1 = B2
         B2 = ITMP
      END DO

      GO TO 10

 999  CONTINUE
      STOP

 98   FORMAT (I5, 3(1PE13.6,1X))
 99   FORMAT (F13.6, 1X, I5, 1X, 3(1PE13.6,1X))

      END
*********************************************************************************************
*potential_lj_ymp.f
*    Compute potential energy from Lennard Jones interaction.
*        V(r) = (1/r**12 - 2/r**6)

*****************************************************************************************
      DOUBLE PRECISION FUNCTION POTENTIAL_LJ(NATOM, X, Y, Z)

      IMPLICIT NONE
      INTEGER I, J, NATOM
*****************************************************************************************
* I,J - Loop indices.
* NATOM - Number of atoms in chain.
*****************************************************************************************

      DOUBLE PRECISION X(*), Y(*), Z(*), POT, RIJSQ, RIJM6
*****************************************************************************************
* X(),Y(),Z() - Position vectors in x,y,z directions.
* POT - Potential energy.
* RIJSQ - Reciprocal of distance between two atoms, squared.
* RIJM6 - Reciprocal of distance between two atoms,
*               to the sixth power.
*********************************************************************

* Initialize potential.
      POT = 0.0

* Compute potential.
      DO I = 1, NATOM-1
         DO J = I+1, NATOM
            RIJSQ = 1.0/((X(I)-X(J))**2 + (Y(I)-Y(J))**2 +
     $           (Z(I)-Z(J))**2)
            RIJM6 = RIJSQ*RIJSQ*RIJSQ
            POT = POT + RIJM6*(RIJM6 - 2.0)
                    END DO
      END DO

* Potential computed.
      POTENTIAL_LJ = POT
      RETURN
      END
*******************************************************************************************
*   force_lj.f
*    Compute forces from Lennard Jones interaction.
*        V(r) = (1/r**12 - 2/r**6)

*************************************************************************************

      SUBROUTINE FORCE_LJ(NATOM,NATOMX,X,Y,Z,B,FX,FY,FZ)

      IMPLICIT NONE
      INTEGER I, J, NATOM, NATOMX, B
      INTEGER NATOMXX
      PARAMETER (NATOMXX = 50)
*************************************************************************************
* I,J - Loop indices.
* NATOM - Number of atoms in chain.
* NATOMX - Maximum number of atoms permitted in chain.
* B - Buffer index.
*************************************************************************************

      DOUBLE PRECISION X(*), Y(*), Z(*),
     $     FX(NATOMX,*), FY(NATOMX,*), FZ(NATOMX,*),
     $     RIJSQ, RIJM6, CIJ, FIJX(NATOMXX),
     $     FIJY(NATOMXX), FIJZ(NATOMXX)
*************************************************************************************
* X(),Y(),Z() - Position vectors in x,y,z directions.
* FX(),FY(),FZ() - Force vectors in x,y,z directions.
* RIJSQ - Reciprocal of distance between two atoms, squared.
* RIJM6 - Reciprocal of distance between two atoms,
*               to the sixth power.
* CIJ - Common subexpression for forces.
* FIJX,FIJY,FIJZ - Force between two atoms in x,y,z directions.
**************************************************************************************

* Initialize forces.
      DO I = 1, NATOM
         FX(I,B) = 0.0
         FY(I,B) = 0.0
         FZ(I,B) = 0.0
      END DO

* Compute forces.
      DO I = 1, NATOM-1
         DO J = I+1, NATOM
            RIJSQ = 1.0/((X(I)-X(J))**2 + (Y(I)-Y(J))**2 +
     $           (Z(I)-Z(J))**2)
            RIJM6 = RIJSQ*RIJSQ*RIJSQ
            CIJ = 12.0*(RIJM6 - 1.0)*RIJM6*RIJSQ
            FIJX(J) = CIJ*(X(I)-X(J))
            FIJY(J) = CIJ*(Y(I)-Y(J))
            FIJZ(J) = CIJ*(Z(I)-Z(J))
            FX(I,B) = FX(I,B) + FIJX(J)
            FY(I,B) = FY(I,B) + FIJY(J)
            FZ(I,B) = FZ(I,B) + FIJZ(J)
         END DO
         DO J = I+1, NATOM
            FX(J,B) = FX(J,B) - FIJX(J)
            FY(J,B) = FY(J,B) - FIJY(J)
            FZ(J,B) = FZ(J,B) - FIJZ(J)
         END DO
      END DO

* Forces computed.
      RETURN
      END




