*JOB TIME=10M,MEMORY=24K,PAGES=30
*
*	This is a simple example batch job, which runs the ChemE
*	Molecular Dynamics benchmark program (cut-down form).
*
*	Steve Thompson, May 1983.
*
*NOECHO
$FORTRAN/F77/OBJECT=EXAMPLE/NOTRACE  TI:
	PROGRAM MDATOM
C**************************************************************
C
C	MOLECULAR DYNAMICS PROGRAM FOR LENNARD-JONES ATOMS
C	THE EQUATIONS OF MOTION ARE SOLVED USING A PREDICTOR
C	CORRECTOR ALGORITHM. BENCHMARK VERSION (STRIPPED DOWN).
C
C	STEVE THOMPSON
C	SCHOOL OF CHEMICAL ENGINEERING
C	OLIN HALL
C	CORNELL UNIVERSITY
C	ITHACA
C	NY 14853
C
C****************************************************************
C
	COMMON /CMCOM0/ X0(108),Y0(108),Z0(108)
	COMMON /CMCOM1/ X1(108),Y1(108),Z1(108)
	COMMON /CMCOM2/ X2(108),Y2(108),Z2(108)
	COMMON /CMCOM3/ X3(108),Y3(108),Z3(108)
	COMMON /CMFRCE/ FRCEX(108),FRCEY(108),FRCEZ(108)
	COMMON /CMMOVE/ XMOVE(108),YMOVE(108),ZMOVE(108)
	DOUBLE PRECISION SUMTEN,SSQTEN
	COMMON /CMAVER/ SUMTMP,SSQTMP,SUMENR,SSQENR,SUMVIR,
	1		SSQVIR,SUMTEN,SSQTEN,SUMFC2,SSQFC2,
	2		SUMVT
	COMMON /CMBOX/ ULEN
	COMMON /CMCCTR/ CT0,CT1,CT2,CT3
	REAL NAOVRN
	COMMON /CMCNST/ PI,EFACT,AVOGAD,BOLTZ,GASCON,NAOVRN
	COMMON /CMCONF/ UCONF,VIRAA,ECONF,EVIRAA
	COMMON /CMDELT/ TSTEP,TSTEP2,TIMSEC
	COMMON /CMHEAT/ ENTRA
	COMMON /CMLINE/ LINPAG,LINCNT,LINRUN,LINEQU,
	1		LINSCL
	REAL MASS
	COMMON /CMPARM/ MASS,EPSILN,SIGMA,SIGPAR,RCUTS,RCUTU
	COMMON /CMPROP/ TOTE,TEMP
	COMMON /CMSCAL/ NSCALE,DTEMP,VSCALE
	COMMON /CMSTEP/ NSTEP,FNSTEP,NRJECT,MAXSTP,
	1		NPRINT,NPROP
	COMMON /CMSTPT/ RHO,TR,VOLM,T
	COMMON /CMSUMV/ SUMX1,SUMY1,SUMZ1,SUMV2
	COMMON /CMTAIL/ ULTAIL,VLTAIL
C
C ===== PERFORM STARTUP PROCESSING
C
	CALL START
C
C ===== SELECT INITIAL CONFIGURATION FOR THIS RUN
C
	CALL FCC
C
C ===== INITIALISE SOME DERIVATIVES IF THIS IS A NEW RUN
C
	CALL FORCEN
	DO 5 I = 1,108
	X2(I) = FRCEX(I)*TSTEP2
	Y2(I) = FRCEY(I)*TSTEP2
	Z2(I) = FRCEZ(I)*TSTEP2
5	CONTINUE
C
C ===== PRINT OUT COLUMN HEADERS
C
10	CALL HEADER
C
C ===== PREDICT NEW COORDINATES AND DERIVATIVES
C
15	CALL PREDIC
C
C ===== EVALUATE FORCES
C
	CALL FORCEN
C
C ===== APPLY CORRECTIONS TO PREDICTED QUANTITIES
C
	CALL CORR
C
C ===== PRINT CURRENT STEP RESULTS IF REQUIRED
C
	IF(NPRINT.EQ.0)GO TO 20
	IF(MOD(NSTEP,NPRINT).EQ.0)CALL PRINT
C
C ===== TEST FOR PHASE OF RUN. PHASES ARE:
C	(1) NSTEP < NRJECT -- EQUILIBRATION PHASE
C	(2) NSTEP = NRJECT -- END OF EQUILIBRATION, CLEANUP REQUIRED
C	(3) NSTEP > NRJECT (NRJECT=0) -- PRODUCTION RUN
C
20	IF(NSTEP - NRJECT)35,30,25
C
C ===== PRODUCTION PHASE; PRINT RUNNING AVERAGES IF REQUIRED
C
25	IF(NPROP.EQ.0)GO TO 35
	IF(MOD(NSTEP,NPROP).NE.0)GO TO 35
	CALL RUNAVS
	CALL HEADER
	GO TO 35
C
C ===== END OF EQUILIBRATION; RESET AVERAGES IF PRODUCTION PHASE
C	IS TO FOLLOW
C
30	IF(MAXSTP.NE.0)CALL ZERO
C
C ===== END OF CURRENT STEP; ALL PHASES CONTINUE HERE
C
35	CONTINUE
C
C ===== CHECK TO SEE IF ALL REQUESTED STEPS HAVE BEEN DONE
C
	IF(NSTEP.LT.NRJECT)GO TO 15
	IF(NSTEP.LT.MAXSTP)GO TO 15
C
C ===== END OF RUN PROCESSING
C
50	CONTINUE
	STOP
	END
	BLOCK DATA
