LCOV - code coverage report
Current view: top level - optional - stden.f90 (source / functions) Hit Total Coverage
Test: combined.info Lines: 96 170 56.5 %
Date: 2019-09-08 04:53:50 Functions: 1 1 100.0 %

          Line data    Source code
       1             : MODULE m_stden
       2             : 
       3             : USE m_juDFT
       4             : !     ************************************************************
       5             : !     generate flapw starting density by superposition of
       6             : !     atomic densities. the non-spherical terms inside
       7             : !     the spheres are obtained by a least squares fit
       8             : !     and the interstitial and vacuum warping terms are calculated
       9             : !     by a fast fourier transform.
      10             : !     e. wimmer   nov. 1984       c.l.fu diagonized 1987
      11             : !     *************************************************************
      12             : 
      13             : CONTAINS
      14             : 
      15          56 : SUBROUTINE stden(mpi,sphhar,stars,atoms,sym,DIMENSION,vacuum,&
      16             :                  input,cell,xcpot,obsolete,noco,oneD)
      17             : 
      18             :    USE m_constants
      19             :    USE m_qsf
      20             :    USE m_checkdopall
      21             :    USE m_cdnovlp
      22             :    USE m_cdn_io
      23             :    USE m_qfix
      24             :    USE m_atom2
      25             :    USE m_types
      26             :    USE m_types_xcpot_inbuild
      27             :    USE m_juDFT_init
      28             : 
      29             :    IMPLICIT NONE
      30             : 
      31             :    TYPE(t_mpi),INTENT(IN)      :: mpi
      32             :    TYPE(t_atoms),INTENT(IN)    :: atoms
      33             :    TYPE(t_dimension),INTENT(IN):: DIMENSION
      34             :    TYPE(t_sphhar),INTENT(IN)   :: sphhar
      35             :    TYPE(t_obsolete),INTENT(IN) :: obsolete
      36             :    TYPE(t_sym),INTENT(IN)      :: sym
      37             :    TYPE(t_stars),INTENT(IN)    :: stars
      38             :    TYPE(t_noco),INTENT(IN)     :: noco
      39             :    TYPE(t_oneD),INTENT(IN)     :: oneD
      40             :    TYPE(t_input),INTENT(IN)    :: input
      41             :    TYPE(t_vacuum),INTENT(IN)   :: vacuum
      42             :    TYPE(t_cell),INTENT(IN)     :: cell
      43             :    CLASS(t_xcpot),INTENT(IN)   :: xcpot
      44             : 
      45             :    ! Local type instances
      46          56 :    TYPE(t_potden)   :: den
      47          56 :    TYPE(t_enpara)   :: enpara
      48          56 :    TYPE(t_xcpot_inbuild)    :: xcpot_dummy
      49             : 
      50             :    ! Local Scalars
      51             :    REAL d,del,fix,h,r,rnot,z,bm,qdel,va
      52             :    REAL denz1(1,1),vacxpot(1,1),vacpot(1,1) 
      53             :    INTEGER i,ivac,iza,j,jr,k,n,n1,ispin 
      54             :    INTEGER nw,ilo,natot,nat 
      55             : 
      56             :    ! Local Arrays
      57          56 :    REAL,    ALLOCATABLE :: vbar(:,:)
      58         168 :    REAL,    ALLOCATABLE :: rat(:,:),eig(:,:,:),sigm(:)
      59          56 :    REAL,    ALLOCATABLE :: rh(:,:,:),rh1(:,:,:),rhoss(:,:)
      60          56 :    REAL,    ALLOCATABLE :: vacpar(:)
      61         112 :    INTEGER lnum(DIMENSION%nstd,atoms%ntype),nst(atoms%ntype) 
      62         112 :    INTEGER jrc(atoms%ntype)
      63         112 :    LOGICAL l_found(0:3),llo_found(atoms%nlod),l_enpara,l_st
      64         112 :    REAL                 :: occ(SIZE(atoms%corestateoccs,1),2)
      65             :    ! Data statements
      66             :    DATA del/1.e-6/
      67             :    PARAMETER (l_st=.true.)
      68             : 
      69             :    IF (input%jspins > input%jspins) CALL juDFT_error("input%jspins > input%jspins", calledby = "stden")
      70             : 
      71          56 :    CALL den%init(stars,atoms,sphhar,vacuum,noco,input%jspins,POTDEN_TYPE_DEN)
      72             : 
      73          56 :    ALLOCATE ( rat(DIMENSION%msh,atoms%ntype),eig(DIMENSION%nstd,input%jspins,atoms%ntype) )
      74          56 :    ALLOCATE ( rh(DIMENSION%msh,atoms%ntype,input%jspins),rh1(DIMENSION%msh,atoms%ntype,input%jspins) )
      75          56 :    ALLOCATE ( vbar(2,atoms%ntype),sigm(vacuum%nmzd) )
      76          56 :    ALLOCATE ( rhoss(DIMENSION%msh,input%jspins) )
      77             : 
      78         138 :    rh = 0.0
      79         138 :    rhoss = 0.0
      80             : 
      81          56 :    IF (mpi%irank == 0) THEN
      82             :       ! if sigma is not 0.0, then divide this charge among all atoms
      83          28 :       IF ( ABS(input%sigma).LT. 1.e-6) THEN
      84          28 :          qdel = 0.0
      85             :       ELSE
      86           0 :          natot = 0
      87           0 :          DO n = 1, atoms%ntype
      88           0 :             IF (atoms%zatom(n).GE.1.0) natot = natot + atoms%neq(n)
      89             :          END DO
      90           0 :          qdel = 2.*input%sigma/natot
      91             :       END IF
      92             : 
      93          28 :       WRITE (6,FMT=8000)
      94             :       8000 FORMAT (/,/,/,' superposition of atomic densities',/,/,' original atomic densities:',/)
      95          91 :       DO n = 1,atoms%ntype
      96          63 :          r = atoms%rmsh(1,n)
      97          63 :          d = EXP(atoms%dx(n))
      98          63 :          jrc(n) = 0
      99       98649 :          DO WHILE (r < atoms%rmt(n) + 20.0) 
     100       49279 :             IF (jrc(n) > DIMENSION%msh) CALL juDFT_error("increase msh in fl7para!",calledby ="stden")
     101       49279 :             jrc(n) = jrc(n) + 1
     102       49279 :             rat(jrc(n),n) = r
     103       49342 :             r = r*d
     104             :          END DO
     105             :       END DO
     106             : 
     107             :       ! Generate the atomic charge densities
     108          91 :       DO n = 1,atoms%ntype
     109          63 :          z = atoms%zatom(n)
     110          63 :          r = atoms%rmt(n)
     111          63 :          h = atoms%dx(n)
     112          63 :          jr = atoms%jri(n)
     113          63 :          IF (input%jspins.EQ.2) THEN
     114          24 :             bm = atoms%bmu(n)
     115             :          ELSE
     116             :             bm = 0.
     117             :          END IF
     118          63 :          occ=atoms%coreStateOccs(:,:,n)
     119             :          ! check whether this atom has been done already
     120         111 :          DO n1 = 1, n - 1
     121          47 :             IF (ABS(z-atoms%zatom(n1)).GT.del) CYCLE
     122          26 :             IF (ABS(r-atoms%rmt(n1)).GT.del) CYCLE
     123          26 :             IF (ABS(h-atoms%dx(n1)).GT.del) CYCLE
     124          26 :             IF (ABS(bm-atoms%bmu(n1)).GT.del) CYCLE
     125          23 :             IF (ANY(ABS(occ(:,:)-atoms%coreStateOccs(:,:,n1))>del)) CYCLE
     126          23 :             IF (jr.NE.atoms%jri(n1)) CYCLE
     127          57 :             DO ispin = 1, input%jspins
     128       25543 :                DO i = 1,jrc(n) ! dimension%msh
     129       25520 :                   rh(i,n,ispin) = rh(i,n1,ispin)
     130             :                END DO
     131             :             END DO
     132          23 :             nst(n) = nst(n1)
     133          23 :             vbar(1,n) = vbar(1,n1)
     134          23 :             vbar(input%jspins,n) = vbar(input%jspins,n1)
     135         213 :             DO i = 1, nst(n1)
     136         190 :                lnum(i,n)  = lnum(i,n1)
     137         190 :                eig(i,1,n) = eig(i,1,n1)
     138         213 :                eig(i,input%jspins,n) = eig(i,input%jspins,n1)
     139             :             END DO
     140          64 :             GO TO 70
     141             :          END DO
     142             :          !--->    new atom
     143          40 :          rnot = atoms%rmsh(1,n)
     144          40 :          IF (z.LT.1.0) THEN
     145           0 :             va = max(z,1.e-8)/(input%jspins*sfp_const*atoms%volmts(n))
     146           0 :             DO ispin = 1, input%jspins
     147           0 :                DO i = 1,jrc(n) ! dimension%msh
     148           0 :                   rh(i,n,ispin) = va/rat(i,n)**2
     149             :                END DO
     150             :             END DO
     151             :          ELSE
     152             :             CALL atom2(DIMENSION,atoms,xcpot,input,n,jrc(n),rnot,qdel,&
     153          40 :                        rhoss,nst(n),lnum(1,n),eig(1,1,n),vbar(1,n))
     154          93 :             DO ispin = 1, input%jspins
     155       41738 :                DO i = 1, jrc(n) ! dimension%msh
     156       41751 :                   rh(i,n,ispin) = rhoss(i,ispin)
     157             :                END DO
     158             :             END DO
     159             :          END IF
     160             :          ! list atomic density
     161          40 :          iza = atoms%zatom(n) + 0.0001
     162          40 :          WRITE (6,FMT=8030) namat_const(iza)
     163             :          8030 FORMAT (/,/,' atom: ',a2,/)
     164             :          8040 FORMAT (4 (3x,i5,f8.5,f12.6))
     165          28 :          70 CONTINUE
     166             :       END DO
     167             : 
     168             :       !roa+
     169             :       ! use cdnovlp to generate total density out of atom densities
     170          69 :       DO ispin = 1, input%jspins
     171          41 :          nat = 1
     172         156 :          DO n = 1,atoms%ntype
     173       67271 :             DO i = 1, jrc(n)
     174       67271 :                rh1(i,n,ispin) = rh(i,n,ispin)*fpi_const*rat(i,n)**2
     175             :             END DO
     176          87 :             rh1(jrc(n):DIMENSION%msh,n,ispin) = 0.0
     177             :             ! prepare spherical mt charge
     178       56454 :             DO i = 1,atoms%jri(n)
     179       56454 :                den%mt(i,0,n,ispin) = rh(i,n,ispin)*sfp_const*atoms%rmsh(i,n)**2
     180             :             END DO
     181             :             ! reset nonspherical mt charge
     182        1784 :             DO k = 1,sphhar%nlh(atoms%ntypsy(nat))
     183     2198906 :                DO j = 1,atoms%jri(n)
     184     1100258 :                   den%mt(j,k,n,ispin) = 0.e0
     185             :                END DO
     186             :             END DO
     187         128 :             nat = nat + atoms%neq(n)
     188             :          END DO
     189             :       END DO ! ispin
     190             :    END IF ! mpi%irank == 0
     191             : 
     192         138 :    DO ispin = 1, input%jspins
     193             :       CALL cdnovlp(mpi,sphhar,stars,atoms,sym,DIMENSION,vacuum,&
     194             :                    cell,input,oneD,l_st,ispin,rh1(:,:,ispin),&
     195         138 :                    den%pw,den%vacxy,den%mt,den%vacz)
     196             :       !roa-
     197             :    END DO
     198             : 
     199          56 :    IF (mpi%irank == 0) THEN
     200             : 
     201             :       ! Check the normalization of total density
     202          28 :       CALL qfix(mpi,stars,atoms,sym,vacuum,sphhar,input,cell,oneD,den,.FALSE.,.FALSE.,.true.,fix)
     203          28 :       z=SUM(atoms%neq(:)*atoms%zatom(:))
     204          28 :       IF (ABS(fix*z-z)>0.5) THEN
     205             :          CALL judft_warn("Starting density not charge neutral",hint= &
     206           0 :                          "Your electronic configuration might be broken",calledby="stden.f90")
     207             :       END IF
     208             : 
     209             :       ! Write superposed density onto density file
     210          28 :       den%iter = 0
     211             :       CALL writeDensity(stars,vacuum,atoms,cell,sphhar,input,sym,oneD,CDN_ARCHIVE_TYPE_CDN1_const,CDN_INPUT_DEN_const,&
     212          28 :                         1,-1.0,0.0,.TRUE.,den)
     213             : 
     214             :       ! Check continuity
     215          28 :       IF (input%vchk) THEN
     216           6 :          DO ispin = 1, input%jspins
     217           3 :             WRITE (6,'(a8,i2)') 'spin No.',ispin
     218             :             CALL checkDOPAll(input,dimension,sphhar,stars,atoms,sym,vacuum,oneD,&
     219           6 :                            cell,den,ispin)
     220             :          END DO ! ispin = 1, input%jspins
     221             :       END IF ! input%vchk
     222             : 
     223          28 :       l_enpara = .FALSE.
     224          28 :       INQUIRE (file='enpara',exist=l_enpara)
     225          28 :       l_enpara = l_enpara.OR.input%l_inpXML
     226             : 
     227             :       ! set up parameters for enpara-file
     228          28 :       IF ((juDFT_was_argument("-genEnpara")).AND..NOT.l_enpara) THEN
     229           0 :          CALL enpara%init(atoms,input%jspins)
     230             : 
     231           0 :          enpara%lchange = .TRUE.
     232           0 :          enpara%llochg = .TRUE.
     233             :                 
     234           0 :          DO ispin = 1, input%jspins
     235             :             ! vacpar is taken as highest occupied level from atomic eigenvalues
     236             :             ! vacpar (+0.3)  serves as fermi level for film or bulk
     237           0 :             vacpar(1) = -999.9
     238           0 :             DO n = 1,atoms%ntype
     239           0 :                vacpar(1) = MAX( vacpar(1),eig(nst(n),ispin,n) )
     240             :             END DO
     241           0 :             IF (.NOT.input%film) vacpar(1) = vacpar(1) + 0.4
     242           0 :             vacpar(2) = vacpar(1)
     243           0 :             DO n = 1,atoms%ntype
     244           0 :                enpara%skiplo(n,ispin) = 0
     245           0 :                DO i = 0, 3
     246           0 :                   l_found(i) = .FALSE.
     247           0 :                   enpara%el0(i,n,ispin) = vacpar(1)
     248             :                END DO
     249           0 :                DO i = 1, atoms%nlo(n)
     250           0 :                   llo_found(i) = .FALSE.
     251           0 :                   enpara%ello0(i,n,ispin) = vacpar(1) - 1.5
     252             :                END DO
     253             : 
     254             :                ! take energy-parameters from atomic calculation
     255           0 :                DO i = nst(n), 1, -1 
     256           0 :                   IF (.NOT.input%film) eig(i,ispin,n) = eig(i,ispin,n) + 0.4
     257           0 :                   IF (.NOT.l_found(lnum(i,n)).AND.(lnum(i,n).LE.3)) THEN
     258           0 :                      enpara%el0(lnum(i,n),n,ispin) = eig(i,ispin,n)
     259           0 :                      IF (enpara%el0(lnum(i,n),n,ispin).LT.input%ellow) THEN
     260           0 :                         enpara%el0(lnum(i,n),n,ispin) = vacpar(1)
     261           0 :                         l_found(lnum(i,n))  = .TRUE.
     262             :                      END IF
     263           0 :                      IF (enpara%el0(lnum(i,n),n,ispin).LT.input%elup) THEN
     264           0 :                         l_found(lnum(i,n))  = .TRUE.
     265             :                      END IF
     266             :                   ELSE
     267           0 :                      IF (l_found(lnum(i,n)).AND.(atoms%nlo(n).GT.0)) THEN
     268           0 :                         DO ilo = 1, atoms%nlo(n)
     269           0 :                            IF (atoms%llo(ilo,n).EQ.lnum(i,n)) THEN
     270           0 :                               IF (.NOT.llo_found(ilo)) THEN
     271           0 :                                  enpara%ello0(ilo,n,ispin) = eig(i,ispin,n)
     272           0 :                                  IF ((enpara%ello0(ilo,n,ispin).GT.input%elup).OR.&
     273             :                                      (enpara%ello0(ilo,n,ispin).LT.input%ellow)) THEN
     274           0 :                                     enpara%ello0(ilo,n,ispin)= vacpar(1)
     275           0 :                                  ELSE IF (atoms%l_dulo(ilo,n)) THEN
     276           0 :                                     enpara%ello0(ilo,n,ispin)= enpara%el0(atoms%llo(ilo,n),n,ispin)
     277             :                                  ELSE
     278           0 :                                     enpara%skiplo(n,ispin) = enpara%skiplo(n,ispin) + 2*atoms%llo(ilo,n)+1
     279             :                                  END IF
     280           0 :                                  llo_found(ilo) = .TRUE.
     281             :                                  IF ((enpara%el0(atoms%llo(ilo,n),n,ispin)-enpara%ello0(ilo,n,ispin).LT.0.5)&
     282           0 :                                      .AND.(atoms%llo(ilo,n).GE.0)) THEN
     283           0 :                                     enpara%el0(atoms%llo(ilo,n),n,ispin) = vacpar(1)
     284           0 :                                     IF (atoms%l_dulo(ilo,n)) enpara%ello0(ilo,n,ispin)= enpara%el0(atoms%llo(ilo,n),n,ispin)
     285             :                                  END IF
     286             :                               END IF
     287             :                            END IF
     288             :                         END DO ! ilo = 1, atoms%nlo(n)
     289             :                      END IF
     290             :                   END IF ! .NOT.l_found(lnum(i,n)).AND.(lnum(i,n).LE.3)
     291             :                END DO ! i = nst(n), 1, -1 
     292           0 :                IF (obsolete%lepr.EQ.1) THEN
     293           0 :                   DO i = 0, 3
     294           0 :                      enpara%el0(i,n,ispin) = enpara%el0(i,n,ispin) - vbar(ispin,n)
     295             :                   END DO
     296           0 :                   DO ilo = 1,atoms%nlo(n)
     297           0 :                      enpara%ello0(ilo,n,ispin) = enpara%ello0(ilo,n,ispin) - vbar(ispin,n)
     298             :                   END DO
     299             :                END IF
     300             :             END DO ! atom types
     301             : 
     302           0 :             IF (input%film) THEN
     303             :                ! get guess for vacuum parameters
     304             :                ! YM : in 1D case should be modified, but also not that bad like this
     305             :                !
     306             :                ! generate coulomb potential by integrating inward to z1
     307             : 
     308           0 :                DO ivac = 1, vacuum%nvac
     309           0 :                   CALL xcpot_dummy%init("vwn",.FALSE.,atoms%ntype)
     310           0 :                   DO i=1,vacuum%nmz
     311           0 :                      sigm(i) = (i-1)*vacuum%delz*den%vacz(i,ivac,ispin)
     312             :                   END DO
     313           0 :                   CALL qsf(vacuum%delz,sigm,vacpar(ivac),vacuum%nmz,0)
     314           0 :                   denz1 = den%vacz(1,ivac,ispin)          ! get estimate for potential at vacuum boundary
     315           0 :                   CALL xcpot%get_vxc(1,denz1,vacpot,vacxpot)
     316             :                   ! seems to be the best choice for 1D not to substract vacpar
     317           0 :                   IF (.NOT.oneD%odi%d1) THEN
     318           0 :                      vacpot = vacpot - fpi_const*vacpar(ivac)
     319             :                   END IF
     320           0 :                   IF (obsolete%lepr.EQ.1) THEN
     321           0 :                      vacpar(ivac) = -0.2 - vacpot(1,1)
     322           0 :                      WRITE (6,'(" vacuum",i2," reference energy =",f12.6)') ivac,vacpot
     323             :                   ELSE
     324           0 :                      vacpar(ivac) = vacpot(1,1)
     325             :                   END IF
     326             :                END DO
     327           0 :                IF (vacuum%nvac.EQ.1) vacpar(2) = vacpar(1)
     328             :             END IF
     329             : 
     330           0 :             IF (obsolete%lepr.EQ.1) THEN
     331           0 :                enpara%enmix = 0.3
     332             :             ELSE
     333           0 :                enpara%enmix = 1.0
     334             :             END IF
     335             : 
     336             :             
     337           0 :             enpara%evac0(:,ispin)=vacpar(:SIZE(enpara%evac0,1))
     338             :            
     339             :          END DO ! ispin
     340           0 :          CALL enpara%WRITE(atoms,input%jspins,input%film)
     341             :       END IF
     342             :    END IF ! mpi%irank == 0
     343             : 
     344         112 : END SUBROUTINE stden
     345             : 
     346             : END MODULE m_stden

Generated by: LCOV version 1.13