VYukawaFilm.f90 17.2 KB
Newer Older
1 2
module m_VYukawaFilm

3 4 5 6 7 8 9 10 11 12 13
  ! Computation of the film-case Yukawa potential for the preconditioning of the charge density residual

  ! 1. part: The pseudo-charge density is used for the vacuum and interstitial potentials.
  ! 4. part: The MT potential is calculated at the end.
  ! These two parts do not change compared to the bulk case.

  ! 2. and 3. part: The vacuum and interstitial potentials are calculated by integration of the product of the charge density and
  ! the 1D Green function with respect to z (from  minus infinity to infinity). The integrals for both potentials are split in 4:
  ! the seperating points are the boundaries of the slab and the point where the integration variable z' equals the point z where we
  ! want to evaluate the potential for. For the contribution from the vacuum charge density we do a numerical integration. For the
  ! contribution from the slab charge density we have analytical expressions of the integrals.
14 15 16 17

  contains 


18

19 20 21 22 23 24
  subroutine VYukawaFilm( stars, vacuum, cell, sym, input, mpi, atoms, sphhar, dimension, oneD, den, &
                          VYukawa )

    use m_constants
    use m_types
    use m_psqpw
25
    use m_vmts
26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42
    implicit none

    type(t_stars),      intent(in)    :: stars
    type(t_vacuum),     intent(in)    :: vacuum
    type(t_cell),       intent(in)    :: cell
    type(t_sym),        intent(in)    :: sym
    type(t_input),      intent(in)    :: input
    type(t_mpi),        intent(in)    :: mpi
    type(t_atoms),      intent(in)    :: atoms 
    type(t_sphhar),     intent(in)    :: sphhar
    type(t_dimension),  intent(in)    :: dimension
    type(t_oneD),       intent(in)    :: oneD
    type(t_potden),     intent(in)    :: den

    type(t_potden),     intent(inout) :: VYukawa

    complex                           :: psq(stars%ng3)
43
    complex                           :: alphm(stars%ng2,2)
44 45 46 47 48 49 50 51 52


    ! PSEUDO-CHARGE DENSITY

    call psqpw( mpi, atoms, sphhar, stars, vacuum, dimension, cell, input, sym, oneD, den%pw(:,1), den%mt(:,:,:,1), den%vacz(:,:,1), .false., VYukawa%potdenType, psq )


    ! VACUUM POTENTIAL

53 54 55
    call VYukawaFilmVacuum( stars, vacuum, cell, sym, input, &
                            psq, den%vacxy(:,:,:,1), den%vacz(:,:,1), &
                            VYukawa%vacxy, VYukawa%vacz, alphm )
56 57 58 59


    ! INTERSTITIAL POTENTIAL

60 61 62
    call VYukawaFilmInterstitial( stars, vacuum, cell, sym, input, &
                                  psq, VYukawa%vacxy, VYukawa%vacz, alphm, &
                                  VYukawa%pw(:,1) )
63 64 65 66


    ! MUFFIN-TIN POTENTIAL

67
    call Vmts( input, mpi, stars, sphhar, atoms, sym, cell, oneD, VYukawa%pw(:,1), den%mt(:,0:,:,1), VYukawa%potdenType, VYukawa%mt(:,0:,:,1) )
68 69 70 71 72


  end subroutine VYukawaFilm


73

74 75 76 77
  subroutine VYukawaFilmVacuum( stars, vacuum, cell, sym, input, &
                                psq, rhtxy, rht, &
                                VVxy, VVz, alphm )

78 79
    ! 1. part: Compute the contribution from the interstitial charge density to the vacuum potential as a function of q_xy and z (analytic expression for integral)
    ! 2. part: Compute the contribution from the vacuum charge density to the vacuum potential as a function of q_xy and z by numerical integration
80

