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

          Line data    Source code
       1             : MODULE m_atom2
       2             :    use m_juDFT
       3             : !     *************************************************************
       4             : !     fully relativistic atomic program based on the subroutines
       5             : !     differ, outint and inwint by d.d.koelling
       6             : !     erich wimmer     august 1981
       7             : !     modified for use in start-up generator.  m.w. april 1982
       8             : !     modified by adding stabilizing well. m. weinert may 1990
       9             : !     *************************************************************
      10             : CONTAINS
      11          40 :    SUBROUTINE atom2(&
      12             :   &                 dimension, atoms, xcpot, input, ntyp, jrc, rnot1,&
      13             :   &                 qdel,&
      14          40 :   &                 rhoss, nst, lnum, eig, vbar)
      15             : 
      16             :       USE m_intgr, ONLY: intgr1, intgr0
      17             :       USE m_constants
      18             :       USE m_potl0
      19             :       USE m_stpot1
      20             :       USE m_setcor
      21             :       USE m_differ
      22             :       USE m_types
      23             :       IMPLICIT NONE
      24             : !     ..
      25             : !     .. Scalar Arguments ..
      26             :       TYPE(t_dimension), INTENT(IN)  :: dimension
      27             :       TYPE(t_atoms), INTENT(IN)      :: atoms
      28             :       CLASS(t_xcpot), INTENT(IN)     :: xcpot
      29             :       TYPE(t_input), INTENT(IN)      :: input
      30             :       INTEGER, INTENT(IN)  :: jrc, ntyp
      31             :       REAL, INTENT(IN)  :: rnot1, qdel
      32             :       REAL, INTENT(OUT) :: rhoss(:, :) !(mshd,input%jspins),
      33             :       REAL, INTENT(OUT) :: eig(dimension%nstd, input%jspins), vbar(input%jspins)
      34             :       INTEGER, INTENT(OUT) :: nst, lnum(dimension%nstd)
      35             : !     ..
      36             : !     .. Local Scalars ..
      37             :       REAL c, d, delrv, dist, distol, e, fisr, fj, fl, fn, h,&
      38             :      &     p, p1, pmax, pmin, r, r3, rn, rnot, z, zero, bmu_l, rho
      39             :       INTEGER i, inr0, it, itmax, k, l, n, ispin, kk, ierr, msh_l
      40             :       LOGICAL conv, lastit, l_start
      41             : !     ..
      42             : !     .. Local Arrays ..
      43         240 :       REAL a(jrc), b(jrc), dens(jrc), occ(dimension%nstd, input%jspins)
      44         160 :       REAL rad(jrc), rev(dimension%nstd, input%jspins), ahelp(jrc), ain(jrc),&
      45          80 :      &     rh(jrc), vr(jrc), f(0:3),&
      46          80 :      &     vr1(jrc, input%jspins), vr2(jrc, input%jspins), vx(dimension%msh, input%jspins), vxc(dimension%msh, input%jspins)
      47          80 :       INTEGER kappa(dimension%nstd), nprnc(dimension%nstd)
      48             : !     ..
      49             : !     ..
      50             : !     .. Data statements ..
      51             : !---->     distol set from 1.0e-6 to 1.0e-3
      52             :       DATA zero, distol/0.0e0, 1.0e-3/
      53             : !     ..
      54          40 :       c = c_light(1.0)
      55          93 :       vxc(:, :) = 0.0
      56          93 :       vx(:, :) = 0.0
      57             : !
      58          40 :       WRITE (6, FMT=8000)
      59             : 8000  FORMAT(' subroutine atom2 entered')
      60          40 :       z = atoms%zatom(ntyp)
      61          40 :       n = jrc
      62          40 :       rnot = rnot1
      63          40 :       itmax = 100
      64          40 :       pmin = 0.01e0
      65          40 :       pmax = 0.2e0
      66          40 :       h = atoms%dx(ntyp)
      67          40 :       d = exp(h)
      68          40 :       r = rnot
      69       32358 :       DO i = 1, n
      70       32318 :          rad(i) = r
      71       32358 :          r = r*d
      72             :       enddo
      73          40 :       rn = rad(n)
      74          40 :       bmu_l = atoms%bmu(ntyp)
      75          40 :       IF (bmu_l > 0.001 .AND. atoms%numStatesProvided(ntyp) .NE. 0) CALL &
      76             :          judft_warn("You specified both: inital moment and occupation numbers.", &
      77           0 :                     hint="The inital moment will be ignored, set magMom=0.0", calledby="atom2.f90")
      78          40 :       CALL setcor(ntyp, input%jspins, atoms, input, bmu_l, nst, kappa, nprnc, occ)
      79             : 
      80             : !
      81             : !--->   for electric field case (sigma.ne.0), add the extra charge
      82             : !--->   to the uppermost level; ignore the possible problem that
      83             : !--->   the occupations may not be between 0 and 2
      84          40 :       IF (input%jspins == 1) THEN
      85          27 :          occ(nst, 1) = occ(nst, 1) + qdel
      86             :       ELSE
      87          13 :          occ(nst, 1) = occ(nst, 1) + qdel/2.
      88          13 :          occ(nst, input%jspins) = occ(nst, input%jspins) + qdel/2.
      89             :       ENDIF
      90             : !
      91             :       CALL stpot1(&
      92             :      &           jrc, n, z, rad,&
      93          40 :      &           vr1)
      94       32358 :       DO i = 1, n
      95       32358 :          vr1(i, input%jspins) = vr1(i, 1)
      96             :       ENDDO
      97             : !
      98             : !     start iterating
      99             : !
     100          40 :       lastit = .false.
     101          40 :       conv = .true.
     102          40 :       delrv = 0.100
     103          40 :       inr0 = log(5.0/rnot)/h + 1.5
     104         924 :       DO it = 1, itmax
     105        2171 :          DO ispin = 1, input%jspins
     106             : !
     107             : !---->     load potential
     108      978425 :             DO i = 1, n
     109      978425 :                vr(i) = vr1(i, ispin)
     110             :             ENDDO
     111             : !----> adding stabilizing well: m. weinert
     112      103505 :             DO i = inr0, n
     113      103505 :                vr(i) = vr(i) + rad(i)*delrv*(rad(i) - rad(inr0))
     114             :             ENDDO
     115             : !---->     note that vr contains r*v(r) in hartree units
     116      978425 :             DO i = 1, n
     117      978425 :                rhoss(i, ispin) = zero
     118             :             ENDDO
     119        1247 :             if (lastit) THEN
     120          53 :                inquire (file="startcharges", exist=l_start)
     121          53 :                if (l_start) then
     122           0 :                   OPEN (61, file="startcharges")
     123             :                   DO WHILE (.true.)
     124           0 :                      read (61, *, end=888, err=888) i, rho
     125           0 :                      if (i == z) then
     126           0 :                         occ(nst, 1) = occ(nst, 1) + rho
     127           0 :                         goto 888
     128             :                      endif
     129             :                   enddo
     130             : 888               continue
     131           0 :                   close (61)
     132             :                endif
     133             :             endif
     134       15074 :             DO 90 k = 1, nst
     135       12903 :                fn = nprnc(k)
     136       12903 :                fj = iabs(kappa(k)) - 0.5e0
     137       12903 :                fl = fj + 0.5e0*isign(1, kappa(k))
     138       12903 :                e = -2*(z/(fn + fl))**2
     139       12903 :                ierr = -1
     140       12903 :                msh_l = jrc
     141       25806 :                DO WHILE (ierr .NE. 0)
     142             :                   CALL differ(&
     143             :        &                      fn, fl, fj, c, z, h, rnot, rn, d, msh_l, vr,&
     144             :        &                      e,&
     145       12903 :        &                      a, b, ierr)!keep
     146       12903 :                   msh_l = msh_l - 1
     147       12903 :                   IF (jrc - msh_l > 100) CALL juDFT_error(&
     148           0 :        &               "atom2", calledby="atom2")
     149             :                ENDDO
     150       25806 :                DO i = msh_l + 1, jrc
     151       12903 :                   a(i) = a(msh_l)
     152       25806 :                   b(i) = b(msh_l)
     153             :                ENDDO
     154             : 
     155    10652317 :                DO i = 1, n
     156    10652317 :                   rh(i) = occ(k, ispin)*(a(i)**2 + b(i)**2)
     157             :                ENDDO
     158             : !+ldau
     159       12903 :                IF (lastit) THEN                         ! calculate slater interals
     160         545 :                   l = int(fl)
     161             : !                 write(*,*) nprnc(k),l
     162        1546 :                   DO kk = 0, 2*l, 2                      ! F0 for s, F0 + F2 for p etc.
     163        1001 :                      r = rnot
     164      840463 :                      DO i = 1, n
     165      839462 :                         ain(i) = a(i)**2*r**(-kk - 1)      ! prepare inner integrand
     166      840463 :                         r = r*d
     167             :                      ENDDO
     168             :                      CALL intgr1(ain, rnot, h, n, &          ! integrate&
     169        1001 :                                 &ahelp)
     170        1001 :                      r = rnot
     171      839462 :                      DO i = 1, n - 1
     172      838461 :                         ain(i) = a(i)**2*r**kk*(ahelp(n) - ahelp(i))
     173      839462 :                         r = r*d
     174             :                      ENDDO
     175             :                      CALL intgr0(ain, rnot, h, n - 1,      &     ! integrate 2nd r&
     176        1546 :                                 &f(kk/2))
     177             : 
     178             :                   ENDDO
     179             : !                 write(*,*) (hartree_to_ev_const*2*f(kk),kk=0,l)
     180             :                ENDIF
     181             : !-ldau
     182       12903 :                eig(k, ispin) = e
     183             : !---->       calculate <r>
     184    10652317 :                DO i = 1, n
     185    10652317 :                   a(i) = (a(i)**2 + b(i)**2)*rad(i)
     186             :                ENDDO
     187       12903 :                CALL intgr1(a, rnot, h, n, b)
     188       12903 :                rev(k, ispin) = b(n)
     189    10652317 :                DO i = 1, n
     190    10652317 :                   rhoss(i, ispin) = rhoss(i, ispin) + rh(i)
     191             :                ENDDO
     192        1247 : 90          ENDDO
     193             :          ENDDO
     194             : !
     195             : !     solve poisson's equation
     196             : !
     197      745398 :          DO i = 1, n
     198      745398 :             dens(i) = rhoss(i, 1)
     199             :          ENDDO
     200         924 :          IF (input%jspins == 2) THEN
     201      233027 :             DO i = 1, n
     202      233027 :                dens(i) = dens(i) + rhoss(i, input%jspins)
     203             :             ENDDO
     204             :          ENDIF
     205         924 :          CALL intgr1(dens, rnot, h, n, a)
     206      745398 :          DO i = 1, n
     207      745398 :             rh(i) = dens(i)/rad(i)
     208             :          ENDDO
     209         924 :          CALL intgr1(rh, rnot, h, n, b)
     210         924 :          fisr = b(n)
     211      745398 :          DO i = 1, n
     212      745398 :             vr(i) = (a(i) + rad(i)*(fisr - b(i)) - z)
     213             :          ENDDO
     214             : !+ta
     215        2171 :          DO ispin = 1, input%jspins
     216      979349 :             DO i = 1, n
     217      978425 :                rhoss(i, ispin) = rhoss(i, ispin)/(fpi_const*rad(i)**2)
     218             :             ENDDO
     219             :          ENDDO
     220         924 :          IF (xcpot%needs_grad()) THEN
     221         704 :             CALL potl0(xcpot, input%jspins, atoms%dx(ntyp), rad, rhoss, vxc)
     222             :          ELSE
     223         220 :             CALL xcpot%get_vxc(input%jspins, rhoss, vxc, vx)
     224             :          ENDIF
     225        2171 :          DO ispin = 1, input%jspins
     226      979349 :             DO i = 1, n
     227      978425 :                vr2(i, ispin) = vr(i) + vxc(i, ispin)*rad(i)
     228             :             ENDDO
     229             :          ENDDO
     230             : !-ta
     231             : !        determine distance of potentials
     232             : !
     233         924 :          r3 = rn**3
     234         924 :          dist = 0.0
     235        2171 :          DO ispin = 1, input%jspins
     236      978425 :             DO i = 1, n
     237      978425 :                a(i) = (vr2(i, ispin) - vr1(i, ispin))**2
     238             :             ENDDO
     239        1247 :             CALL intgr1(a, rnot, h, n, b)
     240        2171 :             dist = dist + sqrt((3.0e0/r3)*b(n))
     241             :          ENDDO
     242         924 :          IF (lastit) GO TO 190
     243         884 :          IF (dist < distol) lastit = .true.
     244             : !     mix new input potential
     245         884 :          p1 = 1.0e0/dist
     246         884 :          p = min(pmax, p1)
     247         884 :          p = max(p, pmin)
     248         884 :          WRITE (6, FMT=8060) it, dist, p
     249         884 :          p1 = 1.0e0 - p
     250        2078 :          DO ispin = 1, input%jspins
     251      937558 :             DO i = 1, n
     252      936674 :                vr1(i, ispin) = p1*vr1(i, ispin) + p*vr2(i, ispin)
     253             :             ENDDO
     254             :          ENDDO
     255             :       ENDDO
     256             : !
     257             : ! output
     258             : !
     259           0 :       WRITE (6, FMT=8030) dist
     260           0 :       conv = .false.
     261             : !     list eigenvalues
     262          40 : 190   IF (conv) WRITE (6, FMT=8040) it, dist
     263          93 :       DO ispin = 1, input%jspins
     264          53 :          WRITE (6, '(a8,i2)') 'spin No.', ispin
     265         598 :          DO k = 1, nst
     266         545 :             fj = iabs(kappa(k)) - 0.5e0
     267         545 :             l = fj + 0.5e0*isign(1, kappa(k)) + 0.01e0
     268         545 :             lnum(k) = l
     269         545 :             WRITE (6, FMT=8050) nprnc(k), kappa(k), l, fj,&
     270        1143 :       &                        occ(k, ispin), eig(k, ispin), rev(k, ispin)
     271             :          ENDDO
     272             : !
     273             : !--->   guess enpara if it doesn't exist, using floating energy parameters
     274             : !
     275          53 :          i = atoms%jri(ntyp) - (log(4.0)/atoms%dx(ntyp) + 1.51)
     276          53 :          vbar(ispin) = vr1(i, ispin)/(rnot*exp(atoms%dx(ntyp)*(i - 1)))
     277          93 :          WRITE (6, '(/,'' reference energy = '',2f12.6,/)') vbar(ispin)
     278             :       ENDDO
     279             : 
     280             : 8030  FORMAT(/, /, /, ' $$$ error: not converged, dist=', f10.6,/)
     281             : 8040  FORMAT(/, /, 3x, 'converged in', i4, ' iterations to a distance of',&
     282             :       &       e12.5, ' har', /, /, 3x, 'n  kappa  l    j  ', 5x,&
     283             :       &       'occ.   eigenvalue (har)  <r>  ',/)
     284             : 8050  FORMAT(3x, i1, i5, i5, f6.1, 2(3x, f7.2, 1x, 2f12.6))
     285             : 8060  FORMAT('it,dist,p=', i4, 2f12.5)
     286             : 
     287          40 :    END SUBROUTINE atom2
     288             : END MODULE m_atom2

Generated by: LCOV version 1.13