C ////////////////////////////////////////////////////////////
C // EWALD SUMMATION FOR TOTAL ENERGY , FORCE (VEC VERSION) //
C // 7 OCTOBER, 1986 ---> 3/15,1990 //
C ////////////////////////////////////////////////////////////
SUBROUTINE EWVMD(ETOT2)
IMPLICIT REAL(A-H,O-Y)
IMPLICIT COMPLEX(Z)
INCLUDE 'PACVPP'
PARAMETER(KEWA=KATM*(KATM+1)/2)
DIMENSION CR(KLM),CI(KLM)
DIMENSION CRX(KLM),CRY(KLM),CRZ(KLM),CLX(KLM),CLY(KLM),CLZ(KLM)
DIMENSION SSS(KLM),SSI(KLM)
DIMENSION ITYP(KATM)
DIMENSION FORC1(KATM,KATM,3)
C
CON=4.D0*PAI/UNIVOL
CC=2.D0/SQRT(PAI)
ETOT2=0.D0
IRL=KLX2*KLY2*KLZ2
IKM=KMX2*KMY2*KMZ2
MM=0
DO 2 II=1,KTYP
I1=MM+1
I2=MM+IATOM(II)
MM=I2
DO 1 JJ=I1,I2
C <------ ATTN]]] KFTYPE <------- ]]]]]
ITYP(JJ)=II
1 CONTINUE
2 CONTINUE
IF(MM.NE.KATM)
& WRITE(6,*) '****WARNING "MM.NE.KATM" IN EWALD DO-LOOP 1****'
C
DO 110 J=1,KATM
DO 111 I=1,KATM
FORC1(I,J,1)=0.0D0
FORC1(I,J,2)=0.0D0
FORC1(I,J,3)=0.0D0
111 CONTINUE
110 CONTINUE
C
DO 100 I=1,KATM
DO 101 J=I,KATM
ITE=1
ALFF=ALF
77 CONTINUE
AL1=SQRT(ALFF)
AL2=0.25D0/ALFF
SUM =0.0D0
SF1X =0.0D0
SF1Y =0.0D0
SF1Z =0.0D0
SF1XI=0.0D0
SF1YI=0.0D0
SF1ZI=0.0D0
C
R1=CATX(I)-CATX(J)
R2=CATY(I)-CATY(J)
R3=CATZ(I)-CATZ(J)
C ///////////////////////////////////
C // STEP 1; REAL-SPACE SUMMATION //
C ///////////////////////////////////
DO 11 L=1,IRL
SSS(L)=0.0D0
SSI(L)=0.0D0
CRX(L)=CX(L)+R1
CRY(L)=CY(L)+R2
CRZ(L)=CZ(L)+R3
CLX(L)=CX(L)-R1
CLY(L)=CY(L)-R2
CLZ(L)=CZ(L)-R3
CR(L)=SQRT((CX(L)+R1)**2 + (CY(L)+R2)**2 + (CZ(L)+R3)**2 )
CI(L)=SQRT((CX(L)-R1)**2 + (CY(L)-R2)**2 + (CZ(L)-R3)**2 )
11 CONTINUE
C
DO 12 L=1,IRL
IF(CR(L).EQ.0.0D0) THEN
SUM=SUM-CC*AL1
ELSE
SUM=SUM+ (1.D0-ERF(AL1*CR(L)))/CR(L)
END IF
12 CONTINUE
C-----FORCE CALCULATION (R-SPACE) ---------------------------------
IF (I.NE.J) THEN
DO 14 L=1,IRL
IF (CR(L).NE.0.0D0)
& SSS(L)=(1.0D0-ERF(AL1*CR(L)))/(CR(L)**3)+(CC*AL1
& *EXP(-ALFF*CR(L)*CR(L)))/(CR(L)*CR(L))
IF (CI(L).NE.0.0D0)
& SSI(L)=(1.0D0-ERF(AL1*CI(L)))/(CI(L)**3)+(CC*AL1
& *EXP(-ALFF*CI(L)*CI(L)))/(CI(L)*CI(L))
14 CONTINUE
C
DO 15 L=1,IRL
SF1X =SF1X +CRX(L)*SSS(L)
SF1Y =SF1Y +CRY(L)*SSS(L)
SF1Z =SF1Z +CRZ(L)*SSS(L)
SF1XI=SF1XI+CLX(L)*SSI(L)
SF1YI=SF1YI+CLY(L)*SSI(L)
SF1ZI=SF1ZI+CLZ(L)*SSI(L)
15 CONTINUE
C FORCE CALCULATION (K-SPACE) FOR OPTICAL VIBRATION
DO 24 L=1,IKM
IF (CKR(L).EQ.0.0D0) GOTO 24
SSK =(CON*SIN(CKX(L)*R1+CKY(L)*R2+CKZ(L)*R3)
& *EXP(-AL2*CKR(L)*CKR(L)))/(CKR(L)*CKR(L))
SF1X =SF1X +CKX(L)*SSK
SF1Y =SF1Y +CKY(L)*SSK
SF1Z =SF1Z +CKZ(L)*SSK
SF1XI=SF1XI-CKX(L)*SSK
SF1YI=SF1YI-CKY(L)*SSK
SF1ZI=SF1ZI-CKZ(L)*SSK
24 CONTINUE
C
FORC1(I,J,1)=
& ACHG(ITYP(I))*ACHG(ITYP(J))*SF1X
FORC1(I,J,2)=
& ACHG(ITYP(I))*ACHG(ITYP(J))*SF1Y
FORC1(I,J,3)=
& ACHG(ITYP(I))*ACHG(ITYP(J))*SF1Z
FORC1(J,I,1)=
& ACHG(ITYP(I))*ACHG(ITYP(J))*SF1XI
FORC1(J,I,2)=
& ACHG(ITYP(I))*ACHG(ITYP(J))*SF1YI
FORC1(J,I,3)=
& ACHG(ITYP(I))*ACHG(ITYP(J))*SF1ZI
END IF
C /////////////////////////////////////////
C // STEP 2; RECIPROCAL-SPACE SUMMATION //
C /////////////////////////////////////////
DO 21 L=1,IKM
C CX = CY = CZ =0
IF (CKR(L).EQ.0.0D0) GOTO 21
SUM=SUM+ CON*EXP(-AL2*CKR(L)*CKR(L))
& *COS(CKX(L)*R1+CKY(L)*R2+CKZ(L)*R3)/(CKR(L)*CKR(L))
21 CONTINUE
SUM= SUM -PAI/ALFF/UNIVOL
C ////////////
C // STEP 3 //
C ////////////
IF(ITE.EQ.1) THEN
SUM1 =SUM
C------------SF2X =SF1X
C SF2Y =SF1Y
C SF2Z =SF1Z
C SF2XI=SF1XI
C SF2YI=SF1YI
C------------SF2ZI=SF1ZI
ALFF=2.D0*ALFF
ITE=2
ELSE
IF(ABS(SUM-SUM1).GT.1.D-5) WRITE(6,*)
& '******WARNING SUM.NE.SUM1***********'
END IF
C-----WRITE (6,22) I,J,SUM1 ,SUM
C WRITE (6,22) I,J,SF2X ,SF1X
C WRITE (6,22) I,J,SF2Y ,SF1Y
C WRITE (6,22) I,J,SF2Z ,SF1Z
C WRITE (6,22) I,J,SF2XI,SF1XI
C WRITE (6,22) I,J,SF2YI,SF1YI
C WRITE (6,22) I,J,SF2ZI,SF1ZI
FAC=1.D0
IF(I.EQ.J) FAC=0.5D0
ETOT2= ETOT2 + FAC*SUM*ACHG(ITYP(I))*ACHG(ITYP(J))
22 FORMAT(1H ,2(I5,1X),' EWALD-SUMM=',D16.8,2X,D16.8)
101 CONTINUE
100 CONTINUE
C
DO 150 I=1,KATM
FFF1(I,1)=0.0D0
FFF1(I,2)=0.0D0
FFF1(I,3)=0.0D0
150 CONTINUE
C
DO 200 I=1,KATM
DO 201 J=1,KATM
FFF1(I,1)=FFF1(I,1)+FORC1(I,J,1)
FFF1(I,2)=FFF1(I,2)+FORC1(I,J,2)
FFF1(I,3)=FFF1(I,3)+FORC1(I,J,3)
201 CONTINUE
200 CONTINUE
C
WRITE(6,*) 'ETOT2(EWALD) = ',ETOT2
RETURN
END