atom2.f90 9.58 KB
Newer Older
Matthias Redies's avatar
Matthias Redies committed
1 2
MODULE m_atom2
   use m_juDFT
3 4 5 6 7 8 9
!     *************************************************************
!     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
!     *************************************************************
Matthias Redies's avatar
Matthias Redies committed
10 11
CONTAINS
   SUBROUTINE atom2(&
Matthias Redies's avatar
Matthias Redies committed
12
  &                 dimension, atoms, xcpot, input, ntyp, jrc, rnot1,&
Matthias Redies's avatar
Matthias Redies committed
13
  &                 qdel,&
Matthias Redies's avatar
Matthias Redies committed
14
  &                 rhoss, nst, lnum, eig, vbar)
15

Matthias Redies's avatar
Matthias Redies committed
16
      USE m_intgr, ONLY: intgr1, intgr0
17 18 19 20 21 22 23 24 25
      USE m_constants
      USE m_potl0
      USE m_stpot1
      USE m_setcor
      USE m_differ
      USE m_types
      IMPLICIT NONE
!     ..
!     .. Scalar Arguments ..
Matthias Redies's avatar
Matthias Redies committed
26 27 28 29 30 31 32 33 34
      TYPE(t_dimension), INTENT(IN)  :: dimension
      TYPE(t_atoms), INTENT(IN)      :: atoms
      CLASS(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,input%jspins),
      REAL, INTENT(OUT) :: eig(dimension%nstd, input%jspins), vbar(input%jspins)
      INTEGER, INTENT(OUT) :: nst, lnum(dimension%nstd)
35 36
!     ..
!     .. Local Scalars ..
Matthias Redies's avatar
Matthias Redies committed
37 38 39 40
      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
41 42
!     ..
!     .. Local Arrays ..
Matthias Redies's avatar
Matthias Redies committed
43 44 45 46 47
      REAL a(jrc), b(jrc), dens(jrc), occ(dimension%nstd, input%jspins)
      REAL rad(jrc), rev(dimension%nstd, input%jspins), ahelp(jrc), ain(jrc),&
     &     rh(jrc), vr(jrc), f(0:3),&
     &     vr1(jrc, input%jspins), vr2(jrc, input%jspins), vx(dimension%msh, input%jspins), vxc(dimension%msh, input%jspins)
      INTEGER kappa(dimension%nstd), nprnc(dimension%nstd)
48 49 50 51
!     ..
!     ..
!     .. Data statements ..
!---->     distol set from 1.0e-6 to 1.0e-3
Matthias Redies's avatar
Matthias Redies committed
52
      DATA zero, distol/0.0e0, 1.0e-3/
53 54
!     ..
      c = c_light(1.0)
Matthias Redies's avatar
Matthias Redies committed
55 56
      vxc(:, :) = 0.0
      vx(:, :) = 0.0
57
!
Matthias Redies's avatar
Matthias Redies committed
58 59
      WRITE (6, FMT=8000)
8000  FORMAT(' subroutine atom2 entered')
60 61 62 63 64 65 66 67 68
      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
Matthias Redies's avatar
Matthias Redies committed
69
      DO i = 1, n
70 71 72 73 74
         rad(i) = r
         r = r*d
      enddo
      rn = rad(n)
      bmu_l = atoms%bmu(ntyp)
Matthias Redies's avatar
Matthias Redies committed
75
      IF (bmu_l > 0.001 .AND. atoms%numStatesProvided(ntyp) .NE. 0) CALL &
Matthias Redies's avatar
Matthias Redies committed
76
         judft_warn("You specified both: inital moment and occupation numbers.", &
Matthias Redies's avatar
Matthias Redies committed
77 78
                    hint="The inital moment will be ignored, set magMom=0.0", calledby="atom2.f90")
      CALL setcor(ntyp, input%jspins, atoms, input, bmu_l, nst, kappa, nprnc, occ)
79 80 81 82 83

!
!--->   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
Matthias Redies's avatar
Matthias Redies committed
84 85
      IF (input%jspins == 1) THEN
         occ(nst, 1) = occ(nst, 1) + qdel
86
      ELSE
Matthias Redies's avatar
Matthias Redies committed
87 88
         occ(nst, 1) = occ(nst, 1) + qdel/2.
         occ(nst, input%jspins) = occ(nst, input%jspins) + qdel/2.
89 90 91
      ENDIF
!
      CALL stpot1(&
Matthias Redies's avatar
Matthias Redies committed
92
     &           jrc, n, z, rad,&
93
     &           vr1)
Matthias Redies's avatar
Matthias Redies committed
94 95
      DO i = 1, n
         vr1(i, input%jspins) = vr1(i, 1)
96 97 98 99 100 101 102 103
      ENDDO
