Line data    Source code

       1             : MODULE m_ccdnup
       2             :   !     *******************************************************
       3             :   !     *****   set up the core densities for compounds.  *****
       4             :   !     *****   in accordanse to d.d.koelling's cored     *****
       5             :   !     *******************************************************
       6             : CONTAINS
       7           2 :   SUBROUTINE ccdnup(atoms,sphhar,input,jatom,rho,sume,vrs,rhochr,rhospn,tecs,qints)
       8             : 
       9             :     USE m_constants
      10             :     USE m_intgr, ONLY : intgr3
      11             :     USE m_types
      12             :     IMPLICIT NONE
      13             :     TYPE(t_input),INTENT(IN)   :: input
      14             :     TYPE(t_sphhar),INTENT(IN)  :: sphhar
      15             :     TYPE(t_atoms),INTENT(IN)   :: atoms
      16             :     !     ..
      17             :     !     .. Scalar Arguments ..
      18             :     INTEGER, INTENT (IN) :: jatom 
      19             :     REAL,    INTENT (IN) :: sume
      20             :     !     ..
      21             :     !     .. Array Arguments ..
      22             :     REAL,    INTENT (IN) :: rhochr(:),rhospn(:)!(atoms%msh)
      23             :     REAL,    INTENT (IN) :: vrs(:,:,:)!(atoms%jmtd,atoms%ntype,input%jspins)
      24             :     REAL,    INTENT (OUT) :: tecs(:,:)!(atoms%ntype,input%jspins)
      25             :     REAL,    INTENT (OUT) :: qints(:,:)!(atoms%ntype,input%jspins)
      26             :     REAL,    INTENT (INOUT) :: rho(:,0:,:,:)!(atoms%jmtd,0:sphhar%nlhd,atoms%ntype,input%jspins)
      27             :     !     ..
      28             :     !     .. Local Scalars ..
      29             :     REAL d,dxx,q,rad,rhs
      30             :     INTEGER i,j,jspin,nm,nm1
      31             :     !     ..
      32             :     !     .. Local Arrays ..
      33           2 :     REAL rhoc(SIZE(rhochr)),rhoss(SIZE(rhochr))
      34             :     !     ..
      35           2 :     nm = atoms%jri(jatom)
      36             :     !     ---->update spherical charge density rho with the core density.
      37             :     !     ---->for spin-polarized (jspins=2), take only half the density
      38           6 :     DO jspin = 1,input%jspins
      39           4 :        IF (input%jspins.EQ.2 .AND. jspin.EQ.1) THEN
      40        1348 :           DO j = 1,SIZE(rhochr)
      41        1348 :              rhoss(j) = rhochr(j) - rhospn(j)
      42             :           END DO
      43           2 :        ELSE IF (input%jspins.EQ.2 .AND. jspin.EQ.2) THEN
      44        1348 :           DO j = 1,SIZE(rhochr)
      45        1348 :              rhoss(j) = rhochr(j) + rhospn(j)
      46             :           END DO
      47             :           ! jspins=1
      48             :        ELSE
      49           0 :           DO j = 1,SIZE(rhochr)
      50           0 :              rhoss(j) = rhochr(j)
      51             :           END DO
      52             :           !
      53             :        END IF
      54             :        !
      55        2264 :        DO  j = 1,nm
      56        2260 :           rhoc(j) = rhoss(j)/input%jspins
      57        2264 :           rho(j,0,jatom,jspin) = rho(j,0,jatom,jspin) + rhoc(j)/sfp_const
      58             :        ENDDO
      59        2264 :        DO  i = 1,nm
      60        2264 :           rhoc(i) = rhoc(i)*vrs(i,jatom,jspin)/atoms%rmsh(i,jatom)
      61             :        ENDDO
      62           4 :        CALL intgr3(rhoc,atoms%rmsh(1,jatom),atoms%dx(jatom),nm,rhs)
      63           4 :        tecs(jatom,jspin) = sume/input%jspins - rhs
      64           4 :        WRITE (oUnit,FMT=8010) jatom,jspin,tecs(jatom,jspin),sume/input%jspins
      65             :   
      66             :        !     ---> simpson integration
      67           4 :        dxx = atoms%dx(jatom)
      68           4 :        d = EXP(atoms%dx(jatom))
      69           4 :        rad = atoms%rmt(jatom)
      70           4 :        q = rad*rhoss(nm)/2.
      71           4 :        DO  nm1 = nm + 1,SIZE(rhochr) - 1,2
      72         216 :           rad = d*rad
      73         216 :           q = q + 2*rad*rhoss(nm1)
      74         216 :           rad = d*rad
      75         216 :           q = q + rad*rhoss(nm1+1)
      76             :        ENDDO
      77           4 :        q = 2*q*dxx/3
      78           4 :        WRITE (oUnit,FMT=8000) q/input%jspins
      79          10 :        qints(jatom,jspin) = q*atoms%neq(jatom)
      80             : 
      81             :     END DO ! end-do-loop input%jspins
      82             : 
      83             : 8000 FORMAT (f20.8,' electrons lost from core.')
      84             : 8010 FORMAT (10x,'atom type',i5,'  (spin',i2,') ',/,10x,&
      85             :          &       'kinetic energy=',e20.12,5x,'sum of the eigenvalues=',&
      86             :          &       e20.12)
      87             : 
      88           2 :   END SUBROUTINE ccdnup
      89             : END MODULE m_ccdnup