C
C
C ===== CERTAIN DERIVATIVES AND RUNNING AVERAGE ACCUMULATORS MUST
C	HAVE ZERO VALUES AT THE START OF A  RUN
C
	COMMON /CMCOM0/ X0(108),Y0(108),Z0(108)
	COMMON /CMCOM1/ X1(108),Y1(108),Z1(108)
	COMMON /CMCOM2/ X2(108),Y2(108),Z2(108)
	COMMON /CMCOM3/ X3(108),Y3(108),Z3(108)
	COMMON /CMMOVE/ XMOVE(108),YMOVE(108),ZMOVE(108)
	DOUBLE PRECISION SUMTEN,SSQTEN
	COMMON /CMAVER/ SUMTMP,SSQTMP,SUMENR,SSQENR,SUMVIR,
	1		SSQVIR,SUMTEN,SSQTEN,SUMFC2,SSQFC2,
	2		SUMVT
	COMMON /CMBOX/ ULEN
	COMMON /CMCCTR/ CT0,CT1,CT2,CT3
	REAL NAOVRN
	COMMON /CMCNST/ PI,EFACT,AVOGAD,BOLTZ,GASCON,NAOVRN
	COMMON /CMCONF/ UCONF,VIRAA,ECONF,EVIRAA
	COMMON /CMDELT/ TSTEP,TSTEP2,TIMSEC
	COMMON /CMHEAT/ ENTRA
	COMMON /CMLINE/ LINPAG,LINCNT,LINRUN,LINEQU,
	1		LINSCL
	REAL MASS
	COMMON /CMPARM/ MASS,EPSILN,SIGMA,SIGPAR,RCUTS,RCUTU
	COMMON /CMPROP/ TOTE,TEMP
	COMMON /CMSCAL/ NSCALE,DTEMP,VSCALE
	COMMON /CMSTEP/ NSTEP,FNSTEP,NRJECT,MAXSTP,
	1		NPRINT,NPROP
	COMMON /CMSTPT/ RHO,TR,VOLM,T
	COMMON /CMSUMV/ SUMX1,SUMY1,SUMZ1,SUMV2
	COMMON /CMTAIL/ ULTAIL,VLTAIL
C
	DATA X3/108*0.0/, Y3/108*0.0/, Z3/108*0.0/
	DATA XMOVE/108*0.0/, YMOVE/108*0.0/, ZMOVE/108*0.0/
	DATA PI /3.1415927/
	DATA AVOGAD,BOLTZ,GASCON /6.0225E23,1.38054E-23,8.3143/
C
	DATA SUMTMP,SSQTMP,SUMENR,SSQENR,SUMVIR,SSQVIR,SUMTEN /7*0.0/
	DATA SSQTEN,SUMVT /2*0.0/
	DATA CT0,CT1,CT2,CT3 /
	1  0.16666666666667,0.83333333333333,1.0,
	2  0.33333333333333 /
	DATA LINPAG,LINRUN,LINEQU,LINSCL /55,2,3,5/
	DATA NSCALE,DTEMP /0,2.5/
	END
	SUBROUTINE HEADER
	COMMON /CMAVER/ SUMTMP,SSQTMP,SUMENR,SSQENR,SUMVIR,
	1		SSQVIR,SUMTEN,SSQTEN,SUMFC2,SSQFC2,
	2		SUMVT
	COMMON /CMBOX/ ULEN
	COMMON /CMCCTR/ CT0,CT1,CT2,CT3
	REAL NAOVRN
	COMMON /CMCNST/ PI,EFACT,AVOGAD,BOLTZ,GASCON,NAOVRN
	COMMON /CMCONF/ UCONF,VIRAA,ECONF,EVIRAA
	COMMON /CMDELT/ TSTEP,TSTEP2,TIMSEC
	COMMON /CMHEAT/ ENTRA
	COMMON /CMLINE/ LINPAG,LINCNT,LINRUN,LINEQU,
	1		LINSCL
	REAL MASS
	COMMON /CMPARM/ MASS,EPSILN,SIGMA,SIGPAR,RCUTS,RCUTU
	COMMON /CMPROP/ TOTE,TEMP
	COMMON /CMSCAL/ NSCALE,DTEMP,VSCALE
	COMMON /CMSTEP/ NSTEP,FNSTEP,NRJECT,MAXSTP,
	1		NPRINT,NPROP
	COMMON /CMSTPT/ RHO,TR,VOLM,T
	COMMON /CMSUMV/ SUMX1,SUMY1,SUMZ1,SUMV2
	COMMON /CMTAIL/ ULTAIL,VLTAIL
C
C ===== DETERMINE IF VELOCITY SCALING DATA PRINTOUT IS REQUIRED
C
	IF(NPRINT.EQ.0)RETURN
	LINCNT = LINPAG
	IF(NRJECT.NE.0)LINCNT = LINCNT - LINSCL
	IF((NRJECT.EQ.0).OR.(NSTEP.EQ.0))GO TO 5
	TEMPX = SUMTMP/FLOAT(NSTEP)
	WRITE(5,2) NSCALE,NSTEP,TEMPX
2	FORMAT(/,T5,'*** Velocities have been scaled',I6,
	1  ' times in',I6,' steps ***',//,
	2  T5,'*** Average Temperature =',F8.2,' K ***')
5	CONTINUE
C
C ===== PRINT COLUMN HEADER ON A NEW PAGE
C
	WRITE(5,10)
10	FORMAT(1H1,//,3X,'Steps',4X,'E(Total)',6X,'Uconf',8X,'Virial',
	1  6X,'Temp',4X,'Mean Square',/,
	2  6X,3(7X,'Kj/Mol'),6X,'(K)',7X,'Displ')
	RETURN
	END
	SUBROUTINE CORR
	COMMON /CMCOM0/ X0(108),Y0(108),Z0(108)
	COMMON /CMCOM1/ X1(108),Y1(108),Z1(108)
	COMMON /CMCOM2/ X2(108),Y2(108),Z2(108)
	COMMON /CMCOM3/ X3(108),Y3(108),Z3(108)
	COMMON /CMMOVE/ XMOVE(108),YMOVE(108),ZMOVE(108)
	COMMON /CMFRCE/ FRCEX(108),FRCEY(108),FRCEZ(108)
	DOUBLE PRECISION SUMTEN,SSQTEN
	COMMON /CMAVER/ SUMTMP,SSQTMP,SUMENR,SSQENR,SUMVIR,
	1		SSQVIR,SUMTEN,SSQTEN,SUMFC2,SSQFC2,
	2		SUMVT
	COMMON /CMBOX/ ULEN
	COMMON /CMCCTR/ CT0,CT1,CT2,CT3
	REAL NAOVRN
	COMMON /CMCNST/ PI,EFACT,AVOGAD,BOLTZ,GASCON,NAOVRN
	COMMON /CMCONF/ UCONF,VIRAA,ECONF,EVIRAA
	COMMON /CMDELT/ TSTEP,TSTEP2,TIMSEC
	COMMON /CMHEAT/ ENTRA
	COMMON /CMLINE/ LINPAG,LINCNT,LINRUN,LINEQU,
	1		LINSCL
	REAL MASS
	COMMON /CMPARM/ MASS,EPSILN,SIGMA,SIGPAR,RCUTS,RCUTU
	COMMON /CMPROP/ TOTE,TEMP
	COMMON /CMSCAL/ NSCALE,DTEMP,VSCALE
	COMMON /CMSTEP/ NSTEP,FNSTEP,NRJECT,MAXSTP,
	1		NPRINT,NPROP
	COMMON /CMSTPT/ RHO,TR,VOLM,T
	COMMON /CMSUMV/ SUMX1,SUMY1,SUMZ1,SUMV2
	COMMON /CMTAIL/ ULTAIL,VLTAIL
