Octopus
xc_photons.F90
Go to the documentation of this file.
1!! Copyright (C) 2022 I.-T Lu
2!!
3!! This program is free software; you can redistribute it and/or modify
4!! it under the terms of the GNU General Public License as published by
5!! the Free Software Foundation; either version 2, or (at your option)
6!! any later version.
7!!
8!! This program is distributed in the hope that it will be useful,
9!! but WITHOUT ANY WARRANTY; without even the implied warranty of
10!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
11!! GNU General Public License for more details.
12!!
13!! You should have received a copy of the GNU General Public License
14!! along with this program; if not, write to the Free Software
15!! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
16!! 02110-1301, USA.
17!!
18
19#include "global.h"
20
27
29 use debug_oct_m
31 use epot_oct_m
32 use global_oct_m
33 use grid_oct_m
34 use iso_c_binding
37 use mesh_oct_m
44 use space_oct_m
46 use parser_oct_m
49
50 implicit none
51
52 private
53 public :: xc_photons_t
54
63 !
64 type xc_photons_t
65 private
66 integer :: method = 0
67 real(real64), allocatable, public :: vpx(:)
68 real(real64), public :: ex
69 type(photon_mode_t) :: pt
70 real(real64) :: pxlda_kappa
71 real(real64) :: eta_c
73 integer :: energy_method = 0
74 logical :: lcorrelations = .false.
75 logical :: llamb_re_mass = .false.
76 logical :: llamb_freespace =.false.
77 real(real64) :: lamb_omega
78
79 logical, public :: lpfmf = .false.
80 real(real64), allocatable, public :: mf_vector_potential(:)
81 real(real64), allocatable :: jp_proj_eo(:,:)
84
85 contains
86 procedure :: init => xc_photons_init
87 procedure :: end => xc_photons_end
88 procedure :: wants_to_renormalize_mass => xc_photons_wants_to_renormalize_mass
89 procedure :: get_renormalized_mass => xc_photons_get_renormalized_emass
90 procedure :: mf_dump => xc_photons_mf_dump
91 procedure :: mf_load => xc_photons_mf_load
92 procedure :: v_ks => xc_photons_v_ks
93 procedure :: add_mean_field => xc_photons_add_mean_field
94 end type xc_photons_t
95
96 ! the PhotonXCXCMethod
97 integer, private, parameter :: &
98 XC_PHOTONS_NONE = 0, &
99 xc_photons_lda = 1, &
101
102contains
103
104 ! ---------------------------------------------------------
106 !
107 subroutine xc_photons_init(xc_photons, namespace, xc_photon, space, gr, st)
108 class(xc_photons_t), intent(out) :: xc_photons
109 type(namespace_t), intent(in) :: namespace
110 integer, intent(in) :: xc_photon
111 class(space_t), intent(in) :: space
112 type(grid_t), intent(in) :: gr
113 type(states_elec_t), intent(in) :: st
114
115 push_sub(xc_photons_init)
116
117 xc_photons%lpfmf = .false.
118
119 call messages_experimental("XCPhotonFunctional /= none")
120
121 call photon_mode_init(xc_photons%pt, namespace, space%dim, .true.)
122 call photon_mode_set_n_electrons(xc_photons%pt, st%qtot)
124 select case(xc_photon)
125 case(option__xcphotonfunctional__photon_x_lda)
126 xc_photons%method = xc_photons_lda
127 xc_photons%lcorrelations = .false.
128 case(option__xcphotonfunctional__photon_xc_lda)
129 xc_photons%method = xc_photons_lda
130 xc_photons%lcorrelations = .true.
131 case(option__xcphotonfunctional__photon_x_wfn)
132 xc_photons%method = xc_photons_wfs
133 xc_photons%lcorrelations = .false.
134 case(option__xcphotonfunctional__photon_xc_wfn)
135 xc_photons%method = xc_photons_wfs
136 xc_photons%lcorrelations = .true.
137 case default
138 xc_photons%method = xc_photons_none
139 return
140 end select
141
142
143 if (xc_photons%method == xc_photons_lda) then
144
145 !%Variable PhotonXCLDAKappa
146 !%Type float
147 !%Default 1.0
148 !%Section Hamiltonian::XC
149 !%Description
150 !% the scaling factor for px-LDA potential
151 !%End
152 call parse_variable(namespace, 'PhotonXCLDAKappa', m_one, xc_photons%pxlda_kappa)
153
154 end if
155
156 !%Variable PhotonXCEnergyMethod
157 !%Type integer
158 !%Default 1
159 !%Section Hamiltonian::XC
160 !%Description
161 !% There are different ways to calculate the energy,
162 !%Option virial 1
163 !% (modified) virial approach</br>
164 !% <math>
165 !% (E_{\rm{px}}^{\rm{virial}} = \frac{1}{2}\int d\mathbf{r}\ \mathbf{r}\cdot[
166 !% -\rho(\mathbf{r})\nabla v_{\rm{px}}(\mathbf{r})])
167 !% </math></br>
168 !%Option expectation_value 2
169 !% expectation value w.tr.t. the wave functions (valid only for 1 electron)</br>
170 !% <math>
171 !% E_{\rm{px}}[\rho] = -\sum_{\alpha=1}^{M_{p}}\frac{\tilde{\lambda}_{\alpha}^{2}}{2\tilde{\omega}_{\alpha}^{2}}
172 !% \langle (\tilde{\mathbf{{\varepsilon}}}_{\alpha}\cdot\hat{\mathbf{J}}_{\rm{p}})\Phi[\rho]
173 !% | (\tilde{\mathbf{{\varepsilon}}}_{\alpha}\cdot\hat{\mathbf{J}}_{\rm{p}})\Phi[\rho] \rangle
174 !% </math></br>
175 !% This option only works for the wave function based electron-photon functionals
176 !%Option LDA 3
177 !% energy from electron density</br>
178 !% <math>
179 !% E_{\rm pxLDA}[\rho] = \frac{-2\pi^{2}}{(d+2)({2V_{d}})^{\frac{2}{d}}}
180 !% \sum_{\alpha=1}^{M_{p}}\frac{\tilde{\lambda}_{\alpha}^{2}}{\tilde{\omega}_{\alpha}^{2}}
181 !% \int d\mathbf{r}\ \rho^{\frac{2+d}{d}}(\mathbf{r})
182 !% </math></br>
183 !% This option only works with LDA electron-photon functionals.
184 !%End
186 call parse_variable(namespace, 'PhotonXCEnergyMethod', 1, xc_photons%energy_method)
188 if( xc_photons%method == xc_photons_wfs .and. xc_photons%energy_method == option__photonxcenergymethod__lda ) then
189 message(1) = "Calculating the electron-photon energy from the LDA expression"
190 message(2) = "is not implemented for wave function based electron-photon functionals"
191 call messages_fatal(2, namespace=namespace)
192 end if
193
194
195 if (xc_photons%lcorrelations) then
196
197 !%Variable PhotonXCEtaC
198 !%Type float
199 !%Default 1.0
200 !%Section Hamiltonian::XC
201 !%Description
202 !% The scaling factor for the px potential to reduce the weak coupling perturbation regime
203 !%End
204
205 if (parse_is_defined(namespace, 'PhotonXCEtaC')) then
206 call parse_variable(namespace, 'PhotonXCEtaC', m_one, xc_photons%eta_c)
207 else
208 message(1) = "Defining PhotonXCEtaC is required for photon functionals containing correlation."
209 call messages_fatal(1, namespace=namespace)
210 end if
211
212 else
213
214 xc_photons%eta_c = m_one
215
216 end if
217
218 ! This variable will keep vpx across iterations
219 safe_allocate(xc_photons%vpx(1:gr%np_part))
220
221 xc_photons%vpx = m_zero
222
223 !%Variable PhotonXCLambShift
224 !%Type logical
225 !%Default .false.
226 !%Section Hamiltonian::XC
227 !%Description
228 !% to deal with the photon free exchange potential for continuum mode in free space
229 !%End
230
231 call parse_variable(namespace, 'PhotonXCLambShift', .false., xc_photons%llamb_freespace)
232 call messages_experimental("PhotonXCLambShift", namespace=namespace)
233
234 if (xc_photons%llamb_freespace) then
235
236 !%Variable PhotonXCLambShiftOmegaCutoff
237 !%Type float
238 !%Default 0.0
239 !%Section Hamiltonian::XC
240 !%Description
241 !% the cutoff frequency (Ha) for Lamb shift
242 !%End
243
244 call parse_variable(namespace, 'PhotonXCLambShiftOmegaCutoff', m_zero, xc_photons%lamb_omega)
245
246 !%Variable PhotonXCLambShiftRenormalizeMass
247 !%Type logical
248 !%Default .false.
249 !%Section Hamiltonian::XC
250 !%Description
251 !% to deal with the photon free exchange potential for continuum mode in free space
252 !%End
253
254 call parse_variable(namespace, 'PhotonXCLambShiftRenormalizeMass', .false., xc_photons%llamb_re_mass)
255
256 end if
257
258 ! compute the dressed photon modes
259 call photon_mode_dressed(xc_photons%pt)
260
261
262 pop_sub(xc_photons_init)
263
264 end subroutine xc_photons_init
265
266 ! ---------------------------------------------------------
267 subroutine xc_photons_end(this)
268 class(xc_photons_t), intent(inout) :: this
269
270 push_sub(xc_photons_end)
271
272 call photon_mode_end(this%pt)
273 safe_deallocate_a(this%vpx)
274
275 if (allocated(this%mf_vector_potential)) then
276 safe_deallocate_a(this%mf_vector_potential)
277 end if
278 if (allocated(this%jp_proj_eo)) then
279 safe_deallocate_a(this%jp_proj_eo)
280 end if
281
282 pop_sub(xc_photons_end)
283 end subroutine xc_photons_end
284
290 !
291 subroutine xc_photons_v_ks(xc_photons, namespace, total_density, gr, space, psolver, ep, st)
292 class(xc_photons_t), intent(inout) :: xc_photons
293 type(namespace_t), intent(in) :: namespace
294 real(real64), pointer, contiguous, intent(in) :: total_density(:)
295 class(grid_t), intent(in) :: gr
296 type(space_t), intent(in) :: space
297 type(poisson_t), intent(in) :: psolver
298 type(epot_t), intent(in) :: ep
299 type(states_elec_t), intent(inout) :: st
300
301 integer :: ia
302
303 push_sub(xc_photons_v_ks)
304
305 xc_photons%lpfmf = xc_photons%method > 0
306
307 xc_photons%vpx = m_zero
308 xc_photons%ex = m_zero
309
310 if ( .not. allocated(xc_photons%mf_vector_potential) ) then
311 safe_allocate(xc_photons%mf_vector_potential(1:space%dim))
312 xc_photons%mf_vector_potential = m_zero
313 end if
314 if ( .not. allocated(xc_photons%jp_proj_eo)) then
315 safe_allocate(xc_photons%jp_proj_eo(1:xc_photons%pt%nmodes,1:2))
316 xc_photons%jp_proj_eo = m_zero
317 end if
318
319
320 select case(xc_photons%method)
321 case(xc_photons_lda) ! LDA approximation for px potential
322 call photon_free_vpx_lda(namespace, xc_photons, total_density, gr, space, psolver)
323 case(xc_photons_wfs) ! wave function approxmation for px potential
324 call photon_free_vpx_wfc(namespace, xc_photons, total_density, gr, space, st)
325 case(xc_photons_none) ! no photon-exchange potential
326 call messages_write('Photon-free px potential is not computed', new_line = .true.)
327 call messages_info()
328 case default
329 assert(.false.)
330 end select
331
332
333 if (.not. xc_photons%llamb_freespace) then
334 ! add the constant energy shift
335 do ia = 1, xc_photons%pt%nmodes
336 xc_photons%ex = xc_photons%ex + 0.5_real64 * (xc_photons%pt%dressed_omega(ia)-xc_photons%pt%omega(ia))
337 end do
338 end if
339
340 pop_sub(xc_photons_v_ks)
341 end subroutine xc_photons_v_ks
342 ! ---------------------------------------------------------
343
344 ! ---------------------------------------------------------
345 ! ---------------------------------------------------------
346 !
378 ! The Lamb shift code is experimental and untested
379 subroutine photon_free_vpx_lda(namespace, xc_photons, total_density, gr, space, psolver)
380 type(namespace_t), intent(in) :: namespace
381 type(xc_photons_t), intent(inout) :: xc_photons
382 real(real64), pointer, contiguous, intent(in) :: total_density(:)
383 type(grid_t), target, intent(in) :: gr
384 type(space_t), intent(in) :: space
385 type(poisson_t), intent(in) :: psolver
387 integer :: ia, ip, iter
388 real(real64) :: unit_volume, r, res, presum, prefact
389 real(real64) :: xx(space%dim), prefactor_lamb
390 real(real64), allocatable :: prefactor(:)
391 real(real64), allocatable :: rho_aux(:)
392 real(real64), allocatable :: grad_rho_aux(:,:)
393 real(real64), allocatable :: px_source(:)
394 real(real64), allocatable :: tmp1(:)
395 real(real64), allocatable :: tmp2(:,:)
396 real(real64), allocatable :: tmp3(:)
397 real(real64), allocatable :: grad_vpx(:,:)
398 real(real64), allocatable :: epsgrad_epsgrad_rho_aux(:)
399 real(real64), allocatable :: epx_force_module(:)
400
401 real(real64), parameter :: threshold = 1e-7_real64
402
403 push_sub(photon_free_vpx_lda)
404
405 if (xc_photons%energy_method == 2 .and. xc_photons%pt%n_electrons >1) then
406 call messages_not_implemented("expectation value for energy for pxLDA more than 1 electron", namespace=namespace)
407 end if
408
409 xc_photons%vpx = m_zero
410 xc_photons%ex = m_zero
411
412 safe_allocate(prefactor(1:xc_photons%pt%nmodes))
413 prefactor = m_zero
414 ! here we will use only one spin channel
415 safe_allocate(rho_aux(1:gr%np_part))
416 safe_allocate(grad_rho_aux(1:gr%np, 1:xc_photons%pt%dim))
417 safe_allocate(px_source(1:gr%np_part))
418 safe_allocate(tmp1(1:gr%np_part))
419 safe_allocate(tmp2(1:gr%np, 1:xc_photons%pt%dim))
420 safe_allocate(tmp3(1:gr%np_part))
421 safe_allocate(grad_vpx(1:gr%np, 1:xc_photons%pt%dim))
422 grad_vpx = m_zero
423 safe_allocate(epx_force_module(1:gr%np_part))
424 epx_force_module = m_zero
425
426 select case(xc_photons%pt%dim)
427 case(1)
428 unit_volume = m_two
429 case(2)
430 unit_volume = m_pi
431 case(3)
432 unit_volume = m_four*m_pi/m_three
433 case default
434 call messages_not_implemented("LDA px more than 3 dimension", namespace=namespace)
435 end select
436
437 !$omp parallel do
438 do ip=1, gr%np
439 rho_aux(ip) = ( abs(total_density(ip))/(m_two*unit_volume) )**(m_two/(xc_photons%pt%dim*m_one))
440 end do
441 !$omp end parallel do
442
443
444 ! compute the electron-photon exchange potential
445
446 if (xc_photons%llamb_freespace) then
447
448 ! Note: The Lamb shift part is currently experimental and untested!
449 ! compute the electron-photon exchange potential
450
451 prefactor_lamb = -(8.0_real64*m_pi*m_third) * xc_photons%lamb_omega / (p_c**3)
452
453 !$OMP parallel do
454 do ip=1,gr%np
455 xc_photons%vpx(ip) = prefactor_lamb*rho_aux(ip)
456 end do
457 !$OMP end parallel do
458
459 else
460
461 do ia = 1, xc_photons%pt%nmodes
462 prefactor(ia) = -m_two*(m_pi * xc_photons%pt%dressed_lambda(ia) / xc_photons%pt%dressed_omega(ia))**2
463 end do
464
465 select case(xc_photons%pt%dim)
466 case(1)
467 ! solve the pxLDA potential using the analytical form in 1D
468 px_source = m_zero
469 do ia = 1, xc_photons%pt%nmodes
470 !$OMP parallel do
471 do ip=1,gr%np_part
472 px_source(ip) = px_source(ip) + prefactor(ia)
473 end do
474 !$OMP end parallel do
475 end do
476
477 !$OMP parallel do
478 do ip=1,gr%np
479 xc_photons%vpx(ip) = px_source(ip)*rho_aux(ip)
480 end do
481 !$OMP end parallel do
482 case(2)
483 call get_px_source(px_source)
484 ! for 2D we solve the Poisson equation using the conjugate gradient method
485
486 ! Note that we need to solve -\Delta f = v, as CG requires SPD operator
487 call lalg_scal(gr%np, -m_one, px_source)
488
489 mesh_aux => gr%der%mesh
490 iter = 1000
491 call dconjugate_gradients(gr%np, xc_photons%vpx(:), px_source, laplacian_op, dmf_dotp_aux, iter, res, threshold, &
492 userdata=[c_loc(gr)])
493 write(message(1),'(a,i6,a)') "Info: CG converged with ", iter, " iterations."
494 write(message(2),'(a,e14.6)') "Info: The residue is ", res
495 call messages_info(2, namespace=namespace)
496
497 case(3)
498 ! for 3D we use thepoisson solver including the prefactor (-4*pi)
499 ! therefore, we need to remove the factor in advance
500 call get_px_source(px_source)
501
502 call lalg_scal(gr%np, m_one/(-m_four*m_pi), px_source)
503
504 ! solve the Poisson equation
505 call dpoisson_solve(psolver, namespace, xc_photons%vpx(:), px_source)
506 case default
507 assert(.false.)
508 end select
509
510 end if
511
512 ! scaling the potential
513 call lalg_scal(gr%np, (xc_photons%eta_c * xc_photons%pxlda_kappa), xc_photons%vpx)
514
515 ! compute electron-photon energy
516
517 select case (xc_photons%energy_method)
518 case(1) ! compute the epx energy from the virial relation
519
520 do ia = 1, xc_photons%pt%nmodes
521
522 ! compute the electron-photon force
523 !$omp parallel do
524 do ip = 1, gr%np
525 epx_force_module(ip) = -prefactor(ia)*m_two*abs(total_density(ip))*rho_aux(ip)/(xc_photons%pt%dim*m_one+m_two)
526 end do
527 !$omp end parallel do
528
529 call dderivatives_grad(gr%der, epx_force_module(:), tmp2)
530 call lalg_gemv(gr%np, xc_photons%pt%dim, m_one, tmp2, xc_photons%pt%dressed_pol(1:xc_photons%pt%dim, ia), m_zero, tmp1)
531
532 !$omp parallel do private(r, xx)
533 do ip = 1, gr%np
534 call mesh_r(gr, ip, r, coords=xx)
535 tmp3(ip) = tmp1(ip)*dot_product(xx(1:xc_photons%pt%dim), xc_photons%pt%dressed_pol(1:xc_photons%pt%dim, ia))
536 end do
537 !$omp end parallel do
538
539 xc_photons%ex = xc_photons%ex + m_half*dmf_integrate(gr, tmp3)
540 end do
541
542 xc_photons%ex = xc_photons%eta_c * xc_photons%pxlda_kappa * xc_photons%ex
543
544 case(2) ! compute the energy as expetation value wrt to the wave functions
545
546 rho_aux(1:gr%np) = sqrt(abs( total_density(1:gr%np)))
547
548 safe_allocate(epsgrad_epsgrad_rho_aux(1:gr%np))
549
550 call dderivatives_grad(gr%der, rho_aux(1:gr%np_part), grad_rho_aux)
551
552 do ia = 1, xc_photons%pt%nmodes
553 prefact = (xc_photons%pt%dressed_lambda(ia)**2) / (m_two*xc_photons%pt%dressed_omega(ia)**2)
554
555 call lalg_gemv(gr%np, xc_photons%pt%dim, m_one, grad_rho_aux, xc_photons%pt%dressed_pol(:, ia), m_zero, tmp1)
556 call dderivatives_grad(gr%der, tmp1, tmp2)
557 call lalg_gemv(gr%np, xc_photons%pt%dim, m_one, tmp2, xc_photons%pt%dressed_pol(1:xc_photons%pt%dim, ia), m_zero, tmp1)
558
559 !$OMP parallel do
560 do ip=1, gr%np
561 epsgrad_epsgrad_rho_aux(ip) = epsgrad_epsgrad_rho_aux(ip) + prefact*tmp1(ip)
562 end do
563 !$OMP end parallel do
564
565 end do
566
567 xc_photons%ex = xc_photons%eta_c * dmf_dotp(gr, rho_aux(1:gr%np), epsgrad_epsgrad_rho_aux(1:gr%np))
568
569 safe_deallocate_a(epsgrad_epsgrad_rho_aux)
570
571 case(3) ! integrate the aux electron density over the volume
572
573 !$OMP parallel do
574 do ip=1,gr%np
575 tmp1(ip) = abs( total_density(ip))**((m_one*xc_photons%pt%dim+m_two)/(xc_photons%pt%dim*m_one))
576 end do
577 !$OMP end parallel do
578
579 ! sum over the prefactors
580 presum = m_zero
581 do ia = 1, xc_photons%pt%nmodes
582 presum = presum + prefactor(ia)
583 end do
584 presum = presum * (m_one/(m_two*unit_volume))**(m_two/(xc_photons%pt%dim*m_one)) / (xc_photons%pt%dim*m_one+m_two)
585 presum = presum * xc_photons%eta_c * xc_photons%pxlda_kappa
586
587 xc_photons%ex = xc_photons%eta_c * xc_photons%pxlda_kappa * presum * dmf_integrate(gr, tmp1)
588
589 end select
590
591 safe_deallocate_a(prefactor)
592 safe_deallocate_a(rho_aux)
593 safe_deallocate_a(grad_rho_aux)
594 safe_deallocate_a(px_source)
595 safe_deallocate_a(tmp1)
596 safe_deallocate_a(tmp2)
597 safe_deallocate_a(tmp3)
598 safe_deallocate_a(grad_vpx)
599 safe_deallocate_a(epx_force_module)
600
601 pop_sub(photon_free_vpx_lda)
602
603 contains
604
605 subroutine get_px_source(px_source)
606 real(real64), contiguous, intent(out) :: px_source(:)
607
608 call dderivatives_grad(gr%der, rho_aux(1:gr%np_part), grad_rho_aux)
609
610 px_source = m_zero
611 do ia = 1, xc_photons%pt%nmodes
612 call lalg_gemv(gr%np, xc_photons%pt%dim, m_one, grad_rho_aux, xc_photons%pt%dressed_pol(:, ia), m_zero, tmp1)
613 call dderivatives_grad(gr%der, tmp1, tmp2)
614 call lalg_gemv(gr%np, xc_photons%pt%dim, m_one, tmp2, xc_photons%pt%dressed_pol(1:xc_photons%pt%dim, ia), m_zero, tmp1)
615 call lalg_axpy(gr%np, prefactor(ia), tmp1, px_source)
616 end do
617
618 end subroutine get_px_source
619
620 end subroutine photon_free_vpx_lda
621 ! ---------------------------------------------------------
622
623 ! ---------------------------------------------------------
625 subroutine laplacian_op(x, hx, userdata)
626 real(real64), contiguous, intent(in) :: x(:)
627 real(real64), contiguous, intent(out) :: Hx(:)
628 type(c_ptr), intent(in) :: userdata(:)
629
630 real(real64), allocatable :: tmpx(:)
631 type(grid_t), pointer :: gr
632
633 assert(size(userdata) == 1)
634 assert(c_associated(userdata(1)))
635 call c_f_pointer(userdata(1), gr)
636
637 safe_allocate(tmpx(1:gr%np_part))
638 call lalg_copy(gr%np, x, tmpx)
639 call dderivatives_lapl(gr%der, tmpx, hx, factor=-m_one)
640 safe_deallocate_a(tmpx)
641
642 end subroutine laplacian_op
643 ! ---------------------------------------------------------
644
645 ! ---------------------------------------------------------
646 !
672 !
673 ! The Lamb shift code is experimental and untested
674
675 subroutine photon_free_vpx_wfc(namespace, xc_photons, total_density, gr, space, st)
676 type(namespace_t), intent(in) :: namespace
677 type(xc_photons_t), intent(inout) :: xc_photons
678 real(real64), pointer, contiguous, intent(in) :: total_density(:)
679 class(grid_t), intent(in) :: gr
680 type(space_t), intent(in) :: space
681 type(states_elec_t), intent(in) :: st
682
683 integer :: ia, ip
684 real(real64) :: prefactor_lamb
685 real(real64) :: xx(space%dim), r
686 real(real64), allocatable :: prefactor(:)
687 real(real64), allocatable :: rho_aux(:)
688 real(real64), allocatable :: grad_rho_aux(:,:)
689 real(real64), allocatable :: grad_vpx(:,:)
690 real(real64), allocatable :: epsgrad_epsgrad_rho_aux(:)
691 real(real64), allocatable :: tmp1(:)
692 real(real64), allocatable :: tmp2(:,:)
693 real(real64) :: shift
694
695 push_sub(photon_free_vpx_wfc)
696
697 if (st%d%nspin >1) then
698 call messages_not_implemented("PhotonXCXCMethod = wavefunction for polarized and spinor cases", namespace=namespace)
699 end if
701 if (xc_photons%pt%n_electrons >1) then
702 call messages_not_implemented("PhotonXCXCMethod = wavefunction for more than 1 electron", namespace=namespace)
703 end if
704
705 xc_photons%vpx = m_zero
706 xc_photons%ex = m_zero
707
708 safe_allocate(prefactor(1:xc_photons%pt%nmodes))
709 prefactor = m_zero
710 safe_allocate(rho_aux(1:gr%np_part))
711 rho_aux = m_zero
712 safe_allocate(grad_rho_aux(1:gr%np, 1:xc_photons%pt%dim))
713 grad_rho_aux = m_zero
714 safe_allocate(epsgrad_epsgrad_rho_aux(1:gr%np))
715 safe_allocate(grad_vpx(1:gr%np, 1:xc_photons%pt%dim))
716 safe_allocate(tmp1(1:gr%np_part))
717 safe_allocate(tmp2(1:gr%np_part, 1:xc_photons%pt%dim))
718
719
720 rho_aux(1:gr%np) = sqrt(abs( total_density(1:gr%np)))
721 !$OMP parallel do
722 do ip = 1, gr%np
723 rho_aux(ip) = safe_tol(rho_aux(ip),1e-18_real64)
724 end do
725 !$OMP end parallel do
726
727 if (xc_photons%llamb_freespace) then
728
729 ! experimental Lamb shift mode
730
731 prefactor_lamb = ( m_two/(m_three*m_pi) ) * xc_photons%lamb_omega / p_c**3
732 call dderivatives_lapl(gr%der, rho_aux(1:gr%np_part), epsgrad_epsgrad_rho_aux)
733 call lalg_scal(gr%np, prefactor_lamb, epsgrad_epsgrad_rho_aux)
734
735 else
736
737 do ia = 1, xc_photons%pt%nmodes
738 prefactor(ia) = (xc_photons%pt%dressed_lambda(ia)**2) / (m_two*xc_photons%pt%dressed_omega(ia)**2)
739 end do
740
741 call dderivatives_grad(gr%der, rho_aux, grad_rho_aux)
742
743 epsgrad_epsgrad_rho_aux = m_zero
744 do ia = 1, xc_photons%pt%nmodes
745 tmp1 = m_zero
746 call lalg_gemv(gr%np, xc_photons%pt%dim, m_one, grad_rho_aux, xc_photons%pt%dressed_pol(:, ia), m_zero, tmp1)
747 call dderivatives_grad(gr%der, tmp1, tmp2)
748 call lalg_gemv(gr%np, xc_photons%pt%dim, prefactor(ia), tmp2, xc_photons%pt%dressed_pol(:, ia), &
749 m_one, epsgrad_epsgrad_rho_aux)
750 end do
751
752 end if
753
754 !$OMP parallel do
755 do ip = 1, gr%np
756 xc_photons%vpx(ip) = epsgrad_epsgrad_rho_aux(ip)/rho_aux(ip)
757 end do
758 !$OMP end parallel do
759
760 if(st%eigenval(1,1) < m_huge .and. .not. space%is_periodic()) then
761 shift = m_two * st%eigenval(1,1) * prefactor(1)
762 !$OMP parallel do
763 do ip=1, gr%np
764 xc_photons%vpx(ip) = xc_photons%vpx(ip) + shift
765 end do
766 !$OMP end parallel do
767 end if
768
769 call lalg_scal(gr%np, xc_photons%eta_c, xc_photons%vpx(:))
771 select case (xc_photons%energy_method)
772 case(1) ! virial
773
774 call dderivatives_grad(gr%der, xc_photons%vpx(:), grad_vpx)
775 !$OMP parallel do private(r, xx)
776 do ip = 1, gr%np
777 call mesh_r(gr, ip, r, coords=xx)
778 tmp1(ip) = - total_density(ip)*dot_product(xx(1:xc_photons%pt%dim), grad_vpx(ip,1:xc_photons%pt%dim))
779 end do
780 !$OMP end parallel do
781
782 xc_photons%ex = m_half*dmf_integrate(gr, tmp1)
783
784 case(2) ! expectation_value
785 xc_photons%ex = xc_photons%eta_c * dmf_dotp(gr, rho_aux(1:gr%np), epsgrad_epsgrad_rho_aux)
786
787 case default
788 assert(.false.)
789 end select
790
791 safe_deallocate_a(prefactor)
792 safe_deallocate_a(rho_aux)
793 safe_deallocate_a(grad_rho_aux)
794 safe_deallocate_a(grad_vpx)
795 safe_deallocate_a(epsgrad_epsgrad_rho_aux)
796 safe_deallocate_a(tmp1)
797 safe_deallocate_a(tmp2)
798
799 pop_sub(photon_free_vpx_wfc)
800
801 end subroutine photon_free_vpx_wfc
802 ! ---------------------------------------------------------
803
804
805 ! ---------------------------------------------------------
815 !
816 subroutine xc_photons_add_mean_field(xc_photons, gr, space, kpoints, st, time, dt)
817 class(xc_photons_t), intent(inout) :: xc_photons
818 class(grid_t), intent(in) :: gr
819 class(space_t), intent(in) :: space
820 type(kpoints_t), intent(in) :: kpoints
821 type(states_elec_t), intent(inout) :: st
822 real(real64), intent(in) :: time
823 real(real64), intent(in) :: dt
824
825
826 integer :: ia, idir, ispin
827 real(real64) :: pol_dot_jp
828 real(real64), allocatable :: current(:,:,:)
829 real(real64), allocatable :: jp(:)
830
832
833 ! compute the dressed photon-free vector potential
834 xc_photons%mf_vector_potential = m_zero
835
836 safe_allocate(current(1:gr%np_part, 1:space%dim, 1:st%d%nspin))
837 current = m_zero
838 ! here we use the paramagnetic current; note that the physical current here
839 ! only contains the paramagnetic current.
840 call states_elec_calc_quantities(gr, st, kpoints, .false., paramagnetic_current = current)
841
842 ! compute the paramagnetic current
843 safe_allocate(jp(1:space%dim))
844 do idir = 1, space%dim
845 jp(idir) = m_zero
846 do ispin = 1, st%d%spin_channels
847 jp(idir) = jp(idir) + dmf_integrate(gr%der%mesh, current(:,idir,ispin))
848 end do
849 end do
850
851 ! update the 'projected' current
852 do ia = 1, xc_photons%pt%nmodes
853 pol_dot_jp = dot_product(xc_photons%pt%dressed_pol(1:space%dim, ia),jp(1:space%dim))
854 xc_photons%jp_proj_eo(ia,1) = xc_photons%jp_proj_eo(ia,1) + &
855 cos(xc_photons%pt%dressed_omega(ia)*( time-dt))*pol_dot_jp*dt
856 xc_photons%jp_proj_eo(ia,2) = xc_photons%jp_proj_eo(ia,2) + &
857 sin(xc_photons%pt%dressed_omega(ia)*( time-dt))*pol_dot_jp*dt
858 end do
859
860 do ia = 1, xc_photons%pt%nmodes
861 xc_photons%mf_vector_potential(1:xc_photons%pt%dim) = xc_photons%mf_vector_potential(1:xc_photons%pt%dim) &
862 + (-p_c*(xc_photons%pt%dressed_lambda(ia)**2) / xc_photons%pt%dressed_omega(ia)) &
863 * (xc_photons%jp_proj_eo(ia,1)*sin(xc_photons%pt%dressed_omega(ia)* time) &
864 - xc_photons%jp_proj_eo(ia,2)*cos(xc_photons%pt%dressed_omega(ia)* time)) &
865 * xc_photons%pt%dressed_pol(1:xc_photons%pt%dim, ia)
866 end do
867
868 safe_deallocate_a(current)
869 safe_deallocate_a(jp)
870
872
873 end subroutine xc_photons_add_mean_field
874 ! ---------------------------------------------------------
875
876
880 !
881 logical pure function xc_photons_wants_to_renormalize_mass(xc_photons) result (renorm)
882 class(xc_photons_t), intent(in) :: xc_photons
883
884 renorm = (xc_photons%method>0) .and. xc_photons%llamb_freespace .and. xc_photons%llamb_re_mass
885
887
889 real(real64) pure function xc_photons_get_renormalized_emass(xc_photons) result(mass)
890 class(xc_photons_t), intent(in) :: xc_photons
891
892 mass = m_one - (m_four*xc_photons%lamb_omega) / (3.0*m_pi * p_c**3)
893
895
897 !
898 subroutine xc_photons_mf_dump(xc_photons, restart, ierr)
899 class(xc_photons_t), intent(in) :: xc_photons
900 type(restart_t), intent(in) :: restart
901 integer, intent(out) :: ierr
902
903 character(len=80), allocatable :: lines(:)
904 integer :: iunit, err, jj, nmodes
905 integer :: pt_dim
906
907 push_sub(photon_free_mf_dump)
908 nmodes = xc_photons%pt%nmodes
909 pt_dim = xc_photons%pt%dim
910
911 safe_allocate(lines(1:nmodes+pt_dim))
912
913 ierr = 0
914
915 iunit = restart%open('photon_free_mf')
916
917 do jj = 1, pt_dim
918 write(lines(jj), '(2x, es19.12)') xc_photons%mf_vector_potential(jj)
919 end do
920
921 do jj = 1, nmodes
922 write(lines(jj+pt_dim), '(a10,1x,I8,a1,2x,2(es19.12,2x))') 'Mode ', jj, ":", xc_photons%jp_proj_eo(jj,1:2)
923 end do
924
925 call restart%write(iunit, lines, nmodes+pt_dim, err)
926 if (err /= 0) ierr = ierr + 1
927 call restart%close(iunit)
928
929 safe_deallocate_a(lines)
930
931 pop_sub(photon_free_mf_dump)
932 end subroutine xc_photons_mf_dump
933
934! ---------------------------------------------------------
935
937 !
938 subroutine xc_photons_mf_load(xc_photons, restart, space, ierr)
939 class(xc_photons_t), intent(inout) :: xc_photons
940 type(restart_t), intent(in) :: restart
941 class(space_t), intent(in) :: space
942 integer, intent(out) :: ierr
943
944 character(len=80), allocatable :: lines(:)
945 character(len=7) :: sdummy
946 integer :: idummy
947 integer :: iunit, err, jj, nmodes
948 integer :: pt_dim
949
950 push_sub(photon_free_mf_load)
951
952 ierr = 0
953 nmodes = xc_photons%pt%nmodes
954 pt_dim = xc_photons%pt%dim
955
956 if (restart%skip()) then
957 ierr = -1
958 pop_sub(photon_free_mf_load)
959 return
960 end if
961
962 message(1) = "Debug: Reading Photon-Free Photons restart."
963 call messages_info(1, namespace=restart%namespace, debug_only=.true.)
964
965 if ( .not. allocated(xc_photons%jp_proj_eo)) then
966 safe_allocate(xc_photons%jp_proj_eo(1:xc_photons%pt%nmodes, 1:2))
967 xc_photons%jp_proj_eo = m_zero
968 end if
969 if ( .not. allocated(xc_photons%mf_vector_potential)) then
970 safe_allocate(xc_photons%mf_vector_potential(1:space%dim))
971 xc_photons%mf_vector_potential = m_zero
972 end if
973
974 safe_allocate(lines(1:nmodes+pt_dim))
975 iunit = restart%open('photon_free_mf')
976 call restart%read(iunit, lines, nmodes+pt_dim, err)
977 if (err /= 0) then
978 ierr = ierr + 1
979 else
980 do jj = 1, pt_dim
981 read(lines(jj),'(2x, es19.12)') xc_photons%mf_vector_potential(jj)
982 end do
983
984 do jj = 1, nmodes
985 read(lines(jj+pt_dim), '(a10,1x,I8,a1,2x,2(es19.12,2x))') sdummy, idummy, sdummy, xc_photons%jp_proj_eo(jj,1:2)
986 end do
987 end if
988 call restart%close(iunit)
989
990 message(1) = "Debug: Reading Photons restart done."
991 call messages_info(1, namespace=restart%namespace, debug_only=.true.)
992
993 safe_deallocate_a(lines)
994
995 pop_sub(photon_free_mf_load)
996 end subroutine xc_photons_mf_load
997
998end module xc_photons_oct_m
999
1000!! Local Variables:
1001!! mode: f90
1002!! coding: utf-8
1003!! End:
constant times a vector plus a vector
Definition: lalg_basic.F90:173
Copies a vector x, to a vector y.
Definition: lalg_basic.F90:188
scales a vector by a constant
Definition: lalg_basic.F90:159
double sin(double __x) __attribute__((__nothrow__
double sqrt(double __x) __attribute__((__nothrow__
double cos(double __x) __attribute__((__nothrow__
This module calculates the derivatives (gradients, Laplacians, etc.) of a function.
subroutine, public dderivatives_grad(der, ff, op_ff, ghost_update, set_bc, to_cartesian)
apply the gradient to a mesh function
subroutine, public dderivatives_lapl(der, ff, op_ff, ghost_update, set_bc, factor)
apply the Laplacian to a mesh function
real(real64), parameter, public m_two
Definition: global.F90:192
real(real64), parameter, public m_huge
Definition: global.F90:208
real(real64), parameter, public m_zero
Definition: global.F90:190
real(real64), parameter, public m_four
Definition: global.F90:194
real(real64), parameter, public m_third
Definition: global.F90:197
real(real64), parameter, public m_pi
some mathematical constants
Definition: global.F90:188
real(real64), parameter, public m_half
Definition: global.F90:196
real(real64), parameter, public p_c
Electron gyromagnetic ratio, see Phys. Rev. Lett. 130, 071801 (2023)
Definition: global.F90:228
real(real64), parameter, public m_one
Definition: global.F90:191
real(real64), parameter, public m_three
Definition: global.F90:193
This module implements the underlying real-space grid.
Definition: grid.F90:119
This module defines various routines, operating on mesh functions.
class(mesh_t), pointer, public mesh_aux
Globally-scoped pointer to the mesh instance.
real(real64) function, public dmf_dotp_aux(f1, f2)
dot product between two vectors (mesh functions)
This module defines the meshes, which are used in Octopus.
Definition: mesh.F90:120
pure subroutine, public mesh_r(mesh, ip, rr, origin, coords)
return the distance to the origin for a given grid point
Definition: mesh.F90:339
subroutine, public messages_not_implemented(feature, namespace)
Definition: messages.F90:1097
character(len=256), dimension(max_lines), public message
to be output by fatal, warning
Definition: messages.F90:162
subroutine, public messages_fatal(no_lines, only_root_writes, namespace)
Definition: messages.F90:416
subroutine, public messages_experimental(name, namespace)
Definition: messages.F90:1069
subroutine, public messages_info(no_lines, iunit, debug_only, stress, all_nodes, namespace)
Definition: messages.F90:600
logical function, public parse_is_defined(namespace, name)
Definition: parser.F90:505
subroutine, public photon_mode_end(this)
subroutine, public photon_mode_dressed(this)
subroutine, public photon_mode_set_n_electrons(this, qtot)
subroutine, public photon_mode_init(this, namespace, dim, photon_free)
subroutine, public dpoisson_solve(this, namespace, pot, rho, all_nodes, kernel, reset)
Calculates the Poisson equation. Given the density returns the corresponding potential.
Definition: poisson.F90:871
This module is intended to contain "only mathematical" functions and procedures.
Definition: solvers.F90:117
subroutine, public states_elec_calc_quantities(gr, st, kpoints, nlcc, kinetic_energy_density, paramagnetic_current, density_gradient, density_laplacian, gi_kinetic_energy_density, st_end)
calculated selected quantities
This module implements the "photon-free" electron-photon exchange-correlation functional.
Definition: xc_photons.F90:123
subroutine laplacian_op(x, hx, userdata)
Computes Hx = (\Laplacian) x.
Definition: xc_photons.F90:721
subroutine photon_free_vpx_wfc(namespace, xc_photons, total_density, gr, space, st)
compute the electron-photon exchange potential based on wave functions
Definition: xc_photons.F90:771
logical pure function xc_photons_wants_to_renormalize_mass(xc_photons)
indicate whether the photon-exchange requires a renormalized electron mass
Definition: xc_photons.F90:977
subroutine xc_photons_init(xc_photons, namespace, xc_photon, space, gr, st)
initialize the photon-exchange functional
Definition: xc_photons.F90:203
subroutine xc_photons_add_mean_field(xc_photons, gr, space, kpoints, st, time, dt)
accumulate the results of time integral the paramagnetic current.
Definition: xc_photons.F90:912
subroutine xc_photons_mf_dump(xc_photons, restart, ierr)
write restart information
Definition: xc_photons.F90:994
subroutine xc_photons_v_ks(xc_photons, namespace, total_density, gr, space, psolver, ep, st)
evaluate the KS potential and energy for the given functional
Definition: xc_photons.F90:387
real(real64) pure function xc_photons_get_renormalized_emass(xc_photons)
return the renormalized electron mass for the electron-photon exhange
Definition: xc_photons.F90:985
integer, parameter, private xc_photons_lda
Definition: xc_photons.F90:192
subroutine photon_free_vpx_lda(namespace, xc_photons, total_density, gr, space, psolver)
compute the electron-photon exchange potential within the LDA
Definition: xc_photons.F90:475
subroutine xc_photons_end(this)
Definition: xc_photons.F90:363
integer, parameter, private xc_photons_wfs
Definition: xc_photons.F90:192
subroutine xc_photons_mf_load(xc_photons, restart, space, ierr)
load restart information
Description of the grid, containing information on derivatives, stencil, and symmetries.
Definition: grid.F90:171
The states_elec_t class contains all electronic wave functions.
This class described the 'photon-exchange' electron-photon xc functionals, based on QEDFT.
Definition: xc_photons.F90:159
int true(void)
subroutine get_px_source(px_source)
Definition: xc_photons.F90:701