Line data Source code
1 : module m_VYukawaFilm
2 :
3 : ! Computation of the film-case Yukawa potential for the preconditioning of the
4 : ! residual charge density in 5 steps:
5 : ! 1. pseudo-charge density generation
6 : ! 2. vacuum potential generation
7 : ! 3. interstitial potential generation
8 : ! 4. muffin-tin potential generation
9 : ! 5. modification for charge neutrality
10 :
11 : ! The Yukawa potential is the solution to the modified Helmholtz equation
12 : ! ( Delta - lambda^2 ) V_lambda = -4 pi ( rho_out - rho_in )
13 : ! subject to some conditions.
14 : ! The general scheme (steps 1 to 4) is the same as for the Poisson equation --
15 : ! we use Green function methods for the z-dependent vacuum and interstitial
16 : ! potentials as well as for the muffin-tin potential and apply Weinert's
17 : ! method.
18 : ! You can choose between two variants:
19 : ! 1. variant:
20 : ! zero Dirichlet boundary conditions at +/- infinity;
21 : ! multiplication with a decaying exponential in vacuum;
22 : ! modification in the film for charge neutrality (step 5)
23 : ! 2. variant:
24 : ! zero Dirichlet boundary conditions near the film boundary in vacuum (D/2+2R);
25 : ! modification in the film for charge neutrality (step 5)
26 : ! In both cases charge neutrality is broken.
27 : ! To restore charge neutrality, we need the integral over the potential to be
28 : ! zero.
29 : ! In step 5 we therefore solve the modified Helmholtz equation again with
30 : ! constant right-hand side, for an additive correction to the potential.
31 : ! The constant is chosen such that the integral over the final potential is
32 : ! zero.
33 :
34 : contains
35 :
36 :
37 :
38 0 : subroutine VYukawaFilm( stars, vacuum, cell, sym, input, mpi, atoms, sphhar, oneD, noco, den, &
39 : VYukawa )
40 :
41 : use m_constants
42 : use m_types
43 : use m_psqpw
44 : use m_vmts
45 : implicit none
46 :
47 : type(t_stars), intent(in) :: stars
48 : type(t_vacuum), intent(in) :: vacuum
49 : type(t_cell), intent(in) :: cell
50 : type(t_sym), intent(in) :: sym
51 : type(t_input), intent(in) :: input
52 : type(t_mpi), intent(in) :: mpi
53 : type(t_atoms), intent(in) :: atoms
54 : type(t_sphhar), intent(in) :: sphhar
55 : type(t_oneD), intent(in) :: oneD
56 : type(t_noco), intent(in) :: noco
57 : type(t_potden), intent(inout) :: den
58 :
59 : type(t_potden), intent(inout) :: VYukawa
60 :
61 0 : complex :: psq(stars%ng3)
62 0 : complex :: alphm(stars%ng2,2)
63 : real :: dh_prec
64 0 : real :: coshdh(stars%ng2)
65 :
66 :
67 : ! PSEUDO-CHARGE DENSITY
68 :
69 : call psqpw( mpi, atoms, sphhar, stars, vacuum, cell, input, sym, oneD, &
70 : den%pw(:,1), den%mt(:,:,:,1), den%vacz(:,:,1), .false., VYukawa%potdenType, &
71 0 : psq )
72 :
73 :
74 : ChooseVariant: if ( .true. ) then
75 :
76 : ! VACUUM POTENTIAL
77 :
78 : call VYukawaFilmVacuumVariant1( &
79 : stars, vacuum, cell, sym, input, atoms%rmt(1), &
80 : psq, den%vacxy(:,:,:,1), den%vacz(:,:,1), &
81 0 : VYukawa%vacxy, VYukawa%vacz, alphm )
82 :
83 :
84 : ! INTERSTITIAL POTENTIAL
85 :
86 : call VYukawaFilmInterstitialVariant1( &
87 : stars, vacuum, cell, sym, input, &
88 : psq, VYukawa%vacxy, VYukawa%vacz, alphm, &
89 0 : VYukawa%pw(:,1) )
90 :
91 : else ChooseVariant
92 :
93 : ! VACUUM POTENTIAL
94 :
95 : call VYukawaFilmVacuumVariant2( &
96 : stars, vacuum, cell, sym, input, 2*atoms%rmt(1), &
97 : psq, den%vacxy(:,:,:,1), den%vacz(:,:,1), &
98 : VYukawa%vacxy, VYukawa%vacz, alphm, dh_prec, coshdh )
99 :
100 :
101 : ! INTERSTITIAL POTENTIAL
102 :
103 : call VYukawaFilmInterstitialVariant2( &
104 : stars, vacuum, cell, sym, input, &
105 : psq, VYukawa%vacxy, VYukawa%vacz, alphm, dh_prec, coshdh, &
106 : VYukawa%pw(:,1) )
107 :
108 : end if ChooseVariant
109 :
110 :
111 : ! MUFFIN-TIN POTENTIAL
112 :
113 : call Vmts( input, mpi, stars, sphhar, atoms, sym, cell, oneD, &
114 : VYukawa%pw(:,1), den%mt(:,0:,:,1), VYukawa%potdenType, &
115 0 : VYukawa%mt(:,0:,:,1) )
116 :
117 :
118 : ! MODIFICATION FOR CHARGE NEUTRALITY
119 :
120 : call VYukawaModify( stars, vacuum, cell, sym, input, mpi, atoms, sphhar, oneD, noco, &
121 : den, &
122 0 : VYukawa )
123 :
124 :
125 0 : end subroutine VYukawaFilm
126 :
127 :
128 :
129 0 : subroutine VYukawaFilmVacuumVariant1( &
130 : stars, vacuum, cell, sym, input, rmt, &
131 0 : psq, rhtxy, rht, &
132 0 : VVxy, VVz, alphm )
133 :
134 : ! 1. part: Compute the contribution from the interstitial charge density to the vacuum potential as a function of q_xy and z (analytic expression for integral)
135 : ! 2. part: Compute the contribution from the vacuum charge density to the vacuum potential as a function of q_xy and z by numerical integration
136 :
137 : use m_ExpSave
138 : use m_constants
139 : use m_types
140 : use m_intgr, only: intgz1Reverse
141 : use m_qsf
142 : implicit none
143 :
144 : type(t_stars), intent(in) :: stars
145 : type(t_vacuum), intent(in) :: vacuum
146 : type(t_cell), intent(in) :: cell
147 : type(t_sym), intent(in) :: sym
148 : type(t_input), intent(in) :: input
149 : real, intent(in) :: rmt
150 : complex, intent(in) :: psq(stars%ng3)
151 : complex, intent(in) :: rhtxy(vacuum%nmzxyd,stars%ng2-1,2)
152 : real, intent(in) :: rht(vacuum%nmzd,2)
153 :
154 : complex, intent(out) :: VVxy(vacuum%nmzxyd,2:stars%ng2,2) ! this is the qxy /= 0 part of the vacuum potential
155 : real, intent(out) :: VVz(vacuum%nmzd,2) ! this is the qxy = 0 part of the vacuum potential
156 : complex, intent(out) :: alphm(stars%ng2,2) ! these are the integrals in upper and lower vacuum, now including the first star---integral for star ig2 is in alphm(ig2,ivac)
157 :
158 0 : complex :: sum_qz(2,stars%ng2)
159 0 : complex :: c_ph(-stars%mx3:stars%mx3,stars%ng2)
160 : complex :: signIqz
161 0 : real :: g_damped(stars%ng2), qz, sign, vcons(stars%ng2)
162 0 : real :: exp_m(vacuum%nmzd,stars%ng2), exp_p(vacuum%nmzd,stars%ng2)
163 0 : real :: expDhg(stars%ng2), expDg(stars%ng2)
164 0 : real :: z(vacuum%nmzd)
165 : integer :: iz, irec2, irec3, ivac, iqz
166 0 : complex :: fa(vacuum%nmzxyd,2:stars%ng2), fb(vacuum%nmzxyd,2:stars%ng2)
167 0 : complex :: alpha(vacuum%nmzxyd,2:stars%ng2,2), beta(vacuum%nmzxyd,2:stars%ng2,2), gamma(vacuum%nmzxyd,2:stars%ng2)
168 0 : real :: ga(vacuum%nmzd), gb(vacuum%nmzd)
169 0 : real :: delta(vacuum%nmzd,2), epsilon(vacuum%nmzd,2), zeta(vacuum%nmzd)
170 :
171 :
172 : ! DEFINITIONS / ALLOCATIONS / INITIALISATIONS
173 :
174 0 : do iz = 1, vacuum%nmz
175 0 : z(iz) = ( iz - 1 ) * vacuum%delz
176 : end do
177 0 : do irec2 = 1, stars%ng2
178 0 : g_damped(irec2) = sqrt( stars%sk2(irec2) ** 2 + input%preconditioning_param ** 2 )
179 0 : vcons(irec2) = tpi_const / g_damped(irec2)
180 0 : do iz = 1, vacuum%nmz
181 0 : exp_m(iz,irec2) = exp_save( - g_damped(irec2) * z(iz) )
182 0 : exp_p(iz,irec2) = exp_save( g_damped(irec2) * z(iz) )
183 : end do
184 0 : expDhg(irec2) = exp_save( - cell%z1 * g_damped(irec2) )
185 0 : expDg(irec2) = exp_save( -2 * cell%z1 * g_damped(irec2) )
186 : end do
187 0 : sum_qz = (0.,0.)
188 0 : VVxy = (0.,0.)
189 0 : VVz = 0.
190 :
191 :
192 : ! CONTRIBUTION FROM THE INTERSTITIAL CHARGE DENSITY
193 :
194 0 : do irec2 = 1, stars%ng2
195 0 : do ivac = 1, vacuum%nvac
196 0 : sign = 3. - 2. * ivac
197 0 : do iqz = -stars%mx3, stars%mx3
198 0 : irec3 = stars%ig(stars%kv2(1,irec2),stars%kv2(2,irec2),iqz)
199 : ! use only stars within the g_max sphere -> stars outside the sphere have per definition index ig3n = 0
200 0 : if ( irec3 /= 0 ) then
201 0 : c_ph(iqz,irec2) = stars%rgphs(stars%kv2(1,irec2),stars%kv2(2,irec2),iqz)
202 0 : qz = iqz * cell%bmat(3,3)
203 0 : signIqz = sign * ImagUnit * qz
204 0 : sum_qz(ivac,irec2) = sum_qz(ivac,irec2) + c_ph(iqz,irec2) * psq(irec3) * ( exp( signIqz * cell%z1 ) - exp( - signIqz * cell%z1 ) * expDg(irec2) ) / ( signIqz + g_damped(irec2) )
205 : endif
206 : enddo
207 0 : if( irec2 /= 1 ) then
208 0 : VVxy(1:vacuum%nmzxy,irec2,ivac) = vcons(irec2) * sum_qz(ivac,irec2) * exp_m(1:vacuum%nmzxy,irec2)
209 : else
210 0 : VVz(1:vacuum%nmz,ivac) = vcons(1) * sum_qz(ivac,1) * exp_m(1:vacuum%nmz,1)
211 : end if
212 : enddo
213 : enddo
214 :
215 :
216 : ! CONTRIBUTION FROM THE VACUUM CHARGE DENSITY
217 :
218 : ! shifting z:
219 0 : do irec2 = 1, stars%ng2
220 0 : exp_m(1:vacuum%nmz,irec2) = exp_m(1:vacuum%nmz,irec2) * expDhg(irec2)
221 0 : exp_p(1:vacuum%nmz,irec2) = exp_p(1:vacuum%nmz,irec2) / expDhg(irec2)
222 : end do
223 :
224 : ! case irec2 > 1:
225 0 : do irec2 = 2, stars%ng2
226 0 : do ivac = 1, vacuum%nvac
227 : ! integrands:
228 0 : fa(1:vacuum%nmzxy,irec2) = rhtxy(1:vacuum%nmzxy,irec2-1,ivac) * exp_m(1:vacuum%nmzxy,irec2)
229 0 : fb(1:vacuum%nmzxy,irec2) = rhtxy(1:vacuum%nmzxy,irec2-1,ivac) * exp_p(1:vacuum%nmzxy,irec2)
230 : ! integrals:
231 : ! alpha(z,q_xy,ivac) = int_z^infty rho(z',q_xy,ivac) exp(-sqrt(q_xy**2+prec_param**2)*z') dz'
232 : ! beta (z,q_xy,ivac) = int_{D/2}^z rho(z',q_xy,ivac) exp(+sqrt(q_xy**2+prec_param**2)*z') dz'
233 : ! where for z < 0 the lower vacuum charge density (ivac=2) is defined by rho(q_xy,z,ivac=2) := rho(q_xy,-z,ivac=2)
234 0 : call intgz1Reverse( fa(:,irec2), vacuum%delz, vacuum%nmzxy, alpha(:,irec2,ivac), .true. )
235 0 : call qsfComplex( vacuum%delz, fb(:,irec2), beta(:,irec2,ivac), vacuum%nmzxy, 1 )
236 : ! alphm(q_xy,ivac) = alpha(D/2,q_xy,ivac) --- these integrals are also needed for the interstitial potential
237 0 : alphm(irec2,ivac) = alpha(1,irec2,ivac)
238 : end do
239 0 : if ( vacuum%nvac == 1 ) then
240 0 : if ( sym%invs ) then
241 0 : alphm(irec2,2) = conjg( alphm(irec2,1) )
242 : else
243 0 : alphm(irec2,2) = alphm(irec2,1)
244 : end if
245 : end if
246 0 : do ivac = 1, vacuum%nvac
247 : gamma(1:vacuum%nmzxy,irec2) = exp_m(1:vacuum%nmzxy,irec2) * ( alphm(irec2,mod(ivac,2)+1) + beta(1:vacuum%nmzxy,irec2,ivac) ) &
248 0 : + exp_p(1:vacuum%nmzxy,irec2) * alpha(1:vacuum%nmzxy,irec2,ivac) ! mod(ivac,2)+1 outputs the other ivac value
249 0 : where ( 2. * gamma(:,irec2) == gamma(:,irec2) ) gamma(:,irec2) = cmplx( 0., 0. )
250 0 : VVxy(1:vacuum%nmzxy,irec2,ivac) = VVxy(1:vacuum%nmzxy,irec2,ivac) + vcons(irec2) * gamma(1:vacuum%nmzxy,irec2)
251 : end do
252 : end do
253 :
254 : ! case irec2 = 1:
255 0 : do ivac = 1, vacuum%nvac
256 0 : ga(1:vacuum%nmz) = rht(1:vacuum%nmz,ivac) * exp_m(1:vacuum%nmz,1)
257 0 : gb(1:vacuum%nmz) = rht(1:vacuum%nmz,ivac) * exp_p(1:vacuum%nmz,1)
258 0 : call intgz1Reverse( ga(:), vacuum%delz, vacuum%nmz, delta(:,ivac), .true. ) ! integrals
259 0 : call qsf( vacuum%delz, gb(:), epsilon(:,ivac), vacuum%nmz, 1 )
260 0 : alphm(1,ivac) = delta(1,ivac)
261 : end do
262 0 : if ( vacuum%nvac == 1 ) alphm(1,2) = alphm(1,1)
263 0 : do ivac = 1, vacuum%nvac
264 : zeta(1:vacuum%nmz) = exp_m(1:vacuum%nmz,1) * ( alphm(1,mod(ivac,2)+1) + epsilon(1:vacuum%nmz,ivac) ) &
265 0 : + exp_p(1:vacuum%nmz,1) * delta(1:vacuum%nmz,ivac)
266 0 : where ( 2. * zeta == zeta ) zeta = 0.
267 0 : VVz(1:vacuum%nmz,ivac) = VVz(1:vacuum%nmz,ivac) + vcons(1) * zeta(1:vacuum%nmz)
268 : end do
269 :
270 : ! damping in vacuum:
271 0 : do ivac = 1, vacuum%nvac
272 0 : VVz(1:vacuum%nmz,ivac) = VVz(1:vacuum%nmz,ivac) * exp( -0.1 / rmt * z(1:vacuum%nmz) )
273 0 : do irec2 = 2, stars%ng2
274 0 : VVxy(1:vacuum%nmzxy,irec2,ivac) = VVxy(1:vacuum%nmzxy,irec2,ivac) * exp( -0.1 / rmt * z(1:vacuum%nmzxy) )
275 : end do
276 : end do
277 :
278 :
279 0 : end subroutine VYukawaFilmVacuumVariant1
280 :
281 :
282 :
283 0 : subroutine VYukawaFilmInterstitialVariant1( &
284 : stars, vacuum, cell, sym, input, &
285 0 : psq, VVxy, VVz, alphm, &
286 0 : VIq )
287 :
288 : ! main parts:
289 : ! 1. part: Compute the contribution from the interstitial charge density to the interstitial potential as a function of q_xy and z (analytic expression for integral)
290 : ! 2. part: Add the contribution from the vacuum charge density to the interstitial potential, which had already been computed earlier for the vacuum potential
291 :
292 : ! 4. part: Compute the coefficients V^I(q_xy,q_z) from the function V^I(q_xy,z) by a 1D Fourier transform:
293 : ! V^I(q_xy,z) = sum_{q_z} V^I(q_xy,q_z) * exp( ImagUnit * q_z * z )
294 : ! In order to be able to match the interstitial and vacuum potentials smoothly at the interface, the Fourier transform is done
295 : ! in a slightly larger region. -> 3. part
296 : ! 3. part: Interpolate the vacuum potential in a small region surrounding the slab
297 :
298 : use m_ExpSave
299 : use m_constants
300 : use m_types
301 : use m_cfft
302 : implicit none
303 :
304 : type(t_stars), intent(in) :: stars
305 : type(t_vacuum), intent(in) :: vacuum
306 : type(t_cell), intent(in) :: cell
307 : type(t_sym), intent(in) :: sym
308 : type(t_input), intent(in) :: input
309 : complex, intent(in) :: psq(stars%ng3)
310 : complex, intent(in) :: VVxy(vacuum%nmzxyd,2:stars%ng2,2)
311 : real, intent(in) :: VVz(vacuum%nmzd,2)
312 : complex, intent(in) :: alphm(stars%ng2,2)
313 :
314 : complex, intent(out) :: VIq(stars%ng3)
315 :
316 : real :: partitioning, rz, qz, q
317 : integer :: irec2, irec3, iz, jz, ivac, iqz, nfft, nzmax, nzmin, nzdh, nLower, nUpper, jvac
318 0 : complex, allocatable :: VIz(:,:), eta(:,:)
319 0 : complex :: VIqz(-stars%mx3:stars%mx3,stars%ng2), c_ph(-stars%mx3:stars%mx3,stars%ng2)
320 0 : complex :: vcons1(stars%ng3)
321 0 : real :: VIzReal(3*stars%mx3,stars%ng2), VIzImag(3*stars%mx3,stars%ng2)
322 0 : real, allocatable :: exp_m(:,:), exp_p(:,:)
323 0 : real, allocatable :: z(:)
324 0 : real :: g_damped(stars%ng2), vcons2(stars%ng2), expDhg(stars%ng2)
325 :
326 :
327 : ! DEFINITIONS / ALLOCATIONS / INITIALISATIONS
328 :
329 : ! grid points z_i:
330 0 : nfft = 3 * stars%mx3 ! number of grid points for Fourier transform
331 0 : partitioning = 1. / real( nfft )
332 0 : nzmax = nfft / 2 ! index of maximal z below D~/2
333 0 : nzmin = nzmax - nfft + 1 ! index of minimal z above -D~/2
334 0 : nzdh = ceiling( cell%z1 / cell%amat(3,3) * nfft ) - 1 ! index of maximal z below D/2
335 0 : allocate( z(nzmin:nzmax) )
336 : ! definition of z_i: ! indexing: z_0 = 0; positive indices for z > 0; negative indices for z < 0
337 0 : do iz = nzmin, nzmax
338 0 : z(iz) = cell%amat(3,3) * iz * partitioning
339 : end do
340 : ! other variables:
341 0 : allocate( VIz(nzmin:nzmax,stars%ng2), eta(nzmin:nzmax,stars%ng2) )
342 0 : allocate( exp_m(nzmin:nzmax,stars%ng2), exp_p(nzmin:nzmax,stars%ng2) )
343 0 : do irec2 = 1, stars%ng2
344 0 : g_damped(irec2) = sqrt( stars%sk2(irec2) ** 2 + input%preconditioning_param ** 2 )
345 0 : vcons2(irec2) = -1. / ( 2. * g_damped(irec2) )
346 0 : do iz = nzmin, nzmax
347 0 : exp_m(iz,irec2) = exp_save( - g_damped(irec2) * z(iz) )
348 0 : exp_p(iz,irec2) = exp_save( g_damped(irec2) * z(iz) )
349 : end do
350 0 : expDhg(irec2) = exp_save( - cell%z1 * g_damped(irec2) )
351 : end do
352 0 : do irec3 = 1, stars%ng3
353 0 : vcons1(irec3) = fpi_const * psq(irec3) / ( stars%sk3(irec3) ** 2 + input%preconditioning_param ** 2 )
354 : end do
355 0 : VIz = (0.,0.)
356 0 : VIq = (0.,0.)
357 :
358 :
359 : ! CONTRIBUTION FROM THE INTERSTITIAL CHARGE DENSITY
360 :
361 : ! compute V^I(q_xy,z) as a function of q_xy and z
362 0 : do irec2 = 1, stars%ng2
363 0 : do iz = -nzdh, nzdh
364 0 : do iqz = -stars%mx3, stars%mx3
365 0 : irec3 = stars%ig(stars%kv2(1,irec2),stars%kv2(2,irec2),iqz)
366 0 : if ( irec3 /= 0 ) then ! use only stars within the g_max sphere
367 0 : c_ph(iqz,irec2) = stars%rgphs(stars%kv2(1,irec2),stars%kv2(2,irec2),iqz)
368 0 : qz = iqz * cell%bmat(3,3)
369 : VIz(iz,irec2) = VIz(iz,irec2) &
370 : + vcons1(irec3) * c_ph(iqz,irec2) * &
371 : ( exp( ImagUnit * qz * z(iz) ) &
372 : + vcons2(irec2) * expDhg(irec2) * &
373 : ( ( g_damped(irec2) + ImagUnit * qz ) * exp_p(iz,irec2) * exp( ImagUnit * qz * cell%z1 ) &
374 0 : + ( g_damped(irec2) - ImagUnit * qz ) * exp_m(iz,irec2) * exp( - ImagUnit * qz * cell%z1 ) ) )
375 : end if
376 : enddo
377 : end do
378 : ! irec2 loop continues
379 :
380 :
381 : ! CONTRIBUTION FROM THE VACUUM CHARGE DENSITY
382 :
383 : ! irec2 loop continues
384 0 : eta(-nzdh:nzdh,irec2) = exp_m(-nzdh:nzdh,irec2) * alphm(irec2,2) + exp_p(-nzdh:nzdh,irec2) * alphm(irec2,1)
385 0 : where ( 2.0 * eta(:,irec2) == eta(:,irec2) ) eta(:,irec2) = cmplx( 0.0, 0.0 )
386 0 : VIz(-nzdh:nzdh,irec2) = VIz(-nzdh:nzdh,irec2) + tpi_const / g_damped(irec2) * eta(-nzdh:nzdh,irec2)
387 : ! irec2 loop continues
388 :
389 :
390 : ! INTERPOLATION IN VACUUM REGION
391 :
392 : ! use Lagrange polynomials of order 3 to interpolate the vacuum potential outside I
393 : ! q, q-1 and q-2 (scaled with +/-1 or +/-0.5) are the factors of the Lagrange basis polynomials
394 : ! irec2 loop continues
395 0 : do ivac = 1, 2
396 0 : select case( ivac )
397 : case( 1 )
398 0 : nUpper = nzmax; nLower = nzdh + 1; jvac = 1
399 : case( 2 )
400 0 : nLower = nzmin; nUpper = -nzdh - 1; jvac = 2; if ( sym%invs .or. sym%zrfs ) jvac = 1
401 : end select
402 0 : do iz = nLower, nUpper
403 0 : rz = ( abs( z(iz) ) - cell%z1 ) / vacuum%delz + 1.0
404 0 : jz = rz ! index of maximal vacuum grid point below z_i
405 0 : q = rz - jz ! factor in Lagrange basis polynomials
406 0 : if ( irec2 == 1 ) then
407 : VIz(iz,irec2) = 0.5 * ( q - 1. ) * ( q - 2. ) * VVz(jz, jvac) &
408 : - q * ( q - 2. ) * VVz(jz+1,jvac) &
409 0 : + 0.5 * q * ( q - 1. ) * VVz(jz+2,jvac)
410 0 : else if ( jz + 2 <= vacuum%nmzxy ) then
411 : VIz(iz,irec2) = 0.5 * ( q - 1. ) * ( q - 2. ) * VVxy(jz, irec2,jvac) &
412 : - q * ( q - 2. ) * VVxy(jz+1,irec2,jvac) &
413 0 : + 0.5 * q * ( q - 1. ) * VVxy(jz+2,irec2,jvac)
414 0 : if ( ( sym%invs .and. .not. sym%zrfs ) .and. ivac == 2 ) VIz(iz,irec2) = conjg( VIz(iz,irec2) )
415 : end if
416 : end do
417 : end do
418 : end do ! irec2
419 :
420 :
421 : ! 1D FOURIER TRANSFORM TO FIND THE COEFFICIENTS V^I(q_xy,q_z)
422 :
423 : ! change the indexing for the subroutine cfft, and split real and imaginary parts
424 0 : VIzReal(1:nzmax+1,:) = real( VIz(0:nzmax,:) ); VIzReal(nzmax+2:nfft,:) = real( VIz(nzmin:-1,:) )
425 0 : VIzImag(1:nzmax+1,:) = aimag( VIz(0:nzmax,:) ); VIzImag(nzmax+2:nfft,:) = aimag( VIz(nzmin:-1,:) )
426 :
427 : ! V^I(q_xy,z) = sum_{q_z} V^I(q_xy,q_z) * exp( ImagUnit * q_z * z )
428 0 : do irec2 = 1, stars%ng2
429 0 : call cfft( VIzReal(:,irec2), VIzImag(:,irec2), nfft, nfft, nfft, -1 )
430 : ! irec2 loop continues
431 :
432 : ! reorder
433 0 : VIqz(0,irec2) = cmplx( VIzReal(1,irec2), VIzImag(1,irec2) )
434 0 : do iqz = 1, stars%mx3
435 0 : VIqz( iqz,irec2) = cmplx( VIzReal(iqz+1, irec2), VIzImag(iqz+1, irec2) )
436 0 : VIqz(-iqz,irec2) = cmplx( VIzReal(nfft+1-iqz,irec2), VIzImag(nfft+1-iqz,irec2) )
437 : end do
438 :
439 : ! add the computed components to V^I(q_xy,q_z):
440 0 : do iqz= -stars%mx3, stars%mx3
441 0 : irec3 = stars%ig(stars%kv2(1,irec2),stars%kv2(2,irec2),iqz)
442 0 : if ( irec3 /= 0 ) VIq(irec3) = VIq(irec3) + VIqz(iqz,irec2) * partitioning / ( stars%nstr(irec3) / stars%nstr2(irec2) )
443 : end do
444 : end do
445 :
446 :
447 0 : end subroutine VYukawaFilmInterstitialVariant1
448 :
449 :
450 :
451 0 : subroutine VYukawaFilmVacuumVariant2( &
452 : stars, vacuum, cell, sym, input, rmt, &
453 0 : psq, rhtxy, rht, &
454 0 : VVxy, VVz, alphm, dh_prec, coshdh )
455 :
456 : ! 1. part: Compute the contribution from the interstitial charge density to the vacuum potential as a function of q_xy and z (analytic expression for integral)
457 : ! 2. part: Compute the contribution from the vacuum charge density to the vacuum potential as a function of q_xy and z by numerical integration
458 :
459 : use m_ExpSave
460 : use m_constants
461 : use m_types
462 : use m_intgr, only: intgz1Reverse
463 : use m_qsf
464 : implicit none
465 :
466 : type(t_stars), intent(in) :: stars
467 : type(t_vacuum), intent(in) :: vacuum
468 : type(t_cell), intent(in) :: cell
469 : type(t_sym), intent(in) :: sym
470 : type(t_input), intent(in) :: input
471 : real, intent(in) :: rmt
472 : complex, intent(in) :: psq(stars%ng3)
473 : complex, intent(in) :: rhtxy(vacuum%nmzxyd,stars%ng2-1,2)
474 : real, intent(in) :: rht(vacuum%nmzd,2)
475 :
476 : complex, intent(out) :: VVxy(vacuum%nmzxyd,2:stars%ng2,2) ! this is the qxy /= 0 part of the vacuum potential
477 : real, intent(out) :: VVz(vacuum%nmzd,2) ! this is the qxy = 0 part of the vacuum potential
478 : complex, intent(out) :: alphm(stars%ng2,2) ! these are the integrals in upper and lower vacuum, now including the first star---integral for star ig2 is in alphm(ig2,ivac)
479 : real, intent(out) :: dh_prec
480 : real, intent(out) :: coshdh(stars%ng2)
481 :
482 : integer :: iz, irec2, irec3, ivac, iqz, nzdhprec, sign
483 : real :: dh, qz, qxy_numerics
484 0 : real, allocatable :: z(:)
485 0 : complex, dimension(stars%ng2) :: sum_qz
486 0 : complex, dimension(stars%ng3) :: vcons3
487 0 : real, dimension(stars%ng2) :: g_damped, vcons2
488 0 : complex, dimension(-stars%mx3:stars%mx3,stars%ng2) :: c_ph, quotq
489 0 : real, allocatable :: quotz(:,:), sinhz(:,:)
490 0 : real, dimension(-1:1,stars%ng2) :: quotvardh
491 0 : complex, dimension(-stars%mx3:stars%mx3) :: expdhqz
492 0 : complex, allocatable :: fa(:,:), fb(:,:)
493 0 : complex, allocatable :: alpha(:,:,:), beta(:,:,:), gamma(:,:)
494 0 : real, allocatable :: ga(:), gb(:)
495 0 : real, allocatable :: delta(:,:), epsilon(:,:), zeta(:)
496 :
497 :
498 : ! DEFINITIONS / ALLOCATIONS / INITIALISATIONS
499 :
500 0 : dh = cell%z1 ! half the thickness of the film
501 0 : nzdhprec = ceiling( rmt / vacuum%delz ) ! index of dh_prec, see below
502 0 : dh_prec = dh + ( nzdhprec - 1 ) * vacuum%delz ! dh_prec is about dh + R; preconditioning boundary
503 0 : qxy_numerics = sqrt( ( 100 / dh_prec ) ** 2 - input%preconditioning_param ** 2 )
504 0 : allocate( z(nzdhprec) )
505 0 : allocate( quotz(-nzdhprec:nzdhprec,stars%ng2) )
506 0 : allocate( sinhz(nzdhprec,stars%ng2) )
507 0 : do iz = 1, nzdhprec ! new boundary
508 0 : z(iz) = ( iz - 1 ) * vacuum%delz + dh
509 : end do
510 0 : do irec2 = 1, stars%ng2
511 0 : g_damped(irec2) = sqrt( stars%sk2(irec2) ** 2 + input%preconditioning_param ** 2 )
512 0 : vcons2(irec2) = fpi_const / g_damped(irec2)
513 0 : coshdh(irec2) = cosh( g_damped(irec2) * ( dh_prec - dh ) )
514 0 : if( stars%sk2(irec2) < qxy_numerics ) then ! numerics ok
515 0 : do iz = 1, nzdhprec
516 : quotz( iz,irec2) = ( cosh( g_damped(irec2) * z(iz) ) / cosh( g_damped(irec2) * dh_prec ) &
517 0 : + sinh( g_damped(irec2) * z(iz) ) / sinh( g_damped(irec2) * dh_prec ) ) / 2
518 : quotz(-iz,irec2) = ( cosh( g_damped(irec2) * z(iz) ) / cosh( g_damped(irec2) * dh_prec ) &
519 0 : - sinh( g_damped(irec2) * z(iz) ) / sinh( g_damped(irec2) * dh_prec ) ) / 2
520 : end do
521 : quotvardh( 1,irec2) = ( cosh( g_damped(irec2) * dh ) / sinh( g_damped(irec2) * dh_prec ) &
522 0 : + sinh( g_damped(irec2) * dh ) / cosh( g_damped(irec2) * dh_prec ) ) / 2
523 : quotvardh(-1,irec2) = ( cosh( g_damped(irec2) * dh ) / sinh( g_damped(irec2) * dh_prec ) &
524 0 : - sinh( g_damped(irec2) * dh ) / cosh( g_damped(irec2) * dh_prec ) ) / 2
525 : else ! numerical treatment necessary
526 0 : do iz = 1, nzdhprec
527 0 : quotz( iz,irec2) = exp_save( g_damped(irec2) * ( z(iz) - dh_prec ) )
528 0 : quotz(-iz,irec2) = exp_save( g_damped(irec2) * ( -z(iz) - dh_prec ) )
529 : end do
530 0 : quotvardh( 1,irec2) = quotz( 1,irec2)
531 0 : quotvardh(-1,irec2) = quotz(-1,irec2)
532 : end if
533 0 : do iz = 1, nzdhprec
534 0 : sinhz(iz,irec2) = sinh( g_damped(irec2) * ( dh_prec - z(iz) ) )
535 : end do
536 0 : do iqz = -stars%mx3, stars%mx3
537 0 : quotq(iqz,irec2) = ImagUnit * iqz * cell%bmat(3,3) / g_damped(irec2)
538 : end do
539 : end do
540 0 : sum_qz = (0.,0.)
541 0 : VVxy = (0.,0.)
542 0 : VVz = 0.
543 0 : do irec3 = 1, stars%ng3
544 0 : vcons3(irec3) = fpi_const * psq(irec3) / ( stars%sk3(irec3) ** 2 + input%preconditioning_param ** 2 )
545 : end do
546 0 : do iqz = -stars%mx3, stars%mx3
547 0 : expdhqz(iqz) = exp( ImagUnit * iqz * cell%bmat(3,3) * dh )
548 : end do
549 :
550 :
551 : ! CONTRIBUTION FROM THE INTERSTITIAL CHARGE DENSITY
552 :
553 0 : do ivac = 1, vacuum%nvac
554 0 : sign = 3 - 2 * ivac
555 0 : do irec2 = 1, stars%ng2
556 0 : do iqz = -stars%mx3, stars%mx3
557 0 : irec3 = stars%ig(stars%kv2(1,irec2),stars%kv2(2,irec2),iqz)
558 : ! use only stars within the g_max sphere -> stars outside the sphere have per definition index ig3n = 0
559 0 : if ( irec3 /= 0 ) then
560 0 : c_ph(iqz,irec2) = stars%rgphs(stars%kv2(1,irec2),stars%kv2(2,irec2),iqz)
561 : sum_qz(irec2) = sum_qz(irec2) + &
562 : c_ph(iqz,irec2) * vcons3(irec3) * &
563 : ( ( quotvardh( 1,irec2) - sign * quotq(iqz,irec2) * quotz( 1,irec2) ) * expdhqz( sign*iqz) &
564 0 : - ( quotvardh(-1,irec2) - sign * quotq(iqz,irec2) * quotz(-1,irec2) ) * expdhqz(-sign*iqz) )
565 : endif
566 : enddo
567 0 : if( irec2 /= 1 ) then
568 0 : VVxy(1:nzdhprec,irec2,ivac) = sum_qz(irec2) * sign * sinhz(1:nzdhprec,irec2)
569 : else
570 0 : VVz( 1:nzdhprec, ivac) = sum_qz(1) * sign * sinhz(1:nzdhprec,1)
571 : end if
572 : enddo
573 : enddo
574 :
575 :
576 : ! CONTRIBUTION FROM THE VACUUM CHARGE DENSITY
577 :
578 0 : allocate( fa(nzdhprec,2:stars%ng2), fb(nzdhprec,2:stars%ng2) )
579 0 : allocate( alpha(nzdhprec,2:stars%ng2,2), beta(nzdhprec,2:stars%ng2,2), gamma(nzdhprec,2:stars%ng2) )
580 : ! case irec2 > 1:
581 0 : do irec2 = 2, stars%ng2
582 0 : do ivac = 1, vacuum%nvac
583 : ! integrands:
584 0 : fa(1:nzdhprec,irec2) = rhtxy(1:nzdhprec,irec2-1,ivac) * sinhz(1:nzdhprec,irec2)
585 0 : fb(1:nzdhprec,irec2) = rhtxy(1:nzdhprec,irec2-1,ivac) * quotz(1:nzdhprec,irec2)
586 : ! integrals:
587 : ! alpha(z,q_xy,ivac) = int_z^infty rho(z',q_xy,ivac) sinhz dz'
588 : ! beta (z,q_xy,ivac) = int_{D/2}^z rho(z',q_xy,ivac) quotz dz'
589 : ! where for z < 0 the lower vacuum charge density (ivac=2) is defined by rho(q_xy,z,ivac=2) := rho(q_xy,-z,ivac=2)
590 0 : call intgz1Reverse( fa(:,irec2), vacuum%delz, nzdhprec, alpha(:,irec2,ivac), .false. )
591 0 : call qsfComplex( vacuum%delz, fb(:,irec2), beta(:,irec2,ivac), nzdhprec, 1 )
592 : ! alphm(q_xy,ivac) = alpha(D/2,q_xy,ivac) --- these integrals are also needed for the interstitial potential
593 0 : alphm(irec2,ivac) = alpha(1,irec2,ivac)
594 : end do
595 0 : if ( vacuum%nvac == 1 ) then
596 0 : if ( sym%invs ) then
597 0 : alphm(irec2,2) = conjg( alphm(irec2,1) )
598 : else
599 0 : alphm(irec2,2) = alphm(irec2,1)
600 : end if
601 : end if
602 0 : do ivac = 1, vacuum%nvac
603 : gamma(1:nzdhprec,irec2) = quotz(-1:-nzdhprec:-1,irec2) * alphm(irec2,mod(ivac,2)+1) &
604 : + quotz(1:nzdhprec,irec2) * alpha(1:nzdhprec,irec2,ivac) &
605 0 : + sinhz(1:nzdhprec,irec2) * beta(1:nzdhprec,irec2,ivac)
606 0 : VVxy(1:nzdhprec,irec2,ivac) = VVxy(1:nzdhprec,irec2,ivac) + vcons2(irec2) * gamma(1:nzdhprec,irec2)
607 : end do
608 : end do
609 :
610 :
611 0 : allocate( ga(nzdhprec), gb(nzdhprec) )
612 0 : allocate( delta(nzdhprec,2), epsilon(nzdhprec,2), zeta(nzdhprec) )
613 : ! case irec2 = 1:
614 0 : do ivac = 1, vacuum%nvac
615 0 : ga(1:nzdhprec) = rht(1:nzdhprec,ivac) * sinhz(1:nzdhprec,1)
616 0 : gb(1:nzdhprec) = rht(1:nzdhprec,ivac) * quotz(1:nzdhprec,1)
617 0 : call intgz1Reverse( ga(:), vacuum%delz, nzdhprec, delta(:,ivac), .false. )
618 0 : call qsf( vacuum%delz, gb(:), epsilon(:,ivac), nzdhprec, 1 )
619 0 : alphm(1,ivac) = delta(1,ivac)
620 : end do
621 0 : if ( vacuum%nvac == 1 ) alphm(1,2) = alphm(1,1)
622 0 : do ivac = 1, vacuum%nvac
623 : zeta(1:nzdhprec) = quotz(-1:-nzdhprec:-1,1) * alphm(1,mod(ivac,2)+1) &
624 : + quotz(1:nzdhprec,1) * delta(1:nzdhprec,ivac) &
625 0 : + sinhz(1:nzdhprec,1) * epsilon(1:nzdhprec,ivac)
626 0 : VVz(1:nzdhprec,ivac) = VVz(1:nzdhprec,ivac) + vcons2(1) * zeta(1:nzdhprec)
627 : end do
628 :
629 :
630 0 : end subroutine VYukawaFilmVacuumVariant2
631 :
632 :
633 :
634 0 : subroutine VYukawaFilmInterstitialVariant2( &
635 : stars, vacuum, cell, sym, input, &
636 0 : psq, VVxy, VVz, alphm, dh_prec, coshdh, &
637 0 : VIq )
638 :
639 : ! main parts:
640 : ! 1. part: Compute the contribution from the interstitial charge density to the interstitial potential as a function of q_xy and z (analytic expression for integral)
641 : ! 2. part: Add the contribution from the vacuum charge density to the interstitial potential, which had already been computed earlier for the vacuum potential
642 :
643 : ! 4. part: Compute the coefficients V^I(q_xy,q_z) from the function V^I(q_xy,z) by a 1D Fourier transform:
644 : ! V^I(q_xy,z) = sum_{q_z} V^I(q_xy,q_z) * exp( ImagUnit * q_z * z )
645 : ! In order to be able to match the interstitial and vacuum potentials smoothly at the interface, the Fourier transform is done
646 : ! in a slightly larger region. -> 3. part
647 : ! 3. part: Interpolate the vacuum potential in a small region surrounding the slab
648 :
649 : use m_ExpSave
650 : use m_constants
651 : use m_types
652 : use m_cfft
653 : implicit none
654 :
655 : type(t_stars), intent(in) :: stars
656 : type(t_vacuum), intent(in) :: vacuum
657 : type(t_cell), intent(in) :: cell
658 : type(t_sym), intent(in) :: sym
659 : type(t_input), intent(in) :: input
660 : complex, intent(in) :: psq(stars%ng3)
661 : complex, intent(in) :: VVxy(vacuum%nmzxyd,2:stars%ng2,2)
662 : real, intent(in) :: VVz(vacuum%nmzd,2)
663 : complex, intent(in) :: alphm(stars%ng2,2)
664 : real, intent(in) :: dh_prec
665 : real, intent(in) :: coshdh(stars%ng2)
666 :
667 : complex, intent(out) :: VIq(stars%ng3)
668 :
669 : real :: partitioning, rz, qz, q, qxy_numerics
670 : integer :: irec2, irec3, iz, jz, ivac, iqz, jvac
671 : integer :: nfft, nzmax, nzmin, nzdh, nLower, nUpper
672 0 : complex, allocatable :: VIz(:,:), eta(:,:)
673 0 : complex, allocatable :: expzqz(:,:)
674 0 : complex, dimension(-stars%mx3:stars%mx3,stars%ng2) :: VIqz, c_ph
675 0 : complex, dimension(-stars%mx3:stars%mx3,stars%ng2) :: qquot, qquottrigp, qquottrigm
676 0 : complex, dimension(stars%ng3) :: vcons3
677 0 : real, dimension(3*stars%mx3,stars%ng2) :: VIzReal, VIzImag
678 0 : real, allocatable :: quotz(:,:), sinhz(:,:)
679 0 : real, allocatable :: z(:)
680 0 : real, dimension(stars%ng2) :: g_damped, vcons2
681 :
682 :
683 : ! DEFINITIONS / ALLOCATIONS / INITIALISATIONS
684 :
685 : ! grid points z_i:
686 0 : qxy_numerics = sqrt( ( 100 / dh_prec ) ** 2 - input%preconditioning_param ** 2 )
687 0 : nfft = 3 * stars%mx3 ! number of grid points for Fourier transform
688 0 : partitioning = 1. / real( nfft )
689 0 : nzmax = nfft / 2 ! index of maximal z below D~/2
690 0 : nzmin = nzmax - nfft + 1 ! index of minimal z above -D~/2
691 0 : nzdh = ceiling( cell%z1 / cell%amat(3,3) * nfft ) - 1 ! index of maximal z below D/2
692 0 : allocate( z(nzmin:nzmax) )
693 : ! definition of z_i: ! indexing: z_0 = 0; positive indices for z > 0; negative indices for z < 0
694 0 : do iz = nzmin, nzmax
695 0 : z(iz) = cell%amat(3,3) * iz * partitioning
696 : end do
697 : ! other variables:
698 0 : allocate( VIz(nzmin:nzmax,stars%ng2) )
699 0 : allocate( eta(-nzdh:nzdh,stars%ng2) )
700 0 : allocate( sinhz(-nzdh:nzdh,stars%ng2) )
701 0 : allocate( quotz(-nzdh:nzdh,stars%ng2) )
702 0 : allocate( expzqz(-nzdh:nzdh,-stars%mx3:stars%mx3) )
703 0 : do irec2 = 1, stars%ng2
704 0 : g_damped(irec2) = sqrt( stars%sk2(irec2) ** 2 + input%preconditioning_param ** 2 )
705 0 : vcons2(irec2) = fpi_const / g_damped(irec2)
706 0 : if( stars%sk2(irec2) < qxy_numerics ) then ! numerics ok
707 0 : do iz = -nzdh, nzdh
708 : quotz( iz,irec2) = ( cosh( g_damped(irec2) * z(iz) ) / cosh( g_damped(irec2) * dh_prec ) &
709 0 : + sinh( g_damped(irec2) * z(iz) ) / sinh( g_damped(irec2) * dh_prec ) ) / 2
710 : end do
711 : else ! numerical treatment necessary
712 0 : do iz = -nzdh, nzdh
713 0 : quotz( iz,irec2) = exp_save( g_damped(irec2) * ( z(iz) - dh_prec ) )
714 : end do
715 : end if
716 0 : do iz = -nzdh, nzdh
717 0 : sinhz(iz,irec2) = sinh( g_damped(irec2) * ( dh_prec - z(iz) ) )
718 0 : do iqz = -stars%mx3, stars%mx3
719 0 : expzqz(iz,iqz) = exp( ImagUnit * iqz * cell%bmat(3,3) * z(iz) )
720 : end do
721 : end do
722 0 : do iqz = -stars%mx3, stars%mx3
723 0 : qquot(iqz,irec2) = ImagUnit * iqz * cell%bmat(3,3) / g_damped(irec2)
724 0 : qquottrigp(iqz,irec2) = coshdh(irec2) + sinhz(nzdh,irec2) * qquot(iqz,irec2)
725 0 : qquottrigm(iqz,irec2) = coshdh(irec2) - sinhz(nzdh,irec2) * qquot(iqz,irec2)
726 : end do
727 : end do
728 0 : do irec3 = 1, stars%ng3
729 0 : vcons3(irec3) = fpi_const * psq(irec3) / ( stars%sk3(irec3) ** 2 + input%preconditioning_param ** 2 )
730 : end do
731 0 : VIz = (0.,0.)
732 0 : VIq = (0.,0.)
733 :
734 :
735 : ! CONTRIBUTION FROM THE INTERSTITIAL CHARGE DENSITY
736 :
737 : ! compute V^I(q_xy,z) as a function of q_xy and z
738 0 : do irec2 = 1, stars%ng2
739 0 : do iz = -nzdh, nzdh
740 0 : do iqz = -stars%mx3, stars%mx3
741 0 : irec3 = stars%ig(stars%kv2(1,irec2),stars%kv2(2,irec2),iqz)
742 0 : if ( irec3 /= 0 ) then ! use only stars within the g_max sphere
743 0 : c_ph(iqz,irec2) = stars%rgphs(stars%kv2(1,irec2),stars%kv2(2,irec2),iqz)
744 0 : qz = iqz * cell%bmat(3,3)
745 : VIz(iz,irec2) = VIz(iz,irec2) &
746 : + vcons3(irec3) * c_ph(iqz,irec2) * &
747 : ( expzqz(iz,iqz) &
748 : - qquottrigp(iqz,irec2) * expzqz( nzdh,iqz) * quotz( iz,irec2) &
749 0 : - qquottrigm(iqz,irec2) * expzqz(-nzdh,iqz) * quotz(-iz,irec2) )
750 : end if
751 : enddo
752 : end do
753 : ! irec2 loop continues
754 :
755 :
756 : ! CONTRIBUTION FROM THE VACUUM CHARGE DENSITY
757 :
758 : ! irec2 loop continues
759 : eta(-nzdh:nzdh,irec2) = quotz(nzdh:-nzdh:-1,irec2) * alphm(irec2,2) &
760 0 : + quotz(-nzdh:nzdh, irec2) * alphm(irec2,1)
761 0 : VIz(-nzdh:nzdh,irec2) = VIz(-nzdh:nzdh,irec2) + vcons2(irec2) * eta(-nzdh:nzdh,irec2)
762 : ! irec2 loop continues
763 :
764 :
765 : ! INTERPOLATION IN VACUUM REGION
766 :
767 : ! use Lagrange polynomials of order 3 to interpolate the vacuum potential outside I
768 : ! q, q-1 and q-2 (scaled with +/-1 or +/-0.5) are the factors of the Lagrange basis polynomials
769 : ! irec2 loop continues
770 0 : do ivac = 1, 2
771 0 : select case( ivac )
772 : case( 1 )
773 0 : nUpper = nzmax; nLower = nzdh + 1; jvac = 1
774 : case( 2 )
775 0 : nLower = nzmin; nUpper = -nzdh - 1; jvac = 2; if ( sym%invs .or. sym%zrfs ) jvac = 1
776 : end select
777 0 : do iz = nLower, nUpper
778 0 : rz = ( abs( z(iz) ) - cell%z1 ) / vacuum%delz + 1.0
779 0 : jz = rz ! index of maximal vacuum grid point below z_i
780 0 : q = rz - jz ! factor in Lagrange basis polynomials
781 0 : if ( irec2 == 1 ) then
782 : VIz(iz,irec2) = 0.5 * ( q - 1. ) * ( q - 2. ) * VVz(jz, jvac) &
783 : - q * ( q - 2. ) * VVz(jz+1,jvac) &
784 0 : + 0.5 * q * ( q - 1. ) * VVz(jz+2,jvac)
785 0 : else if ( jz + 2 <= vacuum%nmzxy ) then
786 : VIz(iz,irec2) = 0.5 * ( q - 1. ) * ( q - 2. ) * VVxy(jz, irec2,jvac) &
787 : - q * ( q - 2. ) * VVxy(jz+1,irec2,jvac) &
788 0 : + 0.5 * q * ( q - 1. ) * VVxy(jz+2,irec2,jvac)
789 0 : if ( ( sym%invs .and. .not. sym%zrfs ) .and. ivac == 2 ) VIz(iz,irec2) = conjg( VIz(iz,irec2) )
790 : end if
791 : end do
792 : end do
793 : end do ! irec2
794 :
795 :
796 : ! 1D FOURIER TRANSFORM TO FIND THE COEFFICIENTS V^I(q_xy,q_z)
797 :
798 : ! change the indexing for the subroutine cfft, and split real and imaginary parts
799 0 : VIzReal(1:nzmax+1,:) = real( VIz(0:nzmax,:) ); VIzReal(nzmax+2:nfft,:) = real( VIz(nzmin:-1,:) )
800 0 : VIzImag(1:nzmax+1,:) = aimag( VIz(0:nzmax,:) ); VIzImag(nzmax+2:nfft,:) = aimag( VIz(nzmin:-1,:) )
801 :
802 : ! V^I(q_xy,z) = sum_{q_z} V^I(q_xy,q_z) * exp( ImagUnit * q_z * z )
803 0 : do irec2 = 1, stars%ng2
804 0 : call cfft( VIzReal(:,irec2), VIzImag(:,irec2), nfft, nfft, nfft, -1 )
805 : ! irec2 loop continues
806 :
807 : ! reorder
808 0 : VIqz(0,irec2) = cmplx( VIzReal(1,irec2), VIzImag(1,irec2) )
809 0 : do iqz = 1, stars%mx3
810 0 : VIqz( iqz,irec2) = cmplx( VIzReal(iqz+1, irec2), VIzImag(iqz+1, irec2) )
811 0 : VIqz(-iqz,irec2) = cmplx( VIzReal(nfft+1-iqz,irec2), VIzImag(nfft+1-iqz,irec2) )
812 : end do
813 :
814 : ! add the computed components to V^I(q_xy,q_z):
815 0 : do iqz= -stars%mx3, stars%mx3
816 0 : irec3 = stars%ig(stars%kv2(1,irec2),stars%kv2(2,irec2),iqz)
817 0 : if ( irec3 /= 0 ) VIq(irec3) = VIq(irec3) + VIqz(iqz,irec2) * partitioning / ( stars%nstr(irec3) / stars%nstr2(irec2) )
818 : end do
819 : end do
820 :
821 :
822 0 : end subroutine VYukawaFilmInterstitialVariant2
823 :
824 :
825 :
826 0 : subroutine VYukawaModify( stars, vacuum, cell, sym, input, mpi, atoms, sphhar, oneD, noco, den, &
827 : VYukawa )
828 :
829 : ! This subroutine adds a potential to the previously computed Yukawa
830 : ! potential to ensure charge neutrality.
831 : ! The added potential itself is a solution to the modified Helmholtz
832 : ! equation with constant right-hand side, where the constant is chosen such
833 : ! that charge neutrality is obtained.
834 : ! The charge is distributed only over the film region, and we therefore
835 : ! solve the differential equation subject to a boundary condition on the
836 : ! film surface, in contrast to the basic Yukawa potential above.
837 :
838 : use m_constants
839 : use m_types
840 : use m_vmts
841 : use m_constants
842 : use m_cdntot
843 : use m_cfft
844 : implicit none
845 :
846 : type(t_stars), intent(in) :: stars
847 : type(t_vacuum), intent(in) :: vacuum
848 : type(t_cell), intent(in) :: cell
849 : type(t_sym), intent(in) :: sym
850 : type(t_input), intent(in) :: input
851 : type(t_mpi), intent(in) :: mpi
852 : type(t_atoms), intent(in) :: atoms
853 : type(t_sphhar), intent(in) :: sphhar
854 : type(t_oneD), intent(in) :: oneD
855 : type(t_noco), intent(in) :: noco
856 : type(t_potden), intent(inout) :: den
857 : type(t_potden), intent(inout) :: VYukawa
858 :
859 : integer :: n, lh, irec3, iz, iqz, nfft, nzmax, nzmin, nzdh
860 : real :: q0, qhat, qbar, ldh, partitioning, dh
861 0 : real :: q(input%jspins), qis(input%jspins), qmt(atoms%ntype,input%jspins),qvac(2,input%jspins), qtot, qistot
862 0 : complex :: psq(stars%ng3)
863 0 : type(t_potden) :: VYukawaModification
864 0 : real, allocatable :: z(:)
865 0 : real :: VIzReal(3*stars%mx3), VIzImag(3*stars%mx3)
866 0 : complex, allocatable :: VIz(:)
867 0 : complex :: VIqz(-stars%mx3:stars%mx3)
868 :
869 :
870 : ! DEFINITIONS / ALLOCATIONS / INITIALISATIONS
871 :
872 : ! constants:
873 0 : dh = cell%z1 ! half the width of the film
874 : ! indexing of grid points z_i:
875 0 : nfft = 3 * stars%mx3 ! number of grid points for Fourier transform
876 0 : partitioning = 1. / real( nfft )
877 0 : nzmax = nfft / 2 ! index of maximal z below D~/2
878 0 : nzmin = nzmax - nfft + 1 ! index of minimal z above -D~/2
879 0 : nzdh = ceiling( dh / cell%amat(3,3) * nfft ) - 1 ! index of maximal z below D/2
880 0 : allocate( z(nzmin:nzmax) )
881 : ! definition of grid points z_i:
882 : ! indexing: z_0 = 0; positive indices for z > 0; negative indices for z < 0
883 0 : do iz = nzmin, nzmax
884 0 : z(iz) = cell%amat(3,3) * iz * partitioning
885 : end do
886 :
887 :
888 : ! INTEGRATION OF THE PREVIOUSLY COMPUTED YUKAWA POTENTIAL
889 :
890 : ! initialise VYukawaModification with in-going VYukawa and prepare for integration
891 0 : call VYukawaModification%init( stars, atoms, sphhar, vacuum, noco, input%jspins, 4 )
892 0 : call VYukawaModification%copyPotDen( VYukawa )
893 0 : do n = 1, atoms%ntype
894 0 : do lh = 0, sphhar%nlhd
895 0 : VYukawaModification%mt(1:atoms%jri(n),lh,n,1) = VYukawaModification%mt(1:atoms%jri(n),lh,n,1) * atoms%rmsh(1:atoms%jri(n),n) ** 2
896 : end do
897 : end do
898 :
899 : ! integrate the potential over the film region
900 0 : call integrate_cdn( stars, atoms, sym, vacuum, input, cell, oneD, VYukawaModification, q, qis, qmt, qvac, qtot, qistot )
901 0 : q0 = qtot / cell%area
902 0 : ldh = input%preconditioning_param * dh
903 0 : qhat = ( q0 / ( 2 * dh ) ) / ( sinh(ldh) / ( ldh * cosh( ldh ) ) - 1 )
904 0 : qbar = input%preconditioning_param ** 2 / fpi_const * qhat
905 :
906 :
907 : ! SET UP CONSTANT CHARGE DENSITY
908 :
909 : ! instead of den%pw(1,1) = qbar we directly set the pseudo charge density
910 0 : den%mt = 0; den%pw = 0; den%vacxy = 0; den%vacz = 0
911 0 : do n = 1, atoms%ntype
912 0 : den%mt(1:atoms%jri(n),0,n,1) = sfp_const * qbar * atoms%rmsh(1:atoms%jri(n),n) ** 2
913 : end do
914 0 : psq = cmplx(0.0,0.0); psq(1) = qbar
915 :
916 :
917 : ! CALCULATE THE INTERSTITIAL POTENTIAL AS A FUNCTION OF z
918 :
919 : ! initialise and calculate out-going modification potential; reuse VYukawaModification
920 0 : VYukawaModification%mt = 0; VYukawaModification%pw = 0; VYukawaModification%vacxy = 0; VYukawaModification%vacz = 0
921 :
922 0 : allocate( VIz(nzmin:nzmax) )
923 0 : VIz = (0.,0.)
924 0 : do iz = -nzdh+1, nzdh-1
925 0 : VIz(iz) = qhat * ( 1 - cosh( input%preconditioning_param * z(iz) ) / cosh( ldh ) )
926 : end do
927 :
928 :
929 : ! 1D FOURIER TRANSFORM TO FIND THE 3D-FOURIER COEFFICIENTS
930 :
931 : ! change the indexing for the subroutine cfft, and split real and imaginary parts
932 0 : VIzReal(1:nzmax+1) = real( VIz(0:nzmax) ); VIzReal(nzmax+2:nfft) = real( VIz(nzmin:-1) )
933 0 : VIzImag(1:nzmax+1) = aimag( VIz(0:nzmax) ); VIzImag(nzmax+2:nfft) = aimag( VIz(nzmin:-1) )
934 :
935 0 : call cfft( VIzReal(:), VIzImag(:), nfft, nfft, nfft, -1 )
936 :
937 : ! reorder
938 0 : VIqz = 0
939 0 : VIqz(0) = cmplx( VIzReal(1), VIzImag(1) )
940 0 : do iqz = 1, stars%mx3
941 0 : VIqz( iqz) = cmplx( VIzReal(iqz+1 ), VIzImag(iqz+1 ) )
942 0 : VIqz(-iqz) = cmplx( VIzReal(nfft+1-iqz), VIzImag(nfft+1-iqz) )
943 : end do
944 :
945 : ! add the computed components
946 0 : do iqz= -stars%mx3, stars%mx3
947 0 : irec3 = stars%ig(stars%kv2(1,1),stars%kv2(2,1),iqz)
948 0 : if ( irec3 /= 0 ) VYukawaModification%pw(irec3,1) = VYukawaModification%pw(irec3,1) + VIqz(iqz) * partitioning / ( stars%nstr(irec3) / stars%nstr2(1) )
949 : end do
950 :
951 :
952 : ! MUFFIN-TIN POTENTIAL
953 :
954 : call Vmts( input, mpi, stars, sphhar, atoms, sym, cell, oneD, &
955 : VYukawaModification%pw(:,1), den%mt(:,0:,:,1), VYukawaModification%potdenType, &
956 0 : VYukawaModification%mt(:,0:,:,1) )
957 :
958 :
959 : ! APPLYING THE MODIFICATION TO THE YUKAWA POTENTIAL
960 :
961 0 : call VYukawa%AddPotDen( VYukawa, VYukawaModification )
962 :
963 :
964 0 : end subroutine VYukawaModify
965 :
966 :
967 :
968 : end module m_VYukawaFilm
|
