#include <hubbard.hpp>
Public Member Functions | |
| HubbardForceCalculator (const Crystal &crystal_structure, const PlanewaveBasis &planewave_basis, const KPointsGrid &kpoint_grid, const NormConservingPlanewaveKohnShamSCFState &scf_state) | |
| std::vector< Vector3D > | compute () const |
Definition at line 43 of file hubbard.hpp.
|
inline |
Definition at line 45 of file hubbard.hpp.
|
inline |
Ortho-atomic projectors are a joint, system-wide transform of every atomic orbital, so their gradient carries a term from orbitals that belong to neither the manifold nor the displaced atom, and is assembled separately.
Ultrasoft and PAW project through the overlap operator, so their gradient carries dS/dR as well as dphi/dR and is assembled separately.
-i(k+G) per direction – identical construction to NonlocalPseudopotentialForceCalculator's own minus_i_k_plus_g.
── Batched <phi_m|psi> and <phi_m|-i(k+G)|psi> ─────────────────────────────── These four overlap sets are the whole cost of this calculator. Each is one ZGEMM over every Hubbard orbital and occupied band at once – see compute_hubbard_projector_overlaps()'s docstring – structurally identical to NonlocalPseudopotentialForceCalculator's own "4 overlaps per band" batching (../internal/nonlocal.hpp), which this calculator already mirrors for the physics. Both operands are packed once per k-point and shared across manifolds, so neither the psi reads nor the orbital packing repeat per Hubbard atom.
Only occupied bands are contracted (as before – the scalar loop skipped the rest), and the results are scattered back into the same per-manifold overlaps matrices, leaving the MPI reduction and the gradient accumulation below untouched.
Hubbard orbitals packed column-major, (number_of_plane_waves x total_hubbard_orbitals), column row_offset + m holding manifold m's orbital.
batched_overlaps[c] is column-major (total_hubbard_orbitals x number_of_occupied_bands): component 0 is <phi|psi>, components 1..3 are <phi|-i(k+G)_direction|psi>.
Gamma-only stores one G of each +-G pair, so each of these sums is half of what it stands for. Component 0 is <phi_m|psi> and needs the G = 0 term the doubling would count twice removed; components 1-3 carry a factor of (k+G)_alpha = G_alpha, which vanishes at G = 0, so they take the doubling alone. Applying the correction to one and not the other is precisely the asymmetry that made the KB nonlocal force skew its components unevenly (KNOWN-ISSUES.md, 2026-08-08).
Scatter this manifold's row block out of the batched ZGEMM results. Rows of the unoccupied bands stay zero, exactly as the scalar loop's continue left them.
Band-parallel (G-vector-distributed) MPI: the g_index loop above only covers this rank's local plane-wave subset – reduce across the band-parallel group before use, mirroring both the Hamiltonian applicator's <phi_m|psi> reduction for H|psi> and this same fix applied to the occupation-matrix accumulation (iterator/kohn_sham/planewave/pseudopotential/init.hpp). No-op when band_size == 1.
Extract this manifold's manifold_size x manifold_size block (row-major, matching accumulate_hubbard_force_gradient's expectation) from the combined block-diagonal buffer (column-major, matching the applicator's BLAS convention – see accumulate_hubbard_occupation_and_potential's own comment).
Definition at line 54 of file hubbard.hpp.
References q3::accumulate_hubbard_force_gradient(), q3::build_atomic_orbitals_g_space(), q3::Crystal::cell_volume(), q3::compute_hubbard_projector_overlaps(), q3::convert_gamma_only_projector_overlaps_to_full_sphere(), q3::Matrix< T >::data(), q3::PlanewaveKohnShamSCFState::execution_policy, q3::find_g_zero_plane_wave_index(), q3::PlanewaveKohnShamSCFState::g_vector_indices_at_kpoint(), q3::PlanewaveBasis::g_vectors(), q3::NormConservingPlanewaveKohnShamQuantumSystem::has_ultrasoft(), q3::NormConservingPlanewaveKohnShamSCFState::hubbard_d_matrix_per_spin, q3::NormConservingPlanewaveKohnShamSCFState::hubbard_manifolds, q3::NormConservingPlanewaveKohnShamSCFState::hubbard_projector_type, q3::PlanewaveBasis::is_gamma_only(), q3::Crystal::number_of_atoms(), q3::PlanewaveKohnShamSCFState::number_of_bands(), q3::KPointsGrid::number_of_kpoints(), q3::PlanewaveKohnShamSCFState::number_of_plane_waves(), q3::PlanewaveKohnShamSCFState::number_of_spin_components(), q3::PlanewaveKohnShamSCFState::occupation(), q3::math::tolerances::occupation_threshold, q3::NormConservingPlanewaveKohnShamQuantumSystem::pseudopotential_basis(), q3::NormConservingPlanewaveBasis::pseudopotential_set(), q3::NormConservingPlanewaveKohnShamSCFState::quantum_system(), q3::Crystal::reciprocal_lattice_vectors(), q3::execution::ExecutionPolicy::reduce_within_group(), q3::tabulate_atomic_wavefunctions_in_q_space(), q3::Crystal::two_pi_over_lattice_parameter(), q3::PlanewaveKohnShamSCFState::wavefunction_coefficient(), q3::Vector3D::x, q3::Vector3D::y, and q3::Vector3D::z.
Referenced by q3::NormConservingPlanewaveForceCalculator::compute().
1.9.8 