cored.F90 8.79 KB
Newer Older
1 2 3 4 5 6 7
MODULE m_cored
CONTAINS
  SUBROUTINE cored(&
       &                 input,jspin,atoms,&
       &                 rho,DIMENSION,&
       &                 sphhar,&
       &                 vr,&
8
       &                 qint,rhc,tec,seig)
9 10 11 12 13 14 15 16 17 18 19

    !     *******************************************************
    !     *****   set up the core densities for compounds.  *****
    !     *****                      d.d.koelling           *****
    !     *******************************************************
    USE m_juDFT
    USE m_intgr, ONLY : intgr3,intgr0,intgr1
    USE m_constants, ONLY : c_light,sfp_const
    USE m_setcor
    USE m_differ
    USE m_types
20
    USE m_xmlOutput
21 22 23 24 25 26 27 28 29 30 31
    IMPLICIT NONE
    TYPE(t_dimension),INTENT(IN)   :: DIMENSION
    TYPE(t_input),INTENT(IN)       :: input
    TYPE(t_sphhar),INTENT(IN)      :: sphhar
    TYPE(t_atoms),INTENT(IN)       :: atoms
    !
    !     .. Scalar Arguments ..
    INTEGER, INTENT (IN) :: jspin   
    REAL,    INTENT (OUT) :: seig
    !     ..
    !     .. Array Arguments ..
32 33 34 35 36
    REAL,    INTENT(IN)    :: vr(atoms%jmtd,atoms%ntype)
    REAL,    INTENT(INOUT) :: rho(atoms%jmtd,0:sphhar%nlhd,atoms%ntype,DIMENSION%jspd)
    REAL,    INTENT(INOUT) :: rhc(DIMENSION%msh,atoms%ntype,DIMENSION%jspd)
    REAL,    INTENT(INOUT) :: qint(atoms%ntype,DIMENSION%jspd)
    REAL,    INTENT(INOUT) :: tec(atoms%ntype,DIMENSION%jspd)
37 38
    !     ..
    !     .. Local Scalars ..
39
    REAL e,fj,fl,fn,q,rad,rhos,rhs,sea,sume,t2
40 41 42 43
    REAL d,dxx,rn,rnot,z,t1,rr,r,lambd,c,bmu,weight
    INTEGER i,j,jatom,jm,korb,n,ncmsh,nm,nm1,nst ,l,ierr
    !     ..
    !     .. Local Arrays ..
44
    
45 46 47 48
    REAL rhcs(DIMENSION%msh),rhoc(DIMENSION%msh),rhoss(DIMENSION%msh),vrd(DIMENSION%msh),f(0:3)
    REAL occ(DIMENSION%nstd),a(DIMENSION%msh),b(DIMENSION%msh),ain(DIMENSION%msh),ahelp(DIMENSION%msh)
    REAL occ_h(DIMENSION%nstd,2)
    INTEGER kappa(DIMENSION%nstd),nprnc(DIMENSION%nstd)
49 50
    CHARACTER(LEN=20) :: attributes(6)
    REAL stateEnergies(29)
51 52 53 54 55 56 57 58 59 60 61 62
    !     ..
    c = c_light(1.0)
    seig = 0.
    !
    IF (input%frcor) THEN
       DO  n = 1,atoms%ntype
          jm = atoms%jri(n)
          rnot = atoms%rmsh(1,n) ; dxx = atoms%dx(n)
          ncmsh = NINT( LOG( (atoms%rmt(n)+10.0)/rnot ) / dxx + 1 )
          ncmsh = MIN( ncmsh, DIMENSION%msh )
          !     --->    update spherical charge density
          DO  i = 1,atoms%jri(n)
63
             rhoc(i) = rhc(i,n,jspin)
64 65 66 67 68 69 70 71 72 73
             rho(i,0,n,jspin) = rho(i,0,n,jspin) + rhoc(i)/sfp_const
          ENDDO
          !     ---> for total energy calculations, determine the sum of the
          !     ---> eigenvalues by requiring that the core kinetic energy
          !     ---> remains constant.
          DO  i = 1,atoms%jri(n)
             rhoc(i) = rhoc(i)*vr(i,n)/atoms%rmsh(i,n)
          ENDDO
          nm = atoms%jri(n)
          CALL intgr3(rhoc,atoms%rmsh(1,n),atoms%dx(n),nm,rhos)
74 75 76
          sea = tec(n,jspin) + rhos
          WRITE (16,FMT=8030) n,jspin,tec(n,jspin),sea
          WRITE (6,FMT=8030) n,jspin,tec(n,jspin),sea
