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, 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(states_elec_t), intent(inout) :: st
299
300 integer :: ia
301
302 push_sub(xc_photons_v_ks)
303
304 xc_photons%lpfmf = xc_photons%method > 0
305
306 xc_photons%vpx = m_zero
307 xc_photons%ex = m_zero
308
309 if ( .not. allocated(xc_photons%mf_vector_potential) ) then
310 safe_allocate(xc_photons%mf_vector_potential(1:space%dim))
311 xc_photons%mf_vector_potential = m_zero
312 end if
313 if ( .not. allocated(xc_photons%jp_proj_eo)) then
314 safe_allocate(xc_photons%jp_proj_eo(1:xc_photons%pt%nmodes,1:2))
315 xc_photons%jp_proj_eo = m_zero
316 end if
317
318
319 select case(xc_photons%method)
320 case(xc_photons_lda) ! LDA approximation for px potential
321 call photon_free_vpx_lda(namespace, xc_photons, total_density, gr, space, psolver)
322 case(xc_photons_wfs) ! wave function approxmation for px potential
323 call photon_free_vpx_wfc(namespace, xc_photons, total_density, gr, space, st)
324 case(xc_photons_none) ! no photon-exchange potential
325 call messages_write('Photon-free px potential is not computed', new_line = .true.)
326 call messages_info()
327 case default
328 assert(.false.)
329 end select
330
331
332 if (.not. xc_photons%llamb_freespace) then
333 ! add the constant energy shift
334 do ia = 1, xc_photons%pt%nmodes
335 xc_photons%ex = xc_photons%ex + 0.5_real64 * (xc_photons%pt%dressed_omega(ia)-xc_photons%pt%omega(ia))
336 end do
337 end if
338
339 pop_sub(xc_photons_v_ks)
340 end subroutine xc_photons_v_ks
341 ! ---------------------------------------------------------
342
343 ! ---------------------------------------------------------
344 ! ---------------------------------------------------------
345 !
377 ! The Lamb shift code is experimental and untested
378 subroutine photon_free_vpx_lda(namespace, xc_photons, total_density, gr, space, psolver)
379 type(namespace_t), intent(in) :: namespace
380 type(xc_photons_t), intent(inout) :: xc_photons
381 real(real64), pointer, contiguous, intent(in) :: total_density(:)
382 type(grid_t), target, intent(in) :: gr
383 type(space_t), intent(in) :: space
384 type(poisson_t), intent(in) :: psolver
385
386 integer :: ia, ip, iter
387 real(real64) :: unit_volume, r, res, presum, prefact
388 real(real64) :: xx(space%dim), prefactor_lamb
389 real(real64), allocatable :: prefactor(:)
390 real(real64), allocatable :: rho_aux(:)
391 real(real64), allocatable :: grad_rho_aux(:,:)
392 real(real64), allocatable :: px_source(:)
393 real(real64), allocatable :: tmp1(:)
394 real(real64), allocatable :: tmp2(:,:)
395 real(real64), allocatable :: tmp3(:)
396 real(real64), allocatable :: grad_vpx(:,:)
397 real(real64), allocatable :: epsgrad_epsgrad_rho_aux(:)
398 real(real64), allocatable :: epx_force_module(:)
399
400 real(real64), parameter :: threshold = 1e-7_real64
401
402 push_sub(photon_free_vpx_lda)
403
404 if (xc_photons%energy_method == 2 .and. xc_photons%pt%n_electrons >1) then
405 call messages_not_implemented("expectation value for energy for pxLDA more than 1 electron", namespace=namespace)
406 end if
407
408 xc_photons%vpx = m_zero
409 xc_photons%ex = m_zero
410
411 safe_allocate(prefactor(1:xc_photons%pt%nmodes))
412 prefactor = m_zero
413 ! here we will use only one spin channel
414 safe_allocate(rho_aux(1:gr%np_part))
415 safe_allocate(grad_rho_aux(1:gr%np, 1:xc_photons%pt%dim))
416 safe_allocate(px_source(1:gr%np_part))
417 safe_allocate(tmp1(1:gr%np_part))
418 safe_allocate(tmp2(1:gr%np, 1:xc_photons%pt%dim))
419 safe_allocate(tmp3(1:gr%np_part))
420 safe_allocate(grad_vpx(1:gr%np, 1:xc_photons%pt%dim))
421 grad_vpx = m_zero
422 safe_allocate(epx_force_module(1:gr%np_part))
423 epx_force_module = m_zero
424
425 select case(xc_photons%pt%dim)
426 case(1)
427 unit_volume = m_two
428 case(2)
429 unit_volume = m_pi
430 case(3)
431 unit_volume = m_four*m_pi/m_three
432 case default
433 call messages_not_implemented("LDA px more than 3 dimension", namespace=namespace)
434 end select
435
436 !$omp parallel do
437 do ip=1, gr%np
438 rho_aux(ip) = ( abs(total_density(ip))/(m_two*unit_volume) )**(m_two/(xc_photons%pt%dim*m_one))
439 end do
440 !$omp end parallel do
441
442
443 ! compute the electron-photon exchange potential
444
445 if (xc_photons%llamb_freespace) then
446
447 ! Note: The Lamb shift part is currently experimental and untested!
448 ! compute the electron-photon exchange potential
449
450 prefactor_lamb = -(8.0_real64*m_pi*m_third) * xc_photons%lamb_omega / (p_c**3)
451
452 !$OMP parallel do
453 do ip=1,gr%np
454 xc_photons%vpx(ip) = prefactor_lamb*rho_aux(ip)
455 end do
456 !$OMP end parallel do
457
458 else
459
460 do ia = 1, xc_photons%pt%nmodes
461 prefactor(ia) = -m_two*(m_pi * xc_photons%pt%dressed_lambda(ia) / xc_photons%pt%dressed_omega(ia))**2
462 end do
463
464 select case(xc_photons%pt%dim)
465 case(1)
466 ! solve the pxLDA potential using the analytical form in 1D
467 px_source = m_zero
468 do ia = 1, xc_photons%pt%nmodes
469 !$OMP parallel do
470 do ip=1,gr%np_part
471 px_source(ip) = px_source(ip) + prefactor(ia)
472 end do
473 !$OMP end parallel do
474 end do
475
476 !$OMP parallel do
477 do ip=1,gr%np
478 xc_photons%vpx(ip) = px_source(ip)*rho_aux(ip)
479 end do
480 !$OMP end parallel do
481 case(2)
482 call get_px_source(px_source)
483 ! for 2D we solve the Poisson equation using the conjugate gradient method
484
485 ! Note that we need to solve -\Delta f = v, as CG requires SPD operator
486 call lalg_scal(gr%np, -m_one, px_source)
487
488 mesh_aux => gr%der%mesh
489 iter = 1000
490 call dconjugate_gradients(gr%np, xc_photons%vpx(:), px_source, laplacian_op, dmf_dotp_aux, iter, res, threshold, &
491 userdata=[c_loc(gr)])
492 write(message(1),'(a,i6,a)') "Info: CG converged with ", iter, " iterations."
493 write(message(2),'(a,e14.6)') "Info: The residue is ", res
494 call messages_info(2, namespace=namespace)
495
496 case(3)
497 ! for 3D we use thepoisson solver including the prefactor (-4*pi)
498 ! therefore, we need to remove the factor in advance
499 call get_px_source(px_source)
500
501 call lalg_scal(gr%np, m_one/(-m_four*m_pi), px_source)
502
503 ! solve the Poisson equation
504 call dpoisson_solve(psolver, namespace, xc_photons%vpx(:), px_source)
505 case default
506 assert(.false.)
507 end select
508
509 end if
510
511 ! scaling the potential
512 call lalg_scal(gr%np, (xc_photons%eta_c * xc_photons%pxlda_kappa), xc_photons%vpx)
513
514 ! compute electron-photon energy
515
516 select case (xc_photons%energy_method)
517 case(1) ! compute the epx energy from the virial relation
518
519 do ia = 1, xc_photons%pt%nmodes
520
521 ! compute the electron-photon force
522 !$omp parallel do
523 do ip = 1, gr%np
524 epx_force_module(ip) = -prefactor(ia)*m_two*abs(total_density(ip))*rho_aux(ip)/(xc_photons%pt%dim*m_one+m_two)
525 end do
526 !$omp end parallel do
527
528 call dderivatives_grad(gr%der, epx_force_module(:), tmp2)
529 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)
530
531 !$omp parallel do private(r, xx)
532 do ip = 1, gr%np
533 call mesh_r(gr, ip, r, coords=xx)
534 tmp3(ip) = tmp1(ip)*dot_product(xx(1:xc_photons%pt%dim), xc_photons%pt%dressed_pol(1:xc_photons%pt%dim, ia))
535 end do
536 !$omp end parallel do
537
538 xc_photons%ex = xc_photons%ex + m_half*dmf_integrate(gr, tmp3)
539 end do
540
541 xc_photons%ex = xc_photons%eta_c * xc_photons%pxlda_kappa * xc_photons%ex
542
543 case(2) ! compute the energy as expetation value wrt to the wave functions
544
545 rho_aux(1:gr%np) = sqrt(abs( total_density(1:gr%np)))
546
547 safe_allocate(epsgrad_epsgrad_rho_aux(1:gr%np))
548
549 call dderivatives_grad(gr%der, rho_aux(1:gr%np_part), grad_rho_aux)
550
551 do ia = 1, xc_photons%pt%nmodes
552 prefact = (xc_photons%pt%dressed_lambda(ia)**2) / (m_two*xc_photons%pt%dressed_omega(ia)**2)
553
554 call lalg_gemv(gr%np, xc_photons%pt%dim, m_one, grad_rho_aux, xc_photons%pt%dressed_pol(:, ia), m_zero, tmp1)
555 call dderivatives_grad(gr%der, tmp1, tmp2)
556 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)
557
558 !$OMP parallel do
559 do ip=1, gr%np
560 epsgrad_epsgrad_rho_aux(ip) = epsgrad_epsgrad_rho_aux(ip) + prefact*tmp1(ip)
561 end do
562 !$OMP end parallel do
563
564 end do
565
566 xc_photons%ex = xc_photons%eta_c * dmf_dotp(gr, rho_aux(1:gr%np), epsgrad_epsgrad_rho_aux(1:gr%np))
567
568 safe_deallocate_a(epsgrad_epsgrad_rho_aux)
569
570 case(3) ! integrate the aux electron density over the volume
571
572 !$OMP parallel do
573 do ip=1,gr%np
574 tmp1(ip) = abs( total_density(ip))**((m_one*xc_photons%pt%dim+m_two)/(xc_photons%pt%dim*m_one))
575 end do
576 !$OMP end parallel do
577
578 ! sum over the prefactors
579 presum = m_zero
580 do ia = 1, xc_photons%pt%nmodes
581 presum = presum + prefactor(ia)
582 end do
583 presum = presum * (m_one/(m_two*unit_volume))**(m_two/(xc_photons%pt%dim*m_one)) / (xc_photons%pt%dim*m_one+m_two)
584 presum = presum * xc_photons%eta_c * xc_photons%pxlda_kappa
585
586 xc_photons%ex = xc_photons%eta_c * xc_photons%pxlda_kappa * presum * dmf_integrate(gr, tmp1)
587
588 end select
589
590 safe_deallocate_a(prefactor)
591 safe_deallocate_a(rho_aux)
592 safe_deallocate_a(grad_rho_aux)
593 safe_deallocate_a(px_source)
594 safe_deallocate_a(tmp1)
595 safe_deallocate_a(tmp2)
596 safe_deallocate_a(tmp3)
597 safe_deallocate_a(grad_vpx)
598 safe_deallocate_a(epx_force_module)
599
600 pop_sub(photon_free_vpx_lda)
601
602 contains
603
604 subroutine get_px_source(px_source)
605 real(real64), contiguous, intent(out) :: px_source(:)
606
607 call dderivatives_grad(gr%der, rho_aux(1:gr%np_part), grad_rho_aux)
608
609 px_source = m_zero
610 do ia = 1, xc_photons%pt%nmodes
611 call lalg_gemv(gr%np, xc_photons%pt%dim, m_one, grad_rho_aux, xc_photons%pt%dressed_pol(:, ia), m_zero, tmp1)
612 call dderivatives_grad(gr%der, tmp1, tmp2)
613 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)
614 call lalg_axpy(gr%np, prefactor(ia), tmp1, px_source)
615 end do
616
617 end subroutine get_px_source
618
619 end subroutine photon_free_vpx_lda
620 ! ---------------------------------------------------------
621
622 ! ---------------------------------------------------------
624 subroutine laplacian_op(x, hx, userdata)
625 real(real64), contiguous, intent(in) :: x(:)
626 real(real64), contiguous, intent(out) :: Hx(:)
627 type(c_ptr), intent(in) :: userdata(:)
628
629 real(real64), allocatable :: tmpx(:)
630 type(grid_t), pointer :: gr
631
632 assert(size(userdata) == 1)
633 assert(c_associated(userdata(1)))
634 call c_f_pointer(userdata(1), gr)
635
636 safe_allocate(tmpx(1:gr%np_part))
637 call lalg_copy(gr%np, x, tmpx)
638 call dderivatives_lapl(gr%der, tmpx, hx, factor=-m_one)
639 safe_deallocate_a(tmpx)
640
641 end subroutine laplacian_op
642 ! ---------------------------------------------------------
643
644 ! ---------------------------------------------------------
645 !
671 !
672 ! The Lamb shift code is experimental and untested
673
674 subroutine photon_free_vpx_wfc(namespace, xc_photons, total_density, gr, space, st)
675 type(namespace_t), intent(in) :: namespace
676 type(xc_photons_t), intent(inout) :: xc_photons
677 real(real64), pointer, contiguous, intent(in) :: total_density(:)
678 class(grid_t), intent(in) :: gr
679 type(space_t), intent(in) :: space
680 type(states_elec_t), intent(in) :: st
681
682 integer :: ia, ip
683 real(real64) :: prefactor_lamb
684 real(real64) :: xx(space%dim), r
685 real(real64), allocatable :: prefactor(:)
686 real(real64), allocatable :: rho_aux(:)
687 real(real64), allocatable :: grad_rho_aux(:,:)
688 real(real64), allocatable :: grad_vpx(:,:)
689 real(real64), allocatable :: epsgrad_epsgrad_rho_aux(:)
690 real(real64), allocatable :: tmp1(:)
691 real(real64), allocatable :: tmp2(:,:)
692 real(real64) :: shift
693
694 push_sub(photon_free_vpx_wfc)
695
696 if (st%d%nspin >1) then
697 call messages_not_implemented("PhotonXCXCMethod = wavefunction for polarized and spinor cases", namespace=namespace)
698 end if
700 if (xc_photons%pt%n_electrons >1) then
701 call messages_not_implemented("PhotonXCXCMethod = wavefunction for more than 1 electron", namespace=namespace)
702 end if
703
704 xc_photons%vpx = m_zero
705 xc_photons%ex = m_zero
706
707 safe_allocate(prefactor(1:xc_photons%pt%nmodes))
708 prefactor = m_zero
709 safe_allocate(rho_aux(1:gr%np_part))
710 rho_aux = m_zero
711 safe_allocate(grad_rho_aux(1:gr%np, 1:xc_photons%pt%dim))
712 grad_rho_aux = m_zero
713 safe_allocate(epsgrad_epsgrad_rho_aux(1:gr%np))
714 safe_allocate(grad_vpx(1:gr%np, 1:xc_photons%pt%dim))
715 safe_allocate(tmp1(1:gr%np_part))
716 safe_allocate(tmp2(1:gr%np_part, 1:xc_photons%pt%dim))
717
718
719 rho_aux(1:gr%np) = sqrt(abs( total_density(1:gr%np)))
720 !$OMP parallel do
721 do ip = 1, gr%np
722 rho_aux(ip) = safe_tol(rho_aux(ip),1e-18_real64)
723 end do
724 !$OMP end parallel do
725
726 if (xc_photons%llamb_freespace) then
727
728 ! experimental Lamb shift mode
729
730 prefactor_lamb = ( m_two/(m_three*m_pi) ) * xc_photons%lamb_omega / p_c**3
731 call dderivatives_lapl(gr%der, rho_aux(1:gr%np_part), epsgrad_epsgrad_rho_aux)
732 call lalg_scal(gr%np, prefactor_lamb, epsgrad_epsgrad_rho_aux)
733
734 else
735
736 do ia = 1, xc_photons%pt%nmodes
737 prefactor(ia) = (xc_photons%pt%dressed_lambda(ia)**2) / (m_two*xc_photons%pt%dressed_omega(ia)**2)
738 end do
739
740 call dderivatives_grad(gr%der, rho_aux, grad_rho_aux)
741
742 epsgrad_epsgrad_rho_aux = m_zero
743 do ia = 1, xc_photons%pt%nmodes
744 tmp1 = m_zero
745 call lalg_gemv(gr%np, xc_photons%pt%dim, m_one, grad_rho_aux, xc_photons%pt%dressed_pol(:, ia), m_zero, tmp1)
746 call dderivatives_grad(gr%der, tmp1, tmp2)
747 call lalg_gemv(gr%np, xc_photons%pt%dim, prefactor(ia), tmp2, xc_photons%pt%dressed_pol(:, ia), &
748 m_one, epsgrad_epsgrad_rho_aux)
749 end do
750
751 end if
752
753 !$OMP parallel do
754 do ip = 1, gr%np
755 xc_photons%vpx(ip) = epsgrad_epsgrad_rho_aux(ip)/rho_aux(ip)
756 end do
757 !$OMP end parallel do
758
759 if(st%eigenval(1,1) < m_huge .and. .not. space%is_periodic()) then
760 shift = m_two * st%eigenval(1,1) * prefactor(1)
761 !$OMP parallel do
762 do ip=1, gr%np
763 xc_photons%vpx(ip) = xc_photons%vpx(ip) + shift
764 end do
765 !$OMP end parallel do
766 end if
767
768 call lalg_scal(gr%np, xc_photons%eta_c, xc_photons%vpx(:))
770 select case (xc_photons%energy_method)
771 case(1) ! virial
772
773 call dderivatives_grad(gr%der, xc_photons%vpx(:), grad_vpx)
774 !$OMP parallel do private(r, xx)
775 do ip = 1, gr%np
776 call mesh_r(gr, ip, r, coords=xx)
777 tmp1(ip) = - total_density(ip)*dot_product(xx(1:xc_photons%pt%dim), grad_vpx(ip,1:xc_photons%pt%dim))
778 end do
779 !$OMP end parallel do
780
781 xc_photons%ex = m_half*dmf_integrate(gr, tmp1)
782
783 case(2) ! expectation_value
784 xc_photons%ex = xc_photons%eta_c * dmf_dotp(gr, rho_aux(1:gr%np), epsgrad_epsgrad_rho_aux)
785
786 case default
787 assert(.false.)
788 end select
789
790 safe_deallocate_a(prefactor)
791 safe_deallocate_a(rho_aux)
792 safe_deallocate_a(grad_rho_aux)
793 safe_deallocate_a(grad_vpx)
794 safe_deallocate_a(epsgrad_epsgrad_rho_aux)
795 safe_deallocate_a(tmp1)
796 safe_deallocate_a(tmp2)
797
798 pop_sub(photon_free_vpx_wfc)
799
800 end subroutine photon_free_vpx_wfc
801 ! ---------------------------------------------------------
802
803
804 ! ---------------------------------------------------------
814 !
815 subroutine xc_photons_add_mean_field(xc_photons, gr, space, kpoints, st, time, dt)
816 class(xc_photons_t), intent(inout) :: xc_photons
817 class(grid_t), intent(in) :: gr
818 class(space_t), intent(in) :: space
819 type(kpoints_t), intent(in) :: kpoints
820 type(states_elec_t), intent(inout) :: st
821 real(real64), intent(in) :: time
822 real(real64), intent(in) :: dt
823
824
825 integer :: ia, idir, ispin
826 real(real64) :: pol_dot_jp
827 real(real64), allocatable :: current(:,:,:)
828 real(real64), allocatable :: jp(:)
829
831
832 ! compute the dressed photon-free vector potential
833 xc_photons%mf_vector_potential = m_zero
834
835 safe_allocate(current(1:gr%np_part, 1:space%dim, 1:st%d%nspin))
836 current = m_zero
837 ! here we use the paramagnetic current; note that the physical current here
838 ! only contains the paramagnetic current.
839 call states_elec_calc_quantities(gr, st, kpoints, .false., paramagnetic_current = current)
840
841 ! compute the paramagnetic current
842 safe_allocate(jp(1:space%dim))
843 do idir = 1, space%dim
844 jp(idir) = m_zero
845 do ispin = 1, st%d%spin_channels
846 jp(idir) = jp(idir) + dmf_integrate(gr%der%mesh, current(:,idir,ispin))
847 end do
848 end do
849
850 ! update the 'projected' current
851 do ia = 1, xc_photons%pt%nmodes
852 pol_dot_jp = dot_product(xc_photons%pt%dressed_pol(1:space%dim, ia),jp(1:space%dim))
853 xc_photons%jp_proj_eo(ia,1) = xc_photons%jp_proj_eo(ia,1) + &
854 cos(xc_photons%pt%dressed_omega(ia)*( time-dt))*pol_dot_jp*dt
855 xc_photons%jp_proj_eo(ia,2) = xc_photons%jp_proj_eo(ia,2) + &
856 sin(xc_photons%pt%dressed_omega(ia)*( time-dt))*pol_dot_jp*dt
857 end do
858
859 do ia = 1, xc_photons%pt%nmodes
860 xc_photons%mf_vector_potential(1:xc_photons%pt%dim) = xc_photons%mf_vector_potential(1:xc_photons%pt%dim) &
861 + (-p_c*(xc_photons%pt%dressed_lambda(ia)**2) / xc_photons%pt%dressed_omega(ia)) &
862 * (xc_photons%jp_proj_eo(ia,1)*sin(xc_photons%pt%dressed_omega(ia)* time) &
863 - xc_photons%jp_proj_eo(ia,2)*cos(xc_photons%pt%dressed_omega(ia)* time)) &
864 * xc_photons%pt%dressed_pol(1:xc_photons%pt%dim, ia)
865 end do
866
867 safe_deallocate_a(current)
868 safe_deallocate_a(jp)
869
871
872 end subroutine xc_photons_add_mean_field
873 ! ---------------------------------------------------------
874
875
879 !
880 logical pure function xc_photons_wants_to_renormalize_mass(xc_photons) result (renorm)
881 class(xc_photons_t), intent(in) :: xc_photons
882
883 renorm = (xc_photons%method>0) .and. xc_photons%llamb_freespace .and. xc_photons%llamb_re_mass
884
886
888 real(real64) pure function xc_photons_get_renormalized_emass(xc_photons) result(mass)
889 class(xc_photons_t), intent(in) :: xc_photons
890
891 mass = m_one - (m_four*xc_photons%lamb_omega) / (3.0*m_pi * p_c**3)
892
894
896 !
897 subroutine xc_photons_mf_dump(xc_photons, restart, ierr)
898 class(xc_photons_t), intent(in) :: xc_photons
899 type(restart_t), intent(in) :: restart
900 integer, intent(out) :: ierr
901
902 character(len=80), allocatable :: lines(:)
903 integer :: iunit, err, jj, nmodes
904 integer :: pt_dim
905
906 push_sub(photon_free_mf_dump)
907 nmodes = xc_photons%pt%nmodes
908 pt_dim = xc_photons%pt%dim
909
910 safe_allocate(lines(1:nmodes+pt_dim))
911
912 ierr = 0
913
914 iunit = restart%open('photon_free_mf')
915
916 do jj = 1, pt_dim
917 write(lines(jj), '(2x, es19.12)') xc_photons%mf_vector_potential(jj)
918 end do
919
920 do jj = 1, nmodes
921 write(lines(jj+pt_dim), '(a10,1x,I8,a1,2x,2(es19.12,2x))') 'Mode ', jj, ":", xc_photons%jp_proj_eo(jj,1:2)
922 end do
923
924 call restart%write(iunit, lines, nmodes+pt_dim, err)
925 if (err /= 0) ierr = ierr + 1
926 call restart%close(iunit)
927
928 safe_deallocate_a(lines)
929
930 pop_sub(photon_free_mf_dump)
931 end subroutine xc_photons_mf_dump
932
933! ---------------------------------------------------------
934
936 !
937 subroutine xc_photons_mf_load(xc_photons, restart, space, ierr)
938 class(xc_photons_t), intent(inout) :: xc_photons
939 type(restart_t), intent(in) :: restart
940 class(space_t), intent(in) :: space
941 integer, intent(out) :: ierr
942
943 character(len=80), allocatable :: lines(:)
944 character(len=7) :: sdummy
945 integer :: idummy
946 integer :: iunit, err, jj, nmodes
947 integer :: pt_dim
948
949 push_sub(photon_free_mf_load)
950
951 ierr = 0
952 nmodes = xc_photons%pt%nmodes
953 pt_dim = xc_photons%pt%dim
954
955 if (restart%skip()) then
956 ierr = -1
957 pop_sub(photon_free_mf_load)
958 return
959 end if
960
961 message(1) = "Debug: Reading Photon-Free Photons restart."
962 call messages_info(1, namespace=restart%namespace, debug_only=.true.)
963
964 if ( .not. allocated(xc_photons%jp_proj_eo)) then
965 safe_allocate(xc_photons%jp_proj_eo(1:xc_photons%pt%nmodes, 1:2))
966 xc_photons%jp_proj_eo = m_zero
967 end if
968 if ( .not. allocated(xc_photons%mf_vector_potential)) then
969 safe_allocate(xc_photons%mf_vector_potential(1:space%dim))
970 xc_photons%mf_vector_potential = m_zero
971 end if
972
973 safe_allocate(lines(1:nmodes+pt_dim))
974 iunit = restart%open('photon_free_mf')
975 call restart%read(iunit, lines, nmodes+pt_dim, err)
976 if (err /= 0) then
977 ierr = ierr + 1
978 else
979 do jj = 1, pt_dim
980 read(lines(jj),'(2x, es19.12)') xc_photons%mf_vector_potential(jj)
981 end do
982
983 do jj = 1, nmodes
984 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)
985 end do
986 end if
987 call restart%close(iunit)
988
989 message(1) = "Debug: Reading Photons restart done."
990 call messages_info(1, namespace=restart%namespace, debug_only=.true.)
991
992 safe_deallocate_a(lines)
993
994 pop_sub(photon_free_mf_load)
995 end subroutine xc_photons_mf_load
996
997end module xc_photons_oct_m
998
999!! Local Variables:
1000!! mode: f90
1001!! coding: utf-8
1002!! 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:193
real(real64), parameter, public m_huge
Definition: global.F90:209
real(real64), parameter, public m_zero
Definition: global.F90:191
real(real64), parameter, public m_four
Definition: global.F90:195
real(real64), parameter, public m_third
Definition: global.F90:198
real(real64), parameter, public m_pi
some mathematical constants
Definition: global.F90:189
real(real64), parameter, public m_half
Definition: global.F90:197
real(real64), parameter, public p_c
Electron gyromagnetic ratio, see Phys. Rev. Lett. 130, 071801 (2023)
Definition: global.F90:229
real(real64), parameter, public m_one
Definition: global.F90:192
real(real64), parameter, public m_three
Definition: global.F90:194
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:341
subroutine, public messages_not_implemented(feature, namespace)
Definition: messages.F90:1091
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:410
subroutine, public messages_experimental(name, namespace)
Definition: messages.F90:1063
subroutine, public messages_info(no_lines, iunit, debug_only, stress, all_nodes, namespace)
Definition: messages.F90:594
logical function, public parse_is_defined(namespace, name)
Definition: parser.F90:455
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:720
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:770
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:976
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:911
subroutine xc_photons_mf_dump(xc_photons, restart, ierr)
write restart information
Definition: xc_photons.F90:993
subroutine xc_photons_v_ks(xc_photons, namespace, total_density, gr, space, psolver, 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:984
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:474
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:700