C ////////////////////////////////////////////////////////////
C // EWALD SUMMATION FOR TOTAL ENERGY , FORCE (VEC VERSION) //
C // 7 OCTOBER, 1986 ---> 3/15,1990 //
C ////////////////////////////////////////////////////////////
SUBROUTINE EWVEC(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 ************************************************************
C ** THE RESULTS SHOULD BE INDEPENDENT OF THE VALUE OF ALF. **
C ************************************************************
A1X=ALTV(1,1)
A1Y=ALTV(2,1)
A1Z=ALTV(3,1)
A2X=ALTV(1,2)
A2Y=ALTV(2,2)
A2Z=ALTV(3,2)
A3X=ALTV(1,3)
A3Y=ALTV(2,3)
A3Z=ALTV(3,3)
B1X=RLTV(1,1)
B1Y=RLTV(2,1)
B1Z=RLTV(3,1)
B2X=RLTV(1,2)
B2Y=RLTV(2,2)
B2Z=RLTV(3,2)
B3X=RLTV(1,3)
B3Y=RLTV(2,3)
B3Z=RLTV(3,3)
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
MM=0
DO 10 I1=1,KLX2
F1=FLOAT(I1-KLX-1)
DO 10 I2=1,KLY2
F2=FLOAT(I2-KLY-1)
DO 10 I3=1,KLZ2
F3=FLOAT(I3-KLZ-1)
MM=MM+1
CX(MM)= F1*A1X + F2*A2X + F3*A3X
CY(MM)= F1*A1Y + F2*A2Y + F3*A3Y
CZ(MM)= F1*A1Z + F2*A2Z + F3*A3Z
10 CONTINUE
C
MK=0
DO 20 I1=1,KMX2
F1=FLOAT(I1-KMX-1)
DO 20 I2=1,KMY2
F2=FLOAT(I2-KMY-1)
DO 20 I3=1,KMZ2
F3=FLOAT(I3-KMZ-1)
MK=MK+1
CKX(MK)= F1*B1X + F2*B2X + F3*B3X
CKY(MK)= F1*B1Y + F2*B2Y + F3*B3Y
CKZ(MK)= F1*B1Z + F2*B2Z + F3*B3Z
CKR(MK)= SQRT( CKX(MK)**2 + CKY(MK)**2 + CKZ(MK)**2 )
20 CONTINUE
C
IF (IRL.NE.MM.OR.IKM.NE.MK)
& WRITE (6,*) 'WORNING] IRL .NE. MM OR IKM .NE. MK]'
C STRESS 1
DO 1000 IS=1,6
SIGEWA(IS)=0.0D0
1000 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
SIG1 =0.0D0
SIG2 =0.0D0
SIG3 =0.0D0
SIG4 =0.0D0
SIG5 =0.0D0
SIG6 =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,MM
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,MM
IF(CR(L).EQ.0.0D0) THEN
SUM=SUM-CC*AL1
ELSE
SUM=SUM+ (1.D0-ERF(AL1*CR(L)))/CR(L)
C-----STRESS CALCULATION
STMP = - (1.D0-ERF(AL1*CR(L)))/(CR(L)*CR(L)*CR(L))
& - CC*AL1*EXP(-ALFF*CR(L)*CR(L))/(CR(L)*CR(L))
SIG1=SIG1+STMP*CRX(L)*CRX(L)
SIG2=SIG2+STMP*CRX(L)*CRY(L)
SIG3=SIG3+STMP*CRX(L)*CRZ(L)
SIG4=SIG4+STMP*CRY(L)*CRY(L)
SIG5=SIG5+STMP*CRY(L)*CRZ(L)
SIG6=SIG6+STMP*CRZ(L)*CRZ(L)
END IF
12 CONTINUE
C-----FORCE CALCULATION (R & K-SPACE) ---------------------------------
IF (I.NE.J) THEN
DO 14 L=1,MM
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,MM
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,MK
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,MK
C CX = CY = CZ =0
IF (CKR(L).EQ.0.0D0) GOTO 21
STMP=CON*EXP(-AL2*CKR(L)*CKR(L))
& *COS(CKX(L)*R1+CKY(L)*R2+CKZ(L)*R3)/(CKR(L)*CKR(L))
SUM=SUM+ STMP
SIG1=SIG1+STMP*( 2.0D0*CKX(L)*CKX(L)*(
& 1.0D0/(CKR(L)*CKR(L)) + AL2 ) - 1.0D0 )
SIG2=SIG2+STMP*( 2.0D0*CKX(L)*CKY(L)*(
& 1.0D0/(CKR(L)*CKR(L)) + AL2 ) )
SIG3=SIG3+STMP*( 2.0D0*CKX(L)*CKZ(L)*(
& 1.0D0/(CKR(L)*CKR(L)) + AL2 ) )
SIG4=SIG4+STMP*( 2.0D0*CKY(L)*CKY(L)*(
& 1.0D0/(CKR(L)*CKR(L)) + AL2 ) - 1.0D0 )
SIG5=SIG5+STMP*( 2.0D0*CKY(L)*CKZ(L)*(
& 1.0D0/(CKR(L)*CKR(L)) + AL2 ) )
SIG6=SIG6+STMP*( 2.0D0*CKZ(L)*CKZ(L)*(
& 1.0D0/(CKR(L)*CKR(L)) + AL2 ) - 1.0D0 )
21 CONTINUE
SUM= SUM -PAI/ALFF/UNIVOL
SIG1=SIG1 +PAI/ALFF/UNIVOL
SIG4=SIG4 +PAI/ALFF/UNIVOL
SIG6=SIG6 +PAI/ALFF/UNIVOL
C ////////////
C // STEP 3 //
C ////////////
IF(ITE.EQ.1) THEN
SUM1 =SUM
SI21 =SIG1
C------------SF2X =SF1X
C SF2Y =SF1Y
C SF2Z =SF1Z
C SF2XI=SF1XI
C SF2YI=SF1YI
C------------SF2ZI=SF1ZI
ALFF=2.D0*ALFF
ITE=2
C
C IF (MOD(ITER,ILAM).EQ.1) GOTO 77
GOTO 77
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,SI21 ,SIG1
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))
SIGEWA(1)=SIGEWA(1)+FAC*SIG1*ACHG(ITYP(I))*ACHG(ITYP(J))
SIGEWA(2)=SIGEWA(2)+FAC*SIG2*ACHG(ITYP(I))*ACHG(ITYP(J))
SIGEWA(3)=SIGEWA(3)+FAC*SIG3*ACHG(ITYP(I))*ACHG(ITYP(J))
SIGEWA(4)=SIGEWA(4)+FAC*SIG4*ACHG(ITYP(I))*ACHG(ITYP(J))
SIGEWA(5)=SIGEWA(5)+FAC*SIG5*ACHG(ITYP(I))*ACHG(ITYP(J))
SIGEWA(6)=SIGEWA(6)+FAC*SIG6*ACHG(ITYP(I))*ACHG(ITYP(J))
22 FORMAT(1H ,2(I5,1X),' EWALD-SUMM=',D16.8,2X,D16.8)
101 CONTINUE
100 CONTINUE
SIGEWA(1)=SIGEWA(1)/UNIVOL
SIGEWA(2)=SIGEWA(2)/UNIVOL
SIGEWA(3)=SIGEWA(3)/UNIVOL
SIGEWA(4)=SIGEWA(4)/UNIVOL
SIGEWA(5)=SIGEWA(5)/UNIVOL
SIGEWA(6)=SIGEWA(6)/UNIVOL
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
IF (ITER.EQ.0) THEN
WRITE(6,*) 'TOTAL ENERGY CONTRIBUTION FROM THE EWALD SUMM.....'
WRITE(6,*) ETOT2
END IF
RETURN
END