56 integer,
parameter,
public :: &
57 SIC_NONE = 1, & !< no self-interaction correction
65 integer,
public :: level = sic_none
66 real(real64),
public :: amaldi_factor
67 type(xc_oep_t),
public :: oep
75 subroutine xc_sic_init(sic, namespace, gr, st, mc, space)
76 type(xc_sic_t),
intent(out) :: sic
77 type(namespace_t),
intent(in) :: namespace
78 type(grid_t),
intent(inout) :: gr
79 type(states_elec_t),
intent(in) :: st
80 type(multicomm_t),
intent(in) :: mc
81 class(space_t),
intent(in) :: space
107 call parse_variable(namespace,
'SICCorrection', sic_none, sic%level)
111 sic%amaldi_factor =
m_one
113 sic%amaldi_factor = (st%qtot -
m_one)/st%qtot
122 if(st%nik > st%d%spin_channels)
then
127 if (
allocated(st%rho_core))
then
135 if (space%is_periodic() .and. sic%level /= sic_none)
then
145 type(xc_sic_t),
intent(inout) :: sic
147 if (sic%level == sic_none)
return
159 type(xc_sic_t),
intent(in) :: sic
160 integer,
optional,
intent(in) :: iunit
163 if (sic%level == sic_none)
return
189 subroutine xc_sic_calc_adsic(sic, namespace, space, gr, st, hm, xc, density, vxc, ex, ec)
192 class(
space_t),
intent(in) :: space
193 type(
grid_t),
intent(in) :: gr
196 type(
xc_t),
intent(inout) :: xc
197 real(real64),
contiguous,
intent(in) :: density(:,:)
198 real(real64),
contiguous,
intent(inout) :: vxc(:,:)
199 real(real64),
optional,
intent(inout) :: ex, ec
201 integer :: ispin, ist, ik, ip
202 real(real64),
allocatable :: vxc_sic(:,:), vh_sic(:), rho(:, :)
203 real(real64) :: ex_sic, ec_sic, qsp(2)
204 real(real64) :: dtot, dpol, vpol
211 if (st%d%ispin ==
spinors .and. .not.
in_family(hm%xc%family, [xc_family_lda, xc_family_gga]))
then
212 write(
message(1),
'(a)')
'ADSIC with non-collinear spin is currently only possible'
213 write(
message(2),
'(a)')
'with LDA and GGA functionals.'
223 if( .not.
allocated(st%frozen_rho))
then
224 select case (st%d%ispin)
228 ispin = st%d%get_spin_index(ik)
229 qsp(ispin) = qsp(ispin) + st%occ(ist, ik) * st%kweights(ik)
239 safe_allocate(vxc_sic(1:gr%np, 1:2))
240 safe_allocate(vh_sic(1:gr%np))
241 safe_allocate(rho(1:gr%np, 1:2))
243 select case (st%d%ispin)
245 do ispin = 1, st%d%nspin
251 rho(:, ispin) = density(:, ispin) / qsp(ispin)
252 if(
present(ex) .and.
present(ec))
then
256 call xc_get_vxc(gr, xc, st, hm%kpoints, hm%psolver, namespace, space, &
257 rho,
spin_polarized, hm%ions%latt%rcell_volume, vxc_sic, ex = ex_sic, ec = ec_sic)
258 ex = ex - ex_sic * qsp(ispin)
259 ec = ec - ec_sic * qsp(ispin)
262 call xc_get_vxc(gr, xc, st, hm%kpoints, hm%psolver, namespace, space, &
271 call dpoisson_solve(hm%psolver, namespace, vh_sic, rho(:, ispin), all_nodes=.false.)
285 assert(
in_family(hm%xc%family, [xc_family_lda, xc_family_gga]))
292 dtot = density(ip, 1) + density(ip, 2)
293 dpol =
sqrt((density(ip, 1) - density(ip, 2))**2 + &
294 m_four*(density(ip, 3)**2 + density(ip, 4)**2))
302 if (nup <= 1e-14_real64) cycle
306 if(
present(ex) .and.
present(ec))
then
309 call xc_get_vxc(gr, xc, st, hm%kpoints, hm%psolver, namespace, space, &
310 rho,
spin_polarized, hm%ions%latt%rcell_volume, vxc_sic, ex = ex_sic, ec = ec_sic)
311 ex = ex - ex_sic * nup
312 ec = ec - ec_sic * nup
314 call xc_get_vxc(gr, xc, st, hm%kpoints, hm%psolver, namespace, space, &
326 call dpoisson_solve(hm%psolver, namespace, vh_sic, rho(:, ispin), all_nodes=.false.)
334 dpol =
sqrt((density(ip, 1) - density(ip, 2))**2 + &
335 m_four*(density(ip, 3)**2 + density(ip, 4)**2))
336 vpol = (vxc_sic(ip, 1) - vxc_sic(ip, 2))*(density(ip, 1) - density(ip, 2))/(safe_tol(dpol,
xc_tiny))
338 vxc(ip, 1) = vxc(ip, 1) -
m_half*(vxc_sic(ip, 1) + vxc_sic(ip, 2) + vpol)
339 vxc(ip, 2) = vxc(ip, 2) -
m_half*(vxc_sic(ip, 1) + vxc_sic(ip, 2) - vpol)
340 vxc(ip, 3) = vxc(ip, 3) - (vxc_sic(ip, 1) - vxc_sic(ip, 2))*density(ip, 3)/(safe_tol(dpol,
xc_tiny))
341 vxc(ip, 4) = vxc(ip, 4) - (vxc_sic(ip, 1) - vxc_sic(ip, 2))*density(ip, 4)/(safe_tol(dpol,
xc_tiny))
348 safe_deallocate_a(vxc_sic)
349 safe_deallocate_a(vh_sic)
350 safe_deallocate_a(rho)
361 type(
xc_t),
intent(in) :: xc
363 type(
grid_t),
intent(in) :: gr
364 real(real64),
intent(in) :: rho(:,:)
365 real(real64),
contiguous,
intent(inout) :: fxc(:,:,:)
367 real(real64),
allocatable :: rho_averaged(:, :)
368 real(real64),
allocatable :: fxc_sic(:,:,:)
376 assert(.not.
allocated(st%frozen_rho))
378 if (
bitand(xc%kernel_family, xc_family_lda) == 0)
then
379 message(1) =
"fxc calculation with ADSIC not implemented beyond LDAs."
387 if (st%d%ispin ==
spinors)
then
392 safe_allocate(rho_averaged(1:gr%np, 1:2))
393 safe_allocate(fxc_sic(1:gr%np, 1:2, 1:2))
395 do ispin = 1, st%d%nspin
400 call lalg_copy(gr%np, rho(:, ispin), rho_averaged(:, ispin))
406 call lalg_axpy(gr%np, -
m_one/qtot, fxc_sic(:, ispin, ispin), fxc(:, ispin, ispin))
409 safe_deallocate_a(rho_averaged)
410 safe_deallocate_a(fxc_sic)
constant times a vector plus a vector
Copies a vector x, to a vector y.
scales a vector by a constant
double sqrt(double __x) __attribute__((__nothrow__
integer, parameter, public unpolarized
Parameters...
integer, parameter, public spinors
integer, parameter, public spin_polarized
real(real64), parameter, public m_zero
real(real64), parameter, public m_four
real(real64), parameter, public m_half
real(real64), parameter, public m_one
real(real64), parameter, public m_min_occ
Minimal occupation that is considered to be non-zero.
This module implements the underlying real-space grid.
This module defines various routines, operating on mesh functions.
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_input_error(namespace, var, details, row, column)
This module handles the communicators for the various parallelization strategies.
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 handles spin dimensions of the states and the k-point distribution.
real(real64), parameter, public xc_tiny
Arbitrary definition of tiny, for use in XC context.
subroutine, public xc_get_vxc(gr, xcs, st, kpoints, psolver, namespace, space, rho, ispin, rcell_volume, vxc, ex, ec, deltaxc, vtau, ex_density, ec_density, stress_xc, force_orbitalfree)
logical function, public xc_is_not_size_consistent(xcs, namespace)
Is one of the x or c functional is not size consistent.
subroutine, public xc_get_fxc(xcs, gr, namespace, rho, ispin, fxc, fxc_grad, fxc_grad_spin)
Returns the exchange-correlation kernel.
pure logical function, public in_family(family, xc_families)
subroutine, public xc_oep_end(oep)
subroutine, public xc_oep_init(oep, namespace, gr, st, mc, space, oep_type)
integer, parameter, public oep_type_sic
subroutine, public xc_sic_write_info(sic, iunit, namespace)
integer, parameter, public sic_adsic
Averaged density SIC.
subroutine, public xc_sic_init(sic, namespace, gr, st, mc, space)
initialize the SIC object
subroutine, public xc_sic_end(sic)
finalize the SIC and, if needed, the included OEP
integer, parameter, public sic_pz_oep
Perdew-Zunger SIC (OEP way)
integer, parameter, public sic_amaldi
Amaldi correction term.
subroutine, public xc_sic_calc_adsic(sic, namespace, space, gr, st, hm, xc, density, vxc, ex, ec)
Computes the ADSIC potential and energy.
subroutine, public xc_sic_add_fxc_adsic(namespace, xc, st, gr, rho, fxc)
Adds to fxc the ADSIC contribution.
Description of the grid, containing information on derivatives, stencil, and symmetries.
The states_elec_t class contains all electronic wave functions.
This class contains information about the self-interaction correction.