C ---------------------------------------------------------------- C PROGRAM ICE (FLEXIBLE MODEL) FEB.28,1995 LIJU. C 1. USE EWALD SUMMATION. C 2. ABLE TO CONTINUE FROM LAST CONFIGURATION C 3. WITH HYDROGEN ATOM VELOCITY AUTOCORRELATION FUNCTIONAL C 4. START WITH ZERO-DIPOLE SELF-CONSISTENT CONFIGURATION C 5. NEW AHARMONIC INTRAMOLECULAR POTENTIAL C 6. PARRINELLO & RAHMAN METHOD WITH SYMMETRIC H C ---------------------------------------------------------------- C USE THE FOLLOWING UNITS: C ENERGY UNIT = THE ELECTROSTATIC ENERGY BETWEEN TWO C ELECTRONS WITH 1 ANGSTROMS APART C = 2.3070795E-18 JOULE C MASS UNIT = 1 AMU = 1.66054E-27 KG C LENGTH UNIT = 1 ANGSTROM C TIME UNIT = 2.6828315E-15 SECOND = 2.6828315E-3 PS C ---------------------------------------------------------------- C HERE THE X,Y,Z,X1,Y1,Z1 ARE ACTUALLY SX,SY,SZ,SX1,SY1,SZ1 C THE RANGE OF X,Y,Z ARE FROM 0 TO CUBE IMPLICIT REAL(A-H,O-Z) PARAMETER (NPD = 648) PARAMETER (ID = 20000) PARAMETER (IV = 200) COMMON /EWALD/ ALPHA,ERR(ID),CLI(ID),DERR,A(9,9,9) COMMON /POS/ X(NPD),Y(NPD),Z(NPD),X1(NPD),Y1(NPD),Z1(NPD) DOUBLE PRECISION X1,Y1,Z1 COMMON /FORCE/ FX(NPD),FY(NPD),FZ(NPD),A12,A6 COMMON /PARM/ RDEL,IGOFR(3,1500),NM,NPO,NPH,NP,PMASS(NPD), 1 ITYPE(NPD),Q(NPD),VOLUME,CUBE,FROG COMMON /ENERGY/ TOTE,POTE,KINE,PLJ12,PLJ6,ELEC,VINTRA,VIRIAL COMMON /CONST/ QO,QH,QOO,QOH,QHH,MO,MH,BOND,HHDIST,ANGLE,PI,P2 COMMON /CONTROL/ KB,KSTART,SCALER,SEED,DELTA,BON,ANG,PM2 COMMON /ADDENDUM/ PMX,PMY,PMZ COMMON /PR/ H(3,3),V(3,3),B(3,3) DOUBLE PRECISION H REAL KINE,MO,MH REAL KINEM,KINESUM DIMENSION V2(IV,216),VX(IV,216),VY(IV,216),VZ(IV,216),CORR(5000), 1 X2(NPD),Y2(NPD),Z2(NPD),VMX(216),VMY(216),VMZ(216), 2 RX1(NPD),RY1(NPD),RZ1(NPD) DOUBLE PRECISION HH(3,3),U(3,3),U1(3,3),U2(3,3),UU2(3,3),P(3,3), 1 WMASS DIMENSION HSUM(3,3),BSUM(3,3),BBSUM(3,3,3,3) INTEGER SCALE DATA CORR(5000) /0./ DATA IGOFR(3,1500) /0/ DATA LP,LP_GR,LP_CORR,LP_POSI,LP_CONFIG,LP_BAK /25,26,27,28,29,30/ C SELECT 200 TIME ORIGINS, EACH WITH 10 STEPS APART C TOTAL CORRELATION LENGTH = 5000 STEPS = 2 ps KVSTEP = 2000/IV AVO=6.0221367E+23 BOLZ=1.380658E-23 AMU=1.6605E-24 WMASS = 200 C ONE TIME STEP = 0.0004ps DELTA = 0.149096206 NM = 216 NPO = 1 NPH = 2 NP = NM*(NPO+NPH) C A12 = 1895.294007 C A6 = -1.883540782 C QO = -0.8476 READ *,A12,A6,QO QH = -QO/2. QOO = QO*QO QOH = QO*QH QHH = QH*QH MO = 15.9994 MH = 1.00797 AMASS = MO+2.*MH HPART = MH/AMASS OPART = MO/AMASS BOND = 0.9572 ANGLE = 1.824218103 HHDIST= 2.*SIN(ANGLE/2.)*BOND PI = 3.141592654 P2 = 2.*PI C RANDOM NUMBER GENERATOR SEED SEED = 0.2573 C WATER OXYGEN AND HYDROGEN DO 1100 I = 1,NP-2,3 PMASS(I) = MO PMASS(I+1)=MH PMASS(I+2)=MH ITYPE(I) = 0 ITYPE(I+1)=1 ITYPE(I+2)=1 Q(I) = QO Q(I+1) = QH Q(I+2) = QH 1100 CONTINUE READ *,T,MAXKB,KSCALE,KSHAKE,KSTART READ *,LAST,K_CORR,K_POSI READ *,DENSITY READ *,P(1,1),P(2,2),P(3,3),P(1,2),P(1,3),P(2,3) P(2,1) = P(1,2) P(3,1) = P(1,3) P(3,2) = P(2,3) TOTALMASS = NM*AMASS VOLUME = TOTALMASS*1.66054/DENSITY CUBE = VOLUME**(1./3.) C SCALER IS THE IDEAL TOTAL KINETIC ENERGY OF CM SCALER = 1.5*BOLZ*T*NM/2.3070795E-18 IF (LAST.EQ.0) THEN DO 388 I=1,3 DO 388 J=1,3 HH(I,J) = 0. IF (I.EQ.J) HH(I,J) = 1. 388 CONTINUE CALL DIAMOND CALL INTVEL ELSE READ *,CUBE1,VOLUME1,T1 READ *,HH(1,1),HH(2,2),HH(3,3),HH(1,2),HH(1,3),HH(2,3) HH(2,1) = HH(1,2) HH(3,1) = HH(1,3) HH(3,2) = HH(2,3) CUBE = CUBE1*((VOLUME/VOLUME1)**(1./3.)) DO 377 I=1,NP READ *,X(I),Y(I),Z(I),X1(I),Y1(I),Z1(I) X(I) = X(I)*CUBE/CUBE1 Y(I) = Y(I)*CUBE/CUBE1 Z(I) = Z(I)*CUBE/CUBE1 X1(I) = X1(I)*SQRT(T/T1) Y1(I) = Y1(I)*SQRT(T/T1) Z1(I) = Z1(I)*SQRT(T/T1) 377 CONTINUE ENDIF FROG = P2/CUBE DO 833 I=1,3 DO 833 J=1,3 U(I,J) = 0. U1(I,J) = 0. U2(I,J) = 0. 833 CONTINUE RDEL = CUBE/2./500. OPEN (UNIT=LP, FORM='FORMATTED',FILE="temp.out") OPEN (UNIT=LP_CONFIG, FORM='FORMATTED',FILE="config") OPEN (UNIT=LP_BAK, FORM='FORMATTED',FILE="config.bak") OPEN (UNIT=LP_GR, FORM='FORMATTED',FILE="gr.out") IF (K_CORR.EQ.1) 1 OPEN (UNIT=LP_CORR, FORM='FORMATTED',FILE="corr.out") IF (K_POSI.EQ.1) 1 OPEN (UNIT=LP_POSI, FORM='FORMATTED',FILE="posi.out") C ------------------------------------------------------------- C KWRITE IS THE FREQUENCY OF WRITING ON TEMP.OUT C SCALE IS THE QUANTA OF TIME STEPS TO SCALE TEMPERATURE C KSHAKE IS THE MAXIMUM TIMES OF MORE FREQUENT SCALING C KLOVE IS THE FREQUENCY OF MORE SLOWLY SCALING IN UNIT OF SCALE C KWGR IS THE FREQUENCY OF WRITING ON GR.OUT C -------------------------------------------------------------- KWRITE = 20 SCALE = 20 KLOVE = 20 KSHAKE = KSHAKE/(SCALE*KLOVE)*KLOVE KWGR = 200000 KWBAK = 11677 C SET CONSTS AND TABLES FOR EWALD SUMMATION ALPHA = 5.714/CUBE DERR = ALPHA*CUBE/2.*1.05 DERR = DERR/ID C = DERR WW = 2./SQRT(PI) DO 500 I=1,ID ERR(I) = ERFC(C) CLI(I) = WW*EXP(-C*C) C = C+DERR 500 CONTINUE C ---------------------------------------------------------- C PRINT PARAMETERS. 5 WRITE(LP,10) 10 FORMAT(///) WRITE(LP,20) 20 FORMAT(7X,50('*')) WRITE(LP,30) 30 FORMAT(7X,'*',T57,'*') WRITE(LP,40) 40 FORMAT(7X,'*',3X,'MOLECULAR DYNAMICS SIMULATION USING SPC/F', A T57,'*') WRITE(LP,50) NM 50 FORMAT(7X,'*',3X,'FLEXIBLE MODEL FOR',I4, A ' WATER MOLECULES,',T57,'*') WRITE(LP,70) 70 FORMAT(7X,'*',3X,'WITH DIAMOND CUBIC STRUCTURE.', A T57,'*') WRITE(LP,30) WRITE(LP,20) WRITE(LP,30) WRITE(LP,81) P(1,1),P(1,2),P(1,3) 81 FORMAT(7X,'*',2X,'STRESS =',3(F11.2),T56,'*') WRITE(LP,82) P(1,2),P(2,2),P(2,3) 82 FORMAT(7X,'*',2X,'(BAR) ',3(F11.2),T56,'*') WRITE(LP,83) P(1,3),P(2,3),P(3,3) 83 FORMAT(7X,'*',2X,' ',3(F11.2),T56,'*') WRITE(LP,30) WRITE(LP,80) T,DENSITY 80 FORMAT(7X,'*',2X,'T(K) =',F8.3,T30,'DENSITY(g/cm) = ', 1 F7.4,T57,'*') WRITE(LP,90) MAXKB,DELTA*2.68279E-3 90 FORMAT(7X,'*',2X,'MAXKB =',I6,T30, 'TIME STEP(ps) = ', 1 F7.4,T57,'*') WRITE(LP,30) WRITE(LP,111) QO,QH 111 FORMAT(7X,'*',2X,'QO = ',F7.4,T30,'QH = ',F7.4,T57,'*') WRITE(LP,112) A12,A6 112 FORMAT(7X,'*',2X,'A12 = ',F11.6,T30,'A6 = ',F11.7,T57,'*') WRITE(LP,30) WRITE(LP,110) 7*7*7-1 110 FORMAT(7X,'*',2X,'EWALD SUMMATION OVER ',I3, 1 ' RECIPROCAL LATTICES',T57,'*') WRITE(LP,30) WRITE(LP,20) WRITE(LP,190) 190 FORMAT(///,4X,'KB',4X,'TIME',5X,'TEMP',5X,'POTE',6X,'PRES', 1 4X,'H11',4X,'H22',4X,'H33',4X,'H12',4X,'H13',4X, 2 'H23'/) DO 55 I=1,NP X2(I) = 0. Y2(I) = 0. Z2(I) = 0. 55 CONTINUE DO 534 I=1,3 DO 534 J=1,3 P(I,J) = P(I,J)/2.3070795E+7 534 CONTINUE TSUM = 0. VINSUM = 0. KINESUM = 0. POTESUM = 0. PRESUM = 0. BONSUM = 0. ANGSUM = 0. PM2SUM = 0. PMXSUM = 0. PMYSUM = 0. PMZSUM = 0. DO 930 I=1,3 DO 930 J=1,3 HSUM(I,J) = 0. BSUM(I,J) = 0. DO 930 K=1,3 DO 930 L=1,3 BBSUM(I,J,K,L) = 0. 930 CONTINUE C ----------------------------------------------------------- C ENTER MAIN LOOP OF SIMULATION DO 310 KB=1,MAXKB C ------------ PARRINELLO & RAHMAN FUNCTIONAL ---------------- DO 380 I=1,3 DO 380 J=1,3 U(I,J) = DELTA*(U1(I,J)+0.5*DELTA*U2(I,J)) IF (I.EQ.J) U(I,J)=U(I,J)+1. 380 CONTINUE DO 381 I=1,3 DO 381 J=1,3 H(I,J) = U(I,1)*HH(1,J)+U(I,2)*HH(2,J)+U(I,3)*HH(3,J) 381 CONTINUE V(1,1) = H(2,2)*H(3,3)-H(2,3)*H(2,3) V(2,2) = H(1,1)*H(3,3)-H(1,3)*H(1,3) V(3,3) = H(1,1)*H(2,2)-H(1,2)*H(1,2) V(1,2) = H(1,3)*H(2,3)-H(1,2)*H(3,3) V(1,3) = H(1,2)*H(2,3)-H(1,3)*H(2,2) V(2,3) = H(1,3)*H(1,2)-H(2,3)*H(1,1) V(2,1) = V(1,2) V(3,1) = V(1,3) V(3,2) = V(2,3) VOLUME = H(1,1)*V(1,1)+H(1,2)*V(1,2)+H(1,3)*V(1,3) DO 633 I=1,3 DO 633 J=1,3 HH(I,J) = H(I,J) V(I,J) = V(I,J)/VOLUME 633 CONTINUE VOLUME = VOLUME*CUBE*CUBE*CUBE DO 600 I = -4,4 DO 600 J = -4,4 DO 600 K = -4,4 QX = (I*V(1,1)+J*V(2,1)+K*V(3,1))*FROG QY = (I*V(1,2)+J*V(2,2)+K*V(3,2))*FROG QZ = (I*V(1,3)+J*V(2,3)+K*V(3,3))*FROG Q2 = QX*QX+QY*QY+QZ*QZ IF (Q2.EQ.0) GOTO 600 C1 = EXP(-Q2/ALPHA/ALPHA/4.) A(I+5,J+5,K+5) = 4.*PI*C1/Q2/VOLUME 600 CONTINUE C ----------- PARRINELLO FUNCTIOAL ENDS ---------------- DO 320 I=1,NP DXX = DELTA*(X1(I)+0.5*DELTA*X2(I)) DYY = DELTA*(Y1(I)+0.5*DELTA*Y2(I)) DZZ = DELTA*(Z1(I)+0.5*DELTA*Z2(I)) X(I) = X(I)+DXX Y(I) = Y(I)+DYY Z(I) = Z(I)+DZZ IF (X(I).GT.CUBE) X(I)=X(I)-CUBE IF (X(I).LT.0) X(I)=X(I)+CUBE IF (Y(I).GT.CUBE) Y(I)=Y(I)-CUBE IF (Y(I).LT.0) Y(I)=Y(I)+CUBE IF (Z(I).GT.CUBE) Z(I)=Z(I)-CUBE IF (Z(I).LT.0) Z(I)=Z(I)+CUBE 320 CONTINUE CALL INTERACT KINE = 0. DO 330 I=1,NP XX2 = FX(I)/PMASS(I) YY2 = FY(I)/PMASS(I) ZZ2 = FZ(I)/PMASS(I) X1(I) = X1(I)+DELTA*(X2(I)+XX2)/2. Y1(I) = Y1(I)+DELTA*(Y2(I)+YY2)/2. Z1(I) = Z1(I)+DELTA*(Z2(I)+ZZ2)/2. X2(I) = XX2 Y2(I) = YY2 Z2(I) = ZZ2 RX1(I) = H(1,1)*X1(I)+H(1,2)*Y1(I)+H(1,3)*Z1(I) RY1(I) = H(2,1)*X1(I)+H(2,2)*Y1(I)+H(2,3)*Z1(I) RZ1(I) = H(3,1)*X1(I)+H(3,2)*Y1(I)+H(3,3)*Z1(I) KINE = KINE+PMASS(I)*(RX1(I)**2+RY1(I)**2+RZ1(I)**2) B(1,1) = B(1,1)+PMASS(I)*RX1(I)*RX1(I) B(2,2) = B(2,2)+PMASS(I)*RY1(I)*RY1(I) B(3,3) = B(3,3)+PMASS(I)*RZ1(I)*RZ1(I) B(1,2) = B(1,2)+PMASS(I)*RX1(I)*RY1(I) B(1,3) = B(1,3)+PMASS(I)*RX1(I)*RZ1(I) B(2,3) = B(2,3)+PMASS(I)*RY1(I)*RZ1(I) 330 CONTINUE KINE = KINE/2. KINEM = 0. IP = 1 DO 66 I=1,NM VMX(I) = RX1(IP)*OPART+(RX1(IP+1)+RX1(IP+2))*HPART VMY(I) = RY1(IP)*OPART+(RY1(IP+1)+RY1(IP+2))*HPART VMZ(I) = RZ1(IP)*OPART+(RZ1(IP+1)+RZ1(IP+2))*HPART KINEM = KINEM+AMASS*(VMX(I)**2+VMY(I)**2+VMZ(I)**2) IP = IP+3 66 CONTINUE KINEM = KINEM/2. IF (MOD(KB,KWGR).EQ.0.OR.KB.EQ.MAXKB) THEN C ------------- PAIR DISTRIBUTION FUCTIONAL ------------- IF (KB.EQ.MAXKB) KWGR = MOD(MAXKB-1,KWGR)+1 C CALC PAIR DISTRIBUTION FUNCTIONS G(R)s DO 475 I = 1,500 DV = 4.*PI*RDEL**3*I*I OP = FLOAT(NM*NPO) HP = FLOAT(NM*NPH) GROO = 2.*VOLUME*FLOAT(IGOFR(1,I))/(DV*KWGR*OP*OP) GRHH = 2.*VOLUME*FLOAT(IGOFR(3,I))/(DV*KWGR*HP*HP) GROH = VOLUME*FLOAT(IGOFR(2,I))/(DV*KWGR*OP*HP) WRITE(LP_GR,4100) I,RDEL*I,GROO,GRHH,GROH 4100 FORMAT(2X,I4,4(2X,F7.4)) 475 CONTINUE IF (KB.NE.MAXKB) THEN WRITE(LP_GR,4200) 4200 FORMAT('B') DO 476 I = 1,1500 IGOFR(1,I) = 0 IGOFR(2,I) = 0 IGOFR(3,I) = 0 476 CONTINUE ENDIF ENDIF C ------------- PAIR DISTRIBUTION FUCTIONAL ------------- IF (.NOT.((LAST.EQ.0).AND.KB.LE.1000)) THEN C -------------- PARRINELLO FUNCTIONAL ------------------ B(2,1) = B(1,2) B(3,1) = B(1,3) B(3,2) = B(2,3) DO 433 I=1,3 DO 433 J=1,3 UU2(I,J) = (B(I,J)/VOLUME-P(I,J))/WMASS U1(I,J) = U1(I,J)+DELTA*(U2(I,J)+UU2(I,J))/2. U2(I,J) = UU2(I,J) 433 CONTINUE C -------------- PARRINELLO FUNCTIONAL ENDS -------------- ENDIF C PRESSURE WILL BE IN BARS & ENERGIES WILL BE IN KJOULE/MOL. PRES = (VIRIAL+2.*KINE)/3./VOLUME*2.3070795E+7 DO 731 I=1,3 DO 731 J=1,3 B(I,J) = B(I,J)/VOLUME*2.3070795E+7 731 CONTINUE TNOW = KINEM*T/SCALER KINE = KINE *2.3070795E-18*AVO/NM/1000. POTE = POTE *2.3070795E-18*AVO/NM/1000. VINTRA = VINTRA*2.3070795E-18*AVO/NM/1000. C ACCUMMULATE SUMS FOR PROPERTY AVERAGES. IF (KB.GT.KSTART) THEN BONSUM = BONSUM + BON ANGSUM = ANGSUM + ANG PRESUM = PRESUM + PRES TSUM = TSUM + TNOW POTESUM = POTESUM + POTE KINESUM = KINESUM + KINE VINSUM = VINSUM + VINTRA PMXSUM = PMXSUM + PMX PMYSUM = PMYSUM + PMY PMZSUM = PMZSUM + PMZ PM2SUM = PM2SUM + PM2 DO 499 I=1,3 DO 499 J=1,3 HSUM(I,J) = HSUM(I,J) + H(I,J) BSUM(I,J) = BSUM(I,J) + B(I,J) DO 499 K=1,3 DO 499 L=1,3 BBSUM(I,J,K,L) = BBSUM(I,J,K,L)+B(I,J)*B(K,L) 499 CONTINUE ENDIF IF(MOD(KB,KWRITE).EQ.0) THEN C -------------------------------------------------------- C PROPERTY CALCULATION & PRINT-OUT AT INTERVALS OF KWRITE C REAL TIME IS IN ps REALTIME = DELTA*KB*2.6828315E-3 44 WRITE(LP,270) KB,REALTIME,TNOW,POTE,PRES,H(1,1), 1 H(2,2),H(3,3),H(1,2),H(1,3),H(2,3) 270 FORMAT(I5,F8.3,2(F9.3),F10.1,6(F8.3)) ENDIF C -------------------------------------------------------- TALL = TALL+TNOW IF (KB.LE.KSCALE.AND.MOD(KB,SCALE).EQ.0) THEN TALL=TALL/SCALE FACTOR=SQRT(T/TALL)-1. DO 383 I=1,NP II = (I-1)/3+1 X1(I)=X1(I)+VMX(II)*FACTOR Y1(I)=Y1(I)+VMY(II)*FACTOR Z1(I)=Z1(I)+VMZ(II)*FACTOR 383 CONTINUE TALL=0. KSHAKE = KSHAKE-1 IF (KSHAKE.EQ.0) SCALE=SCALE*KLOVE ENDIF IF (K_CORR.EQ.1) THEN C --------- VELOCITY CORRELATION FUNCTIONAL -------------- C CALCULATE VELOCITY AUTOCORRELATION FUNCTION I = KB+7000-MAXKB KPT = I/KVSTEP IF (MOD(I,KVSTEP).EQ.0.AND.KPT.GE.1.AND.KPT.LE.IV) THEN IP = 2 DO 1300 IM = 1,NM VX(KPT,IM) = RX1(IP) VY(KPT,IM) = RY1(IP) VZ(KPT,IM) = RZ1(IP) V2(KPT,IM) = RX1(IP)*RX1(IP)+RY1(IP)*RY1(IP)+RZ1(IP)*RZ1(IP) IP = IP+3 1300 CONTINUE ENDIF IF (KPT.GT.IV) KPT = IV DO 1600 KP = 1,KPT J = I-KP*KVSTEP IF (J.GE.1.AND.J.LE.5000) THEN IP = 2 DO 1400 IM = 1,NM CORR(J) = CORR(J)+(RX1(IP)*VX(KP,IM)+RY1(IP)*VY(KP,IM)+ 1 RZ1(IP)*VZ(KP,IM))/V2(KP,IM) IP = IP+3 1400 CONTINUE ENDIF 1600 CONTINUE ENDIF C --------- VELOCITY CORRELATION FUNCTIONAL -------------- IF (MOD(KB,KWBAK).EQ.0) THEN WRITE (LP_BAK,221) CUBE,VOLUME,TNOW WRITE (LP_BAK,222) H(1,1),H(2,2),H(3,3),H(1,2),H(1,3),H(2,3) DO 284 I=1,NP WRITE (LP_BAK,222) X(I),Y(I),Z(I),X1(I),Y1(I),Z1(I) 284 CONTINUE WRITE (LP_BAK,4200) ENDIF print *,'pote =',pote print *,'tote =',pote+kine print *,'pres =',pres print *,'B =',(b(1,1)+b(2,2)+b(3,3))/3. print *,h(1,1),h(1,2),h(1,3) print *,h(2,1),h(2,2),h(2,3) print *,h(3,1),h(3,2),h(3,3) print *,' ' 310 CONTINUE C END OF MAIN LOOP C ----------------------------------------------------------- IF (K_CORR.EQ.1) THEN DO 1500 I = 1,5000 CORR(I) = CORR(I)/NM/IV WRITE(LP_CORR,113) CORR(I) 113 FORMAT(F9.6) 1500 CONTINUE ENDIF IF (K_POSI.EQ.1) THEN DO 1700 I = 1,NP-2,3 WRITE(LP_POSI,233) X(I),Y(I),Z(I),X(I+1),Y(I+1),Z(I+1), 1 X(I+2),Y(I+2),Z(I+2) 233 FORMAT(9(F11.5)) 1700 CONTINUE ENDIF TSUM = TSUM/(MAXKB-KSTART) ANGSUM = ANGSUM/(MAXKB-KSTART) BONSUM = BONSUM/(MAXKB-KSTART) PRESUM = PRESUM/(MAXKB-KSTART) VINSUM = VINSUM/(MAXKB-KSTART) KINESUM = KINESUM/(MAXKB-KSTART) POTESUM = POTESUM/(MAXKB-KSTART) PMXSUM = PMXSUM/(MAXKB-KSTART) PMYSUM = PMYSUM/(MAXKB-KSTART) PMZSUM = PMZSUM/(MAXKB-KSTART) PM2SUM = PM2SUM/(MAXKB-KSTART) DO 444 I=1,3 DO 444 J=1,3 BSUM(I,J) = BSUM(I,J)/(MAXKB-KSTART) DO 444 K=1,3 DO 444 L=1,3 BBSUM(I,J,K,L) = BBSUM(I,J,K,L)/(MAXKB-KSTART) 444 CONTINUE DIE = (PM2SUM-PMXSUM*PMXSUM-PMYSUM*PMYSUM-PMZSUM*PMZSUM)*QH*QH DIE = 1.+DIE*4.*PI*2.3070795E-18/(3*BOLZ*TSUM)/VOLUME WRITE (LP,501) KSTART 501 FORMAT(//,5X,' THE RUNNING AVERAGE BEGINS FROM STEP',I6) WRITE (LP,601) MAXKB 601 FORMAT(5X, ' TO STEP',I6,/) WRITE (LP,701) TSUM 701 FORMAT(5X,' THE AVERAGED TEMPERATURE IS ',F7.2,' KELVIN') WRITE (LP,801) BONSUM 801 FORMAT(5X,' THE AVERAGED BOND LENGTH IS ',F7.5,' ANGSTROM') WRITE (LP,901) ANGSUM*180/PI 901 FORMAT(5X,' THE AVERAGED BOND ANGLE IS ',F7.3,' DEGREE',/) WRITE (LP,1001) PRESUM 1001 FORMAT(5X,' THE AVERAGED PRESSURE IS',F10.2,' BAR') WRITE (LP,1101) KINESUM 1101 FORMAT(5X,' THE AVERAGED KINETIC ENERGY IS',F7.2,' KJOULE/MOL') WRITE (LP,1301) VINSUM 1301 FORMAT(5X,' THE AVERAGED VINTRA IS',F7.4,' KJOULE/MOL') WRITE (LP,1201) POTESUM 1201 FORMAT(5X,' THE AVERAGED POTENTIAL ENERGY IS',F7.2,' KJOULE/MOL') WRITE (LP,1501) DIE 1501 FORMAT(5X,' THE STATIC DIELECTRIC CONSTANT IS',F7.2,/) WRITE (LP,1800) 1800 FORMAT(5X,' THE AVERAGED H MATRIX IS:',/) WRITE (LP,1801) HSUM(1,1),HSUM(1,2),HSUM(1,3) WRITE (LP,1801) HSUM(2,1),HSUM(2,2),HSUM(2,3) WRITE (LP,1801) HSUM(3,1),HSUM(3,2),HSUM(3,3) 1801 FORMAT(7X,'|',3(2X,F10.6,2X),' |') WRITE (LP,1401) 1.5*BOLZ*TSUM*AVO/10./KINESUM 1401 FORMAT(/,5X,' THE CM MOTION ACCOUNT FOR ',F7.4,' PERCENT OF THE', 1 /,5X,' TOTAL KINETIC ENERGY.') 221 FORMAT(3(1X,E13.7)) WRITE (LP_CONFIG,221) CUBE,VOLUME,TNOW WRITE (LP_CONFIG,222) H(1,1),H(2,2),H(3,3),H(1,2),H(1,3),H(2,3) 222 FORMAT(6(1X,E13.7)) DO 244 I=1,NP WRITE (LP_CONFIG,222) X(I),Y(I),Z(I),X1(I),Y1(I),Z1(I) 244 CONTINUE STOP END C ------------------- PROGRAM END HERE ------------------- C -------------------------------------------------------- SUBROUTINE RANDOM IMPLICIT REAL(A-H,O-Z) COMMON /CONTROL/ KB,KSTART,SCALER,SEED,DELTA,BON,ANG,PM2 DATA K,J,M,RM / 5701,3612,566927,566927.0 / IX=INT(SEED*RM) KR = J*IX + K IRAND = MOD(KR,M) SEED = (FLOAT(IRAND) + 0.5) / RM RETURN END C --------------------------------------------------------- C --------------------------------------------------------- SUBROUTINE INTVEL C ASSIGN INITIAL VELOCITIES TO MOLECULES IMPLICIT REAL(A-H,O-Z) PARAMETER (NPD = 648) COMMON /POS/ X(NPD),Y(NPD),Z(NPD),X1(NPD),Y1(NPD),Z1(NPD) DOUBLE PRECISION X1,Y1,Z1 COMMON /FORCE/ FX(NPD),FY(NPD),FZ(NPD),A12,A6 COMMON /PARM/ RDEL,IGOFR(3,1500),NM,NPO,NPH,NP,PMASS(NPD), 1 ITYPE(NPD),Q(NPD),VOLUME,CUBE,FROG COMMON /CONST/ QO,QH,QOO,QOH,QHH,MO,MH,BOND,HHDIST,ANGLE,PI,P2 COMMON /CONTROL/ KB,KSTART,SCALER,SEED,DELTA,BON,ANG,PM2 REAL MO,MH PX=0. PY=0. PZ=0. DO 360 I=1,NP IF (MOD(I,3).EQ.1) THEN CALL RANDOM SEED1 = SEED CALL RANDOM SEED2 = SEED CALL RANDOM SEED3 = SEED ENDIF X1(I) = SEED1 Y1(I) = SEED2 Z1(I) = SEED3 PX = PX + PMASS(I)*X1(I) PY = PY + PMASS(I)*Y1(I) PZ = PZ + PMASS(I)*Z1(I) 360 CONTINUE C SCALE VELOCITIES SO THAT TOTAL MOMEMTUM=0. SUMKINE = 0. C = 1./FLOAT(NM)/(NPO*MO+NPH*MH) DO 370 I=1,NP X1(I)=X1(I)-PX*C Y1(I)=Y1(I)-PY*C Z1(I)=Z1(I)-PZ*C SUMKINE = SUMKINE+PMASS(I)*(X1(I)*X1(I)+Y1(I)*Y1(I)+Z1(I)*Z1(I)) 370 CONTINUE SUMKINE = SUMKINE/2. C SCALE VELOCITIES TO DESIRED TEMPERATURE FACTOR=SQRT(SCALER/SUMKINE) DO 380 I=1,NP X1(I)=X1(I)*FACTOR Y1(I)=Y1(I)*FACTOR Z1(I)=Z1(I)*FACTOR 380 CONTINUE RETURN END C -------------------------------------------------------------- C -------------------------------------------------------------- SUBROUTINE INTERACT IMPLICIT REAL(A-H,O-Z) PARAMETER (NPD = 648) PARAMETER (ID = 20000) COMMON /EWALD/ ALPHA,ERR(ID),CLI(ID),DERR,A(9,9,9) COMMON /POS/ X(NPD),Y(NPD),Z(NPD),X1(NPD),Y1(NPD),Z1(NPD) DOUBLE PRECISION X1,Y1,Z1 COMMON /FORCE/ FX(NPD),FY(NPD),FZ(NPD),A12,A6 COMMON /PARM/ RDEL,IGOFR(3,1500),NM,NPO,NPH,NP,PMASS(NPD), 1 ITYPE(NPD),Q(NPD),VOLUME,CUBE,FROG COMMON /ENERGY/ TOTE,POTE,KINE,PLJ12,PLJ6,ELEC,VINTRA,VIRIAL COMMON /CONST/ QO,QH,QOO,QOH,QHH,MO,MH,BOND,HHDIST,ANGLE,PI,P2 COMMON /CONTROL/ KB,KSTART,SCALER,SEED,DELTA,BON,ANG,PM2 COMMON /ADDENDUM/ PMX,PMY,PMZ COMMON /PR/ H(3,3),V(3,3),B(3,3) DOUBLE PRECISION H REAL KINE,MO,MH DIMENSION SINN(NPD,9,9,9),COSN(NPD,9,9,9) DIMENSION XR(3,3),YR(3,3),ZR(3,3),DX(3,3),DY(3,3),DZ(3,3),DR(3,3) BON = 0. ANG = 0. PMX = 0. PMY = 0. PMZ = 0. ELEC = 0. PLJ6 = 0. PLJ12 = 0. VINTRA = 0. VIRIAL = 0. W = 2.*ALPHA/SQRT(PI) RCUT = CUBE/2. TUCR = -RCUT DO 333 I=1,3 DO 333 J=1,3 B(I,J) = 0. 333 CONTINUE C RECIPROCAL SPACE SUMMATION DO 805 I = 1,NP C THE INERT V2 TERM ELEC = ELEC - Q(I)*Q(I)*W/2. FX(I) = 0. FY(I) = 0. FZ(I) = 0. 805 CONTINUE DO 1000 I=1,NP SX = X(I)*FROG SY = Y(I)*FROG SZ = Z(I)*FROG SINX = SIN(SX) COSX = COS(SX) SINY = SIN(SY) COSY = COS(SY) SINZ = SIN(SZ) COSZ = COS(SZ) SINN(I,1,1,1) = SIN(-4.*SX-4.*SY-4.*SZ) COSN(I,1,1,1) = COS(-4.*SX-4.*SY-4.*SZ) DO 1100 J=2,9 SINN(I,J,1,1) = SINN(I,J-1,1,1)*COSX+COSN(I,J-1,1,1)*SINX COSN(I,J,1,1) = COSN(I,J-1,1,1)*COSX-SINN(I,J-1,1,1)*SINX 1100 CONTINUE DO 1200 J=1,9 DO 1200 K=2,9 SINN(I,J,K,1) = SINN(I,J,K-1,1)*COSY+COSN(I,J,K-1,1)*SINY COSN(I,J,K,1) = COSN(I,J,K-1,1)*COSY-SINN(I,J,K-1,1)*SINY 1200 CONTINUE DO 1300 J=1,9 DO 1300 K=1,9 DO 1300 L=2,9 SINN(I,J,K,L) = SINN(I,J,K,L-1)*COSZ+COSN(I,J,K,L-1)*SINZ COSN(I,J,K,L) = COSN(I,J,K,L-1)*COSZ-SINN(I,J,K,L-1)*SINZ 1300 CONTINUE 1000 CONTINUE DO 1400 J=1,9 DO 1400 K=1,9 DO 1400 L=1,9 XX = (J-5)*FROG YY = (K-5)*FROG ZZ = (L-5)*FROG QX = V(1,1)*XX+V(1,2)*YY+V(1,3)*ZZ QY = V(2,1)*XX+V(2,2)*YY+V(2,3)*ZZ QZ = V(3,1)*XX+V(3,2)*YY+V(3,3)*ZZ Q2 = QX*QX+QY*QY+QZ*QZ IF (Q2.EQ.0.) GOTO 1400 QK = 2./Q2+1./2./ALPHA/ALPHA C (QX,QY,QZ)+Q2+QK ARE THE TRUE FORCE DIRECTION XX = V(1,1)*QX+V(1,2)*QY+V(1,3)*QZ YY = V(2,1)*QX+V(2,2)*QY+V(2,3)*QZ ZZ = V(3,1)*QX+V(3,2)*QY+V(3,3)*QZ C (XX,YY,ZZ) ARE THE FORCE DIRECTION IN S SPACE. COSSUM = 0. SINSUM = 0. DO 1500 I=1,NP COSSUM = COSSUM + Q(I)*COSN(I,J,K,L) SINSUM = SINSUM + Q(I)*SINN(I,J,K,L) 1500 CONTINUE DO 1600 I=1,NP CC = A(J,K,L)*Q(I)*(COSSUM*SINN(I,J,K,L)-SINSUM*COSN(I,J,K,L)) FX(I) = FX(I)+XX*CC FY(I) = FY(I)+YY*CC FZ(I) = FZ(I)+ZZ*CC 1600 CONTINUE CC = A(J,K,L)*(COSSUM*COSSUM+SINSUM*SINSUM)/2. B(1,1) = B(1,1)+(1-QK*QX*QX)*CC B(2,2) = B(2,2)+(1-QK*QY*QY)*CC B(3,3) = B(3,3)+(1-QK*QZ*QZ)*CC B(1,2) = B(1,2)-QK*QX*QY*CC B(1,3) = B(1,3)-QK*QX*QZ*CC B(2,3) = B(2,3)-QK*QY*QZ*CC ELEC = ELEC + CC 1400 CONTINUE DO 522 I=1,NP IMI = (I-1)/3 II = I-IMI*3 DO 5520 J=I+1,NP IMJ = (J-1)/3 IJ = J-IMJ*3 XX = X(I)-X(J) YY = Y(I)-Y(J) ZZ = Z(I)-Z(J) IF (XX.GT.RCUT) XX = XX - CUBE IF (XX.LT.TUCR) XX = XX + CUBE IF (YY.GT.RCUT) YY = YY - CUBE IF (YY.LT.TUCR) YY = YY + CUBE IF (ZZ.GT.RCUT) ZZ = ZZ - CUBE IF (ZZ.LT.TUCR) ZZ = ZZ + CUBE RX = H(1,1)*XX+H(1,2)*YY+H(1,3)*ZZ RY = H(2,1)*XX+H(2,2)*YY+H(2,3)*ZZ RZ = H(3,1)*XX+H(3,2)*YY+H(3,3)*ZZ R2 = RX*RX+RY*RY+RZ*RZ R1 = SQRT(R2) IF (R1.GT.RCUT) GOTO 5520 R3 = R2*R1 IJT = ITYPE(I)+ITYPE(J)+1 IF (IMI.NE.IMJ) THEN IR = R1/RDEL+1 IGOFR(IJT,IR) = IGOFR(IJT,IR) + 1 ENDIF C CALCULATE EWALD FACTOR IN REAL SPACE SUMMATION PRODUCT = ALPHA*R1 I0 = PRODUCT/DERR CC = I0*DERR ARGIN = PRODUCT-CC GACTOR = ERR(I0)-CLI(I0)*ARGIN*(1-CC*ARGIN) C IF THEY ARE THE SAME MOLECULE IF (IMI.EQ.IMJ) THEN GACTOR = GACTOR-1. XR(II,IJ) = XX YR(II,IJ) = YY ZR(II,IJ) = ZZ DX(II,IJ) = RX DY(II,IJ) = RY DZ(II,IJ) = RZ DR(II,IJ) = R1 C XR,YR,ZR ARE S SPACE QUANTITIES, DR ARE REAL SPACE QUANTITY. ENDIF ARGIN = PRODUCT*PRODUCT-CC*CC E1 = 1.-ARGIN*(1-ARGIN*0.5) FACTOR = GACTOR + E1*PRODUCT*CLI(I0) IF ( IJT .EQ. 1 ) THEN R6 = R3*R3 R12 = R6*R6 E1 = A12/R12 E2 = A6/R6 PLJ12 = PLJ12 + E1 PLJ6 = PLJ6 + E2 ELEC = ELEC + GACTOR*QOO/R1 FF = 12.*E1 + 6.*E2 FF = FF/R2 + QOO*FACTOR/R3 ELSEIF ( IJT .EQ. 2 ) THEN IF (R1.LT.0.25) GOTO 5520 ELEC = ELEC + GACTOR*QOH/R1 FF = FACTOR*QOH/R3 ELSEIF ( IJT .EQ. 3 ) THEN ELEC = ELEC + GACTOR*QHH/R1 FF = FACTOR*QHH/R3 ENDIF FFX= FF*XX FFY= FF*YY FFZ= FF*ZZ FX(I) = FX(I) + FFX FY(I) = FY(I) + FFY FZ(I) = FZ(I) + FFZ FX(J) = FX(J) - FFX FY(J) = FY(J) - FFY FZ(J) = FZ(J) - FFZ C THAT'S FORCE IN S SPACE C WE DO STRESS SUMMATION IN REAL SPACE B(1,1) = B(1,1)+RX*RX*FF B(2,2) = B(2,2)+RY*RY*FF B(3,3) = B(3,3)+RZ*RZ*FF B(1,2) = B(1,2)+RX*RY*FF B(1,3) = B(1,3)+RX*RZ*FF B(2,3) = B(2,3)+RY*RZ*FF 5520 CONTINUE IF (II.EQ.3) THEN PMX = PMX+DX(1,2)+DX(1,3) PMY = PMY+DY(1,2)+DY(1,3) PMZ = PMZ+DZ(1,2)+DZ(1,3) BON = BON+DR(1,2)+DR(1,3) ANG = ANG+ACOS((DR(1,2)**2+DR(1,3)**2-DR(2,3)**2)/ 1 2./DR(1,2)/DR(1,3)) DR1 = DR(1,2)-BOND DR2 = DR(1,3)-BOND DR3 = DR(2,3)-HHDIST F1 = 3.959608852*DR1+0.2514577891*DR2-0.5348000472*DR3 G1 = -13.52143988*DR1*DR1-0.768247786*DR1*DR2 1 +1.233829205*DR1*DR3-0.3841238925*DR2*DR2 2 +0.944571334*DR3*DR2-0.6852521196*DR3*DR3 F2 = 3.959608852*DR2+0.2514577891*DR1-0.5348000472*DR3 G2 = -13.52143988*DR2*DR2-0.768247786*DR2*DR1 1 +1.233829205*DR2*DR3-0.3841238925*DR1*DR1 2 +0.944571334*DR3*DR1-0.6852521196*DR3*DR3 F3 = -0.5348000472*(DR1+DR2)+0.8804540624*DR3 G3 = 0.6169146026*(DR1*DR1+DR2*DR2)+0.6267843126*DR3*DR3+ 1 0.944571334*DR1*DR2-1.370504239*(DR1+DR2)*DR3 VINTRA = VINTRA + (F1/2.+G1/3.)*DR1+(F2/2.+G2/3.)*DR2 1 +(F3/2.+G3/3.)*DR3 F1 = (F1+G1)/DR(1,2) F2 = (F2+G2)/DR(1,3) F3 = (F3+G3)/DR(2,3) FX(I-2) = FX(I-2)-F1*XR(1,2)-F2*XR(1,3) FY(I-2) = FY(I-2)-F1*YR(1,2)-F2*YR(1,3) FZ(I-2) = FZ(I-2)-F1*ZR(1,2)-F2*ZR(1,3) FX(I-1) = FX(I-1)+F1*XR(1,2)-F3*XR(2,3) FY(I-1) = FY(I-1)+F1*YR(1,2)-F3*YR(2,3) FZ(I-1) = FZ(I-1)+F1*ZR(1,2)-F3*ZR(2,3) FX(I) = FX(I)+F2*XR(1,3)+F3*XR(2,3) FY(I) = FY(I)+F2*YR(1,3)+F3*YR(2,3) FZ(I) = FZ(I)+F2*ZR(1,3)+F3*ZR(2,3) C WE DO STRESS SUMMATION IN REAL SPACE VIRIAL = VIRIAL-F1*DR(1,2)*DR(1,2)- 1 F2*DR(1,3)*DR(1,3)-F3*DR(2,3)*DR(2,3) B(1,1) = B(1,1)-F1*DX(1,2)*DX(1,2)- 1 F2*DX(1,3)*DX(1,3)-F3*DX(2,3)*DX(2,3) B(2,2) = B(2,2)-F1*DY(1,2)*DY(1,2)- 1 F2*DY(1,3)*DY(1,3)-F3*DY(2,3)*DY(2,3) B(3,3) = B(3,3)-F1*DZ(1,2)*DZ(1,2)- 1 F2*DZ(1,3)*DZ(1,3)-F3*DZ(2,3)*DZ(2,3) B(1,2) = B(1,2)-F1*DX(1,2)*DY(1,2)- 1 F2*DX(1,3)*DY(1,3)-F3*DX(2,3)*DY(2,3) B(1,3) = B(1,3)-F1*DX(1,2)*DZ(1,2)- 1 F2*DX(1,3)*DZ(1,3)-F3*DX(2,3)*DZ(2,3) B(2,3) = B(2,3)-F1*DY(1,2)*DZ(1,2)- 1 F2*DY(1,3)*DZ(1,3)-F3*DY(2,3)*DZ(2,3) ENDIF 522 CONTINUE C CALCULATE TOTAL POTENTIAL ENERGIES BON = BON/2./NM ANG = ANG/NM PM2 = PMX*PMX+PMY*PMY+PMZ*PMZ C THE CUTOFF REMEDY PLJ6 = PLJ6+A6*NM/VOLUME*PI*2./3./RCUT/RCUT/RCUT*NM POTE = PLJ12+PLJ6+ELEC+VINTRA VIRIAL = VIRIAL+12.*PLJ12+6.*PLJ6+ELEC RETURN END C ---------------------------------------------------------------- C ---------------------------------------------------------------- SUBROUTINE DIAMOND IMPLICIT REAL(A-H,O-Z) PARAMETER (NPD = 648) COMMON /POS/ X(NPD),Y(NPD),Z(NPD),X1(NPD),Y1(NPD),Z1(NPD) DOUBLE PRECISION X1,Y1,Z1 COMMON /PARM/ RDEL,IGOFR(3,1500),NM,NPO,NPH,NP,PMASS(NPD), 1 ITYPE(NPD),Q(NPD),VOLUME,CUBE,FROG COMMON /CONST/ QO,QH,QOO,QOH,QHH,MO,MH,BOND,HHDIST,ANGLE,PI,P2 DIMENSION SX(24),SY(24),SZ(24) DIMENSION VECTOR(4,3),MTYPE(16) DATA VECTOR/ 1,1,-1,-1, Y 1,-1,-1,1, Z 1,-1,1,-1/ DATA MTYPE /4,2,4,2,3,1,3,1,-2,-1,-2,-1,-3,-4,-3,-4/ N = NINT((NM/8)**(1./3.)) A = CUBE/N SX(1) = 0. SY(1) = 0. SZ(1) = 0. SX(7) = 1./2. SY(7) = 1./2. SZ(7) = 0. SX(10) = 0. SY(10) = 1./2. SZ(10) = 1./2. SX(4) = 1./2. SY(4) = 0. SZ(4) = 1./2. DO 100 I = 13,22,3 SX(I) = SX(I-12)+0.25 SY(I) = SY(I-12)+0.25 SZ(I) = SZ(I-12)+0.25 100 CONTINUE IO = -2 IB = 1 DO 500 I = 1,24 IF (MOD(I,3).EQ.1) THEN SX(I) = SX(I)*A SY(I) = SY(I)*A SZ(I) = SZ(I)*A IO = IO+3 ELSE IC = MTYPE(IB) IB = IB+1 IF (IC.LT.0) THEN IC = -IC SX(I) = SX(IO)-VECTOR(IC,1)/SQRT(3.)*BOND SY(I) = SY(IO)-VECTOR(IC,2)/SQRT(3.)*BOND SZ(I) = SZ(IO)-VECTOR(IC,3)/SQRT(3.)*BOND ELSE SX(I) = SX(IO)+VECTOR(IC,1)/SQRT(3.)*BOND SY(I) = SY(IO)+VECTOR(IC,2)/SQRT(3.)*BOND SZ(I) = SZ(IO)+VECTOR(IC,3)/SQRT(3.)*BOND ENDIF ENDIF 500 CONTINUE M = 0 DO 200 I = 0,N-1 DO 200 J = 0,N-1 DO 200 K = 0,N-1 DO 300 L = 1,24 X(L+M) = SX(L)+I*A Y(L+M) = SY(L)+J*A Z(L+M) = SZ(L)+K*A 300 CONTINUE M = M + 24 200 CONTINUE DO 700 I= 1,NP IF(X(I).LT.0) X(I) = X(I)+CUBE IF(X(I).GT.CUBE) X(I) = X(I)-CUBE IF(Y(I).LT.0) Y(I) = Y(I)+CUBE IF(Y(I).GT.CUBE) Y(I) = Y(I)-CUBE IF(Z(I).LT.0) Z(I) = Z(I)+CUBE IF(Z(I).GT.CUBE) Z(I) = Z(I)-CUBE 700 CONTINUE RETURN END C ---------------------------------------------------