C
	REAL KETRAN
C
C ===== ZERO COUNTERS
C
	SUMV2 = 0.0
	SUMF2 = 0.0
	SUMX1 = 0.0
	SUMY1 = 0.0
	SUMZ1 = 0.0
C
C ===== CORRECTOR LOOP FOR CENTRE OF MASS MOTION
C
	DO 100 I = 1,108
	XCORR = X2(I) - FRCEX(I)*TSTEP2
	YCORR = Y2(I) - FRCEY(I)*TSTEP2
	ZCORR = Z2(I) - FRCEZ(I)*TSTEP2
	X0(I) = X0(I) - XCORR*CT0
	Y0(I) = Y0(I) - YCORR*CT0
	Z0(I) = Z0(I) - ZCORR*CT0
	X1(I) = X1(I) - XCORR*CT1
	Y1(I) = Y1(I) - YCORR*CT1
	Z1(I) = Z1(I) - ZCORR*CT1
	X2(I) = X2(I) - XCORR
	Y2(I) = Y2(I) - YCORR
	Z2(I) = Z2(I) - ZCORR
	X3(I) = X3(I) - XCORR*CT3
	Y3(I) = Y3(I) - YCORR*CT3
	Z3(I) = Z3(I) - ZCORR*CT3
C
C ===== UPDATE CENTRE OF MASS DISPLACEMENT VECTORS
C
	XMOVE(I) = XMOVE(I) + X0(I)
	YMOVE(I) = YMOVE(I) + Y0(I)
	ZMOVE(I) = ZMOVE(I) + Z0(I)
C
C ===== INCREMENT TRANSLATIONAL KINETIC ENERGY
C
	SUMV2 = SUMV2 + X1(I)**2 + Y1(I)**2 + Z1(I)**2
C
C ===== ACCUMULATED C.O.M VELOCITY OF COMPLETE SYSTEM
C
	SUMX1 = SUMX1 + X1(I)
	SUMY1 = SUMY1 + Y1(I)
	SUMZ1 = SUMZ1 + Z1(I)
C
C ===== APPLY PERIODIC BOUNDARY CONDITIONS
C
	X0(I) = X0(I) - 2*INT(X0(I))
	Y0(I) = Y0(I) - 2*INT(Y0(I))
	Z0(I) = Z0(I) - 2*INT(Z0(I))
C
C ===== INCREMENT MEAN SQUARE FORCE ACCUMULATOR
C
	SUMF2 = SUMF2 + FRCEX(I)**2 + FRCEY(I)**2 + FRCEZ(I)**2
100	CONTINUE
C
C ===== CALCULATE MOLAR PROPERTIES AND TEMPERATURE
C
	SUMX1 = SUMX1/108.0
	SUMY1 = SUMY1/108.0
	SUMZ1 = SUMZ1/108.0
	KETRAN = 0.25*EFACT*SUMV2/TSTEP2
	TEMP = 2.0*KETRAN/3.0/GASCON
	ECONF = EFACT*UCONF + ULTAIL
	EVIRAA = EFACT*VIRAA + VLTAIL
	TOTE = KETRAN + ECONF
C
C ===== UPDATE RUNNING AVERAGES
C
	SUMTMP = SUMTMP + TEMP
	SSQTMP = SSQTMP + TEMP**2
	SUMFC2 = SUMFC2 + SUMF2
	SSQFC2 = SSQFC2 + SUMF2**2
	SUMTEN = SUMTEN + TOTE
	SSQTEN = SSQTEN + TOTE**2
	SUMENR = SUMENR + ECONF
	SSQENR = SSQENR + ECONF*ECONF
	SUMVIR = SUMVIR + EVIRAA
	SSQVIR = SSQVIR + EVIRAA**2
	SUMVT = SUMVT + TEMP*EVIRAA
	IF(NRJECT.EQ.0)RETURN
C
C ===== SCALE VELOCITIES TO MINTAIN TEMPERATURE
C
C
	IF(ABS(T - TEMP).LT.DTEMP)RETURN
	NSCALE = NSCALE + 1
	SUMV2 = SUMV2 - 108.0*(SUMX1**2 + SUMY1**2 + SUMZ1**2)
	VSCALE = SQRT(ENTRA/SUMV2)
	DO 400 I = 1,108
	X1(I) = (X1(I) - SUMX1)*VSCALE
	Y1(I) = (Y1(I) - SUMY1)*VSCALE
	Z1(I) = (Z1(I) - SUMZ1)*VSCALE
