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