mpmom.F90 9.31 KB
Newer Older
1
module m_mpmom
2
  !     ***********************************************************
3 4 5 6 7 8 9 10
  !     calculation of the multipole moments of the original charge 
  !     density minus the interstitial charge density 
  !     for each atom type 
  !     
  !     For yukawa_residual = .true. the subroutines calculate the 
  !     multipole moments for the Yukawa potential instead of the
  !     Coulomb potential. This is used in the preconditioning of
  !     the SCF iteration for metallic systems.
11 12 13 14 15
  !
  !     qlmo(m,l,n) : mult.mom. of the mufftn-tin charge density
  !     qlmp(m,l,n) : mult.mom. of the plane-wave charge density
  !     qlm (m,l,n) : (output) difference of the former quantities
  !     
16 17 18 19 20
  !     references:
  !     for both the Coulomb and the Yukawa pseudo charge density:
  !     F. Tran, P. Blaha: Phys. Rev. B 83, 235118 (2011)
  !     or see the original paper for the normal Coulomb case only:
  !     M. Weinert: J.Math.Phys. 22(11) (1981) p.2434 eq. (10)-(15) 
21
  !     ***********************************************************
22

23 24
contains

25
  subroutine mpmom( input, mpi, atoms, sphhar, stars, sym, cell, oneD, qpw, rho, potdenType, qlm )
26 27 28 29 30 31 32 33 34 35 36 37 38 39

    use m_types
    implicit none

    type(t_input),   intent(in)   :: input
    type(t_mpi),     intent(in)   :: mpi
    type(t_oneD),    intent(in)   :: oneD
    type(t_sym),     intent(in)   :: sym
    type(t_stars),   intent(in)   :: stars
    type(t_cell),    intent(in)   :: cell
    type(t_sphhar),  intent(in)   :: sphhar
    type(t_atoms),   intent(in)   :: atoms
    real,            intent(in)   :: rho(:,0:,:) !(atoms%jmtd,0:sphhar%nlhd,atoms%ntype)
    complex,         intent(in)   :: qpw(:)      !(stars%ng3)
40
    integer,         intent(in)   :: potdenType
41 42 43 44 45 46 47 48
    complex,         intent(out)  :: qlm(-atoms%lmaxd:atoms%lmaxd,0:atoms%lmaxd,atoms%ntype)

    integer                       :: j, jm, lh, mb, mem, mems, n, nd, l, nat, m
    complex                       :: qlmo(-atoms%lmaxd:atoms%lmaxd,0:atoms%lmaxd,atoms%ntype)
    complex                       :: qlmp(-atoms%lmaxd:atoms%lmaxd,0:atoms%lmaxd,atoms%ntype)

    ! multipole moments of original charge density
    if ( mpi%irank == 0 ) then
49
      call mt_moments( input, atoms, sphhar, rho(:,:,:), potdenType, qlmo )
50 51 52
    end if

    ! multipole moments of the interstitial charge density in the spheres
53
    call pw_moments( input, mpi, stars, atoms, cell, sym, oneD, qpw(:), potdenType, qlmp )
54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73

    if ( mpi%irank == 0 ) then
      ! see (A14)
      qlm = qlmo - qlmp
      ! output section
      nat = 1
      do n = 1, atoms%ntype
        write( 6, fmt=8000 ) n
        nd = atoms%ntypsy(nat)
        do lh = 0, sphhar%nlh(nd)
          l = sphhar%llh(lh,nd)
          mems = sphhar%nmem(lh,nd)
          do mem = 1, mems
            m = sphhar%mlh(mem,lh,nd)
            write( 6, fmt=8010 ) l, m, qlmo(m,l,n), qlmp(m,l,n)
          end do
        end do
        nat = nat + atoms%neq(n)
      end do

74 75 76 77
8000   FORMAT (/,10x,'multipole moments for atom type=',i5,/,/,t3,'l',t7,&
            &       'm',t27,'original',t57,'plane wave')
8010   FORMAT (1x,i2,2x,i3,2x,2 (5x,2e15.5))
       !
78
    end if ! mpi%irank == 0
79

80
  end subroutine mpmom
81 82


83
  subroutine mt_moments( input, atoms, sphhar, rho, potdenType, qlmo )