400	CONTINUE
	RETURN
	END
	SUBROUTINE PREDIC
	COMMON /CMCOM0/ X0(108),Y0(108),Z0(108)
	COMMON /CMCOM1/ X1(108),Y1(108),Z1(108)
	COMMON /CMCOM2/ X2(108),Y2(108),Z2(108)
	COMMON /CMCOM3/ X3(108),Y3(108),Z3(108)
	COMMON /CMMOVE/ XMOVE(108),YMOVE(108),ZMOVE(108)
	DOUBLE PRECISION SUMTEN,SSQTEN
	COMMON /CMAVER/ SUMTMP,SSQTMP,SUMENR,SSQENR,SUMVIR,
	1		SSQVIR,SUMTEN,SSQTEN,SUMFC2,SSQFC2,
	2		SUMVT
	COMMON /CMBOX/ ULEN
	COMMON /CMCCTR/ CT0,CT1,CT2,CT3
	REAL NAOVRN
	COMMON /CMCNST/ PI,EFACT,AVOGAD,BOLTZ,GASCON,NAOVRN
	COMMON /CMCONF/ UCONF,VIRAA,ECONF,EVIRAA
	COMMON /CMDELT/ TSTEP,TSTEP2,TIMSEC
	COMMON /CMHEAT/ ENTRA
	COMMON /CMLINE/ LINPAG,LINCNT,LINRUN,LINEQU,
	1		LINSCL
	REAL MASS
	COMMON /CMPARM/ MASS,EPSILN,SIGMA,SIGPAR,RCUTS,RCUTU
	COMMON /CMPROP/ TOTE,TEMP
	COMMON /CMSCAL/ NSCALE,DTEMP,VSCALE
	COMMON /CMSTEP/ NSTEP,FNSTEP,NRJECT,MAXSTP,
	1		NPRINT,NPROP
	COMMON /CMSTPT/ RHO,TR,VOLM,T
	COMMON /CMSUMV/ SUMX1,SUMY1,SUMZ1,SUMV2
	COMMON /CMTAIL/ ULTAIL,VLTAIL
C
C ===== UPDATE TIME STEP COUNTERS
C
	NSTEP = NSTEP + 1
	FNSTEP = FLOAT(NSTEP)
C
C ===== INITIALISE C.O.M. DISPLACEMENT VECTORS
C
	DO 100 I = 1,108
	XMOVE(I) = XMOVE(I) - X0(I)
	YMOVE(I) = YMOVE(I) - Y0(I)
	ZMOVE(I) = ZMOVE(I) - Z0(I)
C
C ===== PREDICTOR LOOP FOR CENTRE OF MASS MOTION
C
	X0(I) = X0(I) + X1(I) + X2(I) + X3(I)
	Y0(I) = Y0(I) + Y1(I) + Y2(I) + Y3(I)
	Z0(I) = Z0(I) + Z1(I) + Z2(I) + Z3(I)
	X1(I) = X1(I) + 2.0*X2(I) + 3.0*X3(I)
	Y1(I) = Y1(I) + 2.0*Y2(I) + 3.0*Y3(I)
	Z1(I) = Z1(I) + 2.0*Z2(I) + 3.0*Z3(I)
	X2(I) = X2(I) + 3.0*X3(I)
	Y2(I) = Y2(I) + 3.0*Y3(I)
	Z2(I) = Z2(I) + 3.0*Z3(I)
100	CONTINUE
	RETURN
	END
	SUBROUTINE FCC
	COMMON /CMCOM0/ X0(108),Y0(108),Z0(108)
	COMMON /CMCOM1/ X1(108),Y1(108),Z1(108)
	COMMON /CMCOM2/ X2(108),Y2(108),Z2(108)
	COMMON /CMCOM3/ X3(108),Y3(108),Z3(108)
	DOUBLE PRECISION SUMTEN,SSQTEN
	COMMON /CMAVER/ SUMTMP,SSQTMP,SUMENR,SSQENR,SUMVIR,
	1		SSQVIR,SUMTEN,SSQTEN,SUMFC2,SSQFC2,
	2		SUMVT
	COMMON /CMBOX/ ULEN
	COMMON /CMCCTR/ CT0,CT1,CT2,CT3
	REAL NAOVRN
	COMMON /CMCNST/ PI,EFACT,AVOGAD,BOLTZ,GASCON,NAOVRN
	COMMON /CMCONF/ UCONF,VIRAA,ECONF,EVIRAA
	COMMON /CMDELT/ TSTEP,TSTEP2,TIMSEC
	COMMON /CMHEAT/ ENTRA
	COMMON /CMLINE/ LINPAG,LINCNT,LINRUN,LINEQU,
	1		LINSCL
	REAL MASS
	COMMON /CMPARM/ MASS,EPSILN,SIGMA,SIGPAR,RCUTS,RCUTU
	COMMON /CMPROP/ TOTE,TEMP
	COMMON /CMSCAL/ NSCALE,DTEMP,VSCALE
	COMMON /CMSTEP/ NSTEP,FNSTEP,NRJECT,MAXSTP,
	1		NPRINT,NPROP
	COMMON /CMSTPT/ RHO,TR,VOLM,T
	COMMON /CMSUMV/ SUMX1,SUMY1,SUMZ1,SUMV2
	COMMON /CMTAIL/ ULTAIL,VLTAIL
C
	INTEGER*4 JSEED
C
	DATA DSPLCE /0.05/
	DATA JSEED /123456789/
C
C ===== SET UP F.C.C. LATTICE AND ASSIGN INITIAL VELOCITIES
C
	NSTEP = 0
	FNSTEP = 0.0
	VT = 2.0*SQRT(ENTRA/108.0/3.0)
	NC = (108.0/4.0)**(1.0/3.0) + 0.1
	CUBE = 2.0/FLOAT(NC)
	DSPLCE = CUBE*DSPLCE
	X0(1) =  - 1.0
	X0(2) = 0.5*CUBE - 1.0
	X0(3) =  - 1.0
	X0(4) = 0.5*CUBE - 1.0
	Y0(1) =  - 1.0
	Y0(2) = 0.5*CUBE - 1.0
	Y0(3) = 0.5*CUBE - 1.0
	Y0(4) =  - 1.0
	Z0(1) =  - 1.0
	Z0(2) =  - 1.0
	Z0(3) = 0.5*CUBE - 1.0
	Z0(4) = 0.5*CUBE - 1.0
	M = 0
	DO 20 I = 1,NC
	DO 20 J = 1,NC
	DO 20 K = 1,NC
	DO 10 IJ = 1,4
	IJM = IJ + M
	IF(IJM.GT.108)GO TO 10
	X0(IJM) = X0(IJ) + CUBE*FLOAT(K - 1)
	Y0(IJM) = Y0(IJ) + CUBE*FLOAT(J - 1)
	Z0(IJM) = Z0(IJ) + CUBE*FLOAT(I - 1)
10	CONTINUE
	M = M + 4
