atom2.f90 9.5 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75
      MODULE m_atom2
      use m_juDFT
!     *************************************************************
!     fully relativistic atomic program based on the subroutines
!     differ, outint and inwint by d.d.koelling
!     erich wimmer     august 1981
!     modified for use in start-up generator.  m.w. april 1982
!     modified by adding stabilizing well. m. weinert may 1990
!     *************************************************************
      CONTAINS
      SUBROUTINE atom2(&
     &                 dimension,atoms,xcpot,input,ntyp,jrc,rnot1,&
     &                 qdel,&
     &                 rhoss,nst,lnum,eig,vbar)

      USE m_intgr, ONLY : intgr1,intgr0
      USE m_constants
      USE m_xcall, ONLY : vxcall
      USE m_potl0
      USE m_stpot1
      USE m_setcor
      USE m_differ
      USE m_types
      IMPLICIT NONE
!     ..
!     .. Scalar Arguments ..
      TYPE(t_dimension),INTENT(IN)  :: dimension
      TYPE(t_atoms),INTENT(IN)      :: atoms
      TYPE(t_xcpot),INTENT(IN)      :: xcpot
      TYPE(t_input),INTENT(IN)      :: input
      INTEGER,INTENT (IN)  :: jrc ,ntyp    
      REAL,   INTENT (IN)  :: rnot1  ,qdel
      REAL,   INTENT (OUT) :: rhoss(:,:) !(mshd,dimension%jspd),
      REAL,   INTENT (OUT) :: eig(dimension%nstd,dimension%jspd),vbar(dimension%jspd)
      INTEGER,INTENT (OUT) :: nst,lnum(dimension%nstd)
!     ..
!     .. Local Scalars ..
      REAL c,d,delrv,dist,distol,e,fisr,fj,fl,fn,h,&
     &     p,p1,pmax,pmin,r,r3,rn,rnot,z,zero,bmu_l,rho
      INTEGER i,inr0,it,itmax,k,l,n,ispin,kk,ierr,msh_l
      LOGICAL conv,lastit,l_start
!     ..
!     .. Local Arrays ..
      REAL a(jrc),b(jrc),dens(jrc),occ(dimension%nstd,dimension%jspd)
      REAL rad(jrc),rev(dimension%nstd,dimension%jspd),ahelp(jrc),ain(jrc),&
     &     rh(jrc),vr(jrc),f(0:3),&
     &     vr1(jrc,dimension%jspd),vr2(jrc,dimension%jspd),vx(dimension%msh,dimension%jspd),vxc(dimension%msh,dimension%jspd)
      INTEGER kappa(dimension%nstd),nprnc(dimension%nstd)
!     ..
!     ..
!     .. Data statements ..
!---->     distol set from 1.0e-6 to 1.0e-3
      DATA zero,distol/0.0e0,1.0e-3/
!     ..
      c = c_light(1.0)
      vxc(:,:) = 0.0
      vx (:,:) = 0.0
!
      WRITE (6,FMT=8000)
 8000 FORMAT (' subroutine atom2 entered')
      z = atoms%zatom(ntyp)
      n = jrc
      rnot = rnot1
      itmax = 100
      pmin = 0.01e0
      pmax = 0.2e0
      h = atoms%dx(ntyp)
      d = exp(h)
      r = rnot
      DO  i = 1,n
         rad(i) = r
         r = r*d
      enddo
      rn = rad(n)
      bmu_l = atoms%bmu(ntyp)
76 77 78
      IF (bmu_l>0.001.AND.atoms%numStatesProvided(ntyp).NE.0) CALL &
           judft_warn("You specified both: inital moment and occupation numbers.", &
           hint="The inital moment will be ignored, set magMom=0.0",calledby="atom2.f90")
79
      CALL setcor(ntyp,input%jspins,atoms,input,bmu_l, nst,kappa,nprnc,occ)
80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180

