LCOV - code coverage report
Current view: top level - core - cored.F90 (source / functions) Hit Total Coverage
Test: combined.info Lines: 105 128 82.0 %
Date: 2019-09-08 04:53:50 Functions: 1 1 100.0 %

          Line data    Source code
       1             : MODULE m_cored
       2             : CONTAINS
       3         295 :    SUBROUTINE cored(input, jspin, atoms, rho, DIMENSION, sphhar, vr, qint, rhc, tec, seig, EnergyDen)
       4             :       !     *******************************************************
       5             :       !     *****   set up the core densities for compounds.  *****
       6             :       !     *****                      d.d.koelling           *****
       7             :       !     *******************************************************
       8             :       USE m_juDFT
       9             :       USE m_intgr, ONLY : intgr3,intgr0,intgr1
      10             :       USE m_constants, ONLY : c_light,sfp_const
      11             :       USE m_setcor
      12             :       USE m_differ
      13             :       USE m_types
      14             :       USE m_xmlOutput
      15             :       IMPLICIT NONE
      16             :       TYPE(t_dimension),INTENT(IN)   :: DIMENSION
      17             :       TYPE(t_input),INTENT(IN)       :: input
      18             :       TYPE(t_sphhar),INTENT(IN)      :: sphhar
      19             :       TYPE(t_atoms),INTENT(IN)       :: atoms
      20             :       !
      21             :       !     .. Scalar Arguments ..
      22             :       INTEGER, INTENT (IN) :: jspin
      23             :       REAL,    INTENT (OUT) :: seig
      24             :       !     ..
      25             :       !     .. Array Arguments ..
      26             :       REAL, INTENT(IN)              :: vr(atoms%jmtd,atoms%ntype)
      27             :       REAL, INTENT(INOUT)           :: rho(atoms%jmtd,0:sphhar%nlhd,atoms%ntype,input%jspins)
      28             :       REAL, INTENT(INOUT)           :: rhc(DIMENSION%msh,atoms%ntype,input%jspins)
      29             :       REAL, INTENT(INOUT)           :: qint(atoms%ntype,input%jspins)
      30             :       REAL, INTENT(INOUT)           :: tec(atoms%ntype,input%jspins)
      31             :       REAL, INTENT(INOUT), OPTIONAL :: EnergyDen(atoms%jmtd,0:sphhar%nlhd,atoms%ntype,input%jspins)
      32             :       !     ..
      33             :       !     .. Local Scalars ..
      34             :       REAL eig,fj,fl,fn,q,rad,rhos,rhs,sea,sume,t2
      35             :       REAL d,dxx,rn,rnot,z,t1,rr,r,lambd,c,bmu,weight, aux_weight
      36             :       INTEGER i,j,jatom,korb,n,ncmsh,nm,nm1,nst ,l,ierr
      37             :       !     ..
      38             :       !     .. Local Arrays ..
      39             : 
      40        1180 :       REAL rhcs(DIMENSION%msh),rhoc(DIMENSION%msh),rhoss(DIMENSION%msh),vrd(DIMENSION%msh),f(0:3)
      41         590 :       REAL rhcs_aux(DIMENSION%msh), rhoss_aux(DIMENSION%msh) !> quantities for energy density calculations
      42        1180 :       REAL occ(DIMENSION%nstd),a(DIMENSION%msh),b(DIMENSION%msh),ain(DIMENSION%msh),ahelp(DIMENSION%msh)
      43         590 :       REAL occ_h(DIMENSION%nstd,2)
      44         590 :       INTEGER kappa(DIMENSION%nstd),nprnc(DIMENSION%nstd)
      45             :       CHARACTER(LEN=20) :: attributes(6)
      46             :       REAL stateEnergies(29)
      47             :       !     ..
      48             :       
      49         295 :       c = c_light(1.0)
      50         295 :       seig = 0.
      51             :       !
      52         295 :       IF (input%frcor) THEN
      53           0 :          DO  n = 1,atoms%ntype
      54           0 :             rnot = atoms%rmsh(1,n) ; dxx = atoms%dx(n)
      55           0 :             ncmsh = NINT( LOG( (atoms%rmt(n)+10.0)/rnot ) / dxx + 1 )
      56           0 :             ncmsh = MIN( ncmsh, DIMENSION%msh )
      57             :             !     --->    update spherical charge density
      58           0 :             DO  i = 1,atoms%jri(n)
      59           0 :                rhoc(i) = rhc(i,n,jspin)
      60           0 :                rho(i,0,n,jspin) = rho(i,0,n,jspin) + rhoc(i)/sfp_const
      61             :             ENDDO
      62             :             !     ---> for total energy calculations, determine the sum of the
      63             :             !     ---> eigenvalues by requiring that the core kinetic energy
      64             :             !     ---> remains constant.
      65           0 :             DO  i = 1,atoms%jri(n)
      66           0 :                rhoc(i) = rhoc(i)*vr(i,n)/atoms%rmsh(i,n)
      67             :             ENDDO
      68           0 :             nm = atoms%jri(n)
      69           0 :             CALL intgr3(rhoc,atoms%rmsh(1,n),atoms%dx(n),nm,rhos)
      70           0 :             sea = tec(n,jspin) + rhos
      71           0 :             WRITE (6,FMT=8030) n,jspin,tec(n,jspin),sea
      72           0 :             seig = seig + atoms%neq(n)*sea
      73             :          ENDDO
      74             :          RETURN
      75             :       END IF
      76             : 
      77             :       !     ---> set up densities
      78         858 :       DO  jatom = 1,atoms%ntype
      79         563 :          sume = 0.
      80         563 :          z = atoms%zatom(jatom)
      81             :          !         rn = rmt(jatom)
      82         563 :          dxx = atoms%dx(jatom)
      83         563 :          bmu = 0.0
      84         563 :          CALL setcor(jatom,input%jspins,atoms,input,bmu,nst,kappa,nprnc,occ_h)
      85         563 :          IF ((bmu > 99.)) THEN
      86           0 :             occ(1:nst) = input%jspins *  occ_h(1:nst,jspin)
      87             :          ELSE
      88         563 :             occ(1:nst) = occ_h(1:nst,1)
      89             :          ENDIF
      90         563 :          rnot = atoms%rmsh(1,jatom)
      91         563 :          d = EXP(atoms%dx(jatom))
      92         563 :          ncmsh = NINT( LOG( (atoms%rmt(jatom)+10.0)/rnot ) / dxx + 1 )
      93         563 :          ncmsh = MIN( ncmsh, DIMENSION%msh )
      94         563 :          rn = rnot* (d** (ncmsh-1))
      95         563 :          WRITE (6,FMT=8000) z,rnot,dxx,atoms%jri(jatom)
      96      333050 :          DO  j = 1,atoms%jri(jatom)
      97      332487 :             rhoss(j)     = 0.0
      98      332487 :             if(present(EnergyDen)) rhoss_aux(j) = 0.0
      99      333050 :             vrd(j) = vr(j,jatom)
     100             :          ENDDO
     101             :          !
     102         563 :          IF (input%l_core_confpot) THEN
     103             :             !--->    linear extension of the potential with slope t1 / a.u.
     104         563 :             t1=0.125
     105             :             t1 = MAX( (vrd(atoms%jri(jatom)) - vrd(atoms%jri(jatom)-1)*d)*&
     106         563 :                      d / (atoms%rmt(jatom)**2 * (d-1) ) , t1)
     107         563 :             t2=vrd(atoms%jri(jatom))/atoms%rmt(jatom)-atoms%rmt(jatom)*t1
     108         563 :             rr = atoms%rmt(jatom)
     109             :          ELSE
     110           0 :             t2 = vrd(atoms%jri(jatom)) / ( atoms%jri(jatom) - ncmsh )
     111             :          ENDIF
     112         563 :          IF ( atoms%jri(jatom) < ncmsh) THEN
     113       47180 :             DO  i = atoms%jri(jatom) + 1,ncmsh
     114       46617 :                rhoss(i) = 0.
     115       46617 :                if(present(EnergyDen)) rhoss_aux(i) = 0.0
     116       47180 :                IF (input%l_core_confpot) THEN
     117       46617 :                   rr = d*rr
     118       46617 :                   vrd(i) = rr*( t2 + rr*t1 )
     119             :                   !               vrd(i) = 2*vrd(jri(jatom)) - rr*( t2 + rr*t1 )
     120             :                ELSE
     121           0 :                   vrd(i) = vrd(atoms%jri(jatom)) + t2* (i-atoms%jri(jatom))
     122             :                ENDIF
     123             :                !
     124             :             ENDDO
     125             :          END IF
     126             : 
     127         563 :          nst = atoms%ncst(jatom)        ! for lda+U
     128             : 
     129         563 :          IF (input%gw==1 .OR. input%gw==3)&
     130           0 :               &                      WRITE(15) nst,atoms%rmsh(1:atoms%jri(jatom),jatom)
     131             : 
     132         563 :          stateEnergies = 0.0
     133        4325 :          DO  korb = 1,nst
     134        4325 :             IF (occ(korb) /= 0.0) THEN
     135        3762 :                fn = nprnc(korb)
     136        3762 :                fj = iabs(kappa(korb)) - .5e0
     137             : 
     138        3762 :                weight = 2*fj + 1.e0
     139        3762 :                IF (bmu > 99.) weight = occ(korb)
     140             : 
     141        3762 :                fl = fj + (.5e0)*isign(1,kappa(korb))
     142             : 
     143        3762 :                eig        = -2* (z/ (fn+fl))**2
     144             : 
     145        3762 :                CALL differ(fn,fl,fj,c,z,dxx,rnot,rn,d,ncmsh,vrd, eig, a,b,ierr)
     146        3762 :                stateEnergies(korb) = eig
     147        3762 :                WRITE (6,FMT=8010) fn,fl,fj,eig,weight
     148             : 
     149        3762 :                IF (ierr/=0)  CALL juDFT_error("error in core-level routine" ,calledby ="cored")
     150        3762 :                IF (input%gw==1 .OR. input%gw==3) WRITE (15) NINT(fl),weight,eig,&
     151           0 :                   a(1:atoms%jri(jatom)),b(1:atoms%jri(jatom))
     152             : 
     153        3762 :                sume = sume + weight*eig/input%jspins
     154     2546042 :                DO j = 1,ncmsh
     155     2542280 :                   rhcs(j)  = weight* (a(j)**2+b(j)**2)
     156     2546042 :                   rhoss(j) = rhoss(j) + rhcs(j)
     157             :                ENDDO
     158             : 
     159        3762 :                IF(present(EnergyDen)) THEN
     160             :                   !rhoss_aux = rhoss
     161           0 :                   DO j = 1,ncmsh
     162             :                      ! for energy density we want to multiply the weights
     163             :                      ! with the eigenenergies
     164           0 :                      rhoss_aux(j) = rhoss_aux(j) + (rhcs(j) * eig)
     165             :                   ENDDO
     166             :                ENDIF
     167             :             ENDIF
     168             :          ENDDO
     169             : 
     170             :          !     ---->update spherical charge density rho with the core density.
     171             :          !     ---->for spin-polarized (jspins=2), take only half the density
     172         563 :          nm = atoms%jri(jatom)
     173      333050 :          DO  j = 1,nm
     174      332487 :             rhoc(j) = rhoss(j)/input%jspins
     175      333050 :             rho(j,0,jatom,jspin) = rho(j,0,jatom,jspin) + rhoc(j)/sfp_const
     176             :          ENDDO
     177             : 
     178         563 :          IF(present(EnergyDen)) then
     179           0 :             DO  j = 1,nm
     180             :                EnergyDen(j,0,jatom,jspin) = EnergyDen(j,0,jatom,jspin) &
     181           0 :                                             + rhoss_aux(j) /(input%jspins * sfp_const)
     182             :             ENDDO
     183             :          ENDIF
     184             : 
     185         563 :          rhc(1:ncmsh,jatom,jspin)   = rhoss(1:ncmsh) / input%jspins
     186         563 :          rhc(ncmsh+1:DIMENSION%msh,jatom,jspin) = 0.0
     187             : 
     188         563 :          seig = seig + atoms%neq(jatom)*sume
     189      333050 :          DO  i = 1,nm
     190      333050 :             rhoc(i) = rhoc(i)*vr(i,jatom)/atoms%rmsh(i,jatom)
     191             :          ENDDO
     192         563 :          CALL intgr3(rhoc,atoms%rmsh(1,jatom),atoms%dx(jatom),nm,rhs)
     193         563 :          tec(jatom,jspin) = sume - rhs
     194         563 :          WRITE (6,FMT=8030) jatom,jspin,tec(jatom,jspin),sume
     195             :          
     196             :          !     ---> simpson integration
     197         563 :          rad = atoms%rmt(jatom)
     198         563 :          q = rad*rhoss(nm)/2.
     199       23596 :          DO  nm1 = nm + 1,ncmsh - 1,2
     200       23033 :             rad = d*rad
     201       23033 :             q = q + 2*rad*rhoss(nm1)
     202       23033 :             rad = d*rad
     203       23596 :             q = q + rad*rhoss(nm1+1)
     204             :          ENDDO
     205         563 :          q = 2*q*dxx/3
     206             :          !+sb
     207         563 :          WRITE (6,FMT=8020) q/input%jspins
     208             :          !-sb
     209         563 :          qint(jatom,jspin) = q*atoms%neq(jatom)
     210        3941 :          attributes = ''
     211         563 :          WRITE(attributes(1),'(i0)') jatom
     212         563 :          WRITE(attributes(2),'(i0)') NINT(z)
     213         563 :          WRITE(attributes(3),'(i0)') jspin
     214         563 :          WRITE(attributes(4),'(f18.10)') tec(jatom,jspin)
     215         563 :          WRITE(attributes(5),'(f18.10)') sume
     216         563 :          WRITE(attributes(6),'(f9.6)') q/input%jspins
     217             :          CALL openXMLElementForm('coreStates',(/'atomType     ','atomicNumber ','spin         ','kinEnergy    ',&
     218             :                                                 'eigValSum    ','lostElectrons'/),&
     219         563 :                                  attributes,RESHAPE((/8,12,4,9,9,13,6,3,1,18,18,9/),(/6,2/)))
     220        4325 :          DO korb = 1, atoms%ncst(jatom)
     221        3762 :             fj = iabs(kappa(korb)) - .5e0
     222        3762 :             weight = 2*fj + 1.e0
     223        3762 :             IF (bmu > 99.) weight = occ(korb)
     224        3762 :             fl = fj + (.5e0)*isign(1,kappa(korb))
     225       26334 :             attributes = ''
     226        3762 :             WRITE(attributes(1),'(i0)') nprnc(korb)
     227        3762 :             WRITE(attributes(2),'(i0)') NINT(fl)
     228        3762 :             WRITE(attributes(3),'(f4.1)') fj
     229        3762 :             WRITE(attributes(4),'(f20.10)') stateEnergies(korb)
     230        3762 :             WRITE(attributes(5),'(f15.10)') weight
     231             :             CALL writeXMLElementForm('state',(/'n     ','l     ','j     ','energy','weight'/),&
     232        4325 :                                      attributes(1:5),RESHAPE((/1,1,1,6,6,2,2,4,20,15/),(/5,2/)))
     233             :          END DO
     234         858 :          CALL closeXMLElement('coreStates')
     235             :       ENDDO
     236             : 
     237             :       RETURN
     238             : 
     239             : 8000  FORMAT (/,/,10x,'z=',f4.0,5x,'r(1)=',e14.6,5x,'dx=',f9.6,5x,&
     240             :           &       'm.t.index=',i4,/,15x,'n',4x,'l',5x,'j',4x,'energy',7x,&
     241             :           &       'weight')
     242             : 8010  FORMAT (12x,2f5.0,f6.1,f10.4,f12.4)
     243             : 8020  FORMAT (f20.8,'  electrons lost from core.')
     244             : 8030  FORMAT (10x,'atom type',i3,'  (spin',i2,') ',/,10x,&
     245             :           &       'kinetic energy=',e20.12,5x,'sum of the eigenvalues=',&
     246             :           &       e20.12)
     247         295 :    END SUBROUTINE cored
     248             : END MODULE m_cored
     249             : 

Generated by: LCOV version 1.13