20	CONTINUE
C
C ===== DISPLACE LATTICE A LITTLE TO HELP EQUILIBRATION
C
	SUMX1 = 0.0
	SUMY1 = 0.0
	SUMZ1 = 0.0
C
	XCOM = 0.0
	YCOM = 0.0
	ZCOM = 0.0
C
	SUMV2 = 0.0
C
	DO 40 I = 1,108
	X0(I) = X0(I) + DSPLCE*(2.0*RAN(JSEED) - 1.0)
	Y0(I) = Y0(I) + DSPLCE*(2.0*RAN(JSEED) - 1.0)
	Z0(I) = Z0(I) + DSPLCE*(2.0*RAN(JSEED) - 1.0)
	X1(I) = VT*(2.0*RAN(JSEED) - 1.0)
	Y1(I) = VT*(2.0*RAN(JSEED) - 1.0)
	Z1(I) = VT*(2.0*RAN(JSEED) - 1.0)
	SUMX1 = SUMX1 + X1(I)
	SUMY1 = SUMY1 + Y1(I)
	SUMZ1 = SUMZ1 + Z1(I)
40	CONTINUE
C
	XCOM = XCOM/108.0
	YCOM = YCOM/108.0
	ZCOM = ZCOM/108.0
C
	SUMX1 = SUMX1/108.0
	SUMY1 = SUMY1/108.0
	SUMZ1 = SUMZ1/108.0
C
	DO 50 I = 1,108
C
C ===== ENSURE ALL CENTRES OF MASS ARE WITHIN CUBE, AND THAT THE
C	CENTRE OF MASS OF THE WHOLE SYSTEM IS IN THE RIGHT PLACE.
C
	X0(I) = X0(I) - XCOM
	Y0(I) = Y0(I) - YCOM
	Z0(I) = Z0(I) - ZCOM
	X0(I) = X0(I) - 2*INT(X0(I))
	Y0(I) = Y0(I) - 2*INT(Y0(I))
	Z0(I) = Z0(I) - 2*INT(Z0(I))
C
C ===== ADJUST VELOCITIES SO THAT C.O.M HAS NO RESULTANT VELOCITY
C
	X1(I) = X1(I) - SUMX1
	Y1(I) = Y1(I) - SUMY1
	Z1(I) = Z1(I) - SUMZ1
	SUMV2 = SUMV2 + X1(I)**2 + Y1(I)**2 + Z1(I)**2
50	CONTINUE
C
C ===== SCALE VELOCITIES TO DESIRED TEMPERATURE.
C
	VSCALE = SQRT(ENTRA/SUMV2)
	DO 60 I = 1,108
	X1(I) = X1(I)*VSCALE
	Y1(I) = Y1(I)*VSCALE
	Z1(I) = Z1(I)*VSCALE
60	CONTINUE
	RETURN
	END
	SUBROUTINE START
	COMMON /CMAVER/ SUMTMP,SSQTMP,SUMENR,SSQENR,SUMVIR,
	1		SSQVIR,SUMTEN,SSQTEN,SUMFC2,SSQFC2,
	2		SUMVT
	COMMON /CMBOX/ ULEN
	COMMON /CMCCTR/ CT0,CT1,CT2,CT3
	REAL NAOVRN
	COMMON /CMCNST/ PI,EFACT,AVOGAD,BOLTZ,GASCON,NAOVRN
	COMMON /CMCONF/ UCONF,VIRAA,ECONF,EVIRAA
	COMMON /CMDELT/ TSTEP,TSTEP2,TIMSEC
	COMMON /CMHEAT/ ENTRA
	COMMON /CMLINE/ LINPAG,LINCNT,LINRUN,LINEQU,
	1		LINSCL
	REAL MASS
	COMMON /CMPARM/ MASS,EPSILN,SIGMA,SIGPAR,RCUTS,RCUTU
	COMMON /CMPROP/ TOTE,TEMP
	COMMON /CMSCAL/ NSCALE,DTEMP,VSCALE
	COMMON /CMSTEP/ NSTEP,FNSTEP,NRJECT,MAXSTP,
	1		NPRINT,NPROP
	COMMON /CMSTPT/ RHO,TR,VOLM,T
	COMMON /CMSUMV/ SUMX1,SUMY1,SUMZ1,SUMV2
	COMMON /CMTAIL/ ULTAIL,VLTAIL
C
	NAOVRN = AVOGAD/108.0
C
C ===== SET RUN PARAMETERS DATA
C
	VOLM = 35.0E-6
	T = 85.0
	MAXSTP = 100
	NRJECT = 100
	EPSILN = 119.8
	SIGMA = 0.3405
	MASS = 39.95
	RCUTS = 2.5
	TIMSEC = 1.0E-14
	NPROP = 500
	NPRINT = 10
C
C ===== VALIDATE CERTAIN INPUT DATA ITEMS
C
	IF(EPSILN.GT.1.0E-5)EPSILN = EPSILN*BOLTZ
	IF(SIGMA.GT.1.0E-3)SIGMA = SIGMA*1.0E-9
C
C ===== PRINT RUN IDENTIFICATION
C
	WRITE(5,503)
503	FORMAT(///,T5,'MOLECULAR DYNAMICS SIMULATION FOR 108',
	1 ' LENNARD-JONES ATOMS')
C
C ===== EVALUATE VARIOUS BOX AND MODEL PARAMETERS
C
	XL = (VOLM/NAOVRN)**(1.0/3.0)
	RHO = AVOGAD*SIGMA**3/VOLM
	ULEN = XL/2.0
	XLS = XL/SIGMA
	TR = BOLTZ*T/EPSILN
	IF(RCUTS.GT.XLS/2.0)RCUTS = XLS/2.0
	RCUTX = RCUTS*SIGMA
	RCUTU = (RCUTS*SIGMA/ULEN)**2
	SIGPAR = (SIGMA/ULEN)**2
	XMASS = MASS/AVOGAD/1000.0
	ULENX = 1.0E10*ULEN
	STEP = SQRT(XMASS*ULENX**2/EPSILN)*1.0E-10
	TSTEP = TIMSEC/STEP
	TSTEP2 = TSTEP**2/2.0
	ENTRA = 6.0*TSTEP2*108.0*TR
	EPSK = EPSILN/BOLTZ
