C ---------------------------------------------------------------- C PROGRAM ICE (FLEXIBLE MODEL) JUL.12,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 ---------------------------------------------------------------- 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 = 50000) PARAMETER (IV = 300) 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 REAL KINE,MO,MH REAL MSD,KINESUM DIMENSION V2(IV,216),VX(IV,216),VY(IV,216),VZ(IV,216),CORR(5000), 1 X2(NPD),Y2(NPD),Z2(NPD),DX(NPD),DY(NPD),DZ(NPD) C 2 VMX(216),VMY(216),VMZ(216) DATA LP,LP_GR,LP_CORR,LP_POSI,LP_CONFIG,LP_BAK /25,26,27,28,29,30/ C SELECT 300 TIME ORIGINS, EACH WITH 10 STEPS APART C TOTAL CORRELATION LENGTH = 5000 STEPS = 2 ps KVSTEP = 3000/IV AVO=6.0221367E+23 BOLZ=1.380658E-23 AMU=1.6605E-24 C ONE TIME STEP = 0.0004ps DELTA = 0.149096206 NM = 216 NPO = 1 NPH = 2 NP = NM*(NPO+NPH) C BERENDSEN T,P COUPLING CONSTANT TAUT = 0.05/2.6828315E-3 TAUP = 0.5/2.6828315E-3 BETA = 4.9E-5 FACTORT = 1. FACTORP = 1. C A12 = 1895.294007 C A6 = -1.883540782 C QO = -0.8476 PI = 3.141592654 P2 = 2.*PI 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 = (104+52./100.)/180.*PI HHDIST= 2.*SIN(ANGLE/2.)*BOND 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,P,MAXKB,KSTART READ *,LAST,K_CORR,K_POSI READ *,DENSITY TOTALMASS = NM*AMASS VOLUME = TOTALMASS*1.66054/DENSITY CUBE = VOLUME**(1./3.) FROG = P2/CUBE C SCALER IS THE IDEAL TOTAL KINETIC ENERGY OF CM SCALER = 1.5*BOLZ*T*NM/2.3070795E-18 IF (LAST.EQ.0) THEN CALL DIAMOND CALL INTVEL ELSE DO 377 I=1,NP READ *,X(I),Y(I),Z(I),X1(I),Y1(I),Z1(I) 377 CONTINUE ENDIF RDEL = CUBE/2./500. DO 153 I=1,1500 IGOFR(1,I) = 0 IGOFR(2,I) = 0 IGOFR(3,I) = 0 153 CONTINUE DO 154 I = 1,5000 CORR(I) = 0. 154 CONTINUE 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 KWBAK IS THE FREQUENCY OF WRITING ON CONFIG.BAK C KWGR IS THE FREQUENCY OF WRITING ON GR.OUT C -------------------------------------------------------------- KWRITE = 20 KWGR = 20000 KWBAK = 10160 C SET CONSTS AND TABLES FOR EWALD SUMMATION ALPHA = 5.714/CUBE DERR = ALPHA*CUBE/2.*SQRT(3.)*1.2 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,80) T,P 80 FORMAT(7X,'*',2X,'T(K) =',F8.3,T30,'P(BAR) = ', 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) 9*9*9-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',5X,'TIME',4X,'TEMP',7X,'PRES',4X,'DENSITY', 1 4X,'POTE',5X,'TOTE',7X,'MSD') DO 55 I=1,NP DX(I) = 0. DY(I) = 0. DZ(I) = 0. X2(I) = 0. Y2(I) = 0. Z2(I) = 0. 55 CONTINUE TSUM = 0. VINSUM = 0. KINESUM = 0. POTESUM = 0. PRESUM = 0. BONSUM = 0. ANGSUM = 0. PM2SUM = 0. PMXSUM = 0. PMYSUM = 0. PMZSUM = 0. DENSUM = 0. C ----------------------------------------------------------- C ENTER MAIN LOOP OF SIMULATION DO 310 KB=1,MAXKB CUBE = CUBE*FACTORP FROG = FROG/FACTORP ALPHA = ALPHA/FACTORP VOLUME = VOLUME*FACTORP**3 DENSITY = DENSITY/FACTORP**3 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)) DX(I) = DX(I) + DXX DY(I) = DY(I) + DYY DZ(I) = DZ(I) + DZZ 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 X(I) = X(I)*FACTORP Y(I) = Y(I)*FACTORP Z(I) = Z(I)*FACTORP 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) = FACTORT*X1(I)+DELTA*(X2(I)+XX2)/2. Y1(I) = FACTORT*Y1(I)+DELTA*(Y2(I)+YY2)/2. Z1(I) = FACTORT*Z1(I)+DELTA*(Z2(I)+ZZ2)/2. X2(I) = XX2 Y2(I) = YY2 Z2(I) = ZZ2 KINE = KINE+PMASS(I)*(X1(I)*X1(I)+Y1(I)*Y1(I)+Z1(I)*Z1(I)) 330 CONTINUE KINE = KINE/2. MSD = 0. IP = 1 DO 66 I=1,NM C VMX(I) = X1(IP)*OPART+(X1(IP+1)+X1(IP+2))*HPART C VMY(I) = Y1(IP)*OPART+(Y1(IP+1)+Y1(IP+2))*HPART C VMZ(I) = Z1(IP)*OPART+(Z1(IP+1)+Z1(IP+2))*HPART XX = DX(IP)*OPART+(DX(IP+1)+DX(IP+2))*HPART YY = DY(IP)*OPART+(DY(IP+1)+DY(IP+2))*HPART ZZ = DZ(IP)*OPART+(DZ(IP+1)+DZ(IP+2))*HPART MSD = MSD+XX*XX+YY*YY+ZZ*ZZ IP = IP+3 66 CONTINUE MSD = MSD/NM 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 ------------- C PRESSURE WILL BE IN BARS & ENERGIES WILL BE IN KJOULE/MOL. PRES = (VIRIAL+2.*KINE)/3./VOLUME*2.3070795E+7 TNOW = KINE *T/(SCALER*3) 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 DENSUM = DENSUM + DENSITY 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,PRES,DENSITY,POTE,KINE+POTE,MSD 270 FORMAT(I6,2(F9.3),F11.3,F9.3,2(F10.3),F9.2) ENDIF C -------------------------------------------------------- FACTORT = 1+DELTA/2./TAUT*(T/TNOW-1) FACTORP = 1+BETA*DELTA/3./TAUP*(PRES-P) IF (K_CORR.EQ.1) THEN C --------- VELOCITY CORRELATION FUNCTIONAL -------------- C CALCULATE VELOCITY AUTOCORRELATION FUNCTION I = KB+8000-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 -------------- IF (MOD(KB,KWBAK).EQ.0) THEN WRITE (LP_BAK,221) DENSITY 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 *,'kb = ',kb print *,'pote =',pote print *,'tote =',pote+kine print *,'T = ',tnow print *,'pres = ',pres print *,'density = ',density 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) DENSUM = DENSUM/(MAXKB-KSTART) PMXSUM = PMXSUM/(MAXKB-KSTART) PMYSUM = PMYSUM/(MAXKB-KSTART) PMZSUM = PMZSUM/(MAXKB-KSTART) PM2SUM = PM2SUM/(MAXKB-KSTART) 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,1001) PRESUM 1001 FORMAT(5X,' THE AVERAGED PRESSURE IS',F10.2,' BAR') WRITE (LP,1011) DENSUM 1011 FORMAT(5X,' THE AVERAGED DENSITY IS ',F7.5,' G/CM^3') 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,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) 221 FORMAT(1X,E13.7) WRITE (LP_CONFIG,221) DENSITY 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 C MULTIPLY BY 6, FOR WATER ONLY. JUN.14 FACTOR=SQRT(6*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 = 50000) 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 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),DR(3,3) DIMENSION CMX(216),CMY(216),CMZ(216) QUOTAH = MH/(MO+2*MH) QUOTAO = MO/(MO+2*MH) 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 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 QX = (J-5)*FROG QY = (K-5)*FROG QZ = (L-5)*FROG Q2 = QX*QX+QY*QY+QZ*QZ IF (Q2.EQ.0.) GOTO 1400 C1 = EXP(-Q2/ALPHA/ALPHA/4.) A(J,K,L) = 4.*PI*C1/Q2/VOLUME 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)) 200 FX(I) = FX(I)+QX*CC FY(I) = FY(I)+QY*CC FZ(I) = FZ(I)+QZ*CC 1600 CONTINUE ELEC = ELEC+A(J,K,L)*(COSSUM*COSSUM+SINSUM*SINSUM)/2. 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 RX = X(I)-X(J) RY = Y(I)-Y(J) RZ = Z(I)-Z(J) IF (RX.GT.RCUT) RX = RX - CUBE IF (RX.LT.TUCR) RX = RX + CUBE IF (RY.GT.RCUT) RY = RY - CUBE IF (RY.LT.TUCR) RY = RY + CUBE IF (RZ.GT.RCUT) RZ = RZ - CUBE IF (RZ.LT.TUCR) RZ = RZ + CUBE R2 = RX*RX+RY*RY+RZ*RZ R1 = SQRT(R2) R3 = R2*R1 IJT = ITYPE(I)+ITYPE(J)+1 IF (IMI.NE.IMJ) THEN IR = R1/RDEL+1 IF (IR.LE.500) IGOFR(IJT,IR) = IGOFR(IJT,IR)+1 ENDIF C CALCULATE EWALD FACTOR IN REAL SPACE SUMMATION PRODUCT = ALPHA*R1 I0 = PRODUCT/DERR XX = I0*DERR ARGIN = PRODUCT-XX IF (I0.GE.ID) GOTO 5520 GACTOR = ERR(I0)-CLI(I0)*ARGIN*(1-XX*ARGIN) C IF THEY ARE THE SAME MOLECULE IF (IMI.EQ.IMJ) THEN GACTOR = GACTOR-1. XR(II,IJ) = RX YR(II,IJ) = RY ZR(II,IJ) = RZ DR(II,IJ) = R1 ENDIF 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 = 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.3) 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*RX FFY= FF*RY FFZ= FF*RZ 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 IF (II.EQ.3) THEN PMX = PMX+XR(1,2)+XR(1,3) PMY = PMY+YR(1,2)+YR(1,3) PMZ = PMZ+ZR(1,2)+ZR(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 DR11 = DR1*DR1 DR22 = DR2*DR2 DR33 = DR3*DR3 DR12 = DR1*DR2 DR13 = DR1*DR3 DR23 = DR2*DR3 DR111 = DR11*DR1 DR222 = DR22*DR2 DR333 = DR33*DR3 DR112 = DR11*DR2 DR113 = DR11*DR3 DR122 = DR22*DR1 DR123 = DR12*DR3 DR133 = DR33*DR1 DR223 = DR22*DR3 DR233 = DR33*DR2 F1 = 3.959608852*DR1+0.2514577891*DR2-0.5348000472*DR3 G1 = -13.52143988*DR11-0.768247786*DR12 1 +1.233829205*DR13-0.3841238925*DR22 2 +0.944571334*DR23-0.6852521196*DR33 F2 = 3.959608852*DR2+0.2514577891*DR1-0.5348000472*DR3 G2 = -13.52143988*DR22-0.768247786*DR12 1 +1.233829205*DR23-0.3841238925*DR11 2 +0.944571334*DR13-0.6852521196*DR33 F3 = -0.5348000472*(DR1+DR2)+0.8804540624*DR3 G3 = 0.6169146026*(DR11+DR22)+0.6267843126*DR33+ 1 0.944571334*DR12-1.370504239*(DR13+DR23) H1 = 24.83326248*DR111-1.017614014*DR333 1 +6.876821916*DR113-6.767337447*DR112 2 -1.100316147*DR133-5.500673460*DR122 3 +5.872225848*DR123-2.255779147*DR222 4 +2.936112915*DR223+1.223369784*DR233 H2 = -1.017614014*DR333+24.83326249*DR222 1 +6.876821910*DR223-2.255779149*DR111 2 -5.500673460*DR112+2.936112924*DR113 3 -6.767337441*DR122+5.872225830*DR123 4 +1.223369784*DR133-1.100316147*DR233 H3 = -3.052842042*(DR133+DR233)+2.292273970*(DR111+DR222) 1 -1.100316147*(DR113+DR223)+2.936112924*(DR112+DR122) 2 +2.446739568*DR123+1.274486629*DR333 C SET CUTOFF SO THAT WON'T COLLAPSE DUE TO 3RD AND 4TH ORDER TERMS. IF (DR1.GT.0.07.OR.DR1.LT.-0.07.OR.DR2.GT.0.07.OR. 1 DR2.LT.-0.07.OR.DR3.GT.0.14.OR.DR3.LT.-0.14) THEN G1 = 0. H1 = 0. G2 = 0. H2 = 0. G3 = 0. H3 = 0. ENDIF VINTRA = VINTRA+(F1/2.+G1/3.+H1/4.)*DR1+(F2/2.+G2/3.+H2/4.)*DR2 1 +(F3/2.+G3/3.+H3/4.)*DR3 F1 = F1+G1+H1 F2 = F2+G2+H2 F3 = F3+G3+H3 FX(I-2) = FX(I-2)-F1*XR(1,2)/DR(1,2)-F2*XR(1,3)/DR(1,3) FY(I-2) = FY(I-2)-F1*YR(1,2)/DR(1,2)-F2*YR(1,3)/DR(1,3) FZ(I-2) = FZ(I-2)-F1*ZR(1,2)/DR(1,2)-F2*ZR(1,3)/DR(1,3) FX(I-1) = FX(I-1)+F1*XR(1,2)/DR(1,2)-F3*XR(2,3)/DR(2,3) FY(I-1) = FY(I-1)+F1*YR(1,2)/DR(1,2)-F3*YR(2,3)/DR(2,3) FZ(I-1) = FZ(I-1)+F1*ZR(1,2)/DR(1,2)-F3*ZR(2,3)/DR(2,3) FX(I) = FX(I)+F2*XR(1,3)/DR(1,3)+F3*XR(2,3)/DR(2,3) FY(I) = FY(I)+F2*YR(1,3)/DR(1,3)+F3*YR(2,3)/DR(2,3) FZ(I) = FZ(I)+F2*ZR(1,3)/DR(1,3)+F3*ZR(2,3)/DR(2,3) VIRIAL = VIRIAL-F1*DR(1,2)-F2*DR(1,3)-F3*DR(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 ---------------------------------------------------