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

Generated by: LCOV version 2.0-1