C
C ===== PRINT SYSTEM PROPERTIES
C
	WRITE(5,20) EPSK,EPSILN,SIGMA,MASS,XMASS
	WRITE(5,30) TR,T,RHO,VOLM
	WRITE(5,40) XLS,XL,RCUTS,RCUTX
	WRITE(5,50) STEP,TSTEP,TIMSEC
20	FORMAT(/,T5,'Lennard-Jones Well Depth',T40,F11.3,
	1  ' Kelvin',T65,1PE11.4,' Joules',//,
	2  T5,'Lennard-Jones Collision Diameter',T40,9PF11.4,
	3  ' nm',//,
	6  T5,'Atomic Mass',T40,0PF11.3,' a.m.u.',
	7  T65,1PE11.4,' kg')
30	FORMAT(/,T5,'Desired Temperature',T40,F11.4,
	1  ' (Reduced)',T65,F11.2,' Kelvin',//,
	2  T5,'Reduced Density',T40,F11.4,//,
	3  T5,'Molar Volume',T40,2PE11.4,' m**3')
40	FORMAT(/,T5,'Box Width',T40,F11.4,' Sigma',
	1  T65,9PF11.4,' nm',//,
	1  T5,'Potential Truncation Distance',T40,0PF11.4,
	2  ' Sigma',T65,9PF11.4,' nm')
50	FORMAT(/,T5,'Uunit of Time',T40,1PE11.4,' Seconds',//,
	1  T5,'Reduced Time Step',T40,0PF11.5,//,
	2  T5,'Real Time Step',T40,1PE11.3,' Seconds')
C
C ===== CALCULATE CONVERSION FACTORS FOR CALCULATED PROPERTIES
C
	EFACT = EPSILN*NAOVRN
C
C ===== CALCULATE POTENTIAL TAIL TRUNCATION CORRECTIONS
C
	CLTAIL = 8.0*PI*RHO*AVOGAD*EPSILN
	RCUT3 = 1.0/RCUTS**3
	RCUT9 = RCUT3**3
	ULTAIL = CLTAIL*(RCUT9/9.0 - RCUT3/3.0)
	VLTAIL = CLTAIL*( - 4.0*RCUT9/3.0 + 2.0*RCUT3)
	WRITE(5,60) ULTAIL,VLTAIL
60	FORMAT(/,T5,'Tail Corrections (Kj/Mol)',
	1  T40,-3PF11.4,' (Energy)',
	2  T65,-3PF11.4,' (Virial)')
	RETURN
	END
	SUBROUTINE PRINT
	COMMON /CMCOM0/ X0(108),Y0(108),Z0(108)
	COMMON /CMCOM1/ X1(108),Y1(108),Z1(108)
	COMMON /CMCOM2/ X2(108),Y2(108),Z2(108)
	COMMON /CMCOM3/ X3(108),Y3(108),Z3(108)
	COMMON /CMMOVE/ XMOVE(108),YMOVE(108),ZMOVE(108)
	DOUBLE PRECISION SUMTEN,SSQTEN
	COMMON /CMAVER/ SUMTMP,SSQTMP,SUMENR,SSQENR,SUMVIR,
	1		SSQVIR,SUMTEN,SSQTEN,SUMFC2,SSQFC2,
	2		SUMVT
	COMMON /CMBOX/ ULEN
	COMMON /CMCCTR/ CT0,CT1,CT2,CT3
	REAL NAOVRN
	COMMON /CMCNST/ PI,EFACT,AVOGAD,BOLTZ,GASCON,NAOVRN
	COMMON /CMCONF/ UCONF,VIRAA,ECONF,EVIRAA
	COMMON /CMDELT/ TSTEP,TSTEP2,TIMSEC
	COMMON /CMHEAT/ ENTRA
	COMMON /CMLINE/ LINPAG,LINCNT,LINRUN,LINEQU,
	1		LINSCL
	REAL MASS
	COMMON /CMPARM/ MASS,EPSILN,SIGMA,SIGPAR,RCUTS,RCUTU
	COMMON /CMPROP/ TOTE,TEMP
	COMMON /CMSCAL/ NSCALE,DTEMP,VSCALE
	COMMON /CMSTEP/ NSTEP,FNSTEP,NRJECT,MAXSTP,
	1		NPRINT,NPROP
	COMMON /CMSTPT/ RHO,TR,VOLM,T
	COMMON /CMSUMV/ SUMX1,SUMY1,SUMZ1,SUMV2
	COMMON /CMTAIL/ ULTAIL,VLTAIL
C
C ===== DETERMINE WHETHER TO START A NEW PAGE OR NOT
C
	IF(NRJECT.EQ.0)LINCNT = LINCNT - LINRUN
	IF(NRJECT.NE.0)LINCNT = LINCNT - LINEQU
	IF(LINCNT.LE.0)CALL HEADER
C
C ===== CALCULATE MEAN SQUARE DISPLACEMENT
C
	RMOVE = 0.0
	DO 100 I = 1,108
	RMOVE = RMOVE + XMOVE(I)**2 + YMOVE(I)**2 + ZMOVE(I)**2
100	CONTINUE
	RMOVE = RMOVE/108.0*ULEN**2
	WRITE(5,10) NSTEP,TOTE,ECONF,EVIRAA,TEMP,RMOVE
10	FORMAT(/,I7,3(-3PF13.4),0PF10.3,1PE14.4)
C
C ===== CALCULATE MELTING FACTOR
C
	IF(NRJECT.EQ.0)RETURN
	CONST = PI*(2.0*108.0)**(1.0/3.0)
C
	CX = 0.0
	CY = 0.0
	CZ = 0.0
	SX = 0.0
	SY = 0.0
	SZ = 0.0
	DO 200 I = 1,108
	CX = CX + COS(CONST*X0(I))
	SX = SX + SIN(CONST*X0(I))
	CY = CY + COS(CONST*Y0(I))
	SY = SY + COS(CONST*Y0(I))
	CZ = SZ + COS(CONST*Z0(I))
	SZ = SZ + SIN(CONST*Z0(I))
200	CONTINUE
	FMELT = (CX**2 + CY**2 + CZ**2 + SX**2 + SY**2 + SZ**2)/
	1 108.0/3.0
	WRITE(5,20) FMELT
