vmts.F90 7.44 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

84
    
85 86 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
    ! 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
127
    !$ vtl = vtl + vtl_loc
128 129
    !$omp end critical
    !$omp end parallel
130
#ifdef CPP_MPI
131 132 133 134
    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/) )
135
    deallocate( c_b )
136 137
#endif

138 139 140


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

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

148
    nat = 1
149 150 151 152
    do n = 1, atoms%ntype
      nd = atoms%ntypsy(nat)
      imax = atoms%jri(n)
      lmax = sphhar%llh(sphhar%nlh(nd), nd)
153
      if ( potdenType == POTDEN_TYPE_POTYUK ) then
154 155
        !do concurrent (i = 1:imax)
        do i = 1,imax
156 157 158 159 160
          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)
161
        if ( potdenType == POTDEN_TYPE_POTYUK ) then
162 163 164 165 166 167 168 169 170 171 172 173 174
          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)
175
        vr(1:imax,lh,n) = green_factor * (   green_1(1:imax) * ( termsR - integral_2(1:imax) ) &
176 177 178 179
                                           + green_2(1:imax) *            integral_1(1:imax)   )
      end do
      nat = nat + atoms%neq(n)
    end do
180
    if ( potdenType == POTDEN_TYPE_POTYUK ) then
181 182 183
      deallocate( il, kl )
    end if

184 185 186
    

    if ( potdenType /= POTDEN_TYPE_POTYUK ) then
187 188 189 190
      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
191 192
    end if

193 194 195
  end subroutine vmts

end module m_vmts