C ---------------------------------------------------------------- C PROGRAM ICE. JAN.19,1995 LIJU. C 1. USE EWARD SUMMATION. C 2. USE CONSTRAINT DYNAMICS "RATTLE" IN VELOCITY C VERLET ALGORITHM. C 3. USE EXTENDED SPC MODEL IN J.PHYS.CHEM(91,6269,1987) C 4. USE PARRINELLO & RAHMAN METHOD (SYMMETRIC H VERSION) C 5. ABLE TO CONTINUE FROM LAST CONFIGURATION C 6. WITH HYDROGEN ATOM VELOCITY AUTOCORRELATION FUNCTIONAL 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 ---------------------------------------------------------------- IMPLICIT REAL(A-H,O-Z) PARAMETER (NPD = 648) PARAMETER (ID = 20000) PARAMETER (IV = 200) COMMON /EWALD/ ALPHA,ERR(ID),CLI(ID),DERR,A(7,7,7) COMMON /POS/ X(NPD),Y(NPD),Z(NPD),SX(NPD),SY(NPD),SZ(NPD), 1 X1(NPD),Y1(NPD),Z1(NPD),SX1(NPD),SY1(NPD),SZ1(NPD) DOUBLE PRECISION X1(NPD),Y1(NPD),Z1(NPD) COMMON /FORCE/ FX(NPD),FY(NPD),FZ(NPD) COMMON /PARM/ RDEL,IGOFR(3,1500),NM,NPO,NPH,NP,PMASS(NPD), 1 ITYPE(NPD),Q(NPD) COMMON /ENERGY/ TOTE,POTE,KINE,PLJ12,PLJ6,ELEC,VIRIAL COMMON /CONST/ QO,QH,QOO,QOH,QHH,MO,MH,BOND,HHDIST,ANGLE,PI,P2 COMMON /CONTROL/ KB,KSTART,SCALER,SEED,DELTA,RX,RY,RZ,WMASS COMMON /H/ H(3,3),HH(3,3),H1(3,3),H2(3,3),V(3,3),P(3,3), 1 B(3,3),CUBE,VOLUME REAL KINE,MO,MH REAL POTESUM,KINESUM,H11,H22,H33,H12,H13,H23 DIMENSION V2(IV,216),VX(IV,216),VY(IV,216),VZ(IV,216),CORR(5000) INTEGER SCALE C SELECT 200 TIME ORIGINS, EACH WITH 10 STEPS APART C TOTAL CORRELATION LENGTH = 5000 STEPS = 5 ps KVSTEP = 2000/IV AVO=6.0221367E+23 BOLZ=1.380658E-23 AMU=1.6605E-24 C ONE TIME STEP = 0.001ps DELTA = 0.372740516 NM = 216 NPO = 1 NPH = 2 NP = NM*(NPO+NPH) C DEFINE CONSTS QO = -0.8476 QH = 0.4238 QOO = QO * QO QOH = QO * QH QHH = QH * QH MO = 15.9994 MH = 1.00797 AMASS = MO+2*MH BOND = 1.0 ANGLE = 1.910633236 HHDIST= 1.632993162 PI = 3.141592654 P2 = 2.*PI WMASS = 500000 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 *,P(1,1),P(2,2),P(3,3),P(1,2),P(1,3),P(2,3) READ *,LAST,K_CORR,K_POSI C SCALER IS THE IDEAL TOTAL KINETIC ENERGY IN REDUCED UNITS SCALER = 3.*BOLZ*T*NM/2.3070795E-18 IF (LAST.EQ.0) THEN READ *,DENSITY TOTALMASS = NM*AMASS VOLUME = TOTALMASS*1.66054/DENSITY CUBE = VOLUME**(1./3.) H(1,1) = CUBE H(2,2) = CUBE H(3,3) = CUBE V(1,1) = P2/CUBE V(2,2) = P2/CUBE V(3,3) = P2/CUBE CALL DIAMOND CALL INTVEL ELSE READ *,VOLUME CUBE = VOLUME**(1./3.) READ *,H(1,1),H(2,2),H(3,3),H(1,2),H(1,3),H(2,3) READ *,V(1,1),V(2,2),V(3,3),V(1,2),V(1,3),V(2,3) READ *,H1(1,1),H1(2,2),H1(3,3),H1(1,2),H1(1,3),H1(2,3) READ *,H2(1,1),H2(2,2),H2(3,3),H2(1,2),H2(1,3),H2(2,3) ENDIF DO 733 I=1,3 DO 733 J=I,3 HH(I,J) = H(I,J)/P2 H(J,I) = H(I,J) V(J,I) = V(I,J) H1(J,I) = H1(I,J) H2(J,I) = H2(I,J) HH(J,I) = HH(I,J) 733 CONTINUE IF (LAST.NE.0) THEN DO 377 I=1,NP READ *,X(I),Y(I),Z(I),X1(I),Y1(I),Z1(I) SX(I) = V(1,1)*X(I)+V(1,2)*Y(I)+V(1,3)*Z(I) SY(I) = V(2,1)*X(I)+V(2,2)*Y(I)+V(2,3)*Z(I) SZ(I) = V(3,1)*X(I)+V(3,2)*Y(I)+V(3,3)*Z(I) 377 CONTINUE ENDIF RDEL = CUBE/2./500. OPEN (UNIT=25, FORM='FORMATTED',FILE="temp.out") OPEN (UNIT=32, FORM='FORMATTED',FILE="config") OPEN (UNIT=26, FORM='FORMATTED',FILE="gr.out") C OPEN (UNIT=27, FORM='FORMATTED',FILE="gr01.out") IF (K_CORR.EQ.1) OPEN (UNIT=30, FORM='FORMATTED',FILE="corr.out") IF (K_POSI.EQ.1) OPEN (UNIT=31, FORM='FORMATTED',FILE="posi.out") LP=25 LP2=26 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 GR0*.OUT C -------------------------------------------------------------- KWRITE = 20 SCALE = 20 KLOVE = 10 KSHAKE = KSHAKE/SCALE/KLOVE*KLOVE KWGR = 10000 C SET CONSTS AND TABLES FOR EWALD SUMMATION ALPHA = 5.6/CUBE DERR = ALPHA*CUBE/2.*1.5 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 DO 600 I = -3,3 DO 600 J = -3,3 DO 600 K = -3,3 QX = I*V(1,1)+J*V(2,1)+K*V(3,1) QY = I*V(1,2)+J*V(2,2)+K*V(3,2) QZ = I*V(1,3)+J*V(2,3)+K*V(3,3) Q2 = QX*QX+QY*QY+QZ*QZ IF (Q2.EQ.0) GOTO 600 C1 = EXP(-Q2/ALPHA/ALPHA/4.) A(I+4,J+4,K+4) = 4.*PI*C1/Q2/VOLUME 600 CONTINUE C ---------------------------------------------------------- C PRINT PARAMETERS. 5 WRITE(LP,10) 10 FORMAT(///) WRITE(LP,20) 20 FORMAT(7X,49('*')) WRITE(LP,30) 30 FORMAT(7X,'*',T56,'*') WRITE(LP,40) 40 FORMAT(7X,'*',2X,'MOLECULAR DYNAMICS SIMULATION USING ', A T56,'*') WRITE(LP,50) NM 50 FORMAT(7X,'*',2X,'CONSTRAINT DYNAMICS "RATTLE" FOR ',I4, A ' WATER',T56,'*') WRITE(LP,60) 60 FORMAT(7X,'*',2X,'MOLECULES.',T56,'*') WRITE(LP,30) WRITE(LP,70) 70 FORMAT(7X,'*',2X,'THIS VERSION IS DIAMOND CUBIC ICE UNDER', A T56,'*') WRITE(LP,99) 99 FORMAT(7X,'*',2X,'CONSTANT STRESS',T56,'*') 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 80 FORMAT(7X,'*',2X,'T(K) =',F9.3,T56,'*') WRITE(LP,90) MAXKB,KWRITE 90 FORMAT(7X,'*',2X,'MAXKB =',I6,T30, 'KWRITE = ', 1 I3,T56,'*') WRITE(LP,101) KSCALE,DELTA*2.68279E-3 101 FORMAT(7X,'*',2X,'KSCALE =',I6,T30, 'TIME STEP(ps) = ', 1 F6.4,T56,'*') WRITE(LP,30) WRITE(LP,110) 7*7*7-1 110 FORMAT(7X,'*',2X,'EWARD SUMMATION OVER ',I3, 1 ' RECIPROCAL',T56,'*') WRITE(LP,115) 115 FORMAT(7X,'*',2X,'LATTICES',T56,'*') WRITE(LP,30) WRITE(LP,20) WRITE(LP,190) 190 FORMAT(///,4X,'KB',4X,'TIME',5X,'TEMP',5X,'POTE',5X,'PRES', 1 4X,'H11',4X,'H22',4X,'H33',4X,'H12',4X,'H13',4X, 2 'H23'/) CALL INTERACT C ----------------------------------------------------------- C ENTER MAIN LOOP OF SIMULATION DO 310 KB=1,MAXKB CALL RATTLE 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,1500 DV = 4.*PI*RDEL**3*I*I OPART = FLOAT(NM*NPO) HPART = FLOAT(NM*NPH) GROO = 2.*VOLUME*FLOAT(IGOFR(1,I))/(DV*KWGR*OPART*OPART) GRHH = 2.*VOLUME*FLOAT(IGOFR(3,I))/(DV*KWGR*HPART*HPART) GROH = VOLUME*FLOAT(IGOFR(2,I))/(DV*KWGR*OPART*HPART) WRITE(LP2,4100) I,RDEL*I,GROO,GRHH,GROH 4100 FORMAT(2X,I4,4(2X,F7.4)) 475 CONTINUE LP2 = LP2+1 DO 476 I = 1,1500 IGOFR(1,I) = 0 IGOFR(2,I) = 0 IGOFR(3,I) = 0 476 CONTINUE ENDIF C ------------- PAIR DISTRIBUTION FUCTIONAL ------------- C STRESS WILL BE IN BARS & ENERGIES WILL BE IN KJOULE/MOL. WW = 0. PRES = (VIRIAL+2.*KINE)/3./VOLUME*2.3070795E+7 DO 463 I=1,3 DO 463 J=1,3 B(I,J) = B(I,J)/VOLUME*2.3070795E+7 WW = WW+H1(I,J)*H1(I,J) 463 CONTINUE TNOW = KINE*T/SCALER KINE = KINE*2.3070795E-18*AVO/NM/1000. POTE = POTE*2.3070795E-18*AVO/NM/1000. WW = WW*WMASS/2.*2.3070795E-18*AVO/NM/1000. C -------------------------------------------------------- C THE POLARIZATION ENERGY OF WATER MOLECULES IN LIQUIDS: C 5.22 KJOULE/MOL. AND THE LENNARD-JONES CUTOFF ENERGY C CORRECTION. C -------------------------------------------------------- POTE=POTE+5.22-2./3.*PI*2616.91*NM/VOLUME/(CUBE/2.)**3 C ACCUMMULATE SUMS FOR PROPERTY AVERAGES. IF (KB.GT.KSTART) THEN POTESUM = POTESUM + POTE KINESUM = KINESUM + KINE H11 = H11 + H(1,1) H22 = H22 + H(2,2) H33 = H33 + H(3,3) H12 = H12 + H(1,2) H13 = H13 + H(1,3) H23 = H23 + H(2,3) PRESUM = PRESUM + PRES B11 = B11 + B(1,1) B22 = B22 + B(2,2) B33 = B33 + B(3,3) B12 = B12 + B(1,2) B13 = B13 + B(1,3) B23 = B23 + B(2,3) 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,3(F9.3),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) DO 380 I=1,NP X1(I)=X1(I)*FACTOR Y1(I)=Y1(I)*FACTOR Z1(I)=Z1(I)*FACTOR 380 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) = X1(IP) VY(KPT,IM) = Y1(IP) VZ(KPT,IM) = Z1(IP) V2(KPT,IM) = X1(IP)*X1(IP)+Y1(IP)*Y1(IP)+Z1(IP)*Z1(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)+(X1(IP)*VX(KP,IM)+Y1(IP)*VY(KP,IM)+ 1 Z1(IP)*VZ(KP,IM))/V2(KP,IM) IP = IP+3 1400 CONTINUE ENDIF 1600 CONTINUE ENDIF C --------- VELOCITY CORRELATION FUNCTIONAL -------------- print *,pote print *,pote+kine 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(30,113) CORR(I) 113 FORMAT(F9.6) 1500 CONTINUE ENDIF IF (K_POSI.EQ.1) THEN DO 1700 I = 1,NP-2,3 WRITE(31,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 KINESUM = KINESUM/(MAXKB-KSTART) POTESUM = POTESUM/(MAXKB-KSTART) PRESUM = PRESUM/(MAXKB-KSTART) H11 = H11/(MAXKB-KSTART) H22 = H22/(MAXKB-KSTART) H33 = H33/(MAXKB-KSTART) H12 = H12/(MAXKB-KSTART) H13 = H13/(MAXKB-KSTART) H23 = H23/(MAXKB-KSTART) B11 = B11/(MAXKB-KSTART) B22 = B22/(MAXKB-KSTART) B33 = B33/(MAXKB-KSTART) B12 = B12/(MAXKB-KSTART) B13 = B13/(MAXKB-KSTART) B23 = B23/(MAXKB-KSTART) 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) KINESUM*1000/AVO/(3.*BOLZ) 701 FORMAT(5X,' THE AVERAGED TEMPERATURE IS ',F7.2,' KELVIN') WRITE (LP,901) POTESUM 901 FORMAT(5X,' THE AVERAGED POTENTIAL ENERGY IS',F7.2,' KJOULE/MOL') WRITE (LP,1001) PRESUM 1001 FORMAT(5X,' THE AVERAGED PRESSURE IS',F7.2,' BAR',/) WRITE (LP,800) 800 FORMAT(5X,' THE AVERAGED H MATRIX IS:',/) WRITE (LP,801) H11,H12,H13 WRITE (LP,801) H12,H22,H23 WRITE (LP,801) H13,H23,H33 801 FORMAT(7X,'|',3(2X,F10.6,2X),' |') WRITE (LP,1800) 1800 FORMAT(/,5X,' THE AVERAGED STRESS TENSOR IS:',/) WRITE (LP,1801) B11,B12,B13 WRITE (LP,1801) B12,B22,B23 WRITE (LP,1801) B13,B23,B33 1801 FORMAT(7X,'|',3(2X,F10.6,2X),' |') 221 FORMAT(1X,E13.7) 222 FORMAT(6(1X,E13.7)) WRITE (32,221) VOLUME WRITE (32,222) H(1,1),H(2,2),H(3,3),H(1,2),H(1,3),H(2,3) WRITE (32,222) V(1,1),V(2,2),V(3,3),V(1,2),V(1,3),V(2,3) WRITE (32,222) H1(1,1),H1(2,2),H1(3,3),H1(1,2),H1(1,3),H1(2,3) WRITE (32,222) H2(1,1),H2(2,2),H2(3,3),H2(1,2),H2(1,3),H2(2,3) DO 244 I=1,NP WRITE (32,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,RX,RY,RZ,WMASS 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 RELATIVE(I,J) C GIVE THE RELATIVE DISTANCE BETWEEN PARTICLES IMPLICIT REAL(A-H,O-Z) PARAMETER (NPD = 648) COMMON /POS/ X(NPD),Y(NPD),Z(NPD),SX(NPD),SY(NPD),SZ(NPD), 1 X1(NPD),Y1(NPD),Z1(NPD),SX1(NPD),SY1(NPD),SZ1(NPD) DOUBLE PRECISION X1(NPD),Y1(NPD),Z1(NPD) COMMON /CONTROL/ KB,KSTART,SCALER,SEED,DELTA,RX,RY,RZ,WMASS COMMON /H/ H(3,3),HH(3,3),H1(3,3),H2(3,3),V(3,3),P(3,3), 1 B(3,3),CUBE,VOLUME COMMON /CONST/ QO,QH,QOO,QOH,QHH,MO,MH,BOND,HHDIST,ANGLE,PI,P2 DSX = SX(I)-SX(J) DSY = SY(I)-SY(J) DSZ = SZ(I)-SZ(J) IF (DSX.LT.-PI) DSX = DSX+P2 IF (DSX.GT.PI) DSX = DSX-P2 IF (DSY.LT.-PI) DSY = DSY+P2 IF (DSY.GT.PI) DSY = DSY-P2 IF (DSZ.LT.-PI) DSZ = DSZ+P2 IF (DSZ.GT.PI) DSZ = DSZ-P2 RX = HH(1,1)*DSX+HH(1,2)*DSY+HH(1,3)*DSZ RY = HH(2,1)*DSX+HH(2,2)*DSY+HH(2,3)*DSZ RZ = HH(3,1)*DSX+HH(3,2)*DSY+HH(3,3)*DSZ RETURN END C --------------------------------------------------------- C --------------------------------------------------------- SUBROUTINE RATTLE C INTEGRATE MOTION WITH CONSTRAINTS IMPLICIT REAL(A-H,O-Z) PARAMETER (NPD = 648) PARAMETER (ID = 20000) COMMON /EWALD/ ALPHA,ERR(ID),CLI(ID),DERR,A(7,7,7) COMMON /POS/ X(NPD),Y(NPD),Z(NPD),SX(NPD),SY(NPD),SZ(NPD), 1 X1(NPD),Y1(NPD),Z1(NPD),SX1(NPD),SY1(NPD),SZ1(NPD) DOUBLE PRECISION X1(NPD),Y1(NPD),Z1(NPD) COMMON /FORCE/ FX(NPD),FY(NPD),FZ(NPD) COMMON /PARM/ RDEL,IGOFR(3,1500),NM,NPO,NPH,NP,PMASS(NPD), 1 ITYPE(NPD),Q(NPD) COMMON /ENERGY/ TOTE,POTE,KINE,PLJ12,PLJ6,ELEC,VIRIAL COMMON /CONST/ QO,QH,QOO,QOH,QHH,MO,MH,BOND,HHDIST,ANGLE,PI,P2 COMMON /CONTROL/ KB,KSTART,SCALER,SEED,DELTA,RX,RY,RZ,WMASS COMMON /H/ H(3,3),HH(3,3),H1(3,3),H2(3,3),V(3,3),P(3,3), 1 B(3,3),CUBE,VOLUME REAL KINE,MO,MH,KX,KY,KZ,K2 DOUBLE PRECISION QX(NPD),QY(NPD),QZ(NPD),X12(216),Y12(216), 1 Z12(216),X13(216),Y13(216),Z13(216),X23(216), 2 Y23(216),Z23(216),G12,G13,G23,S12X,S12Y, 3 S12Z,S13X,S13Y,S13Z,S23X,S23Y,S23Z,B1,B2, 4 B3,BX,BY,BZ,VX12,VY12,VZ12,VX13,VY13, 5 VZ13,VX23,VY23,VZ23,K12,K13,K23,F1,F2,F3, 6 G12S,G13S,G23S,W12,W13,W23,P11,P22,P33,P12, 7 P13,P23 DIMENSION D(3,3),E(3,3) goto 123 IF (KB.GT.1000.OR.LAST.EQ.1) THEN C ------------ PARRINELLO & RAHMAN FUNCTIONAL ------------ DO 555 I=1,NP XX = X1(I)-H1(1,1)*SX(I)-H1(1,2)*SY(I)-H1(1,3)*SZ(I) YY = Y1(I)-H1(2,1)*SX(I)-H1(2,2)*SY(I)-H1(2,3)*SZ(I) ZZ = Z1(I)-H1(3,1)*SX(I)-H1(3,2)*SY(I)-H1(3,3)*SZ(I) SX1(I) = V(1,1)*XX+V(1,2)*YY+V(1,3)*ZZ SY1(I) = V(2,1)*XX+V(2,2)*YY+V(2,3)*ZZ SZ1(I) = V(3,1)*XX+V(3,2)*YY+V(3,3)*ZZ 555 CONTINUE DO 556 I=1,3 DO 556 J=1,3 D(I,J) = HH(1,I)*H1(1,J)+HH(2,I)*H1(2,J)+HH(3,I)*H1(3,J) 556 CONTINUE DO 557 I=1,3 DO 557 J=1,3 E(I,J) = D(I,J)+D(J,I) 557 CONTINUE DO 558 I=1,3 DO 558 J=1,3 D(I,J) = 2*H1(I,J)-V(1,I)*E(1,J)-V(2,I)*E(2,J)-V(3,I)*E(3,J) 558 CONTINUE DO 559 I=1,NP C1X = D(1,1)*SX1(I)+D(1,2)*SY1(I)+D(1,3)*SZ1(I) C1Y = D(2,1)*SX1(I)+D(2,2)*SY1(I)+D(2,3)*SZ1(I) C1Z = D(3,1)*SX1(I)+D(3,2)*SY1(I)+D(3,3)*SZ1(I) C2X = H2(1,1)*SX(I)+H2(1,2)*SY(I)+H2(1,3)*SZ(I) C2Y = H2(2,1)*SX(I)+H2(2,2)*SY(I)+H2(2,3)*SZ(I) C2Z = H2(3,1)*SX(I)+H2(3,2)*SY(I)+H2(3,3)*SZ(I) FX(I) = FX(I)+PMASS(I)*(C1X+C2X) FY(I) = FY(I)+PMASS(I)*(C1Y+C2Y) FZ(I) = FZ(I)+PMASS(I)*(C1Z+C2Z) 559 CONTINUE DO 745 I=1,3 DO 745 J=1,3 HH(I,J) = HH(I,J)+DELTA*H1(I,J)+DELTA*DELTA*H2(I,J)/2. H(I,J) = HH(I,J)*P2 745 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) VOLUME = H(1,1)*V(1,1)+H(1,2)*V(1,2)+H(1,3)*V(1,3) CUBE = VOLUME**(1./3.) DO 633 I=1,3 DO 633 J=I,3 V(I,J) = V(I,J)*P2/VOLUME V(J,I) = V(I,J) 633 CONTINUE DO 634 I = -3,3 DO 634 J = -3,3 DO 634 K = -3,3 KX = I*V(1,1)+J*V(2,1)+K*V(3,1) KY = I*V(1,2)+J*V(2,2)+K*V(3,2) KZ = I*V(1,3)+J*V(2,3)+K*V(3,3) K2 = KX*KX+KY*KY+KZ*KZ IF (K2.EQ.0) GOTO 634 C1 = EXP(-K2/ALPHA/ALPHA/4.) A(I+4,J+4,K+4) = 4*PI*C1/K2/VOLUME 634 CONTINUE ENDIF C ------------ PARRINELLO & RAHMAN FUNCTIONAL ------------ 123 D12 = BOND*BOND D13 = BOND*BOND D23 = HHDIST*HHDIST C12 = 1./(1./MO+1./MH)/DELTA/2. C13 = C12 C23 = MH/4./DELTA HM1 = DELTA/2./MO HM2 = DELTA/2./MH HM3 = DELTA/2./MH I1 = 1 I2 = 2 I3 = 3 DO 100 I = 1,NM CALL RELATIVE(I1,I2) X12(I) = RX Y12(I) = RY Z12(I) = RZ CALL RELATIVE(I1,I3) X13(I) = RX Y13(I) = RY Z13(I) = RZ CALL RELATIVE(I2,I3) X23(I) = RX Y23(I) = RY Z23(I) = RZ QX(I1) = X1(I1)+HM1*FX(I1) QY(I1) = Y1(I1)+HM1*FY(I1) QZ(I1) = Z1(I1)+HM1*FZ(I1) QX(I2) = X1(I2)+HM2*FX(I2) QY(I2) = Y1(I2)+HM2*FY(I2) QZ(I2) = Z1(I2)+HM2*FZ(I2) QX(I3) = X1(I3)+HM3*FX(I3) QY(I3) = Y1(I3)+HM3*FY(I3) QZ(I3) = Z1(I3)+HM3*FZ(I3) 200 S12X = X12(I)+DELTA*(QX(I1)-QX(I2)) S12Y = Y12(I)+DELTA*(QY(I1)-QY(I2)) S12Z = Z12(I)+DELTA*(QZ(I1)-QZ(I2)) S12 = S12X*S12X+S12Y*S12Y+S12Z*S12Z G12 = (S12-D12)*C12/(S12X*X12(I)+S12Y*Y12(I)+S12Z*Z12(I)) B1 = G12/MO B2 = G12/MH QX(I1) = QX(I1)-B1*X12(I) QY(I1) = QY(I1)-B1*Y12(I) QZ(I1) = QZ(I1)-B1*Z12(I) QX(I2) = QX(I2)+B2*X12(I) QY(I2) = QY(I2)+B2*Y12(I) QZ(I2) = QZ(I2)+B2*Z12(I) S13X = X13(I)+DELTA*(QX(I1)-QX(I3)) S13Y = Y13(I)+DELTA*(QY(I1)-QY(I3)) S13Z = Z13(I)+DELTA*(QZ(I1)-QZ(I3)) S13 = S13X*S13X+S13Y*S13Y+S13Z*S13Z G13 = (S13-D13)*C13/(S13X*X13(I)+S13Y*Y13(I)+S13Z*Z13(I)) B1 = G13/MO B3 = G13/MH QX(I1) = QX(I1)-B1*X13(I) QY(I1) = QY(I1)-B1*Y13(I) QZ(I1) = QZ(I1)-B1*Z13(I) QX(I3) = QX(I3)+B3*X13(I) QY(I3) = QY(I3)+B3*Y13(I) QZ(I3) = QZ(I3)+B3*Z13(I) S23X = X23(I)+DELTA*(QX(I2)-QX(I3)) S23Y = Y23(I)+DELTA*(QY(I2)-QY(I3)) S23Z = Z23(I)+DELTA*(QZ(I2)-QZ(I3)) S23 = S23X*S23X+S23Y*S23Y+S23Z*S23Z G23 = (S23-D23)*C23/(S23X*X23(I)+S23Y*Y23(I)+S23Z*Z23(I)) B2 = G23/MH BX = B2*X23(I) BY = B2*Y23(I) BZ = B2*Z23(I) QX(I2) = QX(I2)-BX QY(I2) = QY(I2)-BY QZ(I2) = QZ(I2)-BZ QX(I3) = QX(I3)+BX QY(I3) = QY(I3)+BY QZ(I3) = QZ(I3)+BZ IF (BX.GT.5D-8.OR.BY.GT.5D-8.OR.BZ.GT.5D-8) GOTO 200 IF (BX.LT.-5D-8.OR.BY.LT.-5D-8.OR.BZ.LT.-5D-8) GOTO 200 DX1 = DELTA*QX(I1) DY1 = DELTA*QY(I1) DZ1 = DELTA*QZ(I1) X(I1) = X(I1)+DX1 Y(I1) = Y(I1)+DY1 Z(I1) = Z(I1)+DZ1 DX2 = DELTA*QX(I2) DY2 = DELTA*QY(I2) DZ2 = DELTA*QZ(I2) X(I2) = X(I2)+DX2 Y(I2) = Y(I2)+DY2 Z(I2) = Z(I2)+DZ2 DX3 = DELTA*QX(I3) DY3 = DELTA*QY(I3) DZ3 = DELTA*QZ(I3) X(I3) = X(I3)+DX3 Y(I3) = Y(I3)+DY3 Z(I3) = Z(I3)+DZ3 X12(I) = X12(I)+DX1-DX2 Y12(I) = Y12(I)+DY1-DY2 Z12(I) = Z12(I)+DZ1-DZ2 X13(I) = X13(I)+DX1-DX3 Y13(I) = Y13(I)+DY1-DY3 Z13(I) = Z13(I)+DZ1-DZ3 X23(I) = X23(I)+DX2-DX3 Y23(I) = Y23(I)+DY2-DY3 Z23(I) = Z23(I)+DZ2-DZ3 I1 = I1+3 I2 = I2+3 I3 = I3+3 100 CONTINUE DO 500 I = 1,NP SX(I) = V(1,1)*X(I)+V(1,2)*Y(I)+V(1,3)*Z(I) SY(I) = V(2,1)*X(I)+V(2,2)*Y(I)+V(2,3)*Z(I) SZ(I) = V(3,1)*X(I)+V(3,2)*Y(I)+V(3,3)*Z(I) IF (SX(I).LT.0) THEN SX(I) = SX(I)+P2 X(I) = X(I)+H(1,1) Y(I) = Y(I)+H(2,1) Z(I) = Z(I)+H(3,1) ENDIF IF (SX(I).GT.P2) THEN SX(I) = SX(I)-P2 X(I) = X(I)-H(1,1) Y(I) = Y(I)-H(2,1) Z(I) = Z(I)-H(3,1) ENDIF IF (SY(I).LT.0) THEN SY(I) = SY(I)+P2 X(I) = X(I)+H(1,2) Y(I) = Y(I)+H(2,2) Z(I) = Z(I)+H(3,2) ENDIF IF (SY(I).GT.P2) THEN SY(I) = SY(I)-P2 X(I) = X(I)-H(1,2) Y(I) = Y(I)-H(2,2) Z(I) = Z(I)-H(3,2) ENDIF IF (SZ(I).LT.0) THEN SZ(I) = SZ(I)+P2 X(I) = X(I)+H(1,3) Y(I) = Y(I)+H(2,3) Z(I) = Z(I)+H(3,3) ENDIF IF (SZ(I).GT.P2) THEN SZ(I) = SZ(I)-P2 X(I) = X(I)-H(1,3) Y(I) = Y(I)-H(2,3) Z(I) = Z(I)-H(3,3) ENDIF 500 CONTINUE CALL INTERACT E12 = 1./D12/(1./MH+1./MO) E13 = E12 E23 = MH/D23/2. G12S = 0. G13S = 0. G23S = 0. P11 = 0. P22 = 0. P33 = 0. P12 = 0. P13 = 0. P23 = 0. I1 = 1 I2 = 2 I3 = 3 DO 300 I = 1,NM W12 = 0. W13 = 0. W23 = 0. X1(I1) = QX(I1)+HM1*FX(I1) Y1(I1) = QY(I1)+HM1*FY(I1) Z1(I1) = QZ(I1)+HM1*FZ(I1) X1(I2) = QX(I2)+HM2*FX(I2) Y1(I2) = QY(I2)+HM2*FY(I2) Z1(I2) = QZ(I2)+HM2*FZ(I2) X1(I3) = QX(I3)+HM3*FX(I3) Y1(I3) = QY(I3)+HM3*FY(I3) Z1(I3) = QZ(I3)+HM3*FZ(I3) 400 VX12 = X1(I1)-X1(I2) VY12 = Y1(I1)-Y1(I2) VZ12 = Z1(I1)-Z1(I2) K12 = (X12(I)*VX12+Y12(I)*VY12+Z12(I)*VZ12)*E12 W12 = W12+K12 F1 = K12/MO X1(I1) = X1(I1)-F1*X12(I) Y1(I1) = Y1(I1)-F1*Y12(I) Z1(I1) = Z1(I1)-F1*Z12(I) F2 = K12/MH X1(I2) = X1(I2)+F2*X12(I) Y1(I2) = Y1(I2)+F2*Y12(I) Z1(I2) = Z1(I2)+F2*Z12(I) VX13 = X1(I1)-X1(I3) VY13 = Y1(I1)-Y1(I3) VZ13 = Z1(I1)-Z1(I3) K13 = (X13(I)*VX13+Y13(I)*VY13+Z13(I)*VZ13)*E13 W13 = W13+K13 F1 = K13/MO X1(I1) = X1(I1)-F1*X13(I) Y1(I1) = Y1(I1)-F1*Y13(I) Z1(I1) = Z1(I1)-F1*Z13(I) F3 = K13/MH X1(I3) = X1(I3)+F3*X13(I) Y1(I3) = Y1(I3)+F3*Y13(I) Z1(I3) = Z1(I3)+F3*Z13(I) VX23 = X1(I2)-X1(I3) VY23 = Y1(I2)-Y1(I3) VZ23 = Z1(I2)-Z1(I3) K23 = (X23(I)*VX23+Y23(I)*VY23+Z23(I)*VZ23)*E23 W23 = W23+K23 F3 = K23/MH BX = F3*X23(I) BY = F3*Y23(I) BZ = F3*Z23(I) X1(I2) = X1(I2)-BX Y1(I2) = Y1(I2)-BY Z1(I2) = Z1(I2)-BZ X1(I3) = X1(I3)+BX Y1(I3) = Y1(I3)+BY Z1(I3) = Z1(I3)+BZ IF (BX.GT.5D-8.OR.BY.GT.5D-8.OR.BZ.GT.5D-8) GOTO 400 IF (BX.LT.-5D-8.OR.BY.LT.-5D-8.OR.BZ.LT.-5D-8) GOTO 400 P11 = P11+W12*X12(I)*X12(I)+W13*X13(I)*X13(I)+W23*X23(I)*X23(I) P22 = P22+W12*Y12(I)*Y12(I)+W13*Y13(I)*Y13(I)+W23*Y23(I)*Y23(I) P33 = P33+W12*Z12(I)*Z12(I)+W13*Z13(I)*Z13(I)+W23*Z23(I)*Z23(I) P12 = P12+W12*X12(I)*Y12(I)+W13*X13(I)*Y13(I)+W23*X23(I)*Y23(I) P13 = P13+W12*X12(I)*Z12(I)+W13*X13(I)*Z13(I)+W23*X23(I)*Z23(I) P23 = P23+W12*Y12(I)*Z12(I)+W13*Y13(I)*Z13(I)+W23*Y23(I)*Z23(I) G12S = G12S+W12 G13S = G13S+W13 G23S = G23S+W23 I1 = I1+3 I2 = I2+3 I3 = I3+3 300 CONTINUE VIRIAL = VIRIAL-2./DELTA*(G12S*D12+G13S*D13+G23S*D23) B(1,1) = B(1,1)-2./DELTA*P11 B(2,2) = B(2,2)-2./DELTA*P22 B(3,3) = B(3,3)-2./DELTA*P33 B(1,2) = B(1,2)-2./DELTA*P12 B(1,3) = B(1,3)-2./DELTA*P13 B(2,3) = B(2,3)-2./DELTA*P23 KINE = 0. DO 260 I=1,NP XX = X1(I)-H1(1,1)*SX(I)-H1(1,2)*SY(I)-H1(1,3)*SZ(I) YY = Y1(I)-H1(2,1)*SX(I)-H1(2,2)*SY(I)-H1(2,3)*SZ(I) ZZ = Z1(I)-H1(3,1)*SX(I)-H1(3,2)*SY(I)-H1(3,3)*SZ(I) KINE = KINE+(XX*XX+YY*YY+ZZ*ZZ)*PMASS(I) B(1,1) = B(1,1)+XX*XX*PMASS(I) B(2,2) = B(2,2)+YY*YY*PMASS(I) B(3,3) = B(3,3)+ZZ*ZZ*PMASS(I) B(1,2) = B(1,2)+XX*YY*PMASS(I) B(1,3) = B(1,3)+XX*ZZ*PMASS(I) B(2,3) = B(2,3)+YY*ZZ*PMASS(I) 260 CONTINUE KINE = KINE/2. DO 388 I=1,3 DO 388 J=I,3 D(I,J) = B(I,J)-P(I,J)*VOLUME/2.3070795E+7 D(J,I) = D(I,J) 388 CONTINUE return IF (KB.GT.1000.OR.LAST.EQ.1) THEN DO 389 I=1,3 DO 389 J=I,3 H2IJ = (D(I,1)*V(J,1)+D(I,2)*V(J,2)+D(I,3)*V(J,3))/WMASS H1(I,J) = H1(I,J)+DELTA*(H2(I,J)+H2IJ)/2. H2(I,J) = H2IJ H1(J,I) = H1(I,J) H2(J,I) = H2(I,J) 389 CONTINUE ENDIF 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),SX(NPD),SY(NPD),SZ(NPD), 1 X1(NPD),Y1(NPD),Z1(NPD),SX1(NPD),SY1(NPD),SZ1(NPD) DOUBLE PRECISION X1(NPD),Y1(NPD),Z1(NPD) COMMON /PARM/ RDEL,IGOFR(3,1500),NM,NPO,NPH,NP,PMASS(NPD), 1 ITYPE(NPD),Q(NPD) COMMON /CONST/ QO,QH,QOO,QOH,QHH,MO,MH,BOND,HHDIST,ANGLE,PI,P2 COMMON /CONTROL/ KB,KSTART,SCALER,SEED,DELTA,RX,RY,RZ,WMASS 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(7,7,7) COMMON /POS/ X(NPD),Y(NPD),Z(NPD),SX(NPD),SY(NPD),SZ(NPD), 1 X1(NPD),Y1(NPD),Z1(NPD),SX1(NPD),SY1(NPD),SZ1(NPD) DOUBLE PRECISION X1(NPD),Y1(NPD),Z1(NPD) COMMON /FORCE/ FX(NPD),FY(NPD),FZ(NPD) COMMON /PARM/ RDEL,IGOFR(3,1500),NM,NPO,NPH,NP,PMASS(NPD), 1 ITYPE(NPD),Q(NPD) COMMON /ENERGY/ TOTE,POTE,KINE,PLJ12,PLJ6,ELEC,VIRIAL COMMON /CONST/ QO,QH,QOO,QOH,QHH,MO,MH,BOND,HHDIST,ANGLE,PI,P2 COMMON /CONTROL/ KB,KSTART,SCALER,SEED,DELTA,RX,RY,RZ,WMASS COMMON /H/ H(3,3),HH(3,3),H1(3,3),H2(3,3),V(3,3),P(3,3), 1 B(3,3),CUBE,VOLUME REAL KINE,MO,MH DIMENSION SINN(NPD,7,7,7),COSN(NPD,7,7,7) ELEC = 0. PLJ6 = 0. PLJ12 = 0. DO 300 I=1,3 DO 300 J=I,3 B(I,J) = 0. 300 CONTINUE W = 2.*ALPHA/SQRT(PI) RCUT = CUBE/2. 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 SINX = SIN(SX(I)) COSX = COS(SX(I)) SINY = SIN(SY(I)) COSY = COS(SY(I)) SINZ = SIN(SZ(I)) COSZ = COS(SZ(I)) SINN(I,1,1,1) = SIN(-3.*SX(I)-3.*SY(I)-3.*SZ(I)) COSN(I,1,1,1) = COS(-3.*SX(I)-3.*SY(I)-3.*SZ(I)) DO 1100 J=2,7 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,7 DO 1200 K=2,7 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,7 DO 1300 K=1,7 DO 1300 L=2,7 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,7 DO 1400 K=1,7 DO 1400 L=1,7 QX = (J-4)*V(1,1)+(K-4)*V(2,1)+(L-4)*V(3,1) QY = (J-4)*V(1,2)+(K-4)*V(2,2)+(L-4)*V(3,2) QZ = (J-4)*V(1,3)+(K-4)*V(2,3)+(L-4)*V(3,3) Q2 = QX*QX+QY*QY+QZ*QZ IF (Q2.EQ.0) GOTO 1400 QK = 2./Q2+1./2./ALPHA/ALPHA 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)+QX*CC FY(I) = FY(I)+QY*CC FZ(I) = FZ(I)+QZ*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-1 DO 5520 J=I+1,NP CALL RELATIVE(I,J) 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 ((J-1)/3.NE.(I-1)/3) THEN IR = R1/RDEL+1 IGOFR(IJT,IR) = IGOFR(IJT,IR) + 1 ENDIF C CALCULATE EWARD FACTOR IN REAL SPACE SUMMATION PRODUCT = ALPHA*R1 I0 = PRODUCT/DERR XX = I0*DERR ARGIN = PRODUCT-XX GACTOR = ERR(I0)-CLI(I0)*ARGIN*(1-XX*ARGIN) C IF THEY ARE THE SAME MOLECULE IF ((J-1)/3.EQ.(I-1)/3) GACTOR = GACTOR-1. ARGIN = PRODUCT*PRODUCT-XX*XX 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 = 1895.2940/R12 E2 = 1.8835408/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 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*RX FFY= FF*RY FFZ= FF*RZ B(1,1) = B(1,1)+RX*FFX B(2,2) = B(2,2)+RY*FFY B(3,3) = B(3,3)+RZ*FFZ B(1,2) = B(1,2)+RX*FFY B(1,3) = B(1,3)+RX*FFZ B(2,3) = B(2,3)+RY*FFZ 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 5520 CONTINUE 522 CONTINUE C CALCULATE TOTAL POTENTIAL ENERGIES POTE = PLJ12+PLJ6+ELEC 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),SX(NPD),SY(NPD),SZ(NPD), 1 X1(NPD),Y1(NPD),Z1(NPD),SX1(NPD),SY1(NPD),SZ1(NPD) DOUBLE PRECISION X1(NPD),Y1(NPD),Z1(NPD) COMMON /PARM/ RDEL,IGOFR(3,1500),NM,NPO,NPH,NP,PMASS(NPD), 1 ITYPE(NPD),Q(NPD) COMMON /CONST/ QO,QH,QOO,QOH,QHH,MO,MH,BOND,HHDIST,ANGLE,PI,P2 COMMON /CONTROL/ KB,KSTART,SCALER,SEED,DELTA,RX,RY,RZ,WMASS COMMON /H/ H(3,3),HH(3,3),H1(3,3),H2(3,3),V(3,3),P(3,3), 1 B(3,3),CUBE,VOLUME 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/ NX = NINT((NM/8)**(1./3.)) NY = NX NZ = NX A = CUBE/NX SX(7) = 1./2. SY(7) = 1./2. SY(10) = 1./2. SZ(10) = 1./2. SZ(4)= 1./2. SX(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,NX-1 DO 200 J = 0,NY-1 DO 200 K = 0,NZ-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 SX(I) = X(I)*P2/CUBE SY(I) = Y(I)*P2/CUBE SZ(I) = Z(I)*P2/CUBE 700 CONTINUE RETURN END C ---------------------------------------------------