Fix Bugs in Film-Case Preconditioner

main/mix.F90

@@ -75,6 +75,7 @@ contains
     type(t_potden)                   :: resDen, vYukawa
     integer                          :: ierr(2)
 
+
     !External functions
     real :: CPP_BLAS_sdot
     external :: CPP_BLAS_sdot
@@ -230,7 +231,7 @@ contains
 
     ! KERKER PRECONDITIONER
     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 )
       MPI0_b: if( mpi%irank == 0 ) then 
         call resDen%subPotDen( outDen, inDen )
@@ -258,6 +259,8 @@ 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 )

vgen/VYukawaFilm.f90

 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 
 
@@ -65,8 +75,8 @@ module m_VYukawaFilm
                                 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
-    ! 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
+    ! 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
@@ -95,26 +105,27 @@ module m_VYukawaFilm
     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, kz
+    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 / initialisations:
+    ! DEFINITIONS / ALLOCATIONS / INITIALISATIONS    
+    
     do iz = 1, vacuum%nmz
-      z(iz) = ( iz - 1 ) * vacuum%delz ! z_i starting at 0
+      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) ) ! 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) )
       end do
       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
     sum_qz = (0.,0.)
     VVxy = (0.,0.)
@@ -126,14 +137,14 @@ module m_VYukawaFilm
     do irec2 = 1, stars%ng2
       do ivac = 1, vacuum%nvac
         sign = 3. - 2. * ivac
-        do kz = -stars%mx3, stars%mx3
-          irec3 = stars%ig(stars%kv2(1,irec2),stars%kv2(2,irec2),kz)
+        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(kz,irec2) = stars%rgphs(stars%kv2(1,irec2),stars%kv2(2,irec2),kz)
-            qz = kz * cell%bmat(3,3)
-            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) )
+            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
@@ -149,12 +160,10 @@ module m_VYukawaFilm
 
     ! case irec2 > 1:
     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
         ! 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)
+        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) * expDhg(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'
@@ -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) ) &
                                     + 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) 
+        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
 
     ! 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
-      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)
+      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) * expDhg(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) ! 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
       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) 
+      VVz(1:vacuum%nmz,ivac) = VVz(1:vacuum%nmz,ivac) + vcons(1) * zeta(1:vacuum%nmz) * expDhg(1)
     end do
 
 
@@ -206,9 +213,15 @@ module m_VYukawaFilm
                                       psq, VVxy, VVz, alphm, &
                                       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
-    ! 2. part: add the contribution from the vacuum charge density to the interstitial potential, which had already been computed earlier for the vacuum potential
-    ! 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 )
+    ! 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
 
     use m_ExpSave
     use m_constants
@@ -228,109 +241,136 @@ module m_VYukawaFilm
 
     complex,        intent(out) :: VIq(stars%ng3)
 
-    real                        :: partitioning, rz, qz, fit, q
-    integer                     :: irec2, irec3, iz, jz, ivac, kz, ifft
-    complex                     :: VIz(3*stars%mx3,stars%ng2), sum_qz(3*stars%mx3,stars%ng2), eta(3*stars%mx3,stars%ng2)
+    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), exp_m(3*stars%mx3,stars%ng2), exp_p(3*stars%mx3,stars%ng2)
-    real                        :: z(3*stars%mx3)
-    real                        :: g_damped(stars%ng2), vcons2(stars%ng2)
-
-
-    ! definitions / initialisations
-    ifft = 3 * stars%mx3
-    partitioning = 1. / real( ifft )
-    do iz = 1, ifft
-      z(iz) = cell%amat(3,3) * ( iz - 1 ) * partitioning               ! z_i are equidistantly distributed along ( 0, cell%amat(3,3) ]
-      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 ]
+    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)
+    
+
+    ! DEFINITIONS / ALLOCATIONS / INITIALISATIONS
+
+    ! 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
+    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
+    ! 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
       g_damped(irec2) = sqrt( stars%sk2(irec2) ** 2 + input%preconditioning_param ** 2 )
       vcons2(irec2) = -1. / ( 2. * g_damped(irec2) )
-      do iz = 1, ifft
-        exp_m(iz,irec2) = exp_save( - g_damped(irec2) * ( cell%z1 + z(iz) ) )
-        exp_p(iz,irec2) = exp_save( - g_damped(irec2) * ( cell%z1 - z(iz) ) )
+      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) )
     end do
     do irec3 = 1, stars%ng3
       vcons1(irec3) = fpi_const * psq(irec3) / ( stars%sk3(irec3) ** 2 + input%preconditioning_param ** 2 )
     end do
