vmts.F90 7.5 KB
Newer Older
1
module m_vmts
2

3
contains
4

5
  subroutine vmts( input, mpi, stars, sphhar, atoms, sym, cell, oneD, vpw, rho, potdenType, vr )
6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26

  !-------------------------------------------------------------------------
  ! This subroutine calculates the lattice harmonics expansion coefficients 
  ! of the Coulomb / Yukawa potential for all atom types.                     
  !
  ! In more detail:
  ! the radius-dependent coefficient of the potential is
  ! V_{lm}(r) = \int_0^R G_l(r,r') \rho_{lm}(r') r'^2 dr' + V_{lm}(R) P_l(r)
  ! where 
  ! [sphere boundary contribution - first part of code:]
  ! V_{lm}(R) is the spherical harmonics expansion of the interstitial 
  ! potential,
  ! P_l(r) is the derivative of the spherical Green's function on the sphere 
  ! boundary (times a constant),
  ! [sphere interior contribution - second part of code:]
  ! G_l(r,r') ) = green_factor u_1(r_<) u_2(r_>) is the spherical Green's 
  ! function with homogeneous solutions u_1 and u_2 and 
  ! r_< = min(r,r') and r_> = max(r,r')
  ! The integral is split in a part where r'=r_< and a part where r'=r_> and
  ! the integral from r to R is split in \int_0^R - \int_0^r. Resulting 
  ! terms depending solely on R (not on r) are summarised in the variable 
27 28 29
  ! termsR.
  !
  ! More information and equations can be found in
30
  ! F. Tran, P. Blaha: Phys. Rev. B 83, 235118(2011) 
31
  !-------------------------------------------------------------------------
32 33

#include"cpp_double.h"
34 35 36 37 38 39
    use m_constants
    use m_types
    use m_intgr, only : intgr2
    use m_phasy1
    use m_sphbes
    use m_od_phasy
40
    use m_SphBessel
41 42 43 44 45 46 47 48 49 50 51 52
    implicit none

    type(t_input),  intent(in)        :: input
    type(t_mpi),    intent(in)        :: mpi
    type(t_stars),  intent(in)        :: stars
    type(t_sphhar), intent(in)        :: sphhar
    type(t_atoms),  intent(in)        :: atoms
    type(t_sym),    intent(in)        :: sym
    type(t_cell),   intent(in)        :: cell
    type(t_oneD),   intent(in)        :: oneD
    complex,        intent(in)        :: vpw(:)!(stars%ng3,input%jspins)
    real,           intent(in)        :: rho(:,0:,:)!(atoms%jmtd,0:sphhar%nlhd,atoms%ntype)
53
    integer,        intent(in)        :: potdenType
54 55 56 57 58 59 60 61 62 63 64 65 66
    real,           intent(out)       :: vr(:,0:,:)!(atoms%jmtd,0:sphhar%nlhd,atoms%ntype)

    complex                           :: cp, sm
    integer                           :: i, jm, k, l, lh, n, nd, lm, n1, nat, m, imax, lmax
    complex                           :: vtl(0:sphhar%nlhd, atoms%ntype)
    complex                           :: pylm(( atoms%lmaxd + 1 ) ** 2, atoms%ntype)
    real                              :: green_factor, termsR
    real                              :: green_1    (1:atoms%jmtd), green_2    (1:atoms%jmtd)
    real                              :: integrand_1(1:atoms%jmtd), integrand_2(1:atoms%jmtd)
    real                              :: integral_1 (1:atoms%jmtd), integral_2 (1:atoms%jmtd)
    real                              :: sbf(0:atoms%lmaxd)
    real, allocatable, dimension(:,:) :: il, kl
    
67
    !$ COMPLEX vtl_loc(0:sphhar%nlhd,atoms%ntype)
68
#ifdef CPP_MPI
69 70 71
    include 'mpif.h'
    integer                       :: ierr(3)
    complex, allocatable          :: c_b(:)
72

73
    external MPI_REDUCE
74
#endif
75
    integer :: OMP_GET_NUM_THREADS, OMP_GET_THREAD_NUM
76 77


78 79

    ! SPHERE BOUNDARY CONTRIBUTION to the coefficients calculated from the values
80 81 82
    ! of the interstitial Coulomb / Yukawa potential on the sphere boundary

    vtl(:,:) = cmplx( 0.0, 0.0 )
83
#ifdef CPP_MPI
84
    call MPI_BCAST( vpw, size(vpw), CPP_MPI_COMPLEX, 0, mpi, ierr )
85 86
#endif

87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128
    ! q=0 component
    if ( mpi%irank == 0 ) then
      vtl(0,1:atoms%ntype) = sfp_const * vpw(1)
    end if

    ! q/=0 components
    !$omp parallel default( none ) &
    !$omp& shared( mpi, stars, vpw, oneD, atoms, sym, cell, sphhar, vtl ) &
    !$omp& private( k, cp, pylm, nat, n, sbf, nd, lh, sm, jm, m, lm, l ) &
    !$omp& private( vtl_loc ) 
    !$ vtl_loc(:,:) = cmplx(0.d0,0.d0)
    !$omp do
    do k = mpi%irank+2, stars%ng3, mpi%isize
      cp = vpw(k) * stars%nstr(k)
      if ( .not. oneD%odi%d1 ) then
        call phasy1( atoms, stars, sym, cell, k, pylm )
      else
        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 )
      end if
      nat = 1
      do n = 1, atoms%ntype
        call sphbes( atoms%lmax(n), stars%sk3(k) * atoms%rmt(n), sbf )
        nd = atoms%ntypsy(nat)
        do lh = 0, sphhar%nlh(nd)
          l = sphhar%llh(lh,nd)
          sm = (0.0,0.0)
          do jm = 1, sphhar%nmem(lh,nd)
            m = sphhar%mlh(jm,lh,nd)
            lm = l * ( l + 1 ) + m + 1 
            sm = sm + conjg( sphhar%clnu(jm,lh,nd) ) * pylm(lm,n)
          end do
        !$ if (.false.) then
        vtl(lh,n) = vtl(lh,n) + cp * sbf(l) * sm
        !$ end if
        !$ vtl_loc(lh,n) = vtl_loc(lh,n) + cp * sbf(l) * sm
        end do
        nat = nat + atoms%neq(n)
      end do
    end do
    !$omp end do
    !$omp critical