84 85
    ! multipole moments of original charge density
    ! 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 91
    use m_DoubleFactorial
    use m_SphBessel
92
    implicit none
93

94 95 96 97
    type(t_input),  intent(in)        :: input
    type(t_sphhar), intent(in)        :: sphhar
    type(t_atoms),  intent(in)        :: atoms
    real,           intent(in)        :: rho(: ,0:, :)
98
    integer,        intent(in)        :: potdenType
99
    complex,        intent(out)       :: qlmo(-atoms%lmaxd:,0:,:)
100

101 102 103 104
    integer                           :: n, ns, jm, nl, l, j, mb, m, nat, i, imax, lmax
    real                              :: fint
    real                              :: f( maxval( atoms%jri ) )
    real, allocatable, dimension(:,:) :: il, kl
105

106
    if ( potdenType == POTDEN_TYPE_POTYUK ) then
107 108 109 110
      allocate( il(0:atoms%lmaxd, 1:atoms%jmtd), kl(0:atoms%lmaxd, 1:atoms%jmtd) )
    end if

    qlmo = 0.0
111
    nat = 1
112 113 114 115 116
    do n = 1, atoms%ntype
      ns = atoms%ntypsy(nat)
      jm = atoms%jri(n)
      imax = atoms%jri(n)
      lmax = sphhar%llh(sphhar%nlh(ns), ns)
117
      if ( potdenType == POTDEN_TYPE_POTYUK ) then
118 119
        !do concurrent (i = 1:imax)
        do i = 1,imax
120 121 122 123 124 125
          call ModSphBessel( il(:, i), kl(:, i), input%preconditioning_param * atoms%rmsh(i, n), lmax )
        end do
      end if
      do nl = 0, sphhar%nlh(ns)
        l = sphhar%llh(nl,ns)
        do j = 1, jm
126
          if ( potdenType /= POTDEN_TYPE_POTYUK ) then
127 128 129 130 131 132
            f(j) = atoms%rmsh(j,n) ** l * rho(j,nl,n)
          else
            f(j) = il(l, j) * rho(j,nl,n)
          end if
        end do
        call intgr3( f, atoms%rmsh(:,n), atoms%dx(n), jm, fint )
133
        if ( potdenType == POTDEN_TYPE_POTYUK ) then
134 135 136 137 138 139 140
          fint = fint * DoubleFactorial( l ) / input%preconditioning_param ** l
        end if
        do mb = 1, sphhar%nmem(nl,ns)
          m = sphhar%mlh(mb,nl,ns)
          qlmo(m,l,n) = qlmo(m,l,n) + sphhar%clnu(mb,nl,ns) * fint
        end do
      end do
141
      if ( potdenType /= POTDEN_TYPE_POTYUK ) then
142 143 144 145 146 147 148 149
        qlmo(0,0,n) = qlmo(0,0,n) - atoms%zatom(n) / sfp_const
      end if
      nat = nat + atoms%neq(n)
    end do

  end subroutine mt_moments


150
  subroutine pw_moments( input, mpi, stars, atoms, cell, sym, oneD, qpw, potdenType, qlmp_out )
151 152 153 154 155
    ! multipole moments of the interstitial charge in the spheres

    use m_phasy1
    use m_sphbes
    use m_od_phasy
156
    use m_constants, only: sfp_const, POTDEN_TYPE_POTYUK
157
    use m_types
158 159
    use m_DoubleFactorial
    use m_SphBessel
160 161 162 163 164 165 166 167 168
    implicit none

    type(t_input),    intent(in)   :: input
    type(t_mpi),      intent(in)   :: mpi
    type(t_oneD),     intent(in)   :: oneD
    type(t_sym),      intent(in)   :: sym
    type(t_stars),    intent(in)   :: stars
    type(t_cell),     intent(in)   :: cell
    type(t_atoms),    intent(in)   :: atoms
169
    complex,          intent(in)   :: qpw(:)
170
    integer,          intent(in)   :: potdenType
