C---*----1----*----2----*----3----*----4----*----5----*----6----*----7
SUBROUTINE PCC(KTPCC,RHPC,
< SCHGPC,EPC,PCM)
C SUBROUTINE FOR PARTIAL CORE CORRECTION
C MODIFIED FROM ATOMIC SEP.1990 MORIKAWA
C REF. S.G.LOUIE,S.FROYEN,AND M.L.COHEN,
C PHYS.REV.B26,1738(1982)
IMPLICIT REAL(A-H,O-Y)
IMPLICIT COMPLEX(Z)
INCLUDE 'PACVPP'
INTEGER I,IA,ITY,N,KTPCC(KTYP)
REAL EPC,PCM,RS
REAL XH,H,H3,S,SCHGPC
REAL RHPC(KIZA,KTYP),CHGPC(KTYP),EPCP(KTYP),EMCP(KTYP)
DATA XH/32.0D0/
H = 1.D0/XH
H3 = H/3.D0
CO3=4.D0/3.D0
C WRITE (6,*) 'UNIVOL IN PCC =',UNIVOL
C///////////// FOURIER TRANSFORMATION OF PARTIAL CHARGE DENSITY /////
IF (ITER.EQ.0) THEN
DO 1110 ITY=1,KTYP
DO 1111 N=1,KIZA
RHPC(N,ITY)=0.0D0
1111 CONTINUE
1110 CONTINUE
C----------------------------------------
REWIND 27
READ(27,1100) (I,RAD(N),N=1,KIZA)
1100 FORMAT(3(I4,D22.14))
IF (I.NE.KIZA) THEN
WRITE (6,*) ' READ IN PCC '
STOP
END IF
WRITE(6,*) I,' RAD(KIZA)= ',RAD(KIZA)
DO 1000 ITY=1,KTYP
IF (KTPCC(ITY).EQ.0) GOTO 1000
READ(27,*) CHGPC(ITY)
READ(27,1100) (I,RHPC(N,ITY),N=1,KIZA)
IF (I.NE.KIZA) THEN
WRITE (6,*) ' READ IN PCC ITY= ',ITY
STOP
END IF
1000 CONTINUE
END IF
C--------------------------------
DO 1001 ITY=1,KTYP
WRITE(6,*) ITY,' RHPC(KIZA,ITY)= ',RHPC(KIZA,ITY)
S = 0.D0
EPCP(ITY) = 0.D0
EMCP(ITY) = 0.D0
DO 1150 N=1,KIZA
IF (RHPC(N,ITY).GE.1.D-40) THEN
RS = (3.D0/RHPC(N,ITY))**(1.D0/3.D0)*RAD(N)
EPCP(ITY) = EPCP(ITY)-(0.458D0/RS+0.44D0/(RS+7.8D0))
& *RHPC(N,ITY)*OMO(N)
EMCP(ITY) = EMCP(ITY)-(0.458D0*CO3/RS+(0.44D0*CO3*RS +
& 3.432D0)/(RS+7.8D0)**2)*RHPC(N,ITY)*OMO(N)
END IF
S = S + RHPC(N,ITY)*OMO(N)
1150 CONTINUE
S = S *H3
EPCP(ITY) = EPCP(ITY)*H3
EMCP(ITY) = EMCP(ITY)*H3
WRITE (6,*) ' CHGPC(ITY) = ',CHGPC(ITY)
IF (ABS(S-CHGPC(ITY)).GE.1.D-7) THEN
WRITE (6,*) ITY,'PARTIAL CORE CHARGE S-CHGPC='
/ ,S-CHGPC(ITY)
STOP
END IF
C------------------------------------
DO 1300 I=1,KG
RHPCG(I,ITY)=0.0D0
RRPCG(I,ITY)=0.0D0
1300 CONTINUE
DO 1200 N=1,KIZA
R=RAD(N)
FAC1=OMO(N)*RHPC(N,ITY)*H3/UNIVOL
FAC2=OMO(N)*RHPC(N,ITY)*R*H3/UNIVOL
DO 1210 I=1,KG
XX(I)=R*GR(I)
1210 CONTINUE
CALL DSJNV(0,KG,XX,Y11)
CALL DSJNV(1,KG,XX,Y22)
DO 1220 I=1,KG
RHPCG(I,ITY)=RHPCG(I,ITY)+FAC1*Y11(I)
RRPCG(I,ITY)=RRPCG(I,ITY)-FAC2*Y22(I)
1220 CONTINUE
1200 CONTINUE
1001 CONTINUE
SCHGPC = 0.D0
EPC = 0.D0
PCM = 0.0D0
DO 1410 IA = 1,KATM
IF (KTPCC(KFTYPE(IA)).EQ.0) GOTO 1410
SCHGPC = SCHGPC+CHGPC(KFTYPE(IA))
EPC = EPC +EPCP (KFTYPE(IA))
PCM = PCM +EMCP (KFTYPE(IA))
1410 CONTINUE
WRITE (6,*) 'NO-OVERLAP P.C.(USED) EPC = ',EPC
WRITE (6,*) 'NO-OVERLAP P.C.(USED) PCM = ',PCM
WRITE (6,*) 'TOTAL PARTIAL CHARGE SCHGPC = ',SCHGPC
C///////////////////////////////////////////////////////////////////
ZRHPC(1)=DCMPLX(0.D0,0.D0)
ZRRPC(1)=DCMPLX(0.D0,0.D0)
RG(1)=1.0D0
DO 1500 I=2,KG
ZRHPC(I)=DCMPLX(0.D0,0.D0)
ZRRPC(I)=DCMPLX(0.D0,0.D0)
RG(I)=1.0D0/GR(I)
1500 CONTINUE
DO 2000 IT=1,KTYP
IF (KTPCC(IT).EQ.0) GOTO 2000
DO 2100 I=1,KG
ZRHPC(I) = ZRHPC(I) + ZFM3(I,IT)*RHPCG(I,IT)
ZRRPC(I) = ZRRPC(I) + ZFM3(I,IT)*RRPCG(I,IT)*RG(I)
2100 CONTINUE
2000 CONTINUE
RETURN
END