VYukawaFilm.f90 15.2 KB
Newer Older
1 2
module m_VYukawaFilm

3
  ! computation of the film-case Yukawa potential for the preconditioning of the charge density residual
4 5 6 7

  contains 


8

9 10 11 12 13 14
  subroutine VYukawaFilm( stars, vacuum, cell, sym, input, mpi, atoms, sphhar, dimension, oneD, den, &
                          VYukawa )

    use m_constants
    use m_types
    use m_psqpw
15
    use m_vmts
16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32
    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)
33
    complex                           :: alphm(stars%ng2,2)
34 35 36 37 38 39 40 41 42


    ! 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

43 44 45
    call VYukawaFilmVacuum( stars, vacuum, cell, sym, input, &
                            psq, den%vacxy(:,:,:,1), den%vacz(:,:,1), &
                            VYukawa%vacxy, VYukawa%vacz, alphm )
46 47 48 49


    ! INTERSTITIAL POTENTIAL

50 51 52
    call VYukawaFilmInterstitial( stars, vacuum, cell, sym, input, &
                                  psq, VYukawa%vacxy, VYukawa%vacz, alphm, &
                                  VYukawa%pw(:,1) )
53 54 55 56


    ! MUFFIN-TIN POTENTIAL

57
    call Vmts( input, mpi, stars, sphhar, atoms, sym, cell, oneD, VYukawa%pw(:,1), den%mt(:,0:,:,1), VYukawa%potdenType, VYukawa%mt(:,0:,:,1) )
58 59 60 61 62


  end subroutine VYukawaFilm


63

64 65 66 67
  subroutine VYukawaFilmVacuum( stars, vacuum, cell, sym, input, &
                                psq, rhtxy, rht, &
                                VVxy, VVz, alphm )

68 69 70
    ! 1. part: compute the contribution from the interstitial charge density to the vacuum potential as a function of q_xy and z
    ! 2. part: compute the contribution from the vacuum charge density to the vacuum potential as a function of q_xy and z by numerical integration

71 72 73 74 75 76 77 78 79 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 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130
    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)
    integer                     :: iz, irec2, irec3, ivac, kz
    complex                     :: fa(vacuum%nmzxyd,2:stars%ng2), fb(vacuum%nmzxyd,2:stars%ng2)
    complex                     :: alpha(vacuum%nmzxyd,2:stars%ng2,2), beta(vacuum%nmzxyd,2:stars%ng2,2), gamma(vacuum%nmzxyd,2:stars%ng2)
    real                        :: ga(vacuum%nmzd), gb(vacuum%nmzd)
    real                        :: delta(vacuum%nmzd,2), epsilon(vacuum%nmzd,2), zeta(vacuum%nmzd)


    ! definitions / initialisations:
    do iz = 1, vacuum%nmz
      z(iz) = ( iz - 1 ) * vacuum%delz ! z_i starting at 0
    end do
    do irec2 = 1, stars%ng2
      g_damped(irec2) = sqrt( stars%sk2(irec2) ** 2 + input%preconditioning_param ** 2 )
      vcons(irec2) = tpi_const / g_damped(irec2)
      do iz = 1, vacuum%nmz
        exp_m(iz,irec2) = exp_save( - g_damped(irec2) * z(iz) ) ! redefined in contribution of vacuum charge density part
        exp_p(iz,irec2) = exp_save(   g_damped(irec2) * z(iz) )
      end do
      expDhg(irec2) = exp_save( - cell%z1 * g_damped(irec2) )
      expDg(irec2)  = expDhg(irec2) ** 2
    end do
    sum_qz = (0.,0.)
    VVxy = (0.,0.)
    VVz = 0.
   

    ! CONTRIBUTION FROM THE INTERSTITIAL CHARGE DENSITY

    do irec2 = 1, stars%ng2
      do ivac = 1, vacuum%nvac
        sign = 3. - 2. * ivac
        do kz = -stars%mx3, stars%mx3
          irec3 = stars%ig(stars%kv2(1,irec2),stars%kv2(2,irec2),kz)
