cored.F90 9.51 KB
Newer Older
1 2
MODULE m_cored
CONTAINS
Matthias Redies's avatar
Matthias Redies committed
3
  SUBROUTINE cored(input, jspin, atoms, rho, DIMENSION, sphhar, vr, qint, rhc, tec, seig, EnergyDen)
4 5 6 7 8 9 10 11 12 13
    !     *******************************************************
    !     *****   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
14
    USE m_xmlOutput
15 16 17 18 19 20 21 22 23 24 25
    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 ..
Matthias Redies's avatar
Matthias Redies committed
26 27 28 29 30 31
    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)
    REAL, INTENT(INOUT), OPTIONAL :: EnergyDen(atoms%jmtd,0:sphhar%nlhd,atoms%ntype,DIMENSION%jspd)
32 33
    !     ..
    !     .. Local Scalars ..
Matthias Redies's avatar
Matthias Redies committed
34 35
    REAL eig,fj,fl,fn,q,rad,rhos,rhs,sea,sume,t2
    REAL d,dxx,rn,rnot,z,t1,rr,r,lambd,c,bmu,weight, aux_weight
Matthias Redies's avatar
Matthias Redies committed
36
    INTEGER i,j,jatom,korb,n,ncmsh,nm,nm1,nst ,l,ierr
37 38
    !     ..
    !     .. Local Arrays ..
39
    
40
    REAL rhcs(DIMENSION%msh),rhoc(DIMENSION%msh),rhoss(DIMENSION%msh),vrd(DIMENSION%msh),f(0:3)
Matthias Redies's avatar
Matthias Redies committed
41
    REAL rhcs_aux(DIMENSION%msh), rhoss_aux(DIMENSION%msh) !> quantities for energy density calculations
42 43 44
    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)
45 46
    CHARACTER(LEN=20) :: attributes(6)
    REAL stateEnergies(29)
47
    !     ..
Matthias Redies's avatar
Matthias Redies committed
48

49 50 51 52 53 54 55 56 57 58
    c = c_light(1.0)
    seig = 0.
    !
    IF (input%frcor) THEN
       DO  n = 1,atoms%ntype
          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)
59
             rhoc(i) = rhc(i,n,jspin)
60 61 62 63 64 65 66 67 68 69
             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)
70 71 72
          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
73 74 75 76 77 78 79 80 81 82 83 84
          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
85
       CALL setcor(jatom,DIMENSION%jspd,atoms,input,bmu,nst,kappa,nprnc,occ_h)
86 87 88 89 90 91 92 93 94 95 96 97 98
       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)
Matthias Redies's avatar
Matthias Redies committed
99
          rhoss(j)     = 0.0
100 101 102
          vrd(j) = vr(j,jatom)
       ENDDO
       !
103 104 105 106 107 108 109 110 111 112
       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
113
       IF ( atoms%jri(jatom) < ncmsh) THEN
114 115
          DO  i = atoms%jri(jatom) + 1,ncmsh
             rhoss(i) = 0.
116 117 118 119 120 121 122
             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
123 124 125 126 127
             !
          ENDDO
       END IF

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

129
       IF (input%gw==1 .OR. input%gw==3)&
130 131
            &                      WRITE(15) nst,atoms%rmsh(1:atoms%jri(jatom),jatom)

132
       stateEnergies = 0.0
133
       DO  korb = 1,nst
134 135 136
          IF (occ(korb) /= 0.0) THEN
            fn = nprnc(korb)
            fj = iabs(kappa(korb)) - .5e0
Matthias Redies's avatar
Matthias Redies committed
137

138 139
            weight = 2*fj + 1.e0
            IF (bmu > 99.) weight = occ(korb)
Matthias Redies's avatar
Matthias Redies committed
140

141
            fl = fj + (.5e0)*isign(1,kappa(korb))
