vvac.f90 5.45 KB
Newer Older
1 2 3 4 5 6 7
module m_vvac
  ! ****************************************************************
  ! calculates the g(2-dim)=0 part of the vacuum coulomb potential *
  ! for general symmetry.          c.l.fu, r.podloucky             *
  ! ****************************************************************
  contains

Daniel Wortmann's avatar
Daniel Wortmann committed
8
  subroutine vvac( vacuum, stars, cell, sym, input, field, psq, rht, vz, rhobar, sig1dh, vz1dh ,vslope)
9 10 11 12 13 14 15 16 17 18 19 20 21 22

    use m_constants
    use m_qsf
    use m_types
    implicit none

    type(t_vacuum), intent(in)    :: vacuum
    type(t_stars),  intent(in)    :: stars
    type(t_cell),   intent(in)    :: cell
    type(t_sym),    intent(in)    :: sym

    complex,        intent(out)   :: rhobar
    real,           intent(out)   :: sig1dh, vz1dh
    type(t_input),  intent(in)    :: input
Daniel Wortmann's avatar
Daniel Wortmann committed
23 24
    type(t_field),  intent(in)    :: field
    COMPLEX,INTENT(OUT)           :: vslope
25

Daniel Wortmann's avatar
Daniel Wortmann committed
26
    real,           intent(in)    :: rht(vacuum%nmzd,2)
27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83
    complex,        intent(in)    :: psq(stars%ng3)
    real,           intent(out)   :: vz(vacuum%nmzd,2)

    complex                       :: sumq, vcons
    real                          :: bj0, bj1, dh, qzh, sigmaa(2)
    integer                       :: ig3, imz, ivac, ncsh
    real                          :: f(vacuum%nmzd), sig(vacuum%nmzd), vtemp(vacuum%nmzd)
    intrinsic cos, sin


    vz(:,1:vacuum%nvac) = 0.0 ! initialize potential

    ! obtain mesh point (ncsh) of charge sheet for external electric field
    ncsh = field%efield%zsigma / vacuum%delz + 1.01
    sigmaa(1) = ( field%efield%sigma + field%efield%sig_b(1) ) / cell%area
    sigmaa(2) = ( field%efield%sigma + field%efield%sig_b(2) ) / cell%area

    ! g=0 vacuum potential due to neutral charge density
    ! inside slab and zero charge density outside

    vcons = - 2. * fpi_const * ImagUnit
    dh = cell%z1
    rhobar = - psq(1)
    sumq = cmplx( 0.0, 0.0 )

    do ig3 = 2, stars%ng3
      if (stars%ig2(ig3) == 1) then           ! select g_|| = 0
        qzh = stars%kv3(3,ig3) * cell%bmat(3,3) * dh
        bj0 = sin(qzh) / qzh
        rhobar = rhobar - psq(ig3) * bj0 * stars%nstr(ig3)
        if ( .not.( sym%zrfs .or. sym%invs ) ) then
          bj1 = ( sin(qzh) - qzh * cos(qzh) ) / ( qzh * qzh )
          sumq = sumq + bj1 * psq(ig3) * dh * dh
        endif
      endif
    enddo

    ivac = 2                        ! lower (ivac=2) vacuum
    if ( vacuum%nvac == 2 ) then
      vz(1:vacuum%nmz,ivac) = vcons * sumq
    endif

    ! g=0 vacuum potential due to
    ! negative of rhobar + vacuum (g=0) charge ----> v2(z)

    ivac = 1 ! upper vacuum

    ! =================== dirichlet ==============================
    if ( field%efield%dirichlet ) then
      vz(ncsh+1:vacuum%nmz,ivac) = field%efield%sig_b(1)
      call qsf( vacuum%delz, rht(1,ivac), sig, ncsh, 1 )
      sig(1:ncsh) = sig(ncsh) - sig(1:ncsh)
      call qsf( vacuum%delz, sig, vtemp, ncsh, 1 )
      do imz = 1, ncsh
        vz(imz,ivac) = - fpi_const * ( vtemp(ncsh) - vtemp(imz) ) + field%efield%sig_b(1)
      enddo
      sig1dh = sig(1)
