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