81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107
    use m_ExpSave
    use m_constants
    use m_types
    use m_intgr, only: intgz1Reverse
    use m_qsf
    implicit none

    type(t_stars),  intent(in)  :: stars
    type(t_vacuum), intent(in)  :: vacuum
    type(t_cell),   intent(in)  :: cell
    type(t_sym),    intent(in)  :: sym
    type(t_input),  intent(in)  :: input
    complex,        intent(in)  :: psq(stars%ng3)
    complex,        intent(in)  :: rhtxy(vacuum%nmzxyd,stars%ng2-1,2)
    real,           intent(in)  :: rht(vacuum%nmzd,2) 

    complex,        intent(out) :: VVxy(vacuum%nmzxyd,2:stars%ng2,2) ! this is the qxy /= 0 part of the vacuum potential
    real,           intent(out) :: VVz(vacuum%nmzd,2)                ! this is the qxy  = 0 part of the vacuum potential
    complex,        intent(out) :: alphm(stars%ng2,2)                ! these are the integrals in upper and lower vacuum, now including the first star---integral for star ig2 is in alphm(ig2,ivac) 

    complex                     :: sum_qz(2,stars%ng2)
    complex                     :: c_ph(-stars%mx3:stars%mx3,stars%ng2)
    complex                     :: signIqz
    real                        :: g_damped(stars%ng2), qz, sign, vcons(stars%ng2) 
    real                        :: exp_m(vacuum%nmzd,stars%ng2), exp_p(vacuum%nmzd,stars%ng2)
    real                        :: expDhg(stars%ng2), expDg(stars%ng2)
    real                        :: z(vacuum%nmzd)
108
    integer                     :: iz, irec2, irec3, ivac, iqz
109 110 111 112 113 114
    complex                     :: fa(vacuum%nmzxyd,2:stars%ng2), fb(vacuum%nmzxyd,2:stars%ng2)
    complex                     :: alpha(vacuum%nmzxyd,2:stars%ng2,2), beta(vacuum%nmzxyd,2:stars%ng2,2), gamma(vacuum%nmzxyd,2:stars%ng2)
    real                        :: ga(vacuum%nmzd), gb(vacuum%nmzd)
    real                        :: delta(vacuum%nmzd,2), epsilon(vacuum%nmzd,2), zeta(vacuum%nmzd)


115 116
    ! DEFINITIONS / ALLOCATIONS / INITIALISATIONS    
    
117
    do iz = 1, vacuum%nmz
118
      z(iz) = ( iz - 1 ) * vacuum%delz
119 120 121 122 123
    end do
    do irec2 = 1, stars%ng2
      g_damped(irec2) = sqrt( stars%sk2(irec2) ** 2 + input%preconditioning_param ** 2 )
      vcons(irec2) = tpi_const / g_damped(irec2)
      do iz = 1, vacuum%nmz
124
        exp_m(iz,irec2) = exp_save( - g_damped(irec2) * z(iz) )
125 126 127
        exp_p(iz,irec2) = exp_save(   g_damped(irec2) * z(iz) )
      end do
      expDhg(irec2) = exp_save( - cell%z1 * g_damped(irec2) )
128
      expDg(irec2)  = exp_save( -2 * cell%z1 * g_damped(irec2) )
129 130 131 132 133 134 135 136 137 138 139
    end do
    sum_qz = (0.,0.)
    VVxy = (0.,0.)
    VVz = 0.
   

    ! CONTRIBUTION FROM THE INTERSTITIAL CHARGE DENSITY

    do irec2 = 1, stars%ng2
      do ivac = 1, vacuum%nvac
        sign = 3. - 2. * ivac
140 141
        do iqz = -stars%mx3, stars%mx3
          irec3 = stars%ig(stars%kv2(1,irec2),stars%kv2(2,irec2),iqz)
142
          ! use only stars within the g_max sphere -> stars outside the sphere have per definition index ig3n = 0
143
          if ( irec3 /= 0 ) then
