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, st)
294 real(real64),
pointer,
contiguous,
intent(in) :: total_density(:)
295 class(
grid_t),
intent(in) :: gr
296 type(
space_t),
intent(in) :: space
304 xc_photons%lpfmf = xc_photons%method > 0
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
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
319 select case(xc_photons%method)
324 case(xc_photons_none)
332 if (.not. xc_photons%llamb_freespace)
then
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))
381 real(real64),
pointer,
contiguous,
intent(in) :: total_density(:)
382 type(
grid_t),
target,
intent(in) :: gr
383 type(
space_t),
intent(in) :: space
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(:)
400 real(real64),
parameter :: threshold = 1e-7_real64
404 if (xc_photons%energy_method == 2 .and. xc_photons%pt%n_electrons >1)
then
411 safe_allocate(prefactor(1:xc_photons%pt%nmodes))
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))
422 safe_allocate(epx_force_module(1:gr%np_part))
425 select case(xc_photons%pt%dim)
438 rho_aux(ip) = ( abs(total_density(ip))/(
m_two*unit_volume) )**(
m_two/(xc_photons%pt%dim*
m_one))
445 if (xc_photons%llamb_freespace)
then
450 prefactor_lamb = -(8.0_real64*
m_pi*
m_third) * xc_photons%lamb_omega / (
p_c**3)
454 xc_photons%vpx(ip) = prefactor_lamb*rho_aux(ip)
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
464 select case(xc_photons%pt%dim)
468 do ia = 1, xc_photons%pt%nmodes
471 px_source(ip) = px_source(ip) + prefactor(ia)
478 xc_photons%vpx(ip) = px_source(ip)*rho_aux(ip)
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
504 call dpoisson_solve(psolver, namespace, xc_photons%vpx(:), px_source)
512 call lalg_scal(gr%np, (xc_photons%eta_c * xc_photons%pxlda_kappa), xc_photons%vpx)
516 select case (xc_photons%energy_method)
519 do ia = 1, xc_photons%pt%nmodes
524 epx_force_module(ip) = -prefactor(ia)*
m_two*abs(total_density(ip))*rho_aux(ip)/(xc_photons%pt%dim*
m_one+
m_two)
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)
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))
541 xc_photons%ex = xc_photons%eta_c * xc_photons%pxlda_kappa * xc_photons%ex
545 rho_aux(1:gr%np) =
sqrt(abs( total_density(1:gr%np)))
547 safe_allocate(epsgrad_epsgrad_rho_aux(1:gr%np))
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)
554 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 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)
560 epsgrad_epsgrad_rho_aux(ip) = epsgrad_epsgrad_rho_aux(ip) + prefact*tmp1(ip)
566 xc_photons%ex = xc_photons%eta_c *
dmf_dotp(gr, rho_aux(1:gr%np), epsgrad_epsgrad_rho_aux(1:gr%np))
568 safe_deallocate_a(epsgrad_epsgrad_rho_aux)
574 tmp1(ip) = abs( total_density(ip))**((
m_one*xc_photons%pt%dim+
m_two)/(xc_photons%pt%dim*
m_one))
580 do ia = 1, xc_photons%pt%nmodes
581 presum = presum + prefactor(ia)
584 presum = presum * xc_photons%eta_c * xc_photons%pxlda_kappa
586 xc_photons%ex = xc_photons%eta_c * xc_photons%pxlda_kappa * presum *
dmf_integrate(gr, tmp1)
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)
605 real(real64),
contiguous,
intent(out) :: px_source(:)
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)
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)
625 real(real64),
contiguous,
intent(in) :: x(:)
626 real(real64),
contiguous,
intent(out) :: Hx(:)
627 type(c_ptr),
intent(in) :: userdata(:)
629 real(real64),
allocatable :: tmpx(:)
630 type(
grid_t),
pointer :: gr
632 assert(
size(userdata) == 1)
633 assert(c_associated(userdata(1)))
634 call c_f_pointer(userdata(1), gr)
636 safe_allocate(tmpx(1:gr%np_part))
639 safe_deallocate_a(tmpx)
677 real(real64),
pointer,
contiguous,
intent(in) :: total_density(:)
678 class(
grid_t),
intent(in) :: gr
679 type(
space_t),
intent(in) :: space
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
696 if (st%d%nspin >1)
then
697 call messages_not_implemented(
"PhotonXCXCMethod = wavefunction for polarized and spinor cases", namespace=namespace)
700 if (xc_photons%pt%n_electrons >1)
then
707 safe_allocate(prefactor(1:xc_photons%pt%nmodes))
709 safe_allocate(rho_aux(1:gr%np_part))
711 safe_allocate(grad_rho_aux(1:gr%np, 1:xc_photons%pt%dim))
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))
719 rho_aux(1:gr%np) =
sqrt(abs( total_density(1:gr%np)))
722 rho_aux(ip) = safe_tol(rho_aux(ip),1e-18_real64)
726 if (xc_photons%llamb_freespace)
then
732 call lalg_scal(gr%np, prefactor_lamb, epsgrad_epsgrad_rho_aux)
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)
742 epsgrad_epsgrad_rho_aux =
m_zero
743 do ia = 1, xc_photons%pt%nmodes
745 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 lalg_gemv(gr%np, xc_photons%pt%dim, prefactor(ia), tmp2, xc_photons%pt%dressed_pol(:, ia), &
748 m_one, epsgrad_epsgrad_rho_aux)
755 xc_photons%vpx(ip) = epsgrad_epsgrad_rho_aux(ip)/rho_aux(ip)
759 if(st%eigenval(1,1) <
m_huge .and. .not. space%is_periodic())
then
760 shift =
m_two * st%eigenval(1,1) * prefactor(1)
763 xc_photons%vpx(ip) = xc_photons%vpx(ip) + shift
768 call lalg_scal(gr%np, xc_photons%eta_c, xc_photons%vpx(:))
770 select case (xc_photons%energy_method)
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))
784 xc_photons%ex = xc_photons%eta_c *
dmf_dotp(gr, rho_aux(1:gr%np), epsgrad_epsgrad_rho_aux)
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)
817 class(
grid_t),
intent(in) :: gr
818 class(
space_t),
intent(in) :: space
821 real(real64),
intent(in) :: time
822 real(real64),
intent(in) :: dt
825 integer :: ia, idir, ispin
826 real(real64) :: pol_dot_jp
827 real(real64),
allocatable :: current(:,:,:)
828 real(real64),
allocatable :: jp(:)
833 xc_photons%mf_vector_potential =
m_zero
835 safe_allocate(current(1:gr%np_part, 1:space%dim, 1:st%d%nspin))
842 safe_allocate(jp(1:space%dim))
843 do idir = 1, space%dim
845 do ispin = 1, st%d%spin_channels
846 jp(idir) = jp(idir) +
dmf_integrate(gr%der%mesh, current(:,idir,ispin))
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
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)
867 safe_deallocate_a(current)
868 safe_deallocate_a(jp)
880 logical pure function xc_photons_wants_to_renormalize_mass(xc_photons) result (renorm)
883 renorm = (xc_photons%method>0) .and. xc_photons%llamb_freespace .and. xc_photons%llamb_re_mass
888 real(real64)
pure function xc_photons_get_renormalized_emass(xc_photons) result(mass)
891 mass = m_one - (m_four*xc_photons%lamb_omega) / (3.0*m_pi * p_c**3)
899 type(restart_t),
intent(in) :: restart
900 integer,
intent(out) :: ierr
902 character(len=80),
allocatable :: lines(:)
903 integer :: iunit, err, jj, nmodes
906 push_sub(photon_free_mf_dump)
907 nmodes = xc_photons%pt%nmodes
908 pt_dim = xc_photons%pt%dim
910 safe_allocate(lines(1:nmodes+pt_dim))
914 iunit = restart%open(
'photon_free_mf')
917 write(lines(jj),
'(2x, es19.12)') xc_photons%mf_vector_potential(jj)
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)
924 call restart%write(iunit, lines, nmodes+pt_dim, err)
925 if (err /= 0) ierr = ierr + 1
926 call restart%close(iunit)
928 safe_deallocate_a(lines)
930 pop_sub(photon_free_mf_dump)
939 type(restart_t),
intent(in) :: restart
940 class(space_t),
intent(in) :: space
941 integer,
intent(out) :: ierr
943 character(len=80),
allocatable :: lines(:)
944 character(len=7) :: sdummy
946 integer :: iunit, err, jj, nmodes
949 push_sub(photon_free_mf_load)
952 nmodes = xc_photons%pt%nmodes
953 pt_dim = xc_photons%pt%dim
955 if (restart%skip())
then
957 pop_sub(photon_free_mf_load)
961 message(1) =
"Debug: Reading Photon-Free Photons restart."
962 call messages_info(1, namespace=restart%namespace, debug_only=.
true.)
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
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
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)
980 read(lines(jj),
'(2x, es19.12)') xc_photons%mf_vector_potential(jj)
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)
987 call restart%close(iunit)
989 message(1) =
"Debug: Reading Photons restart done."
990 call messages_info(1, namespace=restart%namespace, debug_only=.
true.)
992 safe_deallocate_a(lines)
994 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, 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)