77 78 79 80 81 82 83 84 85 86 87 88
          seig = seig + atoms%neq(n)*sea
       ENDDO
       RETURN
    END IF

    !     ---> set up densities
    DO  jatom = 1,atoms%ntype
       sume = 0.
       z = atoms%zatom(jatom)
       !         rn = rmt(jatom)
       dxx = atoms%dx(jatom)
       bmu = 0.0
89
       CALL setcor(jatom,DIMENSION%jspd,atoms,input,bmu,nst,kappa,nprnc,occ_h)
90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106
       IF ((bmu > 99.)) THEN
          occ(1:nst) = input%jspins *  occ_h(1:nst,jspin)
       ELSE
          occ(1:nst) = occ_h(1:nst,1)
       ENDIF
       rnot = atoms%rmsh(1,jatom)
       d = EXP(atoms%dx(jatom))
       ncmsh = NINT( LOG( (atoms%rmt(jatom)+10.0)/rnot ) / dxx + 1 )
       ncmsh = MIN( ncmsh, DIMENSION%msh )
       rn = rnot* (d** (ncmsh-1))
       WRITE (6,FMT=8000) z,rnot,dxx,atoms%jri(jatom)
       WRITE (16,FMT=8000) z,rnot,dxx,atoms%jri(jatom)
       DO  j = 1,atoms%jri(jatom)
          rhoss(j) = 0.
          vrd(j) = vr(j,jatom)
       ENDDO
       !
107 108 109 110 111 112 113 114 115 116
       IF (input%l_core_confpot) THEN
          !--->    linear extension of the potential with slope t1 / a.u.
          t1=0.125
          t1 = MAX( (vrd(atoms%jri(jatom)) - vrd(atoms%jri(jatom)-1)*d)*&
               d / (atoms%rmt(jatom)**2 * (d-1) ) , t1)
          t2=vrd(atoms%jri(jatom))/atoms%rmt(jatom)-atoms%rmt(jatom)*t1
          rr = atoms%rmt(jatom)
       ELSE
          t2 = vrd(atoms%jri(jatom)) / ( atoms%jri(jatom) - ncmsh )
       ENDIF
117
       IF ( atoms%jri(jatom) < ncmsh) THEN
118 119
          DO  i = atoms%jri(jatom) + 1,ncmsh
             rhoss(i) = 0.
120 121 122 123 124 125 126
             IF (input%l_core_confpot) THEN
                rr = d*rr
                vrd(i) = rr*( t2 + rr*t1 )
                !               vrd(i) = 2*vrd(jri(jatom)) - rr*( t2 + rr*t1 )
             ELSE
                vrd(i) = vrd(atoms%jri(jatom)) + t2* (i-atoms%jri(jatom))
             ENDIF
127 128 129 130 131
             !
          ENDDO
       END IF

       nst = atoms%ncst(jatom)        ! for lda+U
132

133
       IF (input%gw==1 .OR. input%gw==3)&
134 135
            &                      WRITE(15) nst,atoms%rmsh(1:atoms%jri(jatom),jatom)

136
       stateEnergies = 0.0
137
       DO  korb = 1,nst
138 139 140 141 142 143 144 145 146 147 148 149 150 151
          IF (occ(korb) /= 0.0) THEN
            fn = nprnc(korb)
            fj = iabs(kappa(korb)) - .5e0
            weight = 2*fj + 1.e0
            IF (bmu > 99.) weight = occ(korb)
            fl = fj + (.5e0)*isign(1,kappa(korb))
            e = -2* (z/ (fn+fl))**2
            CALL differ(fn,fl,fj,c,z,dxx,rnot,rn,d,ncmsh,vrd, e, a,b,ierr)
            stateEnergies(korb) = e
            WRITE (6,FMT=8010) fn,fl,fj,e,weight
            WRITE (16,FMT=8010) fn,fl,fj,e,weight
            IF (ierr/=0)  CALL juDFT_error("error in core-level routine" ,calledby ="cored")
            IF (input%gw==1 .OR. input%gw==3) WRITE (15) NINT(fl),weight,e,&
                a(1:atoms%jri(jatom)),b(1:atoms%jri(jatom))
152

153 154 155 156 157 158
            sume = sume + weight*e/input%jspins
            DO j = 1,ncmsh
              rhcs(j) = weight* (a(j)**2+b(j)**2)
              rhoss(j) = rhoss(j) + rhcs(j)
            ENDDO
          ENDIF
