LCOV - code coverage report
Current view: top level - hybrid - read_core.F90 (source / functions) Hit Total Coverage
Test: combined.info Lines: 0 230 0.0 %
Date: 2019-09-08 04:53:50 Functions: 0 4 0.0 %

          Line data    Source code
       1             : MODULE m_read_core
       2             : 
       3             :    ! read core radial wavefunctions from corebas
       4             :    ! corebas is written in cored.F
       5             :    ! (core basis functions can be read in once during an iteration)
       6             : 
       7             : CONTAINS
       8             : 
       9           0 :    SUBROUTINE read_core(atoms, hybdat)
      10             : 
      11             :       USE m_types
      12             :       IMPLICIT NONE
      13             : 
      14             :       TYPE(t_hybdat), INTENT(INOUT)   :: hybdat
      15             :       TYPE(t_atoms), INTENT(IN)       :: atoms
      16             : 
      17             :       ! - local scalars -
      18             :       INTEGER                            ::  ncst
      19             :       INTEGER                            ::  ok, itype, i, idum, l
      20             : 
      21             :       REAL                               ::  rdum
      22             :       REAL                               ::  weight1, weight2
      23             :       ! - local arrays -
      24           0 :       INTEGER, ALLOCATABLE                ::  nindxcr(:, :)
      25           0 :       INTEGER, ALLOCATABLE                ::  l_qn(:, :), n_qn(:, :)
      26             : 
      27           0 :       REAL, ALLOCATABLE                   ::  j_qn(:, :)
      28           0 :       REAL, ALLOCATABLE                   ::  core1r(:, :, :, :),&
      29           0 :      &                                       core2r(:, :, :, :)
      30             : 
      31           0 :       OPEN (UNIT=77, FILE='corebas', FORM='unformatted')
      32           0 :       READ (77) ncst
      33             : 
      34             :       ALLOCATE (n_qn(0:ncst, atoms%ntype), l_qn(0:ncst, atoms%ntype), &
      35           0 :      &          j_qn(0:ncst, atoms%ntype), stat=ok)
      36           0 :       IF (ok /= 0) STOP 'mhsfock: failure allocation n_qn,l_qn,j_qn'
      37           0 :       ALLOCATE (nindxcr(0:ncst, atoms%ntype), stat=ok)
      38           0 :       IF (ok /= 0) STOP 'mhsfock: failure allocation nindxcr'
      39             : 
      40           0 :       nindxcr = 0
      41           0 :       hybdat%lmaxc = 0
      42           0 :       l_qn = 0
      43           0 :       ncst = 0
      44           0 :       DO itype = 1, atoms%ntype
      45           0 :          READ (77) ncst
      46           0 :          IF (ncst == 0) CYCLE
      47           0 :          DO i = 1, ncst
      48           0 :             READ (77) n_qn(i, itype), l_qn(i, itype), j_qn(i, itype)
      49           0 :             nindxcr(l_qn(i, itype), itype) = nindxcr(l_qn(i, itype), itype) + 1
      50             :          END DO
      51           0 :          hybdat%lmaxc(itype) = maxval(l_qn(:, itype))
      52             :       END DO
      53             : 
      54           0 :       ALLOCATE (core1r(atoms%jmtd, 0:maxval(hybdat%lmaxc), maxval(nindxcr), atoms%ntype))
      55           0 :       ALLOCATE (core2r(atoms%jmtd, 0:maxval(hybdat%lmaxc), maxval(nindxcr), atoms%ntype))
      56           0 :       core1r = 0
      57           0 :       core2r = 0
      58           0 :       REWIND (77)
      59             : 
      60           0 :       DEALLOCATE (nindxcr)
      61           0 :       ALLOCATE (nindxcr(0:maxval(l_qn), atoms%ntype))
      62           0 :       nindxcr = 0
      63           0 :       READ (77) idum
      64           0 :       DO itype = 1, atoms%ntype
      65           0 :          READ (77) ncst
      66           0 :          DO i = 1, ncst
      67           0 :             nindxcr(l_qn(i, itype), itype) = nindxcr(l_qn(i, itype), itype) + 1
      68           0 :             READ (77) idum, idum, rdum, core1r(:atoms%jri(itype), l_qn(i, itype),&
      69           0 :        &                            nindxcr(l_qn(i, itype), itype), itype),&
      70             :        &                            core2r(:atoms%jri(itype), l_qn(i, itype),&
      71           0 :        &                            nindxcr(l_qn(i, itype), itype), itype)
      72             :          END DO
      73             :       END DO
      74             : 
      75           0 :       ALLOCATE (hybdat%nindxc(0:maxval(hybdat%lmaxc), atoms%ntype), stat=ok)
      76           0 :       IF (ok /= 0) STOP 'mhsfock: failure allocation nindxc'
      77             :       ALLOCATE (hybdat%core1(atoms%jmtd, 0:maxval(hybdat%lmaxc), maxval(nindxcr(0, :),&
      78           0 :      &                nint((maxval(nindxcr)/2.0))), atoms%ntype), stat=ok)
      79           0 :       IF (ok /= 0) STOP 'mhsfock: failure allocation core1'
      80             :       ALLOCATE (hybdat%core2(atoms%jmtd, 0:maxval(hybdat%lmaxc), maxval(nindxcr(0, :),&
      81           0 :      &                nint((maxval(nindxcr)/2.0))), atoms%ntype), stat=ok)
      82           0 :       IF (ok /= 0) STOP 'mhsfock: failure allocation core2'
      83           0 :       hybdat%nindxc = 0; hybdat%core1 = 0; hybdat%core2 = 0
      84             : 
      85             :       ! average over core states that only differ in j
      86             :       ! core functions with l-qn equal 0 doesnot change during averaging
      87             : 
      88           0 :       hybdat%nindxc(0, :) = nindxcr(0, :)
      89           0 :       DO itype = 1, atoms%ntype
      90             :          hybdat%core1(:, 0, :hybdat%nindxc(0, itype), itype)&
      91           0 :       &    = core1r(:, 0, :hybdat%nindxc(0, itype), itype)
      92             :          hybdat%core2(:, 0, :hybdat%nindxc(0, itype), itype)&
      93           0 :       &    = core2r(:, 0, :hybdat%nindxc(0, itype), itype)
      94             :       END DO
      95             : 
      96           0 :       DO itype = 1, atoms%ntype
      97           0 :          DO l = 1, hybdat%lmaxc(itype)
      98           0 :             weight1 = 2*(l - 0.5) + 1
      99           0 :             weight2 = 2*(l + 0.5) + 1
     100           0 :             IF (modulo(nindxcr(l, itype), 2) == 0) THEN
     101           0 :                DO i = 1, nindxcr(l, itype), 2
     102           0 :                   hybdat%nindxc(l, itype) = hybdat%nindxc(l, itype) + 1
     103             :                   hybdat%core1(:atoms%jri(itype), l, hybdat%nindxc(l, itype), itype) =&
     104             :          &                        (weight1*core1r(:atoms%jri(itype), l, i, itype) +&
     105             :          &                         weight2*core1r(:atoms%jri(itype), l, i + 1, itype))&
     106           0 :          &                        /(weight1 + weight2)
     107             :                   hybdat%core2(:atoms%jri(itype), l, hybdat%nindxc(l, itype), itype) =&
     108             :          &                        (weight1*core2r(:atoms%jri(itype), l, i, itype) +&
     109             :          &                         weight2*core2r(:atoms%jri(itype), l, i + 1, itype))&
     110           0 :          &                        /(weight1 + weight2)
     111             :                END DO
     112             :             ELSE
     113           0 :                DO i = 1, nindxcr(l, itype) - 1, 2
     114           0 :                   hybdat%nindxc(l, itype) = hybdat%nindxc(l, itype) + 1
     115             :                   hybdat%core1(:atoms%jri(itype), l, hybdat%nindxc(l, itype), itype) =&
     116             :          &                        (weight1*core1r(:atoms%jri(itype), l, i, itype) +&
     117             :          &                         weight2*core1r(:atoms%jri(itype), l, i + 1, itype))&
     118           0 :          &                        /(weight1 + weight2)
     119             :                   hybdat%core2(:atoms%jri(itype), l, hybdat%nindxc(l, itype), itype) =&
     120             :          &                        (weight1*core2r(:atoms%jri(itype), l, i, itype) +&
     121             :          &                         weight2*core2r(:atoms%jri(itype), l, i + 1, itype))&
     122           0 :          &                        /(weight1 + weight2)
     123             :                END DO
     124           0 :                hybdat%nindxc(l, itype) = hybdat%nindxc(l, itype) + 1
     125             :                hybdat%core1(:atoms%jri(itype), l, hybdat%nindxc(l, itype), itype)&
     126           0 :       &          = core1r(:atoms%jri(itype), l, nindxcr(l, itype), itype)
     127             :                hybdat%core2(:atoms%jri(itype), l, hybdat%nindxc(l, itype), itype)&
     128           0 :       &          = core2r(:atoms%jri(itype), l, nindxcr(l, itype), itype)
     129             :             END IF
     130             : 
     131             :          END DO
     132             :       END DO
     133             : 
     134           0 :       DEALLOCATE (nindxcr, core1r, core2r)
     135             : 
     136           0 :       hybdat%maxindxc = maxval(hybdat%nindxc)
     137           0 :       CLOSE (77)
     138             : 
     139           0 :    END SUBROUTINE read_core
     140             : 
     141           0 :    SUBROUTINE corewf(atoms, jsp, input, dimension,&
     142           0 :   &                   vr, lmaxcd, maxindxc, mpi, lmaxc, nindxc, core1, core2, eig_c)
     143             :       USE m_types
     144             :       IMPLICIT NONE
     145             : 
     146             :       TYPE(t_mpi), INTENT(IN)   :: mpi
     147             :       TYPE(t_dimension), INTENT(IN)   :: dimension
     148             :       TYPE(t_input), INTENT(IN)   :: input
     149             :       TYPE(t_atoms), INTENT(IN)   :: atoms
     150             : 
     151             :       ! - scalars -
     152             :       INTEGER, INTENT(IN)        ::  jsp
     153             :       INTEGER, INTENT(IN)        :: lmaxcd
     154             :       INTEGER, INTENT(INOUT)     ::  maxindxc
     155             : 
     156             :       ! -arrays -
     157             :       INTEGER, INTENT(INOUT)     ::  lmaxc(:)
     158             :       INTEGER, INTENT(INOUT)     ::  nindxc(0:lmaxcd, atoms%ntype)
     159             : 
     160             :       REAL, INTENT(IN)       ::  vr(:, :, :)!(atoms%jmtd,atoms%ntypd,input%jspins)
     161             :       REAL, INTENT(INOUT)    ::  core1(:, :, 0:, :) !(atoms%jmtd,maxindxc,0:lmaxcd,atoms%ntype)
     162             :       REAL, INTENT(INOUT)     ::  core2(:, :, 0:, :) !(jmtd,maxindxc,0:lmaxcd,ntype)
     163             : 
     164             :       REAL, INTENT(INOUT)    ::  eig_c(maxindxc, 0:lmaxcd, atoms%ntype)
     165             : 
     166             :       ! - local scalars -
     167             :       INTEGER                   ::  ncstd
     168             :       INTEGER                   ::  ok, itype, i, j, idum, l
     169             : 
     170             :       REAL                      ::  rdum
     171             :       REAL                      ::  weight1, weight2
     172             :       ! - local arrays -
     173           0 :       INTEGER, ALLOCATABLE       ::  nindxcr(:, :)
     174             :       INTEGER, ALLOCATABLE       ::  l_qn(:, :), n_qn(:, :)
     175             : 
     176             :       REAL, ALLOCATABLE          ::  j_qn(:, :)
     177           0 :       REAL, ALLOCATABLE          ::  core1r(:, :, :, :), core2r(:, :, :, :)
     178             :       REAL, ALLOCATABLE          ::  core11r(:, :, :, :), core22r(:, :, :, :)
     179           0 :       REAL, ALLOCATABLE          ::  eig_cr(:, :, :)
     180             : 
     181           0 :       ncstd = maxval(atoms%ncst)
     182           0 :       ALLOCATE (nindxcr(0:ncstd, atoms%ntype), stat=ok)
     183             : 
     184             :       ! generate relativistic core wave functions( ->core1r,core2r )
     185             :       CALL calcorewf(dimension, input, jsp, atoms,&
     186             :      &                ncstd, vr,&
     187           0 :      &                lmaxc, nindxcr, core1r, core2r, eig_cr, mpi)
     188             : 
     189           0 :       nindxc = 0
     190             : 
     191             :       ! average over core states that only differ in j
     192             :       ! core functions with l-qn equal 0 doesnot change during the average process
     193             : 
     194           0 :       nindxc(0, :) = nindxcr(0, :)
     195           0 :       DO itype = 1, atoms%ntype
     196             :          core1(:, :nindxc(0, itype), 0, itype)&
     197           0 :       &    = core1r(:, 0, :nindxc(0, itype), itype)
     198             :          core2(:, :nindxc(0, itype), 0, itype)&
     199           0 :       &    = core2r(:, 0, :nindxc(0, itype), itype)
     200             :          eig_c(:nindxc(0, itype), 0, itype)&
     201           0 :       &    = eig_cr(0, :nindxc(0, itype), itype)
     202             :       END DO
     203             : 
     204           0 :       DO itype = 1, atoms%ntype
     205           0 :          DO l = 1, lmaxc(itype)
     206           0 :             weight1 = 2*(l - 0.5) + 1
     207           0 :             weight2 = 2*(l + 0.5) + 1
     208           0 :             IF (modulo(nindxcr(l, itype), 2) == 0) THEN
     209           0 :                DO i = 1, nindxcr(l, itype), 2
     210           0 :                   nindxc(l, itype) = nindxc(l, itype) + 1
     211             :                   core1(:atoms%jri(itype), nindxc(l, itype), l, itype) =&
     212             :          &                       (weight1*core1r(:atoms%jri(itype), l, i, itype) +&
     213             :          &                        weight2*core1r(:atoms%jri(itype), l, i + 1, itype))&
     214           0 :          &                       /(weight1 + weight2)
     215             :                   core2(:atoms%jri(itype), nindxc(l, itype), l, itype) =&
     216             :          &                       (weight1*core2r(:atoms%jri(itype), l, i, itype) + &
     217             :          &                        weight2*core2r(:atoms%jri(itype), l, i + 1, itype))&
     218           0 :          &                       /(weight1 + weight2)
     219             : 
     220             :                   eig_c(nindxc(l, itype), l, itype) =&
     221             :          &                       (weight1*eig_cr(l, i, itype) +&
     222             :          &                        weight2*eig_cr(l, i + 1, itype))&
     223           0 :          &                       /(weight1 + weight2)
     224             :                END DO
     225             :             ELSE
     226           0 :                DO i = 1, nindxcr(l, itype) - 1, 2
     227           0 :                   nindxc(l, itype) = nindxc(l, itype) + 1
     228             :                   core1(:atoms%jri(itype), nindxc(l, itype), l, itype) =&
     229             :          &                       (weight1*core1r(:atoms%jri(itype), l, i, itype) +&
     230             :          &                        weight2*core1r(:atoms%jri(itype), l, i + 1, itype))&
     231           0 :          &                       /(weight1 + weight2)
     232             :                   core2(:atoms%jri(itype), nindxc(l, itype), l, itype) =&
     233             :          &                       (weight1*core2r(:atoms%jri(itype), l, i, itype) + &
     234             :          &                        weight2*core2r(:atoms%jri(itype), l, i + 1, itype))&
     235           0 :          &                       /(weight1 + weight2)
     236             : 
     237             :                   eig_c(nindxc(l, itype), l, itype) =&
     238             :          &                       (weight1*eig_cr(l, i, itype) +&
     239             :          &                        weight2*eig_cr(l, i + 1, itype))&
     240           0 :          &                       /(weight1 + weight2)
     241             :                END DO
     242           0 :                nindxc(l, itype) = nindxc(l, itype) + 1
     243             :                core1(:atoms%jri(itype), nindxc(l, itype), l, itype)&
     244           0 :         &        = core1r(:atoms%jri(itype), l, nindxcr(l, itype), itype)
     245             :                core2(:atoms%jri(itype), nindxc(l, itype), l, itype)&
     246           0 :         &        = core2r(:atoms%jri(itype), l, nindxcr(l, itype), itype)
     247             :                eig_c(nindxc(l, itype), l, itype)&
     248           0 :         &        = eig_cr(l, nindxcr(l, itype), itype)
     249             :             END IF
     250             : 
     251             :          END DO
     252             :       END DO
     253             : 
     254           0 :       DEALLOCATE (nindxcr, core1r, core2r, eig_cr)
     255             : 
     256           0 :       IF (maxindxc /= maxval(nindxc))&
     257           0 :      &   STOP 'corewf: counting error nindxc'
     258             : 
     259           0 :    END SUBROUTINE corewf
     260             : 
     261           0 :    SUBROUTINE calcorewf(dimension, input, jspin, atoms,&
     262           0 :   &                     ncstd, vr,&
     263           0 :   &                     lmaxc, nindxcr, core1, core2, eig_c, mpi)
     264             : 
     265             :       USE m_intgr, ONLY: intgr3, intgr0, intgr1
     266             :       USE m_constants, ONLY: c_light
     267             :       USE m_setcor
     268             :       USE m_differ
     269             :       USE m_types
     270             :       IMPLICIT NONE
     271             : 
     272             :       TYPE(t_mpi), INTENT(IN)   :: mpi
     273             :       TYPE(t_dimension), INTENT(IN)   :: dimension
     274             :       TYPE(t_input), INTENT(IN)   :: input
     275             :       TYPE(t_atoms), INTENT(IN)   :: atoms
     276             : 
     277             :       !  - scalars -
     278             :       INTEGER, INTENT(IN) :: ncstd
     279             :       INTEGER, INTENT(IN) :: jspin
     280             :       INTEGER, INTENT(OUT):: lmaxc(:)
     281             : 
     282             :       !  - arrays -
     283             :       INTEGER, INTENT(OUT):: nindxcr(0:ncstd, atoms%ntype)
     284             :       REAL, INTENT(IN) :: vr(:, :, :)!(atoms%jmtd,atoms%ntypd,input%jspins)
     285             :       REAL, ALLOCATABLE :: core1(:, :, :, :), core2(:, :, :, :)
     286             :       REAL, ALLOCATABLE :: eig_c(:, :, :)
     287             : 
     288             :       !  - local scalars -
     289             :       INTEGER              :: i, j, itype, korb, ncmsh, nst, ierr
     290             :       REAL                 :: e, fj, fl, fn, t2, c, bmu, weight
     291             :       REAL                 :: d, dxx, rn, rnot, z, t1, rr
     292             :       LOGICAL, SAVE        :: first = .true.
     293             : 
     294             :       !  - local arrays -
     295           0 :       INTEGER              :: kappa(dimension%nstd), nprnc(dimension%nstd)
     296           0 :       REAL                 :: vrd(dimension%msh)
     297           0 :       REAL                 :: occ(dimension%nstd), occ_h(dimension%nstd, 2), a(dimension%msh), b(dimension%msh)
     298             :       REAL, ALLOCATABLE, SAVE:: vr0(:, :, :)
     299             : 
     300             :       !   - intrinsic functions -
     301             :       INTRINSIC exp, iabs, isign
     302             : 
     303           0 :       c = c_light(1.0)
     304             : 
     305           0 :       IF (first) THEN
     306           0 :          ALLOCATE (vr0(atoms%jmtd, atoms%ntype, input%jspins))
     307             :       END IF
     308             : 
     309           0 :       IF (input%frcor) THEN
     310           0 :          IF (first) THEN
     311           0 :             vr0 = vr
     312           0 :             first = .false.
     313             :          END IF
     314             :       ELSE
     315           0 :          vr0 = vr
     316           0 :          first = .false.
     317             :       END IF
     318             : 
     319             :       ! this loop determines the dimensions
     320             : 
     321           0 :       lmaxc = 0; nindxcr = 0
     322           0 :       DO itype = 1, atoms%ntype
     323           0 :          z = atoms%zatom(itype)
     324           0 :          dxx = atoms%dx(itype)
     325           0 :          bmu = 0.0
     326             :          CALL setcor(itype, input%jspins, atoms, input, bmu,&
     327           0 :       &               nst, kappa, nprnc, occ_h)
     328             : 
     329           0 :          IF ((bmu > 99.)) THEN
     330           0 :             occ(1:nst) = input%jspins*occ_h(1:nst, jspin)
     331             :          ELSE
     332           0 :             occ(1:nst) = occ_h(1:nst, 1)
     333             :          END IF
     334           0 :          rnot = atoms%rmsh(1, itype)
     335           0 :          d = exp(atoms%dx(itype))
     336           0 :          ncmsh = nint(log((atoms%rmt(itype) + 10.0)/rnot)/dxx + 1)
     337           0 :          ncmsh = min(ncmsh, dimension%msh)
     338           0 :          rn = rnot*(d**(ncmsh - 1))
     339             : 
     340           0 :          nst = atoms%ncst(itype)
     341             : 
     342           0 :          DO 80 korb = 1, nst
     343           0 :             IF (occ(korb) == 0) GOTO 80
     344           0 :             fn = nprnc(korb)
     345           0 :             fj = iabs(kappa(korb)) - .5e0
     346           0 :             weight = 2*fj + 1.e0
     347           0 :             IF (bmu > 99.) weight = occ(korb)
     348           0 :             fl = fj + (.5e0)*isign(1, kappa(korb))
     349           0 :             e = -2*(z/(fn + fl))**2
     350             : 
     351           0 :             nindxcr(NINT(fl), itype) = nindxcr(NINT(fl), itype) + 1
     352           0 :             lmaxc(itype) = max(lmaxc(itype), NINT(fl))
     353           0 : 80       END DO
     354             : 
     355             :       END DO
     356             : 
     357           0 :       ALLOCATE (core1(atoms%jmtd, 0:maxval(lmaxc), maxval(nindxcr), atoms%ntype))
     358           0 :       ALLOCATE (core2(atoms%jmtd, 0:maxval(lmaxc), maxval(nindxcr), atoms%ntype))
     359           0 :       ALLOCATE (eig_c(0:maxval(lmaxc), maxval(nindxcr), atoms%ntype))
     360             : 
     361           0 :       core1 = 0; core2 = 0
     362           0 :       nindxcr = 0
     363           0 :       DO itype = 1, atoms%ntype
     364           0 :          z = atoms%zatom(itype)
     365           0 :          dxx = atoms%dx(itype)
     366           0 :          bmu = 0.0
     367           0 :          CALL setcor(itype, input%jspins, atoms, input, bmu, nst, kappa, nprnc, occ_h)
     368             : 
     369           0 :          IF ((bmu > 99.)) THEN
     370           0 :             occ(1:nst) = input%jspins*occ_h(1:nst, jspin)
     371             :          ELSE
     372           0 :             occ(1:nst) = occ_h(1:nst, 1)
     373             :          END IF
     374           0 :          rnot = atoms%rmsh(1, itype)
     375           0 :          d = exp(atoms%dx(itype))
     376           0 :          ncmsh = nint(log((atoms%rmt(itype) + 10.0)/rnot)/dxx + 1)
     377           0 :          ncmsh = min(ncmsh, dimension%msh)
     378           0 :          rn = rnot*(d**(ncmsh - 1))
     379           0 :          IF (mpi%irank == 0) THEN
     380           0 :             WRITE (6, FMT=8000) z, rnot, dxx, atoms%jri(itype)
     381             :          END IF
     382           0 :          DO j = 1, atoms%jri(itype)
     383           0 :             vrd(j) = vr0(j, itype, jspin)
     384             :          END DO
     385             : 
     386           0 :          if (input%l_core_confpot) THEN
     387             : 
     388             :             ! linear extension of the potential with slope t1 / a.u.
     389           0 :             t1 = 0.125
     390           0 :             t2 = vrd(atoms%jri(itype))/atoms%rmt(itype) - atoms%rmt(itype)*t1
     391           0 :             rr = atoms%rmt(itype)
     392             :          else
     393           0 :             t2 = vrd(atoms%jri(itype))/(atoms%jri(itype) - ncmsh)
     394             :          endif
     395           0 :          IF (atoms%jri(itype) < ncmsh) THEN
     396           0 :             DO i = atoms%jri(itype) + 1, ncmsh
     397           0 :                if (input%l_core_confpot) THEN
     398           0 :                   rr = d*rr
     399           0 :                   vrd(i) = rr*(t2 + rr*t1)
     400             :                else
     401           0 :                   vrd(i) = vrd(atoms%jri(itype)) + t2*(i - atoms%jri(itype))
     402             :                endif
     403             : !
     404             :             END DO
     405             :          END IF
     406             : 
     407           0 :          nst = atoms%ncst(itype)
     408             : 
     409           0 :          DO 90 korb = 1, nst
     410           0 :             IF (occ(korb) == 0) GOTO 90
     411           0 :             fn = nprnc(korb)
     412           0 :             fj = iabs(kappa(korb)) - .5e0
     413           0 :             weight = 2*fj + 1.e0
     414           0 :             IF (bmu > 99.) weight = occ(korb)
     415           0 :             fl = fj + (.5e0)*isign(1, kappa(korb))
     416           0 :             e = -2*(z/(fn + fl))**2
     417           0 :             CALL differ(fn, fl, fj, c, z, dxx, rnot, rn, d, ncmsh, vrd, e, a, b, ierr)
     418             : 
     419           0 :             nindxcr(NINT(fl), itype) = nindxcr(NINT(fl), itype) + 1
     420             : 
     421             :             core1(:atoms%jri(itype), NINT(fl), nindxcr(NINT(fl), itype), itype)&
     422           0 :      &                                                 = a(:atoms%jri(itype))
     423             :             core2(:atoms%jri(itype), NINT(fl), nindxcr(NINT(fl), itype), itype)&
     424           0 :      &                                                 = b(:atoms%jri(itype))
     425             : 
     426           0 :             eig_c(NINT(fl), nindxcr(NINT(fl), itype), itype) = e
     427             : 
     428           0 :             IF (mpi%irank == 0) THEN
     429           0 :                WRITE (6, FMT=8010) fn, fl, fj, e, weight
     430             :             END IF
     431           0 :             IF (ierr /= 0) STOP 'error in core-level routine'
     432             : 
     433           0 : 90       END DO
     434             : 
     435             :       END DO
     436             : 
     437             : 8000  FORMAT(/, /, 10x, 'z=', f4.0, 5x, 'r(1)=', e14.6, 5x, 'dx=', f8.6, 5x,&
     438             :       &       'm.t.index=', i4, /, 15x, 'n', 4x, 'l', 5x, 'j', 4x, 'energy', 7x,&
     439             :       &       'weight')
     440             : 8010  FORMAT(12x, 2f5.0, f6.1, f10.4, f10.0)
     441             : 
     442           0 :    END SUBROUTINE calcorewf
     443             : 
     444           0 :    SUBROUTINE core_init(dimension, input, atoms, lmaxcd, maxindxc)
     445             : 
     446             :       USE m_intgr, ONLY: intgr3, intgr0, intgr1
     447             :       USE m_constants, ONLY: c_light
     448             :       USE m_setcor
     449             :       USE m_differ
     450             :       USE m_types
     451             :       IMPLICIT NONE
     452             : 
     453             :       TYPE(t_dimension), INTENT(IN)   :: dimension
     454             :       TYPE(t_input), INTENT(IN)       :: input
     455             :       TYPE(t_atoms), INTENT(IN)       :: atoms
     456             :       INTEGER, INTENT(OUT)            :: maxindxc, lmaxcd
     457             : 
     458             :       !  - local scalars -
     459             :       INTEGER              :: i, j, itype, korb, ncmsh, nst, ierr
     460             :       REAL                 :: e, fj, fl, fn, t, bmu, c
     461             :       REAL                 :: d, dxx, rn, rnot, z, t1, rr
     462             : 
     463             :       !  - local arrays -
     464           0 :       INTEGER              :: kappa(dimension%nstd), nprnc(dimension%nstd)
     465           0 :       INTEGER              :: nindxcr(0:dimension%nstd, atoms%ntype)
     466           0 :       REAL                 :: occ(dimension%nstd), occ_h(dimension%nstd, 2), a(dimension%msh), b(dimension%msh)
     467           0 :       INTEGER              :: lmaxc(atoms%ntype)
     468             : 
     469             :       !   - intrinsic functions -
     470             :       INTRINSIC exp, iabs, isign
     471             : 
     472           0 :       c = c_light(1.0)
     473             : 
     474             :       ! this loop determines the dimensions
     475             : 
     476           0 :       lmaxc = 0; nindxcr = 0
     477           0 :       DO itype = 1, atoms%ntype
     478           0 :          z = atoms%zatom(itype)
     479           0 :          dxx = atoms%dx(itype)
     480           0 :          bmu = 0.0
     481           0 :          CALL setcor(itype, input%jspins, atoms, input, bmu, nst, kappa, nprnc, occ_h)
     482             : 
     483           0 :          occ(1:nst) = occ_h(1:nst, 1)
     484             : 
     485           0 :          rnot = atoms%rmsh(1, itype)
     486           0 :          d = exp(atoms%dx(itype))
     487           0 :          ncmsh = nint(log((atoms%rmt(itype) + 10.0)/rnot)/dxx + 1)
     488           0 :          ncmsh = min(ncmsh, dimension%msh)
     489           0 :          rn = rnot*(d**(ncmsh - 1))
     490             : 
     491           0 :          nst = atoms%ncst(itype)
     492             : 
     493           0 :          DO korb = 1, nst
     494           0 :             IF (occ(korb) == 0) CYCLE
     495           0 :             fn = nprnc(korb)
     496           0 :             fj = iabs(kappa(korb)) - .5e0
     497             : 
     498           0 :             fl = fj + (.5e0)*isign(1, kappa(korb))
     499           0 :             e = -2*(z/(fn + fl))**2
     500             : 
     501           0 :             nindxcr(NINT(fl), itype) = nindxcr(NINT(fl), itype) + 1
     502           0 :             lmaxc(itype) = max(lmaxc(itype), NINT(fl))
     503             :          END DO
     504             : 
     505             :       END DO
     506             : 
     507           0 :       lmaxcd = maxval(lmaxc)
     508           0 :       maxindxc = maxval(nindxcr(0, :), nint((maxval(nindxcr)/2.0)))
     509             : 
     510           0 :    END SUBROUTINE core_init
     511             : 
     512             : END MODULE m_read_core

Generated by: LCOV version 1.13