20	FORMAT(T15,' Melting Factor =',F10.4)
	RETURN
	END
	SUBROUTINE RUNAVS
	DOUBLE PRECISION SUMTEN,SSQTEN
	COMMON /CMAVER/ SUMTMP,SSQTMP,SUMENR,SSQENR,SUMVIR,
	1		SSQVIR,SUMTEN,SSQTEN,SUMFC2,SSQFC2,
	2		SUMVT
	COMMON /CMBOX/ ULEN
	COMMON /CMCCTR/ CT0,CT1,CT2,CT3
	REAL NAOVRN
	COMMON /CMCNST/ PI,EFACT,AVOGAD,BOLTZ,GASCON,NAOVRN
	COMMON /CMCONF/ UCONF,VIRAA,ECONF,EVIRAA
	COMMON /CMDELT/ TSTEP,TSTEP2,TIMSEC
	COMMON /CMHEAT/ ENTRA
	COMMON /CMLINE/ LINPAG,LINCNT,LINRUN,LINEQU,
	1		LINSCL
	REAL MASS
	COMMON /CMPARM/ MASS,EPSILN,SIGMA,SIGPAR,RCUTS,RCUTU
	COMMON /CMPROP/ TOTE,TEMP
	COMMON /CMSCAL/ NSCALE,DTEMP,VSCALE
	COMMON /CMSTEP/ NSTEP,FNSTEP,NRJECT,MAXSTP,
	1		NPRINT,NPROP
	COMMON /CMSTPT/ RHO,TR,VOLM,T
	COMMON /CMSUMV/ SUMX1,SUMY1,SUMZ1,SUMV2
	COMMON /CMTAIL/ ULTAIL,VLTAIL
C
C ===== PRINT OUT PAGE HEADER
C
	WRITE(5,5) NSTEP
