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 ---------------------------------------------------