Commit 934aafd6 authored by Miriam Hinzen's avatar Miriam Hinzen

Extend functionality of vmts to cover the Yukawa case

parent 6085af3a
......@@ -178,7 +178,7 @@ CONTAINS
#ifdef CPP_MPI
CALL MPI_BCAST( den%mt, atoms%jmtd*(1+sphhar%nlhd) * atoms%ntype * dimension%jspd, MPI_DOUBLE_PRECISION, 0, mpi%mpi_comm, ierr )
#endif
CALL vmts( mpi, stars, sphhar, atoms, sym, cell, oneD, vCoul%pw(:,ispin), den%mt(:,0:,:,ispin), &
CALL vmts( input, mpi, stars, sphhar, atoms, sym, cell, oneD, vCoul%pw(:,ispin), den%mt(:,0:,:,ispin), &
yukawa_residual, vCoul%mt(:,0:,:,ispin) )
CALL timestop( "MT-spheres" )
......
MODULE m_vmts
! *******************************************************************
! this subroutine calculates the lattice harmonics expansion coeffi-*
! cients of the coulomb potential for all atom types *
! c.l.fu, r.podloucky *
! *******************************************************************
module m_vmts
CONTAINS
contains
SUBROUTINE vmts( mpi, stars, sphhar, atoms, sym, cell, oneD, vpw, rho, yukawa_residual, vr )
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.
!-------------------------------------------------------------------------
#include"cpp_double.h"
USE m_constants
USE m_types
USE m_intgr, ONLY : intgr2
USE m_phasy1
USE m_sphbes
USE m_od_phasy
use m_constants
use m_types
use m_intgr, only : intgr2
use m_phasy1
use m_sphbes
use m_od_phasy
use m_SphBessel
IMPLICIT NONE
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
REAL :: rmt2l, rmtl, ror, rr, rrlr, fpl21
INTEGER :: i, jm, k, l, l21, lh, n, nd, lm, n1, nat, m
COMPLEX :: vtl(0:sphhar%nlhd,atoms%ntype)
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
!$ COMPLEX vtl_loc(0:sphhar%nlhd,atoms%ntype)
COMPLEX :: pylm(( atoms%lmaxd + 1 ) ** 2, atoms%ntype)
REAL :: f1r(atoms%jmtd),f2r(atoms%jmtd),x1r(atoms%jmtd),x2r(atoms%jmtd)
REAL :: sbf(0:atoms%lmaxd),rrl(atoms%jmtd),rrl1(atoms%jmtd)
#ifdef CPP_MPI
INCLUDE 'mpif.h'
INTEGER :: ierr(3)
COMPLEX, ALLOCATABLE :: c_b(:)
include 'mpif.h'
integer :: ierr(3)
complex, allocatable :: c_b(:)
! .. External Subroutines
EXTERNAL MPI_REDUCE
external MPI_REDUCE
#endif
integer :: OMP_GET_NUM_THREADS, OMP_GET_THREAD_NUM
! calculate lattice harmonics expansion coefficients of the
! interstitial coulomb potential on the sphere boundaries
vtl(:,:) = CMPLX(0.0,0.0)
!-------------------------------------------------------------------------
! 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 )
#ifdef CPP_MPI
CALL MPI_BCAST(vpw,SIZE(vpw),CPP_MPI_COMPLEX,0,mpi,ierr)
call MPI_BCAST( vpw, size(vpw), CPP_MPI_COMPLEX, 0, mpi, ierr )
#endif
! g=0 component
IF (mpi%irank == 0) THEN
DO n = 1,atoms%ntype
vtl(0,n) = sfp_const*vpw(1)
ENDDO
ENDIF
! g.ne.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
!-odim
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)
!+odim
END IF
! 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.)
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)
ENDDO
!$ IF (.false.) THEN
vtl(lh,n) = vtl(lh,n) + cp*sbf(l)*sm
!$ ENDIF
!$ vtl_loc(lh,n) = vtl_loc(lh,n) + cp*sbf(l)*sm
ENDDO
nat = nat + atoms%neq(n)
ENDDO
ENDDO
!$OMP END DO
!$OMP CRITICAL
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
!$ vtl = vtl + vtl_loc
!$OMP END CRITICAL
!$OMP END PARALLEL
!$omp end critical
!$omp end parallel
#ifdef CPP_MPI
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.EQ.0) vtl=reshape(c_b,(/sphhar%nlhd+1,atoms%ntype/))
DEALLOCATE (c_b)
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)
#endif
! ----> solution of the poisson's equation
!-------------------------------------------------------------------------
! 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
nat = 1
DO n = 1,atoms%ntype
nd = atoms%ntypsy(nat)
DO lh = 0,sphhar%nlh(nd)
l = sphhar%llh(lh,nd)
l21 = 2*l + 1
fpl21 = fpi_const/l21
DO i = 1,atoms%jri(n)
rrl(i) = atoms%rmsh(i,n)**l
rrl1(i) = 1./( rrl(i) * atoms%rmsh(i,n) )
x1r(i) = rrl(i)*rho(i,lh,n)
x2r(i) = rrl1(i)*rho(i,lh,n)
ENDDO
CALL intgr2(x1r,atoms%rmsh(1,n),atoms%dx(n),atoms%jri(n),f1r)
CALL intgr2(x2r,atoms%rmsh(1,n),atoms%dx(n),atoms%jri(n),f2r)
rmtl = 1./atoms%rmt(n)**l
rmt2l = 1./atoms%rmt(n)**l21
DO i = 1,atoms%jri(n)
rrlr = rrl(i)*rmt2l
ror = rrl(i)*rmtl
vr(i,lh,n) = fpl21 * (rrl1(i)*f1r(i)-rrlr*f1r(atoms%jri(n))+&
rrl(i) * (f2r(atoms%jri(n))-f2r(i))) + ror*vtl(lh,n)
ENDDO
ENDDO
nat = nat + atoms%neq(n)
ENDDO
DO n = 1,atoms%ntype
DO i = 1,atoms%jri(n)
rr = atoms%rmsh(i,n)/atoms%rmt(n)
vr(i,0,n) = vr(i,0,n)-sfp_const*(1.-rr)/atoms%rmsh(i,n)*atoms%zatom(n)
ENDDO
ENDDO
END SUBROUTINE vmts
END MODULE m_vmts
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
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment