67 real(real64),
allocatable,
public :: vpx(:)
68 real(real64),
public :: ex
69 type(photon_mode_t) :: pt
70 real(real64) :: pxlda_kappa
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
79 logical,
public :: lpfmf = .false.
80 real(real64),
allocatable,
public :: mf_vector_potential(:)
81 real(real64),
allocatable :: jp_proj_eo(:,:)
97 integer,
private,
parameter :: &
98 XC_PHOTONS_NONE = 0, &
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
117 xc_photons%lpfmf = .false.
124 select case(xc_photon)
125 case(option__xcphotonfunctional__photon_x_lda)
127 xc_photons%lcorrelations = .false.
128 case(option__xcphotonfunctional__photon_xc_lda)
130 xc_photons%lcorrelations = .
true.
131 case(option__xcphotonfunctional__photon_x_wfn)
133 xc_photons%lcorrelations = .false.
134 case(option__xcphotonfunctional__photon_xc_wfn)
136 xc_photons%lcorrelations = .
true.
138 xc_photons%method = xc_photons_none
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"
195 if (xc_photons%lcorrelations)
then
208 message(1) =
"Defining PhotonXCEtaC is required for photon functionals containing correlation."
214 xc_photons%eta_c =
m_one
219 safe_allocate(xc_photons%vpx(1:gr%np_part))
231 call parse_variable(namespace,
'PhotonXCLambShift', .false., xc_photons%llamb_freespace)
234 if (xc_photons%llamb_freespace)
then
254 call parse_variable(namespace,
'PhotonXCLambShiftRenormalizeMass', .false., xc_photons%llamb_re_mass)
273 safe_deallocate_a(this%vpx)
275 if (
allocated(this%mf_vector_potential))
then
276 safe_deallocate_a(this%mf_vector_potential)
278 if (
allocated(this%jp_proj_eo))
then
279 safe_deallocate_a(this%jp_proj_eo)
291 subroutine xc_photons_v_ks(xc_photons, namespace, total_density, gr, space, psolver, ep, st)
294 real(real64),
pointer,
contiguous,
intent(in) :: total_density(:)
295 class(
grid_t),
intent(in) :: gr
296 type(
space_t),
intent(in) :: space
298 type(
epot_t),
intent(in) :: ep
305 xc_photons%lpfmf = xc_photons%method > 0
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
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
320 select case(xc_photons%method)
325 case(xc_photons_none)
333 if (.not. xc_photons%llamb_freespace)
then
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))
382 real(real64),
pointer,
contiguous,
intent(in) :: total_density(:)
383 type(
grid_t),
target,
intent(in) :: gr
384 type(
space_t),
intent(in) :: space
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(:)
401 real(real64),
parameter :: threshold = 1e-7_real64
405 if (xc_photons%energy_method == 2 .and. xc_photons%pt%n_electrons >1)
then
412 safe_allocate(prefactor(1:xc_photons%pt%nmodes))
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))
423 safe_allocate(epx_force_module(1:gr%np_part))
426 select case(xc_photons%pt%dim)
439 rho_aux(ip) = ( abs(total_density(ip))/(
m_two*unit_volume) )**(
m_two/(xc_photons%pt%dim*
m_one))
446 if (xc_photons%llamb_freespace)
then
451 prefactor_lamb = -(8.0_real64*
m_pi*
m_third) * xc_photons%lamb_omega / (
p_c**3)
455 xc_photons%vpx(ip) = prefactor_lamb*rho_aux(ip)
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
465 select case(xc_photons%pt%dim)
469 do ia = 1, xc_photons%pt%nmodes
472 px_source(ip) = px_source(ip) + prefactor(ia)
479 xc_photons%vpx(ip) = px_source(ip)*rho_aux(ip)
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
505 call dpoisson_solve(psolver, namespace, xc_photons%vpx(:), px_source)
513 call lalg_scal(gr%np, (xc_photons%eta_c * xc_photons%pxlda_kappa), xc_photons%vpx)
517 select case (xc_photons%energy_method)
520 do ia = 1, xc_photons%pt%nmodes
525 epx_force_module(ip) = -prefactor(ia)*
m_two*abs(total_density(ip))*rho_aux(ip)/(xc_photons%pt%dim*
m_one+
m_two)
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)
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))
542 xc_photons%ex = xc_photons%eta_c * xc_photons%pxlda_kappa * xc_photons%ex
546 rho_aux(1:gr%np) =
sqrt(abs( total_density(1:gr%np)))
548 safe_allocate(epsgrad_epsgrad_rho_aux(1:gr%np))
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)
555 call lalg_gemv(gr%np, xc_photons%pt%dim,
m_one, grad_rho_aux, xc_photons%pt%dressed_pol(:, ia),
m_zero, tmp1)
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)
561 epsgrad_epsgrad_rho_aux(ip) = epsgrad_epsgrad_rho_aux(ip) + prefact*tmp1(ip)
567 xc_photons%ex = xc_photons%eta_c *
dmf_dotp(gr, rho_aux(1:gr%np), epsgrad_epsgrad_rho_aux(1:gr%np))
569 safe_deallocate_a(epsgrad_epsgrad_rho_aux)
575 tmp1(ip) = abs( total_density(ip))**((
m_one*xc_photons%pt%dim+
m_two)/(xc_photons%pt%dim*
m_one))
581 do ia = 1, xc_photons%pt%nmodes
582 presum = presum + prefactor(ia)
585 presum = presum * xc_photons%eta_c * xc_photons%pxlda_kappa
587 xc_photons%ex = xc_photons%eta_c * xc_photons%pxlda_kappa * presum *
dmf_integrate(gr, tmp1)
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)
606 real(real64),
contiguous,
intent(out) :: px_source(:)
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)
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)
626 real(real64),
contiguous,
intent(in) :: x(:)
627 real(real64),
contiguous,
intent(out) :: Hx(:)
628 type(c_ptr),
intent(in) :: userdata(:)
630 real(real64),
allocatable :: tmpx(:)
631 type(
grid_t),
pointer :: gr
633 assert(
size(userdata) == 1)
634 assert(c_associated(userdata(1)))
635 call c_f_pointer(userdata(1), gr)
637 safe_allocate(tmpx(1:gr%np_part))
640 safe_deallocate_a(tmpx)
678 real(real64),
pointer,
contiguous,
intent(in) :: total_density(:)
679 class(
grid_t),
intent(in) :: gr
680 type(
space_t),
intent(in) :: space
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
697 if (st%d%nspin >1)
then
698 call messages_not_implemented(
"PhotonXCXCMethod = wavefunction for polarized and spinor cases", namespace=namespace)
701 if (xc_photons%pt%n_electrons >1)
then
708 safe_allocate(prefactor(1:xc_photons%pt%nmodes))
710 safe_allocate(rho_aux(1:gr%np_part))
712 safe_allocate(grad_rho_aux(1:gr%np, 1:xc_photons%pt%dim))
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))
720 rho_aux(1:gr%np) =
sqrt(abs( total_density(1:gr%np)))
723 rho_aux(ip) = safe_tol(rho_aux(ip),1e-18_real64)
727 if (xc_photons%llamb_freespace)
then
733 call lalg_scal(gr%np, prefactor_lamb, epsgrad_epsgrad_rho_aux)
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)
743 epsgrad_epsgrad_rho_aux =
m_zero
744 do ia = 1, xc_photons%pt%nmodes
746 call lalg_gemv(gr%np, xc_photons%pt%dim,
m_one, grad_rho_aux, xc_photons%pt%dressed_pol(:, ia),
m_zero, tmp1)
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)
756 xc_photons%vpx(ip) = epsgrad_epsgrad_rho_aux(ip)/rho_aux(ip)
760 if(st%eigenval(1,1) <
m_huge .and. .not. space%is_periodic())
then
761 shift =
m_two * st%eigenval(1,1) * prefactor(1)
764 xc_photons%vpx(ip) = xc_photons%vpx(ip) + shift
769 call lalg_scal(gr%np, xc_photons%eta_c, xc_photons%vpx(:))
771 select case (xc_photons%energy_method)
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))
785 xc_photons%ex = xc_photons%eta_c *
dmf_dotp(gr, rho_aux(1:gr%np), epsgrad_epsgrad_rho_aux)
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)
818 class(
grid_t),
intent(in) :: gr
819 class(
space_t),
intent(in) :: space
822 real(real64),
intent(in) :: time
823 real(real64),
intent(in) :: dt
826 integer :: ia, idir, ispin
827 real(real64) :: pol_dot_jp
828 real(real64),
allocatable :: current(:,:,:)
829 real(real64),
allocatable :: jp(:)
834 xc_photons%mf_vector_potential =
m_zero
836 safe_allocate(current(1:gr%np_part, 1:space%dim, 1:st%d%nspin))
843 safe_allocate(jp(1:space%dim))
844 do idir = 1, space%dim
846 do ispin = 1, st%d%spin_channels
847 jp(idir) = jp(idir) +
dmf_integrate(gr%der%mesh, current(:,idir,ispin))
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
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)
868 safe_deallocate_a(current)
869 safe_deallocate_a(jp)
881 logical pure function xc_photons_wants_to_renormalize_mass(xc_photons) result (renorm)
884 renorm = (xc_photons%method>0) .and. xc_photons%llamb_freespace .and. xc_photons%llamb_re_mass
889 real(real64)
pure function xc_photons_get_renormalized_emass(xc_photons) result(mass)
892 mass = m_one - (m_four*xc_photons%lamb_omega) / (3.0*m_pi * p_c**3)
900 type(restart_t),
intent(in) :: restart
901 integer,
intent(out) :: ierr
903 character(len=80),
allocatable :: lines(:)
904 integer :: iunit, err, jj, nmodes
907 push_sub(photon_free_mf_dump)
908 nmodes = xc_photons%pt%nmodes
909 pt_dim = xc_photons%pt%dim
911 safe_allocate(lines(1:nmodes+pt_dim))
915 iunit = restart%open(
'photon_free_mf')
918 write(lines(jj),
'(2x, es19.12)') xc_photons%mf_vector_potential(jj)
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)
925 call restart%write(iunit, lines, nmodes+pt_dim, err)
926 if (err /= 0) ierr = ierr + 1
927 call restart%close(iunit)
929 safe_deallocate_a(lines)
931 pop_sub(photon_free_mf_dump)
940 type(restart_t),
intent(in) :: restart
941 class(space_t),
intent(in) :: space
942 integer,
intent(out) :: ierr
944 character(len=80),
allocatable :: lines(:)
945 character(len=7) :: sdummy
947 integer :: iunit, err, jj, nmodes
950 push_sub(photon_free_mf_load)
953 nmodes = xc_photons%pt%nmodes
954 pt_dim = xc_photons%pt%dim
956 if (restart%skip())
then
958 pop_sub(photon_free_mf_load)
962 message(1) =
"Debug: Reading Photon-Free Photons restart."
963 call messages_info(1, namespace=restart%namespace, debug_only=.
true.)
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
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
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)
981 read(lines(jj),
'(2x, es19.12)') xc_photons%mf_vector_potential(jj)
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)
988 call restart%close(iunit)
990 message(1) =
"Debug: Reading Photons restart done."
991 call messages_info(1, namespace=restart%namespace, debug_only=.
true.)
993 safe_deallocate_a(lines)
995 pop_sub(photon_free_mf_load)
constant times a vector plus a vector
Copies a vector x, to a vector y.
scales a vector by a constant
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
real(real64), parameter, public m_huge
real(real64), parameter, public m_zero
real(real64), parameter, public m_four
real(real64), parameter, public m_third
real(real64), parameter, public m_pi
some mathematical constants
real(real64), parameter, public m_half
real(real64), parameter, public p_c
Electron gyromagnetic ratio, see Phys. Rev. Lett. 130, 071801 (2023)
real(real64), parameter, public m_one
real(real64), parameter, public m_three
This module implements the underlying real-space grid.
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.
pure subroutine, public mesh_r(mesh, ip, rr, origin, coords)
return the distance to the origin for a given grid point
subroutine, public messages_not_implemented(feature, namespace)
character(len=256), dimension(max_lines), public message
to be output by fatal, warning
subroutine, public messages_fatal(no_lines, only_root_writes, namespace)
subroutine, public messages_experimental(name, namespace)
subroutine, public messages_info(no_lines, iunit, debug_only, stress, all_nodes, namespace)
logical function, public parse_is_defined(namespace, name)
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.
This module is intended to contain "only mathematical" functions and procedures.
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.
subroutine laplacian_op(x, hx, userdata)
Computes Hx = (\Laplacian) x.
subroutine photon_free_vpx_wfc(namespace, xc_photons, total_density, gr, space, st)
compute the electron-photon exchange potential based on wave functions
logical pure function xc_photons_wants_to_renormalize_mass(xc_photons)
indicate whether the photon-exchange requires a renormalized electron mass
subroutine xc_photons_init(xc_photons, namespace, xc_photon, space, gr, st)
initialize the photon-exchange functional
subroutine xc_photons_add_mean_field(xc_photons, gr, space, kpoints, st, time, dt)
accumulate the results of time integral the paramagnetic current.
subroutine xc_photons_mf_dump(xc_photons, restart, ierr)
write restart information
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
real(real64) pure function xc_photons_get_renormalized_emass(xc_photons)
return the renormalized electron mass for the electron-photon exhange
integer, parameter, private xc_photons_lda
subroutine photon_free_vpx_lda(namespace, xc_photons, total_density, gr, space, psolver)
compute the electron-photon exchange potential within the LDA
subroutine xc_photons_end(this)
integer, parameter, private xc_photons_wfs
subroutine xc_photons_mf_load(xc_photons, restart, space, ierr)
load restart information
Description of the grid, containing information on derivatives, stencil, and symmetries.
The states_elec_t class contains all electronic wave functions.
This class described the 'photon-exchange' electron-photon xc functionals, based on QEDFT.
subroutine get_px_source(px_source)