144 145 146 147
            c_ph(iqz,irec2) = stars%rgphs(stars%kv2(1,irec2),stars%kv2(2,irec2),iqz)
            qz = iqz * cell%bmat(3,3)
            signIqz = sign * ImagUnit * qz 
            sum_qz(ivac,irec2) = sum_qz(ivac,irec2) + c_ph(iqz,irec2) * psq(irec3) * ( exp( signIqz * cell%z1 ) - exp( - signIqz * cell%z1 ) * expDg(irec2) ) / ( signIqz + g_damped(irec2) )
148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163
          endif
        enddo
        if( irec2 /= 1 ) then
          VVxy(1:vacuum%nmzxy,irec2,ivac) = vcons(irec2) * sum_qz(ivac,irec2) * exp_m(1:vacuum%nmzxy,irec2)
        else
          VVz(1:vacuum%nmz,ivac)          = vcons(1)     * sum_qz(ivac,1)     * exp_m(1:vacuum%nmz,1)
        end if
      enddo
    enddo


    ! CONTRIBUTION FROM THE VACUUM CHARGE DENSITY

    ! case irec2 > 1:
    do irec2 = 2, stars%ng2
      do ivac = 1, vacuum%nvac
164
        ! integrands:
165 166
        fa(1:vacuum%nmzxy,irec2) = rhtxy(1:vacuum%nmzxy,irec2-1,ivac) * exp_m(1:vacuum%nmzxy,irec2) * expDhg(irec2)
        fb(1:vacuum%nmzxy,irec2) = rhtxy(1:vacuum%nmzxy,irec2-1,ivac) * exp_p(1:vacuum%nmzxy,irec2) * expDhg(irec2)
167 168 169 170 171 172 173
        ! integrals:
        ! alpha(z,q_xy,ivac) = int_z^infty rho(z',q_xy,ivac) exp(-sqrt(q_xy**2+prec_param**2)*z') dz'
        ! beta (z,q_xy,ivac) = int_{D/2}^z rho(z',q_xy,ivac) exp(+sqrt(q_xy**2+prec_param**2)*z') dz'
        ! where for z < 0 the lower vacuum charge density (ivac=2) is defined by rho(q_xy,z,ivac=2) := rho(q_xy,-z,ivac=2)
        call intgz1Reverse( fa(:,irec2), vacuum%delz, vacuum%nmzxy, alpha(:,irec2,ivac), .true. ) 
        call qsfComplex( vacuum%delz, fb(:,irec2), beta(:,irec2,ivac), vacuum%nmzxy, 1 )
        ! alphm(q_xy,ivac) = alpha(D/2,q_xy,ivac) --- these integrals are also needed for the interstitial potential
174 175
        alphm(irec2,ivac) = alpha(1,irec2,ivac)
      end do
176
      if ( vacuum%nvac == 1 ) then
177 178 179 180 181 182 183 184 185 186
        if ( sym%invs ) then
          alphm(irec2,2) = conjg( alphm(irec2,1) )
        else
          alphm(irec2,2) = alphm(irec2,1)
        end if
      end if
      do ivac = 1, vacuum%nvac
        gamma(1:vacuum%nmzxy,irec2) = exp_m(1:vacuum%nmzxy,irec2) * ( alphm(irec2,mod(ivac,2)+1) + beta(1:vacuum%nmzxy,irec2,ivac) ) &
                                    + exp_p(1:vacuum%nmzxy,irec2) *                               alpha(1:vacuum%nmzxy,irec2,ivac) ! mod(ivac,2)+1 outputs the other ivac value 
        where ( 2. * gamma(:,irec2) == gamma(:,irec2) ) gamma(:,irec2) = cmplx( 0., 0. )
187
        VVxy(1:vacuum%nmzxy,irec2,ivac) = VVxy(1:vacuum%nmzxy,irec2,ivac) + vcons(irec2) * gamma(1:vacuum%nmzxy,irec2) * expDhg(irec2) 
188 189 190 191 192
      end do
    end do

    ! case irec2 = 1:
    do ivac = 1, vacuum%nvac
193 194
      ga(1:vacuum%nmz) = rht(1:vacuum%nmz,ivac) * exp_m(1:vacuum%nmz,1) * expDhg(1)
      gb(1:vacuum%nmz) = rht(1:vacuum%nmz,ivac) * exp_p(1:vacuum%nmz,1) * expDhg(1)
