C ******************************************************************* C ** FICHE F.22 ** C ** REAL-SPACE AND RECIPROCAL-SPACE PARTS OF EWALD SUM FOR IONS. ** C ******************************************************************* C ******************************************************************* C ** REAL-SPACE AND RECIPROCAL-SPACE PARTS OF EWALD SUM FOR IONS. ** C ** ** C ** REFERENCES: ** C ** ** C ** WOODCOCK AND SINGER, TRANS. FARADAY SOC. 67, 12, 1971. ** C ** DE LEEUW ET AL., PROC. ROY. SOC. A 373, 27, 1980. ** C ** HEYES, J. CHEM. PHYS. 74, 1924, 1981. ** C ** SEE ALSO FINCHAM, MDIONS, CCP5 PROGRAM LIBRARY. ** C ** ** C ** ROUTINES SUPPLIED: ** C ** ** C ** SUBROUTINE SETUP ( KAPPA ) ** C ** SETS UP THE WAVEVECTORS FOR USE IN THE EWALD SUM ** C ** SUBROUTINE RWALD ( KAPPA, VR ) ** C ** CALCULATES THE R-SPACE PART OF THE SUM ** C ** SUBROUTINE KWALD ( KAPPA, VK ) ** C ** CALCULATES THE K-SPACE PART OF THE SUM ** C ** REAL FUNCTION ERFC ( X ) ** C ** RETURNS THE COMPLEMENTARY ERROR FUNCTION ** C ** ** C ** PRINCIPAL VARIABLES: ** C ** ** C ** INTEGER TOTK THE TOTAL NUMBER OF K-VECTORS STORED ** C ** INTEGER MAXK MAXIMUM POSSIBLE NUMBER OF K-VECTORS ** C ** INTEGER KMAX MAX INTEGER COMPONENT OF THE K-VECTOR ** C ** INTEGER KSQMAX MAX SQUARE MOD OF THE K-VECTOR REQUIRED ** C ** REAL VR ENERGY FROM R-SPACE SUM ** C ** REAL VK ENERGY FROM K-SPACE SUM ** C ** REAL KVEC(MAXK) ARRAY USED TO STORE K-VECTORS ** C ** REAL KAPPA WIDTH OF CANCELLING DISTRIBUTION ** C ** ** C ** USAGE: ** C ** ** C ** SETUP IS CALLED ONCE AT THE BEGINNING OF THE SIMULATION ** C ** TO CALCULATE ALL THE K-VECTORS REQUIRED IN THE EWALD SUM. ** C ** THESE VECTORS ARE USED THROUGHOUT THE SIMULATION IN THE ** C ** SUBROUTINE KWALD TO CALCULATE THE K-SPACE CONTRIBUTION TO THE ** C ** POTENTIAL ENERGY AT EACH CONFIGURATION. THE SELF TERM IS ** C ** SUBTRACTED FROM THE K-SPACE CONTRIBUTION IN KWALD. ** C ** THE SURFACE TERM FOR SIMULATIONS IN VACUUM IS NOT INCLUDED. ** C ** ROUTINE RWALD RETURNS THE R-SPACE CONTRIBUTION TO THE EWALD ** C ** SUM AND IS CALLED FOR EACH CONFIGURATION IN THE SIMULATION. ** C ** A CUBIC BOX AND UNIT BOX LENGTH ARE ASSUMED THROUGHOUT. ** C ******************************************************************* SUBROUTINE SETUP ( KAPPA ) COMMON / BLOCK2 / KVEC C ******************************************************************* C ** ROUTINE TO SET UP THE WAVE-VECTORS FOR THE EWALD SUM. ** C ** ** C ** THE WAVEVECTORS MUST FIT INTO A BOX OF UNIT LENGTH. ** C ** IN THIS EXAMPLE WE ALLOW A MAXIMUM OF 1000 WAVEVECTORS. ** C ******************************************************************* INTEGER MAXK PARAMETER ( MAXK = 1000 ) REAL KVEC(MAXK), KAPPA INTEGER KMAX, KSQMAX, KSQ, KX, KY, KZ, TOTK REAL TWOPI, B, RKX, RKY, RKZ, RKSQ PARAMETER ( KMAX = 5, KSQMAX = 27 , TWOPI = 6.2831853 ) C ******************************************************************* B = 1.0 / 4.0 / KAPPA / KAPPA C ** LOOP OVER K-VECTORS. NOTE KX IS NON-NEGATIVE ** TOTK = 0 DO 100 KX = 0, KMAX RKX = TWOPI * REAL ( KX ) DO 99 KY = -KMAX, KMAX RKY = TWOPI * REAL ( KY ) DO 98 KZ = -KMAX, KMAX RKZ = TWOPI * REAL ( KZ ) KSQ = KX * KX + KY * KY + KZ * KZ IF ( ( KSQ .LT. KSQMAX ) .AND. ( KSQ .NE. 0 ) ) THEN TOTK = TOTK + 1 IF ( TOTK .GT. MAXK ) STOP 'KVEC IS TOO SMALL' RKSQ = RKX * RKX + RKY * RKY + RKZ * RKZ KVEC(TOTK) = TWOPI * EXP ( -B * RKSQ ) / RKSQ ENDIF 98 CONTINUE 99 CONTINUE 100 CONTINUE WRITE( *, ' ( '' EWALD SUM SETUP COMPLETE '' ) ' ) WRITE( *, ' ( '' NUMBER OF WAVEVECTORS IS '', I5 ) ' ) TOTK RETURN END SUBROUTINE RWALD ( KAPPA, VR ) COMMON / BLOCK1 / RX, RY, RZ, Z C ******************************************************************* C ** CALCULATES R-SPACE PART OF POTENTIAL ENERGY BY EWALD METHOD. ** C ** ** C ** PRINCIPAL VARIABLES: ** C ** ** C ** INTEGER N NUMBER OF IONS ** C ** REAL RX(N),RY(N),RZ(N) POSITIONS OF IONS ** C ** REAL Z(N) IONIC CHARGES ** C ** REAL VR R-SPACE POTENTIAL ENERGY ** C ** ** C ** ROUTINE REFERENCED: ** C ** ** C ** REAL FUNCTION ERFC ( X ) ** C ** RETURNS THE COMPLEMENTARY ERROR FUNCTION ** C ******************************************************************* INTEGER N PARAMETER ( N = 216 ) REAL RX(N), RY(N), RZ(N), Z(N) REAL KAPPA, VR REAL RXI, RYI, RZI, ZI, RXIJ, RYIJ, RZIJ REAL RIJSQ, RIJ, KRIJ, ERFC, VIJ INTEGER I, J C ******************************************************************* VR = 0.0 DO 100 I = 1, N - 1 RXI = RX(I) RYI = RY(I) RZI = RZ(I) ZI = Z(I) DO 99 J = I + 1, N RXIJ = RXI - RX(J) RYIJ = RYI - RY(J) RZIJ = RZI - RZ(J) RXIJ = RXIJ - ANINT ( RXIJ ) RYIJ = RYIJ - ANINT ( RYIJ ) RZIJ = RZIJ - ANINT ( RZIJ ) RIJSQ = RXIJ * RXIJ + RYIJ * RYIJ + RZIJ * RZIJ RIJ = SQRT ( RIJSQ ) KRIJ = KAPPA * RIJ VIJ = ZI * Z(J) * ERFC ( KRIJ ) / RIJ VR = VR + VIJ 99 CONTINUE 100 CONTINUE RETURN END SUBROUTINE KWALD ( KAPPA, VK ) COMMON / BLOCK1 / RX, RY, RZ, Z COMMON / BLOCK2 / KVEC C ******************************************************************* C ** CALCULATES K-SPACE PART OF POTENTIAL ENERGY BY EWALD METHOD. ** C ** ** C ** THE SELF TERM IS SUBTRACTED. ** C ** IN ONE COORDINATE DIRECTION (X), SYMMETRY IS USED TO REDUCE ** C ** THE SUM TO INCLUDE ONLY POSITIVE K-VECTORS. ** C ** THE NEGATIVE VECTORS IN THIS DIRECTION ARE INCLUDED BY USE ** C ** OF THE MULTIPLICATIVE VARIABLE 'FACTOR'. ** C ** ** C ** PRINCIPAL VARIABLES: ** C ** ** C ** INTEGER N NUMBER OF IONS ** C ** REAL RX(N),RY(N),RZ(N) POSITIONS OF IONS ** C ** REAL Z(N) IONIC CHARGES ** C ** REAL VK K-SPACE POTENTIAL ENERGY ** C ** REAL VKS SELF PART OF K-SPACE SUM ** C ******************************************************************* INTEGER MAXK, N PARAMETER ( MAXK = 1000, N = 216 ) REAL KVEC(MAXK), RX(N), RY(N), RZ(N), Z(N) REAL KAPPA, VK INTEGER TOTK INTEGER KMAX, KX, KY, KZ, I, KSQMAX, KSQ REAL TWOPI, FACTOR, VD, VS, RSQPI PARAMETER ( KMAX = 5, KSQMAX = 27 ) PARAMETER ( TWOPI = 6.2831853, RSQPI = 0.5641896 ) COMPLEX EIKX(1:N, 0:KMAX) COMPLEX EIKY(1:N, -KMAX:KMAX) COMPLEX EIKZ(1:N, -KMAX:KMAX) COMPLEX EIKR(N), SUM C ******************************************************************* C ** CONSTRUCT EXP(IK.R) FOR ALL IONS AND K-VECTORS ** C ** CALCULATE KX, KY, KZ = 0 , -1 AND 1 EXPLICITLY ** DO 10 I = 1, N EIKX(I, 0) = (1.0, 0.0) EIKY(I, 0) = (1.0, 0.0) EIKZ(I, 0) = (1.0, 0.0) EIKX(I, 1) = CMPLX ( COS ( TWOPI * RX(I) ) , : SIN ( TWOPI * RX(I) ) ) EIKY(I, 1) = CMPLX ( COS ( TWOPI * RY(I) ) , : SIN ( TWOPI * RY(I) ) ) EIKZ(I, 1) = CMPLX ( COS ( TWOPI * RZ(I) ) , : SIN ( TWOPI * RZ(I) ) ) EIKY(I, -1) = CONJG ( EIKY(I, 1) ) EIKZ(I, -1) = CONJG ( EIKZ(I, 1) ) 10 CONTINUE C ** CALCULATE REMAINING KX, KY AND KZ BY RECURRENCE ** DO 12 KX = 2, KMAX DO 11 I = 1, N EIKX(I, KX) = EIKX(I, KX-1) * EIKX(I, 1) 11 CONTINUE 12 CONTINUE DO 14 KY = 2, KMAX DO 13 I = 1, N EIKY(I, KY) = EIKY(I, KY-1) * EIKY(I, 1) EIKY(I, -KY) = CONJG ( EIKY(I, KY) ) 13 CONTINUE 14 CONTINUE DO 16 KZ = 2, KMAX DO 15 I = 1, N EIKZ(I, KZ) = EIKZ(I, KZ-1) * EIKZ(I, 1) EIKZ(I, -KZ) = CONJG ( EIKZ(I, KZ) ) 15 CONTINUE 16 CONTINUE C ** SUM OVER ALL VECTORS ** VD = 0.0 TOTK = 0 DO 24 KX = 0, KMAX IF ( KX .EQ. 0 ) THEN FACTOR = 1.0 ELSE FACTOR = 2.0 ENDIF DO 23 KY = -KMAX, KMAX DO 22 KZ = -KMAX, KMAX KSQ = KX * KX + KY * KY + KZ * KZ IF ( ( KSQ .LT. KSQMAX ) .AND. ( KSQ .NE. 0 ) ) THEN TOTK = TOTK + 1 SUM = (0.0, 0.0) DO 21 I = 1, N EIKR(I) = EIKX(I, KX) * EIKY(I, KY) * EIKZ(I, KZ) SUM = SUM + Z(I) * EIKR(I) 21 CONTINUE VD = VD + FACTOR * KVEC(TOTK) * CONJG ( SUM ) * SUM ENDIF 22 CONTINUE 23 CONTINUE 24 CONTINUE C ** CALCULATES SELF PART OF K-SPACE SUM ** VS = 0.0 DO 25 I = 1, N VS = VS + Z(I) * Z(I) 25 CONTINUE VS = RSQPI * KAPPA * VS C ** CALCULATE THE TOTAL K-SPACE POTENTIAL ** VK = VD - VS RETURN END REAL FUNCTION ERFC ( X ) C ******************************************************************* C ** APPROXIMATION TO THE COMPLEMENTARY ERROR FUNCTION ** C ** ** C ** REFERENCE: ** C ** ** C ** ABRAMOWITZ AND STEGUN, HANDBOOK OF MATHEMATICAL FUNCTIONS, ** C ** NATIONAL BUREAU OF STANDARDS, FORMULA 7.1.26 ** C ******************************************************************* REAL A1, A2, A3, A4, A5, P PARAMETER ( A1 = 0.254829592, A2 = -0.284496736 ) PARAMETER ( A3 = 1.421413741, A4 = -1.453152027 ) PARAMETER ( A5 = 1.061405429, P = 0.3275911 ) REAL T, X, XSQ, TP C ******************************************************************* T = 1.0 / ( 1.0 + P * X ) XSQ = X * X TP = T * ( A1 + T * ( A2 + T * ( A3 + T * ( A4 + T * A5 ) ) ) ) ERFC = TP * EXP ( -XSQ ) RETURN END