129
    !$ vtl = vtl + vtl_loc
130 131
    !$omp end critical
    !$omp end parallel
132
#ifdef CPP_MPI
133 134 135 136
    n1 = ( sphhar%nlhd + 1 ) * atoms%ntype
    allocate( c_b(n1) )
    call MPI_REDUCE( vtl, c_b, n1, CPP_MPI_COMPLEX, MPI_SUM, 0, mpi%mpi_comm, ierr )
    if ( mpi%irank == 0 ) vtl = reshape( c_b, (/sphhar%nlhd+1,atoms%ntype/) )
137
    deallocate( c_b )
138 139
#endif

140 141 142


    ! SPHERE INTERIOR CONTRIBUTION to the coefficients calculated from the 
143
    ! values of the sphere Coulomb/Yukawa potential on the sphere boundary
144

145
    if( mpi%irank == 0 ) then
146
    if ( potdenType == POTDEN_TYPE_POTYUK ) then
147 148 149
      allocate( il(0:atoms%lmaxd, 1:atoms%jmtd), kl(0:atoms%lmaxd, 1:atoms%jmtd) )
    end if

150
    nat = 1
151 152 153 154
    do n = 1, atoms%ntype
      nd = atoms%ntypsy(nat)
      imax = atoms%jri(n)
      lmax = sphhar%llh(sphhar%nlh(nd), nd)
155
      if ( potdenType == POTDEN_TYPE_POTYUK ) then
156 157 158 159 160 161
        do concurrent (i = 1:imax)
          call ModSphBessel( il(0:,i), kl(0:,i), input%preconditioning_param * atoms%rmsh(i,n), lmax )
        end do
      end if
      do lh = 0, sphhar%nlh(nd)
        l = sphhar%llh(lh,nd)
162
        if ( potdenType == POTDEN_TYPE_POTYUK ) then
163 164 165 166 167 168 169 170 171 172 173 174 175
          green_1(1:imax) = il(l,1:imax)
          green_2(1:imax) = kl(l,1:imax)
          green_factor    = fpi_const * input%preconditioning_param
        else
          green_1(1:imax) = atoms%rmsh(1:imax,n) ** l
          green_2(1:imax) = 1.0 / ( green_1(1:imax) * atoms%rmsh(1:imax,n) )
          green_factor    = fpi_const / ( 2 * l + 1 )
        end if
        integrand_1(1:imax) = green_1(1:imax) * rho(1:imax,lh,n)
        integrand_2(1:imax) = green_2(1:imax) * rho(1:imax,lh,n)
        call intgr2( integrand_1(1:imax), atoms%rmsh(1,n), atoms%dx(n), imax, integral_1(1:imax) )
        call intgr2( integrand_2(1:imax), atoms%rmsh(1,n), atoms%dx(n), imax, integral_2(1:imax) )
        termsR = integral_2(imax) + ( vtl(lh,n) / green_factor - integral_1(imax) * green_2(imax) ) / green_1(imax)
176
        vr(1:imax,lh,n) = green_factor * (   green_1(1:imax) * ( termsR - integral_2(1:imax) ) &
177 178 179 180
                                           + green_2(1:imax) *            integral_1(1:imax)   )
      end do
      nat = nat + atoms%neq(n)
    end do
181
    if ( potdenType == POTDEN_TYPE_POTYUK ) then
182 183 184
      deallocate( il, kl )
    end if

185 186 187
    

    if ( potdenType /= POTDEN_TYPE_POTYUK ) then
188 189 190 191 192
      do n = 1, atoms%ntype
        vr(1:atoms%jri(n),0,n) = vr(1:atoms%jri(n),0,n) - sfp_const * ( 1.0 / atoms%rmsh(1:atoms%jri(n),n) - 1.0 / atoms%rmt(n) ) * atoms%zatom(n)
      end do
    end if

193 194
    end if

195 196 197
  end subroutine vmts

end module m_vmts