C PROGRAM MDNEW.F
C
C 10/23/90 Switch the EQBRAT subroutine (remove the old one which runs
C NEQ+100 steps instead of NEQ during equilibration and replace
C with the one from 1989 version (md.f)) Tien
C 10/12/90 Adding INPUT DATA to the heading of the output file (in PRNT
C subroutine); change output file to mdout.new.
C Tien Nguyen
C
C DATE OF REVISION: OCT 13, 1988 Kin S. Cheung
C
C THIS IS A MODIFIED VERSION OF HAILE.F:
c
c - The maximum number of particle allowed is now 1000.
c - The equilibration criterion is changed from MSD to number of timesteps
c NEQ such that the program can accommodate solids.
c - The bug of np=108 in CORR is fixed.
c - Interactive input parameters include: NP, TR, DR, NEQ, MAXKB.
c - If NEQ = 0, IFLG is set to 1.
c - All output information is stored in file 'mdout.new'.
C
C DATE OF REVISION: Fall, 1991 Don Puffer
C
c - Additional output files have been added which contain only numbers.
C These can be read into gaphics packages or used to calculate
C parameters. The main file of interest is 'mdnew2.matlab.out'.
C This file contains all the numbers found in 'mdout.new' but
C without the text. This is the file use for the tutorial.
C
C
PARAMETER ( NPD = 1000, NP6 = NPD*6, NP9 = NPD*9 )
PARAMETER ( KBSKP = 10, ISTEP = 10 )
COMMON /POS/ X0(NPD),Y0(NPD),Z0(NPD)
COMMON /VEL/ X1(NPD),Y1(NPD),Z1(NPD)
COMMON /DER/ X2(NPD),Y2(NPD),Z2(NPD),X3(NPD),Y3(NPD),Z3(NPD),
A X4(NPD),Y4(NPD),Z4(NPD),X5(NPD),Y5(NPD),Z5(NPD)
COMMON /FOR/ FX(NPD),FY(NPD),FZ(NPD)
COMMON /DISP/ DAX(NPD),DAY(NPD),DAZ(NPD),X0L(NPD),Y0L(NPD),
A Z0L(NPD)
COMMON /NABLST/ LIST(NPD*50),NABORS(NPD)
COMMON /PARM/ F02,F12,F32,F42,F52
COMMON /PROP/ IDIST(1000),SUME,SUMV,XSUM,ISUM
COMMON /PROPA/ RDEL,RMAX,RDMAX,RLIST,FSHFT,ESHFT,ESHFTA,NP1,
A NP2,KSORT
COMMON /ISF/ XX0(NPD),YY0(NPD),ZZ0(NPD)
COMMON /ISCFT/ TMP1(NPD),TMP2(NPD),TMP3(NPD),TMP4(NPD),TMP5(NPD),
A TMP6(NPD)
COMMON /CORTMP/ TX0(NPD,ISTEP+1),TY0(NPD,ISTEP+1),TZ0(NPD,ISTEP+1)
DIMENSION IVALS(50),RVALS(50)
C
C ASSIGN PRINT UNIT
C -----------------
c
OPEN (UNIT=23, FORM='FORMATTED', FILE='mdout2.new')
LP = 23
C
c*****************************Puffer*********************************
C
OPEN (UNIT=24, FORM='FORMATTED', FILE='mdnew2.matlab.out')
matlab = 24
OPEN (UNIT=25, FORM='FORMATTED', FILE='md.force.matlab.out')
OPEN (UNIT=22, FORM='FORMATTED', FILE='md.rdf.matlab.out')
OPEN (UNIT=26, FORM='FORMATTED', FILE='md.position.matlab.out')
c OPEN (UNIT=26, FORM='FORMATTED', FILE='md.nbrdis.matlab.out')
C
c*****************************Puffer*********************************
C
C
DO 30 I=1,50
IVALS(I) = 0
RVALS(I) = 0.0
30 CONTINUE
C
C ASSIGN OUTPUT FILE LOGICAL UNIT NUMBER
C --------------------------------------
LU = 9
C
C NUMBER OF PARTICLES
C
PRINT*, "ENTER NUMBER OF PARTICLES (<=1000) : NP = "
READ(*,*)NP
PRINT*, "ENTER NUMBER OF EQUILIBRATING TIMESTEPS : NEQ = "
READ(*,*)NEQ
PRINT*, "ENTER NUMBER OF TIMESTEPS : MAXKB = "
READ(*,*)MAXKB
C
c31 FORMAT(I4)
c NP=108
C
C ******DUMMY VARIABLES FOR FORCE OUTPUT FILE TO SEPARATE DATA BLOCKS***
zeroI=0
zeroF=0.0
C
C TEST NUMBER OF PARTICLES VS. DIMENSION
C
IF ( NP .GT. NPD ) THEN
WRITE ( LP,70 ) NP,NPD
70 FORMAT ( /,' ERROR - NUMBER OF PARTICLES (',I4,')',/,
1 ' EXCEED PROGRAM DIMENSION (',I4,')' )
CALL EXIT
ENDIF
PART=NP
NP1=NP-1
NP2=.5*PART+.01
NP22=NP-2
C
C READ IN FLUID STATE CONDITIONS FOR KRYPTON
C
PRINT*, "ENTER DESIRED REDUCED TEMPERATURE : TR = "
READ(*,1)TR
1 FORMAT(E10.3)
C
PRINT*, "ENTER REDUCED DENSITY : DR ="
READ(*,2)DR
2 FORMAT(E10.3)
C
C
RC = 2.5
C
C SET RUN FLAGS AND PARAMETERS
IFLG=0
KB=0
KSAVE=10
KSORT=10
KWRITE=10
SEED =0.597
c
IF(NEQ.EQ.0) IFLG=1
C
NT = (MAXKB-KBSKP)/ISTEP
C IF ( NT .GT. NTD ) THEN
C WRITE ( LP,90 ) NT,NTD
C90 FORMAT ( /,' ERROR - NUMBER OF SAVED TIME STEPS (',I4,')',/,
C 1 ' EXCEED PROGRAM DIMENSION (',I4,')' )
C CALL EXIT
C ENDIF
TDIST=0.
XDIST=0.36
C
C SET VALUES OF PHYSICAL CONSTANTS FOR KRYPTON
SIGMA=3.405
EPSI=120.
AVO=6.0225E+23
BOLZ=1.38054
WTMOL=39.948
THIRD=1./3.
PI=3.14159265
C
C CALCULATE DIMENSIONAL TEMPERATURE AND VOLUME
T=TR*EPSI
VOL=PART/DR
VSCALE=THIRD/TR
VOLCC=AVO*(SIGMA*1.E-8)**3/DR
C
C SET TIME-STEP AND ITS MULTIPLES
DELTA=0.005
DELSQ=DELTA*DELTA
DELTSQ=.5*DELSQ
TSTEP= SQRT(WTMOL*SIGMA**2/AVO/EPSI/BOLZ)*1.E+12
TT1=TSTEP*DELTA*1.E-12
C
C SET PARAMETERS IN PREDICTOR-CORRECTOR METHOD
F02=3./16.
F12=251./360.
F32=11./18.
F42=1./6.
F52=1./60.
C
C SET DISTANCES FOR POTENTIAL CUT-OFF, VERLET LIST, ETC.
CUBE=VOL**THIRD
CUBE2=.5*CUBE
RLIST=(RC+.3)**2
RDMAX=CUBE2*CUBE2
IF(RLIST.GT.RDMAX) RLIST=RDMAX
IF(RC.GT.CUBE2) RC=CUBE2-0.35
RMAX=RC*RC
C
C SCALE FACTORS FOR VELOCITIES DURING EQUILIBRATION
AHEAT=DELSQ*PART*3.*TR
RTR=DELTA* SQRT(3.*TR)
C
C INCREMENT FOR SAMPLING FOR G(R)
RDEL=0.025
NY=CUBE2/RDEL-1
C
C
C SHIFTED-FORCE CONSTANTS
RRMAX=1./ SQRT(RMAX)
RRMAX6=RRMAX**6
ESHFT=RRMAX6*(28.-52.*RRMAX6)
ESHFTA=48.*RRMAX*RRMAX6*(RRMAX6-.5)
FSHFT=ESHFTA
C
C CORRECTIONS FOR LONG-RANGE INTERACTIONS
RC3=RC**3
RC9=RC3**3
CORE=8.*PI*DR*(1./9./RC9-THIRD/RC3)
CORV=96.*PI*DR*(.5*THIRD/RC3-1./9./RC9)
DE=CORE
DP=-VSCALE*CORV
C
C INITIALIZE SUM ACCUMMULATORS
ISUM=0
XSUM=0.
SUME=0.
SUMV=0.
DO 510 K=1,NY
510 IDIST(K)=0
C
C PRINT PARAMETERS
CALL PRNT(LP,NP,MAXKB,EPSI,SIGMA,TR,T,DR,VOL,CUBE,RMAX,RC,
A RDEL,RLIST,RDMAX,DELTA,NY,TSTEP,TT1,DE,DP,SEED,KBSKP,ISTEP,NEQ)
C
C PRINT RUN-TABLE HEADING
WRITE(LP,930)
930 FORMAT(1H1///7X,'KB',2X,'TIME',5X,'PRES',7X,'P1',5X,'ENRG'
A ,6X,'E1',7X,'DIST',4X,'TEMP',3X,'ETOTAL'/)
C
C
C LOAD INITIAL POSITIONS OF ATOMS
CALL FCC(CUBE,NP)
C
C LOAD INITIAL VELOCITIES OF ATOMS
CALL INTVEL(AHEAT,RTR,PART,NP,SEED)
C
C ASSIGN INITIAL ACCELERATIONS BASED ON INITIAL POSITIONS
CALL EVAL(TOTV,TOTE,CUBE,CUBE2,KB,NP,NABTOT)
C
C SCALE ACCELERATIONS AND STORE STARTING POSITIONS
DO 530 I=1,NP
X2(I)=FX(I)*DELTSQ
Y2(I)=FY(I)*DELTSQ
Z2(I)=FZ(I)*DELTSQ
X0L(I)=X0(I)
Y0L(I)=Y0(I)
Z0L(I)=Z0(I)
530 CONTINUE
C
C
C
C
C OUTPUT HEADER INFO TO UNIT 9 FOR FUTURE USE
C -------------------------------------------
C
C
IVALS(01) = NP
IVALS(02) = NT
IVALS(03) = NY
IVALS(04) = MAXKB
IVALS(05) = KBSKP
IVALS(06) = ISTEP
C
RVALS(01) = EPSI
RVALS(02) = SIGMA
RVALS(03) = TR
RVALS(04) = T
RVALS(05) = DR
RVALS(06) = VOL
RVALS(07) = CUBE
RVALS(08) = RMAX
RVALS(09) = RC
RVALS(10) = RDEL
RVALS(11) = RLIST
RVALS(12) = RDMAX
RVALS(13) = DELTA
RVALS(14) = TSTEP
RVALS(15) = TT1
RVALS(16) = DE
RVALS(17) = DP
RVALS(18) = SEED
C
CXXX REWIND ( LU )
CXXX WRITE ( LU ) IVALS,RVALS
C
C
C
C
C ENTER MAIN LOOP OF SIMULATION
C
NTSTP =0
NSTART=1
533 DO 599 NTIMES=NSTART,MAXKB
KB=KB+1
CALL PREDCT(NP)
CALL EVAL(TOTV,TOTE,CUBE,CUBE2,KB,NP,NABTOT)
CALL CORR(LP,DELTSQ,CUBE,KB,IFLG,NTSTP,LU,NP)
C
C CALCULATE MEAN SQUARE DISPLACEMENT & KINETIC ENERGY
SUMVEL=0.
TDIST=0.
DO 540 I=1,NP
TDIST=TDIST+DAX(I)**2+DAY(I)**2+DAZ(I)**2
SUMVEL=SUMVEL+X1(I)**2+Y1(I)**2+Z1(I)**2
540 CONTINUE
TDIST=TDIST/PART
EK=SUMVEL/(2.*PART*DELSQ)
C
C ACCUMMULATE SUMS FOR PROPERTY AVERAGES
XSUM=XSUM+SUMVEL
SUME=SUME+TOTE
SUMV=SUMV+TOTV
C
C PROPERTY CALCULATION & PRINT-OUT AT INTERVALS
IF(MOD(KB,KWRITE).NE.0) GO TO 550
FKB= FLOAT(KB)*PART
TMP=XSUM/(3.*DELSQ*FKB)
ENR=(SUME/FKB+CORE)
VIR=VSCALE*(SUMV/FKB+CORV)
PRES=TMP*DR*(1.-VIR)
E1=TOTE/PART+CORE
ETOT=E1+EK
P1=TMP*DR*(1.-VSCALE*(TOTV/PART+CORV))
RLTIM=DELTA* FLOAT(KB)*TSTEP
WRITE(LP,940) KB,RLTIM,PRES,P1,ENR,E1,TDIST,TMP,ETOT
940 FORMAT(1H ,3X,I5,F7.3,5F9.3,F9.4,2X,F8.5)
WRITE(matlab,941) KB,RLTIM,PRES,P1,ENR,E1,TDIST,TMP,ETOT
941 FORMAT(1H ,3X,I5,F7.3,5F9.3,F9.4,2X,F8.5)
DO 005 II=1,NP
C chkforce=FX(II)+FY(II)+FZ(II)
c if (KB.GE.100) then
c if(KB.LE.150) then
c WRITE(25,222) II, FX(II),FY(II),FZ(II)
C WRITE(26,222) II, X0(II),Y0(II),Z0(II)
222 FORMAT(I5,F11.6,F11.6,F11.6)
c endif
c endif
005 CONTINUE
C WRITE(25,222) zeroI,zeroF,zeroF,zeroF
C
C PRINT G(R) AT INTERVALS
IF(MOD(KB,1000).EQ.0.OR.KB.EQ.MAXKB) CALL RDF(LP,DR,TMP,
A RDEL,NP,NY,KB)
C
C DURING FIRST OF RUN, SCALE VELOCITIES FOR TEMPERATURE
550 IF(IFLG.LT.1) CALL EQBRAT(LP,SUMVEL,AHEAT,TDIST,XDIST,NP,
A IFLG,NY,KB,NEQ)
599 CONTINUE
IF(KB.GE.MAXKB) GO TO 600
NSTART=KB+1
GO TO 533
600 CONTINUE
C
C NOTIFY USER OF NUMBER OF TIME STEPS ACTUALLY SAVED
C --------------------------------------------------
WRITE ( LP,950 ) MAXKB,KBSKP,ISTEP,NTSTP
950 FORMAT ( ///,' SUMMARY :',/,' *********',//,
1 ' STARTING TIME STEP : 1',/,
2 ' ENDING TIME STEP :',I6,/,
3 ' SKIP THROUGH TIME STEP :',I6,/,
4 ' SAVE TIME STEP INCREMENT :',I6,/,
5 ' ACTUAL NUMBER OF TIME STEPS SAVED :',I6 )
C
C === END OF EXECUTION
C
END
SUBROUTINE FCC(CUBE,NP)
C
C CALCULATE POSITIONS OF SITES ON FACE-CENTERED CUBIC LATTICE
C BASED ON A BOX OF SIDE = 100
C
PARAMETER ( NPD = 1000, NP6 = NPD*6, NP9 = NPD*9 )
COMMON /POS/ X0(NPD),Y0(NPD),Z0(NPD)
NC=( FLOAT(NP)/4.)**(1./3.)+.1
XL=100./ FLOAT(NC)
Y=.5*XL
X0(1)=0.
Y0(1)=0.
Z0(1)=0.
X0(2)=0.
Y0(2)=Y
Z0(2)=Y
X0(3)=Y
Y0(3)=0.
Z0(3)=Y
X0(4)=Y
Y0(4)=Y
Z0(4)=0.
M=0
DO 10 I=1,NC
XLI = XL*(I-1)
DO 10 J=1,NC
XLJ = XL*(J-1)
DO 10 K=1,NC
XLK = XL*(K-1)
DO 11 IJ=1,4
X0(IJ+M)=X0(IJ)+XLK
Y0(IJ+M)=Y0(IJ)+XLJ
Z0(IJ+M)=Z0(IJ)+XLI
11 CONTINUE
M=M+4
10 CONTINUE
C
C SCALE POSITIONS TO BOX OF SIDE = CUBE
CUB01 = CUBE*.01
DO 12 I=1,NP
X0(I)=X0(I)*CUB01
Y0(I)=Y0(I)*CUB01
Z0(I)=Z0(I)*CUB01
12 CONTINUE
RETURN
END
SUBROUTINE INTVEL(AHEAT,RTR,PART,NP,SEED)
C
C ASSIGN INITIAL VELOCITIES TO ATOMS
C
PARAMETER ( NPD = 1000, NP6 = NPD*6, NP9 = NPD*9 )
COMMON/VEL/X1(NPD),Y1(NPD),Z1(NPD)
CXXX IRR=34567
X=SEED
SUMX=0.
SUMY=0.
SUMZ=0.
DO 200 I=1,NP
CXXX IRR=IRR*65539
CXXX XX= FLOAT(IRR)*.4656613D-9
CALL RANDOM(X,XX)
X=XX
XX=2.0*(XX-0.5)
CXXX IRR=IRR*65539
CXXX YY= FLOAT(IRR)*.4656613D-9
CALL RANDOM(X,YY)
X=YY
YY=2.0*(YY-0.5)
CXXX IRR=IRR*65539
CXXX ZZ= FLOAT(IRR)*.4656613D-9
CALL RANDOM(X,ZZ)
X=ZZ
ZZ=2.0*(ZZ-0.5)
XYZ=1./ SQRT(XX*XX+YY*YY+ZZ*ZZ)
X1(I)=XX*XYZ*RTR
Y1(I)=YY*XYZ*RTR
Z1(I)=ZZ*XYZ*RTR
SUMX=SUMX+X1(I)
SUMY=SUMY+Y1(I)
SUMZ=SUMZ+Z1(I)
200 CONTINUE
C
C SCALE VELOCITIES SO THAT TOTAL MOMENTUM=ZERO
X=0
DO 210 I=1,NP
X1(I)=X1(I)-SUMX/PART
Y1(I)=Y1(I)-SUMY/PART
Z1(I)=Z1(I)-SUMZ/PART
X=X+X1(I)**2+Y1(I)**2+Z1(I)**2
210 CONTINUE
C
C SCALE VELOCITIES TO DESIRED TEMPERATURE
HEAT= SQRT(AHEAT/X)
Heat=1.
DO 220 I=1,NP
X1(I)=X1(I)*HEAT
Y1(I)=Y1(I)*HEAT
Z1(I)=Z1(I)*HEAT
220 CONTINUE
RETURN
END
SUBROUTINE PREDCT(NP)
C
C USE FIFTH-ORDER TAYLOR SERIES TO PREDICT POSITIONS & THEIR
C DERIVATIVES AT NEXT TIME-STEP
C
PARAMETER ( NPD = 1000, NP6 = NPD*6, NP9 = NPD*9 )
COMMON/POS/X0(NPD),Y0(NPD),Z0(NPD)
COMMON/VEL/X1(NPD),Y1(NPD),Z1(NPD)
COMMON/DER/X2(NPD),Y2(NPD),Z2(NPD),X3(NPD),Y3(NPD),Z3(NPD)
A ,X4(NPD),Y4(NPD),Z4(NPD),X5(NPD),Y5(NPD),Z5(NPD)
COMMON/FOR/FX(NPD),FY(NPD),FZ(NPD)
C
C
DO 300 I=1,NP
X0(I)=X0(I)+X1(I)+X2(I)+X3(I)+X4(I)+X5(I)
Y0(I)=Y0(I)+Y1(I)+Y2(I)+Y3(I)+Y4(I)+Y5(I)
Z0(I)=Z0(I)+Z1(I)+Z2(I)+Z3(I)+Z4(I)+Z5(I)
X1(I)=X1(I)+2.*X2(I)+3.*X3(I)+4.*X4(I)+5.*X5(I)
Y1(I)=Y1(I)+2.*Y2(I)+3.*Y3(I)+4.*Y4(I)+5.*Y5(I)
Z1(I)=Z1(I)+2.*Z2(I)+3.*Z3(I)+4.*Z4(I)+5.*Z5(I)
X2(I)=X2(I)+3.*X3(I)+6.*X4(I)+10.*X5(I)
Y2(I)=Y2(I)+3.*Y3(I)+6.*Y4(I)+10.*Y5(I)
Z2(I)=Z2(I)+3.*Z3(I)+6.*Z4(I)+10.*Z5(I)
X3(I)=X3(I)+4.*X4(I)+10.*X5(I)
Y3(I)=Y3(I)+4.*Y4(I)+10.*Y5(I)
Z3(I)=Z3(I)+4.*Z4(I)+10.*Z5(I)
X4(I)=X4(I)+5.*X5(I)
Y4(I)=Y4(I)+5.*Y5(I)
Z4(I)=Z4(I)+5.*Z5(I)
C
C INITIALIZE FORCE ARRAYS
FX(I)=0.
FY(I)=0.
FZ(I)=0.
300 CONTINUE
RETURN
END
SUBROUTINE EVAL(TOTV,TOTE,CUBE,CUBE2,KB,NP,NABTOT)
C
C EVALUATE FORCES ON ATOMS USING PAIRWISE ADDITIVE LENNARD-
C JONES (6-12) INTERMOLECULAR POTENTIAL
C
PARAMETER ( NPD = 1000, NP6 = NPD*6, NP9 = NPD*9 )
COMMON/POS/X0(NPD),Y0(NPD),Z0(NPD)
COMMON/FOR/FX(NPD),FY(NPD),FZ(NPD)
COMMON/NABLST/LIST(NPD*50),NABORS(NPD)
COMMON/PROP/IDIST(1000),SUME,SUMV,XSUM,ISUM
COMMON/PROPA/RDEL,RMAX,RDMAX,RLIST,FSHFT,ESHFT,ESHFTA,NP1,
A NP2,KSORT
C
C
C STATEMENT FUNCTIONS
C -------------------
RPLFUN(R,RC) = 48.*R*(R-0.5)-RC*FSHFT
TOTFUN(R,RC) = 4.*R*(R-1.0)+RC*ESHFTA+ESHFT
XYZIMG(W,WR) = (WR-CUBE)*SIGN(1.00,W)
C
K=0
TOTE=0.
TOTV=0.
C
C SET FLAG FOR LIST UP-DATE
LCHK=2
IF(MOD(KB,KSORT).EQ.0) LCHK=1
C
C*********Puffer mod**********************
IF(MOD(KB,KSORT).EQ.0) then
WRITE(26,220) kb
220 FORMAT(I5)
endif
C******************************************
C OUTER LOOP OVER ATOMS
DO 490 I=1,NP1
IF (LCHK .EQ. 1) GO TO 100
C
C USE NEIGHBOR-LIST
JBEGIN = NABORS(I)
JEND = NABORS(I+1) - 1
IF (JBEGIN .GT. JEND) GO TO 490
GO TO 150
C
C LOOP OVER ALL ATOMS AT KSORT INTERVALS
100 NABORS(I)=K+1
JBEGIN=I+1
JEND=NP
C
C STORE POSITION OF ATOM I
150 XI=X0(I)
YI=Y0(I)
ZI=Z0(I)
C
C INNER LOOPS OVER ATOMS
C
DO 495 JX = JBEGIN,JEND
J = JX
IF (LCHK .EQ. 2) J = LIST(JX)
C
C DISTANCE BETWEEN ATOM I AND ATOM J
X = XI - X0(J)
Y = YI - Y0(J)
Z = ZI - Z0(J)
C
C FIND IMAGE OF ATOM J CLOSEST TO ATOM I
C
XR = ABS ( X )
YR = ABS ( Y )
ZR = ABS ( Z )
IF ( XR .GT. CUBE2 ) X=XYZIMG(X,XR)
IF ( YR .GT. CUBE2 ) Y=XYZIMG(Y,YR)
IF ( ZR .GT. CUBE2 ) Z=XYZIMG(Z,ZR)
C
C CALCULATE R(IJ)
C
RSQ = X**2 + Y**2 + Z**2
RCTR = SQRT ( RSQ )
C
C INCREMENT COUNTER FOR G(R) & UPDATE LIST AT INTERVALS
IF (LCHK .NE. 1) GO TO 175
IF ( RSQ .GT. RDMAX ) GO TO 495
IJ = RCTR / RDEL + .5
IDIST(IJ)=IDIST(IJ)+1
C
IF ( RSQ .GT. RLIST ) GO TO 495
K=K+1
LIST(K)=J
C
175 IF ( RSQ .GT. RMAX ) GO TO 495
C
C EVALUATE FORCES ON ATOMS
RSI=1./RSQ
R6=RSI**3
RPL = RPLFUN(R6,RCTR)
RP1=RPL*RSI*X
RP2=RPL*RSI*Y
RP3=RPL*RSI*Z
FX(I)=FX(I)+RP1
FX(J)=FX(J)-RP1
FY(I)=FY(I)+RP2
FY(J)=FY(J)-RP2
FZ(I)=FZ(I)+RP3
FZ(J)=FZ(J)-RP3
C
C ACCUMMULATE ENERGY AND VIRIAL
C******** puffer modification********
energyfun=totfun(r6,rctr)
TOTE = TOTE + energyfun
C TOTE = TOTE + TOTFUN(R6,RCTR)
C if (KB.GE.10) then
IF(MOD(KB,KSORT).EQ.0) then
C IF(MOD(KB,KWRITE).EQ.0) then
if (abs(energyfun).GE.3) then
c WRITE(26,220) kb
c220 FORMAT(I5)
WRITE(26,221) I,J,energyfun,RCTR,TOTE
221 FORMAT(I5,I5,f11.6,f11.6,f11.6)
endif
endif
C endif
C*************************************
TOTV=TOTV-RPL
495 CONTINUE
490 CONTINUE
C
IF(LCHK.NE.1) RETURN
NABORS(NP)=K+1
ISUM=ISUM+NP2
NABTOT=K
C
RETURN
END
SUBROUTINE CORR(LP,DELTSQ,CUBE,KB,IFLG,NTSTP,LU,NP)
C
C CORRECT PREDICTED POSITIONS AND THEIR DERIVATIVES
C
PARAMETER ( NPD = 1000, NP6 = NPD*6, NP9 = NPD*9 )
PARAMETER ( KBSKP = 10, ISTEP = 10 )
COMMON/POS/X0(NPD),Y0(NPD),Z0(NPD)
COMMON/VEL/X1(NPD),Y1(NPD),Z1(NPD)
COMMON/DER/X2(NPD),Y2(NPD),Z2(NPD),X3(NPD),Y3(NPD),Z3(NPD)
A ,X4(NPD),Y4(NPD),Z4(NPD),X5(NPD),Y5(NPD),Z5(NPD)
COMMON/FOR/FX(NPD),FY(NPD),FZ(NPD)
COMMON/DISP/DAX(NPD),DAY(NPD),DAZ(NPD),X0L(NPD),Y0L(NPD),
A Z0L(NPD)
COMMON/PARM/F02,F12,F32,F42,F52
COMMON/ISF/XX0(NPD),YY0(NPD),ZZ0(NPD)
COMMON/ISCFT/TMP1(NPD),TMP2(NPD),TMP3(NPD),TMP4(NPD),TMP5(NPD),
A TMP6(NPD)
COMMON/CORTMP/TX0(NPD,ISTEP+1),TY0(NPD,ISTEP+1),TZ0(NPD,ISTEP+1)
C
DO 670 I=1,NP
XERR=X2(I)-DELTSQ*FX(I)
YERR=Y2(I)-DELTSQ*FY(I)
ZERR=Z2(I)-DELTSQ*FZ(I)
X0(I)=X0(I)-XERR*F02
X1(I)=X1(I)-XERR*F12
X2(I)=X2(I)-XERR
X3(I)=X3(I)-XERR*F32
X4(I)=X4(I)-XERR*F42
X5(I)=X5(I)-XERR*F52
Y0(I)=Y0(I)-YERR*F02
Y1(I)=Y1(I)-YERR*F12
Y2(I)=Y2(I)-YERR
Y3(I)=Y3(I)-YERR*F32
Y4(I)=Y4(I)-YERR*F42
Y5(I)=Y5(I)-YERR*F52
Z0(I)=Z0(I)-ZERR*F02
Z1(I)=Z1(I)-ZERR*F12
Z2(I)=Z2(I)-ZERR
Z3(I)=Z3(I)-ZERR*F32
Z4(I)=Z4(I)-ZERR*F42
Z5(I)=Z5(I)-ZERR*F52
C
C DISPLACEMENTS
DAX(I)=DAX(I)-X0(I)+X0L(I)
DAY(I)=DAY(I)-Y0(I)+Y0L(I)
DAZ(I)=DAZ(I)-Z0(I)+Z0L(I)
670 CONTINUE
C
C CHECK FOR CONSERVATION OF PARTICLES IN PRIMARY CELL
C
IND = KB - KC*ISTEP + 1
IF ( KB .LE. ISTEP ) IND = KB
C WRITE ( LP,777 ) KB,IND
777 FORMAT ( ' KB:',I4,' IND:',I3 )
C
DO 690 I=1,NP
XX=X0(I)
YY=Y0(I)
ZZ=Z0(I)
C
IF(IFLG .LT. 1) GO TO 1000
IF(IFLG .EQ. 1 .AND. KB .EQ. 1) GO TO 800
IF(IFLG .EQ. 1 .AND. IND .GE. 2) GO TO 900
GO TO 1000
800 TX0(I,1)=X0(I)
TY0(I,1)=Y0(I)
TZ0(I,1)=Z0(I)
IND=1
GO TO 1000
C
900 TX0(I,IND)=TX0(I,IND-1)+(X0(I)-X0L(I))
TY0(I,IND)=TY0(I,IND-1)+(Y0(I)-Y0L(I))
TZ0(I,IND)=TZ0(I,IND-1)+(Z0(I)-Z0L(I))
C
IF(MOD(KB,ISTEP) .NE. 0 ) GO TO 1000
KC=KB/ISTEP
IF ( KB .LE. KBSKP ) GO TO 950
XX0(I)=TX0(I,IND)
YY0(I)=TY0(I,IND)
ZZ0(I)=TZ0(I,IND)
950 TX0(I,1) =TX0(I,IND)
TY0(I,1) =TY0(I,IND)
TZ0(I,1) =TZ0(I,IND)
C
1000 CONTINUE
IF(XX.LE.0.) XX=XX+CUBE
IF(XX.GE.CUBE) XX=XX-CUBE
IF(YY.LE.0.) YY=YY+CUBE
IF(YY.GE.CUBE) YY=YY-CUBE
IF(ZZ.LE.0.) ZZ=ZZ+CUBE
IF(ZZ.GE.CUBE) ZZ=ZZ-CUBE
C
C
C STORE NEW POSITIONS
X0(I)=XX
Y0(I)=YY
Z0(I)=ZZ
X0L(I)=XX
Y0L(I)=YY
Z0L(I)=ZZ
690 CONTINUE
C
C
C SAVE POSITIONS AND VELOCITIES, IMPLICIT TIME STEP
C -------------------------------------------------
C
IF(MOD(KB,ISTEP) .NE. 0 .OR. KB .LE. KBSKP ) GO TO 2000
NTSTP = NTSTP + 1
CXXX WRITE ( LU ) XX0,YY0,ZZ0,X1,Y1,Z1
C************************************************************************
CXXX IF ( NTSTP .LE. 10 )
CXXX 1 WRITE ( LP,888 ) NTSTP,XX0(1),YY0(1),ZZ0(1),X1(1),Y1(1),Z1(1),
CXXX 2 XX0(NP),YY0(NP),ZZ0(NP),X1(NP),Y1(NP),Z1(NP)
888 FORMAT ( I8,6F15.9,/,8X,6F15.9 )
C************************************************************************
C
2000 RETURN
END
SUBROUTINE RDF(LP,DR,TMP,RDEL,NP,NY,KB)
C
C NORMALIZE COUNTERS FOR RADIAL DISTRIBUTION FUNCTION
C
COMMON/PROP/IDIST(1000),SUME,SUMV,XSUM,ISUM
C
C PRINT HEADING
WRITE(LP,963) KB
963 FORMAT(1H1////T20,'RADIAL DISTRIBUTION FUNCTION AT TIME-',
A 'STEP',I6,/)
WRITE(LP,965) NP,DR,TMP
965 FORMAT(T25,'NP = ',I4,3X,'DR = ',F6.3,3X,'TR = ',F6.3//)
WRITE(LP,968)
968 FORMAT(29X,'I',6X,'R',6X,'IDIST',3X,'G(R)'/)
C
Y=ISUM
PDEN=DR*(RDEL**3)
DO 780 J=1,NY
IF(IDIST(J).EQ.0) GO TO 780
V=(4.*3.14159/3.)*(3.* FLOAT(J)**2+.25)
X=IDIST(J)
GR=X/V/Y/PDEN
RRR=RDEL* FLOAT(J)
WRITE(LP,972) J,RRR,IDIST(J),GR
972 FORMAT(28X,I3,3X,F5.3,3X,I6,F8.3)
c************************************************puffer modification*******************
WRITE(22,223) J,RRR,IDIST(J),GR
223 FORMAT(I3,3X,F5.3,3X,I6,F8.3)
c************************************************puffer modification*******************
780 CONTINUE
C
WRITE(LP,970)
970 FORMAT(1H1///)
RETURN
END
c
c************************************
SUBROUTINE EQBRAT(LP,SUMVEL,AHEAT,TDIST,XDIST,NP,IFLG,NY,KB,NEQ)
C
C SCALE VELOCITIES DURING INITIAL TIME-STEPS
C
PARAMETER ( NPD = 1000, NP6 = NPD*6, NP9 = NPD*9 )
COMMON/VEL/X1(NPD),Y1(NPD),Z1(NPD)
COMMON/PROP/IDIST(1000),SUME,SUMV,XSUM,ISUM
C
C
c To perform vel rescaling during the first NEQ timesteps:
c
HEAT= SQRT(AHEAT/SUMVEL)
SUMVELN = 0.0
DO 730 I=1,NP
X1(I)=X1(I)*HEAT
Y1(I)=Y1(I)*HEAT
Z1(I)=Z1(I)*HEAT
SUMVELN=SUMVELN+X1(I)**2+Y1(I)**2+Z1(I)**2
730 CONTINUE
c
sumvel = sumveln
c
c Reset all sum counters to zero at the end of equilibration
c
if(kb.eq.neq) then
iflg = 1
kb = 0
isum = 0
SUME=0.
SUMV=0.
XSUM=0.
CTT ************************
DO 100 I=1,1000
100 IDIST(I)=0
CTT ************************
WRITE(lp,777)
777 FORMAT(/10X,'**** EQUILIBRATION COMPLETE ****'//)
WRITE(lp,930)
930 FORMAT(1H1///7X,'KB',2X,'TIME',5X,'PRES',7X,'P1',5X,'ENRG'
A ,6X,'E1',7X,'DIST',4X,'TEMP',3X,'ETOTAL'/)
endif
RETURN
C
END
c
SUBROUTINE RANDOM(X,XR)
C
C SUBROUTINE TO GENRATE RANDOM NUMBER.
C
C
DATA K,J,M,RM / 5701,3612,566927,566927.0 /
C
IX=INT(X*RM)
KR=J*IX+K
IRAND=MOD(KR,M)
XR=(FLOAT(IRAND) + 0.5) / RM
RETURN
END
SUBROUTINE PRNT(LP,NP,MAXKB,EPSI,SIGMA,TR,T,DR,VOL,CUBE,RMAX,RC,
A RDEL,RLIST,RDMAX,DELTA,NY,TSTEP,TT1,DE,DP,SEED,KBSKP,ISTEP,
A NEQ)
C
C PRINT PARAMETERS
C
C
WRITE(LP,900)
900 FORMAT(1H1///)
WRITE(LP,901) NP,NEQ,MAXKB,TR,DR
901 FORMAT(T2,'INPUT DATA:',2X,I5/,T15,I5,2X,I6/,T15,F5.3,2X,F6.3//,
A T6,'# Line 1: NP'/,T6,'# Line 2: NEQ, MAXKB'/,
A T6,'# Line 3: TR, DR'/,40('-')///, T2,'OUTPUT:'/)
WRITE(LP,902)
902 FORMAT(T20,49('*'))
WRITE(LP,904)
904 FORMAT(T20,'*',T68,'*')
WRITE(LP,906) NP
906 FORMAT(T20,'*',6X,'MOLECULAR DYNAMICS FOR',I4,' LJ ATOMS'
A ,T68,'*')
WRITE(LP,904)
WRITE(LP,902)
WRITE(LP,904)
WRITE(LP,908) EPSI,SIGMA
908 FORMAT(T20,'*',2X,'EPSI/K = ',F7.3,T49,'SIGMA = ',F7.3,
A T68,'*')
WRITE(LP,910) TR,T
910 FORMAT(T20,'*',2X,'TR = ',F7.3,T49,'T = ',F7.3,
A T68,'*')
WRITE(LP,912) DR,VOL
912 FORMAT(T20,'*',2X,'DR = ',F7.3,T49,'VOL = ',F7.3,
A T68,'*')
WRITE(LP,914) CUBE,RMAX
914 FORMAT(T20,'*',2X,'CUBE = ',F7.3,T49,'RMAX = ',F7.3,
A T68,'*')
WRITE(LP,916) RC,RDEL
916 FORMAT(T20,'*',2X,'RC = ',F7.3,T49,'DELR = ',F7.3,
A T68,'*')
WRITE(LP,918) RLIST,RDMAX
918 FORMAT(T20,'*',2X,'RLIST = ',F7.3,T49,'RDMAX = ',F7.3,
A T68,'*')
WRITE(LP,920) DELTA,NY
920 FORMAT(T20,'*',2X,'DELTA = ',F7.3,T49,'NY = ',I3,
A T68,'*')
WRITE(LP,904)
WRITE(LP,922) TSTEP
922 FORMAT(T20,'*',2X,'TIME UNIT = ',F6.3,'D-12 SEC',T68,'*')
WRITE(LP,924) TT1
924 FORMAT(T20,'*',2X,'TIME STEP = ',1PD10.3,' SEC',T68, '*')
WRITE(LP,904)
WRITE(LP,926) DE
926 FORMAT(T20,'*',2X,'ENERGY CORRECTION = ',F7.3,T68,'*')
WRITE(LP,928) DP
928 FORMAT(T20,'*',2X,'PRESSURE CORRECTION = ',F7.3,T68,'*')
WRITE(LP,904)
WRITE(LP,902)
C
WRITE(LP,950) SEED
950 FORMAT(/////T20,'RANDOM NUMBER SEED = ',F5.3,/)
C
C PRINT IDENTIFICATION OF OUTPUT COLUMN HEADINGS
WRITE(LP,951)
951 FORMAT(/////T20,'OUTPUT COLUMN HEADING SYMBOLS'/)
WRITE(LP,952)
952 FORMAT(T24,'TIME ELAPSED TIME IN PICOSECONDS')
WRITE(LP,953)
953 FORMAT(T24,'PRES RUNNING AVERAGE PRESSURE, P*')
WRITE(LP,954)
954 FORMAT(T24,'P1 INSTANTANEOUS PRESSURE')
WRITE(LP,955)
955 FORMAT(T24,'ENRG RUNNING AVERAGE CONF INTERNAL ENERGY')
WRITE(LP,956)
956 FORMAT(T24,'E1 INSTANTANEOUS CONF INTERNAL ENERGY')
WRITE(LP,957)
957 FORMAT(T24,'DIST MEAN SQ DISPLACEMENT FROM TIME ZERO')
WRITE(LP,958)
958 FORMAT(T24,'TEMP RUNNING AVERAGE TEMPERATURE')
WRITE(LP,959)
959 FORMAT(T24,'NAB INSTANTANEOUS TOTAL NEIGHBORS IN LIST')
WRITE(LP,960)
960 FORMAT(T24,'ETOTAL INSTANTANEOUS TOTAL ENERGY')
C
WRITE ( LP,970 ) MAXKB,KBSKP,ISTEP
970 FORMAT ( //T20,'MAXKB =',I6,5X,'SKIP =',I5,5X,'SAVE EVERY',I3,
1 'TH STEP')
RETURN
END
BLOCK DATA
PARAMETER ( NPD = 1000, NP6 = NPD*6, NP9 = NPD*9 )
COMMON/DER/X2(NPD),Y2(NPD),Z2(NPD),X3(NPD),Y3(NPD),Z3(NPD)
A ,X4(NPD),Y4(NPD),Z4(NPD),X5(NPD),Y5(NPD),Z5(NPD)
COMMON/FOR/FX(NPD),FY(NPD),FZ(NPD)
COMMON/DISP/DAX(NPD),DAY(NPD),DAZ(NPD),X0L(NPD),Y0L(NPD),
A Z0L(NPD)
COMMON/PROP/IDIST(1000),SUME,SUMV,XSUM,ISUM
DATA X3,Y3,Z3,X4,Y4,Z4,X5,Y5,Z5/NP9*0.0/
DATA DAX,DAY,DAZ,FX,FY,FZ/NP6*0.0/
DATA IDIST/1000*0/
END