Here is the Fortran77 Program for Molecular Dynamics Excersise
Subroutines for Potential Energy and Force are included
<![CDATA[
* 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
***********************************************************************
* This work has been supported by the National Science Foundation
* under an Educational Infrastructure grant, CDA-9017953. It has
* been produced by the HPSC Group, Department of Computer Science,
* University of Colorado, Boulder, CO 80309. Please direct
* comments or queries to Elizabeth Jessup at this address or e-mail
* jessup@cs.colorado.edu.
************************************************************************
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
]]></plaintext></body>