131
          ! use only stars within the g_max sphere -> stars outside the sphere have per definition index ig3n = 0
132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154
          if ( irec3 /= 0 ) then
            c_ph(kz,irec2) = stars%rgphs(stars%kv2(1,irec2),stars%kv2(2,irec2),kz)
            qz = kz * cell%bmat(3,3)
            signIqz = sign * ImagUnit * qz
            sum_qz(ivac,irec2) = sum_qz(ivac,irec2) + c_ph(kz,irec2) * psq(irec2) * ( exp( signIqz * cell%z1 ) - exp( - signIqz * cell%z1 ) * expDg(irec2) ) / ( signIqz + g_damped(irec2) )
          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
      exp_m(1:vacuum%nmzxy,irec2) = exp_m(1:vacuum%nmzxy,irec2) * expDhg(irec2) ! z_i can now be thought of as starting from D/2 -> expDhg is the correction
      exp_p(1:vacuum%nmzxy,irec2) = exp_p(1:vacuum%nmzxy,irec2) * expDhg(irec2)
      do ivac = 1, vacuum%nvac
155
        ! integrands:
156
        fa(1:vacuum%nmzxy,irec2) = rhtxy(1:vacuum%nmzxy,irec2-1,ivac) * exp_m(1:vacuum%nmzxy,irec2)
157 158 159 160 161 162 163 164
        fb(1:vacuum%nmzxy,irec2) = rhtxy(1:vacuum%nmzxy,irec2-1,ivac) * exp_p(1:vacuum%nmzxy,irec2)
        ! integrals:
        ! alpha(z,q_xy,ivac) = int_z^infty rho(z',q_xy,ivac) exp(-sqrt(q_xy**2+prec_param**2)*z') dz'
        ! beta (z,q_xy,ivac) = int_{D/2}^z rho(z',q_xy,ivac) exp(+sqrt(q_xy**2+prec_param**2)*z') dz'
        ! where for z < 0 the lower vacuum charge density (ivac=2) is defined by rho(q_xy,z,ivac=2) := rho(q_xy,-z,ivac=2)
        call intgz1Reverse( fa(:,irec2), vacuum%delz, vacuum%nmzxy, alpha(:,irec2,ivac), .true. ) 
        call qsfComplex( vacuum%delz, fb(:,irec2), beta(:,irec2,ivac), vacuum%nmzxy, 1 )
        ! alphm(q_xy,ivac) = alpha(D/2,q_xy,ivac) --- these integrals are also needed for the interstitial potential
165 166
        alphm(irec2,ivac) = alpha(1,irec2,ivac)
      end do
167
      if ( vacuum%nvac == 1 ) then
168 169 170 171 172 173 174 175 176 177
        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. )
178
        VVxy(1:vacuum%nmzxy,irec2,ivac) = VVxy(1:vacuum%nmzxy,irec2,ivac) + vcons(irec2) * gamma(1:vacuum%nmzxy,irec2) 
179 180 181 182 183 184 185 186
      end do
    end do

    ! case irec2 = 1:
    exp_m(1:vacuum%nmz,1) = exp_m(1:vacuum%nmz,1) * expDhg(1)
    exp_p(1:vacuum%nmz,1) = exp_p(1:vacuum%nmz,1) * expDhg(1)
    do ivac = 1, vacuum%nvac
      ga(1:vacuum%nmz) = rht(1:vacuum%nmz,ivac) * exp_m(1:vacuum%nmz,1)
187 188 189 190
      gb(1:vacuum%nmz) = rht(1:vacuum%nmz,ivac) * exp_p(1:vacuum%nmz,1)
      call intgz1Reverse( ga(:), vacuum%delz, vacuum%nmz, delta(:,ivac), .true. ) ! integrals 
      call qsf( vacuum%delz, gb(:), epsilon(:,ivac), vacuum%nmz, 1 )
      alphm(1,ivac) = delta(1,ivac)
191 192 193 194 195 196
    end do
    if ( vacuum%nvac == 1 ) alphm(1,2) = alphm(1,1) ! is real, thus no conjg as in the irec2 > 1 case
    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.