159 160 161 162 163 164 165 166 167 168
       ENDDO

       !     ---->update spherical charge density rho with the core density.
       !     ---->for spin-polarized (jspins=2), take only half the density
       nm = atoms%jri(jatom)
       DO  j = 1,nm
          rhoc(j) = rhoss(j)/input%jspins
          rho(j,0,jatom,jspin) = rho(j,0,jatom,jspin) + rhoc(j)/sfp_const
       ENDDO

169 170
       rhc(1:ncmsh,jatom,jspin)   = rhoss(1:ncmsh) / input%jspins
       rhc(ncmsh+1:DIMENSION%msh,jatom,jspin) = 0.0
171 172 173 174 175 176

       seig = seig + atoms%neq(jatom)*sume
       DO  i = 1,nm
          rhoc(i) = rhoc(i)*vr(i,jatom)/atoms%rmsh(i,jatom)
       ENDDO
       CALL intgr3(rhoc,atoms%rmsh(1,jatom),atoms%dx(jatom),nm,rhs)
177 178 179
       tec(jatom,jspin) = sume - rhs
       WRITE (6,FMT=8030) jatom,jspin,tec(jatom,jspin),sume
       WRITE (16,FMT=8030) jatom,jspin,tec(jatom,jspin),sume
180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195

       !     ---> simpson integration
       rad = atoms%rmt(jatom)
       q = rad*rhoss(nm)/2.
       DO  nm1 = nm + 1,ncmsh - 1,2
          rad = d*rad
          q = q + 2*rad*rhoss(nm1)
          rad = d*rad
          q = q + rad*rhoss(nm1+1)
       ENDDO
       q = 2*q*dxx/3
       !+sb
       WRITE (6,FMT=8020) q/input%jspins
       WRITE (16,FMT=8020) q/input%jspins
       !-sb
       qint(jatom,jspin) = q*atoms%neq(jatom)
196 197 198 199
       attributes = ''
       WRITE(attributes(1),'(i0)') jatom
       WRITE(attributes(2),'(i0)') NINT(z)
       WRITE(attributes(3),'(i0)') jspin
200
       WRITE(attributes(4),'(f18.10)') tec(jatom,jspin)
201 202
       WRITE(attributes(5),'(f18.10)') sume
       WRITE(attributes(6),'(f9.6)') q/input%jspins
203
       CALL openXMLElementForm('coreStates',(/'atomType     ','atomicNumber ','spin         ','kinEnergy    ',&
204 205
            'eigValSum    ','lostElectrons'/),&
            attributes,RESHAPE((/8,12,4,9,9,13,6,3,1,18,18,9/),(/6,2/)))
206 207 208 209 210 211 212
       DO korb = 1, atoms%ncst(jatom)
          fj = iabs(kappa(korb)) - .5e0
          weight = 2*fj + 1.e0
          IF (bmu > 99.) weight = occ(korb)
          fl = fj + (.5e0)*isign(1,kappa(korb))
          attributes = ''
          WRITE(attributes(1),'(i0)') nprnc(korb)
213
          WRITE(attributes(2),'(i0)') NINT(fl)
214 215 216 217
          WRITE(attributes(3),'(f4.1)') fj
          WRITE(attributes(4),'(f20.10)') stateEnergies(korb)
          WRITE(attributes(5),'(f15.10)') weight
          CALL writeXMLElementForm('state',(/'n     ','l     ','j     ','energy','weight'/),&
218
               attributes(1:5),RESHAPE((/1,1,1,6,6,2,2,4,20,15/),(/5,2/)))
219 220
       END DO
       CALL closeXMLElement('coreStates')
221 222 223 224
    ENDDO

    RETURN

225
8000 FORMAT (/,/,10x,'z=',f4.0,5x,'r(1)=',e14.6,5x,'dx=',f9.6,5x,&
226 227
         &       'm.t.index=',i4,/,15x,'n',4x,'l',5x,'j',4x,'energy',7x,&
         &       'weight')
Daniel Wortmann's avatar
Daniel Wortmann committed
228
8010 FORMAT (12x,2f5.0,f6.1,f10.4,f12.4)
229 230 231 232
8020 FORMAT (f20.8,'  electrons lost from core.')
8030 FORMAT (10x,'atom type',i3,'  (spin',i2,') ',/,10x,&
         &       'kinetic energy=',e20.12,5x,'sum of the eigenvalues=',&
         &       e20.12)
233
  END SUBROUTINE cored
234 235
END MODULE m_cored