Line data Source code
1 : !-------------------------------------------------------------------------------
2 : ! Copyright (c) 2022 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
3 : ! This file is part of FLEUR and available as free software under the conditions
4 : ! of the MIT license as expressed in the LICENSE file in more detail.
5 : !-------------------------------------------------------------------------------
6 :
7 : module m_vvac
8 : ! Legacy comment:
9 : ! ****************************************************************
10 : ! calculates the g(2-dim)=0 part of the vacuum coulomb potential *
11 : ! for general symmetry. c.l.fu, r.podloucky *
12 : ! ****************************************************************
13 : contains
14 37 : subroutine vvac(vacuum, stars, cell, input, field, psq, rht, vnew, rhobar, sig1dh, vz1dh, vslope, l_bind, vmz1dh, sigma_disc, sigma_disc2)
15 : !! Calculates the \(\boldsymbol{G}_{||}=0\) part of the vacuum Coulomb potential.
16 : !! There are two possible cases for Dirichlet and von Neumann boundary conditions, respectively.
17 : !! von Neumann case:
18 : !! upper vacuum \( z\ge D/2)\):
19 : !! $$V^{0}(z)=-4\pi[\int_{z}^{\infty}\sigma_{+}^{0}(z')dz'+(z-z_{\sigma})\sigma_{+}]=: V_{+}^{0}(z)$$
20 : !! with
21 : !! $$\sigma_{+}^{0}(z):=\int_{z}^{\infty}n_{V}^{0}(z')dz'$$
22 : use m_constants
23 : use m_qsf
24 : use m_types
25 :
26 : implicit none
27 :
28 : type(t_vacuum), intent(in) :: vacuum
29 : type(t_stars), intent(in) :: stars
30 : type(t_cell), intent(in) :: cell
31 : type(t_input), intent(in) :: input
32 : type(t_field), intent(in) :: field
33 :
34 : complex, intent(in) :: psq(stars%ng3)
35 : complex, intent(in) :: rht(vacuum%nmzd,2)
36 :
37 : complex, intent(out) :: vnew(vacuum%nmzd,2)
38 : complex, intent(out) :: rhobar
39 : complex, intent(out) :: sig1dh, vz1dh
40 : complex, intent(out) :: vslope
41 : logical, intent(in) :: l_bind
42 : complex, intent(out) :: vmz1dh
43 : complex, intent(in) :: sigma_disc(2)
44 :
45 : complex, optional, intent(in) :: sigma_disc2(2)
46 :
47 : complex :: sumq, newdp, newdm, newdp2, newdm2
48 : real :: bj0, bj1, qzh, sigmaa(2)
49 : integer :: ig3, imz, ivac, ncsh
50 37 : real :: f(vacuum%nmzd), sig(vacuum%nmzd), vtemp(vacuum%nmzd)
51 :
52 37 : vmz1dh = cmplx(0.0,0.0)
53 :
54 37 : newdp = cmplx(0.0,0.0)
55 37 : newdm = cmplx(0.0,0.0)
56 37 : newdp2 = cmplx(0.0,0.0)
57 37 : newdm2 = cmplx(0.0,0.0)
58 :
59 11834 : vnew(:,1:vacuum%nvac) = 0.0 ! initialize potential
60 :
61 : ! obtain mesh point (ncsh) of charge sheet for external electric field
62 37 : ncsh = field%efield%zsigma / vacuum%delz + 1.01
63 37 : sigmaa(1) = ( field%efield%sigma + field%efield%sig_b(1) ) / cell%area
64 37 : sigmaa(2) = ( field%efield%sigma + field%efield%sig_b(2) ) / cell%area
65 :
66 : ! g=0 vacuum potential due to neutral charge density
67 : ! inside slab and zero charge density outside
68 :
69 37 : rhobar = - psq(1)
70 37 : sumq = 0.0
71 :
72 37 : IF (l_bind) newdp = -psq(1) * cell%z1
73 : IF (l_bind) newdm = -psq(1) * cell%z1
74 37 : IF (l_bind) newdp2 = -psq(1) * cell%z1 * cell%z1
75 37 : IF (l_bind) newdm2 = psq(1) * cell%z1 * cell%z1
76 :
77 138057 : do ig3 = 2, stars%ng3
78 138057 : if (stars%ig2(ig3) == 1) then ! select g_|| = 0
79 730 : qzh = stars%kv3(3,ig3) * cell%bmat(3,3) * cell%z1
80 730 : bj0 = sin(qzh) / qzh
81 730 : rhobar = rhobar - psq(ig3) * bj0 * stars%nstr(ig3)
82 730 : if ( vacuum%nvac==2 ) then
83 340 : bj1 = ( sin(qzh) - qzh * cos(qzh) ) / ( qzh * qzh )
84 340 : sumq = sumq - 2. * fpi_const * ImagUnit * bj1 * psq(ig3) * cell%z1 *cell%z1
85 :
86 : ! New discontinuity correction
87 340 : IF (l_bind) newdp = newdp + ImagUnit * psq(ig3) * cmplx(cos(qzh), sin(qzh)) * cell%z1 / qzh
88 340 : IF (l_bind) newdm = newdm - ImagUnit * psq(ig3) * cmplx(cos(qzh),-sin(qzh)) * cell%z1 / qzh
89 340 : IF (l_bind) newdp2 = newdp2 + psq(ig3) * cmplx(cos(qzh), sin(qzh)) * cell%z1**2 / qzh**2
90 340 : IF (l_bind) newdm2 = newdm2 - psq(ig3) * cmplx(cos(qzh),-sin(qzh)) * cell%z1**2 / qzh**2
91 :
92 : end if
93 : end if
94 : end do ! --> qzh, bj0, bj1, psq finished; rhobar, sumq passed on and unchanged
95 :
96 : ! lower (nvac=2) vacuum
97 2537 : if ( vacuum%nvac == 2 ) vnew(1:vacuum%nmz,2) = sumq
98 :
99 : ! g=0 vacuum potential due to
100 : ! negative of rhobar + vacuum (g=0) charge ----> v2(z)
101 :
102 37 : ivac = 1 ! upper vacuum
103 :
104 37 : if ( field%efield%dirichlet ) then ! Dirichlet
105 0 : vnew(ncsh+1:vacuum%nmz,ivac) = field%efield%sig_b(1)
106 0 : call qsf( vacuum%delz, REAL(rht(:,ivac)), sig, ncsh, 1 )
107 0 : sig(1:ncsh) = sig(ncsh) - sig(1:ncsh)
108 0 : call qsf( vacuum%delz, sig, vtemp, ncsh, 1 )
109 0 : do imz = 1, ncsh
110 0 : vnew(imz,ivac) = - fpi_const * ( vtemp(ncsh) - vtemp(imz) ) + field%efield%sig_b(1)
111 : end do
112 0 : sig1dh = sig(1)
113 0 : vz1dh = vnew(1,ivac) ! potential on vacuum boundary
114 27 : if ( vacuum%nvac == 1 ) return
115 :
116 0 : ivac = 2 ! lower vacuum
117 0 : call qsf( vacuum%delz, REAL(rht(:,ivac)), sig, ncsh, 1 )
118 0 : f(1:ncsh) = sig(1:ncsh) - rhobar*vacuum%dvac + sig1dh
119 0 : call qsf( vacuum%delz, f, vtemp, ncsh, 1 )
120 0 : do imz = 1,ncsh
121 0 : vnew(imz,ivac) = - fpi_const * ( vtemp(imz) + sig1dh * vacuum%dvac - rhobar * vacuum%dvac * vacuum%dvac / 2. ) + vnew(imz,ivac) + vz1dh
122 : end do
123 :
124 : ! force matching on the other side
125 0 : vslope = ( field%efield%sig_b(2) - vnew(ncsh,1) ) / ( 2 * vacuum%delz * ( ncsh + 1 ) + vacuum%dvac )
126 0 : ivac = 1
127 0 : do imz = 1, ncsh
128 0 : vnew(imz,ivac) = vnew(imz,ivac) + vslope * vacuum%delz * ( ncsh - imz + 1 )
129 : end do
130 0 : ivac = 2
131 0 : do imz = 1, ncsh
132 0 : vnew(imz,ivac) = vnew(imz,ivac) + vslope * ( vacuum%dvac + vacuum%delz * imz + vacuum%delz * ncsh )
133 : end do
134 0 : vnew(ncsh+1:vacuum%nmz,ivac) = field%efield%sig_b(2)
135 : else ! Neumann
136 9324 : call qsf( vacuum%delz, REAL(rht(:,ivac)), sig, vacuum%nmz, 1 )
137 37 : sig1dh = sig(vacuum%nmz) - sigmaa(1) ! need to include contribution from electric field
138 9287 : sig(1:vacuum%nmz) = sig(vacuum%nmz) - sig(1:vacuum%nmz)
139 37 : call qsf( vacuum%delz, sig, vtemp, vacuum%nmz, 1 )
140 : ! external electric field contribution
141 3774 : do imz = 1, ncsh
142 3774 : vnew(imz,ivac) = - fpi_const * ( vtemp(vacuum%nmz) - vtemp(imz) ) + vnew(imz,ivac) - fpi_const * ( imz - ncsh ) * vacuum%delz * sigmaa(1)
143 : end do
144 5550 : do imz = ncsh + 1, vacuum%nmz
145 5550 : vnew(imz,ivac) = - fpi_const * ( vtemp(vacuum%nmz) - vtemp(imz) ) + vnew(imz,ivac)
146 : end do
147 37 : vz1dh = vnew(1,ivac) ! potential on vacuum boundary
148 37 : if ( vacuum%nvac == 1 ) return
149 :
150 10 : ivac = 2 ! lower vacuum
151 2510 : call qsf( vacuum%delz, REAL(rht(:,ivac)), sig, vacuum%nmz, 1 )
152 2510 : f(1:vacuum%nmz) = sig(1:vacuum%nmz) - rhobar * vacuum%dvac + sig1dh
153 10 : call qsf( vacuum%delz, f, vtemp, vacuum%nmz, 1 )
154 : ! external electric field contribution
155 1020 : do imz = 1, ncsh
156 : vnew(imz,ivac) = - fpi_const * ( vtemp(imz) + sig1dh * vacuum%dvac - rhobar * vacuum%dvac * vacuum%dvac / 2. ) + vz1dh + vnew(imz,ivac) &
157 : - fpi_const * (sigma_disc(1) * ( vacuum%dvac / 2. - (imz-1) * vacuum%delz ) - sigma_disc(2) * ( vacuum%dvac / 2. + (imz-1) * vacuum%delz )) &! Discontinuity correction
158 1020 : - fpi_const * (-newdp2-newdm2+newdp * ( vacuum%dvac / 2. - (imz-1) * vacuum%delz ) - newdm * ( vacuum%dvac / 2. + (imz-1) * vacuum%delz ))! New discontinuity correction
159 : !if (present(sigma_disc2)) vnew(imz,ivac) = vnew(imz,ivac) - fpi_const * (sigma_disc2(1)+sigma_disc2(2))
160 : end do
161 1500 : do imz = ncsh + 1, vacuum%nmz
162 : vnew(imz,ivac) = - fpi_const * ( vtemp(imz) + sig1dh * vacuum%dvac - rhobar * vacuum%dvac * vacuum%dvac / 2. ) + vz1dh + vnew(imz,ivac) &
163 : + fpi_const * ( imz - ncsh ) * vacuum%delz * sigmaa(2) & ! Discontinuity correction
164 : - fpi_const * (sigma_disc(1) * ( vacuum%dvac / 2. - (imz-1) * vacuum%delz ) - sigma_disc(2) * ( vacuum%dvac / 2. + (imz-1) * vacuum%delz )) &! Discontinuity correction
165 1500 : - fpi_const * (-newdp2-newdm2+newdp * ( vacuum%dvac / 2. - (imz-1) * vacuum%delz ) - newdm * ( vacuum%dvac / 2. + (imz-1) * vacuum%delz )) ! New discontinuity correction
166 : !if (present(sigma_disc2)) vnew(imz,ivac) = vnew(imz,ivac) - fpi_const * (sigma_disc2(1)+sigma_disc2(2))
167 : end do
168 : !if (l_bind) then
169 : ! ! Fix the potential to 0 at -infinity and save the resulting value at the vacuum border -D/2
170 : ! vnew(:,2) = vnew(:,2) - vnew(vacuum%nmz,2)
171 : ! vmz1dh = vnew(1,2)
172 : !end if
173 : ! Discontinuity correction
174 : !if (l_bind) then
175 : ! Fix the potential to 0 at -infinity and save the resulting value at the vacuum border -D/2
176 : !vnew(:,2) = cmplx(0.0,0.0)
177 : !vmz1dh = vnew(1,2)
178 : !end if
179 : end if ! Dirichlet/Neumann
180 : end subroutine vvac
181 : end module m_vvac
|
