vmts.F90 7.5 KB
 Miriam Hinzen committed Jun 06, 2018 1 module m_vmts  Miriam Hinzen committed Jun 04, 2018 2   Miriam Hinzen committed Jun 06, 2018 3 contains  Miriam Hinzen committed Jun 04, 2018 4   Miriam Hinzen committed Jun 07, 2018 5  subroutine vmts( input, mpi, stars, sphhar, atoms, sym, cell, oneD, vpw, rho, potdenType, vr )  Miriam Hinzen committed Jun 06, 2018 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  Miriam Hinzen committed Jun 06, 2018 27 28 29  ! termsR. ! ! More information and equations can be found in  Miriam Hinzen committed Jun 07, 2018 30  ! F. Tran, P. Blaha: Phys. Rev. B 83, 235118(2011)  Miriam Hinzen committed Jun 06, 2018 31  !-------------------------------------------------------------------------  Markus Betzinger committed Apr 26, 2016 32 33  #include"cpp_double.h"  Miriam Hinzen committed Jun 06, 2018 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  Miriam Hinzen committed Jun 04, 2018 40  use m_SphBessel  Miriam Hinzen committed Jun 06, 2018 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)  Miriam Hinzen committed Jun 07, 2018 53  integer, intent(in) :: potdenType  Miriam Hinzen committed Jun 06, 2018 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  Uliana Alekseeva committed Mar 02, 2018 67  !$COMPLEX vtl_loc(0:sphhar%nlhd,atoms%ntype)  Markus Betzinger committed Apr 26, 2016 68 #ifdef CPP_MPI  Miriam Hinzen committed Jun 06, 2018 69 70 71  include 'mpif.h' integer :: ierr(3) complex, allocatable :: c_b(:)  Gregor Michalicek committed May 29, 2018 72   Miriam Hinzen committed Jun 06, 2018 73  external MPI_REDUCE  Markus Betzinger committed Apr 26, 2016 74 #endif  Uliana Alekseeva committed Mar 02, 2018 75  integer :: OMP_GET_NUM_THREADS, OMP_GET_THREAD_NUM  Gregor Michalicek committed May 29, 2018 76 77   Miriam Hinzen committed Jun 07, 2018 78 79  ! SPHERE BOUNDARY CONTRIBUTION to the coefficients calculated from the values  Miriam Hinzen committed Jun 06, 2018 80 81 82  ! of the interstitial Coulomb / Yukawa potential on the sphere boundary vtl(:,:) = cmplx( 0.0, 0.0 )  Gregor Michalicek committed May 29, 2018 83 #ifdef CPP_MPI  Miriam Hinzen committed Jun 06, 2018 84  call MPI_BCAST( vpw, size(vpw), CPP_MPI_COMPLEX, 0, mpi, ierr )  Markus Betzinger committed Apr 26, 2016 85 86 #endif  Miriam Hinzen committed Jun 06, 2018 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  Uliana Alekseeva committed Mar 02, 2018 129  !$vtl = vtl + vtl_loc  Miriam Hinzen committed Jun 06, 2018 130 131  !$omp end critical !\$omp end parallel  Markus Betzinger committed Apr 26, 2016 132 #ifdef CPP_MPI  Miriam Hinzen committed Jun 06, 2018 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/) )  Miriam Hinzen committed Jun 07, 2018 137  deallocate( c_b )  Markus Betzinger committed Apr 26, 2016 138 139 #endif  Miriam Hinzen committed Jun 07, 2018 140 141 142  ! SPHERE INTERIOR CONTRIBUTION to the coefficients calculated from the  Miriam Hinzen committed Jun 06, 2018 143  ! values of the sphere Coulomb/Yukawa potential on the sphere boundary  Miriam Hinzen committed Jun 07, 2018 144   Miriam Hinzen committed Jun 11, 2018 145  if( mpi%irank == 0 ) then  Miriam Hinzen committed Jun 07, 2018 146  if ( potdenType == POTDEN_TYPE_POTYUK ) then  Miriam Hinzen committed Jun 06, 2018 147 148 149  allocate( il(0:atoms%lmaxd, 1:atoms%jmtd), kl(0:atoms%lmaxd, 1:atoms%jmtd) ) end if  Markus Betzinger committed Apr 26, 2016 150  nat = 1  Miriam Hinzen committed Jun 06, 2018 151 152 153 154  do n = 1, atoms%ntype nd = atoms%ntypsy(nat) imax = atoms%jri(n) lmax = sphhar%llh(sphhar%nlh(nd), nd)  Miriam Hinzen committed Jun 07, 2018 155  if ( potdenType == POTDEN_TYPE_POTYUK ) then  Miriam Hinzen committed Jun 06, 2018 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)  Miriam Hinzen committed Jun 07, 2018 162  if ( potdenType == POTDEN_TYPE_POTYUK ) then  Miriam Hinzen committed Jun 06, 2018 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)  Miriam Hinzen committed Jun 07, 2018 176  vr(1:imax,lh,n) = green_factor * ( green_1(1:imax) * ( termsR - integral_2(1:imax) ) &  Miriam Hinzen committed Jun 06, 2018 177 178 179 180  + green_2(1:imax) * integral_1(1:imax) ) end do nat = nat + atoms%neq(n) end do  Miriam Hinzen committed Jun 07, 2018 181  if ( potdenType == POTDEN_TYPE_POTYUK ) then  Miriam Hinzen committed Jun 06, 2018 182 183 184  deallocate( il, kl ) end if  Miriam Hinzen committed Jun 07, 2018 185 186 187  if ( potdenType /= POTDEN_TYPE_POTYUK ) then  Miriam Hinzen committed Jun 06, 2018 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  Miriam Hinzen committed Jun 11, 2018 193 194  end if  Miriam Hinzen committed Jun 06, 2018 195 196 197  end subroutine vmts end module m_vmts