Matthias Redies's avatar
Matthias Redies committed
142 143 144 145 146 147
            eig        = -2* (z/ (fn+fl))**2

            CALL differ(fn,fl,fj,c,z,dxx,rnot,rn,d,ncmsh,vrd, eig, a,b,ierr)
            stateEnergies(korb) = eig
            WRITE (6,FMT=8010) fn,fl,fj,eig,weight
            WRITE (16,FMT=8010) fn,fl,fj,eig,weight
148
            IF (ierr/=0)  CALL juDFT_error("error in core-level routine" ,calledby ="cored")
Matthias Redies's avatar
Matthias Redies committed
149
            IF (input%gw==1 .OR. input%gw==3) WRITE (15) NINT(fl),weight,eig,&
150
                a(1:atoms%jri(jatom)),b(1:atoms%jri(jatom))
151

Matthias Redies's avatar
Matthias Redies committed
152
            sume = sume + weight*eig/input%jspins
153
            DO j = 1,ncmsh
Matthias Redies's avatar
Matthias Redies committed
154
              rhcs(j)  = weight* (a(j)**2+b(j)**2)
155 156
              rhoss(j) = rhoss(j) + rhcs(j)
            ENDDO
Matthias Redies's avatar
Matthias Redies committed
157 158 159 160 161 162 163 164 165

            IF(present(EnergyDen)) THEN
              rhoss_aux = rhoss
              DO j = 1,ncmsh
                ! for energy density we want to multiply the weights 
                ! with the eigenenergies
                rhoss_aux(j) = rhoss_aux(j) + (rhcs(j) * eig)
              ENDDO
            ENDIF
166
          ENDIF
167 168 169 170 171 172 173 174 175 176
       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

Matthias Redies's avatar
Matthias Redies committed
177 178 179 180 181 182 183 184
       IF(present(EnergyDen)) then
        DO  j = 1,nm
            rhoc(j) = rhoss(j)/input%jspins
            EnergyDen(j,0,jatom,jspin) = EnergyDen(j,0,jatom,jspin) &
                                       + rhoss_aux(j) /(input%jspins * sfp_const)
        ENDDO
      ENDIF

185 186
       rhc(1:ncmsh,jatom,jspin)   = rhoss(1:ncmsh) / input%jspins
       rhc(ncmsh+1:DIMENSION%msh,jatom,jspin) = 0.0
187 188 189 190 191 192

       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)
193 194 195
       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
196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211

       !     ---> 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)
212 213 214 215
       attributes = ''
       WRITE(attributes(1),'(i0)') jatom
       WRITE(attributes(2),'(i0)') NINT(z)
       WRITE(attributes(3),'(i0)') jspin
216
       WRITE(attributes(4),'(f18.10)') tec(jatom,jspin)
217 218
       WRITE(attributes(5),'(f18.10)') sume
       WRITE(attributes(6),'(f9.6)') q/input%jspins
219
       CALL openXMLElementForm('coreStates',(/'atomType     ','atomicNumber ','spin         ','kinEnergy    ',&
220 221
            'eigValSum    ','lostElectrons'/),&
            attributes,RESHAPE((/8,12,4,9,9,13,6,3,1,18,18,9/),(/6,2/)))
222 223 224 225 226 227 228
       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)
229
          WRITE(attributes(2),'(i0)') NINT(fl)
230 231 232 233
          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'/),&
234
               attributes(1:5),RESHAPE((/1,1,1,6,6,2,2,4,20,15/),(/5,2/)))
235 236
       END DO
       CALL closeXMLElement('coreStates')
237 238 239 240
    ENDDO

    RETURN

241
8000 FORMAT (/,/,10x,'z=',f4.0,5x,'r(1)=',e14.6,5x,'dx=',f9.6,5x,&
242 243
         &       'm.t.index=',i4,/,15x,'n',4x,'l',5x,'j',4x,'energy',7x,&
         &       'weight')
Daniel Wortmann's avatar
Daniel Wortmann committed
244
8010 FORMAT (12x,2f5.0,f6.1,f10.4,f12.4)
245 246 247 248
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)
249
  END SUBROUTINE cored
250 251
END MODULE m_cored