Commit 2cc2779b authored by S.Rost's avatar S.Rost

Merge branch 'develop' of iffgit.fz-juelich.de:fleur/fleur into develop

parents 3c90dc6c eb35c2e9
...@@ -17,7 +17,7 @@ CONTAINS ...@@ -17,7 +17,7 @@ CONTAINS
DO iintsp=1,2 DO iintsp=1,2
DO jintsp=1,2 DO jintsp=1,2
mat(jintsp,iintsp)%data_c(:,:)=chi(jintsp,iintsp)*mat_tmp%data_c(:,:)+mat(jintsp,iintsp)%data_c(:,:) mat(jintsp,iintsp)%data_c(:,:)=chi(iintsp,jintsp)*mat_tmp%data_c(:,:)+mat(jintsp,iintsp)%data_c(:,:)
ENDDO ENDDO
ENDDO ENDDO
END SUBROUTINE hsmt_distspins END SUBROUTINE hsmt_distspins
......
...@@ -75,6 +75,7 @@ contains ...@@ -75,6 +75,7 @@ contains
type(t_potden) :: resDen, vYukawa type(t_potden) :: resDen, vYukawa
integer :: ierr(2) integer :: ierr(2)
!External functions !External functions
real :: CPP_BLAS_sdot real :: CPP_BLAS_sdot
external :: CPP_BLAS_sdot external :: CPP_BLAS_sdot
...@@ -230,7 +231,7 @@ contains ...@@ -230,7 +231,7 @@ contains
! KERKER PRECONDITIONER ! KERKER PRECONDITIONER
if( input%preconditioning_param /= 0 ) then if( input%preconditioning_param /= 0 ) then
call resDen%init( stars, atoms, sphhar, vacuum, input%jspins, noco%l_noco, 1001 ) call resDen%init( stars, atoms, sphhar, vacuum, input%jspins, noco%l_noco, POTDEN_TYPE_DEN )
call vYukawa%init( stars, atoms, sphhar, vacuum, input%jspins, noco%l_noco, 4 ) call vYukawa%init( stars, atoms, sphhar, vacuum, input%jspins, noco%l_noco, 4 )
MPI0_b: if( mpi%irank == 0 ) then MPI0_b: if( mpi%irank == 0 ) then
call resDen%subPotDen( outDen, inDen ) call resDen%subPotDen( outDen, inDen )
...@@ -258,6 +259,8 @@ contains ...@@ -258,6 +259,8 @@ contains
* vYukawa%mt(1:atoms%jri(n),lh,n,1) * atoms%rmsh(1:atoms%jri(n),n) ** 2 * vYukawa%mt(1:atoms%jri(n),lh,n,1) * atoms%rmsh(1:atoms%jri(n),n) ** 2
end do end do
end do end do
resDen%vacz = resDen%vacz - input%preconditioning_param ** 2 / fpi_const * vYukawa%vacz
resDen%vacxy = resDen%vacxy - input%preconditioning_param ** 2 / fpi_const * vYukawa%vacxy
if( input%jspins == 2 ) call resDen%ChargeAndMagnetisationToSpins() if( input%jspins == 2 ) call resDen%ChargeAndMagnetisationToSpins()
! fix the preconditioned density ! fix the preconditioned density
call outDen%addPotDen( resDen, inDen ) call outDen%addPotDen( resDen, inDen )
......
...@@ -55,10 +55,9 @@ module m_SphBessel ...@@ -55,10 +55,9 @@ module m_SphBessel
! modernised and extended by M. Hinzen (2016) ! modernised and extended by M. Hinzen (2016)
!------------------------------------------------------------------------- !-------------------------------------------------------------------------
use m_constants, only: ImagUnit
implicit none implicit none
complex, parameter :: CI = (0.0, 1.0)
interface SphBessel interface SphBessel
module procedure SphBesselComplex, SphBesselReal module procedure SphBesselComplex, SphBesselReal
...@@ -105,7 +104,7 @@ contains ...@@ -105,7 +104,7 @@ contains
end do EXPANSION end do EXPANSION
jl(l) = jl(l) * zj jl(l) = jl(l) * zj
nl(l) = -nl(l) * zn nl(l) = -nl(l) * zn
hl(l) = jl(l) + nl(l) * CI hl(l) = jl(l) + nl(l) * ImagUnit
zj = zj * z / ( rl + 3.0 ) zj = zj * z / ( rl + 3.0 )
zn = zn / z * ( rl + 1.0 ) zn = zn / z * ( rl + 1.0 )
end do SERIAL_L_LOOP end do SERIAL_L_LOOP
...@@ -120,12 +119,12 @@ contains ...@@ -120,12 +119,12 @@ contains
SERIAL_M_LOOP: do m = 0, l SERIAL_M_LOOP: do m = 0, l
hl(l) = hl(l) + (-1) ** m * rnm(l) hl(l) = hl(l) + (-1) ** m * rnm(l)
nl(l) = nl(l) + rnm(l) nl(l) = nl(l) + rnm(l)
rnm(l) = rnm(l) / ( m + 1.0 ) * ( l * ( l + 1 ) - m * ( m + 1 ) ) / ( CI * ( z + z ) ) rnm(l) = rnm(l) / ( m + 1.0 ) * ( l * ( l + 1 ) - m * ( m + 1 ) ) / ( ImagUnit * ( z + z ) )
end do SERIAL_M_LOOP end do SERIAL_M_LOOP
hl(l) = hl(l) * (-CI) ** l * exp( CI * z ) / ( CI * z ) hl(l) = hl(l) * (-ImagUnit) ** l * exp( ImagUnit * z ) / ( ImagUnit * z )
nl(l) = nl(l) * CI ** l * exp( -CI * z ) / ( -CI * z ) nl(l) = nl(l) * ImagUnit ** l * exp( -ImagUnit * z ) / ( -ImagUnit * z )
jl(l) = ( hl(l) + nl(l) ) / 2.0 jl(l) = ( hl(l) + nl(l) ) / 2.0
nl(l) = ( hl(l) - jl(l) ) * (-CI) nl(l) = ( hl(l) - jl(l) ) * (-ImagUnit)
end if end if
end do PARALLEL_L_LOOP end do PARALLEL_L_LOOP
...@@ -165,11 +164,11 @@ contains ...@@ -165,11 +164,11 @@ contains
integer :: l integer :: l
call SphBesselComplex( il, nl, kl, CI * z, lmax ) call SphBesselComplex( il, nl, kl, ImagUnit * z, lmax )
do l = 0, lmax do l = 0, lmax
il(l) = (-CI) ** l * il(l) il(l) = (-ImagUnit) ** l * il(l)
kl(l) = - CI ** l * kl(l) kl(l) = - ImagUnit ** l * kl(l)
end do end do
end subroutine ModSphBesselComplex end subroutine ModSphBesselComplex
...@@ -185,11 +184,11 @@ contains ...@@ -185,11 +184,11 @@ contains
integer :: l integer :: l
call SphBesselComplex( il, nl, kl, CI * z, lmax ) call SphBesselComplex( il, nl, kl, ImagUnit * z, lmax )
do l = 0, lmax do l = 0, lmax
il(l) = (-CI) ** l * il(l) il(l) = (-ImagUnit) ** l * il(l)
kl(l) = - CI ** l * kl(l) kl(l) = - ImagUnit ** l * kl(l)
end do end do
end subroutine ModSphBesselComplex2 end subroutine ModSphBesselComplex2
...@@ -207,12 +206,12 @@ contains ...@@ -207,12 +206,12 @@ contains
complex :: z complex :: z
z = CI * x z = ImagUnit * x
call SphBesselComplex( jl, nl, hl, z, lmax ) call SphBesselComplex( jl, nl, hl, z, lmax )
do l = 0, lmax do l = 0, lmax
il(l) = (-CI) ** l * jl(l) il(l) = (-ImagUnit) ** l * jl(l)
kl(l) = - CI ** l * hl(l) kl(l) = - ImagUnit ** l * hl(l)
end do end do
end subroutine ModSphBesselReal end subroutine ModSphBesselReal
......
module m_VYukawaFilm module m_VYukawaFilm
! computation of the film-case Yukawa potential for the preconditioning of the charge density residual ! Computation of the film-case Yukawa potential for the preconditioning of the charge density residual
! 1. part: The pseudo-charge density is used for the vacuum and interstitial potentials.
! 4. part: The MT potential is calculated at the end.
! These two parts do not change compared to the bulk case.
! 2. and 3. part: The vacuum and interstitial potentials are calculated by integration of the product of the charge density and
! the 1D Green function with respect to z (from minus infinity to infinity). The integrals for both potentials are split in 4:
! the seperating points are the boundaries of the slab and the point where the integration variable z' equals the point z where we
! want to evaluate the potential for. For the contribution from the vacuum charge density we do a numerical integration. For the
! contribution from the slab charge density we have analytical expressions of the integrals.
contains contains
...@@ -65,8 +75,8 @@ module m_VYukawaFilm ...@@ -65,8 +75,8 @@ module m_VYukawaFilm
psq, rhtxy, rht, & psq, rhtxy, rht, &
VVxy, VVz, alphm ) VVxy, VVz, alphm )
! 1. part: compute the contribution from the interstitial charge density to the vacuum potential as a function of q_xy and z ! 1. part: Compute the contribution from the interstitial charge density to the vacuum potential as a function of q_xy and z (analytic expression for integral)
! 2. part: compute the contribution from the vacuum charge density to the vacuum potential as a function of q_xy and z by numerical integration ! 2. part: Compute the contribution from the vacuum charge density to the vacuum potential as a function of q_xy and z by numerical integration
use m_ExpSave use m_ExpSave
use m_constants use m_constants
...@@ -95,26 +105,27 @@ module m_VYukawaFilm ...@@ -95,26 +105,27 @@ module m_VYukawaFilm
real :: exp_m(vacuum%nmzd,stars%ng2), exp_p(vacuum%nmzd,stars%ng2) real :: exp_m(vacuum%nmzd,stars%ng2), exp_p(vacuum%nmzd,stars%ng2)
real :: expDhg(stars%ng2), expDg(stars%ng2) real :: expDhg(stars%ng2), expDg(stars%ng2)
real :: z(vacuum%nmzd) real :: z(vacuum%nmzd)
integer :: iz, irec2, irec3, ivac, kz integer :: iz, irec2, irec3, ivac, iqz
complex :: fa(vacuum%nmzxyd,2:stars%ng2), fb(vacuum%nmzxyd,2:stars%ng2) complex :: fa(vacuum%nmzxyd,2:stars%ng2), fb(vacuum%nmzxyd,2:stars%ng2)
complex :: alpha(vacuum%nmzxyd,2:stars%ng2,2), beta(vacuum%nmzxyd,2:stars%ng2,2), gamma(vacuum%nmzxyd,2:stars%ng2) complex :: alpha(vacuum%nmzxyd,2:stars%ng2,2), beta(vacuum%nmzxyd,2:stars%ng2,2), gamma(vacuum%nmzxyd,2:stars%ng2)
real :: ga(vacuum%nmzd), gb(vacuum%nmzd) real :: ga(vacuum%nmzd), gb(vacuum%nmzd)
real :: delta(vacuum%nmzd,2), epsilon(vacuum%nmzd,2), zeta(vacuum%nmzd) real :: delta(vacuum%nmzd,2), epsilon(vacuum%nmzd,2), zeta(vacuum%nmzd)
! definitions / initialisations: ! DEFINITIONS / ALLOCATIONS / INITIALISATIONS
do iz = 1, vacuum%nmz do iz = 1, vacuum%nmz
z(iz) = ( iz - 1 ) * vacuum%delz ! z_i starting at 0 z(iz) = ( iz - 1 ) * vacuum%delz
end do end do
do irec2 = 1, stars%ng2 do irec2 = 1, stars%ng2
g_damped(irec2) = sqrt( stars%sk2(irec2) ** 2 + input%preconditioning_param ** 2 ) g_damped(irec2) = sqrt( stars%sk2(irec2) ** 2 + input%preconditioning_param ** 2 )
vcons(irec2) = tpi_const / g_damped(irec2) vcons(irec2) = tpi_const / g_damped(irec2)
do iz = 1, vacuum%nmz do iz = 1, vacuum%nmz
exp_m(iz,irec2) = exp_save( - g_damped(irec2) * z(iz) ) ! redefined in contribution of vacuum charge density part exp_m(iz,irec2) = exp_save( - g_damped(irec2) * z(iz) )
exp_p(iz,irec2) = exp_save( g_damped(irec2) * z(iz) ) exp_p(iz,irec2) = exp_save( g_damped(irec2) * z(iz) )
end do end do
expDhg(irec2) = exp_save( - cell%z1 * g_damped(irec2) ) expDhg(irec2) = exp_save( - cell%z1 * g_damped(irec2) )
expDg(irec2) = expDhg(irec2) ** 2 expDg(irec2) = exp_save( -2 * cell%z1 * g_damped(irec2) )
end do end do
sum_qz = (0.,0.) sum_qz = (0.,0.)
VVxy = (0.,0.) VVxy = (0.,0.)
...@@ -126,14 +137,14 @@ module m_VYukawaFilm ...@@ -126,14 +137,14 @@ module m_VYukawaFilm
do irec2 = 1, stars%ng2 do irec2 = 1, stars%ng2
do ivac = 1, vacuum%nvac do ivac = 1, vacuum%nvac
sign = 3. - 2. * ivac sign = 3. - 2. * ivac
do kz = -stars%mx3, stars%mx3 do iqz = -stars%mx3, stars%mx3
irec3 = stars%ig(stars%kv2(1,irec2),stars%kv2(2,irec2),kz) irec3 = stars%ig(stars%kv2(1,irec2),stars%kv2(2,irec2),iqz)
! use only stars within the g_max sphere -> stars outside the sphere have per definition index ig3n = 0 ! use only stars within the g_max sphere -> stars outside the sphere have per definition index ig3n = 0
if ( irec3 /= 0 ) then if ( irec3 /= 0 ) then
c_ph(kz,irec2) = stars%rgphs(stars%kv2(1,irec2),stars%kv2(2,irec2),kz) c_ph(iqz,irec2) = stars%rgphs(stars%kv2(1,irec2),stars%kv2(2,irec2),iqz)
qz = kz * cell%bmat(3,3) qz = iqz * cell%bmat(3,3)
signIqz = sign * ImagUnit * qz signIqz = sign * ImagUnit * qz
sum_qz(ivac,irec2) = sum_qz(ivac,irec2) + c_ph(kz,irec2) * psq(irec2) * ( exp( signIqz * cell%z1 ) - exp( - signIqz * cell%z1 ) * expDg(irec2) ) / ( signIqz + g_damped(irec2) ) sum_qz(ivac,irec2) = sum_qz(ivac,irec2) + c_ph(iqz,irec2) * psq(irec3) * ( exp( signIqz * cell%z1 ) - exp( - signIqz * cell%z1 ) * expDg(irec2) ) / ( signIqz + g_damped(irec2) )
endif endif
enddo enddo
if( irec2 /= 1 ) then if( irec2 /= 1 ) then
...@@ -149,12 +160,10 @@ module m_VYukawaFilm ...@@ -149,12 +160,10 @@ module m_VYukawaFilm
! case irec2 > 1: ! case irec2 > 1:
do irec2 = 2, stars%ng2 do irec2 = 2, stars%ng2
exp_m(1:vacuum%nmzxy,irec2) = exp_m(1:vacuum%nmzxy,irec2) * expDhg(irec2) ! z_i can now be thought of as starting from D/2 -> expDhg is the correction
exp_p(1:vacuum%nmzxy,irec2) = exp_p(1:vacuum%nmzxy,irec2) * expDhg(irec2)
do ivac = 1, vacuum%nvac do ivac = 1, vacuum%nvac
! integrands: ! integrands:
fa(1:vacuum%nmzxy,irec2) = rhtxy(1:vacuum%nmzxy,irec2-1,ivac) * exp_m(1:vacuum%nmzxy,irec2) fa(1:vacuum%nmzxy,irec2) = rhtxy(1:vacuum%nmzxy,irec2-1,ivac) * exp_m(1:vacuum%nmzxy,irec2) * expDhg(irec2)
fb(1:vacuum%nmzxy,irec2) = rhtxy(1:vacuum%nmzxy,irec2-1,ivac) * exp_p(1:vacuum%nmzxy,irec2) fb(1:vacuum%nmzxy,irec2) = rhtxy(1:vacuum%nmzxy,irec2-1,ivac) * exp_p(1:vacuum%nmzxy,irec2) * expDhg(irec2)
! integrals: ! integrals:
! alpha(z,q_xy,ivac) = int_z^infty rho(z',q_xy,ivac) exp(-sqrt(q_xy**2+prec_param**2)*z') dz' ! alpha(z,q_xy,ivac) = int_z^infty rho(z',q_xy,ivac) exp(-sqrt(q_xy**2+prec_param**2)*z') dz'
! beta (z,q_xy,ivac) = int_{D/2}^z rho(z',q_xy,ivac) exp(+sqrt(q_xy**2+prec_param**2)*z') dz' ! beta (z,q_xy,ivac) = int_{D/2}^z rho(z',q_xy,ivac) exp(+sqrt(q_xy**2+prec_param**2)*z') dz'
...@@ -175,26 +184,24 @@ module m_VYukawaFilm ...@@ -175,26 +184,24 @@ module m_VYukawaFilm
gamma(1:vacuum%nmzxy,irec2) = exp_m(1:vacuum%nmzxy,irec2) * ( alphm(irec2,mod(ivac,2)+1) + beta(1:vacuum%nmzxy,irec2,ivac) ) & gamma(1:vacuum%nmzxy,irec2) = exp_m(1:vacuum%nmzxy,irec2) * ( alphm(irec2,mod(ivac,2)+1) + beta(1:vacuum%nmzxy,irec2,ivac) ) &
+ exp_p(1:vacuum%nmzxy,irec2) * alpha(1:vacuum%nmzxy,irec2,ivac) ! mod(ivac,2)+1 outputs the other ivac value + exp_p(1:vacuum%nmzxy,irec2) * alpha(1:vacuum%nmzxy,irec2,ivac) ! mod(ivac,2)+1 outputs the other ivac value
where ( 2. * gamma(:,irec2) == gamma(:,irec2) ) gamma(:,irec2) = cmplx( 0., 0. ) where ( 2. * gamma(:,irec2) == gamma(:,irec2) ) gamma(:,irec2) = cmplx( 0., 0. )
VVxy(1:vacuum%nmzxy,irec2,ivac) = VVxy(1:vacuum%nmzxy,irec2,ivac) + vcons(irec2) * gamma(1:vacuum%nmzxy,irec2) VVxy(1:vacuum%nmzxy,irec2,ivac) = VVxy(1:vacuum%nmzxy,irec2,ivac) + vcons(irec2) * gamma(1:vacuum%nmzxy,irec2) * expDhg(irec2)
end do end do
end do end do
! case irec2 = 1: ! case irec2 = 1:
exp_m(1:vacuum%nmz,1) = exp_m(1:vacuum%nmz,1) * expDhg(1)
exp_p(1:vacuum%nmz,1) = exp_p(1:vacuum%nmz,1) * expDhg(1)
do ivac = 1, vacuum%nvac do ivac = 1, vacuum%nvac
ga(1:vacuum%nmz) = rht(1:vacuum%nmz,ivac) * exp_m(1:vacuum%nmz,1) ga(1:vacuum%nmz) = rht(1:vacuum%nmz,ivac) * exp_m(1:vacuum%nmz,1) * expDhg(1)
gb(1:vacuum%nmz) = rht(1:vacuum%nmz,ivac) * exp_p(1:vacuum%nmz,1) gb(1:vacuum%nmz) = rht(1:vacuum%nmz,ivac) * exp_p(1:vacuum%nmz,1) * expDhg(1)
call intgz1Reverse( ga(:), vacuum%delz, vacuum%nmz, delta(:,ivac), .true. ) ! integrals call intgz1Reverse( ga(:), vacuum%delz, vacuum%nmz, delta(:,ivac), .true. ) ! integrals
call qsf( vacuum%delz, gb(:), epsilon(:,ivac), vacuum%nmz, 1 ) call qsf( vacuum%delz, gb(:), epsilon(:,ivac), vacuum%nmz, 1 )
alphm(1,ivac) = delta(1,ivac) alphm(1,ivac) = delta(1,ivac)
end do end do
if ( vacuum%nvac == 1 ) alphm(1,2) = alphm(1,1) ! is real, thus no conjg as in the irec2 > 1 case if ( vacuum%nvac == 1 ) alphm(1,2) = alphm(1,1)
do ivac = 1, vacuum%nvac do ivac = 1, vacuum%nvac
zeta(1:vacuum%nmz) = exp_m(1:vacuum%nmz,1) * ( alphm(1,mod(ivac,2)+1) + epsilon(1:vacuum%nmz,ivac) ) & zeta(1:vacuum%nmz) = exp_m(1:vacuum%nmz,1) * ( alphm(1,mod(ivac,2)+1) + epsilon(1:vacuum%nmz,ivac) ) &
+ exp_p(1:vacuum%nmz,1) * delta(1:vacuum%nmz,ivac) + exp_p(1:vacuum%nmz,1) * delta(1:vacuum%nmz,ivac)
where ( 2. * zeta == zeta ) zeta = 0. where ( 2. * zeta == zeta ) zeta = 0.
VVz(1:vacuum%nmz,ivac) = VVz(1:vacuum%nmz,ivac) + vcons(1) * zeta(1:vacuum%nmz) VVz(1:vacuum%nmz,ivac) = VVz(1:vacuum%nmz,ivac) + vcons(1) * zeta(1:vacuum%nmz) * expDhg(1)
end do end do
...@@ -206,9 +213,15 @@ module m_VYukawaFilm ...@@ -206,9 +213,15 @@ module m_VYukawaFilm
psq, VVxy, VVz, alphm, & psq, VVxy, VVz, alphm, &
VIq ) VIq )
! 1. part: compute the contribution from the interstitial charge density to the interstitial potential (largest part) as a function of q_xy and z ! main parts:
! 2. part: add the contribution from the vacuum charge density to the interstitial potential, which had already been computed earlier for the vacuum potential ! 1. part: Compute the contribution from the interstitial charge density to the interstitial potential as a function of q_xy and z (analytic expression for integral)
! 3. part: compute the coefficients V^I(q_xy,q_z) from the function V^I(q_xy,z) by a 1D Fourier transform: V^I(q_xy,z) = sum_{q_z} V^I(q_xy,q_z) * exp( ImagUnit * q_z * z ) ! 2. part: Add the contribution from the vacuum charge density to the interstitial potential, which had already been computed earlier for the vacuum potential
! 4. part: Compute the coefficients V^I(q_xy,q_z) from the function V^I(q_xy,z) by a 1D Fourier transform:
! V^I(q_xy,z) = sum_{q_z} V^I(q_xy,q_z) * exp( ImagUnit * q_z * z )
! In order to be able to match the interstitial and vacuum potentials smoothly at the interface, the Fourier transform is done
! in a slightly larger region. -> 3. part
! 3. part: Interpolate the vacuum potential in a small region surrounding the slab
use m_ExpSave use m_ExpSave
use m_constants use m_constants
...@@ -228,109 +241,136 @@ module m_VYukawaFilm ...@@ -228,109 +241,136 @@ module m_VYukawaFilm
complex, intent(out) :: VIq(stars%ng3) complex, intent(out) :: VIq(stars%ng3)
real :: partitioning, rz, qz, fit, q real :: partitioning, rz, qz, q
integer :: irec2, irec3, iz, jz, ivac, kz, ifft integer :: irec2, irec3, iz, jz, ivac, iqz, nfft, nzmax, nzmin, nzdh, nLower, nUpper, jvac
complex :: VIz(3*stars%mx3,stars%ng2), sum_qz(3*stars%mx3,stars%ng2), eta(3*stars%mx3,stars%ng2) complex, allocatable :: VIz(:,:), eta(:,:)
complex :: VIqz(-stars%mx3:stars%mx3,stars%ng2), c_ph(-stars%mx3:stars%mx3,stars%ng2) complex :: VIqz(-stars%mx3:stars%mx3,stars%ng2), c_ph(-stars%mx3:stars%mx3,stars%ng2)
complex :: vcons1(stars%ng3) complex :: vcons1(stars%ng3)
real :: VIzReal(3*stars%mx3,stars%ng2), VIzImag(3*stars%mx3,stars%ng2), exp_m(3*stars%mx3,stars%ng2), exp_p(3*stars%mx3,stars%ng2) real :: VIzReal(3*stars%mx3,stars%ng2), VIzImag(3*stars%mx3,stars%ng2)
real :: z(3*stars%mx3) real, allocatable :: exp_m(:,:), exp_p(:,:)
real :: g_damped(stars%ng2), vcons2(stars%ng2) real, allocatable :: z(:)
real :: g_damped(stars%ng2), vcons2(stars%ng2), expDhg(stars%ng2)
! definitions / initialisations
ifft = 3 * stars%mx3 ! DEFINITIONS / ALLOCATIONS / INITIALISATIONS
partitioning = 1. / real( ifft )
do iz = 1, ifft ! grid points z_i:
z(iz) = cell%amat(3,3) * ( iz - 1 ) * partitioning ! z_i are equidistantly distributed along ( 0, cell%amat(3,3) ] nfft = 3 * stars%mx3 ! number of grid points for Fourier transform
if( z(iz) > cell%amat(3,3) / 2. ) z(iz) = z(iz) - cell%amat(3,3) ! z_i are equidistantly distributed along ( -cell%amat(3,3) / 2, cell%amat(3,3) / 2 ] partitioning = 1. / real( nfft )
nzmax = nfft / 2 ! index of maximal z below D~/2
nzmin = nzmax - nfft + 1 ! index of minimal z above -D~/2
nzdh = ceiling( cell%z1 / cell%amat(3,3) * nfft ) - 1 ! index of maximal z below D/2
allocate( z(nzmin:nzmax) )
! definition of z_i: ! indexing: z_0 = 0; positive indices for z > 0; negative indices for z < 0
do iz = nzmin, nzmax
z(iz) = cell%amat(3,3) * iz * partitioning
end do end do
! other variables:
allocate( VIz(nzmin:nzmax,stars%ng2), eta(nzmin:nzmax,stars%ng2) )
allocate( exp_m(nzmin:nzmax,stars%ng2), exp_p(nzmin:nzmax,stars%ng2) )
do irec2 = 1, stars%ng2 do irec2 = 1, stars%ng2
g_damped(irec2) = sqrt( stars%sk2(irec2) ** 2 + input%preconditioning_param ** 2 ) g_damped(irec2) = sqrt( stars%sk2(irec2) ** 2 + input%preconditioning_param ** 2 )
vcons2(irec2) = -1. / ( 2. * g_damped(irec2) ) vcons2(irec2) = -1. / ( 2. * g_damped(irec2) )
do iz = 1, ifft do iz = nzmin, nzmax
exp_m(iz,irec2) = exp_save( - g_damped(irec2) * ( cell%z1 + z(iz) ) ) exp_m(iz,irec2) = exp_save( - g_damped(irec2) * z(iz) )
exp_p(iz,irec2) = exp_save( - g_damped(irec2) * ( cell%z1 - z(iz) ) ) exp_p(iz,irec2) = exp_save( g_damped(irec2) * z(iz) )
end do end do
expDhg(irec2) = exp_save( - cell%z1 * g_damped(irec2) )
end do end do
do irec3 = 1, stars%ng3 do irec3 = 1, stars%ng3
vcons1(irec3) = fpi_const * psq(irec3) / ( stars%sk3(irec3) ** 2 + input%preconditioning_param ** 2 ) vcons1(irec3) = fpi_const * psq(irec3) / ( stars%sk3(irec3) ** 2 + input%preconditioning_param ** 2 )
end do end do
sum_qz = (0.,0.)
VIz = (0.,0.) VIz = (0.,0.)
VIq = (0.,0.) VIq = (0.,0.)
! CONTRIBUTION FROM THE INTERSTITIAL CHARGE DENSITY ! CONTRIBUTION FROM THE INTERSTITIAL CHARGE DENSITY
! compute V^I(q_xy,z) as a function of q_xy and z ! compute V^I(q_xy,z) as a function of q_xy and z
do irec2 = 1, stars%ng2 do irec2 = 1, stars%ng2
do iz = 1, ifft do iz = -nzdh, nzdh
do iqz = -stars%mx3, stars%mx3
! in the transition region ( D/2, D~/2 ), smooth out the potential numerically irec3 = stars%ig(stars%kv2(1,irec2),stars%kv2(2,irec2),iqz)
if ( abs( z(iz) ) >= cell%z1 ) then if ( irec3 /= 0 ) then ! use only stars within the g_max sphere
ivac = 1 ! upper vacuum c_ph(iqz,irec2) = stars%rgphs(stars%kv2(1,irec2),stars%kv2(2,irec2),iqz)
if ( z(iz) < 0.0 .and. .not. ( sym%invs .or. sym%zrfs ) ) ivac = 2 ! lower vacuum qz = iqz * cell%bmat(3,3)
rz = ( abs( z(iz) ) - cell%z1 ) / vacuum%delz + 1.0 ! numbering the grid points in the vacuum region VIz(iz,irec2) = VIz(iz,irec2) &
jz = rz + vcons1(irec3) * c_ph(iqz,irec2) * &
q = rz - jz ( exp( ImagUnit * qz * z(iz) ) &
if ( irec2 == 1 ) then + vcons2(irec2) * expDhg(irec2) * &
fit = 0.5 * ( q - 1. ) * ( q - 2. ) * VVz(jz, ivac) & ( ( g_damped(irec2) + ImagUnit * qz ) * exp_p(iz,irec2) * exp( ImagUnit * qz * cell%z1 ) &
- q * ( q - 2. ) * VVz(jz+1,ivac) & + ( g_damped(irec2) - ImagUnit * qz ) * exp_m(iz,irec2) * exp( - ImagUnit * qz * cell%z1 ) ) )
+ 0.5 * q * ( q - 1. ) * VVz(jz+2,ivac)
VIz(iz,irec2) = cmplx( fit, 0.0 )
else if ( jz + 2 <= vacuum%nmzxy ) then
VIz(iz,irec2) = 0.5 * ( q - 1. ) * ( q - 2. ) * VVxy(jz, irec2-1,ivac) &
- q * ( q - 2. ) * VVxy(jz+1,irec2-1,ivac) &
+ 0.5 * q * ( q - 1. ) * VVxy(jz+2,irec2-1,ivac)
if ( ( sym%invs .and. .not. sym%zrfs ) .and. z(iz) < 0 ) VIz(iz,irec2) = conjg( VIz(iz,irec2) )
end if end if
enddo
! z in ( -D/2, D/2 )
else
do kz = -stars%mx3, stars%mx3
irec3 = stars%ig(stars%kv2(1,irec2),stars%kv2(2,irec2),kz)
if ( irec3 /= 0 ) then ! use only stars within the g_max sphere
c_ph(kz,irec2) = stars%rgphs(stars%kv2(1,irec2),stars%kv2(2,irec2),kz)
qz = kz * cell%bmat(3,3)
sum_qz(iz,irec2) = sum_qz(iz,irec2) + exp( ImagUnit * qz * z(iz) ) &
+ vcons2(irec2) * ( ( g_damped(irec2) + ImagUnit * qz ) * exp_p(iz,irec2) * exp( ImagUnit * qz * cell%z1 ) &
+ ( g_damped(irec2) - ImagUnit * qz ) * exp_m(iz,irec2) * exp( - ImagUnit * qz * cell%z1 ) )
sum_qz(iz,irec2) = vcons1(irec3) * c_ph(kz,irec2) * sum_qz(iz,irec2)
end if
enddo
VIz(iz,irec2) = VIz(iz,irec2) + sum_qz(iz,irec2)
end if
end do end do
! irec2 loop continues
! CONTRIBUTION FROM THE VACUUM CHARGE DENSITY ! CONTRIBUTION FROM THE VACUUM CHARGE DENSITY
eta(:,irec2) = exp_m(:,irec2) * alphm(irec2,2) + exp_p(:,irec2) * alphm(irec2,1) ! irec2 loop continues
eta(-nzdh:nzdh,irec2) = exp_m(-nzdh:nzdh,irec2) * alphm(irec2,2) + exp_p(-nzdh:nzdh,irec2) * alphm(irec2,1)
where ( 2.0 * eta(:,irec2) == eta(:,irec2) ) eta(:,irec2) = cmplx( 0.0, 0.0 ) where ( 2.0 * eta(:,irec2) == eta(:,irec2) ) eta(:,irec2) = cmplx( 0.0, 0.0 )
VIz(:,irec2) = VIz(:,irec2) + tpi_const / g_damped(irec2) * eta(:,irec2) VIz(-nzdh:nzdh,irec2) = VIz(-nzdh:nzdh,irec2) + tpi_const / g_damped(irec2) * eta(-nzdh:nzdh,irec2)
! irec2 loop continues
! INTERPOLATION IN VACUUM REGION
! use Lagrange polynomials of order 3 to interpolate the vacuum potential outside I
! q, q-1 and q-2 (scaled with +/-1 or +/-0.5) are the factors of the Lagrange basis polynomials
! irec2 loop continues
do ivac = 1, 2
select case( ivac )
case( 1 )
nUpper = nzmax; nLower = nzdh + 1; jvac = 1
case( 2 )
nUpper = nzmin; nLower = -nzdh - 1; jvac = 2; if ( sym%invs .or. sym%zrfs ) jvac = 1
end select
do iz = nLower, nUpper
rz = ( abs( z(iz) ) - cell%z1 ) / vacuum%delz + 1.0
jz = rz ! index of maximal vacuum grid point below z_i
q = rz - jz ! factor in Lagrange basis polynomials
if ( irec2 == 1 ) then
VIz(iz,irec2) = 0.5 * ( q - 1. ) * ( q - 2. ) * VVz(jz, jvac) &
- q * ( q - 2. ) * VVz(jz+1,jvac) &
+ 0.5 * q * ( q - 1. ) * VVz(jz+2,jvac)
else if ( jz + 2 <= vacuum%nmzxy ) then
VIz(iz,irec2) = 0.5 * ( q - 1. ) * ( q - 2. ) * VVxy(jz, irec2,jvac) &
- q * ( q - 2. ) * VVxy(jz+1,irec2,jvac) &
+ 0.5 * q * ( q - 1. ) * VVxy(jz+2,irec2,jvac)
if ( ( sym%invs .and. .not. sym%zrfs ) .and. ivac == 2 ) VIz(iz,irec2) = conjg( VIz(iz,irec2) )
end if
end do
end do
end do ! irec2
! 1D FOURIER TRANSFORM TO FIND THE COEFFICIENTS V^I(q_xy,q_z)
! 1D FOURIER TRANSFORM TO FIND THE COEFFICIENTS V^I(q_xy,q_z) ! change the indexing for the subroutine cfft, and split real and imaginary parts
VIzReal(1:nzmax+1,:) = real( VIz(0:nzmax,:) ); VIzReal(nzmax+2:nfft,:) = real( VIz(nzmin:-1,:) )
VIzImag(1:nzmax+1,:) = aimag( VIz(0:nzmax,:) ); VIzImag(nzmax+2:nfft,:) = aimag( VIz(nzmin:-1,:) )
! V^I(q_xy,z) = sum_{q_z} V^I(q_xy,q_z) * exp( ImagUnit * q_z * z )
do irec2 = 1, stars%ng2
call cfft( VIzReal(:,irec2), VIzImag(:,irec2), nfft, nfft, nfft, -1 )
! irec2 loop continues
! V^I(q_xy,z) = sum_{q_z} V^I(q_xy,q_z) * exp( ImagUnit * q_z * z ) ! reorder
VIzReal = real( VIz )
VIzImag = aimag( VIz )
call cfft( VIzReal(:,irec2), VIzImag(:,irec2), ifft, ifft, ifft, -1 )
! reorder:
VIqz(0,irec2) = cmplx( VIzReal(1,irec2), VIzImag(1,irec2) ) VIqz(0,irec2) = cmplx( VIzReal(1,irec2), VIzImag(1,irec2) )
do kz = 1, stars%mx3 do iqz = 1, stars%mx3
VIqz( kz,irec2) = cmplx( VIzReal(kz+1, irec2), VIzImag(kz+1, irec2) ) VIqz( iqz,irec2) = cmplx( VIzReal(iqz+1, irec2), VIzImag(iqz+1, irec2) )
VIqz(-kz,irec2) = cmplx( VIzReal(ifft+1-kz,irec2), VIzImag(ifft+1-kz,irec2) ) VIqz(-iqz,irec2) = cmplx( VIzReal(nfft+1-iqz,irec2), VIzImag(nfft+1-iqz,irec2) )
end do