!
!--->   for electric field case (sigma.ne.0), add the extra charge
!--->   to the uppermost level; ignore the possible problem that
!--->   the occupations may not be between 0 and 2
      IF (input%jspins.EQ.1) THEN
        occ(nst,1) = occ(nst,1) + qdel
      ELSE
        occ(nst,1) = occ(nst,1) + qdel/2.
        occ(nst,input%jspins) = occ(nst,input%jspins) + qdel/2.
      ENDIF
!
      CALL stpot1(&
     &           jrc,n,z,rad,&
     &           vr1)
      DO i = 1,n
         vr1(i,input%jspins) = vr1(i,1)
      ENDDO
!
!     start iterating
!
      lastit = .false.
      conv = .true.
      delrv = 0.100
      inr0 = log(5.0/rnot)/h + 1.5
      DO 180 it = 1,itmax
         DO ispin = 1,input%jspins
!
!---->     load potential
           DO i = 1,n
              vr(i) = vr1(i,ispin)
           ENDDO
!----> adding stabilizing well: m. weinert
           DO i = inr0,n
              vr(i) = vr(i) + rad(i)*delrv* (rad(i)-rad(inr0))
           ENDDO
!---->     note that vr contains r*v(r) in hartree units
           DO i = 1,n
             rhoss(i,ispin) = zero
           ENDDO
           if (lastit) THEN
                inquire(file="startcharges",exist=l_start)
                if (l_start) then
                    OPEN (61,file="startcharges")
                    DO WHILE (.true.)
                        read(61,*,end=888,err=888) i,rho
                        if (i==z) then
                            occ(nst,1)=occ(nst,1)+rho
                            goto 888
                        endif
                    enddo
  888               continue
                    close(61)
                 endif
              endif
           DO 90 k = 1,nst
              fn = nprnc(k)
              fj = iabs(kappa(k)) - 0.5e0
              fl = fj + 0.5e0*isign(1,kappa(k))
              e = -2* (z/ (fn+fl))**2
              ierr = -1 
              msh_l = jrc 
              DO WHILE (ierr.NE.0) 
                CALL differ(&
     &                      fn,fl,fj,c,z,h,rnot,rn,d,msh_l,vr,&
     &                      e,&
     &                      a,b,ierr)!keep
                msh_l = msh_l - 1
                IF (jrc-msh_l > 100)  CALL juDFT_error(&
     &               "atom2",calledby ="atom2")
              ENDDO
              DO i = msh_l+1, jrc
                a(i) = a(msh_l) 
                b(i) = b(msh_l) 
              ENDDO


              DO i = 1,n
                 rh(i) = occ(k,ispin)* (a(i)**2+b(i)**2)
              ENDDO
!+ldau
              IF (lastit) THEN                         ! calculate slater interals
                 l = int(fl)
!                 write(*,*) nprnc(k),l 
                 DO kk = 0, 2*l, 2                      ! F0 for s, F0 + F2 for p etc.
                   r = rnot
                   DO i = 1, n
                     ain(i) = a(i)**2 * r**(-kk-1)      ! prepare inner integrand
                     r = r * d
                   ENDDO
                   CALL intgr1(ain,rnot,h,n, &          ! integrate&
                              &ahelp)
                   r = rnot
                   DO i = 1, n-1
                     ain(i) =  a(i)**2 * r**kk * (ahelp(n) - ahelp(i))
                     r = r * d
                   ENDDO
                   CALL intgr0(ain,rnot,h,n-1,      &     ! integrate 2nd r&
                              &f(kk/2))

                 ENDDO
181
!                 write(*,*) (hartree_to_ev_const*2*f(kk),kk=0,l)
182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221
              ENDIF
!-ldau
              eig(k,ispin) = e
!---->       calculate <r>
              DO i = 1,n
                 a(i) = (a(i)**2+b(i)**2)*rad(i)
              ENDDO
              CALL intgr1(a,rnot,h,n,b)
              rev(k,ispin) = b(n)
              DO i = 1,n
                 rhoss(i,ispin) = rhoss(i,ispin) + rh(i)
              ENDDO
   90      ENDDO
         ENDDO
