Skip to content

GitLab

fleur
fleur
Commit a25ca042 authored by Miriam Hinzen's avatar Miriam Hinzen
Browse files
Options

Modify Pseudo-Charge Generation Subroutines to Cover the Preconditioning Case

parent 5c698212
Show whitespace changes
Inline Side-by-side
Showing with 401 additions and 364 deletions
vgen/mpmom.F90
View file @ a25ca042
MODULE m_mpmom module m_mpmom
! *********************************************************** ! ***********************************************************
! determine the multipole moments of (original charge minus ! calculation of the multipole moments of the original charge
! plane wave charge) for each atom type ! density minus the interstitial charge density
! c.l.fu ! for each atom type
! cf. M.Weinert J.Math.Phys. 22(11) (1981) p.2434 eq. (10)-(15) !
! For yukawa_residual = .true. the subroutines calculate the
! multipole moments for the Yukawa potential instead of the
! Coulomb potential. This is used in the preconditioning of
! the SCF iteration for metallic systems.
! !
! qlmo(m,l,n) : mult.mom. of the mufftn-tin charge density ! qlmo(m,l,n) : mult.mom. of the mufftn-tin charge density
! qlmp(m,l,n) : mult.mom. of the plane-wave charge density ! qlmp(m,l,n) : mult.mom. of the plane-wave charge density
! qlm (m,l,n) : (output) difference of the former quantities ! qlm (m,l,n) : (output) difference of the former quantities
! !
! references:
! for both the Coulomb and the Yukawa pseudo charge density:
! F. Tran, P. Blaha: Phys. Rev. B 83, 235118 (2011)
! or see the original paper for the normal Coulomb case only:
! M. Weinert: J.Math.Phys. 22(11) (1981) p.2434 eq. (10)-(15)
! *********************************************************** ! ***********************************************************
CONTAINS contains
SUBROUTINE mpmom( mpi, atoms, sphhar, stars, sym, cell, oneD, qpw, rho, yukawa_residual, qlm ) subroutine mpmom( input, mpi, atoms, sphhar, stars, sym, cell, oneD, qpw, rho, yukawa_residual, qlm )
USE m_types
IMPLICIT NONE
TYPE(t_mpi), INTENT(IN) :: mpi
TYPE(t_oneD), INTENT(IN) :: oneD
TYPE(t_sym), INTENT(IN) :: sym
TYPE(t_stars), INTENT(IN) :: stars
TYPE(t_cell), INTENT(IN) :: cell
TYPE(t_sphhar), INTENT(IN) :: sphhar
TYPE(t_atoms), INTENT(IN) :: atoms
REAL, INTENT(IN) :: rho(:,0:,:) !(atoms%jmtd,0:sphhar%nlhd,atoms%ntype)
COMPLEX, INTENT(IN) :: qpw(:) !(stars%ng3)
logical, intent(in) :: yukawa_residual
COMPLEX, INTENT (OUT) :: qlm(-atoms%lmaxd:atoms%lmaxd,0:atoms%lmaxd,atoms%ntype)
INTEGER :: j, jm, lh, mb, mem, mems, n, nd, l, nat, m use m_types
COMPLEX :: qlmo(-atoms%lmaxd:atoms%lmaxd,0:atoms%lmaxd,atoms%ntype) implicit none
COMPLEX :: qlmp(-atoms%lmaxd:atoms%lmaxd,0:atoms%lmaxd,atoms%ntype)
! multipole moments of original charge (q_{lm}^i) type(t_input), intent(in) :: input
IF (mpi%irank == 0) THEN type(t_mpi), intent(in) :: mpi
CALL mt_moments( atoms, sphhar, rho(:,:,:), yukawa_residual, qlmo ) type(t_oneD), intent(in) :: oneD
ENDIF ! mpi%irank == 0 type(t_sym), intent(in) :: sym
type(t_stars), intent(in) :: stars
type(t_cell), intent(in) :: cell
type(t_sphhar), intent(in) :: sphhar
type(t_atoms), intent(in) :: atoms
real, intent(in) :: rho(:,0:,:) !(atoms%jmtd,0:sphhar%nlhd,atoms%ntype)
complex, intent(in) :: qpw(:) !(stars%ng3)
logical, intent(in) :: yukawa_residual
complex, intent(out) :: qlm(-atoms%lmaxd:atoms%lmaxd,0:atoms%lmaxd,atoms%ntype)
integer :: j, jm, lh, mb, mem, mems, n, nd, l, nat, m
complex :: qlmo(-atoms%lmaxd:atoms%lmaxd,0:atoms%lmaxd,atoms%ntype)
complex :: qlmp(-atoms%lmaxd:atoms%lmaxd,0:atoms%lmaxd,atoms%ntype)
CALL pw_moments( mpi, stars, atoms, cell, sym, oneD, qpw(:), yukawa_residual, qlmp ) ! multipole moments of original charge density
if ( mpi%irank == 0 ) then
call mt_moments( input, atoms, sphhar, rho(:,:,:), yukawa_residual, qlmo )
end if
! eq.(15): \tilde q_(lm}^i = q_{lm}^i - q_{lm}^{Ii} ! multipole moments of the interstitial charge density in the spheres
IF (mpi%irank == 0) THEN call pw_moments( input, mpi, stars, atoms, cell, sym, oneD, qpw(:), yukawa_residual, qlmp )
qlm=qlmo-qlmp
! Output section if ( mpi%irank == 0 ) then
! see (A14)
qlm = qlmo - qlmp
! output section
nat = 1 nat = 1
DO n = 1,atoms%ntype do n = 1, atoms%ntype
WRITE (6,FMT=8000) n write( 6, fmt=8000 ) n
nd = atoms%ntypsy(nat) nd = atoms%ntypsy(nat)
DO lh = 0,sphhar%nlh(nd) do lh = 0, sphhar%nlh(nd)
l = sphhar%llh(lh,nd) l = sphhar%llh(lh,nd)
mems = sphhar%nmem(lh,nd) mems = sphhar%nmem(lh,nd)
DO mem = 1,mems do mem = 1, mems
m = sphhar%mlh(mem,lh,nd) m = sphhar%mlh(mem,lh,nd)
WRITE (6,FMT=8010) l,m,qlmo(m,l,n),qlmp(m,l,n) write( 6, fmt=8010 ) l, m, qlmo(m,l,n), qlmp(m,l,n)
! write(16,1002) l,m,qlmo(m,l,n),qlmp(m,l,n) end do
ENDDO end do
ENDDO
nat = nat + atoms%neq(n) nat = nat + atoms%neq(n)
ENDDO end do
!
8000 FORMAT (/,10x,'multipole moments for atom type=',i5,/,/,t3,'l',t7,& 8000 FORMAT (/,10x,'multipole moments for atom type=',i5,/,/,t3,'l',t7,&
& 'm',t27,'original',t57,'plane wave') & 'm',t27,'original',t57,'plane wave')
8010 FORMAT (1x,i2,2x,i3,2x,2 (5x,2e15.5)) 8010 FORMAT (1x,i2,2x,i3,2x,2 (5x,2e15.5))
! !
ENDIF ! mpi%irank == 0 end if ! mpi%irank == 0
END SUBROUTINE mpmom end subroutine mpmom
SUBROUTINE mt_moments( atoms, sphhar, rho, yukawa_residual, qlmo ) subroutine mt_moments( input, atoms, sphhar, rho, yukawa_residual, qlmo )
!multipole moments of original charge (q_{lm}^i) ! multipole moments of original charge density
! see (A15) (Coulomb case) or (A17) (Yukawa case)
USE m_intgr, ONLY : intgr3 use m_intgr, only: intgr3
USE m_constants,ONLY:sfp_const use m_constants, only: sfp_const
USE m_types use m_types
use m_DoubleFactorial use m_DoubleFactorial
use m_SphBessel use m_SphBessel
IMPLICIT NONE implicit none
TYPE(t_sphhar), INTENT(IN) :: sphhar type(t_input), intent(in) :: input
TYPE(t_atoms), INTENT(IN) :: atoms type(t_sphhar), intent(in) :: sphhar
REAL, INTENT(IN) :: rho(: ,0:, :) type(t_atoms), intent(in) :: atoms
real, intent(in) :: rho(: ,0:, :)
logical, intent(in) :: yukawa_residual logical, intent(in) :: yukawa_residual
COMPLEX, INTENT(OUT) :: qlmo(-atoms%lmaxd:,0:,:) complex, intent(out) :: qlmo(-atoms%lmaxd:,0:,:)
integer :: n, ns, jm, nl, l, j, mb, m, nat, i, imax, lmax
real :: fint
real :: f( maxval( atoms%jri ) )
real, allocatable, dimension(:,:) :: il, kl
INTEGER :: n, ns, jm, nl, l, j, mb, m, nat if ( yukawa_residual ) then
REAL :: fint allocate( il(0:atoms%lmaxd, 1:atoms%jmtd), kl(0:atoms%lmaxd, 1:atoms%jmtd) )
REAL :: f(MAXVAL(atoms%jri)) end if
qlmo=0.0 qlmo = 0.0
nat = 1 nat = 1
DO n = 1, atoms%ntype do n = 1, atoms%ntype
ns = atoms%ntypsy(nat) ns = atoms%ntypsy(nat)
jm = atoms%jri(n) jm = atoms%jri(n)
DO nl = 0, sphhar%nlh(ns) imax = atoms%jri(n)
lmax = sphhar%llh(sphhar%nlh(ns), ns)
if ( yukawa_residual ) then
do concurrent (i = 1:imax)
call ModSphBessel( il(:, i), kl(:, i), input%preconditioning_param * atoms%rmsh(i, n), lmax )
end do
end if
do nl = 0, sphhar%nlh(ns)
l = sphhar%llh(nl,ns) l = sphhar%llh(nl,ns)
DO j = 1, jm do j = 1, jm
f(j) = (atoms%rmsh(j,n)**l)*rho(j,nl,n) if ( .not. yukawa_residual ) then
ENDDO f(j) = atoms%rmsh(j,n) ** l * rho(j,nl,n)
CALL intgr3(f,atoms%rmsh(:,n),atoms%dx(n),jm,fint) else
DO mb = 1, sphhar%nmem(nl,ns) f(j) = il(l, j) * rho(j,nl,n)
end if
end do
call intgr3( f, atoms%rmsh(:,n), atoms%dx(n), jm, fint )
if ( yukawa_residual ) then
fint = fint * DoubleFactorial( l ) / input%preconditioning_param ** l
end if
do mb = 1, sphhar%nmem(nl,ns)
m = sphhar%mlh(mb,nl,ns) m = sphhar%mlh(mb,nl,ns)
qlmo(m,l,n) = qlmo(m,l,n) + sphhar%clnu(mb,nl,ns)*fint qlmo(m,l,n) = qlmo(m,l,n) + sphhar%clnu(mb,nl,ns) * fint
ENDDO end do
ENDDO end do
qlmo(0,0,n) = qlmo(0,0,n) - atoms%zatom(n)/sfp_const if ( .not. yukawa_residual ) then
qlmo(0,0,n) = qlmo(0,0,n) - atoms%zatom(n) / sfp_const
end if
nat = nat + atoms%neq(n) nat = nat + atoms%neq(n)
ENDDO end do
END SUBROUTINE mt_moments end subroutine mt_moments
SUBROUTINE pw_moments( mpi, stars, atoms, cell, sym, oneD, qpw_in, yukawa_residual, qlmp_out ) subroutine pw_moments( input, mpi, stars, atoms, cell, sym, oneD, qpw_in, yukawa_residual, qlmp_out )
!multipole moments of plane wave charge inside the spheres (q_{lm}^{Ii}) ! multipole moments of the interstitial charge in the spheres
USE m_phasy1 use m_phasy1
USE m_sphbes use m_sphbes
USE m_od_phasy use m_od_phasy
USE m_constants,ONLY:sfp_const use m_constants, only: sfp_const
USE m_types use m_types
use m_DoubleFactorial use m_DoubleFactorial
use m_SphBessel use m_SphBessel
IMPLICIT NONE implicit none
TYPE(t_mpi), INTENT(IN) :: mpi type(t_input), intent(in) :: input
TYPE(t_oneD), INTENT(IN) :: oneD type(t_mpi), intent(in) :: mpi
TYPE(t_sym), INTENT(IN) :: sym type(t_oneD), intent(in) :: oneD
TYPE(t_stars), INTENT(IN) :: stars type(t_sym), intent(in) :: sym
TYPE(t_cell), INTENT(IN) :: cell type(t_stars), intent(in) :: stars
TYPE(t_atoms), INTENT(IN) :: atoms type(t_cell), intent(in) :: cell
COMPLEX, INTENT(IN) :: qpw_in(:) type(t_atoms), intent(in) :: atoms
complex, intent(in) :: qpw_in(:)
logical, intent(in) :: yukawa_residual logical, intent(in) :: yukawa_residual
COMPLEX, INTENT(OUT) :: qlmp_out(-atoms%lmaxd:,0:,:) complex, intent(out) :: qlmp_out(-atoms%lmaxd:,0:,:)
INTEGER :: n, k, l, ll1, lm, ierr(3), m integer :: n, k, l, ll1, lm, ierr(3), m
COMPLEX :: sk3i, cil, nqpw complex :: sk3i, cil, nqpw
COMPLEX :: pylm(( MAXVAL(atoms%lmax) + 1 ) ** 2, atoms%ntype) complex :: pylm(( maxval( atoms%lmax ) + 1 ) ** 2, atoms%ntype)
REAL :: sk3r, rl3 real :: sk3r, rl2
REAL :: aj(0:MAXVAL(atoms%lmax)+1) real :: aj(0:maxval( atoms%lmax ) + 1 )
COMPLEX :: qpw(stars%ng3) complex :: qpw(stars%ng3)
LOGICAL :: od logical :: od
real :: il(0:maxval( atoms%lmax ) + 1 )
real :: kl(0:maxval( atoms%lmax ) + 1 )
#ifdef CPP_MPI #ifdef CPP_MPI
INCLUDE 'mpif.h' include 'mpif.h'
#endif #endif
COMPLEX :: qlmp(-atoms%lmaxd:atoms%lmaxd,0:atoms%lmaxd,atoms%ntype) complex :: qlmp(-atoms%lmaxd:atoms%lmaxd,0:atoms%lmaxd,atoms%ntype)
qpw=qpw_in(:stars%ng3) qpw = qpw_in(:stars%ng3)
qlmp= 0.0 qlmp = 0.0
IF (mpi%irank==0) THEN if ( mpi%irank == 0 ) then
!g eq 0 term : \sqrt{4 \pi}/3 R_i^3 \rho_I(0) \delta_{l,0} ! q=0 term: see (A19) (Coulomb case) or (A20) (Yukawa case)
do n = 1, atoms%ntype
DO n = 1,atoms%ntype if ( .not. yukawa_residual ) then
qlmp(0,0,n) = qpw(1)*stars%nstr(1)*atoms%volmts(n)/sfp_const qlmp(0,0,n) = qpw(1) * stars%nstr(1) * atoms%volmts(n) / sfp_const
ENDDO else
ENDIF call ModSphBessel( il(0:1), kl(0:1), input%preconditioning_param * atoms%rmt(n), 1 )
qlmp(0,0,n) = qpw(1) * stars%nstr(1) * sfp_const * atoms%rmt(n) ** 2 * il(1) / input%preconditioning_param
end if
end do
end if
#ifdef CPP_MPI #ifdef CPP_MPI
CALL MPI_BCAST(qpw,SIZE(qpw),MPI_DOUBLE_COMPLEX,0, mpi%mpi_comm,ierr) call MPI_BCAST( qpw, size(qpw), MPI_DOUBLE_COMPLEX, 0, mpi%mpi_comm, ierr )
#endif #endif
! g ne 0 terms : \sum_{K \= 0} 4 \pi i^l \rho_I(K) R_i^{l+3} \times
! j_{l+1} (KR_i) / KR_i \exp{iK\xi_i} Y^*_{lm} (K) ! q/=0 terms: see (A16) (Coulomb case) or (A18) (Yukawa case)
od=oneD%odi%d1 od = oneD%odi%d1
! !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(pylm,nqpw,n,sk3r,aj,rl3,sk3i,& !$omp parallel do default( shared ) private( pylm, nqpw, n, sk3r, aj, rl3, sk3i, &
! !$OMP& l,cil,ll1,m,lm,k) REDUCTION(+:qlmp) !$omp& l, cil, ll1, m, lm, k ) reduction( +:qlmp )
DO k = mpi%irank+2, stars%ng3, mpi%isize do k = mpi%irank+2, stars%ng3, mpi%isize
IF (od) THEN if ( od ) then
CALL od_phasy(atoms%ntype,stars%ng3,atoms%nat,atoms%lmaxd,atoms%ntype,atoms%neq,atoms%lmax,& call od_phasy( atoms%ntype, stars%ng3, atoms%nat, atoms%lmaxd, atoms%ntype, &
atoms%taual,cell%bmat,stars%kv3,k,oneD%odi,oneD%ods,& atoms%neq, atoms%lmax, atoms%taual, cell%bmat, stars%kv3, k, oneD%odi, oneD%ods, pylm)
pylm) else
ELSE call phasy1( atoms, stars, sym, cell, k, pylm )
CALL phasy1(atoms,stars,sym,cell,k,pylm) end if
END IF
! nqpw = qpw(k) * stars%nstr(k)
nqpw = qpw(k)*stars%nstr(k) do n = 1, atoms%ntype
DO n = 1,atoms%ntype sk3r = stars%sk3(k) * atoms%rmt(n)
sk3r = stars%sk3(k)*atoms%rmt(n) call sphbes( atoms%lmax(n) + 1, sk3r, aj )
CALL sphbes(atoms%lmax(n)+1,sk3r,aj) rl2 = atoms%rmt(n) ** 2
rl3 = atoms%rmt(n)**3 if ( yukawa_residual ) then
sk3i = nqpw/sk3r call ModSphBessel( il(0:atoms%lmax(n)+1), kl(0:atoms%lmax(n)+1), input%preconditioning_param * atoms%rmt(n), atoms%lmax(n) + 1 )
DO l = 0,atoms%lmax(n) sk3i = nqpw / ( stars%sk3(k) ** 2 + input%preconditioning_param ** 2 ) * rl2
cil = aj(l+1)*sk3i*rl3 else
ll1 = l*(l+1) + 1 sk3i = nqpw / stars%sk3(k)
DO m = -l,l end if
do l = 0, atoms%lmax(n)
if ( yukawa_residual ) then
cil = ( stars%sk3(k) * il(l) * aj(l+1) + input%preconditioning_param * il(l+1) * aj(l) ) * ( DoubleFactorial( l ) / input%preconditioning_param ** l ) * sk3i
else
cil = aj(l+1) * sk3i * rl2
rl2 = rl2 * atoms%rmt(n)
end if
ll1 = l * ( l + 1 ) + 1
do m = -l, l
lm = ll1 + m lm = ll1 + m
qlmp(m,l,n) = qlmp(m,l,n) + cil*pylm(lm,n) qlmp(m,l,n) = qlmp(m,l,n) + cil * pylm(lm,n)
ENDDO end do
rl3 = rl3*atoms%rmt(n) end do ! l = 0, atoms%lmax(n)
ENDDO ! l = 0, atoms%lmax(n) end do ! n = 1, atoms%ntype
ENDDO ! n = 1, atoms%ntype end do ! k = 2, stars%ng3
ENDDO ! k = 2, stars%ng3 !$omp end parallel do
! !$OMP END PARALLEL DO
#ifdef CPP_MPI #ifdef CPP_MPI
PRINT *,"mpi",mpi%irank,qlmp(0,0,:) print *, "mpi", mpi%irank, qlmp(0,0,:)
CALL MPI_REDUCE(qlmp,qlmp_out,SIZE(qlmp),MPI_DOUBLE_COMPLEX,MPI_SUM,0,mpi%mpi_comm,ierr) call MPI_REDUCE( qlmp, qlmp_out, size(qlmp), MPI_DOUBLE_COMPLEX, MPI_SUM, 0, mpi%mpi_comm, ierr )
#else #else
qlmp_out=qlmp qlmp_out = qlmp
#endif #endif
END SUBROUTINE pw_moments end subroutine pw_moments
END MODULE m_mpmom end module m_mpmom
vgen/psqpw.F90
View file @ a25ca042
MODULE m_psqpw module m_psqpw
! *********************************************************** ! ***********************************************************
! generates the fourier coefficients of pseudo charge density ! generates the fourier coefficients of pseudo charge density
! c.l.fu
! corrected april 1990 m.w.
! !
! cf. M.Weinert J.Math.Phys. 22(11) (1981) p.2434 eq. (10)-(15) ! For yukawa_residual = .true. the subroutines calculate the
! pseudo charge density for the generation of the Yukawa
! potential instead of the Coulomb potential. This is used in
! the preconditioning of the SCF iteration for metallic systems.
! !
! ! references:
! parallelized 04/08 gb ! for both the Coulomb and Yukawa cases:
! F. Tran, P. Blaha: Phys. Rev. B 83, 235118 (2011)
! or see the original paper for the normal Coulomb case only:
! M. Weinert: J. Math. Phys. 22(11) (1981) p.2434 eq. (10)-(15)
! *********************************************************** ! ***********************************************************
CONTAINS contains
SUBROUTINE psqpw( mpi, atoms, sphhar, stars, vacuum, DIMENSION, cell, input, sym, oneD, & subroutine psqpw( mpi, atoms, sphhar, stars, vacuum, dimension, cell, input, sym, oneD, &
& qpw, rho, rht, l_xyav, yukawa_residual, psq ) & qpw, rho, rht, l_xyav, yukawa_residual, psq )
#include"cpp_double.h" #include"cpp_double.h"
USE m_constants use m_constants
USE m_phasy1 use m_phasy1
USE m_mpmom use m_mpmom
USE m_sphbes use m_sphbes
USE m_qsf use m_qsf
USE m_od_phasy use m_od_phasy
USE m_od_cylbes use m_od_cylbes
USE m_types use m_types
use m_DoubleFactorial use m_DoubleFactorial
use m_SphBessel use m_SphBessel
IMPLICIT NONE implicit none
TYPE(t_mpi), INTENT(IN) :: mpi type(t_mpi), intent(in) :: mpi
TYPE(t_atoms), INTENT(IN) :: atoms type(t_atoms), intent(in) :: atoms
TYPE(t_sphhar), INTENT(IN) :: sphhar type(t_sphhar), intent(in) :: sphhar
TYPE(t_stars), INTENT(IN) :: stars type(t_stars), intent(in) :: stars
TYPE(t_vacuum), INTENT(IN) :: vacuum type(t_vacuum), intent(in) :: vacuum
TYPE(t_dimension), INTENT(IN) :: DIMENSION type(t_dimension), intent(in) :: dimension
TYPE(t_cell), INTENT(IN) :: cell type(t_cell), intent(in) :: cell
TYPE(t_input), INTENT(IN) :: input type(t_input), intent(in) :: input
TYPE(t_sym), INTENT(IN) :: sym type(t_sym), intent(in) :: sym
TYPE(t_oneD), INTENT(IN) :: oneD type(t_oneD), intent(in) :: oneD
LOGICAL, INTENT(IN) :: l_xyav logical, intent(in) :: l_xyav
COMPLEX, INTENT(IN) :: qpw(stars%ng3) complex, intent(in) :: qpw(stars%ng3)
REAL, INTENT(IN) :: rho(atoms%jmtd,0:sphhar%nlhd,atoms%ntype) real, intent(in) :: rho(atoms%jmtd,0:sphhar%nlhd,atoms%ntype)
REAL, INTENT(IN) :: rht(vacuum%nmzd,2) real, intent(in) :: rht(vacuum%nmzd,2)
logical, intent(in) :: yukawa_residual logical, intent(in) :: yukawa_residual
COMPLEX, INTENT(OUT) :: psq(stars%ng3) complex, intent(out) :: psq(stars%ng3)
COMPLEX :: psint, sa, sl, sm complex :: psint, sa, sl, sm
REAL :: f, fact, fpo, gz, p, qvac, rmtl, s, fJ, gr, g real :: f, fact, fpo, gz, p, qvac, rmtl, s, fJ, gr, g
INTEGER :: ivac, k, l, n, n1, nc, ncvn, lm, ll1, nd, m, nz integer :: ivac, k, l, n, n1, nc, ncvn, lm, ll1, nd, m, nz
COMPLEX :: pylm(( atoms%lmaxd + 1 ) ** 2, atoms%ntype) complex :: pylm(( atoms%lmaxd + 1 ) ** 2, atoms%ntype)
COMPLEX :: qlm(-atoms%lmaxd:atoms%lmaxd,0:atoms%lmaxd,atoms%ntype) complex :: qlm(-atoms%lmaxd:atoms%lmaxd,0:atoms%lmaxd,atoms%ntype)
REAL :: q2(vacuum%nmzd) real :: q2(vacuum%nmzd)
real :: pn(0:atoms%lmaxd,atoms%ntype) real :: pn(0:atoms%lmaxd,atoms%ntype)
real :: aj(0:atoms%lmaxd+DIMENSION%ncvd+1) real :: aj(0:atoms%lmaxd+DIMENSION%ncvd+1)
REAL :: rht1(vacuum%nmz) real :: rht1(vacuum%nmz)
real, allocatable, dimension(:) :: il, kl
real :: g0(atoms%ntype)
#ifdef CPP_MPI #ifdef CPP_MPI
INCLUDE 'mpif.h' include 'mpif.h'
INTEGER :: ierr(3) integer :: ierr(3)
COMPLEX, ALLOCATABLE :: c_b(:) complex, allocatable :: c_b(:)
#endif #endif
! Calculate multipole moments ! Calculate multipole moments
CALL mpmom( mpi, atoms, sphhar, stars, sym, cell, oneD, qpw, rho, yukawa_residual, qlm ) call mpmom( input, mpi, atoms, sphhar, stars, sym, cell, oneD, qpw, rho, yukawa_residual, qlm )
#ifdef CPP_MPI #ifdef CPP_MPI
psq(:) = CMPLX(0.0,0.0) psq(:) = cmplx( 0.0, 0.0 )
CALL MPI_BCAST(qpw,size(qpw),CPP_MPI_COMPLEX,0,mpi%mpi_comm,ierr) call MPI_BCAST( qpw, size(qpw), CPP_MPI_COMPLEX, 0, mpi%mpi_comm, ierr )
nd = (2*atoms%lmaxd+1)*(atoms%lmaxd+1)*atoms%ntype nd = ( 2 * atoms%lmaxd + 1 ) * ( atoms%lmaxd + 1 ) * atoms%ntype
CALL MPI_BCAST(qlm,nd,CPP_MPI_COMPLEX,0,mpi%MPI_COMM,ierr) call MPI_BCAST( qlm, nd, CPP_MPI_COMPLEX, 0, mpi%MPI_COMM, ierr )
#endif #endif
!
! pn(l,n) = (2l + 2nc(n) + 3)!! / (2l + 1)!! R^l ; ncv(n)=n+l in paper ! prefactor in (A10) (Coulomb case) or (A11) (Yukawa case)
! ! nc(n) is the integer p in the paper; ncv(n) is l + p
DO n = 1,atoms%ntype ! Coulomb case: pn(l,n) = (2 * l + 2 * p + 3)!! / ( (2 * l + 1)!! * R ** (ncv(n) + 1) ),
rmtl = 1.0 ! Yukawa case: pn(l,n) = lambda ** (l + p + 1) / ( i_{l+p+1}(lambda * R) * (2 * l + 1)!! )
DO l = 0,atoms%lmax(n) ! g0 is the prefactor for the q=0 component in (A13)
IF (l.GE.atoms%ncv(n)) THEN pn = 0.
pn(l,n) = 0.0 do n = 1, atoms%ntype
ELSE if ( .not. yukawa_residual ) then
p = 1. do l = 0, min( atoms%ncv(n) - 1, atoms%lmax(n) )
DO nc = l,atoms%ncv(n) pn(l, n) = DoubleFactorial( atoms%ncv(n) + 1, l ) / ( atoms%rmt(n) ** ( atoms%ncv(n) + 1 ) )
p = p* (2*nc+3) end do