197
      VVz(1:vacuum%nmz,ivac) = VVz(1:vacuum%nmz,ivac) + vcons(1) * zeta(1:vacuum%nmz) 
198 199 200 201 202 203
    end do


  end subroutine VYukawaFilmVacuum


204

205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337
  subroutine VYukawaFilmInterstitial( stars, vacuum, cell, sym, input, &
                                      psq, VVxy, VVz, alphm, &
                                      VIq )

    ! 1. part: compute the contribution from the interstitial charge density to the interstitial potential (largest part) as a function of q_xy and z
    ! 2. part: add the contribution from the vacuum charge density to the interstitial potential, which had already been computed earlier for the vacuum potential
    ! 3. part: compute the coefficients V^I(q_xy,q_z) from the function V^I(q_xy,z) by a 1D Fourier transform: V^I(q_xy,z) = sum_{q_z} V^I(q_xy,q_z) * exp( ImagUnit * q_z * z )

    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)

    real                        :: partitioning, rz, qz, fit, q
    integer                     :: irec2, irec3, iz, jz, ivac, kz, ifft
    complex                     :: VIz(3*stars%mx3,stars%ng2), sum_qz(3*stars%mx3,stars%ng2), eta(3*stars%mx3,stars%ng2)
    complex                     :: VIqz(-stars%mx3:stars%mx3,stars%ng2), c_ph(-stars%mx3:stars%mx3,stars%ng2)
    complex                     :: vcons1(stars%ng3)
    real                        :: VIzReal(3*stars%mx3,stars%ng2), VIzImag(3*stars%mx3,stars%ng2), exp_m(3*stars%mx3,stars%ng2), exp_p(3*stars%mx3,stars%ng2)
    real                        :: z(3*stars%mx3)
    real                        :: g_damped(stars%ng2), vcons2(stars%ng2)


    ! definitions / initialisations
    ifft = 3 * stars%mx3
    partitioning = 1. / real( ifft )
    do iz = 1, ifft
      z(iz) = cell%amat(3,3) * ( iz - 1 ) * partitioning               ! z_i are equidistantly distributed along ( 0, cell%amat(3,3) ]
      if( z(iz) > cell%amat(3,3) / 2. ) z(iz) = z(iz) - cell%amat(3,3) ! z_i are equidistantly distributed along ( -cell%amat(3,3) / 2, cell%amat(3,3) / 2 ]
    end do
    do irec2 = 1, stars%ng2
      g_damped(irec2) = sqrt( stars%sk2(irec2) ** 2 + input%preconditioning_param ** 2 )
      vcons2(irec2) = -1. / ( 2. * g_damped(irec2) )
      do iz = 1, ifft
        exp_m(iz,irec2) = exp_save( - g_damped(irec2) * ( cell%z1 + z(iz) ) )
        exp_p(iz,irec2) = exp_save( - g_damped(irec2) * ( cell%z1 - z(iz) ) )
      end do
    end do
    do irec3 = 1, stars%ng3
      vcons1(irec3) = fpi_const * psq(irec3) / ( stars%sk3(irec3) ** 2 + input%preconditioning_param ** 2 )
    end do
    sum_qz = (0.,0.)
    VIz = (0.,0.)
    VIq = (0.,0.)


    ! CONTRIBUTION FROM THE INTERSTITIAL CHARGE DENSITY

    ! compute V^I(q_xy,z) as a function of q_xy and z
    do irec2 = 1, stars%ng2
      do iz = 1, ifft

        ! in the transition region ( D/2, D~/2 ), smooth out the potential numerically
        if ( abs( z(iz) ) >= cell%z1 ) then
          ivac = 1 ! upper vacuum
          if ( z(iz) < 0.0 .and. .not. ( sym%invs .or. sym%zrfs ) ) ivac = 2 ! lower vacuum
          rz = ( abs( z(iz) ) - cell%z1 ) / vacuum%delz + 1.0 ! numbering the grid points in the vacuum region
          jz = rz
          q = rz - jz
          if ( irec2 == 1 ) then
            fit = 0.5     * ( q - 1. ) * ( q - 2. ) * VVz(jz,  ivac) &
                -       q *              ( q - 2. ) * VVz(jz+1,ivac) &
                + 0.5 * q * ( q - 1. )              * VVz(jz+2,ivac)
            VIz(iz,irec2) = cmplx( fit, 0.0 )
          else if ( jz + 2 <= vacuum%nmzxy ) then
            VIz(iz,irec2) = 0.5 *     ( q - 1. ) * ( q - 2. ) * VVxy(jz,  irec2-1,ivac) &
                          -       q              * ( q - 2. ) * VVxy(jz+1,irec2-1,ivac) &
                          + 0.5 * q * ( q - 1. )              * VVxy(jz+2,irec2-1,ivac)
            if ( ( sym%invs .and. .not. sym%zrfs ) .and. z(iz) < 0 ) VIz(iz,irec2) = conjg( VIz(iz,irec2) )
          end if

        ! z in ( -D/2, D/2 )
        else
          do kz = -stars%mx3, stars%mx3
            irec3 = stars%ig(stars%kv2(1,irec2),stars%kv2(2,irec2),kz)
            if ( irec3 /= 0 ) then ! use only stars within the g_max sphere
              c_ph(kz,irec2) = stars%rgphs(stars%kv2(1,irec2),stars%kv2(2,irec2),kz)
              qz = kz * cell%bmat(3,3)
              sum_qz(iz,irec2) = sum_qz(iz,irec2) + exp( ImagUnit * qz * z(iz) ) &
                               + vcons2(irec2) * ( ( g_damped(irec2) + ImagUnit * qz ) * exp_p(iz,irec2) * exp(   ImagUnit * qz * cell%z1 ) &
                                                 + ( g_damped(irec2) - ImagUnit * qz ) * exp_m(iz,irec2) * exp( - ImagUnit * qz * cell%z1 ) )
              sum_qz(iz,irec2) = vcons1(irec3) * c_ph(kz,irec2) * sum_qz(iz,irec2)
            end if
          enddo
          VIz(iz,irec2) = VIz(iz,irec2) + sum_qz(iz,irec2)
        end if
      end do


      ! CONTRIBUTION FROM THE VACUUM CHARGE DENSITY

      eta(:,irec2) = exp_m(:,irec2) * alphm(irec2,2) + exp_p(:,irec2) * alphm(irec2,1)
      where ( 2.0 * eta(:,irec2) == eta(:,irec2) ) eta(:,irec2) = cmplx( 0.0, 0.0 )
      VIz(:,irec2) = VIz(:,irec2) + tpi_const / g_damped(irec2) * eta(:,irec2)


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

      ! V^I(q_xy,z) = sum_{q_z} V^I(q_xy,q_z) * exp( ImagUnit * q_z * z )
      VIzReal =  real( VIz )
      VIzImag = aimag( VIz )
      call cfft( VIzReal(:,irec2), VIzImag(:,irec2), ifft, ifft, ifft, -1 )
      ! reorder:
      VIqz(0,irec2) = cmplx( VIzReal(1,irec2), VIzImag(1,irec2) )
      do kz = 1, stars%mx3
        VIqz( kz,irec2) = cmplx( VIzReal(kz+1,     irec2), VIzImag(kz+1,     irec2) )
        VIqz(-kz,irec2) = cmplx( VIzReal(ifft+1-kz,irec2), VIzImag(ifft+1-kz,irec2) )
      end do
      ! add the computed components to V^I(q_xy,q_z):
      do kz= -stars%mx3, stars%mx3
        irec3 = stars%ig(stars%kv2(1,irec2),stars%kv2(2,irec2),kz)
        VIq(irec3) = VIq(irec3) + VIqz(kz,irec2) / partitioning * stars%nstr(irec3) / stars%nstr2(irec2)
      end do
    end do


  end subroutine VYukawaFilmInterstitial



338
end module m_VYukawaFilm