!
!     solve poisson's equation
!
         DO i = 1,n
            dens(i) = rhoss(i,1)
         ENDDO
         IF (input%jspins.EQ.2) THEN
           DO i = 1,n
              dens(i) = dens(i) + rhoss(i,input%jspins)
           ENDDO
         ENDIF
         CALL intgr1(dens,rnot,h,n,a)
         DO 110 i = 1,n
            rh(i) = dens(i)/rad(i)
  110    CONTINUE
         CALL intgr1(rh,rnot,h,n,b)
         fisr = b(n)
         DO 120 i = 1,n
            vr(i) = (a(i)+rad(i)* (fisr-b(i))-z)
  120    CONTINUE
!+ta
         DO ispin = 1, input%jspins
           DO i = 1,n
             rhoss(i,ispin) = rhoss(i,ispin) / (fpi_const*rad(i)**2)
           ENDDO
         ENDDO
222 223 224 225 226 227 228 229 230
         IF (xcpot%is_gga()) THEN
            CALL potl0(&
                 xcpot,DIMENSION%msh,DIMENSION%jspd,input%jspins,n,&
                 atoms%dx(ntyp),rad,rhoss, vxc)
         ELSE
            CALL  vxcall(6,xcpot,input%jspins,&
                 SIZE(vx,1),jrc,rhoss,&
                 vx,vxc)
            
231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295
         ENDIF
         DO ispin = 1, input%jspins
           DO i = 1,n
             vr2(i,ispin) = vr(i) + vxc(i,ispin)*rad(i)
           ENDDO
         ENDDO
!-ta
!        determine distance of potentials
!
         r3 = rn**3
         dist = 0.0
         DO ispin = 1, input%jspins
           DO i = 1,n
              a(i) = (vr2(i,ispin)-vr1(i,ispin))**2
           ENDDO
           CALL intgr1(a,rnot,h,n,b)
           dist = dist + sqrt((3.0e0/r3)*b(n))
         ENDDO
         IF (lastit) GO TO 190
         IF (dist.LT.distol) lastit = .true.
!     mix new input potential
         p1 = 1.0e0/dist
         p = min(pmax,p1)
         p = max(p,pmin)
         WRITE (6,FMT=8060) it,dist,p
         p1 = 1.0e0 - p
         DO ispin = 1, input%jspins
           DO i = 1,n
             vr1(i,ispin) = p1*vr1(i,ispin) + p*vr2(i,ispin)
           ENDDO
         ENDDO
  180 CONTINUE
!
! output 
!
      WRITE (6,FMT=8030) dist
      conv = .false.
!     list eigenvalues
  190 IF (conv) WRITE (6,FMT=8040) it,dist
      DO ispin = 1,input%jspins
        WRITE (6,'(a8,i2)') 'spin No.',ispin 
        DO k = 1,nst
           fj = iabs(kappa(k)) - 0.5e0
           l = fj + 0.5e0*isign(1,kappa(k)) + 0.01e0
           lnum(k) = l
           WRITE (6,FMT=8050) nprnc(k),kappa(k),l,fj,&
     &                        occ(k,ispin),eig(k,ispin),rev(k,ispin)
        ENDDO  
!
!--->   guess enpara if it doesn't exist, using floating energy parameters
!
        i = atoms%jri(ntyp) - (log(4.0)/atoms%dx(ntyp)+1.51)
        vbar(ispin) = vr1(i,ispin)/( rnot*exp(atoms%dx(ntyp)*(i-1)) )
        WRITE ( 6,'(/,'' reference energy = '',2f12.6,/)') vbar(ispin)
      ENDDO

 8030 FORMAT (/,/,/,' $$$ error: not converged, dist=',f10.6,/)
 8040 FORMAT (/,/,3x,'converged in',i4,' iterations to a distance of',&
     &       e12.5,' har',/,/,3x,'n  kappa  l    j  ',5x,&
     &       'occ.   eigenvalue (har)  <r>  ',/)
 8050 FORMAT (3x,i1,i5,i5,f6.1,2 (3x,f7.2,1x,2f12.6))
 8060 FORMAT ('it,dist,p=',i4,2f12.5)

      END SUBROUTINE atom2
      END MODULE m_atom2