vmts.F90 7.44 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 06, 2018 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28  subroutine vmts( input, mpi, stars, sphhar, atoms, sym, cell, oneD, vpw, rho, yukawa_residual, vr ) !------------------------------------------------------------------------- ! 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 ! termsR. !-------------------------------------------------------------------------  Markus Betzinger committed Apr 26, 2016 29 30  #include"cpp_double.h"  Miriam Hinzen committed Jun 06, 2018 31 32 33 34 35 36  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 37  use m_SphBessel  Miriam Hinzen committed Jun 06, 2018 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63  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) logical, intent(in) :: yukawa_residual 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 64  !$COMPLEX vtl_loc(0:sphhar%nlhd,atoms%ntype)  Markus Betzinger committed Apr 26, 2016 65 #ifdef CPP_MPI  Miriam Hinzen committed Jun 06, 2018 66 67 68  include 'mpif.h' integer :: ierr(3) complex, allocatable :: c_b(:)  Gregor Michalicek committed May 29, 2018 69   Miriam Hinzen committed Jun 06, 2018 70  external MPI_REDUCE  Markus Betzinger committed Apr 26, 2016 71 #endif  Uliana Alekseeva committed Mar 02, 2018 72  integer :: OMP_GET_NUM_THREADS, OMP_GET_THREAD_NUM  Gregor Michalicek committed May 29, 2018 73 74   Miriam Hinzen committed Jun 06, 2018 75 76 77 78 79  !------------------------------------------------------------------------- ! sphere boundary contribution to the coefficients calculated from the values ! of the interstitial Coulomb / Yukawa potential on the sphere boundary vtl(:,:) = cmplx( 0.0, 0.0 )  Gregor Michalicek committed May 29, 2018 80 81  #ifdef CPP_MPI  Miriam Hinzen committed Jun 06, 2018 82  call MPI_BCAST( vpw, size(vpw), CPP_MPI_COMPLEX, 0, mpi, ierr )  Markus Betzinger committed Apr 26, 2016 83 84 #endif  Miriam Hinzen committed Jun 06, 2018 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104  ! 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  Uliana Alekseeva committed Mar 02, 2018 105   Miriam Hinzen committed Jun 06, 2018 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127  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 128  !$vtl = vtl + vtl_loc  Miriam Hinzen committed Jun 06, 2018 129 130  !$omp end critical !\$omp end parallel  Markus Betzinger committed Apr 26, 2016 131 #ifdef CPP_MPI  Miriam Hinzen committed Jun 06, 2018 132 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/) ) deallocate (c_b)  Markus Betzinger committed Apr 26, 2016 137 138 #endif  Miriam Hinzen committed Jun 06, 2018 139 140 141 142 143 144 145 146  !------------------------------------------------------------------------- ! sphere interior contribution to the coefficients calculated from the ! values of the sphere Coulomb/Yukawa potential on the sphere boundary if ( yukawa_residual ) then allocate( il(0:atoms%lmaxd, 1:atoms%jmtd), kl(0:atoms%lmaxd, 1:atoms%jmtd) ) end if  Markus Betzinger committed Apr 26, 2016 147  nat = 1  Miriam Hinzen committed Jun 06, 2018 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190  do n = 1, atoms%ntype nd = atoms%ntypsy(nat) imax = atoms%jri(n) lmax = sphhar%llh(sphhar%nlh(nd), nd) if ( yukawa_residual ) then 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) if ( yukawa_residual ) then 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) vr(1:imax,lh,n) = green_factor * ( green_1(1:imax) * ( termsR - integral_2(1:imax) ) & + green_2(1:imax) * integral_1(1:imax) ) end do nat = nat + atoms%neq(n) end do if ( yukawa_residual ) then deallocate( il, kl ) end if if ( .not. yukawa_residual ) then 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 end subroutine vmts end module m_vmts