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