195 196 197
      call intgz1Reverse( ga(:), vacuum%delz, vacuum%nmz, delta(:,ivac), .true. ) ! integrals 
      call qsf( vacuum%delz, gb(:), epsilon(:,ivac), vacuum%nmz, 1 )
      alphm(1,ivac) = delta(1,ivac)
198
    end do
199
    if ( vacuum%nvac == 1 ) alphm(1,2) = alphm(1,1)
200 201 202 203
    do ivac = 1, vacuum%nvac
      zeta(1:vacuum%nmz) = exp_m(1:vacuum%nmz,1) * ( alphm(1,mod(ivac,2)+1) + epsilon(1:vacuum%nmz,ivac) ) &
                         + exp_p(1:vacuum%nmz,1) *                              delta(1:vacuum%nmz,ivac) 
      where ( 2. * zeta == zeta ) zeta = 0.
204
      VVz(1:vacuum%nmz,ivac) = VVz(1:vacuum%nmz,ivac) + vcons(1) * zeta(1:vacuum%nmz) * expDhg(1)
205 206 207 208 209 210
    end do


  end subroutine VYukawaFilmVacuum


211

212 213 214 215
  subroutine VYukawaFilmInterstitial( stars, vacuum, cell, sym, input, &
                                      psq, VVxy, VVz, alphm, &
                                      VIq )

216 217 218 219 220 221 222 223 224
    ! main parts:
    ! 1. part: Compute the contribution from the interstitial charge density to the interstitial potential as a function of q_xy and z (analytic expression for integral)
    ! 2. part: Add the contribution from the vacuum charge density to the interstitial potential, which had already been computed earlier for the vacuum potential
 
    ! 4. part: Compute the coefficients V^I(q_xy,q_z) from the function V^I(q_xy,z) by a 1D Fourier transform:
    !          V^I(q_xy,z) = sum_{q_z} V^I(q_xy,q_z) * exp( ImagUnit * q_z * z )
    ! In order to be able to match the interstitial and vacuum potentials smoothly at the interface, the Fourier transform is done
    ! in a slightly larger region. -> 3. part
    ! 3. part: Interpolate the vacuum potential in a small region surrounding the slab
225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243

    use m_ExpSave
    use m_constants
    use m_types
    use m_cfft
    implicit none

    type(t_stars),  intent(in)  :: stars
    type(t_vacuum), intent(in)  :: vacuum
    type(t_cell),   intent(in)  :: cell
    type(t_sym),    intent(in)  :: sym
    type(t_input),  intent(in)  :: input
    complex,        intent(in)  :: psq(stars%ng3)
    complex,        intent(in)  :: VVxy(vacuum%nmzxyd,2:stars%ng2,2)
    real,           intent(in)  :: VVz(vacuum%nmzd,2)
    complex,        intent(in)  :: alphm(stars%ng2,2)

    complex,        intent(out) :: VIq(stars%ng3)

244 245 246
    real                        :: partitioning, rz, qz, q
    integer                     :: irec2, irec3, iz, jz, ivac, iqz, nfft, nzmax, nzmin, nzdh, nLower, nUpper, jvac
    complex, allocatable        :: VIz(:,:), eta(:,:)
247 248
    complex                     :: VIqz(-stars%mx3:stars%mx3,stars%ng2), c_ph(-stars%mx3:stars%mx3,stars%ng2)
    complex                     :: vcons1(stars%ng3)
