LCOV - code coverage report
Current view: top level - vgen - mpmom.F90 (source / functions) Hit Total Coverage
Test: combined.info Lines: 88 89 98.9 %
Date: 2019-09-08 04:53:50 Functions: 3 3 100.0 %

          Line data    Source code
       1             : module m_mpmom
       2             :   !     ***********************************************************
       3             :   !     calculation of the multipole moments of the original charge 
       4             :   !     density minus the interstitial charge density 
       5             :   !     for each atom type 
       6             :   !     
       7             :   !     For yukawa_residual = .true. the subroutines calculate the 
       8             :   !     multipole moments for the Yukawa potential instead of the
       9             :   !     Coulomb potential. This is used in the preconditioning of
      10             :   !     the SCF iteration for metallic systems.
      11             :   !
      12             :   !     qlmo(m,l,n) : mult.mom. of the mufftn-tin charge density
      13             :   !     qlmp(m,l,n) : mult.mom. of the plane-wave charge density
      14             :   !     qlm (m,l,n) : (output) difference of the former quantities
      15             :   !     
      16             :   !     references:
      17             :   !     for both the Coulomb and the Yukawa pseudo charge density:
      18             :   !     F. Tran, P. Blaha: Phys. Rev. B 83, 235118 (2011)
      19             :   !     or see the original paper for the normal Coulomb case only:
      20             :   !     M. Weinert: J.Math.Phys. 22(11) (1981) p.2434 eq. (10)-(15) 
      21             :   !     ***********************************************************
      22             : 
      23             : contains
      24             : 
      25         346 :   subroutine mpmom( input, mpi, atoms, sphhar, stars, sym, cell, oneD, qpw, rho, potdenType, qlm )
      26             : 
      27             :     use m_types
      28             :     implicit none
      29             : 
      30             :     type(t_input),   intent(in)   :: input
      31             :     type(t_mpi),     intent(in)   :: mpi
      32             :     type(t_oneD),    intent(in)   :: oneD
      33             :     type(t_sym),     intent(in)   :: sym
      34             :     type(t_stars),   intent(in)   :: stars
      35             :     type(t_cell),    intent(in)   :: cell
      36             :     type(t_sphhar),  intent(in)   :: sphhar
      37             :     type(t_atoms),   intent(in)   :: atoms
      38             :     real,            intent(in)   :: rho(:,0:,:) !(atoms%jmtd,0:sphhar%nlhd,atoms%ntype)
      39             :     complex,         intent(in)   :: qpw(:)      !(stars%ng3)
      40             :     integer,         intent(in)   :: potdenType
      41             :     complex,         intent(out)  :: qlm(-atoms%lmaxd:atoms%lmaxd,0:atoms%lmaxd,atoms%ntype)
      42             : 
      43             :     integer                       :: j, jm, lh, mb, mem, mems, n, nd, l, nat, m
      44         692 :     complex                       :: qlmo(-atoms%lmaxd:atoms%lmaxd,0:atoms%lmaxd,atoms%ntype)
      45         692 :     complex                       :: qlmp(-atoms%lmaxd:atoms%lmaxd,0:atoms%lmaxd,atoms%ntype)
      46             : 
      47             :     ! multipole moments of original charge density
      48         346 :     if ( mpi%irank == 0 ) then
      49         173 :       call mt_moments( input, atoms, sphhar, rho(:,:,:), potdenType, qlmo )
      50             :     end if
      51             : 
      52             :     ! multipole moments of the interstitial charge density in the spheres
      53         346 :     call pw_moments( input, mpi, stars, atoms, cell, sym, oneD, qpw(:), potdenType, qlmp )
      54             : 
      55         346 :     if ( mpi%irank == 0 ) then
      56             :       ! see (A14)
      57         517 :       qlm = qlmo - qlmp
      58             :       ! output section
      59         173 :       nat = 1
      60         517 :       do n = 1, atoms%ntype
      61         344 :         write( 6, fmt=8000 ) n
      62         344 :         nd = atoms%ntypsy(nat)
      63        5867 :         do lh = 0, sphhar%nlh(nd)
      64        5523 :           l = sphhar%llh(lh,nd)
      65        5523 :           mems = sphhar%nmem(lh,nd)
      66       17322 :           do mem = 1, mems
      67       11455 :             m = sphhar%mlh(mem,lh,nd)
      68       16978 :             write( 6, fmt=8010 ) l, m, qlmo(m,l,n), qlmp(m,l,n)
      69             :           end do
      70             :         end do
      71         517 :         nat = nat + atoms%neq(n)
      72             :       end do
      73             : 
      74             : 8000   FORMAT (/,10x,'multipole moments for atom type=',i5,/,/,t3,'l',t7,&
      75             :             &       'm',t27,'original',t57,'plane wave')
      76             : 8010   FORMAT (1x,i2,2x,i3,2x,2 (5x,2e15.5))
      77             :        !
      78             :     end if ! mpi%irank == 0
      79             : 
      80         346 :   end subroutine mpmom
      81             : 
      82             : 
      83         173 :   subroutine mt_moments( input, atoms, sphhar, rho, potdenType, qlmo )
      84             :     ! multipole moments of original charge density
      85             :     ! see (A15) (Coulomb case) or (A17) (Yukawa case)
      86             : 
      87             :     use m_intgr,     only: intgr3
      88             :     use m_constants, only: sfp_const, POTDEN_TYPE_POTYUK
      89             :     use m_types
      90             :     use m_DoubleFactorial
      91             :     use m_SphBessel
      92             :     implicit none
      93             : 
      94             :     type(t_input),  intent(in)        :: input
      95             :     type(t_sphhar), intent(in)        :: sphhar
      96             :     type(t_atoms),  intent(in)        :: atoms
      97             :     real,           intent(in)        :: rho(: ,0:, :)
      98             :     integer,        intent(in)        :: potdenType
      99             :     complex,        intent(out)       :: qlmo(-atoms%lmaxd:,0:,:)
     100             : 
     101             :     integer                           :: n, ns, jm, nl, l, j, mb, m, nat, i, imax, lmax
     102             :     real                              :: fint
     103         346 :     real                              :: f( maxval( atoms%jri ) )
     104         173 :     real, allocatable, dimension(:,:) :: il, kl
     105             : 
     106         173 :     if ( potdenType == POTDEN_TYPE_POTYUK ) then
     107           3 :       allocate( il(0:atoms%lmaxd, 1:atoms%jmtd), kl(0:atoms%lmaxd, 1:atoms%jmtd) )
     108             :     end if
     109             : 
     110         517 :     qlmo = 0.0
     111         173 :     nat = 1
     112         517 :     do n = 1, atoms%ntype
     113         344 :       ns = atoms%ntypsy(nat)
     114         344 :       jm = atoms%jri(n)
     115         344 :       imax = atoms%jri(n)
     116         344 :       lmax = sphhar%llh(sphhar%nlh(ns), ns)
     117         344 :       if ( potdenType == POTDEN_TYPE_POTYUK ) then
     118             :         !do concurrent (i = 1:imax)
     119       22725 :         do i = 1,imax
     120       11370 :           call ModSphBessel( il(:, i), kl(:, i), input%preconditioning_param * atoms%rmsh(i, n), lmax )
     121             :         end do
     122             :       end if
     123        5867 :       do nl = 0, sphhar%nlh(ns)
     124        5523 :         l = sphhar%llh(nl,ns)
     125     3432098 :         do j = 1, jm
     126     3432098 :           if ( potdenType /= POTDEN_TYPE_POTYUK ) then
     127     3269876 :             f(j) = atoms%rmsh(j,n) ** l * rho(j,nl,n)
     128             :           else
     129      156699 :             f(j) = il(l, j) * rho(j,nl,n)
     130             :           end if
     131             :         end do
     132        5523 :         call intgr3( f, atoms%rmsh(:,n), atoms%dx(n), jm, fint )
     133        5523 :         if ( potdenType == POTDEN_TYPE_POTYUK ) then
     134         207 :           fint = fint * DoubleFactorial( l ) / input%preconditioning_param ** l
     135             :         end if
     136       17322 :         do mb = 1, sphhar%nmem(nl,ns)
     137       11455 :           m = sphhar%mlh(mb,nl,ns)
     138       16978 :           qlmo(m,l,n) = qlmo(m,l,n) + sphhar%clnu(mb,nl,ns) * fint
     139             :         end do
     140             :       end do
     141         344 :       if ( potdenType /= POTDEN_TYPE_POTYUK ) then
     142         329 :         qlmo(0,0,n) = qlmo(0,0,n) - atoms%zatom(n) / sfp_const
     143             :       end if
     144         517 :       nat = nat + atoms%neq(n)
     145             :     end do
     146             : 
     147         346 :   end subroutine mt_moments
     148             : 
     149             : 
     150         346 :   subroutine pw_moments( input, mpi, stars, atoms, cell, sym, oneD, qpw, potdenType, qlmp_out )
     151             :     ! multipole moments of the interstitial charge in the spheres
     152             : 
     153             :     use m_phasy1
     154             :     use m_sphbes
     155             :     use m_od_phasy
     156             :     use m_constants, only: sfp_const, POTDEN_TYPE_POTYUK
     157             :     use m_types
     158             :     use m_DoubleFactorial
     159             :     use m_SphBessel
     160             :     implicit none
     161             : 
     162             :     type(t_input),    intent(in)   :: input
     163             :     type(t_mpi),      intent(in)   :: mpi
     164             :     type(t_oneD),     intent(in)   :: oneD
     165             :     type(t_sym),      intent(in)   :: sym
     166             :     type(t_stars),    intent(in)   :: stars
     167             :     type(t_cell),     intent(in)   :: cell
     168             :     type(t_atoms),    intent(in)   :: atoms
     169             :     complex,          intent(in)   :: qpw(:)
     170             :     integer,          intent(in)   :: potdenType
     171             :     complex,          intent(out)  :: qlmp_out(-atoms%lmaxd:,0:,:)
     172             : 
     173             :     integer                        :: n, k, l, ll1, lm, ierr(3), m
     174             :     complex                        :: sk3i, cil, nqpw
     175        1380 :     complex                        :: pylm(( maxval( atoms%lmax ) + 1 ) ** 2, atoms%ntype)
     176             :     real                           :: sk3r, rl2
     177         692 :     real                           :: aj(0:maxval( atoms%lmax ) + 1 )
     178             :     logical                        :: od
     179        1380 :     real                           :: il(0:maxval( atoms%lmax ) + 1 )
     180        1380 :     real                           :: kl(0:maxval( atoms%lmax ) + 1 )
     181             : #ifdef CPP_MPI
     182             :     include 'mpif.h'
     183             : #endif
     184         692 :     complex                        :: qlmp(-atoms%lmaxd:atoms%lmaxd,0:atoms%lmaxd,atoms%ntype)
     185             : 
     186        1034 :     qlmp = 0.0
     187         346 :     if ( mpi%irank == 0 ) then
     188             :       ! q=0 term: see (A19) (Coulomb case) or (A20) (Yukawa case)
     189         861 :       do n = 1, atoms%ntype
     190         517 :         if ( potdenType /= POTDEN_TYPE_POTYUK ) then  
     191         329 :           qlmp(0,0,n) = qpw(1) * stars%nstr(1) * atoms%volmts(n) / sfp_const
     192             :         else
     193          15 :           call ModSphBessel( il(0:1), kl(0:1), input%preconditioning_param * atoms%rmt(n), 1 )
     194          15 :           qlmp(0,0,n) = qpw(1) * stars%nstr(1) * sfp_const * atoms%rmt(n) ** 2 * il(1) / input%preconditioning_param
     195             :         end if
     196             :       end do
     197             :     end if
     198             : 
     199             :     ! q/=0 terms: see (A16) (Coulomb case) or (A18) (Yukawa case)
     200         346 :     od = oneD%odi%d1
     201             : !    !$omp parallel do default( shared ) private( pylm, nqpw, n, sk3r, aj, rl2, sk3i, &
     202             : !    !$omp& l, cil, ll1, m, lm, k ) reduction( +:qlmp )
     203      174821 :     do k = mpi%irank+2, stars%ng3, mpi%isize
     204      174129 :       if ( od ) then
     205             :         call od_phasy( atoms%ntype, stars%ng3, atoms%nat, atoms%lmaxd, atoms%ntype, &
     206           0 :              atoms%neq, atoms%lmax, atoms%taual, cell%bmat, stars%kv3, k, oneD%odi, oneD%ods, pylm)
     207             :       else
     208      174129 :         call phasy1( atoms, stars, sym, cell, k, pylm )
     209             :       end if
     210             :      
     211      174129 :       nqpw = qpw(k) * stars%nstr(k)
     212      464256 :       do n = 1, atoms%ntype
     213      289781 :         sk3r = stars%sk3(k) * atoms%rmt(n)
     214      289781 :         call sphbes( atoms%lmax(n) + 1, sk3r, aj )
     215      289781 :         rl2 = atoms%rmt(n) ** 2
     216      289781 :         if ( potdenType == POTDEN_TYPE_POTYUK ) then
     217       19170 :           call ModSphBessel( il(0:atoms%lmax(n)+1), kl(0:atoms%lmax(n)+1), input%preconditioning_param * atoms%rmt(n), atoms%lmax(n) + 1 )
     218       19170 :           sk3i = nqpw / ( stars%sk3(k) ** 2 + input%preconditioning_param ** 2 ) * rl2
     219             :         else
     220      270611 :           sk3i = nqpw / stars%sk3(k)
     221             :         end if
     222     3143321 :         do l = 0, atoms%lmax(n)
     223     2679411 :           if ( potdenType == POTDEN_TYPE_POTYUK ) then
     224      172530 :             cil = ( stars%sk3(k) * il(l) * aj(l+1) + input%preconditioning_param * il(l+1) * aj(l) ) * ( DoubleFactorial( l ) / input%preconditioning_param ** l ) * sk3i
     225             :           else
     226     2506881 :             cil = aj(l+1) * sk3i * rl2  
     227     2506881 :             rl2 = rl2 * atoms%rmt(n)
     228             :           end if
     229     2679411 :           ll1 = l * ( l + 1 ) + 1
     230    28150373 :           do m = -l, l
     231    25181181 :             lm = ll1 + m
     232    27860592 :             qlmp(m,l,n) = qlmp(m,l,n) + cil * pylm(lm,n)
     233             :           end do
     234             :         end do                  ! l = 0, atoms%lmax(n)
     235             :       end do                    ! n = 1, atoms%ntype
     236             :     end do                      ! k = 2, stars%ng3
     237             : !    !$omp end parallel do
     238             : #ifdef CPP_MPI
     239         346 :     CALL MPI_REDUCE( qlmp, qlmp_out, SIZE(qlmp), MPI_DOUBLE_COMPLEX, MPI_SUM,0, mpi%mpi_comm, ierr )
     240             : #else
     241             :     qlmp_out = qlmp
     242             : #endif
     243             : 
     244         346 :   end subroutine pw_moments
     245             : 
     246             : end module m_mpmom

Generated by: LCOV version 1.13