Commit 0bf2eadd authored by Miriam Hinzen's avatar Miriam Hinzen

Film-Case Kerker Preconditions in Film Region Only

parent 400cfea6
......@@ -28,7 +28,6 @@ CONTAINS
REAL :: q2(vacuum%nmz), w, rht1(vacuum%nmzd,2,input%jspins)
COMPLEX :: x(stars%ng3)
CALL timestart("cdntot")
qtot = 0.0
qistot = 0.0
DO jsp = 1,input%jspins
......@@ -109,7 +108,10 @@ CONTAINS
REAL qmt(atoms%ntype,input%jspins),qvac(2,input%jspins)
INTEGER, ALLOCATABLE :: lengths(:,:)
CHARACTER(LEN=20) :: attributes(6), names(6)
CALL timestart("cdntot")
call cdntot_integrate(stars,atoms,sym,vacuum,input,cell,oneD, den, &
q, qis, qmt, qvac, qtot, qistot)
......
......@@ -87,9 +87,9 @@ contains
! KERKER PRECONDITIONER
IF( input%preconditioning_param /= 0 ) THEN
CALL timestart("Preconditioner")
CALL kerker(field, DIMENSION, mpi, &
stars, atoms, sphhar, vacuum, input, sym, cell, noco, &
oneD, inDen, outDen, fsm(it) )
CALL kerker( field, DIMENSION, mpi, &
stars, atoms, sphhar, vacuum, input, sym, cell, noco, &
oneD, inDen, outDen, fsm(it) )
!Store modified density in history
CALL mixing_history_store(fsm(it))
CALL timestop("Preconditioner")
......
......@@ -4,12 +4,15 @@
! of the MIT license as expressed in the LICENSE file in more detail.
!--------------------------------------------------------------------------------
MODULE m_kerker
CONTAINS
SUBROUTINE kerker(field, DIMENSION, mpi, &
SUBROUTINE kerker( field, DIMENSION, mpi, &
stars, atoms, sphhar, vacuum, input, sym, cell, noco, &
oneD, inDen, outDen, precon_v )
!Implementation of the Kerker preconditioner by M.Hinzen
USE m_vgen_coulomb
USE m_VYukawaFilm
USE m_juDFT
......@@ -19,6 +22,7 @@ CONTAINS
USE m_constants
USE m_mpi_bc_potden
IMPLICIT NONE
TYPE(t_oneD), INTENT(in) :: oneD
TYPE(t_input), INTENT(in) :: input
TYPE(t_vacuum), INTENT(in) :: vacuum
......@@ -35,10 +39,9 @@ CONTAINS
TYPE(t_potden), INTENT(in) :: inDen
TYPE(t_mixvector), INTENT(INOUT) :: precon_v
!Locals
type(t_potden) :: resDen, vYukawa
real :: fix
integer :: lh,n
type(t_potden) :: resDen, vYukawa, resDenMod
real :: fix
integer :: lh,n
CALL resDen%init( stars, atoms, sphhar, vacuum, noco, input%jspins, POTDEN_TYPE_DEN )
CALL vYukawa%init( stars, atoms, sphhar, vacuum, noco, input%jspins, 4 )
......@@ -53,8 +56,12 @@ CONTAINS
CALL vgen_coulomb( 1, mpi, DIMENSION, oneD, input, field, vacuum, sym, stars, cell, &
sphhar, atoms, resDen, vYukawa )
ELSE
if( mpi%irank == 0 ) then
call resDenMod%init( stars, atoms, sphhar, vacuum, noco, input%jspins, POTDEN_TYPE_DEN )
call resDenMod%copyPotDen( resDen )
end if
vYukawa%iter = resDen%iter
CALL VYukawaFilm( stars, vacuum, cell, sym, input, mpi, atoms, sphhar, DIMENSION, oneD, resDen, &
CALL VYukawaFilm( stars, vacuum, cell, sym, input, mpi, atoms, sphhar, DIMENSION, oneD, noco, resDenMod, &
vYukawa )
END IF
......@@ -67,8 +74,6 @@ CONTAINS
* vYukawa%mt(1:atoms%jri(n),lh,n,1) * atoms%rmsh(1:atoms%jri(n),n) ** 2
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()
! fix the preconditioned density
CALL outDen%addPotDen( resDen, inDen )
......@@ -76,6 +81,7 @@ CONTAINS
CALL resDen%subPotDen( outDen, inDen )
END IF MPI0_c
CALL precon_v%from_density(resden)
! end of preconditioner
END SUBROUTINE kerker
END SUBROUTINE kerker
END MODULE m_kerker
......@@ -42,6 +42,7 @@ MODULE m_types_potden
procedure :: ChargeAndMagnetisationToSpins
procedure :: addPotDen
procedure :: subPotDen
procedure :: copyPotDen
procedure :: distribute
procedure :: collect
END TYPE t_potden
......@@ -233,6 +234,21 @@ CONTAINS
end subroutine
subroutine copyPotDen( PotDenCopy, PotDen )
implicit none
class(t_potden), intent(in) :: PotDen
class(t_potden), intent(inout) :: PotDenCopy
PotDenCopy%iter = PotDen%iter
PotDenCopy%potdenType = PotDen%potdenType
PotDenCopy%mt = PotDen%mt
PotDenCopy%pw = PotDen%pw
PotDenCopy%vacz = PotDen%vacz
PotDenCopy%vacxy = PotDen%vacxy
end subroutine copyPotDen
SUBROUTINE init_potden_types(pd,stars,atoms,sphhar,vacuum,noco,jspins,potden_type)
USE m_judft
USE m_types_setup
......
module m_VYukawaFilm
! 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
subroutine VYukawaFilm( stars, vacuum, cell, sym, input, mpi, atoms, sphhar, dimension, oneD, den, &
! In a first step we solve the modified Helmholtz equation with right-hand
! side -4 pi ( rho_out - rho_in ) subject to the boundary condition that the
! potential equals zero on the surface of the film.
! In the interstitial we use a Green function method to solve the
! differential equation (separately for each q_xy) for a z-dependent function.
! The muffin-tin potential generation is the same as for bulk materials.
! Unfortunately, the potential cannot satisfy the condition that its
! integral over the film is zero at the same time, which is needed for
! charge neutrality of the preconditioned density.
! Therefore we add a solution to the differential equation with constant
! right-hand side sufficing the same boundary condition, where the constant
! is chosen such that the integral over the sum of these potentials is zero.
! So we basically added a constant to the charge density in order to make
! the resulting Yukawa potential charge neutral.
contains
subroutine VYukawaFilm( stars, vacuum, cell, sym, input, mpi, atoms, sphhar, dimension, oneD, noco, den, &
VYukawa )
use m_constants
......@@ -35,199 +37,51 @@ module m_VYukawaFilm
type(t_sphhar), intent(in) :: sphhar
type(t_dimension), intent(in) :: dimension
type(t_oneD), intent(in) :: oneD
type(t_potden), intent(in) :: den
type(t_noco), intent(in) :: noco
type(t_potden), intent(inout) :: den
type(t_potden), intent(inout) :: VYukawa
complex :: psq(stars%ng3)
complex :: alphm(stars%ng2,2)
! PSEUDO-CHARGE DENSITY
call psqpw( mpi, atoms, sphhar, stars, vacuum, dimension, cell, input, sym, oneD, den%pw(:,1), den%mt(:,:,:,1), den%vacz(:,:,1), .false., VYukawa%potdenType, psq )
! VACUUM POTENTIAL
call VYukawaFilmVacuum( stars, vacuum, cell, sym, input, &
psq, den%vacxy(:,:,:,1), den%vacz(:,:,1), &
VYukawa%vacxy, VYukawa%vacz, alphm )
call psqpw( mpi, atoms, sphhar, stars, vacuum, dimension, cell, input, sym, oneD, &
den%pw(:,1), den%mt(:,:,:,1), den%vacz(:,:,1), .false., VYukawa%potdenType, &
psq )
! INTERSTITIAL POTENTIAL
call VYukawaFilmInterstitial( stars, vacuum, cell, sym, input, &
psq, VYukawa%vacxy, VYukawa%vacz, alphm, &
psq, &
VYukawa%pw(:,1) )
! MUFFIN-TIN POTENTIAL
call Vmts( input, mpi, stars, sphhar, atoms, sym, cell, oneD, VYukawa%pw(:,1), den%mt(:,0:,:,1), VYukawa%potdenType, VYukawa%mt(:,0:,:,1) )
end subroutine VYukawaFilm
subroutine VYukawaFilmVacuum( stars, vacuum, cell, sym, input, &
psq, rhtxy, rht, &
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 (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
use m_ExpSave
use m_constants
use m_types
use m_intgr, only: intgz1Reverse
use m_qsf
implicit none
call Vmts( input, mpi, stars, sphhar, atoms, sym, cell, oneD, &
VYukawa%pw(:,1), den%mt(:,0:,:,1), VYukawa%potdenType, &
VYukawa%mt(:,0:,:,1) )
type(t_stars), intent(in) :: stars
type(t_vacuum), intent(in) :: vacuum
type(t_cell), intent(in) :: cell
type(t_sym), intent(in) :: sym
type(t_input), intent(in) :: input
complex, intent(in) :: psq(stars%ng3)
complex, intent(in) :: rhtxy(vacuum%nmzxyd,stars%ng2-1,2)
real, intent(in) :: rht(vacuum%nmzd,2)
complex, intent(out) :: VVxy(vacuum%nmzxyd,2:stars%ng2,2) ! this is the qxy /= 0 part of the vacuum potential
real, intent(out) :: VVz(vacuum%nmzd,2) ! this is the qxy = 0 part of the vacuum potential
complex, intent(out) :: alphm(stars%ng2,2) ! these are the integrals in upper and lower vacuum, now including the first star---integral for star ig2 is in alphm(ig2,ivac)
complex :: sum_qz(2,stars%ng2)
complex :: c_ph(-stars%mx3:stars%mx3,stars%ng2)
complex :: signIqz
real :: g_damped(stars%ng2), qz, sign, vcons(stars%ng2)
real :: exp_m(vacuum%nmzd,stars%ng2), exp_p(vacuum%nmzd,stars%ng2)
real :: expDhg(stars%ng2), expDg(stars%ng2)
real :: z(vacuum%nmzd)
integer :: iz, irec2, irec3, ivac, iqz
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)
real :: ga(vacuum%nmzd), gb(vacuum%nmzd)
real :: delta(vacuum%nmzd,2), epsilon(vacuum%nmzd,2), zeta(vacuum%nmzd)
! DEFINITIONS / ALLOCATIONS / INITIALISATIONS
do iz = 1, vacuum%nmz
z(iz) = ( iz - 1 ) * vacuum%delz
end do
do irec2 = 1, stars%ng2
g_damped(irec2) = sqrt( stars%sk2(irec2) ** 2 + input%preconditioning_param ** 2 )
vcons(irec2) = tpi_const / g_damped(irec2)
do iz = 1, vacuum%nmz
exp_m(iz,irec2) = exp_save( - g_damped(irec2) * z(iz) )
exp_p(iz,irec2) = exp_save( g_damped(irec2) * z(iz) )
end do
expDhg(irec2) = exp_save( - cell%z1 * g_damped(irec2) )
expDg(irec2) = exp_save( -2 * cell%z1 * g_damped(irec2) )
end do
sum_qz = (0.,0.)
VVxy = (0.,0.)
VVz = 0.
! CONTRIBUTION FROM THE INTERSTITIAL CHARGE DENSITY
do irec2 = 1, stars%ng2
do ivac = 1, vacuum%nvac
sign = 3. - 2. * ivac
do iqz = -stars%mx3, stars%mx3
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
if ( irec3 /= 0 ) then
c_ph(iqz,irec2) = stars%rgphs(stars%kv2(1,irec2),stars%kv2(2,irec2),iqz)
qz = iqz * cell%bmat(3,3)
signIqz = sign * ImagUnit * qz
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
enddo
if( irec2 /= 1 ) then
VVxy(1:vacuum%nmzxy,irec2,ivac) = vcons(irec2) * sum_qz(ivac,irec2) * exp_m(1:vacuum%nmzxy,irec2)
else
VVz(1:vacuum%nmz,ivac) = vcons(1) * sum_qz(ivac,1) * exp_m(1:vacuum%nmz,1)
end if
enddo
enddo
! CONTRIBUTION FROM THE VACUUM CHARGE DENSITY
! shifting z:
do irec2 = 1, stars%ng2
exp_m(1:vacuum%nmz,irec2) = exp_m(1:vacuum%nmz,irec2) * expDhg(irec2)
exp_p(1:vacuum%nmz,irec2) = exp_p(1:vacuum%nmz,irec2) / expDhg(irec2)
end do
! case irec2 > 1:
do irec2 = 2, stars%ng2
do ivac = 1, vacuum%nvac
! integrands:
fa(1:vacuum%nmzxy,irec2) = rhtxy(1:vacuum%nmzxy,irec2-1,ivac) * exp_m(1:vacuum%nmzxy,irec2)
fb(1:vacuum%nmzxy,irec2) = rhtxy(1:vacuum%nmzxy,irec2-1,ivac) * exp_p(1:vacuum%nmzxy,irec2)
! integrals:
! 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'
! where for z < 0 the lower vacuum charge density (ivac=2) is defined by rho(q_xy,z,ivac=2) := rho(q_xy,-z,ivac=2)
call intgz1Reverse( fa(:,irec2), vacuum%delz, vacuum%nmzxy, alpha(:,irec2,ivac), .true. )
call qsfComplex( vacuum%delz, fb(:,irec2), beta(:,irec2,ivac), vacuum%nmzxy, 1 )
! alphm(q_xy,ivac) = alpha(D/2,q_xy,ivac) --- these integrals are also needed for the interstitial potential
alphm(irec2,ivac) = alpha(1,irec2,ivac)
end do
if ( vacuum%nvac == 1 ) then
if ( sym%invs ) then
alphm(irec2,2) = conjg( alphm(irec2,1) )
else
alphm(irec2,2) = alphm(irec2,1)
end if
end if
do ivac = 1, vacuum%nvac
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
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)
end do
end do
! MODIFICATION FOR CHARGE NEUTRALITY
! case irec2 = 1:
do ivac = 1, vacuum%nvac
ga(1:vacuum%nmz) = rht(1:vacuum%nmz,ivac) * exp_m(1:vacuum%nmz,1)
gb(1:vacuum%nmz) = rht(1:vacuum%nmz,ivac) * exp_p(1:vacuum%nmz,1)
call intgz1Reverse( ga(:), vacuum%delz, vacuum%nmz, delta(:,ivac), .true. ) ! integrals
call qsf( vacuum%delz, gb(:), epsilon(:,ivac), vacuum%nmz, 1 )
alphm(1,ivac) = delta(1,ivac)
end do
if ( vacuum%nvac == 1 ) alphm(1,2) = alphm(1,1)
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) ) &
+ exp_p(1:vacuum%nmz,1) * delta(1:vacuum%nmz,ivac)
where ( 2. * zeta == zeta ) zeta = 0.
VVz(1:vacuum%nmz,ivac) = VVz(1:vacuum%nmz,ivac) + vcons(1) * zeta(1:vacuum%nmz)
end do
call VYukawaModify( stars, vacuum, cell, sym, input, mpi, atoms, sphhar, dimension, oneD, noco, den, &
VYukawa )
end subroutine VYukawaFilmVacuum
end subroutine VYukawaFilm
subroutine VYukawaFilmInterstitial( stars, vacuum, cell, sym, input, &
psq, VVxy, VVz, alphm, &
psq, &
VIq )
! main parts:
! 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)
! 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
! This subroutine computes the basic interstitial Yukawa potential that is
! zero on the film surface. Its integral is not zero.
use m_ExpSave
use m_constants
......@@ -235,53 +89,62 @@ module m_VYukawaFilm
use m_cfft
implicit none
type(t_stars), intent(in) :: stars
type(t_vacuum), intent(in) :: vacuum
type(t_cell), intent(in) :: cell
type(t_sym), intent(in) :: sym
type(t_input), intent(in) :: input
complex, intent(in) :: psq(stars%ng3)
complex, intent(in) :: VVxy(vacuum%nmzxyd,2:stars%ng2,2)
real, intent(in) :: VVz(vacuum%nmzd,2)
complex, intent(in) :: alphm(stars%ng2,2)
complex, intent(out) :: VIq(stars%ng3)
real :: partitioning, rz, qz, q
integer :: irec2, irec3, iz, jz, ivac, iqz, nfft, nzmax, nzmin, nzdh, nLower, nUpper, jvac
complex, allocatable :: VIz(:,:), eta(:,:)
complex :: VIqz(-stars%mx3:stars%mx3,stars%ng2), c_ph(-stars%mx3:stars%mx3,stars%ng2)
complex :: vcons1(stars%ng3)
real :: VIzReal(3*stars%mx3,stars%ng2), VIzImag(3*stars%mx3,stars%ng2)
real, allocatable :: exp_m(:,:), exp_p(:,:)
real, allocatable :: z(:)
real :: g_damped(stars%ng2), vcons2(stars%ng2), expDhg(stars%ng2)
type(t_stars), intent(in) :: stars
type(t_vacuum), intent(in) :: vacuum
type(t_cell), intent(in) :: cell
type(t_sym), intent(in) :: sym
type(t_input), intent(in) :: input
complex, intent(inout) :: psq(stars%ng3)
complex, intent(out) :: VIq(stars%ng3)
real :: partitioning, rz, qz, q, q0, dh, qxy_numerics
integer :: irec2, irec3, iz, ivac, iqz, nfft, nzmax, nzmin, nzdh
complex, allocatable :: VIz(:,:)
complex :: VIqz(-stars%mx3:stars%mx3,stars%ng2), c_ph(-stars%mx3:stars%mx3,stars%ng2)
complex :: vcons1(stars%ng3)
real :: VIzReal(3*stars%mx3,stars%ng2), VIzImag(3*stars%mx3,stars%ng2)
real, allocatable :: cosh_quot(:,:), sinh_quot(:,:), exp_m(:,:), exp_p(:,:)
real, allocatable :: z(:)
real :: g_damped(stars%ng2), cosh_dh(stars%ng2), sinh_dh(stars%ng2)
! DEFINITIONS / ALLOCATIONS / INITIALISATIONS
! grid points z_i:
nfft = 3 * stars%mx3 ! number of grid points for Fourier transform
! constants:
dh = cell%z1
qxy_numerics = sqrt( ( 100 / dh ) ** 2 - input%preconditioning_param ** 2 )
! indexing of the grid points z_i:
nfft = 3 * stars%mx3 ! number of grid points for Fourier transform
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
nzmax = nfft / 2 ! index of maximal z below D~/2
nzmin = nzmax - nfft + 1 ! index of minimal z above -D~/2
nzdh = ceiling( dh / 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
! definition of the grid points 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
! 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) )
allocate( VIz(nzmin:nzmax,stars%ng2) )
allocate( cosh_quot(-nzdh+1:nzdh-1,stars%ng2), sinh_quot(-nzdh+1:nzdh-1,stars%ng2) )
allocate( exp_m(-nzdh+1:nzdh-1,stars%ng2), exp_p(-nzdh+1:nzdh-1,stars%ng2) )
do irec2 = 1, stars%ng2
g_damped(irec2) = sqrt( stars%sk2(irec2) ** 2 + input%preconditioning_param ** 2 )
vcons2(irec2) = -1. / ( 2. * g_damped(irec2) )
do iz = nzmin, nzmax
exp_m(iz,irec2) = exp_save( - g_damped(irec2) * z(iz) )
exp_p(iz,irec2) = exp_save( g_damped(irec2) * z(iz) )
end do
expDhg(irec2) = exp_save( - cell%z1 * g_damped(irec2) )
if( stars%sk2(irec2) < qxy_numerics ) then ! numerics ok
cosh_dh(irec2) = cosh( g_damped(irec2) * dh )
sinh_dh(irec2) = sinh( g_damped(irec2) * dh )
do iz = -nzdh+1, nzdh-1
cosh_quot(iz,irec2) = cosh( g_damped(irec2) * z(iz) ) / cosh_dh(irec2)
sinh_quot(iz,irec2) = sinh( g_damped(irec2) * z(iz) ) / sinh_dh(irec2)
end do
else ! numerical treatment necessary
do iz = -nzdh+1, nzdh-1
exp_m(iz,irec2) = exp_save( - g_damped(irec2) * ( dh - z(iz) ) )
exp_p(iz,irec2) = exp_save( - g_damped(irec2) * ( dh + z(iz) ) )
end do
end if
end do
do irec3 = 1, stars%ng3
vcons1(irec3) = fpi_const * psq(irec3) / ( stars%sk3(irec3) ** 2 + input%preconditioning_param ** 2 )
......@@ -291,65 +154,40 @@ module m_VYukawaFilm
! CONTRIBUTION FROM THE INTERSTITIAL CHARGE DENSITY
! compute V^I(q_xy,z) as a function of q_xy and z
do irec2 = 1, stars%ng2
do iz = -nzdh, nzdh
do iqz = -stars%mx3, stars%mx3
irec3 = stars%ig(stars%kv2(1,irec2),stars%kv2(2,irec2),iqz)
if ( irec3 /= 0 ) then ! use only stars within the g_max sphere
c_ph(iqz,irec2) = stars%rgphs(stars%kv2(1,irec2),stars%kv2(2,irec2),iqz)
qz = iqz * cell%bmat(3,3)
VIz(iz,irec2) = VIz(iz,irec2) &
+ vcons1(irec3) * c_ph(iqz,irec2) * &
( exp( ImagUnit * qz * z(iz) ) &
+ vcons2(irec2) * expDhg(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 ) ) )
end if
enddo
end do
! irec2 loop continues
! CONTRIBUTION FROM THE VACUUM CHARGE DENSITY
! 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 )
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 )
nLower = nzmin; nUpper = -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
if( stars%sk2(irec2) < qxy_numerics ) then ! numerics ok
do iz = -nzdh+1, nzdh-1 ! for z=+-D/2 this is zero anyways
do iqz = -stars%mx3, stars%mx3
irec3 = stars%ig(stars%kv2(1,irec2),stars%kv2(2,irec2),iqz)
if ( irec3 /= 0 ) then ! use only stars within the g_max sphere
c_ph(iqz,irec2) = stars%rgphs(stars%kv2(1,irec2),stars%kv2(2,irec2),iqz)
qz = iqz * cell%bmat(3,3)
VIz(iz,irec2) = VIz(iz,irec2) &
+ vcons1(irec3) * c_ph(iqz,irec2) * &
( exp( ImagUnit * qz * z(iz) ) &
- cosh_quot(iz,irec2) * cosh( ImagUnit * qz * dh ) &
- sinh_quot(iz,irec2) * sinh( ImagUnit * qz * dh ) )
end if
enddo
end do
end do
end do ! irec2
else ! numerical treatment necessary
do iz = -nzdh+1, nzdh-1 ! for z=+-D/2 this is zero anyways
do iqz = -stars%mx3, stars%mx3
irec3 = stars%ig(stars%kv2(1,irec2),stars%kv2(2,irec2),iqz)
if ( irec3 /= 0 ) then ! use only stars within the g_max sphere
c_ph(iqz,irec2) = stars%rgphs(stars%kv2(1,irec2),stars%kv2(2,irec2),iqz)
qz = iqz * cell%bmat(3,3)
VIz(iz,irec2) = VIz(iz,irec2) &
+ vcons1(irec3) * c_ph(iqz,irec2) * &
( exp( ImagUnit * qz * z(iz) ) &
- exp_m(iz,irec2) * exp( ImagUnit * qz * dh ) &
- exp_p(iz,irec2) * exp( - ImagUnit * <