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 266 : subroutine vvac(vacuum, stars, cell, input, field, psq, rht, vnew, rhobar, sig1dh, vz1dh, vslope, vmz1dh,l_dfptvgen)
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 : complex, intent(out) :: vmz1dh
42 : logical, intent(in) :: l_dfptvgen
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 266 : real :: f(vacuum%nmzd), sig(vacuum%nmzd), vtemp(vacuum%nmzd)
51 :
52 : !vmz1dh = cmplx(0.0,0.0)
53 :
54 : ! newdp = cmplx(0.0,0.0)
55 : ! newdm = cmplx(0.0,0.0)
56 : ! newdp2 = cmplx(0.0,0.0)
57 : ! newdm2 = cmplx(0.0,0.0)
58 :
59 124260 : vnew(:,1:vacuum%nvac) = 0.0 ! initialize potential
60 :
61 : ! obtain mesh point (ncsh) of charge sheet for external electric field
62 266 : ncsh = field%efield%zsigma / vacuum%delz + 1.01
63 266 : sigmaa = 0.0
64 266 : if (.not. l_dfptvgen) then
65 56 : sigmaa(1) = ( field%efield%sigma + field%efield%sig_b(1) ) / cell%area
66 56 : sigmaa(2) = ( field%efield%sigma + field%efield%sig_b(2) ) / cell%area
67 : end if
68 : ! g=0 vacuum potential due to neutral charge density
69 : ! inside slab and zero charge density outside
70 :
71 266 : rhobar = - psq(1)
72 266 : sumq = 0.0
73 :
74 : !IF (l_bind) newdp = -psq(1) * cell%z1
75 : !IF (l_bind) newdm = -psq(1) * cell%z1
76 : !IF (l_bind) newdp2 = -psq(1) * cell%z1 * cell%z1
77 : !IF (l_bind) newdm2 = psq(1) * cell%z1 * cell%z1
78 :
79 1368307 : do ig3 = 2, stars%ng3
80 1368307 : if (stars%ig2(ig3) == 1) then ! select g_|| = 0
81 10755 : qzh = stars%kv3(3,ig3) * cell%bmat(3,3) * cell%z1
82 10755 : bj0 = sin(qzh) / qzh
83 10755 : rhobar = rhobar - psq(ig3) * bj0 * stars%nstr(ig3)
84 10755 : if ( vacuum%nvac==2 ) then
85 10196 : bj1 = ( sin(qzh) - qzh * cos(qzh) ) / ( qzh * qzh )
86 10196 : sumq = sumq - 2. * fpi_const * ImagUnit * bj1 * psq(ig3) * cell%z1 *cell%z1
87 :
88 : ! ! New discontinuity correction
89 : ! IF (l_bind) newdp = newdp + ImagUnit * psq(ig3) * cmplx(cos(qzh), sin(qzh)) * cell%z1 / qzh
90 : ! IF (l_bind) newdm = newdm - ImagUnit * psq(ig3) * cmplx(cos(qzh),-sin(qzh)) * cell%z1 / qzh
91 : ! IF (l_bind) newdp2 = newdp2 + psq(ig3) * cmplx(cos(qzh), sin(qzh)) * cell%z1**2 / qzh**2
92 : ! IF (l_bind) newdm2 = newdm2 - psq(ig3) * cmplx(cos(qzh),-sin(qzh)) * cell%z1**2 / qzh**2
93 :
94 : end if
95 : end if
96 : end do ! --> qzh, bj0, bj1, psq finished; rhobar, sumq passed on and unchanged
97 :
98 : ! lower (nvac=2) vacuum
99 57266 : if ( vacuum%nvac == 2 ) vnew(1:vacuum%nmz,2) = sumq
100 :
101 : ! g=0 vacuum potential due to
102 : ! negative of rhobar + vacuum (g=0) charge ----> v2(z)
103 :
104 266 : ivac = 1 ! upper vacuum
105 :
106 266 : if ( field%efield%dirichlet .and. .not. l_dfptvgen ) then ! Dirichlet
107 0 : vnew(ncsh+1:vacuum%nmz,ivac) = field%efield%sig_b(1)
108 0 : call qsf( vacuum%delz, REAL(rht(:,ivac)), sig, ncsh, 1 )
109 0 : sig(1:ncsh) = sig(ncsh) - sig(1:ncsh)
110 0 : call qsf( vacuum%delz, sig, vtemp, ncsh, 1 )
111 0 : do imz = 1, ncsh
112 0 : vnew(imz,ivac) = - fpi_const * ( vtemp(ncsh) - vtemp(imz) ) + field%efield%sig_b(1)
113 : end do
114 0 : sig1dh = sig(1)
115 0 : vz1dh = vnew(1,ivac) ! potential on vacuum boundary
116 38 : if ( vacuum%nvac == 1 ) return
117 :
118 0 : ivac = 2 ! lower vacuum
119 0 : call qsf( vacuum%delz, REAL(rht(:,ivac)), sig, ncsh, 1 )
120 0 : f(1:ncsh) = sig(1:ncsh) - rhobar*vacuum%dvac + sig1dh
121 0 : call qsf( vacuum%delz, f, vtemp, ncsh, 1 )
122 0 : do imz = 1,ncsh
123 0 : vnew(imz,ivac) = - fpi_const * ( vtemp(imz) + sig1dh * vacuum%dvac - rhobar * vacuum%dvac * vacuum%dvac / 2. ) + vnew(imz,ivac) + vz1dh
124 : end do
125 :
126 : ! force matching on the other side
127 0 : vslope = ( field%efield%sig_b(2) - vnew(ncsh,1) ) / ( 2 * vacuum%delz * ( ncsh + 1 ) + vacuum%dvac )
128 0 : ivac = 1
129 0 : do imz = 1, ncsh
130 0 : vnew(imz,ivac) = vnew(imz,ivac) + vslope * vacuum%delz * ( ncsh - imz + 1 )
131 : end do
132 0 : ivac = 2
133 0 : do imz = 1, ncsh
134 0 : vnew(imz,ivac) = vnew(imz,ivac) + vslope * ( vacuum%dvac + vacuum%delz * imz + vacuum%delz * ncsh )
135 : end do
136 0 : vnew(ncsh+1:vacuum%nmz,ivac) = field%efield%sig_b(2)
137 : else ! Neumann
138 67032 : call qsf( vacuum%delz, REAL(rht(:,ivac)), sig, vacuum%nmz, 1 )
139 266 : sig1dh = sig(vacuum%nmz) - sigmaa(1) ! need to include contribution from electric field
140 66766 : sig(1:vacuum%nmz) = sig(vacuum%nmz) - sig(1:vacuum%nmz)
141 266 : call qsf( vacuum%delz, sig, vtemp, vacuum%nmz, 1 )
142 : ! external electric field contribution
143 27132 : do imz = 1, ncsh
144 27132 : vnew(imz,ivac) = - fpi_const * ( vtemp(vacuum%nmz) - vtemp(imz) ) + vnew(imz,ivac) - fpi_const * ( imz - ncsh ) * vacuum%delz * sigmaa(1)
145 : end do
146 39900 : do imz = ncsh + 1, vacuum%nmz
147 39900 : vnew(imz,ivac) = - fpi_const * ( vtemp(vacuum%nmz) - vtemp(imz) ) + vnew(imz,ivac)
148 : end do
149 266 : vz1dh = vnew(1,ivac) ! potential on vacuum boundary
150 266 : if ( vacuum%nvac == 1 ) return
151 :
152 228 : ivac = 2 ! lower vacuum
153 57228 : call qsf( vacuum%delz, REAL(rht(:,ivac)), sig, vacuum%nmz, 1 )
154 57228 : f(1:vacuum%nmz) = sig(1:vacuum%nmz) - rhobar * vacuum%dvac + sig1dh
155 228 : call qsf( vacuum%delz, f, vtemp, vacuum%nmz, 1 )
156 : ! external electric field contribution
157 23256 : do imz = 1, ncsh
158 23256 : vnew(imz,ivac) = - fpi_const * ( vtemp(imz) + sig1dh * vacuum%dvac - rhobar * vacuum%dvac * vacuum%dvac / 2. ) + vz1dh + vnew(imz,ivac) !&
159 : ! - fpi_const * (sigma_disc(1) * ( vacuum%dvac / 2. - (imz-1) * vacuum%delz ) - sigma_disc(2) * ( vacuum%dvac / 2. + (imz-1) * vacuum%delz )) &! Discontinuity correction
160 : ! - fpi_const * (-newdp2-newdm2+newdp * ( vacuum%dvac / 2. - (imz-1) * vacuum%delz ) - newdm * ( vacuum%dvac / 2. + (imz-1) * vacuum%delz ))! New discontinuity correction
161 : !if (present(sigma_disc2)) vnew(imz,ivac) = vnew(imz,ivac) - fpi_const * (sigma_disc2(1)+sigma_disc2(2))
162 : end do
163 34200 : do imz = ncsh + 1, vacuum%nmz
164 : vnew(imz,ivac) = - fpi_const * ( vtemp(imz) + sig1dh * vacuum%dvac - rhobar * vacuum%dvac * vacuum%dvac / 2. ) + vz1dh + vnew(imz,ivac) &
165 34200 : + fpi_const * ( imz - ncsh ) * vacuum%delz * sigmaa(2) !& ! Discontinuity correction
166 : ! - fpi_const * (sigma_disc(1) * ( vacuum%dvac / 2. - (imz-1) * vacuum%delz ) - sigma_disc(2) * ( vacuum%dvac / 2. + (imz-1) * vacuum%delz )) &! Discontinuity correction
167 : ! - fpi_const * (-newdp2-newdm2+newdp * ( vacuum%dvac / 2. - (imz-1) * vacuum%delz ) - newdm * ( vacuum%dvac / 2. + (imz-1) * vacuum%delz )) ! New discontinuity correction
168 : !if (present(sigma_disc2)) vnew(imz,ivac) = vnew(imz,ivac) - fpi_const * (sigma_disc2(1)+sigma_disc2(2))
169 : end do
170 : !if (l_bind) then
171 : ! ! Fix the potential to 0 at -infinity and save the resulting value at the vacuum border -D/2
172 : ! vnew(:,2) = vnew(:,2) - vnew(vacuum%nmz,2)
173 : ! vmz1dh = vnew(1,2)
174 : !end if
175 : ! Discontinuity correction
176 : !if (l_bind) then
177 : ! Fix the potential to 0 at -infinity and save the resulting value at the vacuum border -D/2
178 : !vnew(:,2) = cmplx(0.0,0.0)
179 : !vmz1dh = vnew(1,2)
180 : !end if
181 : end if ! Dirichlet/Neumann
182 : end subroutine vvac
183 : end module m_vvac
|