249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266
    real                        :: VIzReal(3*stars%mx3,stars%ng2), VIzImag(3*stars%mx3,stars%ng2)
    real, allocatable           :: exp_m(:,:), exp_p(:,:)
    real, allocatable           :: z(:)
    real                        :: g_damped(stars%ng2), vcons2(stars%ng2), expDhg(stars%ng2)
    

    ! DEFINITIONS / ALLOCATIONS / INITIALISATIONS

    ! grid points z_i:
    nfft = 3 * stars%mx3                                  ! number of grid points for Fourier transform
    partitioning = 1. / real( nfft )
    nzmax = nfft / 2                                      ! index of maximal z below D~/2 
    nzmin = nzmax - nfft + 1                              ! index of minimal z above -D~/2
    nzdh = ceiling( cell%z1 / cell%amat(3,3) * nfft ) - 1 ! index of maximal z below D/2
    allocate( z(nzmin:nzmax) )
    ! definition of z_i:        ! indexing: z_0 = 0; positive indices for z > 0; negative indices for z < 0
    do iz = nzmin, nzmax
      z(iz) = cell%amat(3,3) * iz * partitioning
267
    end do
268 269 270
    ! other variables:
    allocate( VIz(nzmin:nzmax,stars%ng2), eta(nzmin:nzmax,stars%ng2) )
    allocate( exp_m(nzmin:nzmax,stars%ng2), exp_p(nzmin:nzmax,stars%ng2) )
271 272 273
    do irec2 = 1, stars%ng2
      g_damped(irec2) = sqrt( stars%sk2(irec2) ** 2 + input%preconditioning_param ** 2 )
      vcons2(irec2) = -1. / ( 2. * g_damped(irec2) )
274 275 276
      do iz = nzmin, nzmax
        exp_m(iz,irec2) = exp_save( - g_damped(irec2) * z(iz) )
        exp_p(iz,irec2) = exp_save(   g_damped(irec2) * z(iz) )
277
      end do
278
      expDhg(irec2) = exp_save( - cell%z1 * g_damped(irec2) )
279 280 281 282 283 284 285 286 287
    end do
    do irec3 = 1, stars%ng3
      vcons1(irec3) = fpi_const * psq(irec3) / ( stars%sk3(irec3) ** 2 + input%preconditioning_param ** 2 )
    end do
    VIz = (0.,0.)
    VIq = (0.,0.)


    ! CONTRIBUTION FROM THE INTERSTITIAL CHARGE DENSITY
288
    
289
    ! compute V^I(q_xy,z) as a function of q_xy and z
290 291 292 293 294 295 296 297 298 299 300 301 302
    do irec2 = 1, stars%ng2  
      do iz = -nzdh, nzdh
        do iqz = -stars%mx3, stars%mx3
          irec3 = stars%ig(stars%kv2(1,irec2),stars%kv2(2,irec2),iqz)
          if ( irec3 /= 0 ) then ! use only stars within the g_max sphere
            c_ph(iqz,irec2) = stars%rgphs(stars%kv2(1,irec2),stars%kv2(2,irec2),iqz)
            qz = iqz * cell%bmat(3,3)
            VIz(iz,irec2) = VIz(iz,irec2) &
                          + vcons1(irec3) * c_ph(iqz,irec2) * &
                            ( exp( ImagUnit * qz * z(iz) ) &
                            + vcons2(irec2) * expDhg(irec2) * &
                              ( ( g_damped(irec2) + ImagUnit * qz ) * exp_p(iz,irec2) * exp(   ImagUnit * qz * cell%z1 ) &
                              + ( g_damped(irec2) - ImagUnit * qz ) * exp_m(iz,irec2) * exp( - ImagUnit * qz * cell%z1 ) ) )
303
          end if
304
        enddo
305
      end do
306
    ! irec2 loop continues
307 308


309
    ! CONTRIBUTION FROM THE VACUUM CHARGE DENSITY
310

311 312
    ! irec2 loop continues
      eta(-nzdh:nzdh,irec2) = exp_m(-nzdh:nzdh,irec2) * alphm(irec2,2) + exp_p(-nzdh:nzdh,irec2) * alphm(irec2,1)
313
      where ( 2.0 * eta(:,irec2) == eta(:,irec2) ) eta(:,irec2) = cmplx( 0.0, 0.0 )