5	FORMAT(1H1,////,T5,24('*'),' AVERAGES OVER',I5,' STEPS ',
	1  24('*'),///)
C
C ===== PRINT OUT RUNNING AVERAGES SO FAR
C
	EINT = SUMTEN/FNSTEP
	ENERGY = SUMENR/FNSTEP
	VIRIAL = SUMVIR/FNSTEP
	VTMEAN = SUMVT/FNSTEP
	TEMPK = SUMTMP/FNSTEP
	TEMPA = BOLTZ*TEMPK/EPSILN
	F2 = SUMFC2/FNSTEP
	SDEINT = SQRT(ABS(SSQTEN - FNSTEP*EINT*EINT)/FNSTEP)
	SDFRCE = SQRT(ABS(SSQFC2 - FNSTEP*F2*F2)/FNSTEP)
	SDTMPK = SQRT(ABS(SSQTMP - FNSTEP*TEMPK*TEMPK)/FNSTEP)
	SDTMP = BOLTZ*SDTMPK/EPSILN
	SDUCF = SQRT(ABS(SSQENR - FNSTEP*ENERGY*ENERGY)/FNSTEP)
	SDVRAA = SQRT(ABS(SSQVIR - FNSTEP*VIRIAL*VIRIAL)/FNSTEP)
	PRESS = (GASCON*TEMPK - VIRIAL/3.0)/VOLM
	SDPRSS = SQRT(ABS(GASCON*SDTMPK**2 + SDVRAA**2/9.0
	1   - 2.0*GASCON/3.0*(VTMEAN - TEMPK*VIRIAL)))/VOLM
	WRITE(5,10) TEMPA,SDTMP,TEMPK,SDTMPK
	WRITE(5,20) ENERGY,SDUCF
	WRITE(5,25) VIRIAL,SDVRAA,PRESS,SDPRSS
	WRITE(5,30) EINT,SDEINT
10	FORMAT(/,T5,'Reduced Temperature =',F8.3,' +/-',F8.3,
	1  '  (',F6.2,' +/-',F6.2,' Kelvin)')
20	FORMAT(/,T5,'Mean Configurational Energy =',
	1  -3PF8.3,' +/- ',F8.3,' Kj/Mol')
25	FORMAT(/,T5,'Mean Virial =',-3PF10.3,' +/- ',F10.3,
	1 ' Kj/Mol',//,
	2  T5,'Pressure =',-5PF10.3,' +/-',F10.3,' Bars')
30	FORMAT(/,
	2  T5,'Total Energy =',-3PF9.4,' +/-',F9.4,' Kj/Mol')
	F2 = (EPSILN/ULEN)**2*F2
	SDFRCE = (EPSILN/ULEN)**2*SDFRCE
	WRITE(5,50) F2,SDFRCE
50	FORMAT(/,T5,'Mean-Squared Force =',1PE13.4,' +/-',
	1  E13.4,' (J/M)**2')
	CV = 1.5*GASCON/(1.0 - 2.0/3.0*108.0*(SDTMPK/TEMPK)**2)
	GAMMAV = CV/VOLM*(2.0/3.0 - 108.0/GASCON/TEMPK**2*(GASCON*
	1  SDTMPK**2 - (VTMEAN - TEMPK*VIRIAL)/3.0))
	WRITE(5,70) CV,GAMMAV
70	FORMAT(/,T5,'Constant Volume Specific Heat =',F10.3,
	1  ' J/Mol/K',//,
	2  T5,'Thermal Pressure Coefficient =',-5PF10.3)
	RETURN
	END
	SUBROUTINE ZERO
	DOUBLE PRECISION SUMTEN,SSQTEN
	COMMON /CMAVER/ SUMTMP,SSQTMP,SUMENR,SSQENR,SUMVIR,
	1		SSQVIR,SUMTEN,SSQTEN,SUMFC2,SSQFC2,
	2		SUMVT
	COMMON /CMBOX/ ULEN
	COMMON /CMCCTR/ CT0,CT1,CT2,CT3
	REAL NAOVRN
	COMMON /CMCNST/ PI,EFACT,AVOGAD,BOLTZ,GASCON,NAOVRN
	COMMON /CMCONF/ UCONF,VIRAA,ECONF,EVIRAA
	COMMON /CMDELT/ TSTEP,TSTEP2,TIMSEC
	COMMON /CMHEAT/ ENTRA
	COMMON /CMLINE/ LINPAG,LINCNT,LINRUN,LINEQU,
	1		LINSCL
	REAL MASS
	COMMON /CMPARM/ MASS,EPSILN,SIGMA,SIGPAR,RCUTS,RCUTU
	COMMON /CMPROP/ TOTE,TEMP
	COMMON /CMSCAL/ NSCALE,DTEMP,VSCALE
	COMMON /CMSTEP/ NSTEP,FNSTEP,NRJECT,MAXSTP,
	1		NPRINT,NPROP
	COMMON /CMSTPT/ RHO,TR,VOLM,T
	COMMON /CMSUMV/ SUMX1,SUMY1,SUMZ1,SUMV2
	COMMON /CMTAIL/ ULTAIL,VLTAIL
	COMMON /CMMOVE/ XMOVE(108),YMOVE(108),ZMOVE(108)
C
C ===== RESET TO ZERO ALL PROPERTY COUNTERS AT END OF EQUILIBRATION
C
	WRITE(5,10) NSCALE
10	FORMAT(///,T5,'*** End of Equilibration -- Above Time ',
	1  'steps have been rejected ***',//,
	2  T5,'*** Velocities were rescaled',I6,' times ***')
	SUMTMP = 0.0
	SSQTMP = 0.0
	SUMENR = 0.0
	SSQENR = 0.0
	SUMVIR = 0.0
	SSQVIR = 0.0
	SUMVT = 0.0
	SUMTEN = 0.0
	SSQTEN = 0.0
	SUMFC2 = 0.0
	SSQFC2 = 0.0
C
	NSTEP = 0
	FNSTEP = 0.0
	NRJECT = 0
C
C ===== ZERO C.O.M. DISPLACEMENTS
C
	DO 100 I = 1,108
	XMOVE(I) = 0.0
	YMOVE(I) = 0.0
	ZMOVE(I) = 0.0
100	CONTINUE
C
C ===== START NEW PAGE AND WRITE HEADER
C
	CALL HEADER
	RETURN
	END
	SUBROUTINE FORCEN
C
C ===== EVALUATE FORCES ACTING ON ALL ATOMS
C
	COMMON /CMCOM0/ X0(108),Y0(108),Z0(108)
	COMMON /CMFRCE/ FRCEX(108),FRCEY(108),FRCEZ(108)
	COMMON /CMBOX/ ULEN
	COMMON /CMCCTR/ CT0,CT1,CT2,CT3
	REAL NAOVRN
	COMMON /CMCNST/ PI,EFACT,AVOGAD,BOLTZ,GASCON,NAOVRN
	COMMON /CMCONF/ UCONF,VIRAA,ECONF,EVIRAA
	COMMON /CMDELT/ TSTEP,TSTEP2,TIMSEC
	COMMON /CMHEAT/ ENTRA
	COMMON /CMLINE/ LINPAG,LINCNT,LINRUN,LINEQU,
	1		LINSCL
	REAL MASS
	COMMON /CMPARM/ MASS,EPSILN,SIGMA,SIGPAR,RCUTS,RCUTU
	COMMON /CMPROP/ TOTE,TEMP
	COMMON /CMSCAL/ NSCALE,DTEMP,VSCALE
	COMMON /CMSTEP/ NSTEP,FNSTEP,NRJECT,MAXSTP,
	1		NPRINT,NPROP
	COMMON /CMSTPT/ RHO,TR,VOLM,T
	COMMON /CMSUMV/ SUMX1,SUMY1,SUMZ1,SUMV2
	COMMON /CMTAIL/ ULTAIL,VLTAIL
C
C ===== ZERO CONFIGURATIONAL PROPERTIES
C
	UCONF = 0.0
	VIRAA = 0.0
C
C ===== ZERO FORCE ACCUMULATORS
C
	DO 50 I = 1,108
	FRCEX(I) = 0.0
	FRCEY(I) = 0.0
	FRCEZ(I) = 0.0
50	CONTINUE
C
C ===== LOOP OVER ATOM I
C
	DO 100 I = 1,107
	XI = X0(I)
	YI = Y0(I)
	ZI = Z0(I)
	JBEG = I + 1
C
C ===== LOOP OVER ATOM J
C
	DO 200 J = JBEG,108
C
C ===== ATOM - ATOM SEPARATION
C
	XDIFF = XI - X0(J)
	XDIFF = XDIFF - 2*INT(XDIFF)
	YDIFF = YI - Y0(J)
	YDIFF = YDIFF - 2*INT(YDIFF)
	ZDIFF = ZI - Z0(J)
	ZDIFF = ZDIFF - 2*INT(ZDIFF)
	RSQ = XDIFF*XDIFF + YDIFF*YDIFF + ZDIFF*ZDIFF
	IF(RSQ.GT.RCUTU)GO TO 200
	RSQI = 1.0/RSQ
	R6 = (SIGPAR*RSQI)**3
	R12 = R6*R6
	AAVIR = R12 + R12 - R6
	DUDRR1 = AAVIR*RSQI
	UCONF = UCONF + R12 - R6
	VIRAA = VIRAA + AAVIR
	FXIJ = XDIFF*DUDRR1
	FYIJ = YDIFF*DUDRR1
	FZIJ = ZDIFF*DUDRR1
	FRCEX(I) = FRCEX(I) + FXIJ
	FRCEY(I) = FRCEY(I) + FYIJ
	FRCEZ(I) = FRCEZ(I) + FZIJ
	FRCEX(J) = FRCEX(J) - FXIJ
	FRCEY(J) = FRCEY(J) - FYIJ
	FRCEZ(J) = FRCEZ(J) - FZIJ
200	CONTINUE
100	CONTINUE
C
C ===== ALLOW FOR FACTOR OF 4.0 IN LENNARD-JONES POTENTIAL
C
	UCONF = 4.0*UCONF
	VIRAA = -24.0*VIRAA
	DO 500 I = 1,108
	FRCEX(I) = 24.0*FRCEX(I)
	FRCEY(I) = 24.0*FRCEY(I)
	FRCEZ(I) = 24.0*FRCEZ(I)
500	CONTINUE
	RETURN
	END
$LINK EXAMPLE
$DELETE EXAMPLE.OBJ;*
$RUN EXAMPLE
$DELETE EXAMPLE.TSK;*