Here is the Fortran90 Program for Molecular Dynamics Excersise
Subroutines for Potential Energy and Force are included
<![CDATA[
PROGRAM MOLEC_D
!
! 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, NSTEP, PRNFREQ
INTEGER, 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.
!*******************************************************************
REAL :: H, HO2, HP, PE, KE, TEND
REAL, DIMENSION(NATOMX) :: X, Y, Z, U, V, W
REAL, DIMENSION(NATOMX,2) :: FX, FY, FZ
OPEN (UNIT = 2, FILE = "md90.in")
OPEN (UNIT = 3, FILE = "md90.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 --
!Read computation parameters
10 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,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,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)
99 FORMAT(F13.6, TR1, I5, TR1, 3(TR1,E13.6,TR1))
END DO
END IF
ITMP = B1
B1 = B2
B2 = ITMP
END DO
GO TO 10
999 STOP
CONTAINS
!*****************************************************************
!potential_lj_ymp.f
! Compute potential energy from Lennard Jones interaction.
! V(r) = (1/r**12 - 2/r**6)
!*****************************************************************
FUNCTION POTENTIAL_LJ(NATOM,X,Y,Z)
IMPLICIT NONE
INTEGER :: I,J
INTEGER, INTENT(IN) :: NATOM
!******************************************************************
! I,J - Loop indices.
! NATOM - Number of atoms in chain.
!******************************************************************
REAL :: POT, RIJSQ, RIJM6, POTENTIAL_LJ
REAL, DIMENSION(:), INTENT(IN) :: X, Y, Z
!******************************************************************
! 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 FUNCTION POTENTIAL_LJ
!********************************************************************
! force_lj.f
! Compute forces from Lennard Jones interaction.
! V(r) = (1/r**12 - 2/r**6)
!********************************************************************
SUBROUTINE FORCE_LJ(NATOM,X,Y,Z,B,FX,FY,FZ)
IMPLICIT NONE
INTEGER, PARAMETER :: NATOMXX = 50
INTEGER, INTENT(IN) :: NATOM, B
INTEGER :: I, J
!********************************************************************
! I,J - Loop indices.
! NATOM - Number of atoms in chain.
! NATOMX - Maximum number of atoms permitted in chain.
! B - Buffer index.
!********************************************************************
REAL, DIMENSION(:), INTENT(IN) :: X, Y, Z
REAL, DIMENSION(:,:), INTENT(IN OUT) :: FX, FY, FZ
REAL :: RIJSQ, RIJM6, CIJ
REAL, DIMENSION(NATOMXX) :: FIJX, FIJY, FIJZ
!********************************************************************
! 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 SUBROUTINE FORCE_LJ
END PROGRAM MOLEC_D
]]></plaintext></body>