LCOV - code coverage report
Current view: top level - vgen - VYukawaFilm.f90 (source / functions) Hit Total Coverage
Test: combined.info Lines: 0 371 0.0 %
Date: 2019-09-08 04:53:50 Functions: 0 6 0.0 %

          Line data    Source code
       1             : module m_VYukawaFilm
       2             : 
       3             :   ! Computation of the film-case Yukawa potential for the preconditioning of the
       4             :   ! residual charge density in 5 steps:
       5             :   ! 1. pseudo-charge density generation
       6             :   ! 2. vacuum potential generation
       7             :   ! 3. interstitial potential generation
       8             :   ! 4. muffin-tin potential generation
       9             :   ! 5. modification for charge neutrality
      10             : 
      11             :   ! The Yukawa potential is the solution to the modified Helmholtz equation
      12             :   ! ( Delta - lambda^2 ) V_lambda = -4 pi ( rho_out - rho_in )
      13             :   ! subject to some conditions. 
      14             :   ! The general scheme (steps 1 to 4) is the same as for the Poisson equation --
      15             :   ! we use Green function methods for the z-dependent vacuum and interstitial
      16             :   ! potentials as well as for the muffin-tin potential and apply Weinert's
      17             :   ! method.
      18             :   ! You can choose between two variants:
      19             :   ! 1. variant: 
      20             :   ! zero Dirichlet boundary conditions at +/- infinity; 
      21             :   ! multiplication with a decaying exponential in vacuum;
      22             :   ! modification in the film for charge neutrality (step 5)
      23             :   ! 2. variant:
      24             :   ! zero Dirichlet boundary conditions near the film boundary in vacuum (D/2+2R);
      25             :   ! modification in the film for charge neutrality (step 5) 
      26             :   ! In both cases charge neutrality is broken.
      27             :   ! To restore charge neutrality, we need the integral over the potential to be
      28             :   ! zero.
      29             :   ! In step 5 we therefore solve the modified Helmholtz equation again with
      30             :   ! constant right-hand side, for an additive correction to the potential.
      31             :   ! The constant is chosen such that the integral over the final potential is
      32             :   ! zero.
      33             : 
      34             :   contains 
      35             : 
      36             : 
      37             : 
      38           0 :   subroutine VYukawaFilm( stars, vacuum, cell, sym, input, mpi, atoms, sphhar, oneD, noco, den, &
      39             :                           VYukawa )
      40             : 
      41             :     use m_constants
      42             :     use m_types
      43             :     use m_psqpw
      44             :     use m_vmts
      45             :     implicit none
      46             : 
      47             :     type(t_stars),      intent(in)    :: stars
      48             :     type(t_vacuum),     intent(in)    :: vacuum
      49             :     type(t_cell),       intent(in)    :: cell
      50             :     type(t_sym),        intent(in)    :: sym
      51             :     type(t_input),      intent(in)    :: input
      52             :     type(t_mpi),        intent(in)    :: mpi
      53             :     type(t_atoms),      intent(in)    :: atoms 
      54             :     type(t_sphhar),     intent(in)    :: sphhar
      55             :     type(t_oneD),       intent(in)    :: oneD
      56             :     type(t_noco),       intent(in)    :: noco
      57             :     type(t_potden),     intent(inout) :: den
      58             : 
      59             :     type(t_potden),     intent(inout) :: VYukawa
      60             : 
      61           0 :     complex                           :: psq(stars%ng3)
      62           0 :     complex                           :: alphm(stars%ng2,2)
      63             :     real                              :: dh_prec
      64           0 :     real                              :: coshdh(stars%ng2)
      65             : 
      66             :  
      67             :     ! PSEUDO-CHARGE DENSITY
      68             : 
      69             :     call psqpw( mpi, atoms, sphhar, stars, vacuum, cell, input, sym, oneD, &
      70             :                 den%pw(:,1), den%mt(:,:,:,1), den%vacz(:,:,1), .false., VYukawa%potdenType, &
      71           0 :                 psq )
      72             : 
      73             : 
      74             :     ChooseVariant: if ( .true. ) then
      75             : 
      76             :       ! VACUUM POTENTIAL
      77             : 
      78             :       call VYukawaFilmVacuumVariant1( &
      79             :               stars, vacuum, cell, sym, input, atoms%rmt(1), &
      80             :               psq, den%vacxy(:,:,:,1), den%vacz(:,:,1), &
      81           0 :               VYukawa%vacxy, VYukawa%vacz, alphm )
      82             : 
      83             : 
      84             :       ! INTERSTITIAL POTENTIAL
      85             : 
      86             :       call VYukawaFilmInterstitialVariant1( &
      87             :               stars, vacuum, cell, sym, input, &
      88             :               psq, VYukawa%vacxy, VYukawa%vacz, alphm, &
      89           0 :               VYukawa%pw(:,1) )
      90             : 
      91             :     else ChooseVariant
      92             : 
      93             :       ! VACUUM POTENTIAL
      94             : 
      95             :       call VYukawaFilmVacuumVariant2( &
      96             :               stars, vacuum, cell, sym, input, 2*atoms%rmt(1), &
      97             :               psq, den%vacxy(:,:,:,1), den%vacz(:,:,1), &
      98             :               VYukawa%vacxy, VYukawa%vacz, alphm, dh_prec, coshdh )
      99             : 
     100             : 
     101             :       ! INTERSTITIAL POTENTIAL
     102             : 
     103             :       call VYukawaFilmInterstitialVariant2( &
     104             :               stars, vacuum, cell, sym, input, &
     105             :               psq, VYukawa%vacxy, VYukawa%vacz, alphm, dh_prec, coshdh, &
     106             :               VYukawa%pw(:,1) )
     107             : 
     108             :     end if ChooseVariant
     109             : 
     110             : 
     111             :     ! MUFFIN-TIN POTENTIAL
     112             : 
     113             :     call Vmts( input, mpi, stars, sphhar, atoms, sym, cell, oneD, &
     114             :                VYukawa%pw(:,1), den%mt(:,0:,:,1), VYukawa%potdenType, &
     115           0 :                VYukawa%mt(:,0:,:,1) )
     116             : 
     117             :  
     118             :     ! MODIFICATION FOR CHARGE NEUTRALITY
     119             : 
     120             :     call VYukawaModify( stars, vacuum, cell, sym, input, mpi, atoms, sphhar, oneD, noco, &
     121             :                         den, &
     122           0 :                         VYukawa )
     123             : 
     124             : 
     125           0 :   end subroutine VYukawaFilm
     126             : 
     127             : 
     128             : 
     129           0 :   subroutine VYukawaFilmVacuumVariant1( &
     130             :                 stars, vacuum, cell, sym, input, rmt, &
     131           0 :                 psq, rhtxy, rht, &
     132           0 :                 VVxy, VVz, alphm )
     133             : 
     134             :     ! 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)
     135             :     ! 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
     136             : 
     137             :     use m_ExpSave
     138             :     use m_constants
     139             :     use m_types
     140             :     use m_intgr, only: intgz1Reverse
     141             :     use m_qsf
     142             :     implicit none
     143             : 
     144             :     type(t_stars),  intent(in)  :: stars
     145             :     type(t_vacuum), intent(in)  :: vacuum
     146             :     type(t_cell),   intent(in)  :: cell
     147             :     type(t_sym),    intent(in)  :: sym
     148             :     type(t_input),  intent(in)  :: input
     149             :     real,           intent(in)  :: rmt
     150             :     complex,        intent(in)  :: psq(stars%ng3)
     151             :     complex,        intent(in)  :: rhtxy(vacuum%nmzxyd,stars%ng2-1,2)
     152             :     real,           intent(in)  :: rht(vacuum%nmzd,2) 
     153             : 
     154             :     complex,        intent(out) :: VVxy(vacuum%nmzxyd,2:stars%ng2,2) ! this is the qxy /= 0 part of the vacuum potential
     155             :     real,           intent(out) :: VVz(vacuum%nmzd,2)                ! this is the qxy  = 0 part of the vacuum potential
     156             :     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) 
     157             : 
     158           0 :     complex                     :: sum_qz(2,stars%ng2)
     159           0 :     complex                     :: c_ph(-stars%mx3:stars%mx3,stars%ng2)
     160             :     complex                     :: signIqz
     161           0 :     real                        :: g_damped(stars%ng2), qz, sign, vcons(stars%ng2) 
     162           0 :     real                        :: exp_m(vacuum%nmzd,stars%ng2), exp_p(vacuum%nmzd,stars%ng2)
     163           0 :     real                        :: expDhg(stars%ng2), expDg(stars%ng2)
     164           0 :     real                        :: z(vacuum%nmzd)
     165             :     integer                     :: iz, irec2, irec3, ivac, iqz
     166           0 :     complex                     :: fa(vacuum%nmzxyd,2:stars%ng2), fb(vacuum%nmzxyd,2:stars%ng2)
     167           0 :     complex                     :: alpha(vacuum%nmzxyd,2:stars%ng2,2), beta(vacuum%nmzxyd,2:stars%ng2,2), gamma(vacuum%nmzxyd,2:stars%ng2)
     168           0 :     real                        :: ga(vacuum%nmzd), gb(vacuum%nmzd)
     169           0 :     real                        :: delta(vacuum%nmzd,2), epsilon(vacuum%nmzd,2), zeta(vacuum%nmzd)
     170             : 
     171             : 
     172             :     ! DEFINITIONS / ALLOCATIONS / INITIALISATIONS    
     173             :     
     174           0 :     do iz = 1, vacuum%nmz
     175           0 :       z(iz) = ( iz - 1 ) * vacuum%delz
     176             :     end do
     177           0 :     do irec2 = 1, stars%ng2
     178           0 :       g_damped(irec2) = sqrt( stars%sk2(irec2) ** 2 + input%preconditioning_param ** 2 )
     179           0 :       vcons(irec2) = tpi_const / g_damped(irec2)
     180           0 :       do iz = 1, vacuum%nmz
     181           0 :         exp_m(iz,irec2) = exp_save( - g_damped(irec2) * z(iz) )
     182           0 :         exp_p(iz,irec2) = exp_save(   g_damped(irec2) * z(iz) )
     183             :       end do
     184           0 :       expDhg(irec2) = exp_save( - cell%z1 * g_damped(irec2) )
     185           0 :       expDg(irec2)  = exp_save( -2 * cell%z1 * g_damped(irec2) )
     186             :     end do
     187           0 :     sum_qz = (0.,0.)
     188           0 :     VVxy = (0.,0.)
     189           0 :     VVz = 0.
     190             :    
     191             : 
     192             :     ! CONTRIBUTION FROM THE INTERSTITIAL CHARGE DENSITY
     193             : 
     194           0 :     do irec2 = 1, stars%ng2
     195           0 :       do ivac = 1, vacuum%nvac
     196           0 :         sign = 3. - 2. * ivac
     197           0 :         do iqz = -stars%mx3, stars%mx3
     198           0 :           irec3 = stars%ig(stars%kv2(1,irec2),stars%kv2(2,irec2),iqz)
     199             :           ! use only stars within the g_max sphere -> stars outside the sphere have per definition index ig3n = 0
     200           0 :           if ( irec3 /= 0 ) then
     201           0 :             c_ph(iqz,irec2) = stars%rgphs(stars%kv2(1,irec2),stars%kv2(2,irec2),iqz)
     202           0 :             qz = iqz * cell%bmat(3,3)
     203           0 :             signIqz = sign * ImagUnit * qz 
     204           0 :             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) )
     205             :           endif
     206             :         enddo
     207           0 :         if( irec2 /= 1 ) then
     208           0 :           VVxy(1:vacuum%nmzxy,irec2,ivac) = vcons(irec2) * sum_qz(ivac,irec2) * exp_m(1:vacuum%nmzxy,irec2)
     209             :         else
     210           0 :           VVz(1:vacuum%nmz,ivac)          = vcons(1)     * sum_qz(ivac,1)     * exp_m(1:vacuum%nmz,1)
     211             :         end if
     212             :       enddo
     213             :     enddo
     214             : 
     215             : 
     216             :     ! CONTRIBUTION FROM THE VACUUM CHARGE DENSITY
     217             : 
     218             :     ! shifting z:
     219           0 :     do irec2 = 1, stars%ng2
     220           0 :       exp_m(1:vacuum%nmz,irec2) = exp_m(1:vacuum%nmz,irec2) * expDhg(irec2)
     221           0 :       exp_p(1:vacuum%nmz,irec2) = exp_p(1:vacuum%nmz,irec2) / expDhg(irec2)
     222             :     end do
     223             : 
     224             :     ! case irec2 > 1:
     225           0 :     do irec2 = 2, stars%ng2
     226           0 :       do ivac = 1, vacuum%nvac
     227             :         ! integrands:
     228           0 :         fa(1:vacuum%nmzxy,irec2) = rhtxy(1:vacuum%nmzxy,irec2-1,ivac) * exp_m(1:vacuum%nmzxy,irec2)
     229           0 :         fb(1:vacuum%nmzxy,irec2) = rhtxy(1:vacuum%nmzxy,irec2-1,ivac) * exp_p(1:vacuum%nmzxy,irec2)
     230             :         ! integrals:
     231             :         ! alpha(z,q_xy,ivac) = int_z^infty rho(z',q_xy,ivac) exp(-sqrt(q_xy**2+prec_param**2)*z') dz'
     232             :         ! beta (z,q_xy,ivac) = int_{D/2}^z rho(z',q_xy,ivac) exp(+sqrt(q_xy**2+prec_param**2)*z') dz'
     233             :         ! 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)
     234           0 :         call intgz1Reverse( fa(:,irec2), vacuum%delz, vacuum%nmzxy, alpha(:,irec2,ivac), .true. ) 
     235           0 :         call qsfComplex( vacuum%delz, fb(:,irec2), beta(:,irec2,ivac), vacuum%nmzxy, 1 )
     236             :         ! alphm(q_xy,ivac) = alpha(D/2,q_xy,ivac) --- these integrals are also needed for the interstitial potential
     237           0 :         alphm(irec2,ivac) = alpha(1,irec2,ivac)
     238             :       end do
     239           0 :       if ( vacuum%nvac == 1 ) then
     240           0 :         if ( sym%invs ) then
     241           0 :           alphm(irec2,2) = conjg( alphm(irec2,1) )
     242             :         else
     243           0 :           alphm(irec2,2) = alphm(irec2,1)
     244             :         end if
     245             :       end if
     246           0 :       do ivac = 1, vacuum%nvac
     247             :         gamma(1:vacuum%nmzxy,irec2) = exp_m(1:vacuum%nmzxy,irec2) * ( alphm(irec2,mod(ivac,2)+1) + beta(1:vacuum%nmzxy,irec2,ivac) ) &
     248           0 :                                     + exp_p(1:vacuum%nmzxy,irec2) *                               alpha(1:vacuum%nmzxy,irec2,ivac) ! mod(ivac,2)+1 outputs the other ivac value 
     249           0 :         where ( 2. * gamma(:,irec2) == gamma(:,irec2) ) gamma(:,irec2) = cmplx( 0., 0. )
     250           0 :         VVxy(1:vacuum%nmzxy,irec2,ivac) = VVxy(1:vacuum%nmzxy,irec2,ivac) + vcons(irec2) * gamma(1:vacuum%nmzxy,irec2) 
     251             :       end do
     252             :     end do
     253             : 
     254             :     ! case irec2 = 1:
     255           0 :     do ivac = 1, vacuum%nvac
     256           0 :       ga(1:vacuum%nmz) = rht(1:vacuum%nmz,ivac) * exp_m(1:vacuum%nmz,1)
     257           0 :       gb(1:vacuum%nmz) = rht(1:vacuum%nmz,ivac) * exp_p(1:vacuum%nmz,1)
     258           0 :       call intgz1Reverse( ga(:), vacuum%delz, vacuum%nmz, delta(:,ivac), .true. ) ! integrals 
     259           0 :       call qsf( vacuum%delz, gb(:), epsilon(:,ivac), vacuum%nmz, 1 )
     260           0 :       alphm(1,ivac) = delta(1,ivac)
     261             :     end do
     262           0 :     if ( vacuum%nvac == 1 ) alphm(1,2) = alphm(1,1)
     263           0 :     do ivac = 1, vacuum%nvac
     264             :       zeta(1:vacuum%nmz) = exp_m(1:vacuum%nmz,1) * ( alphm(1,mod(ivac,2)+1) + epsilon(1:vacuum%nmz,ivac) ) &
     265           0 :                          + exp_p(1:vacuum%nmz,1) *                              delta(1:vacuum%nmz,ivac) 
     266           0 :       where ( 2. * zeta == zeta ) zeta = 0.
     267           0 :       VVz(1:vacuum%nmz,ivac) = VVz(1:vacuum%nmz,ivac) + vcons(1) * zeta(1:vacuum%nmz)
     268             :     end do
     269             : 
     270             :     ! damping in vacuum:
     271           0 :     do ivac = 1, vacuum%nvac
     272           0 :       VVz(1:vacuum%nmz,ivac) = VVz(1:vacuum%nmz,ivac) * exp( -0.1 / rmt * z(1:vacuum%nmz) )
     273           0 :       do irec2 = 2, stars%ng2
     274           0 :         VVxy(1:vacuum%nmzxy,irec2,ivac) = VVxy(1:vacuum%nmzxy,irec2,ivac) * exp( -0.1 / rmt * z(1:vacuum%nmzxy) )
     275             :       end do
     276             :     end do
     277             : 
     278             : 
     279           0 :   end subroutine VYukawaFilmVacuumVariant1
     280             : 
     281             : 
     282             : 
     283           0 :   subroutine VYukawaFilmInterstitialVariant1( &
     284             :                 stars, vacuum, cell, sym, input, &
     285           0 :                 psq, VVxy, VVz, alphm, &
     286           0 :                 VIq )
     287             : 
     288             :     ! main parts:
     289             :     ! 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)
     290             :     ! 2. part: Add the contribution from the vacuum charge density to the interstitial potential, which had already been computed earlier for the vacuum potential
     291             :  
     292             :     ! 4. part: Compute the coefficients V^I(q_xy,q_z) from the function V^I(q_xy,z) by a 1D Fourier transform:
     293             :     !          V^I(q_xy,z) = sum_{q_z} V^I(q_xy,q_z) * exp( ImagUnit * q_z * z )
     294             :     ! In order to be able to match the interstitial and vacuum potentials smoothly at the interface, the Fourier transform is done
     295             :     ! in a slightly larger region. -> 3. part
     296             :     ! 3. part: Interpolate the vacuum potential in a small region surrounding the slab
     297             : 
     298             :     use m_ExpSave
     299             :     use m_constants
     300             :     use m_types
     301             :     use m_cfft
     302             :     implicit none
     303             : 
     304             :     type(t_stars),  intent(in)  :: stars
     305             :     type(t_vacuum), intent(in)  :: vacuum
     306             :     type(t_cell),   intent(in)  :: cell
     307             :     type(t_sym),    intent(in)  :: sym
     308             :     type(t_input),  intent(in)  :: input
     309             :     complex,        intent(in)  :: psq(stars%ng3)
     310             :     complex,        intent(in)  :: VVxy(vacuum%nmzxyd,2:stars%ng2,2)
     311             :     real,           intent(in)  :: VVz(vacuum%nmzd,2)
     312             :     complex,        intent(in)  :: alphm(stars%ng2,2)
     313             : 
     314             :     complex,        intent(out) :: VIq(stars%ng3)
     315             : 
     316             :     real                        :: partitioning, rz, qz, q
     317             :     integer                     :: irec2, irec3, iz, jz, ivac, iqz, nfft, nzmax, nzmin, nzdh, nLower, nUpper, jvac
     318           0 :     complex, allocatable        :: VIz(:,:), eta(:,:)
     319           0 :     complex                     :: VIqz(-stars%mx3:stars%mx3,stars%ng2), c_ph(-stars%mx3:stars%mx3,stars%ng2)
     320           0 :     complex                     :: vcons1(stars%ng3)
     321           0 :     real                        :: VIzReal(3*stars%mx3,stars%ng2), VIzImag(3*stars%mx3,stars%ng2)
     322           0 :     real, allocatable           :: exp_m(:,:), exp_p(:,:)
     323           0 :     real, allocatable           :: z(:)
     324           0 :     real                        :: g_damped(stars%ng2), vcons2(stars%ng2), expDhg(stars%ng2)
     325             :     
     326             : 
     327             :     ! DEFINITIONS / ALLOCATIONS / INITIALISATIONS
     328             : 
     329             :     ! grid points z_i:
     330           0 :     nfft = 3 * stars%mx3                                  ! number of grid points for Fourier transform
     331           0 :     partitioning = 1. / real( nfft )
     332           0 :     nzmax = nfft / 2                                      ! index of maximal z below D~/2 
     333           0 :     nzmin = nzmax - nfft + 1                              ! index of minimal z above -D~/2
     334           0 :     nzdh = ceiling( cell%z1 / cell%amat(3,3) * nfft ) - 1 ! index of maximal z below D/2
     335           0 :     allocate( z(nzmin:nzmax) )
     336             :     ! definition of z_i:        ! indexing: z_0 = 0; positive indices for z > 0; negative indices for z < 0
     337           0 :     do iz = nzmin, nzmax
     338           0 :       z(iz) = cell%amat(3,3) * iz * partitioning
     339             :     end do
     340             :     ! other variables:
     341           0 :     allocate( VIz(nzmin:nzmax,stars%ng2), eta(nzmin:nzmax,stars%ng2) )
     342           0 :     allocate( exp_m(nzmin:nzmax,stars%ng2), exp_p(nzmin:nzmax,stars%ng2) )
     343           0 :     do irec2 = 1, stars%ng2
     344           0 :       g_damped(irec2) = sqrt( stars%sk2(irec2) ** 2 + input%preconditioning_param ** 2 )
     345           0 :       vcons2(irec2) = -1. / ( 2. * g_damped(irec2) )
     346           0 :       do iz = nzmin, nzmax
     347           0 :         exp_m(iz,irec2) = exp_save( - g_damped(irec2) * z(iz) )
     348           0 :         exp_p(iz,irec2) = exp_save(   g_damped(irec2) * z(iz) )
     349             :       end do
     350           0 :       expDhg(irec2) = exp_save( - cell%z1 * g_damped(irec2) )
     351             :     end do
     352           0 :     do irec3 = 1, stars%ng3
     353           0 :       vcons1(irec3) = fpi_const * psq(irec3) / ( stars%sk3(irec3) ** 2 + input%preconditioning_param ** 2 )
     354             :     end do
     355           0 :     VIz = (0.,0.)
     356           0 :     VIq = (0.,0.)
     357             : 
     358             : 
     359             :     ! CONTRIBUTION FROM THE INTERSTITIAL CHARGE DENSITY
     360             :     
     361             :     ! compute V^I(q_xy,z) as a function of q_xy and z
     362           0 :     do irec2 = 1, stars%ng2  
     363           0 :       do iz = -nzdh, nzdh
     364           0 :         do iqz = -stars%mx3, stars%mx3
     365           0 :           irec3 = stars%ig(stars%kv2(1,irec2),stars%kv2(2,irec2),iqz)
     366           0 :           if ( irec3 /= 0 ) then ! use only stars within the g_max sphere
     367           0 :             c_ph(iqz,irec2) = stars%rgphs(stars%kv2(1,irec2),stars%kv2(2,irec2),iqz)
     368           0 :             qz = iqz * cell%bmat(3,3)
     369             :             VIz(iz,irec2) = VIz(iz,irec2) &
     370             :                           + vcons1(irec3) * c_ph(iqz,irec2) * &
     371             :                             ( exp( ImagUnit * qz * z(iz) ) &
     372             :                             + vcons2(irec2) * expDhg(irec2) * &
     373             :                               ( ( g_damped(irec2) + ImagUnit * qz ) * exp_p(iz,irec2) * exp(   ImagUnit * qz * cell%z1 ) &
     374           0 :                               + ( g_damped(irec2) - ImagUnit * qz ) * exp_m(iz,irec2) * exp( - ImagUnit * qz * cell%z1 ) ) )
     375             :           end if
     376             :         enddo
     377             :       end do
     378             :     ! irec2 loop continues
     379             : 
     380             : 
     381             :     ! CONTRIBUTION FROM THE VACUUM CHARGE DENSITY
     382             : 
     383             :     ! irec2 loop continues
     384           0 :       eta(-nzdh:nzdh,irec2) = exp_m(-nzdh:nzdh,irec2) * alphm(irec2,2) + exp_p(-nzdh:nzdh,irec2) * alphm(irec2,1)
     385           0 :       where ( 2.0 * eta(:,irec2) == eta(:,irec2) ) eta(:,irec2) = cmplx( 0.0, 0.0 )
     386           0 :       VIz(-nzdh:nzdh,irec2) = VIz(-nzdh:nzdh,irec2) + tpi_const / g_damped(irec2) * eta(-nzdh:nzdh,irec2)
     387             :     ! irec2 loop continues
     388             :     
     389             :     
     390             :     ! INTERPOLATION IN VACUUM REGION
     391             : 
     392             :     ! use Lagrange polynomials of order 3 to interpolate the vacuum potential outside I
     393             :     ! q, q-1 and q-2 (scaled with +/-1 or +/-0.5) are the factors of the Lagrange basis polynomials
     394             :     ! irec2 loop continues
     395           0 :       do ivac = 1, 2
     396           0 :         select case( ivac )
     397             :           case( 1 )
     398           0 :             nUpper = nzmax; nLower =  nzdh + 1; jvac = 1
     399             :           case( 2 )
     400           0 :             nLower = nzmin; nUpper = -nzdh - 1; jvac = 2; if ( sym%invs .or. sym%zrfs ) jvac = 1
     401             :         end select
     402           0 :         do iz = nLower, nUpper
     403           0 :           rz = ( abs( z(iz) ) - cell%z1 ) / vacuum%delz + 1.0
     404           0 :           jz = rz      ! index of maximal vacuum grid point below z_i
     405           0 :           q = rz - jz  ! factor in Lagrange basis polynomials
     406           0 :           if ( irec2 == 1 ) then
     407             :             VIz(iz,irec2) = 0.5     * ( q - 1. ) * ( q - 2. ) * VVz(jz,  jvac) &
     408             :                           -       q *              ( q - 2. ) * VVz(jz+1,jvac) &
     409           0 :                           + 0.5 * q * ( q - 1. )              * VVz(jz+2,jvac)
     410           0 :           else if ( jz + 2 <= vacuum%nmzxy ) then
     411             :             VIz(iz,irec2) = 0.5 *     ( q - 1. ) * ( q - 2. ) * VVxy(jz,  irec2,jvac) &
     412             :                           -       q              * ( q - 2. ) * VVxy(jz+1,irec2,jvac) &
     413           0 :                           + 0.5 * q * ( q - 1. )              * VVxy(jz+2,irec2,jvac)
     414           0 :             if ( ( sym%invs .and. .not. sym%zrfs ) .and. ivac == 2 ) VIz(iz,irec2) = conjg( VIz(iz,irec2) )
     415             :           end if
     416             :         end do
     417             :       end do
     418             :     end do ! irec2
     419             : 
     420             : 
     421             :     ! 1D FOURIER TRANSFORM TO FIND THE COEFFICIENTS V^I(q_xy,q_z)
     422             : 
     423             :     ! change the indexing for the subroutine cfft, and split real and imaginary parts
     424           0 :     VIzReal(1:nzmax+1,:) =  real( VIz(0:nzmax,:) ); VIzReal(nzmax+2:nfft,:) =  real( VIz(nzmin:-1,:) )
     425           0 :     VIzImag(1:nzmax+1,:) = aimag( VIz(0:nzmax,:) ); VIzImag(nzmax+2:nfft,:) = aimag( VIz(nzmin:-1,:) )
     426             :  
     427             :     ! V^I(q_xy,z) = sum_{q_z} V^I(q_xy,q_z) * exp( ImagUnit * q_z * z )
     428           0 :     do irec2 = 1, stars%ng2
     429           0 :       call cfft( VIzReal(:,irec2), VIzImag(:,irec2), nfft, nfft, nfft, -1 )
     430             :     ! irec2 loop continues
     431             : 
     432             :     ! reorder
     433           0 :       VIqz(0,irec2) = cmplx( VIzReal(1,irec2), VIzImag(1,irec2) )
     434           0 :       do iqz = 1, stars%mx3
     435           0 :         VIqz( iqz,irec2) = cmplx( VIzReal(iqz+1,     irec2), VIzImag(iqz+1,     irec2) )
     436           0 :         VIqz(-iqz,irec2) = cmplx( VIzReal(nfft+1-iqz,irec2), VIzImag(nfft+1-iqz,irec2) )
     437             :       end do
     438             :    
     439             :     ! add the computed components to V^I(q_xy,q_z):
     440           0 :       do iqz= -stars%mx3, stars%mx3
     441           0 :         irec3 = stars%ig(stars%kv2(1,irec2),stars%kv2(2,irec2),iqz)
     442           0 :         if ( irec3 /= 0 ) VIq(irec3) = VIq(irec3) + VIqz(iqz,irec2) * partitioning / ( stars%nstr(irec3) / stars%nstr2(irec2) )
     443             :       end do
     444             :     end do
     445             : 
     446             : 
     447           0 :   end subroutine VYukawaFilmInterstitialVariant1
     448             : 
     449             : 
     450             : 
     451           0 :   subroutine VYukawaFilmVacuumVariant2( &
     452             :                 stars, vacuum, cell, sym, input, rmt, &
     453           0 :                 psq, rhtxy, rht, &
     454           0 :                 VVxy, VVz, alphm, dh_prec, coshdh )
     455             : 
     456             :     ! 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)
     457             :     ! 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
     458             : 
     459             :     use m_ExpSave
     460             :     use m_constants
     461             :     use m_types
     462             :     use m_intgr, only: intgz1Reverse
     463             :     use m_qsf
     464             :     implicit none
     465             : 
     466             :     type(t_stars),  intent(in)  :: stars
     467             :     type(t_vacuum), intent(in)  :: vacuum
     468             :     type(t_cell),   intent(in)  :: cell
     469             :     type(t_sym),    intent(in)  :: sym
     470             :     type(t_input),  intent(in)  :: input
     471             :     real,           intent(in)  :: rmt
     472             :     complex,        intent(in)  :: psq(stars%ng3)
     473             :     complex,        intent(in)  :: rhtxy(vacuum%nmzxyd,stars%ng2-1,2)
     474             :     real,           intent(in)  :: rht(vacuum%nmzd,2) 
     475             : 
     476             :     complex,        intent(out) :: VVxy(vacuum%nmzxyd,2:stars%ng2,2) ! this is the qxy /= 0 part of the vacuum potential
     477             :     real,           intent(out) :: VVz(vacuum%nmzd,2)                ! this is the qxy  = 0 part of the vacuum potential
     478             :     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) 
     479             :     real,           intent(out) :: dh_prec
     480             :     real,           intent(out) :: coshdh(stars%ng2)
     481             : 
     482             :     integer                                             :: iz, irec2, irec3, ivac, iqz, nzdhprec, sign
     483             :     real                                                :: dh, qz, qxy_numerics 
     484           0 :     real, allocatable                                   :: z(:)
     485           0 :     complex, dimension(stars%ng2)                       :: sum_qz
     486           0 :     complex, dimension(stars%ng3)                       :: vcons3
     487           0 :     real, dimension(stars%ng2)                          :: g_damped, vcons2
     488           0 :     complex, dimension(-stars%mx3:stars%mx3,stars%ng2)  :: c_ph, quotq
     489           0 :     real, allocatable                                   :: quotz(:,:), sinhz(:,:)
     490           0 :     real, dimension(-1:1,stars%ng2)                     :: quotvardh
     491           0 :     complex, dimension(-stars%mx3:stars%mx3)            :: expdhqz
     492           0 :     complex, allocatable                                :: fa(:,:), fb(:,:)
     493           0 :     complex, allocatable                                :: alpha(:,:,:), beta(:,:,:), gamma(:,:)
     494           0 :     real, allocatable                                   :: ga(:), gb(:)
     495           0 :     real, allocatable                                   :: delta(:,:), epsilon(:,:), zeta(:)
     496             : 
     497             : 
     498             :     ! DEFINITIONS / ALLOCATIONS / INITIALISATIONS    
     499             :     
     500           0 :     dh = cell%z1 ! half the thickness of the film
     501           0 :     nzdhprec = ceiling( rmt / vacuum%delz ) ! index of dh_prec, see below
     502           0 :     dh_prec = dh + ( nzdhprec - 1 ) * vacuum%delz ! dh_prec is about dh + R; preconditioning boundary
     503           0 :     qxy_numerics = sqrt( ( 100 / dh_prec ) ** 2 - input%preconditioning_param ** 2 )
     504           0 :     allocate( z(nzdhprec) )
     505           0 :     allocate( quotz(-nzdhprec:nzdhprec,stars%ng2) )
     506           0 :     allocate( sinhz(nzdhprec,stars%ng2) )
     507           0 :     do iz = 1, nzdhprec ! new boundary
     508           0 :       z(iz) = ( iz - 1 ) * vacuum%delz + dh
     509             :     end do
     510           0 :     do irec2 = 1, stars%ng2
     511           0 :       g_damped(irec2) = sqrt( stars%sk2(irec2) ** 2 + input%preconditioning_param ** 2 )
     512           0 :       vcons2(irec2) = fpi_const / g_damped(irec2)
     513           0 :       coshdh(irec2) = cosh( g_damped(irec2) * ( dh_prec - dh ) )
     514           0 :       if( stars%sk2(irec2) < qxy_numerics ) then ! numerics ok
     515           0 :         do iz = 1, nzdhprec
     516             :           quotz( iz,irec2) = ( cosh( g_damped(irec2) * z(iz) ) / cosh( g_damped(irec2) * dh_prec ) &
     517           0 :                              + sinh( g_damped(irec2) * z(iz) ) / sinh( g_damped(irec2) * dh_prec ) ) / 2
     518             :           quotz(-iz,irec2) = ( cosh( g_damped(irec2) * z(iz) ) / cosh( g_damped(irec2) * dh_prec ) &
     519           0 :                              - sinh( g_damped(irec2) * z(iz) ) / sinh( g_damped(irec2) * dh_prec ) ) / 2
     520             :         end do
     521             :         quotvardh( 1,irec2) = ( cosh( g_damped(irec2) * dh ) / sinh( g_damped(irec2) * dh_prec ) &
     522           0 :                               + sinh( g_damped(irec2) * dh ) / cosh( g_damped(irec2) * dh_prec ) ) / 2
     523             :         quotvardh(-1,irec2) = ( cosh( g_damped(irec2) * dh ) / sinh( g_damped(irec2) * dh_prec ) &
     524           0 :                               - sinh( g_damped(irec2) * dh ) / cosh( g_damped(irec2) * dh_prec ) ) / 2
     525             :       else ! numerical treatment necessary
     526           0 :         do iz = 1, nzdhprec
     527           0 :           quotz( iz,irec2) = exp_save( g_damped(irec2) * (  z(iz) - dh_prec ) )
     528           0 :           quotz(-iz,irec2) = exp_save( g_damped(irec2) * ( -z(iz) - dh_prec ) )
     529             :         end do
     530           0 :         quotvardh( 1,irec2) = quotz( 1,irec2)
     531           0 :         quotvardh(-1,irec2) = quotz(-1,irec2)
     532             :       end if
     533           0 :       do iz = 1, nzdhprec
     534           0 :         sinhz(iz,irec2) = sinh( g_damped(irec2) * ( dh_prec - z(iz) ) )
     535             :       end do
     536           0 :       do iqz = -stars%mx3, stars%mx3
     537           0 :         quotq(iqz,irec2) = ImagUnit * iqz * cell%bmat(3,3) / g_damped(irec2)
     538             :       end do
     539             :     end do
     540           0 :     sum_qz = (0.,0.)
     541           0 :     VVxy = (0.,0.)
     542           0 :     VVz = 0.
     543           0 :     do irec3 = 1, stars%ng3
     544           0 :       vcons3(irec3) = fpi_const * psq(irec3) / ( stars%sk3(irec3) ** 2 + input%preconditioning_param ** 2 )
     545             :     end do
     546           0 :     do iqz = -stars%mx3, stars%mx3
     547           0 :       expdhqz(iqz) = exp( ImagUnit * iqz * cell%bmat(3,3) * dh )
     548             :     end do
     549             : 
     550             : 
     551             :     ! CONTRIBUTION FROM THE INTERSTITIAL CHARGE DENSITY
     552             : 
     553           0 :     do ivac = 1, vacuum%nvac
     554           0 :       sign = 3 - 2 * ivac
     555           0 :       do irec2 = 1, stars%ng2
     556           0 :         do iqz = -stars%mx3, stars%mx3
     557           0 :           irec3 = stars%ig(stars%kv2(1,irec2),stars%kv2(2,irec2),iqz)
     558             :           ! use only stars within the g_max sphere -> stars outside the sphere have per definition index ig3n = 0
     559           0 :           if ( irec3 /= 0 ) then
     560           0 :             c_ph(iqz,irec2) = stars%rgphs(stars%kv2(1,irec2),stars%kv2(2,irec2),iqz)
     561             :             sum_qz(irec2) = sum_qz(irec2) + &
     562             :                             c_ph(iqz,irec2) * vcons3(irec3) * &
     563             :                             ( ( quotvardh( 1,irec2) - sign * quotq(iqz,irec2) * quotz( 1,irec2) ) * expdhqz( sign*iqz) &
     564           0 :                             - ( quotvardh(-1,irec2) - sign * quotq(iqz,irec2) * quotz(-1,irec2) ) * expdhqz(-sign*iqz) )
     565             :           endif
     566             :         enddo
     567           0 :         if( irec2 /= 1 ) then
     568           0 :           VVxy(1:nzdhprec,irec2,ivac) = sum_qz(irec2) * sign * sinhz(1:nzdhprec,irec2)
     569             :         else
     570           0 :           VVz( 1:nzdhprec,      ivac) = sum_qz(1)     * sign * sinhz(1:nzdhprec,1)
     571             :         end if
     572             :       enddo
     573             :     enddo
     574             : 
     575             : 
     576             :     ! CONTRIBUTION FROM THE VACUUM CHARGE DENSITY
     577             : 
     578           0 :     allocate( fa(nzdhprec,2:stars%ng2), fb(nzdhprec,2:stars%ng2) )
     579           0 :     allocate( alpha(nzdhprec,2:stars%ng2,2), beta(nzdhprec,2:stars%ng2,2), gamma(nzdhprec,2:stars%ng2) )
     580             :     ! case irec2 > 1:
     581           0 :     do irec2 = 2, stars%ng2
     582           0 :       do ivac = 1, vacuum%nvac
     583             :         ! integrands:
     584           0 :         fa(1:nzdhprec,irec2) = rhtxy(1:nzdhprec,irec2-1,ivac) * sinhz(1:nzdhprec,irec2)
     585           0 :         fb(1:nzdhprec,irec2) = rhtxy(1:nzdhprec,irec2-1,ivac) * quotz(1:nzdhprec,irec2)
     586             :         ! integrals:
     587             :         ! alpha(z,q_xy,ivac) = int_z^infty rho(z',q_xy,ivac) sinhz dz'
     588             :         ! beta (z,q_xy,ivac) = int_{D/2}^z rho(z',q_xy,ivac) quotz dz'
     589             :         ! 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)
     590           0 :         call intgz1Reverse( fa(:,irec2), vacuum%delz, nzdhprec, alpha(:,irec2,ivac), .false. ) 
     591           0 :         call qsfComplex( vacuum%delz, fb(:,irec2), beta(:,irec2,ivac), nzdhprec, 1 )
     592             :         ! alphm(q_xy,ivac) = alpha(D/2,q_xy,ivac) --- these integrals are also needed for the interstitial potential
     593           0 :         alphm(irec2,ivac) = alpha(1,irec2,ivac)
     594             :       end do
     595           0 :       if ( vacuum%nvac == 1 ) then
     596           0 :         if ( sym%invs ) then
     597           0 :           alphm(irec2,2) = conjg( alphm(irec2,1) )
     598             :         else
     599           0 :           alphm(irec2,2) = alphm(irec2,1)
     600             :         end if
     601             :       end if
     602           0 :       do ivac = 1, vacuum%nvac
     603             :         gamma(1:nzdhprec,irec2) = quotz(-1:-nzdhprec:-1,irec2) * alphm(irec2,mod(ivac,2)+1) &
     604             :                                 + quotz(1:nzdhprec,irec2) * alpha(1:nzdhprec,irec2,ivac) &
     605           0 :                                 + sinhz(1:nzdhprec,irec2) * beta(1:nzdhprec,irec2,ivac)
     606           0 :         VVxy(1:nzdhprec,irec2,ivac) = VVxy(1:nzdhprec,irec2,ivac) + vcons2(irec2) * gamma(1:nzdhprec,irec2) 
     607             :       end do
     608             :     end do
     609             : 
     610             : 
     611           0 :     allocate( ga(nzdhprec), gb(nzdhprec) )
     612           0 :     allocate( delta(nzdhprec,2), epsilon(nzdhprec,2), zeta(nzdhprec) )
     613             :     ! case irec2 = 1:
     614           0 :     do ivac = 1, vacuum%nvac
     615           0 :       ga(1:nzdhprec) = rht(1:nzdhprec,ivac) * sinhz(1:nzdhprec,1)
     616           0 :       gb(1:nzdhprec) = rht(1:nzdhprec,ivac) * quotz(1:nzdhprec,1)
     617           0 :       call intgz1Reverse( ga(:), vacuum%delz, nzdhprec, delta(:,ivac), .false. ) 
     618           0 :       call qsf( vacuum%delz, gb(:), epsilon(:,ivac), nzdhprec, 1 )
     619           0 :       alphm(1,ivac) = delta(1,ivac)
     620             :     end do
     621           0 :     if ( vacuum%nvac == 1 ) alphm(1,2) = alphm(1,1)
     622           0 :     do ivac = 1, vacuum%nvac
     623             :       zeta(1:nzdhprec) = quotz(-1:-nzdhprec:-1,1) * alphm(1,mod(ivac,2)+1) &
     624             :                        + quotz(1:nzdhprec,1) * delta(1:nzdhprec,ivac) &
     625           0 :                        + sinhz(1:nzdhprec,1) * epsilon(1:nzdhprec,ivac) 
     626           0 :       VVz(1:nzdhprec,ivac) = VVz(1:nzdhprec,ivac) + vcons2(1) * zeta(1:nzdhprec)
     627             :     end do
     628             : 
     629             : 
     630           0 :   end subroutine VYukawaFilmVacuumVariant2
     631             : 
     632             : 
     633             : 
     634           0 :   subroutine VYukawaFilmInterstitialVariant2( &
     635             :                 stars, vacuum, cell, sym, input, &
     636           0 :                 psq, VVxy, VVz, alphm, dh_prec, coshdh, &
     637           0 :                 VIq )
     638             : 
     639             :     ! main parts:
     640             :     ! 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)
     641             :     ! 2. part: Add the contribution from the vacuum charge density to the interstitial potential, which had already been computed earlier for the vacuum potential
     642             :  
     643             :     ! 4. part: Compute the coefficients V^I(q_xy,q_z) from the function V^I(q_xy,z) by a 1D Fourier transform:
     644             :     !          V^I(q_xy,z) = sum_{q_z} V^I(q_xy,q_z) * exp( ImagUnit * q_z * z )
     645             :     ! In order to be able to match the interstitial and vacuum potentials smoothly at the interface, the Fourier transform is done
     646             :     ! in a slightly larger region. -> 3. part
     647             :     ! 3. part: Interpolate the vacuum potential in a small region surrounding the slab
     648             : 
     649             :     use m_ExpSave
     650             :     use m_constants
     651             :     use m_types
     652             :     use m_cfft
     653             :     implicit none
     654             : 
     655             :     type(t_stars),  intent(in)  :: stars
     656             :     type(t_vacuum), intent(in)  :: vacuum
     657             :     type(t_cell),   intent(in)  :: cell
     658             :     type(t_sym),    intent(in)  :: sym
     659             :     type(t_input),  intent(in)  :: input
     660             :     complex,        intent(in)  :: psq(stars%ng3)
     661             :     complex,        intent(in)  :: VVxy(vacuum%nmzxyd,2:stars%ng2,2)
     662             :     real,           intent(in)  :: VVz(vacuum%nmzd,2)
     663             :     complex,        intent(in)  :: alphm(stars%ng2,2)
     664             :     real,           intent(in)  :: dh_prec
     665             :     real,           intent(in)  :: coshdh(stars%ng2)
     666             : 
     667             :     complex,        intent(out) :: VIq(stars%ng3)
     668             : 
     669             :     real                                               :: partitioning, rz, qz, q, qxy_numerics
     670             :     integer                                            :: irec2, irec3, iz, jz, ivac, iqz, jvac
     671             :     integer                                            :: nfft, nzmax, nzmin, nzdh, nLower, nUpper
     672           0 :     complex, allocatable                               :: VIz(:,:), eta(:,:)
     673           0 :     complex, allocatable                               :: expzqz(:,:)
     674           0 :     complex, dimension(-stars%mx3:stars%mx3,stars%ng2) :: VIqz, c_ph
     675           0 :     complex, dimension(-stars%mx3:stars%mx3,stars%ng2) :: qquot, qquottrigp, qquottrigm
     676           0 :     complex, dimension(stars%ng3)                      :: vcons3
     677           0 :     real, dimension(3*stars%mx3,stars%ng2)             :: VIzReal, VIzImag
     678           0 :     real, allocatable                                  :: quotz(:,:), sinhz(:,:)
     679           0 :     real, allocatable                                  :: z(:)
     680           0 :     real, dimension(stars%ng2)                         :: g_damped, vcons2
     681             :     
     682             : 
     683             :     ! DEFINITIONS / ALLOCATIONS / INITIALISATIONS
     684             : 
     685             :     ! grid points z_i:
     686           0 :     qxy_numerics = sqrt( ( 100 / dh_prec ) ** 2 - input%preconditioning_param ** 2 )
     687           0 :     nfft = 3 * stars%mx3                                  ! number of grid points for Fourier transform
     688           0 :     partitioning = 1. / real( nfft )
     689           0 :     nzmax = nfft / 2                                      ! index of maximal z below D~/2 
     690           0 :     nzmin = nzmax - nfft + 1                              ! index of minimal z above -D~/2
     691           0 :     nzdh = ceiling( cell%z1 / cell%amat(3,3) * nfft ) - 1 ! index of maximal z below D/2
     692           0 :     allocate( z(nzmin:nzmax) )
     693             :     ! definition of z_i:        ! indexing: z_0 = 0; positive indices for z > 0; negative indices for z < 0
     694           0 :     do iz = nzmin, nzmax
     695           0 :       z(iz) = cell%amat(3,3) * iz * partitioning
     696             :     end do
     697             :     ! other variables:
     698           0 :     allocate( VIz(nzmin:nzmax,stars%ng2) )
     699           0 :     allocate( eta(-nzdh:nzdh,stars%ng2) )
     700           0 :     allocate( sinhz(-nzdh:nzdh,stars%ng2) )
     701           0 :     allocate( quotz(-nzdh:nzdh,stars%ng2) )
     702           0 :     allocate( expzqz(-nzdh:nzdh,-stars%mx3:stars%mx3) )
     703           0 :     do irec2 = 1, stars%ng2
     704           0 :       g_damped(irec2) = sqrt( stars%sk2(irec2) ** 2 + input%preconditioning_param ** 2 )
     705           0 :       vcons2(irec2) = fpi_const / g_damped(irec2)
     706           0 :       if( stars%sk2(irec2) < qxy_numerics ) then ! numerics ok
     707           0 :         do iz = -nzdh, nzdh
     708             :           quotz( iz,irec2) = ( cosh( g_damped(irec2) * z(iz) ) / cosh( g_damped(irec2) * dh_prec ) &
     709           0 :                              + sinh( g_damped(irec2) * z(iz) ) / sinh( g_damped(irec2) * dh_prec ) ) / 2
     710             :         end do
     711             :       else ! numerical treatment necessary
     712           0 :         do iz = -nzdh, nzdh
     713           0 :           quotz( iz,irec2) = exp_save( g_damped(irec2) * (  z(iz) - dh_prec ) )
     714             :         end do
     715             :       end if
     716           0 :       do iz = -nzdh, nzdh
     717           0 :         sinhz(iz,irec2) = sinh( g_damped(irec2) * ( dh_prec - z(iz) ) )
     718           0 :         do iqz = -stars%mx3, stars%mx3
     719           0 :           expzqz(iz,iqz) = exp( ImagUnit * iqz * cell%bmat(3,3) * z(iz) )
     720             :         end do
     721             :       end do
     722           0 :       do iqz = -stars%mx3, stars%mx3
     723           0 :         qquot(iqz,irec2) = ImagUnit * iqz * cell%bmat(3,3) / g_damped(irec2)
     724           0 :         qquottrigp(iqz,irec2) = coshdh(irec2) + sinhz(nzdh,irec2) * qquot(iqz,irec2)
     725           0 :         qquottrigm(iqz,irec2) = coshdh(irec2) - sinhz(nzdh,irec2) * qquot(iqz,irec2)
     726             :       end do
     727             :     end do
     728           0 :     do irec3 = 1, stars%ng3
     729           0 :       vcons3(irec3) = fpi_const * psq(irec3) / ( stars%sk3(irec3) ** 2 + input%preconditioning_param ** 2 )
     730             :     end do
     731           0 :     VIz = (0.,0.)
     732           0 :     VIq = (0.,0.)
     733             : 
     734             : 
     735             :     ! CONTRIBUTION FROM THE INTERSTITIAL CHARGE DENSITY
     736             :     
     737             :     ! compute V^I(q_xy,z) as a function of q_xy and z
     738           0 :     do irec2 = 1, stars%ng2  
     739           0 :       do iz = -nzdh, nzdh
     740           0 :         do iqz = -stars%mx3, stars%mx3
     741           0 :           irec3 = stars%ig(stars%kv2(1,irec2),stars%kv2(2,irec2),iqz)
     742           0 :           if ( irec3 /= 0 ) then ! use only stars within the g_max sphere
     743           0 :             c_ph(iqz,irec2) = stars%rgphs(stars%kv2(1,irec2),stars%kv2(2,irec2),iqz)
     744           0 :             qz = iqz * cell%bmat(3,3)
     745             :             VIz(iz,irec2) = VIz(iz,irec2) &
     746             :                           + vcons3(irec3) * c_ph(iqz,irec2) * &
     747             :                             ( expzqz(iz,iqz) &
     748             :                             - qquottrigp(iqz,irec2) * expzqz( nzdh,iqz) * quotz( iz,irec2) &
     749           0 :                             - qquottrigm(iqz,irec2) * expzqz(-nzdh,iqz) * quotz(-iz,irec2) )
     750             :           end if
     751             :         enddo
     752             :       end do
     753             :     ! irec2 loop continues
     754             : 
     755             : 
     756             :     ! CONTRIBUTION FROM THE VACUUM CHARGE DENSITY
     757             : 
     758             :     ! irec2 loop continues
     759             :       eta(-nzdh:nzdh,irec2) = quotz(nzdh:-nzdh:-1,irec2) * alphm(irec2,2) &
     760           0 :                             + quotz(-nzdh:nzdh,   irec2) * alphm(irec2,1)
     761           0 :       VIz(-nzdh:nzdh,irec2) = VIz(-nzdh:nzdh,irec2) + vcons2(irec2) * eta(-nzdh:nzdh,irec2)
     762             :     ! irec2 loop continues
     763             :    
     764             :     
     765             :     ! INTERPOLATION IN VACUUM REGION
     766             : 
     767             :     ! use Lagrange polynomials of order 3 to interpolate the vacuum potential outside I
     768             :     ! q, q-1 and q-2 (scaled with +/-1 or +/-0.5) are the factors of the Lagrange basis polynomials
     769             :     ! irec2 loop continues
     770           0 :       do ivac = 1, 2
     771           0 :         select case( ivac )
     772             :           case( 1 )
     773           0 :             nUpper = nzmax; nLower =  nzdh + 1; jvac = 1
     774             :           case( 2 )
     775           0 :             nLower = nzmin; nUpper = -nzdh - 1; jvac = 2; if ( sym%invs .or. sym%zrfs ) jvac = 1
     776             :         end select
     777           0 :         do iz = nLower, nUpper
     778           0 :           rz = ( abs( z(iz) ) - cell%z1 ) / vacuum%delz + 1.0
     779           0 :           jz = rz      ! index of maximal vacuum grid point below z_i
     780           0 :           q = rz - jz  ! factor in Lagrange basis polynomials
     781           0 :           if ( irec2 == 1 ) then
     782             :             VIz(iz,irec2) = 0.5     * ( q - 1. ) * ( q - 2. ) * VVz(jz,  jvac) &
     783             :                           -       q *              ( q - 2. ) * VVz(jz+1,jvac) &
     784           0 :                           + 0.5 * q * ( q - 1. )              * VVz(jz+2,jvac)
     785           0 :           else if ( jz + 2 <= vacuum%nmzxy ) then
     786             :             VIz(iz,irec2) = 0.5 *     ( q - 1. ) * ( q - 2. ) * VVxy(jz,  irec2,jvac) &
     787             :                           -       q              * ( q - 2. ) * VVxy(jz+1,irec2,jvac) &
     788           0 :                           + 0.5 * q * ( q - 1. )              * VVxy(jz+2,irec2,jvac)
     789           0 :             if ( ( sym%invs .and. .not. sym%zrfs ) .and. ivac == 2 ) VIz(iz,irec2) = conjg( VIz(iz,irec2) )
     790             :           end if
     791             :         end do
     792             :       end do
     793             :     end do ! irec2
     794             : 
     795             : 
     796             :     ! 1D FOURIER TRANSFORM TO FIND THE COEFFICIENTS V^I(q_xy,q_z)
     797             : 
     798             :     ! change the indexing for the subroutine cfft, and split real and imaginary parts
     799           0 :     VIzReal(1:nzmax+1,:) =  real( VIz(0:nzmax,:) ); VIzReal(nzmax+2:nfft,:) =  real( VIz(nzmin:-1,:) )
     800           0 :     VIzImag(1:nzmax+1,:) = aimag( VIz(0:nzmax,:) ); VIzImag(nzmax+2:nfft,:) = aimag( VIz(nzmin:-1,:) )
     801             :  
     802             :     ! V^I(q_xy,z) = sum_{q_z} V^I(q_xy,q_z) * exp( ImagUnit * q_z * z )
     803           0 :     do irec2 = 1, stars%ng2
     804           0 :       call cfft( VIzReal(:,irec2), VIzImag(:,irec2), nfft, nfft, nfft, -1 )
     805             :     ! irec2 loop continues
     806             : 
     807             :     ! reorder
     808           0 :       VIqz(0,irec2) = cmplx( VIzReal(1,irec2), VIzImag(1,irec2) )
     809           0 :       do iqz = 1, stars%mx3
     810           0 :         VIqz( iqz,irec2) = cmplx( VIzReal(iqz+1,     irec2), VIzImag(iqz+1,     irec2) )
     811           0 :         VIqz(-iqz,irec2) = cmplx( VIzReal(nfft+1-iqz,irec2), VIzImag(nfft+1-iqz,irec2) )
     812             :       end do
     813             :    
     814             :     ! add the computed components to V^I(q_xy,q_z):
     815           0 :       do iqz= -stars%mx3, stars%mx3
     816           0 :         irec3 = stars%ig(stars%kv2(1,irec2),stars%kv2(2,irec2),iqz)
     817           0 :         if ( irec3 /= 0 ) VIq(irec3) = VIq(irec3) + VIqz(iqz,irec2) * partitioning / ( stars%nstr(irec3) / stars%nstr2(irec2) )
     818             :       end do
     819             :     end do
     820             : 
     821             : 
     822           0 :   end subroutine VYukawaFilmInterstitialVariant2
     823             : 
     824             : 
     825             : 
     826           0 :   subroutine VYukawaModify( stars, vacuum, cell, sym, input, mpi, atoms, sphhar, oneD, noco, den, &
     827             :                             VYukawa )
     828             : 
     829             :     ! This subroutine adds a potential to the previously computed Yukawa
     830             :     ! potential to ensure charge neutrality.
     831             :     ! The added potential itself is a solution to the modified Helmholtz
     832             :     ! equation with constant right-hand side, where the constant is chosen such
     833             :     ! that charge neutrality is obtained.
     834             :     ! The charge is distributed only over the film region, and we therefore 
     835             :     ! solve the differential equation subject to a boundary condition on the 
     836             :     ! film surface, in contrast to the basic Yukawa potential above. 
     837             : 
     838             :     use m_constants
     839             :     use m_types
     840             :     use m_vmts
     841             :     use m_constants
     842             :     use m_cdntot
     843             :     use m_cfft
     844             :     implicit none
     845             : 
     846             :     type(t_stars),      intent(in)    :: stars
     847             :     type(t_vacuum),     intent(in)    :: vacuum
     848             :     type(t_cell),       intent(in)    :: cell
     849             :     type(t_sym),        intent(in)    :: sym
     850             :     type(t_input),      intent(in)    :: input
     851             :     type(t_mpi),        intent(in)    :: mpi
     852             :     type(t_atoms),      intent(in)    :: atoms
     853             :     type(t_sphhar),     intent(in)    :: sphhar
     854             :     type(t_oneD),       intent(in)    :: oneD
     855             :     type(t_noco),       intent(in)    :: noco
     856             :     type(t_potden),     intent(inout) :: den
     857             :     type(t_potden),     intent(inout) :: VYukawa
     858             : 
     859             :     integer                           :: n, lh, irec3, iz, iqz, nfft, nzmax, nzmin, nzdh
     860             :     real                              :: q0, qhat, qbar, ldh, partitioning, dh
     861           0 :     real                              :: q(input%jspins), qis(input%jspins), qmt(atoms%ntype,input%jspins),qvac(2,input%jspins), qtot, qistot
     862           0 :     complex                           :: psq(stars%ng3)
     863           0 :     type(t_potden)                    :: VYukawaModification
     864           0 :     real, allocatable                 :: z(:)
     865           0 :     real                              :: VIzReal(3*stars%mx3), VIzImag(3*stars%mx3)
     866           0 :     complex, allocatable              :: VIz(:)
     867           0 :     complex                           :: VIqz(-stars%mx3:stars%mx3)
     868             : 
     869             : 
     870             :     ! DEFINITIONS / ALLOCATIONS / INITIALISATIONS
     871             : 
     872             :     ! constants:
     873           0 :     dh = cell%z1   ! half the width of the film
     874             :     ! indexing of grid points z_i:
     875           0 :     nfft = 3 * stars%mx3   ! number of grid points for Fourier transform
     876           0 :     partitioning = 1. / real( nfft )
     877           0 :     nzmax = nfft / 2           ! index of maximal z below D~/2 
     878           0 :     nzmin = nzmax - nfft + 1   ! index of minimal z above -D~/2
     879           0 :     nzdh = ceiling( dh / cell%amat(3,3) * nfft ) - 1   ! index of maximal z below D/2
     880           0 :     allocate( z(nzmin:nzmax) )
     881             :     ! definition of grid points z_i:
     882             :     ! indexing: z_0 = 0; positive indices for z > 0; negative indices for z < 0
     883           0 :     do iz = nzmin, nzmax
     884           0 :       z(iz) = cell%amat(3,3) * iz * partitioning
     885             :     end do
     886             : 
     887             : 
     888             :     ! INTEGRATION OF THE PREVIOUSLY COMPUTED YUKAWA POTENTIAL
     889             : 
     890             :     ! initialise VYukawaModification with in-going VYukawa and prepare for integration
     891           0 :     call VYukawaModification%init( stars, atoms, sphhar, vacuum, noco, input%jspins, 4 )
     892           0 :     call VYukawaModification%copyPotDen( VYukawa )
     893           0 :     do n = 1, atoms%ntype
     894           0 :       do lh = 0, sphhar%nlhd    
     895           0 :         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
     896             :       end do
     897             :     end do
     898             : 
     899             :     ! integrate the potential over the film region
     900           0 :     call integrate_cdn( stars, atoms, sym, vacuum, input, cell, oneD, VYukawaModification, q, qis, qmt, qvac, qtot, qistot  )
     901           0 :     q0 = qtot / cell%area
     902           0 :     ldh = input%preconditioning_param * dh
     903           0 :     qhat = ( q0 / ( 2 * dh ) ) / ( sinh(ldh) / ( ldh * cosh( ldh ) ) - 1 )
     904           0 :     qbar = input%preconditioning_param ** 2 / fpi_const *  qhat
     905             : 
     906             : 
     907             :     ! SET UP CONSTANT CHARGE DENSITY
     908             : 
     909             :     ! instead of den%pw(1,1) = qbar we directly set the pseudo charge density
     910           0 :     den%mt = 0; den%pw = 0; den%vacxy = 0; den%vacz = 0    
     911           0 :     do n = 1, atoms%ntype
     912           0 :       den%mt(1:atoms%jri(n),0,n,1) = sfp_const * qbar * atoms%rmsh(1:atoms%jri(n),n) ** 2
     913             :     end do
     914           0 :     psq = cmplx(0.0,0.0); psq(1) = qbar
     915             : 
     916             : 
     917             :     ! CALCULATE THE INTERSTITIAL POTENTIAL AS A FUNCTION OF z
     918             : 
     919             :     ! initialise and calculate out-going modification potential; reuse VYukawaModification
     920           0 :     VYukawaModification%mt = 0; VYukawaModification%pw = 0; VYukawaModification%vacxy = 0; VYukawaModification%vacz = 0
     921             : 
     922           0 :     allocate( VIz(nzmin:nzmax) )
     923           0 :     VIz = (0.,0.)
     924           0 :     do iz = -nzdh+1, nzdh-1
     925           0 :       VIz(iz) = qhat * ( 1 - cosh( input%preconditioning_param * z(iz) ) / cosh( ldh ) ) 
     926             :     end do
     927             : 
     928             : 
     929             :     ! 1D FOURIER TRANSFORM TO FIND THE 3D-FOURIER COEFFICIENTS
     930             : 
     931             :     ! change the indexing for the subroutine cfft, and split real and imaginary parts
     932           0 :     VIzReal(1:nzmax+1) =  real( VIz(0:nzmax) ); VIzReal(nzmax+2:nfft) = real( VIz(nzmin:-1) )
     933           0 :     VIzImag(1:nzmax+1) = aimag( VIz(0:nzmax) ); VIzImag(nzmax+2:nfft) = aimag( VIz(nzmin:-1) )
     934             : 
     935           0 :     call cfft( VIzReal(:), VIzImag(:), nfft, nfft, nfft, -1 )
     936             : 
     937             :     ! reorder
     938           0 :     VIqz = 0
     939           0 :     VIqz(0) = cmplx( VIzReal(1), VIzImag(1) )
     940           0 :     do iqz = 1, stars%mx3
     941           0 :       VIqz( iqz) = cmplx( VIzReal(iqz+1     ), VIzImag(iqz+1     ) )
     942           0 :       VIqz(-iqz) = cmplx( VIzReal(nfft+1-iqz), VIzImag(nfft+1-iqz) )
     943             :     end do
     944             : 
     945             :     ! add the computed components
     946           0 :     do iqz= -stars%mx3, stars%mx3
     947           0 :       irec3 = stars%ig(stars%kv2(1,1),stars%kv2(2,1),iqz)
     948           0 :       if ( irec3 /= 0 ) VYukawaModification%pw(irec3,1) = VYukawaModification%pw(irec3,1) + VIqz(iqz) * partitioning / ( stars%nstr(irec3) / stars%nstr2(1) )
     949             :     end do
     950             : 
     951             : 
     952             :     ! MUFFIN-TIN POTENTIAL
     953             : 
     954             :     call Vmts( input, mpi, stars, sphhar, atoms, sym, cell, oneD, &
     955             :                VYukawaModification%pw(:,1), den%mt(:,0:,:,1), VYukawaModification%potdenType, &
     956           0 :                VYukawaModification%mt(:,0:,:,1) )
     957             : 
     958             : 
     959             :     ! APPLYING THE MODIFICATION TO THE YUKAWA POTENTIAL
     960             : 
     961           0 :     call VYukawa%AddPotDen( VYukawa, VYukawaModification )
     962             : 
     963             : 
     964           0 :   end subroutine VYukawaModify
     965             : 
     966             : 
     967             : 
     968             : end module m_VYukawaFilm

Generated by: LCOV version 1.13