diff --git a/cdn/cdntot.f90 b/cdn/cdntot.f90 index 7b9bc92ce1d10440427a415850bfb3e34ddc033b..af7ec2da5748b7b1e56c50ad3d85b5a83654ff25 100644 --- a/cdn/cdntot.f90 +++ b/cdn/cdntot.f90 @@ -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) diff --git a/main/mix.F90 b/main/mix.F90 index 395089bc1fb5889dbb3d6cf227c0ae38f7ee51b7..e3932d063847df49df890511222e48d71341d5f1 100644 --- a/main/mix.F90 +++ b/main/mix.F90 @@ -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") diff --git a/mix/kerker.F90 b/mix/kerker.F90 index 07fdd222e931cd1440cbbab107971ed64f30a24c..26c2678f327152199e4a1c1082a805d04642ede7 100644 --- a/mix/kerker.F90 +++ b/mix/kerker.F90 @@ -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 diff --git a/types/types_potden.F90 b/types/types_potden.F90 index ba6042b53cf4e343d113ea90485ae033b0362df8..1a127ce9a3814ec5b24e8d3fad4d88be859c3a08 100644 --- a/types/types_potden.F90 +++ b/types/types_potden.F90 @@ -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 diff --git a/vgen/VYukawaFilm.f90 b/vgen/VYukawaFilm.f90 index 148b2ab489bd1219d32d83a4bb977d30432302c5..0aa327063cb89af833b95235a0ef9afb0d13af41 100644 --- a/vgen/VYukawaFilm.f90 +++ b/vgen/VYukawaFilm.f90 @@ -1,22 +1,24 @@ 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 * qz * dh ) ) + end if + enddo + end do + end if + end do ! 1D FOURIER TRANSFORM TO FIND THE COEFFICIENTS V^I(q_xy,q_z) @@ -377,8 +215,149 @@ module m_VYukawaFilm end do end do + end subroutine VYukawaFilmInterstitial + subroutine VYukawaModify( stars, vacuum, cell, sym, input, mpi, atoms, sphhar, dimension, oneD, noco, den, & + VYukawa ) + + ! This subroutine adds a potential to the previously computed Yukawa + ! potential to ensure charge neutrality. + ! The added potential itself is a solution to the modified Helmholtz + ! equation with constant right-hand side, where the constant is chosen such + ! that charge neutrality is obtained. + + use m_constants + use m_types + use m_vmts + use m_constants + use m_cdntot + 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 + type(t_mpi), intent(in) :: mpi + type(t_atoms), intent(in) :: atoms + type(t_sphhar), intent(in) :: sphhar + type(t_dimension), intent(in) :: dimension + type(t_oneD), intent(in) :: oneD + type(t_noco), intent(in) :: noco + type(t_potden), intent(inout) :: den + type(t_potden), intent(inout) :: VYukawa + + integer :: n, lh, irec3, iz, iqz, nfft, nzmax, nzmin, nzdh + real :: q0, qhat, qbar, ldh, partitioning, dh + real :: q(input%jspins), qis(input%jspins), qmt(atoms%ntype,input%jspins),qvac(2,input%jspins), qtot, qistot + complex :: psq(stars%ng3) + type(t_potden) :: VYukawaModification + real, allocatable :: z(:) + real :: VIzReal(3*stars%mx3), VIzImag(3*stars%mx3) + complex, allocatable :: VIz(:) + complex :: VIqz(-stars%mx3:stars%mx3) + + + ! DEFINITIONS / ALLOCATIONS / INITIALISATIONS + + ! constants: + dh = cell%z1 ! half the width of the film + ! indexing of 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( dh / cell%amat(3,3) * nfft ) - 1 ! index of maximal z below D/2 + allocate( z(nzmin:nzmax) ) + ! definition of 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 + + + ! INTEGRATION OF THE PREVIOUSLY COMPUTED YUKAWA POTENTIAL + + ! initialise VYukawaModification with in-going VYukawa and prepare for integration + call VYukawaModification%init( stars, atoms, sphhar, vacuum, noco, input%jspins, 4 ) + call VYukawaModification%copyPotDen( VYukawa ) + do n = 1, atoms%ntype + do lh = 0, sphhar%nlhd + VYukawaModification%mt(1:atoms%jri(n),lh,n,1) = VYukawaModification%mt(1:atoms%jri(n),lh,n,1) * atoms%rmsh(1:atoms%jri(n),n) ** 2 + end do + end do + + ! integrate the potential over the film region + call cdntot_integrate( stars, atoms, sym, vacuum, input, cell, oneD, VYukawaModification, q, qis, qmt, qvac, qtot, qistot ) + q0 = qtot / cell%area + ldh = input%preconditioning_param * cell%z1 + qhat = ( q0 / ( 2 * dh ) ) / ( sinh(ldh) / ( ldh * cosh( ldh ) ) - 1 ) + qbar = input%preconditioning_param ** 2 / fpi_const * qhat + + + ! SET UP CONSTANT CHARGE DENSITY + + ! instead of den%pw(1,1) = qbar we directly set the pseudo charge density + den%mt = 0; den%pw = 0; den%vacxy = 0; den%vacz = 0 + do n = 1, atoms%ntype + den%mt(1:atoms%jri(n),0,n,1) = sfp_const * qbar * atoms%rmsh(1:atoms%jri(n),n) ** 2 + end do + psq = cmplx(0.0,0.0); psq(1) = qbar + + + ! CALCULATE THE INTERSTITIAL POTENTIAL AS A FUNCTION OF z + + ! initialise and calculate out-going modification potential; reuse VYukawaModification + VYukawaModification%mt = 0; VYukawaModification%pw = 0; VYukawaModification%vacxy = 0; VYukawaModification%vacz = 0 + + allocate( VIz(nzmin:nzmax) ) + VIz = (0.,0.) + do iz = -nzdh+1, nzdh-1 + VIz(iz) = qhat * ( 1 - cosh( input%preconditioning_param * z(iz) ) / cosh( ldh ) ) + end do + + + ! 1D FOURIER TRANSFORM TO FIND THE 3D-FOURIER COEFFICIENTS + + ! 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) ) + + call cfft( VIzReal(:), VIzImag(:), nfft, nfft, nfft, -1 ) + + ! reorder + VIqz = 0 + VIqz(0) = cmplx( VIzReal(1), VIzImag(1) ) + do iqz = 1, stars%mx3 + VIqz( iqz) = cmplx( VIzReal(iqz+1 ), VIzImag(iqz+1 ) ) + VIqz(-iqz) = cmplx( VIzReal(nfft+1-iqz), VIzImag(nfft+1-iqz) ) + end do + + ! add the computed components + do iqz= -stars%mx3, stars%mx3 + irec3 = stars%ig(stars%kv2(1,1),stars%kv2(2,1),iqz) + if ( irec3 /= 0 ) VYukawaModification%pw(irec3,1) = VYukawaModification%pw(irec3,1) + VIqz(iqz) * partitioning / ( stars%nstr(irec3) / stars%nstr2(1) ) + end do + + + ! MUFFIN-TIN POTENTIAL + + call Vmts( input, mpi, stars, sphhar, atoms, sym, cell, oneD, & + VYukawaModification%pw(:,1), den%mt(:,0:,:,1), VYukawaModification%potdenType, & + VYukawaModification%mt(:,0:,:,1) ) + + + ! APPLYING THE MODIFICATION TO THE YUKAWA POTENTIAL + + call VYukawa%AddPotDen( VYukawa, VYukawaModification ) + + + end subroutine VYukawaModify + + + end module m_VYukawaFilm