C ---------------------------------------------------------------- C PROGRAM MINI (FLEXIBLE MODEL) MARCH.7, 1995 LIJU. C THIS PROGRAM IS FOR THE ENERGY MINIMIZATION OF A C 216 WATER MOLECULES SYSTEM IN DIAMOND CUBIC STRUCTURE C USING CONSTRAINED SIMULATED ANNEALING. 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) 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,BXX,BYY,BZZ COMMON /OTHER/ AMASS,HPART,OPART,DX(NPD),DY(NPD),DZ(NPD) REAL KINE,MO,MH DIMENSION KTLIST(200,2),X2(NPD),Y2(NPD),Z2(NPD) C KTLIST IS THE CONTROL PANEL FOR SIMULATED ANNEALING DATA LP_POSI,LP_CONFIG,LP_BAK /25,26,27/ KWBAK = 10 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 A12 = 1895.294007 C A6 = -1.883540782 C QO = -0.8476 C SPRING IS THE SPRING CONSTANT WHICH ALIGNS CM TO C DIAMOND CUBIC STRUCTURE READ *,A12,A6,QO,SPRING,DENSITY READ *,LAST,KSCALE READ *,KTERMS DO 111 I=1,KTERMS READ *,KTLIST(I,1),KTLIST(I,2) 111 CONTINUE KPT = 1 T = KTLIST(1,1) 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 TOTALMASS = NM*AMASS VOLUME = TOTALMASS*1.66054/DENSITY CUBE = VOLUME**(1./3.) FROG = P2/CUBE SCALER = 4.5*BOLZ*NM/2.3070795E-18 IF (LAST.EQ.0) THEN CALL DIAMOND ELSE READ *,CUBE1,T1 DO 366 I=1,NP READ *,DX(I),DY(I),DZ(I) 366 CONTINUE 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 377 CONTINUE ENDIF OPEN (UNIT=LP_CONFIG, FORM='FORMATTED',FILE="config") OPEN (UNIT=LP_BAK, FORM='FORMATTED',FILE="config.bak") OPEN (UNIT=LP_POSI, FORM='FORMATTED',FILE="posi.out") 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 DO 600 I = -4,4 DO 600 J = -4,4 DO 600 K = -4,4 QX = I*FROG QY = J*FROG QZ = K*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 DO 55 I=1,NP X2(I) = 0. Y2(I) = 0. Z2(I) = 0. 55 CONTINUE KB = 0 BONSUM = 0. ANGSUM = 0. BXXSUM = 0. BYYSUM = 0. BZZSUM = 0. 310 CONTINUE DO 754 MM = 1,KTLIST(KPT,2) KB = KB+1 print *,kb 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 DX(I) = DX(I)+DXX DY(I) = DY(I)+DYY DZ(I) = DZ(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 290 I=1,NP IF (MOD(I,3).EQ.1) THEN XX = OPART*DX(I)+HPART*(DX(I+1)+DX(I+2)) YY = OPART*DY(I)+HPART*(DY(I+1)+DY(I+2)) ZZ = OPART*DZ(I)+HPART*(DZ(I+1)+DZ(I+2)) XX = XX*SPRING/AMASS YY = YY*SPRING/AMASS ZZ = ZZ*SPRING/AMASS ENDIF XX2 = FX(I)/PMASS(I)-XX YY2 = FY(I)/PMASS(I)-YY ZZ2 = FZ(I)/PMASS(I)-ZZ 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 KINE = KINE+PMASS(I)*(X1(I)**2+Y1(I)**2+Z1(I)**2) 290 CONTINUE KINE = KINE/2. TNOW = KINE/SCALER POTE = POTE*2.3070795E-18*AVO/NM/1000. BXX = BXX/VOLUME*2.3070795E+7 BYY = BYY/VOLUME*2.3070795E+7 BZZ = BZZ/VOLUME*2.3070795E+7 BONSUM = BONSUM + BON ANGSUM = ANGSUM + ANG BXXSUM = BXXSUM + BXX BYYSUM = BYYSUM + BYY BZZSUM = BZZSUM + BZZ TALL = TALL+TNOW IF (MOD(KB,KSCALE).EQ.0) THEN TALL=TALL/KSCALE FACTOR=(SQRT(T/TALL)-1.)/2. DO 380 I=1,NP X1(I)=X1(I)*(1+FACTOR) Y1(I)=Y1(I)*(1+FACTOR) Z1(I)=Z1(I)*(1+FACTOR) 380 CONTINUE TALL=0. ENDIF IF (MOD(KB,KWBAK).EQ.0) THEN BONSUM = BONSUM/KWBAK ANGSUM = ANGSUM/KWBAK/PI*180. BXXSUM = BXXSUM/KWBAK BYYSUM = BYYSUM/KWBAK BZZSUM = BZZSUM/KWBAK 190 FORMAT(//,6X,'KB',5X,'TEMP',5X,'POTE',5X,'BOND', 1 4X,'ANGLE',6X,'BXX',6X,'BYY',6X,'BZZ') WRITE(LP_BAK,190) 270 FORMAT(I8,4(F9.3),3(F9.1),/) WRITE(LP_BAK,270) KB,TNOW,POTE,BONSUM,ANGSUM,BXX, 1 BYY,BZZ WRITE (LP_BAK,221) CUBE,TNOW DO 887 I=1,NP WRITE (LP_BAK,232) DX(I),DY(I),DZ(I) 887 CONTINUE 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) 4200 FORMAT('B') BONSUM = 0. ANGSUM = 0. BXXSUM = 0. BYYSUM = 0. BZZSUM = 0. ENDIF T = T-1./KTLIST(KPT,2) 754 CONTINUE IF (T.LE.KTLIST(KPT+1,1)) KPT=KPT+1 IF (KPT.LT.KTERMS) GOTO 310 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 221 FORMAT(2(1X,E13.7)) WRITE (LP_CONFIG,221) CUBE,TNOW 232 FORMAT(3(1X,E13.7)) DO 655 I=1,NP WRITE (LP_CONFIG,232) DX(I),DY(I),DZ(I) 655 CONTINUE 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(300.*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,BXX,BYY,BZZ 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) BXX = 0. BYY = 0. BZZ = 0. 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 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) 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 XX = I0*DERR ARGIN = PRODUCT-XX 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.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*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 BXX = BXX+FF*RX*RX BYY = BYY+FF*RY*RY BZZ = BZZ+FF*RZ*RZ 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 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) VIRIAL = VIRIAL-F1*DR(1,2)*DR(1,2) 1 -F2*DR(1,3)*DR(1,3)-F3*DR(2,3)*DR(2,3) BXX = BXX-F1*XR(1,2)*XR(1,2) 1 -F2*XR(1,3)*XR(1,3)-F3*XR(2,3)*XR(2,3) BYY = BYY-F1*YR(1,2)*YR(1,2) 1 -F2*YR(1,3)*YR(1,3)-F3*YR(2,3)*YR(2,3) BZZ = BZZ-F1*ZR(1,2)*ZR(1,2) 1 -F2*ZR(1,3)*ZR(1,3)-F3*ZR(2,3)*ZR(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 COMMON /OTHER/ AMASS,HPART,OPART,DX(NPD),DY(NPD),DZ(NPD) 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 IF (MOD(L,3).EQ.1) THEN XX=HPART*(SX(L+1)+SX(L+2)-2*SX(L)) YY=HPART*(SY(L+1)+SY(L+2)-2*SY(L)) ZZ=HPART*(SZ(L+1)+SZ(L+2)-2*SZ(L)) ENDIF DX(L+M) = XX DY(L+M) = YY DZ(L+M) = ZZ X(L+M) = SX(L)+(I+1./8.)*A Y(L+M) = SY(L)+(J+1./8.)*A Z(L+M) = SZ(L)+(K+1./8.)*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 -----------------------------------------------------------