vvac.f90 6.05 KB
Newer Older
1 2 3 4 5 6 7 8
      MODULE m_vvac
!     ****************************************************************
!     calculates the g(2-dim)=0 part of the vacuum coulomb potential *
!     for general symmetry.          c.l.fu, r.podloucky             *
!     ****************************************************************
      CONTAINS
      SUBROUTINE vvac(&
     &                vacuum,stars,&
9
     &                cell,sym,input,field,&
10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
     &                psq,rht,&
     &                vz,rhobar,sig1dh,vz1dh)

      USE m_constants
      USE m_qsf
      USE m_types
      IMPLICIT NONE
!     ..
!     .. Scalar Arguments ..
      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
26 27
      TYPE(t_input), INTENT(IN) :: input
      TYPE(t_field),INTENT(INOUT) :: field !efield is modified here
28 29
!     ..
!     .. Array Arguments ..
30
      REAL,    INTENT (IN) :: rht(vacuum%nmzd,2) 
31
      COMPLEX, INTENT (IN) :: psq(stars%ng3) ! pseudo charge density (see psqpw.F)
32
      REAL,    INTENT (OUT):: vz(vacuum%nmzd,2)
33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48
!     ..
!     .. Local Scalar Parameters ..
      COMPLEX, PARAMETER :: ci = (0.0,1.0)
!     ..
!     .. Local Scalars ..
      COMPLEX sumq,vcons
      REAL bj0,bj1,dh,qzh,sigmaa(2)
      INTEGER ig3,imz,ivac,ncsh
!     ..
!     .. Local Arrays ..
      REAL f(vacuum%nmzd),sig(vacuum%nmzd),vtemp(vacuum%nmzd)
!     ..
!     .. Intrinsic Functions ..
      INTRINSIC cos,sin

     
49
      vz(:,1:vacuum%nvac) = 0.0  ! initialize potential
50 51 52

      ! obtain mesh point (ncsh) of charge sheet for external electric field
      !                                                 (X.Nie, IFF, 10/97)
53 54 55
      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
56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78

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

      vcons = -2.*fpi_const*ci
      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.EQ.2) THEN
79
        vz(1:vacuum%nmz,ivac) = vcons*sumq
80 81 82 83 84 85 86 87
      ENDIF

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

      ivac = 1     ! upper vacuum

      ! =================== DIRICHLET ==============================
88 89
      IF (field%efield%dirichlet) THEN
        vz(ncsh+1:vacuum%nmz,ivac) = field%efield%sig_b(1)
90

91
        CALL qsf(vacuum%delz,rht(1,ivac),sig,ncsh,1)
92 93 94 95
        sig(1:ncsh) = sig(ncsh)- sig(1:ncsh)
        CALL qsf(vacuum%delz,sig,vtemp,ncsh,1)

        DO imz = 1,ncsh
96
          vz(imz,ivac) = -fpi_const* (vtemp(ncsh)-vtemp(imz))&
97
     &                     + field%efield%sig_b(1)
98 99 100
        ENDDO

        sig1dh = sig(1)
101
        vz1dh = vz(1,ivac)   ! potential on vacuum boundary
102 103 104 105 106 107

        IF (vacuum%nvac == 1) RETURN


        ivac = 2     ! lower vacuum

108
        CALL qsf(vacuum%delz,rht(1,ivac),sig,ncsh,1)
109 110 111 112 113 114

        f(1:ncsh) = sig(1:ncsh) - rhobar*vacuum%dvac + sig1dh

        CALL qsf(vacuum%delz,f,vtemp,ncsh,1)

        DO imz = 1,ncsh
115 116
          vz(imz,ivac) = -fpi_const* (vtemp(imz)+sig1dh*vacuum%dvac&
     &                     -rhobar*vacuum%dvac*vacuum%dvac/2.) + vz(imz,ivac)&
117 118 119 120
     &                     +vz1dh
        END DO

        ! Force matching on the other side
121
        field%efield%vslope = (field%efield%sig_b(2)-vz(ncsh,1)) / (2*vacuum%delz*(ncsh+1)+vacuum%dvac)
122 123
        ivac = 1
        DO imz = 1,ncsh
124
          vz(imz,ivac) = vz(imz,ivac) + field%efield%vslope*vacuum%delz*(ncsh-imz+1)
125 126 127
        END DO
        ivac = 2
        DO imz = 1,ncsh
128
          vz(imz,ivac) = vz(imz,ivac)&
129
     &                    + field%efield%vslope*(vacuum%dvac+vacuum%delz*imz+vacuum%delz*ncsh)
130 131
        END DO

132
        vz(ncsh+1:vacuum%nmz,ivac) = field%efield%sig_b(2)
133 134 135 136

      ! =================== NEUMANN ==============================
      ELSE ! Neumann

137
      CALL qsf(vacuum%delz,rht(1,ivac),sig,vacuum%nmz,1)
138 139 140 141 142 143 144 145 146 147

      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)

      ! external electric field contribution (X.Nie, IFF, 10/97)
      !                                       corrected 10/99 mw
      DO imz = 1,ncsh
148
         vz(imz,ivac) = -fpi_const* (vtemp(vacuum%nmz)-vtemp(imz)) + vz(imz,ivac)&
149 150 151
     &                    -fpi_const*(imz-ncsh)*vacuum%delz*sigmaa(1)
      ENDDO
      DO imz =ncsh+1,vacuum%nmz
152
         vz(imz,ivac) = -fpi_const* (vtemp(vacuum%nmz)-vtemp(imz)) + vz(imz,ivac)
153 154
      ENDDO

155
      vz1dh = vz(1,ivac)   ! potential on vacuum boundary
156 157 158 159 160

      IF (vacuum%nvac.EQ.1) RETURN

      ivac = 2     ! lower vacuum

161
      CALL qsf(vacuum%delz,rht(1,ivac),sig,vacuum%nmz,1)
162 163 164 165 166 167 168

      f(1:vacuum%nmz) = sig(1:vacuum%nmz) - rhobar*vacuum%dvac + sig1dh

      CALL qsf(vacuum%delz,f,vtemp,vacuum%nmz,1)

      !   external electric field contribution
      DO imz = 1,ncsh
169 170
         vz(imz,ivac) = -fpi_const* (vtemp(imz)+sig1dh*vacuum%dvac-&
     &                    rhobar*vacuum%dvac*vacuum%dvac/2.) + vz1dh + vz(imz,ivac)
171 172
      ENDDO
      DO imz =ncsh+1,vacuum%nmz
173 174
         vz(imz,ivac) = -fpi_const* (vtemp(imz)+sig1dh*vacuum%dvac-&
     &                    rhobar*vacuum%dvac*vacuum%dvac/2.) + vz1dh + vz(imz,ivac)&
175 176 177 178 179 180 181
     &                    +fpi_const*(imz-ncsh)*vacuum%delz*sigmaa(2)
      ENDDO

      END IF ! Dirichlet (vs. Neumann)

      END SUBROUTINE vvac
      END MODULE m_vvac