314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347
      VIz(-nzdh:nzdh,irec2) = VIz(-nzdh:nzdh,irec2) + tpi_const / g_damped(irec2) * eta(-nzdh:nzdh,irec2)
    ! irec2 loop continues
    
    
    ! INTERPOLATION IN VACUUM REGION

    ! use Lagrange polynomials of order 3 to interpolate the vacuum potential outside I
    ! q, q-1 and q-2 (scaled with +/-1 or +/-0.5) are the factors of the Lagrange basis polynomials
    ! irec2 loop continues
      do ivac = 1, 2
        select case( ivac )
          case( 1 )
            nUpper = nzmax; nLower =  nzdh + 1; jvac = 1
          case( 2 )
            nUpper = nzmin; nLower = -nzdh - 1; jvac = 2; if ( sym%invs .or. sym%zrfs ) jvac = 1
        end select
        do iz = nLower, nUpper
          rz = ( abs( z(iz) ) - cell%z1 ) / vacuum%delz + 1.0
          jz = rz      ! index of maximal vacuum grid point below z_i
          q = rz - jz  ! factor in Lagrange basis polynomials
          if ( irec2 == 1 ) then
            VIz(iz,irec2) = 0.5     * ( q - 1. ) * ( q - 2. ) * VVz(jz,  jvac) &
                          -       q *              ( q - 2. ) * VVz(jz+1,jvac) &
                          + 0.5 * q * ( q - 1. )              * VVz(jz+2,jvac)
          else if ( jz + 2 <= vacuum%nmzxy ) then
            VIz(iz,irec2) = 0.5 *     ( q - 1. ) * ( q - 2. ) * VVxy(jz,  irec2,jvac) &
                          -       q              * ( q - 2. ) * VVxy(jz+1,irec2,jvac) &
                          + 0.5 * q * ( q - 1. )              * VVxy(jz+2,irec2,jvac)
            if ( ( sym%invs .and. .not. sym%zrfs ) .and. ivac == 2 ) VIz(iz,irec2) = conjg( VIz(iz,irec2) )
          end if
        end do
      end do
    end do ! irec2

348

349
    ! 1D FOURIER TRANSFORM TO FIND THE COEFFICIENTS V^I(q_xy,q_z)
350

351 352 353 354 355 356 357 358
    ! change the indexing for the subroutine cfft, and split real and imaginary parts
    VIzReal(1:nzmax+1,:) =  real( VIz(0:nzmax,:) ); VIzReal(nzmax+2:nfft,:) =  real( VIz(nzmin:-1,:) )
    VIzImag(1:nzmax+1,:) = aimag( VIz(0:nzmax,:) ); VIzImag(nzmax+2:nfft,:) = aimag( VIz(nzmin:-1,:) )
 
    ! V^I(q_xy,z) = sum_{q_z} V^I(q_xy,q_z) * exp( ImagUnit * q_z * z )
    do irec2 = 1, stars%ng2
      call cfft( VIzReal(:,irec2), VIzImag(:,irec2), nfft, nfft, nfft, -1 )
    ! irec2 loop continues
359

360
    ! reorder
361
      VIqz(0,irec2) = cmplx( VIzReal(1,irec2), VIzImag(1,irec2) )
362 363 364
      do iqz = 1, stars%mx3
        VIqz( iqz,irec2) = cmplx( VIzReal(iqz+1,     irec2), VIzImag(iqz+1,     irec2) )
        VIqz(-iqz,irec2) = cmplx( VIzReal(nfft+1-iqz,irec2), VIzImag(nfft+1-iqz,irec2) )
365
      end do
366 367 368 369 370
   
    ! add the computed components to V^I(q_xy,q_z):
      do iqz= -stars%mx3, stars%mx3
        irec3 = stars%ig(stars%kv2(1,irec2),stars%kv2(2,irec2),iqz)
        if ( irec3 /= 0 ) VIq(irec3) = VIq(irec3) + VIqz(iqz,irec2) * partitioning / ( stars%nstr(irec3) / stars%nstr2(irec2) )
371 372 373 374 375 376 377
      end do
    end do

  end subroutine VYukawaFilmInterstitial



378
end module m_VYukawaFilm