#include <init.hpp>
Inheritance diagram for q3::PawPlanewavePseudopotentialKohnShamSCFIterator:Static Protected Attributes | |
| static constexpr Real | occupancy_damping_ = 0.4 |
Additional Inherited Members | |
| Public Attributes inherited from q3::KohnShamSCFIterator | |
| KohnShamSCFIteratorConfig | config_ |
| std::string | exchange_correlation_functional_ |
| int | methfessel_paxton_order_ |
| Real | mixing_beta_ |
| int | mixing_history_size_ |
| std::string | mixing_mode_ |
| std::string | mixing_preconditioning_ |
| Real | mixing_q0_ |
| std::string | smearing_ |
| Real | smearing_width_ |
|
inlineexplicit |
|
overridedefault |
|
inlineoverrideprotectedvirtual |
Damp ρ_ij between iterations, the way QE mixes rhobec alongside the density.
Without this the on-site occupancies are recomputed from fresh ψ every step and fed straight back through D, a loop with gain 1 that the density mixer cannot damp — lowering the mixing beta does nothing to it, which is how the problem identifies itself. Whether that matters is a property of the dataset: Q3 converges Si and C undamped, but the PseudoDojo Cu PAW dataset has a second p channel carrying an augmentation charge of 25.2 electrons and D_ion of −71.8 Ry, giving a one-center D of ~217 Ry, and the loop oscillates between two states ±212 Ry apart until it is damped.
Linear damping rather than the full mixer: the fixed point is unchanged either way, since at self-consistency the new and old occupancies are equal, so this only sets how fast the loop is allowed to move. Feeding ρ_ij through the density mixer as QE does would converge stiff cases in fewer iterations, and is the upgrade if one proves too slow.
The damped array is what the augmentation charge on the grid sees too, which is deliberate: the density and the on-site terms have to be built from the same occupancies or the energy is not a functional of one state.
The per-spin arrays are damped by the same factor and in the same step. Damping them independently, or damping only the total, would let ρ↑ + ρ↓ drift away from the total the augmentation charge on the grid was built from, and the energy would stop being a functional of one state.
Reimplemented from q3::UltrasoftPlanewavePseudopotentialKohnShamSCFIterator.
Definition at line 218 of file init.hpp.
References q3::UltrasoftPlanewavePseudopotentialKohnShamSCFIterator::accumulate_projector_overlaps(), damp_towards_previous(), occupancy_residual_, output_projector_overlaps_by_spin_, previous_projector_overlaps_, previous_projector_overlaps_by_spin_, q3::UltrasoftPlanewavePseudopotentialKohnShamSCFIterator::projector_overlap_coefficients_, and q3::UltrasoftPlanewavePseudopotentialKohnShamSCFIterator::projector_overlap_coefficients_by_spin_.
|
inlineoverrideprotectedvirtual |
−Σ_a Σ_ij ρ_ij D^oc_ij, on top of the base's −∫ρ·V_Hxc.
The base's term corrects for the grid potential inside the band energy; this corrects for the on-site potential that entered it through D. They are the same correction applied to two different parts of the same sum, which is why this adds rather than replaces.
Three occupancy arrays are in play by the time the energy is assembled, and each term has to use the right one:
ρ^in built the D this iteration's Hamiltonian was diagonalized with ρ^out the projections of the ψ that came out of it ρ^mix the damped array the next iteration will start from, which is also what E¹ and the refreshed D^oc are now built from
The band energy contains Σ ρ^out · D^oc(ρ^in), so that is what is removed, and E¹ is evaluated at ρ^mix with the first-order term that carries it back to ρ^out:
E¹(ρ^out) ≈ E¹(ρ^mix) − Σ (ρ^mix − ρ^out) · D^oc(ρ^mix)
which is the one-centre half of the SCF correction and the same expansion compute_delta_escf_correction makes for the grid. With both, the reported energy is second order in the residual of either channel. Quantum ESPRESSO's delta_escf carries the same PAW term.
What this used to subtract was Σ ρ^mix · D^oc(ρ^mix), mixed against refreshed, while the band energy holds out against input. The difference, Σ ρ^out · (D^oc(ρ^mix) − D^oc(ρ^in)), is first order in the step the damping takes, and the damping takes a fixed fraction of the residual – so the energy carried an error first order in the on-site residual. On silicon PAW at Γ that kept the energy drifting by ~3e-9 Ry per iteration after the density residual had reached 1e-14, and left a run stopped at a 1e-12 threshold 6e-9 Ry from where it settles.
At self-consistency the three arrays coincide and this reduces to the old expression, so a converged energy does not move; what moves is how close an unconverged one already is.
Implements q3::KohnShamSCFIterator.
Definition at line 322 of file init.hpp.
References q3::KohnShamSCFIterator::compute_double_counting_correction(), one_center_band_energy_term_, and one_center_scf_correction_.
|
inlineoverrideprotectedvirtual |
Compute sub-components (e.g. band energy, entropy) needed for total(). Default no-op; override in representations that compute eagerly.
Reimplemented from q3::KohnShamSCFIterator.
Definition at line 362 of file init.hpp.
References q3::KohnShamSCFIterator::compute_energy_components(), q3::PlanewavePseudopotentialKohnShamSCFIterator::energies_, q3::paw_one_center::PawOneCenterTerms::energy(), q3::NormConservingPlanewaveKohnShamEnergies::one_center_energy, and one_center_terms_.
|
inlineoverrideprotectedvirtual |
A run given the converged occupations keeps the one-centre term built from them. The ones just accumulated come from this run's own k-points and occupations – a band path's, say – and rebuilding the term from them would hand the first-iteration re-diagonalization a different D than the pass it repeats.
D still holds the one-centre term built from ρ^in – the one the Hamiltonian just diagonalized used – so this is the moment to contract it with ρ^out.
Implements q3::KohnShamSCFIterator.
Definition at line 164 of file init.hpp.
References q3::UltrasoftPlanewavePseudopotentialKohnShamSCFIterator::construct_density_from_states(), q3::log::info(), one_center_band_energy_term_, output_projector_overlaps_by_spin_, and preloaded_one_center_occupations_.
|
inlineprotected |
One damping step, in place; seeds previous on the first call.
Returns ‖fresh − damped‖² for that array, or the unconverged sentinel on the seeding call, where there is no previous array to have moved away from and zero would be a false "already converged".
Definition at line 286 of file init.hpp.
References occupancy_damping_, and q3::Matrix< T >::raw().
Referenced by accumulate_projector_overlaps().
|
inlineoverrideprotectedvirtual |
The density residual alone cannot see the on-site channel.
ρ_ij is a second self-consistent quantity with its own feedback loop, damped here at a fixed rate rather than driven by the density mixer, and it converges more slowly than the density and with a period-2 oscillation. Measured on Cu PAW at 80 Ry, the density residual runs two to four orders of magnitude below ‖δρ_ij‖² near convergence, and the gap widens as the run proceeds: at dr2 = 1.5e-13 the occupancy residual is still 6.7e-10. ρ_ij reaches the grid only through the compensation charge, which is confined to the augmentation sphere and so sits at the high-|G| end the Hartree metric's 1/|G|² weight suppresses, while the one-center energy integrates it unweighted. Exiting on dr2 alone therefore stops while E¹ − Ẽ¹ is still moving, and because the one-center channel has no first-order SCF correction of its own (the delta_escf analogue that makes the grid part second-order in its error), that residual motion lands in the total energy at first order. On Cu at 80 Ry that left ~0.6 mRy of scatter in the total energy while the stress on the same runs moved only ~0.4 kbar, which is why the defect reads as an energy-only one: the stress is evaluated once from the final state rather than assembled across iterations, so it carries no analogous first-order term.
Same reasoning, same shape as the DFT+U occupation-matrix gate in PlanewavePseudopotentialKohnShamSCFIterator::is_iteration_converged(): the density can be self-consistent while a potential built from a still-relaxing occupancy is not.
occupancy_residual_ is refreshed in accumulate_projector_overlaps(), which perform_iteration() always runs before this is checked, so it is never stale by an iteration. It is computed from arrays the ultrasoft accumulation has already pool- and band-reduced, by pure arithmetic with a rank-invariant damping scalar, so it needs no reduction of its own – unlike density_residual_, whose per-rank FFTW plan noise does.
Reimplemented from q3::KohnShamSCFIterator.
Definition at line 274 of file init.hpp.
References q3::SCFIterator::convergence_threshold_, q3::SCFIterator::is_iteration_converged(), and occupancy_residual_.
|
inlineoverridevirtual |
Hand over ρ_ij from a checkpoint, so D's one-centre term is built from the converged occupations before the first Hamiltonian and stays so for the whole run. Call after construction and before run(), like load_effective_potential(). Nothing to do where D has no such term.
Reimplemented from q3::PlanewavePseudopotentialKohnShamSCFIterator.
Definition at line 76 of file init.hpp.
References preloaded_one_center_occupations_.
|
inlineoverrideprotectedvirtual |
Hook invoked once the initial atomic-superposition charge density is ready.
Default no-op. The ultrasoft subclass overrides it to initialise the smooth density and projector overlap storage.
Reimplemented from q3::PlanewaveKohnShamSCFIterator.
Definition at line 149 of file init.hpp.
References q3::paw_one_center::PawOneCenterTerms::build(), q3::KohnShamSCFIterator::crystal_structure(), q3::KohnShamSCFIterator::exchange_correlation(), q3::log::info(), q3::paw_one_center::PawOneCenterTerms::is_active(), q3::UltrasoftPlanewavePseudopotentialKohnShamSCFIterator::on_charge_density_initialized(), one_center_terms_, preloaded_one_center_occupations_, and q3::PlanewavePseudopotentialKohnShamSCFIterator::quantum_system_snapshot_.
|
inlineoverridevirtual |
The on-site occupations ρ_ij per spin channel, which the checkpoint stores: what D's one-centre term is built from. Empty wherever D has no such term; the PAW iterator overrides this.
Reimplemented from q3::PlanewavePseudopotentialKohnShamSCFIterator.
Definition at line 67 of file init.hpp.
References q3::UltrasoftPlanewavePseudopotentialKohnShamSCFIterator::projector_overlap_coefficients_by_spin_.
|
inlineoverrideprotectedvirtual |
The base wired its own D into the quantum system; the one-center term is added to that same matrix in place, so nothing downstream needs to know a PAW branch exists.
Reimplemented from q3::KohnShamSCFIterator.
Definition at line 304 of file init.hpp.
References q3::UltrasoftPlanewavePseudopotentialKohnShamSCFIterator::prepare_hamiltonian().
|
inlineoverrideprotectedvirtual |
Representation-specific tail of the convergence report (eigenvalues, magnetization, final energy breakdown). Generic banner lives in report_convergence().
Second rung of the failure ladder. Printed in QE's labels and order for readability, but see BranchTotals: QE's four values are on a different footing and are not directly differenceable against these.
Implements q3::KohnShamSCFIterator.
Definition at line 367 of file init.hpp.
References q3::paw_one_center::PawOneCenterTerms::branch_totals(), q3::log::info(), one_center_terms_, and q3::KohnShamSCFIterator::report_convergence_details().
|
inlineoverridevirtual |
Whether a non-self-consistent pass must be handed those occupations before it diagonalizes – true exactly where D carries a one-centre term.
Reimplemented from q3::PlanewavePseudopotentialKohnShamSCFIterator.
|
staticconstexprprotected |
How far ρ_ij is allowed to move per iteration. Chosen to be loose enough that Si and C converge in the same iteration count they did undamped, and tight enough to hold the stiff Cu channel.
Definition at line 147 of file init.hpp.
Referenced by damp_towards_previous().
|
protected |
‖ρ_ij^fresh − ρ_ij^damped‖² — how far the occupancies just accumulated from ψ are from the ones that built D for this iteration, before the damping factor scales the step actually taken. Same squared-magnitude convention density_residual_ uses, and the same role hubbard_occupation_residual_ plays for DFT+U: it is what tells the SCF that the on-site channel is still moving.
Large sentinel until the first real damping step, so it can never make a run look converged before ρ_ij has been compared against anything.
Definition at line 140 of file init.hpp.
Referenced by accumulate_projector_overlaps(), and is_iteration_converged().
|
protected |
The point-group average of ρ_ij, built once from this cell's operations.
Inactive — and symmetrize() a no-op — whenever the density symmetrizer is, which is the point: control.disable_symmetry, a Γ-only basis and a cell with no operation beyond the identity all stand this down as well, so the on-site channel and the grid never disagree about whether the run is imposing symmetry.
|
protected |
Σ ρ_ij^out · D^oc(ρ^in): the part of the band energy that D's one-centre term put there. Band energy is a sum over ψ_out, whose projections are ρ^out, against the Hamiltonian this iteration was diagonalized with, whose one-centre D was built from ρ^in – so this is what the double-counting term has to remove.
Definition at line 118 of file init.hpp.
Referenced by compute_double_counting_correction(), and construct_density_from_states().
|
protected |
−Σ (ρ_ij^mix − ρ_ij^out) · D^oc(ρ^mix): the one-centre half of the SCF correction, the term that makes the reported energy second order in the on-site residual. The grid density has always had its half (compute_delta_escf_correction); see compute_double_counting_correction for why the one-centre channel needs one too.
Definition at line 127 of file init.hpp.
Referenced by compute_double_counting_correction().
|
protected |
Definition at line 81 of file init.hpp.
Referenced by compute_energy_components(), on_charge_density_initialized(), and report_convergence_details().
|
protected |
ρ_ij^out per spin: the occupancies this iteration's ψ produced, symmetrized, before damping replaces them with the next input. The energy needs all three arrays – in, out and mixed – and damping overwrites the out array in place, so this is the copy that survives.
Definition at line 109 of file init.hpp.
Referenced by accumulate_projector_overlaps(), and construct_density_from_states().
|
protected |
ρ_ij handed over from a checkpoint, held until the one-centre engine is built and then kept as what D's one-centre term is built from. Empty for a self-consistent run, which builds the term from its own occupations.
Definition at line 99 of file init.hpp.
Referenced by construct_density_from_states(), load_one_center_occupations(), and on_charge_density_initialized().
|
protected |
Definition at line 100 of file init.hpp.
Referenced by accumulate_projector_overlaps().
|
protected |
Definition at line 101 of file init.hpp.
Referenced by accumulate_projector_overlaps().
1.9.8 