171 172 173 174 175 176 177 178 179 180
    complex,          intent(out)  :: qlmp_out(-atoms%lmaxd:,0:,:)

    integer                        :: n, k, l, ll1, lm, ierr(3), m
    complex                        :: sk3i, cil, nqpw
    complex                        :: pylm(( maxval( atoms%lmax ) + 1 ) ** 2, atoms%ntype)
    real                           :: sk3r, rl2
    real                           :: aj(0:maxval( atoms%lmax ) + 1 )
    logical                        :: od
    real                           :: il(0:maxval( atoms%lmax ) + 1 )
    real                           :: kl(0:maxval( atoms%lmax ) + 1 )
181
#ifdef CPP_MPI
182
    include 'mpif.h'
183
#endif
184 185 186 187 188 189
    complex                        :: qlmp(-atoms%lmaxd:atoms%lmaxd,0:atoms%lmaxd,atoms%ntype)

    qlmp = 0.0
    if ( mpi%irank == 0 ) then
      ! q=0 term: see (A19) (Coulomb case) or (A20) (Yukawa case)
      do n = 1, atoms%ntype
190
        if ( potdenType /= POTDEN_TYPE_POTYUK ) then  
191 192 193 194 195 196 197 198 199 200
          qlmp(0,0,n) = qpw(1) * stars%nstr(1) * atoms%volmts(n) / sfp_const
        else
          call ModSphBessel( il(0:1), kl(0:1), input%preconditioning_param * atoms%rmt(n), 1 )
          qlmp(0,0,n) = qpw(1) * stars%nstr(1) * sfp_const * atoms%rmt(n) ** 2 * il(1) / input%preconditioning_param
        end if
      end do
    end if

    ! q/=0 terms: see (A16) (Coulomb case) or (A18) (Yukawa case)
    od = oneD%odi%d1
201 202
!    !$omp parallel do default( shared ) private( pylm, nqpw, n, sk3r, aj, rl2, sk3i, &
!    !$omp& l, cil, ll1, m, lm, k ) reduction( +:qlmp )
203 204 205 206 207 208 209 210 211 212 213 214 215
    do k = mpi%irank+2, stars%ng3, mpi%isize
      if ( od ) then
        call od_phasy( atoms%ntype, stars%ng3, atoms%nat, atoms%lmaxd, atoms%ntype, &
             atoms%neq, atoms%lmax, atoms%taual, cell%bmat, stars%kv3, k, oneD%odi, oneD%ods, pylm)
      else
        call phasy1( atoms, stars, sym, cell, k, pylm )
      end if
     
      nqpw = qpw(k) * stars%nstr(k)
      do n = 1, atoms%ntype
        sk3r = stars%sk3(k) * atoms%rmt(n)
        call sphbes( atoms%lmax(n) + 1, sk3r, aj )
        rl2 = atoms%rmt(n) ** 2
216
        if ( potdenType == POTDEN_TYPE_POTYUK ) then
217 218 219 220 221 222
          call ModSphBessel( il(0:atoms%lmax(n)+1), kl(0:atoms%lmax(n)+1), input%preconditioning_param * atoms%rmt(n), atoms%lmax(n) + 1 )
          sk3i = nqpw / ( stars%sk3(k) ** 2 + input%preconditioning_param ** 2 ) * rl2
        else
          sk3i = nqpw / stars%sk3(k)
        end if
        do l = 0, atoms%lmax(n)
223
          if ( potdenType == POTDEN_TYPE_POTYUK ) then
224 225 226 227 228 229 230 231 232 233 234 235 236
            cil = ( stars%sk3(k) * il(l) * aj(l+1) + input%preconditioning_param * il(l+1) * aj(l) ) * ( DoubleFactorial( l ) / input%preconditioning_param ** l ) * sk3i
          else
            cil = aj(l+1) * sk3i * rl2  
            rl2 = rl2 * atoms%rmt(n)
          end if
          ll1 = l * ( l + 1 ) + 1
          do m = -l, l
            lm = ll1 + m
            qlmp(m,l,n) = qlmp(m,l,n) + cil * pylm(lm,n)
          end do
        end do                  ! l = 0, atoms%lmax(n)
      end do                    ! n = 1, atoms%ntype
    end do                      ! k = 2, stars%ng3
237
!    !$omp end parallel do
Daniel Wortmann's avatar
Daniel Wortmann committed
238
#ifdef CPP_MPI
239
    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 243
#endif

244
  end subroutine pw_moments
245

246
end module m_mpmom