!
!     start iterating
!
      lastit = .false.
      conv = .true.
      delrv = 0.100
      inr0 = log(5.0/rnot)/h + 1.5
Matthias Redies's avatar
Matthias Redies committed
104
      DO it = 1, itmax
Matthias Redies's avatar
Matthias Redies committed
105
         DO ispin = 1, input%jspins
106 107
!
!---->     load potential
Matthias Redies's avatar
Matthias Redies committed
108 109
            DO i = 1, n
               vr(i) = vr1(i, ispin)
Matthias Redies's avatar
Matthias Redies committed
110
            ENDDO
111
!----> adding stabilizing well: m. weinert
Matthias Redies's avatar
Matthias Redies committed
112 113
            DO i = inr0, n
               vr(i) = vr(i) + rad(i)*delrv*(rad(i) - rad(inr0))
Matthias Redies's avatar
Matthias Redies committed
114
            ENDDO
115
!---->     note that vr contains r*v(r) in hartree units
Matthias Redies's avatar
Matthias Redies committed
116 117
            DO i = 1, n
               rhoss(i, ispin) = zero
Matthias Redies's avatar
Matthias Redies committed
118 119
            ENDDO
            if (lastit) THEN
Matthias Redies's avatar
Matthias Redies committed
120
               inquire (file="startcharges", exist=l_start)
Matthias Redies's avatar
Matthias Redies committed
121
               if (l_start) then
Matthias Redies's avatar
Matthias Redies committed
122
                  OPEN (61, file="startcharges")
Matthias Redies's avatar
Matthias Redies committed
123
                  DO WHILE (.true.)
Matthias Redies's avatar
Matthias Redies committed
124 125 126
                     read (61, *, end=888, err=888) i, rho
                     if (i == z) then
                        occ(nst, 1) = occ(nst, 1) + rho
Matthias Redies's avatar
Matthias Redies committed
127 128 129 130
                        goto 888
                     endif
                  enddo
888               continue
Matthias Redies's avatar
Matthias Redies committed
131
                  close (61)
Matthias Redies's avatar
Matthias Redies committed
132 133
               endif
            endif
Matthias Redies's avatar
Matthias Redies committed
134
            DO 90 k = 1, nst
Matthias Redies's avatar
Matthias Redies committed
135 136
               fn = nprnc(k)
               fj = iabs(kappa(k)) - 0.5e0
Matthias Redies's avatar
Matthias Redies committed
137 138
               fl = fj + 0.5e0*isign(1, kappa(k))
               e = -2*(z/(fn + fl))**2
Matthias Redies's avatar
Matthias Redies committed
139 140
               ierr = -1
               msh_l = jrc
Matthias Redies's avatar
Matthias Redies committed
141
               DO WHILE (ierr .NE. 0)
Matthias Redies's avatar
Matthias Redies committed
142
                  CALL differ(&
Matthias Redies's avatar
Matthias Redies committed
143
       &                      fn, fl, fj, c, z, h, rnot, rn, d, msh_l, vr,&
Matthias Redies's avatar
Matthias Redies committed
144
       &                      e,&
Matthias Redies's avatar
Matthias Redies committed
145
       &                      a, b, ierr)!keep
Matthias Redies's avatar
Matthias Redies committed
146
                  msh_l = msh_l - 1
Matthias Redies's avatar
Matthias Redies committed
147 148
                  IF (jrc - msh_l > 100) CALL juDFT_error(&
       &               "atom2", calledby="atom2")
Matthias Redies's avatar
Matthias Redies committed
149
               ENDDO
Matthias Redies's avatar
Matthias Redies committed
150
               DO i = msh_l + 1, jrc
Matthias Redies's avatar
Matthias Redies committed
151 152 153
                  a(i) = a(msh_l)
                  b(i) = b(msh_l)
               ENDDO
154

Matthias Redies's avatar
Matthias Redies committed
155 156
               DO i = 1, n
                  rh(i) = occ(k, ispin)*(a(i)**2 + b(i)**2)
Matthias Redies's avatar
Matthias Redies committed
157
               ENDDO
158
!+ldau
Matthias Redies's avatar
Matthias Redies committed
159 160 161 162 163 164
               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
Matthias Redies's avatar
Matthias Redies committed
165 166
                        ain(i) = a(i)**2*r**(-kk - 1)      ! prepare inner integrand
                        r = r*d
Matthias Redies's avatar
Matthias Redies committed
167
                     ENDDO
Matthias Redies's avatar
Matthias Redies committed
168
                     CALL intgr1(ain, rnot, h, n, &          ! integrate&
Matthias Redies's avatar
Matthias Redies committed
169 170
                                &ahelp)
                     r = rnot
