LCOV - code coverage report
Current view: top level - vgen - vvac.f90 (source / functions) Hit Total Coverage
Test: combined.info Lines: 30 63 47.6 %
Date: 2019-09-08 04:53:50 Functions: 1 1 100.0 %

          Line data    Source code
       1             : module m_vvac
       2             :   ! ****************************************************************
       3             :   ! calculates the g(2-dim)=0 part of the vacuum coulomb potential *
       4             :   ! for general symmetry.          c.l.fu, r.podloucky             *
       5             :   ! ****************************************************************
       6             :   contains
       7             : 
       8           5 :   subroutine vvac( vacuum, stars, cell, sym, input, field, psq, rht, vz, rhobar, sig1dh, vz1dh )
       9             : 
      10             :     use m_constants
      11             :     use m_qsf
      12             :     use m_types
      13             :     implicit none
      14             : 
      15             :     type(t_vacuum), intent(in)    :: vacuum
      16             :     type(t_stars),  intent(in)    :: stars
      17             :     type(t_cell),   intent(in)    :: cell
      18             :     type(t_sym),    intent(in)    :: sym
      19             : 
      20             :     complex,        intent(out)   :: rhobar
      21             :     real,           intent(out)   :: sig1dh, vz1dh
      22             :     type(t_input),  intent(in)    :: input
      23             :     type(t_field),  intent(inout) :: field ! efield is modified here
      24             : 
      25             :     real,           intent(in)    :: rht(vacuum%nmzd,2) 
      26             :     complex,        intent(in)    :: psq(stars%ng3)
      27             :     real,           intent(out)   :: vz(vacuum%nmzd,2)
      28             : 
      29             :     complex                       :: sumq, vcons
      30             :     real                          :: bj0, bj1, dh, qzh, sigmaa(2)
      31             :     integer                       :: ig3, imz, ivac, ncsh
      32           5 :     real                          :: f(vacuum%nmzd), sig(vacuum%nmzd), vtemp(vacuum%nmzd)
      33             :     intrinsic cos, sin
      34             : 
      35             : 
      36           5 :     vz(:,1:vacuum%nvac) = 0.0 ! initialize potential
      37             : 
      38             :     ! obtain mesh point (ncsh) of charge sheet for external electric field
      39           5 :     ncsh = field%efield%zsigma / vacuum%delz + 1.01
      40           5 :     sigmaa(1) = ( field%efield%sigma + field%efield%sig_b(1) ) / cell%area
      41           5 :     sigmaa(2) = ( field%efield%sigma + field%efield%sig_b(2) ) / cell%area
      42             : 
      43             :     ! g=0 vacuum potential due to neutral charge density
      44             :     ! inside slab and zero charge density outside
      45             : 
      46           5 :     vcons = - 2. * fpi_const * ImagUnit
      47           5 :     dh = cell%z1
      48           5 :     rhobar = - psq(1)
      49           5 :     sumq = cmplx( 0.0, 0.0 )
      50             : 
      51       16294 :     do ig3 = 2, stars%ng3
      52       16294 :       if (stars%ig2(ig3) == 1) then           ! select g_|| = 0
      53         125 :         qzh = stars%kv3(3,ig3) * cell%bmat(3,3) * dh
      54         125 :         bj0 = sin(qzh) / qzh
      55         125 :         rhobar = rhobar - psq(ig3) * bj0 * stars%nstr(ig3)
      56         125 :         if ( .not.( sym%zrfs .or. sym%invs ) ) then
      57           0 :           bj1 = ( sin(qzh) - qzh * cos(qzh) ) / ( qzh * qzh )
      58           0 :           sumq = sumq + bj1 * psq(ig3) * dh * dh
      59             :         endif
      60             :       endif
      61             :     enddo
      62             : 
      63           5 :     ivac = 2                        ! lower (ivac=2) vacuum
      64           5 :     if ( vacuum%nvac == 2 ) then
      65           0 :       vz(1:vacuum%nmz,ivac) = vcons * sumq
      66             :     endif
      67             : 
      68             :     ! g=0 vacuum potential due to
      69             :     ! negative of rhobar + vacuum (g=0) charge ----> v2(z)
      70             : 
      71           5 :     ivac = 1 ! upper vacuum
      72             : 
      73             :     ! =================== dirichlet ==============================
      74           5 :     if ( field%efield%dirichlet ) then
      75           0 :       vz(ncsh+1:vacuum%nmz,ivac) = field%efield%sig_b(1)
      76           0 :       call qsf( vacuum%delz, rht(1,ivac), sig, ncsh, 1 )
      77           0 :       sig(1:ncsh) = sig(ncsh) - sig(1:ncsh)
      78           0 :       call qsf( vacuum%delz, sig, vtemp, ncsh, 1 )
      79           0 :       do imz = 1, ncsh
      80           0 :         vz(imz,ivac) = - fpi_const * ( vtemp(ncsh) - vtemp(imz) ) + field%efield%sig_b(1)
      81             :       enddo
      82           0 :       sig1dh = sig(1)
      83           0 :       vz1dh = vz(1,ivac)   ! potential on vacuum boundary
      84           0 :       if ( vacuum%nvac == 1 ) return
      85             : 
      86           0 :       ivac = 2     ! lower vacuum
      87           0 :       call qsf( vacuum%delz, rht(1,ivac), sig, ncsh, 1 )
      88           0 :       f(1:ncsh) = sig(1:ncsh) - rhobar*vacuum%dvac + sig1dh
      89           0 :       call qsf( vacuum%delz, f, vtemp, ncsh, 1 )
      90           0 :       do imz = 1,ncsh
      91           0 :         vz(imz,ivac) = - fpi_const * ( vtemp(imz) + sig1dh * vacuum%dvac - rhobar * vacuum%dvac * vacuum%dvac / 2. ) + vz(imz,ivac) + vz1dh
      92             :       end do
      93             : 
      94             :       ! force matching on the other side
      95           0 :       field%efield%vslope = ( field%efield%sig_b(2) - vz(ncsh,1) ) / ( 2 * vacuum%delz * ( ncsh + 1 ) + vacuum%dvac )
      96           0 :       ivac = 1
      97           0 :       do imz = 1, ncsh
      98           0 :         vz(imz,ivac) = vz(imz,ivac) + field%efield%vslope * vacuum%delz * ( ncsh - imz + 1 )
      99             :       end do
     100             :       ivac = 2
     101           0 :       do imz = 1, ncsh
     102           0 :         vz(imz,ivac) = vz(imz,ivac) + field%efield%vslope * ( vacuum%dvac + vacuum%delz * imz + vacuum%delz * ncsh )
     103             :       end do
     104           0 :       vz(ncsh+1:vacuum%nmz,ivac) = field%efield%sig_b(2)
     105             : 
     106             :       ! =================== neumann ==============================
     107             :     else ! neumann
     108             : 
     109           5 :       call qsf( vacuum%delz, rht(1,ivac), sig, vacuum%nmz, 1 )
     110           5 :       sig1dh = sig(vacuum%nmz) - sigmaa(1)  ! need to include contribution from electric field
     111        1255 :       sig(1:vacuum%nmz) = sig(vacuum%nmz) - sig(1:vacuum%nmz)
     112           5 :       call qsf( vacuum%delz, sig, vtemp, vacuum%nmz, 1 )
     113             :       ! external electric field contribution 
     114         510 :       do imz = 1, ncsh
     115         510 :         vz(imz,ivac) = - fpi_const * ( vtemp(vacuum%nmz) - vtemp(imz) ) + vz(imz,ivac) - fpi_const * ( imz - ncsh ) * vacuum%delz * sigmaa(1)
     116             :       enddo
     117         750 :       do imz = ncsh + 1, vacuum%nmz
     118         750 :         vz(imz,ivac) = - fpi_const * ( vtemp(vacuum%nmz) - vtemp(imz) ) + vz(imz,ivac)
     119             :       enddo
     120           5 :       vz1dh = vz(1,ivac)   ! potential on vacuum boundary
     121           5 :       if ( vacuum%nvac == 1 ) return
     122             : 
     123           0 :       ivac = 2 ! lower vacuum
     124           0 :       call qsf( vacuum%delz, rht(1,ivac), sig, vacuum%nmz, 1 )
     125           0 :       f(1:vacuum%nmz) = sig(1:vacuum%nmz) - rhobar * vacuum%dvac + sig1dh
     126           0 :       call qsf( vacuum%delz, f, vtemp, vacuum%nmz, 1 )
     127             : 
     128             :       ! external electric field contribution
     129           0 :       do imz = 1, ncsh
     130           0 :         vz(imz,ivac) = - fpi_const * ( vtemp(imz) + sig1dh * vacuum%dvac - rhobar * vacuum%dvac * vacuum%dvac / 2. ) + vz1dh + vz(imz,ivac)
     131             :       enddo
     132           0 :       do imz = ncsh + 1, vacuum%nmz
     133           0 :         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)
     134             :       enddo
     135             : 
     136             :     end if ! dirichlet (vs. neumann)
     137             : 
     138             :   end subroutine vvac
     139             : 
     140             : end module m_vvac

Generated by: LCOV version 1.13