84
      vz1dh = vz(1,ivac)   ! potential on vacuum boundary
85
      if ( vacuum%nvac == 1 ) return
86 87

      ivac = 2     ! lower vacuum
88 89 90 91 92 93 94 95
      call qsf( vacuum%delz, rht(1,ivac), sig, ncsh, 1 )
      f(1:ncsh) = sig(1:ncsh) - rhobar*vacuum%dvac + sig1dh
      call qsf( vacuum%delz, f, vtemp, ncsh, 1 )
      do imz = 1,ncsh
        vz(imz,ivac) = - fpi_const * ( vtemp(imz) + sig1dh * vacuum%dvac - rhobar * vacuum%dvac * vacuum%dvac / 2. ) + vz(imz,ivac) + vz1dh
      end do

      ! force matching on the other side
Daniel Wortmann's avatar
Daniel Wortmann committed
96
      vslope = ( field%efield%sig_b(2) - vz(ncsh,1) ) / ( 2 * vacuum%delz * ( ncsh + 1 ) + vacuum%dvac )
97 98
      ivac = 1
      do imz = 1, ncsh
Daniel Wortmann's avatar
Daniel Wortmann committed
99
        vz(imz,ivac) = vz(imz,ivac) + vslope * vacuum%delz * ( ncsh - imz + 1 )
100 101 102
      end do
      ivac = 2
      do imz = 1, ncsh
Daniel Wortmann's avatar
Daniel Wortmann committed
103
        vz(imz,ivac) = vz(imz,ivac) + vslope * ( vacuum%dvac + vacuum%delz * imz + vacuum%delz * ncsh )
104 105 106 107 108 109 110 111 112 113
      end do
      vz(ncsh+1:vacuum%nmz,ivac) = field%efield%sig_b(2)

      ! =================== neumann ==============================
    else ! neumann

      call qsf( vacuum%delz, rht(1,ivac), sig, vacuum%nmz, 1 )
      sig1dh = sig(vacuum%nmz) - sigmaa(1)  ! need to include contribution from electric field
      sig(1:vacuum%nmz) = sig(vacuum%nmz) - sig(1:vacuum%nmz)
      call qsf( vacuum%delz, sig, vtemp, vacuum%nmz, 1 )
Daniel Wortmann's avatar
Daniel Wortmann committed
114
      ! external electric field contribution
115 116 117 118 119 120 121 122
      do imz = 1, ncsh
        vz(imz,ivac) = - fpi_const * ( vtemp(vacuum%nmz) - vtemp(imz) ) + vz(imz,ivac) - fpi_const * ( imz - ncsh ) * vacuum%delz * sigmaa(1)
      enddo
      do imz = ncsh + 1, vacuum%nmz
        vz(imz,ivac) = - fpi_const * ( vtemp(vacuum%nmz) - vtemp(imz) ) + vz(imz,ivac)
      enddo
      vz1dh = vz(1,ivac)   ! potential on vacuum boundary
      if ( vacuum%nvac == 1 ) return
123

124 125 126 127
      ivac = 2 ! lower vacuum
      call qsf( vacuum%delz, rht(1,ivac), sig, vacuum%nmz, 1 )
      f(1:vacuum%nmz) = sig(1:vacuum%nmz) - rhobar * vacuum%dvac + sig1dh
      call qsf( vacuum%delz, f, vtemp, vacuum%nmz, 1 )
128

129 130 131 132 133 134 135
      ! external electric field contribution
      do imz = 1, ncsh
        vz(imz,ivac) = - fpi_const * ( vtemp(imz) + sig1dh * vacuum%dvac - rhobar * vacuum%dvac * vacuum%dvac / 2. ) + vz1dh + vz(imz,ivac)
      enddo
      do imz = ncsh + 1, vacuum%nmz
        vz(imz,ivac) = - fpi_const * ( vtemp(imz) + sig1dh * vacuum%dvac - rhobar * vacuum%dvac * vacuum%dvac / 2. ) + vz1dh + vz(imz,ivac) + fpi_const * ( imz - ncsh ) * vacuum%delz * sigmaa(2)
      enddo
136

137
    end if ! dirichlet (vs. neumann)
138

139
  end subroutine vvac
140

141
end module m_vvac