-    sum_qz = (0.,0.)
     VIz = (0.,0.)
     VIq = (0.,0.)
 
 
     ! 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 = 1, ifft
-
-        ! in the transition region ( D/2, D~/2 ), smooth out the potential numerically
-        if ( abs( z(iz) ) >= cell%z1 ) then
-          ivac = 1 ! upper vacuum
-          if ( z(iz) < 0.0 .and. .not. ( sym%invs .or. sym%zrfs ) ) ivac = 2 ! lower vacuum
-          rz = ( abs( z(iz) ) - cell%z1 ) / vacuum%delz + 1.0 ! numbering the grid points in the vacuum region
-          jz = rz
-          q = rz - jz
-          if ( irec2 == 1 ) then
-            fit = 0.5     * ( q - 1. ) * ( q - 2. ) * VVz(jz,  ivac) &
-                -       q *              ( q - 2. ) * VVz(jz+1,ivac) &
-                + 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) )
+    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
-
-        ! 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
+        enddo
       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 )
-      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 )
-      VIzReal =  real( VIz )
-      VIzImag = aimag( VIz )
-      call cfft( VIzReal(:,irec2), VIzImag(:,irec2), ifft, ifft, ifft, -1 )
-      ! reorder:
+    ! reorder
       VIqz(0,irec2) = cmplx( VIzReal(1,irec2), VIzImag(1,irec2) )
-      do kz = 1, stars%mx3
-        VIqz( kz,irec2) = cmplx( VIzReal(kz+1,     irec2), VIzImag(kz+1,     irec2) )
-        VIqz(-kz,irec2) = cmplx( VIzReal(ifft+1-kz,irec2), VIzImag(ifft+1-kz,irec2) )
+      do iqz = 1, stars%mx3
+        VIqz( iqz,irec2) = cmplx( VIzReal(iqz+1,     irec2), VIzImag(iqz+1,     irec2) )
+        VIqz(-iqz,irec2) = cmplx( VIzReal(nfft+1-iqz,irec2), VIzImag(nfft+1-iqz,irec2) )
       end do
-      ! add the computed components to V^I(q_xy,q_z):
-      do kz= -stars%mx3, stars%mx3
-        irec3 = stars%ig(stars%kv2(1,irec2),stars%kv2(2,irec2),kz)
-        VIq(irec3) = VIq(irec3) + VIqz(kz,irec2) / partitioning * stars%nstr(irec3) / stars%nstr2(irec2)
+   
+    ! add the computed components to V^I(q_xy,q_z):
+      do iqz= -stars%mx3, stars%mx3
+        irec3 = stars%ig(stars%kv2(1,irec2),stars%kv2(2,irec2),iqz)
+        if ( irec3 /= 0 ) VIq(irec3) = VIq(irec3) + VIqz(iqz,irec2) * partitioning / ( stars%nstr(irec3) / stars%nstr2(irec2) )
       end do
     end do
 
-
   end subroutine VYukawaFilmInterstitial
 
 

vgen/psqpw.F90

@@ -152,6 +152,7 @@ contains
 #endif
 
     if ( mpi%irank == 0 ) then
+      if ( potdenType == POTDEN_TYPE_POTYUK ) return
       ! Check: integral of the pseudo charge density within the slab
       if ( input%film .and. .not. oneD%odi%d1 ) then
         psint = psq(1) * stars%nstr(1) * vacuum%dvac

vgen/vgen_coulomb.F90

@@ -67,7 +67,6 @@ contains
     integer:: ierr
 #endif
  
-
     
     allocate ( alphm(stars%ng2,2), af1(3*stars%mx3), bf1(3*stars%mx3), psq(stars%ng3)  )
     vCoul%iter = den%iter
@@ -135,6 +134,8 @@ contains
           !                bf1(i_sm) = bf1(i_sm) + z * deltb
           !              ENDDO
           !            ENDIF
+
+
           !        --> 1-d fourier transform and store the coefficients in vTot%pw( ,1)
           call cfft( af1, bf1, ivfft, ivfft, ivfft, -1 )
           !            delta = ivfft * delta * 2 / fpi ! * amat(3,3)**2 * ani