Matthias Redies's avatar
Matthias Redies committed
171 172 173
                     DO i = 1, n - 1
                        ain(i) = a(i)**2*r**kk*(ahelp(n) - ahelp(i))
                        r = r*d
Matthias Redies's avatar
Matthias Redies committed
174
                     ENDDO
Matthias Redies's avatar
Matthias Redies committed
175
                     CALL intgr0(ain, rnot, h, n - 1,      &     ! integrate 2nd r&
Matthias Redies's avatar
Matthias Redies committed
176
                                &f(kk/2))
177

Matthias Redies's avatar
Matthias Redies committed
178
                  ENDDO
179
!                 write(*,*) (hartree_to_ev_const*2*f(kk),kk=0,l)
Matthias Redies's avatar
Matthias Redies committed
180
               ENDIF
181
!-ldau
Matthias Redies's avatar
Matthias Redies committed
182
               eig(k, ispin) = e
183
!---->       calculate <r>
Matthias Redies's avatar
Matthias Redies committed
184 185
               DO i = 1, n
                  a(i) = (a(i)**2 + b(i)**2)*rad(i)
Matthias Redies's avatar
Matthias Redies committed
186
               ENDDO
Matthias Redies's avatar
Matthias Redies committed
187 188 189 190
               CALL intgr1(a, rnot, h, n, b)
               rev(k, ispin) = b(n)
               DO i = 1, n
                  rhoss(i, ispin) = rhoss(i, ispin) + rh(i)
Matthias Redies's avatar
Matthias Redies committed
191 192
               ENDDO
90          ENDDO
193 194 195 196
         ENDDO
!
!     solve poisson's equation
!
Matthias Redies's avatar
Matthias Redies committed
197 198
         DO i = 1, n
            dens(i) = rhoss(i, 1)
199
         ENDDO
Matthias Redies's avatar
Matthias Redies committed
200 201 202
         IF (input%jspins == 2) THEN
            DO i = 1, n
               dens(i) = dens(i) + rhoss(i, input%jspins)
Matthias Redies's avatar
Matthias Redies committed
203
            ENDDO
204
         ENDIF
Matthias Redies's avatar
Matthias Redies committed
205 206
         CALL intgr1(dens, rnot, h, n, a)
         DO i = 1, n
207
            rh(i) = dens(i)/rad(i)
Matthias Redies's avatar
Matthias Redies committed
208
         ENDDO
Matthias Redies's avatar
Matthias Redies committed
209 210 211 212 213
         CALL intgr1(rh, rnot, h, n, b)
         fisr = b(n)
         DO i = 1, n
            vr(i) = (a(i) + rad(i)*(fisr - b(i)) - z)
         ENDDO
214
!+ta
Matthias Redies's avatar
Matthias Redies committed
215 216 217 218 219
         DO ispin = 1, input%jspins
            DO i = 1, n
               rhoss(i, ispin) = rhoss(i, ispin)/(fpi_const*rad(i)**2)
            ENDDO
         ENDDO
Matthias Redies's avatar
Matthias Redies committed
220
         IF (xcpot%needs_grad()) THEN
Matthias Redies's avatar
Matthias Redies committed
221 222 223 224 225 226 227 228 229
            CALL potl0(xcpot, input%jspins, atoms%dx(ntyp), rad, rhoss, vxc)
         ELSE
            CALL xcpot%get_vxc(input%jspins, rhoss, vxc, vx)
         ENDIF
         DO ispin = 1, input%jspins
            DO i = 1, n
               vr2(i, ispin) = vr(i) + vxc(i, ispin)*rad(i)
            ENDDO
         ENDDO
230 231 232
!-ta
!        determine distance of potentials
!
Matthias Redies's avatar
Matthias Redies committed
233 234 235 236 237 238 239 240 241 242 243
         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 < distol) lastit = .true.
244
!     mix new input potential
Matthias Redies's avatar
Matthias Redies committed
245 246 247 248 249 250 251 252 253 254
         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
Matthias Redies's avatar
Matthias Redies committed
255
      ENDDO
256
!
Matthias Redies's avatar
Matthias Redies committed
257
! output
258
!
Matthias Redies's avatar
Matthias Redies committed
259 260
      WRITE (6, FMT=8030) dist
      conv = .false.
261
!     list eigenvalues
Matthias Redies's avatar
Matthias Redies committed
262 263 264 265 266 267 268 269 270 271
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
272 273 274
!
!--->   guess enpara if it doesn't exist, using floating energy parameters
!
Matthias Redies's avatar
Matthias Redies committed
275 276 277 278
         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
279

Matthias Redies's avatar
Matthias Redies committed
280 281 282 283 284 285
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)
286

Matthias Redies's avatar
Matthias Redies committed
287 288
   END SUBROUTINE atom2
END MODULE m_atom2