Loading...
Searching...
No Matches
hubbard.hpp
Go to the documentation of this file.
1#pragma once
2
3/**
4 * DFT+U (Hubbard, Dudarev) Hellmann-Feynman force contribution.
5 *
6 * Computes the analytic gradient of E_U = (U/2)Tr[n(1-n)] with respect to atomic positions,
7 * using the same "atomic" (non-orthogonalized) projector convention as the energy/potential
8 * (effective_potential/internal/hubbard/hubbard.hpp). Structurally identical to
9 * NonlocalPseudopotentialForceCalculator (../internal/nonlocal.hpp): both are D-matrix-weighted
10 * sums over pairs of position-dependent G-space projectors (Kleinman-Bylander beta functions
11 * there, atomic orbitals here), differentiated via the same -i(k+G) structure-factor-derivative
12 * trick and collapsed into the same "2*Re[...]" form by the D-matrix's Hermiticity -- see
13 * hubbard.hpp's accumulate_hubbard_force_gradient() for the full derivation.
14 *
15 * The atomic-orbital projectors are rebuilt fresh here (not reused from the SCF iterator's
16 * cache), mirroring both NonlocalPseudopotentialForceCalculator's own KB projector rebuild and
17 * PlanewaveWavefunctionProjectionCalculator's (projections/planewave/init.hpp) identical
18 * pattern for the same atomic orbitals in a different post-SCF context.
19 */
20
21#include <algorithm>
22#include <stdexcept>
23#include <vector>
24
25#include "layer_0_core/logging/init.hpp"
26#include "layer_0_core/math/tolerances.hpp"
27#include "layer_0_core/types/init.hpp"
28#include "layer_0_core/types/matrix.hpp"
29#include "layer_1_structure/crystal/crystal.hpp"
30#include "layer_1_structure/reciprocal_space/kpoints_grid.hpp"
31#include "layer_3_representation/basis/planewave/init.hpp"
32#include "layer_3_representation/pseudopotential/planewave/internal/atomic_orbital_builder.hpp"
33#include "layer_3_representation/pseudopotential/planewave/internal/lowdin_orthogonalization.hpp"
34#include "layer_3_representation/pseudopotential/planewave/internal/projector.hpp"
35#include "layer_3_representation/pseudopotential/planewave/internal/wavefunction_atomic.hpp"
36#include "layer_3_representation/quantum_system/kohn_sham/planewave/ultrasoft/init.hpp"
37#include "layer_4_computation/diagonalizer/planewave/internal/overlap/pseudopotential/factory.hpp"
38#include "layer_4_computation/effective_potential/internal/hubbard/hubbard.hpp"
39#include "layer_4_computation/state/kohn_sham/planewave/norm_conserving/init.hpp"
40
41namespace q3 {
42
43class HubbardForceCalculator {
44 public:
45 HubbardForceCalculator(const Crystal& crystal_structure,
46 const PlanewaveBasis& planewave_basis,
47 const KPointsGrid& kpoint_grid,
48 const NormConservingPlanewaveKohnShamSCFState& scf_state)
49 : crystal_structure_(crystal_structure)
50 , planewave_basis_(planewave_basis)
51 , kpoint_grid_(kpoint_grid)
52 , scf_state_(scf_state) {}
53
54 [[nodiscard]] std::vector<Vector3D> compute() const {
55 int number_of_atoms = crystal_structure_.number_of_atoms();
56 std::vector<Vector3D> hubbard_force(number_of_atoms, {0.0, 0.0, 0.0});
57
58 const auto& manifolds = scf_state_.hubbard_manifolds;
59 if (manifolds.empty()) {
60 return hubbard_force;
61 }
62 /**
63 * Ortho-atomic projectors are a joint, system-wide transform of every
64 * atomic orbital, so their gradient carries a term from orbitals that
65 * belong to neither the manifold nor the displaced atom, and is
66 * assembled separately.
67 */
68 if (scf_state_.hubbard_projector_type == "ortho-atomic") {
69 return compute_with_orthogonal_projectors();
70 }
71
72 /**
73 * Ultrasoft and PAW project through the overlap operator, so their
74 * gradient carries dS/dR as well as dphi/dR and is assembled
75 * separately.
76 */
77 if (scf_state_.quantum_system().has_ultrasoft()) {
78 return compute_with_overlap_operator();
79 }
80
81 const auto& pseudopotential_set =
82 scf_state_.quantum_system().pseudopotential_basis().pseudopotential_set();
83
84 constexpr Real q_spacing = 0.01;
85 constexpr int number_of_q_points = 900;
86 auto atomic_wavefunction_table = tabulate_atomic_wavefunctions_in_q_space(
87 pseudopotential_set, crystal_structure_.cell_volume(), q_spacing, number_of_q_points);
88
89 size_t number_of_manifolds = manifolds.size();
90 std::vector<int> manifold_row_offsets(number_of_manifolds);
91 int total_hubbard_orbitals = 0;
92 for (size_t manifold_index = 0; manifold_index < number_of_manifolds; manifold_index++) {
93 manifold_row_offsets[manifold_index] = total_hubbard_orbitals;
94 total_hubbard_orbitals +=
95 static_cast<int>(manifolds[manifold_index].orbital_indices.size());
96 }
97
98 int number_of_kpoints = kpoint_grid_.number_of_kpoints();
99 int number_of_bands = scf_state_.number_of_bands();
100 int number_of_kpoint_states = number_of_kpoints * scf_state_.number_of_spin_components();
101 const Matrix3x3& reciprocal_lattice_vectors =
102 crystal_structure_.reciprocal_lattice_vectors();
103
104 for (int kpoint_index_state = 0; kpoint_index_state < number_of_kpoint_states;
105 kpoint_index_state++) {
106 int spin_index = kpoint_index_state / number_of_kpoints;
107 int kpoint_index = kpoint_index_state % number_of_kpoints;
108 int number_of_plane_waves = scf_state_.number_of_plane_waves(kpoint_index);
109 const auto& plane_wave_indices = scf_state_.g_vector_indices_at_kpoint(kpoint_index);
110 Real kpoint_weight = kpoint_grid_[kpoint_index].weight;
111 bool is_gamma_only = planewave_basis_.is_gamma_only();
112 int g_zero_index = is_gamma_only ? find_g_zero_plane_wave_index(planewave_basis_,
113 plane_wave_indices,
114 number_of_plane_waves)
115 : -1;
116
117 Vector3D kpoint_vector = kpoint_grid_[kpoint_index].coordinates;
118 Vector3D kpoint_cartesian;
119 for (int direction = 0; direction < 3; direction++) {
120 kpoint_cartesian[direction] =
121 kpoint_vector.x * reciprocal_lattice_vectors[0][direction] +
122 kpoint_vector.y * reciprocal_lattice_vectors[1][direction] +
123 kpoint_vector.z * reciprocal_lattice_vectors[2][direction];
124 }
125
126 std::vector<std::vector<Complex>> all_atomic_orbitals;
127 build_atomic_orbitals_g_space(kpoint_index,
128 number_of_plane_waves,
129 plane_wave_indices,
130 crystal_structure_,
131 kpoint_grid_,
132 planewave_basis_,
133 pseudopotential_set,
134 atomic_wavefunction_table,
135 q_spacing,
136 number_of_q_points,
137 all_atomic_orbitals);
138
139 /**
140 * -i(k+G) per direction -- identical construction to
141 * NonlocalPseudopotentialForceCalculator's own minus_i_k_plus_g.
142 */
143 std::vector<std::vector<Complex>> minus_i_k_plus_g(
144 3, std::vector<Complex>(number_of_plane_waves));
145 for (int g_index = 0; g_index < number_of_plane_waves; g_index++) {
146 const auto& g_vector = planewave_basis_.g_vectors()[plane_wave_indices[g_index]];
147 for (int direction = 0; direction < 3; direction++) {
148 Real k_plus_g_direction =
149 (g_vector.g_cart[direction] + kpoint_cartesian[direction]) *
150 crystal_structure_.two_pi_over_lattice_parameter();
151 minus_i_k_plus_g[direction][g_index] = Complex(0.0, -k_plus_g_direction);
152 }
153 }
154
155 std::vector<Real> band_occupations(number_of_bands);
156 for (int band_index = 0; band_index < number_of_bands; band_index++) {
157 band_occupations[band_index] =
158 scf_state_.occupation(kpoint_index_state, band_index);
159 }
160
161 /**
162 * ── Batched <phi_m|psi> and <phi_m|-i(k+G)|psi> ───────────────────────────────
163 * These four overlap sets are the whole cost of this calculator. Each is one
164 * ZGEMM over every Hubbard orbital and occupied band at once -- see
165 * compute_hubbard_projector_overlaps()'s docstring -- structurally identical to
166 * NonlocalPseudopotentialForceCalculator's own "4 overlaps per band" batching
167 * (../internal/nonlocal.hpp), which this calculator already mirrors for the
168 * physics. Both operands are packed once per k-point and shared across
169 * manifolds, so neither the psi reads nor the orbital packing repeat per
170 * Hubbard atom.
171 *
172 * Only occupied bands are contracted (as before -- the scalar loop skipped the
173 * rest), and the results are scattered back into the same per-manifold
174 * `overlaps` matrices, leaving the MPI reduction and the gradient accumulation
175 * below untouched.
176 */
177 std::vector<int> occupied_band_indices;
178 for (int band_index = 0; band_index < number_of_bands; band_index++) {
179 if (band_occupations[band_index] >= math::tolerances::occupation_threshold) {
180 occupied_band_indices.push_back(band_index);
181 }
182 }
183 int number_of_occupied_bands = static_cast<int>(occupied_band_indices.size());
184
185 /**
186 * Hubbard orbitals packed column-major, (number_of_plane_waves x
187 * total_hubbard_orbitals), column `row_offset + m` holding manifold m's orbital.
188 */
189 std::vector<Complex> packed_orbitals(static_cast<size_t>(number_of_plane_waves) *
190 total_hubbard_orbitals);
191 for (size_t manifold_index = 0; manifold_index < number_of_manifolds;
192 manifold_index++) {
193 const auto& manifold = manifolds[manifold_index];
194 int row_offset = manifold_row_offsets[manifold_index];
195 for (size_t m = 0; m < manifold.orbital_indices.size(); m++) {
196 const auto& orbital =
197 all_atomic_orbitals[static_cast<size_t>(manifold.orbital_indices[m])];
198 std::copy(orbital.begin(),
199 orbital.begin() + number_of_plane_waves,
200 packed_orbitals.begin() +
201 (static_cast<size_t>(row_offset) + m) * number_of_plane_waves);
202 }
203 }
204
205 /**
206 * batched_overlaps[c] is column-major (total_hubbard_orbitals x
207 * number_of_occupied_bands): component 0 is <phi|psi>, components 1..3 are
208 * <phi|-i(k+G)_direction|psi>.
209 */
210 std::vector<std::vector<Complex>> batched_overlaps(4);
211 {
212 std::vector<Complex> modified_wavefunction(
213 static_cast<size_t>(number_of_plane_waves) * number_of_occupied_bands);
214 for (int component = 0; component < 4; component++) {
215 for (int occupied_band_index = 0;
216 occupied_band_index < number_of_occupied_bands;
217 occupied_band_index++) {
218 int band_index = occupied_band_indices[occupied_band_index];
219 size_t offset =
220 static_cast<size_t>(occupied_band_index) * number_of_plane_waves;
221 for (int g_index = 0; g_index < number_of_plane_waves; g_index++) {
222 Complex psi = scf_state_.wavefunction_coefficient(
223 kpoint_index_state, band_index, g_index);
224 modified_wavefunction[offset + g_index] =
225 (component == 0) ? psi
226 : psi * minus_i_k_plus_g[component - 1][g_index];
227 }
228 }
229 compute_hubbard_projector_overlaps(
230 number_of_plane_waves,
231 total_hubbard_orbitals,
232 number_of_occupied_bands,
233 packed_orbitals.data(),
234 modified_wavefunction.data(),
235 batched_overlaps[static_cast<size_t>(component)]);
236
237 /**
238 * Gamma-only stores one G of each +-G pair, so each of
239 * these sums is half of what it stands for. Component 0 is
240 * <phi_m|psi> and needs the G = 0 term the doubling would
241 * count twice removed; components 1-3 carry a factor of
242 * (k+G)_alpha = G_alpha, which vanishes at G = 0, so they
243 * take the doubling alone. Applying the correction to one
244 * and not the other is precisely the asymmetry that made
245 * the KB nonlocal force skew its components unevenly
246 * (KNOWN-ISSUES.md, 2026-08-08).
247 */
248 if (is_gamma_only) {
249 convert_gamma_only_projector_overlaps_to_full_sphere(
250 number_of_plane_waves,
251 total_hubbard_orbitals,
252 number_of_occupied_bands,
253 packed_orbitals.data(),
254 modified_wavefunction.data(),
255 g_zero_index,
256 /*band_functions_carry_g_component=*/component != 0,
257 batched_overlaps[static_cast<size_t>(component)]);
258 }
259 }
260 }
261
262 for (size_t manifold_index = 0; manifold_index < number_of_manifolds;
263 manifold_index++) {
264 const auto& manifold = manifolds[manifold_index];
265 int manifold_size = static_cast<int>(manifold.orbital_indices.size());
266 int row_offset = manifold_row_offsets[manifold_index];
267
268 Matrix<Complex> overlaps(static_cast<size_t>(number_of_bands),
269 static_cast<size_t>(manifold_size));
270 std::vector<Matrix<Complex>> direction_derivative_overlaps(
271 3,
272 Matrix<Complex>(static_cast<size_t>(number_of_bands),
273 static_cast<size_t>(manifold_size)));
274
275 /**
276 * Scatter this manifold's row block out of the batched ZGEMM results.
277 * Rows of the unoccupied bands stay zero, exactly as the scalar loop's
278 * `continue` left them.
279 */
280 for (int occupied_band_index = 0; occupied_band_index < number_of_occupied_bands;
281 occupied_band_index++) {
282 size_t band_index =
283 static_cast<size_t>(occupied_band_indices[occupied_band_index]);
284 size_t column_offset =
285 static_cast<size_t>(occupied_band_index) * total_hubbard_orbitals;
286 for (int m = 0; m < manifold_size; m++) {
287 size_t packed_index = column_offset + static_cast<size_t>(row_offset) +
288 static_cast<size_t>(m);
289 overlaps(band_index, static_cast<size_t>(m)) =
290 batched_overlaps[0][packed_index];
291 for (int direction = 0; direction < 3; direction++) {
292 direction_derivative_overlaps[direction](band_index,
293 static_cast<size_t>(m)) =
294 batched_overlaps[static_cast<size_t>(direction) + 1][packed_index];
295 }
296 }
297 }
298
299 /**
300 * Band-parallel (G-vector-distributed) MPI: the g_index loop above only
301 * covers this rank's local plane-wave subset -- reduce across the
302 * band-parallel group before use, mirroring both the Hamiltonian
303 * applicator's <phi_m|psi> reduction for H|psi> and this same fix applied to
304 * the occupation-matrix accumulation
305 * (iterator/kohn_sham/planewave/pseudopotential/init.hpp). No-op when
306 * band_size == 1.
307 */
308 if (scf_state_.execution_policy != nullptr) {
309 scf_state_.execution_policy->reduce_within_group(
310 overlaps.data(), number_of_bands * manifold_size);
311 for (int direction = 0; direction < 3; direction++) {
312 scf_state_.execution_policy->reduce_within_group(
313 direction_derivative_overlaps[direction].data(),
314 number_of_bands * manifold_size);
315 }
316 }
317
318 /**
319 * Extract this manifold's manifold_size x manifold_size block (row-major,
320 * matching accumulate_hubbard_force_gradient's expectation) from the combined
321 * block-diagonal buffer (column-major, matching the applicator's BLAS
322 * convention -- see accumulate_hubbard_occupation_and_potential's own comment).
323 */
324 const auto& combined_d_matrix = scf_state_.hubbard_d_matrix_per_spin[spin_index];
325 std::vector<Complex> manifold_d_matrix(static_cast<size_t>(manifold_size) *
326 manifold_size);
327 for (int m = 0; m < manifold_size; m++) {
328 for (int m_prime = 0; m_prime < manifold_size; m_prime++) {
329 int global_row = row_offset + m;
330 int global_column = row_offset + m_prime;
331 manifold_d_matrix[static_cast<size_t>(m * manifold_size + m_prime)] =
332 combined_d_matrix[static_cast<size_t>(global_row) +
333 static_cast<size_t>(global_column) *
334 total_hubbard_orbitals];
335 }
336 }
337
338 for (int direction = 0; direction < 3; direction++) {
339 Real force_contribution =
340 accumulate_hubbard_force_gradient(number_of_bands,
341 manifold_size,
342 overlaps,
343 direction_derivative_overlaps[direction],
344 band_occupations,
345 kpoint_weight,
346 manifold_d_matrix);
347 hubbard_force[manifold.atom_index][direction] += force_contribution;
348 }
349 }
350 }
351
352 return hubbard_force;
353 }
354
355 private:
356 /**
357 * The DFT+U force when the projector is the Löwdin-orthogonalized orbital.
358 *
359 * ── What is being differentiated ──────────────────────────────────────
360 *
361 * The ortho-atomic projector is not a property of its own atom. The whole
362 * system-wide atomic-orbital set {phi_J} is orthogonalized jointly,
363 *
364 * phi^ortho_I = sum_J (O^-1/2)_JI phi_J, O_IJ = <phi_I|phi_J>,
365 *
366 * so the projection whose occupation DFT+U is built from,
367 *
368 * P_Iv = <phi^ortho_I|psi_v> = sum_J (O^-1/2)_IJ <phi_J|psi_v>,
369 *
370 * depends on the position of *every* atom in the cell through O. That is
371 * the term the atomic path (compute() above) has no analogue of: there a
372 * Hubbard force appears only on a Hubbard atom, here displacing any atom
373 * exerts one on every manifold. It is Pulay-like in origin -- a
374 * position-dependent basis -- but it is not the usual Pulay term, because
375 * it is the *projector* basis that moves, not the wavefunction basis.
376 *
377 * Differentiating, with alpha the displaced atom and its own orbitals the
378 * only ones with a non-zero dphi/dtau,
379 *
380 * dP_Iv/dtau = sum_J d(O^-1/2)_IJ <phi_J|psi_v>
381 * + sum_{J in alpha} (O^-1/2)_IJ <dphi_J/dtau|psi_v>
382 *
383 * (using that O^-1/2 and its derivative are both Hermitian, so the
384 * conjugate of the (J,I) element is the (I,J) element). dO/dtau follows
385 * from the same product rule,
386 *
387 * dO_IJ/dtau = [J in alpha] T_IJ - [I in alpha] T_IJ,
388 * T_IJ = <phi_I|-i(k+G)|phi_J>,
389 *
390 * which is antisymmetric enough to vanish identically when I and J both
391 * belong to alpha -- rigidly translating an atom cannot change the
392 * overlaps among its own orbitals. d(O^-1/2) then comes from dO and O's
393 * eigendecomposition in closed form (lowdin_inverse_square_root_derivative),
394 * which evaluates Eq. (32) of Phys. Rev. B 102, 235159 (2020).
395 *
396 * ── Sign convention ───────────────────────────────────────────────────
397 *
398 * `accumulate_hubbard_force_gradient` is shared with the norm-conserving
399 * atomic path, which hands it <phi_m|-i(k+G)|psi>, i.e. **-dP/dtau** (see
400 * compute_with_overlap_operator's own note). The same quantity is passed
401 * here:
402 *
403 * -dP_Iv/dtau = -sum_J d(O^-1/2)_IJ <phi_J|psi_v>
404 * + sum_{J in alpha} (O^-1/2)_IJ T^psi_Jv
405 *
406 * with T^psi_Jv = <phi_J|-i(k+G)|psi_v>, using
407 * <dphi_J/dtau|psi> = -T^psi_Jv. Freeze O^-1/2 at the identity and drop
408 * its derivative and this collapses to T^psi_Iv on the manifold's own
409 * atom and to zero everywhere else -- exactly what the atomic path passes.
410 * Pinning the sign by reduction to an already-validated path is safer
411 * than re-deriving it from the energy.
412 *
413 * For ultrasoft and PAW, O is the generalized <phi|S|phi> metric. The
414 * additional branch below differentiates both S occurrences explicitly;
415 * the norm-conserving algebra remains the simpler S=1 specialization.
416 */
417 [[nodiscard]] std::vector<Vector3D> compute_with_orthogonal_projectors() const {
418 int number_of_atoms = crystal_structure_.number_of_atoms();
419 std::vector<Vector3D> hubbard_force(number_of_atoms, {0.0, 0.0, 0.0});
420
421 const auto& manifolds = scf_state_.hubbard_manifolds;
422 const auto& pseudopotential_basis = scf_state_.quantum_system().pseudopotential_basis();
423 const auto& pseudopotential_set = pseudopotential_basis.pseudopotential_set();
424 bool has_augmented_overlap = scf_state_.quantum_system().has_ultrasoft();
425
426 constexpr Real q_spacing = 0.01;
427 constexpr int number_of_q_points = 900;
428 auto atomic_wavefunction_table = tabulate_atomic_wavefunctions_in_q_space(
429 pseudopotential_set, crystal_structure_.cell_volume(), q_spacing, number_of_q_points);
430
431 /**
432 * Where each atom's own atomic orbitals sit in the system-wide list.
433 * enumerate_atomic_orbitals is atom-major, so each atom owns one
434 * contiguous run, and displacing that atom is the only thing that gives
435 * the orbitals in it a non-zero dphi/dtau. An atom contributing no
436 * orbitals is left with first > last, making every loop over its range
437 * empty.
438 */
439 auto orbital_labels = enumerate_atomic_orbitals(crystal_structure_, pseudopotential_set);
440 int number_of_atomic_orbitals = static_cast<int>(orbital_labels.size());
441 std::vector<int> atom_orbital_first(number_of_atoms, number_of_atomic_orbitals);
442 std::vector<int> atom_orbital_last(number_of_atoms, 0);
443 for (int orbital = 0; orbital < number_of_atomic_orbitals; orbital++) {
444 int atom_index = orbital_labels[static_cast<size_t>(orbital)].atom_index;
445 atom_orbital_first[atom_index] = std::min(atom_orbital_first[atom_index], orbital);
446 atom_orbital_last[atom_index] = std::max(atom_orbital_last[atom_index], orbital + 1);
447 }
448
449 int total_projectors = pseudopotential_basis.total_number_of_projectors();
450 std::vector<int> projector_offset_by_atom(number_of_atoms, 0);
451 std::vector<int> projector_channels_by_atom(number_of_atoms, 0);
452 {
453 int projector_offset = 0;
454 for (int atom_index = 0; atom_index < number_of_atoms; atom_index++) {
455 int species_index = crystal_structure_.species_type_index(atom_index);
456 int number_of_projector_channels =
457 pseudopotential_basis.number_of_projector_channels(species_index);
458 projector_offset_by_atom[atom_index] = projector_offset;
459 projector_channels_by_atom[atom_index] = number_of_projector_channels;
460 projector_offset += number_of_projector_channels;
461 }
462 }
463 auto augmentation_overlap = [&](int atom_index, int channel_i, int channel_j) -> Real {
464 return hubbard_augmentation_overlap(
465 static_cast<const UltrasoftPlanewaveKohnShamQuantumSystem&>(
466 scf_state_.quantum_system())
467 .ultrasoft_basis(),
468 pseudopotential_basis,
469 crystal_structure_.species_type_index(atom_index),
470 channel_i,
471 channel_j);
472 };
473
474 size_t number_of_manifolds = manifolds.size();
475 std::vector<int> manifold_row_offsets(number_of_manifolds);
476 int total_hubbard_orbitals = 0;
477 for (size_t manifold_index = 0; manifold_index < number_of_manifolds; manifold_index++) {
478 manifold_row_offsets[manifold_index] = total_hubbard_orbitals;
479 total_hubbard_orbitals +=
480 static_cast<int>(manifolds[manifold_index].orbital_indices.size());
481 }
482
483 int number_of_kpoints = kpoint_grid_.number_of_kpoints();
484 int number_of_bands = scf_state_.number_of_bands();
485 int number_of_kpoint_states = number_of_kpoints * scf_state_.number_of_spin_components();
486 const Matrix3x3& reciprocal_lattice_vectors =
487 crystal_structure_.reciprocal_lattice_vectors();
488
489 for (int kpoint_index_state = 0; kpoint_index_state < number_of_kpoint_states;
490 kpoint_index_state++) {
491 int spin_index = kpoint_index_state / number_of_kpoints;
492 int kpoint_index = kpoint_index_state % number_of_kpoints;
493 int number_of_plane_waves = scf_state_.number_of_plane_waves(kpoint_index);
494 const auto& plane_wave_indices = scf_state_.g_vector_indices_at_kpoint(kpoint_index);
495 Real kpoint_weight = kpoint_grid_[kpoint_index].weight;
496 bool is_gamma_only = planewave_basis_.is_gamma_only();
497 int g_zero_index = is_gamma_only ? find_g_zero_plane_wave_index(planewave_basis_,
498 plane_wave_indices,
499 number_of_plane_waves)
500 : -1;
501
502 Vector3D kpoint_cartesian;
503 {
504 Vector3D kpoint_vector = kpoint_grid_[kpoint_index].coordinates;
505 for (int direction = 0; direction < 3; direction++)
506 kpoint_cartesian[direction] =
507 kpoint_vector.x * reciprocal_lattice_vectors[0][direction] +
508 kpoint_vector.y * reciprocal_lattice_vectors[1][direction] +
509 kpoint_vector.z * reciprocal_lattice_vectors[2][direction];
510 }
511
512 std::vector<std::vector<Complex>> all_atomic_orbitals;
513 build_atomic_orbitals_g_space(kpoint_index,
514 number_of_plane_waves,
515 plane_wave_indices,
516 crystal_structure_,
517 kpoint_grid_,
518 planewave_basis_,
519 pseudopotential_set,
520 atomic_wavefunction_table,
521 q_spacing,
522 number_of_q_points,
523 all_atomic_orbitals);
524
525 std::vector<std::vector<Complex>> minus_i_k_plus_g(
526 3, std::vector<Complex>(number_of_plane_waves));
527 for (int g_index = 0; g_index < number_of_plane_waves; g_index++) {
528 const auto& g_vector = planewave_basis_.g_vectors()[plane_wave_indices[g_index]];
529 for (int direction = 0; direction < 3; direction++)
530 minus_i_k_plus_g[direction][g_index] =
531 Complex(0.0,
532 -(g_vector.g_cart[direction] + kpoint_cartesian[direction]) *
533 crystal_structure_.two_pi_over_lattice_parameter());
534 }
535
536 std::vector<Real> band_occupations(number_of_bands);
537 std::vector<int> occupied_band_indices;
538 for (int band_index = 0; band_index < number_of_bands; band_index++) {
539 band_occupations[band_index] =
540 scf_state_.occupation(kpoint_index_state, band_index);
541 if (band_occupations[band_index] >= math::tolerances::occupation_threshold)
542 occupied_band_indices.push_back(band_index);
543 }
544 int number_of_occupied_bands = static_cast<int>(occupied_band_indices.size());
545
546 /**
547 * Column-major (number_of_plane_waves x N) operands, the layout
548 * compute_hubbard_projector_overlaps' ZGEMM expects, plus their
549 * -i(k+G)-weighted counterparts. Unlike the atomic path this needs
550 * the FULL orbital set, not just the Hubbard manifolds' columns:
551 * O and its derivative run over every atomic orbital in the cell.
552 */
553 std::vector<Complex> packed_orbitals(static_cast<size_t>(number_of_plane_waves) *
554 number_of_atomic_orbitals);
555 for (int orbital = 0; orbital < number_of_atomic_orbitals; orbital++)
556 std::copy(all_atomic_orbitals[static_cast<size_t>(orbital)].begin(),
557 all_atomic_orbitals[static_cast<size_t>(orbital)].begin() +
558 number_of_plane_waves,
559 packed_orbitals.begin() +
560 static_cast<size_t>(orbital) * number_of_plane_waves);
561 std::vector<Complex> packed_bands(static_cast<size_t>(number_of_plane_waves) *
562 number_of_occupied_bands);
563 for (int occupied_band_index = 0; occupied_band_index < number_of_occupied_bands;
564 occupied_band_index++)
565 for (int g_index = 0; g_index < number_of_plane_waves; g_index++)
566 packed_bands[static_cast<size_t>(occupied_band_index) * number_of_plane_waves +
567 g_index] =
568 scf_state_.wavefunction_coefficient(
569 kpoint_index_state,
570 occupied_band_indices[static_cast<size_t>(occupied_band_index)],
571 g_index);
572
573 auto weight_by_direction = [&](const std::vector<Complex>& source, int direction) {
574 std::vector<Complex> weighted(source.size());
575 size_t columns = source.size() / static_cast<size_t>(number_of_plane_waves);
576 for (size_t column = 0; column < columns; column++)
577 for (int g_index = 0; g_index < number_of_plane_waves; g_index++) {
578 size_t offset = column * number_of_plane_waves + g_index;
579 weighted[offset] = source[offset] * minus_i_k_plus_g[direction][g_index];
580 }
581 return weighted;
582 };
583
584 /**
585 * Every contraction is a partial sum under band (G-vector)
586 * parallelism and is reduced before use, the same reason the SCF
587 * reduces the occupation matrix at its own G-sum. Gamma-only
588 * contractions are first reconstructed from the stored half sphere;
589 * derivative operands carry a G component and therefore no G=0 term.
590 */
591 auto contract = [&](const Complex* left,
592 int number_of_left,
593 const Complex* right,
594 int number_of_right,
595 bool right_carries_g_component = false) {
596 std::vector<Complex> result;
597 compute_hubbard_projector_overlaps(
598 number_of_plane_waves, number_of_left, number_of_right, left, right, result);
599 if (is_gamma_only)
600 convert_gamma_only_projector_overlaps_to_full_sphere(number_of_plane_waves,
601 number_of_left,
602 number_of_right,
603 left,
604 right,
605 g_zero_index,
606 right_carries_g_component,
607 result);
608 if (scf_state_.execution_policy != nullptr)
609 scf_state_.execution_policy->reduce_within_group(
610 result.data(), static_cast<int>(result.size()));
611 return result;
612 };
613
614 if (has_augmented_overlap) {
615 if (scf_state_.execution_policy == nullptr) {
616 throw std::runtime_error(
617 "Ortho-atomic ultrasoft Hubbard forces require an execution policy");
618 }
619
620 Matrix<Complex> kleinman_bylander_projectors;
621 compute_kleinman_bylander_projectors(kpoint_grid_[kpoint_index].coordinates,
622 crystal_structure_,
623 planewave_basis_,
624 pseudopotential_set,
625 pseudopotential_basis.beta_function_table(),
626 plane_wave_indices,
627 number_of_plane_waves,
628 total_projectors,
629 kleinman_bylander_projectors);
630 auto overlap_operator = create_overlap_operator(crystal_structure_,
631 scf_state_.quantum_system(),
632 kleinman_bylander_projectors,
633 *scf_state_.execution_policy);
634
635 std::vector<Complex> overlap_applied_orbitals(packed_orbitals.size());
636 overlap_operator->apply_batch(number_of_plane_waves,
637 number_of_atomic_orbitals,
638 packed_orbitals.data(),
639 overlap_applied_orbitals.data());
640 std::vector<Complex> overlap_applied_bands(packed_bands.size());
641 overlap_operator->apply_batch(number_of_plane_waves,
642 number_of_occupied_bands,
643 packed_bands.data(),
644 overlap_applied_bands.data());
645
646 auto overlap_matrix = contract(packed_orbitals.data(),
647 number_of_atomic_orbitals,
648 overlap_applied_orbitals.data(),
649 number_of_atomic_orbitals);
650 auto transform =
651 compute_lowdin_transform(overlap_matrix, number_of_atomic_orbitals);
652 auto atomic_band_overlaps = contract(packed_orbitals.data(),
653 number_of_atomic_orbitals,
654 overlap_applied_bands.data(),
655 number_of_occupied_bands);
656 auto element = [number_of_atomic_orbitals](int row, int column) {
657 return static_cast<size_t>(row + column * number_of_atomic_orbitals);
658 };
659
660 // P_Av = sum_I (O^-1/2)_AI <phi_I|S|psi_v>.
661 std::vector<Complex> projection(static_cast<size_t>(number_of_atomic_orbitals) *
662 number_of_occupied_bands,
663 Complex(0.0, 0.0));
664 for (int occupied_band_index = 0; occupied_band_index < number_of_occupied_bands;
665 occupied_band_index++) {
666 for (int orbital = 0; orbital < number_of_atomic_orbitals; orbital++) {
667 for (int inner_orbital = 0; inner_orbital < number_of_atomic_orbitals;
668 inner_orbital++) {
669 projection[element(orbital, occupied_band_index)] +=
670 transform.inverse_square_root[element(orbital, inner_orbital)] *
671 atomic_band_overlaps[element(inner_orbital, occupied_band_index)];
672 }
673 }
674 }
675
676 /**
677 * Apply dS/dR to a column-major vector block. Both beta overlaps are
678 * completed across the band group before the local G-space result is
679 * assembled, matching the production overlap operator's collective order.
680 */
681 auto apply_overlap_derivative = [&](const Complex* input_vectors,
682 int number_of_vectors,
683 int displaced_atom,
684 int direction) {
685 std::vector<Complex> derivative_projectors(
686 static_cast<size_t>(number_of_plane_waves) * total_projectors,
687 Complex(0.0, 0.0));
688 int displaced_projector_offset = projector_offset_by_atom[displaced_atom];
689 int displaced_projector_channels = projector_channels_by_atom[displaced_atom];
690 for (int channel = 0; channel < displaced_projector_channels; channel++) {
691 int projector_index = displaced_projector_offset + channel;
692 for (int g_index = 0; g_index < number_of_plane_waves; g_index++) {
693 derivative_projectors[static_cast<size_t>(projector_index) *
694 number_of_plane_waves +
695 g_index] =
696 kleinman_bylander_projectors[projector_index][g_index] *
697 minus_i_k_plus_g[direction][g_index];
698 }
699 }
700
701 auto projector_overlaps = contract(kleinman_bylander_projectors.data(),
702 total_projectors,
703 input_vectors,
704 number_of_vectors);
705 auto derivative_projector_overlaps = contract(derivative_projectors.data(),
706 total_projectors,
707 input_vectors,
708 number_of_vectors);
709 std::vector<Complex> result(static_cast<size_t>(number_of_plane_waves) *
710 number_of_vectors,
711 Complex(0.0, 0.0));
712 for (int channel_row = 0; channel_row < displaced_projector_channels;
713 channel_row++) {
714 int projector_row = displaced_projector_offset + channel_row;
715 for (int channel_column = 0; channel_column < displaced_projector_channels;
716 channel_column++) {
717 Real overlap_value =
718 augmentation_overlap(displaced_atom, channel_row, channel_column);
719 if (std::abs(overlap_value) <
720 math::tolerances::d_matrix_element_threshold) {
721 continue;
722 }
723 int projector_column = displaced_projector_offset + channel_column;
724 for (int vector_index = 0; vector_index < number_of_vectors;
725 vector_index++) {
726 size_t overlap_index =
727 static_cast<size_t>(projector_column) +
728 static_cast<size_t>(vector_index) * total_projectors;
729 Complex derivative_projector_coefficient =
730 overlap_value * projector_overlaps[overlap_index];
731 Complex projector_coefficient =
732 overlap_value * derivative_projector_overlaps[overlap_index];
733 for (int g_index = 0; g_index < number_of_plane_waves; g_index++) {
734 size_t result_index =
735 static_cast<size_t>(vector_index) * number_of_plane_waves +
736 g_index;
737 result[result_index] +=
738 derivative_projectors[static_cast<size_t>(projector_row) *
739 number_of_plane_waves +
740 g_index] *
741 derivative_projector_coefficient +
742 kleinman_bylander_projectors[projector_row][g_index] *
743 projector_coefficient;
744 }
745 }
746 }
747 }
748 return result;
749 };
750
751 for (int displaced_atom = 0; displaced_atom < number_of_atoms; displaced_atom++) {
752 int displaced_first = atom_orbital_first[displaced_atom];
753 int displaced_last = atom_orbital_last[displaced_atom];
754
755 for (int direction = 0; direction < 3; direction++) {
756 std::vector<Complex> orbital_derivatives(packed_orbitals.size(),
757 Complex(0.0, 0.0));
758 for (int orbital = displaced_first; orbital < displaced_last; orbital++) {
759 for (int g_index = 0; g_index < number_of_plane_waves; g_index++) {
760 size_t orbital_index =
761 static_cast<size_t>(orbital) * number_of_plane_waves + g_index;
762 orbital_derivatives[orbital_index] =
763 packed_orbitals[orbital_index] *
764 minus_i_k_plus_g[direction][g_index];
765 }
766 }
767
768 std::vector<Complex> overlap_applied_orbital_derivatives(
769 orbital_derivatives.size());
770 overlap_operator->apply_batch(number_of_plane_waves,
771 number_of_atomic_orbitals,
772 orbital_derivatives.data(),
773 overlap_applied_orbital_derivatives.data());
774 auto overlap_derivative_on_orbitals =
775 apply_overlap_derivative(packed_orbitals.data(),
776 number_of_atomic_orbitals,
777 displaced_atom,
778 direction);
779 auto overlap_derivative_on_bands =
780 apply_overlap_derivative(packed_bands.data(),
781 number_of_occupied_bands,
782 displaced_atom,
783 direction);
784
785 auto left_orbital_metric_derivative =
786 contract(orbital_derivatives.data(),
787 number_of_atomic_orbitals,
788 overlap_applied_orbitals.data(),
789 number_of_atomic_orbitals);
790 auto overlap_metric_derivative =
791 contract(packed_orbitals.data(),
792 number_of_atomic_orbitals,
793 overlap_derivative_on_orbitals.data(),
794 number_of_atomic_orbitals);
795 auto right_orbital_metric_derivative =
796 contract(packed_orbitals.data(),
797 number_of_atomic_orbitals,
798 overlap_applied_orbital_derivatives.data(),
799 number_of_atomic_orbitals);
800 std::vector<Complex> metric_derivative(overlap_matrix.size());
801 for (size_t matrix_index = 0; matrix_index < metric_derivative.size();
802 matrix_index++) {
803 metric_derivative[matrix_index] =
804 left_orbital_metric_derivative[matrix_index] +
805 overlap_metric_derivative[matrix_index] +
806 right_orbital_metric_derivative[matrix_index];
807 }
808 for (int column = 0; column < number_of_atomic_orbitals; column++) {
809 for (int row = 0; row <= column; row++) {
810 Complex hermitian_value =
811 0.5 * (metric_derivative[element(row, column)] +
812 std::conj(metric_derivative[element(column, row)]));
813 metric_derivative[element(row, column)] = hermitian_value;
814 metric_derivative[element(column, row)] =
815 std::conj(hermitian_value);
816 }
817 }
818 auto inverse_square_root_derivative =
819 lowdin_inverse_square_root_derivative(transform, metric_derivative);
820
821 auto left_orbital_projection_derivative =
822 contract(orbital_derivatives.data(),
823 number_of_atomic_orbitals,
824 overlap_applied_bands.data(),
825 number_of_occupied_bands);
826 auto overlap_projection_derivative =
827 contract(packed_orbitals.data(),
828 number_of_atomic_orbitals,
829 overlap_derivative_on_bands.data(),
830 number_of_occupied_bands);
831 std::vector<Complex> raw_projection_derivative(atomic_band_overlaps.size());
832 for (size_t overlap_index = 0;
833 overlap_index < raw_projection_derivative.size();
834 overlap_index++) {
835 raw_projection_derivative[overlap_index] =
836 left_orbital_projection_derivative[overlap_index] +
837 overlap_projection_derivative[overlap_index];
838 }
839
840 std::vector<Complex> negative_projection_derivative(projection.size(),
841 Complex(0.0, 0.0));
842 for (int occupied_band_index = 0;
843 occupied_band_index < number_of_occupied_bands;
844 occupied_band_index++) {
845 for (int orbital = 0; orbital < number_of_atomic_orbitals; orbital++) {
846 for (int inner_orbital = 0;
847 inner_orbital < number_of_atomic_orbitals;
848 inner_orbital++) {
849 negative_projection_derivative[element(orbital,
850 occupied_band_index)] -=
851 inverse_square_root_derivative[element(orbital,
852 inner_orbital)] *
853 atomic_band_overlaps[element(inner_orbital,
854 occupied_band_index)] +
855 transform.inverse_square_root[element(orbital,
856 inner_orbital)] *
857 raw_projection_derivative[element(inner_orbital,
858 occupied_band_index)];
859 }
860 }
861 }
862
863 for (size_t manifold_index = 0; manifold_index < number_of_manifolds;
864 manifold_index++) {
865 const auto& manifold = manifolds[manifold_index];
866 int manifold_size = static_cast<int>(manifold.orbital_indices.size());
867 int row_offset = manifold_row_offsets[manifold_index];
868 Matrix<Complex> manifold_overlaps(static_cast<size_t>(number_of_bands),
869 static_cast<size_t>(manifold_size));
870 Matrix<Complex> derivative_overlaps(
871 static_cast<size_t>(number_of_bands),
872 static_cast<size_t>(manifold_size));
873 for (int manifold_orbital = 0; manifold_orbital < manifold_size;
874 manifold_orbital++) {
875 int orbital =
876 manifold.orbital_indices[static_cast<size_t>(manifold_orbital)];
877 for (int occupied_band_index = 0;
878 occupied_band_index < number_of_occupied_bands;
879 occupied_band_index++) {
880 size_t band_index = static_cast<size_t>(
881 occupied_band_indices[static_cast<size_t>(
882 occupied_band_index)]);
883 manifold_overlaps(band_index,
884 static_cast<size_t>(manifold_orbital)) =
885 projection[element(orbital, occupied_band_index)];
886 derivative_overlaps(band_index,
887 static_cast<size_t>(manifold_orbital)) =
888 negative_projection_derivative[element(
889 orbital, occupied_band_index)];
890 }
891 }
892
893 const auto& combined_d_matrix =
894 scf_state_.hubbard_d_matrix_per_spin[spin_index];
895 std::vector<Complex> manifold_d_matrix(
896 static_cast<size_t>(manifold_size) * manifold_size);
897 for (int row = 0; row < manifold_size; row++) {
898 for (int column = 0; column < manifold_size; column++) {
899 manifold_d_matrix[static_cast<size_t>(row * manifold_size +
900 column)] =
901 combined_d_matrix[static_cast<size_t>(row_offset + row) +
902 static_cast<size_t>(row_offset + column) *
903 total_hubbard_orbitals];
904 }
905 }
906
907 hubbard_force[displaced_atom][direction] +=
908 accumulate_hubbard_force_gradient(number_of_bands,
909 manifold_size,
910 manifold_overlaps,
911 derivative_overlaps,
912 band_occupations,
913 kpoint_weight,
914 manifold_d_matrix);
915 }
916 }
917 }
918 continue;
919 }
920
921 auto overlap_matrix = contract(packed_orbitals.data(),
922 number_of_atomic_orbitals,
923 packed_orbitals.data(),
924 number_of_atomic_orbitals);
925 auto transform = compute_lowdin_transform(overlap_matrix, number_of_atomic_orbitals);
926 auto atomic_band_overlaps = contract(packed_orbitals.data(),
927 number_of_atomic_orbitals,
928 packed_bands.data(),
929 number_of_occupied_bands);
930 std::vector<std::vector<Complex>> orbital_derivative_overlaps(3);
931 std::vector<std::vector<Complex>> band_derivative_overlaps(3);
932 for (int direction = 0; direction < 3; direction++) {
933 auto weighted_orbitals = weight_by_direction(packed_orbitals, direction);
934 auto weighted_bands = weight_by_direction(packed_bands, direction);
935 orbital_derivative_overlaps[static_cast<size_t>(direction)] =
936 contract(packed_orbitals.data(),
937 number_of_atomic_orbitals,
938 weighted_orbitals.data(),
939 number_of_atomic_orbitals,
940 /*right_carries_g_component=*/true);
941 band_derivative_overlaps[static_cast<size_t>(direction)] =
942 contract(packed_orbitals.data(),
943 number_of_atomic_orbitals,
944 weighted_bands.data(),
945 number_of_occupied_bands,
946 /*right_carries_g_component=*/true);
947 }
948
949 auto element = [number_of_atomic_orbitals](int row, int column) {
950 return static_cast<size_t>(row + column * number_of_atomic_orbitals);
951 };
952
953 // P_Iv = sum_J (O^-1/2)_IJ <phi_J|psi_v>.
954 std::vector<Complex> projection(static_cast<size_t>(number_of_atomic_orbitals) *
955 number_of_occupied_bands,
956 Complex(0.0, 0.0));
957 for (int occupied_band_index = 0; occupied_band_index < number_of_occupied_bands;
958 occupied_band_index++)
959 for (int orbital = 0; orbital < number_of_atomic_orbitals; orbital++) {
960 Complex accumulated(0.0, 0.0);
961 for (int inner = 0; inner < number_of_atomic_orbitals; inner++)
962 accumulated += transform.inverse_square_root[element(orbital, inner)] *
963 atomic_band_overlaps[element(inner, occupied_band_index)];
964 projection[element(orbital, occupied_band_index)] = accumulated;
965 }
966
967 for (int displaced_atom = 0; displaced_atom < number_of_atoms; displaced_atom++) {
968 int displaced_first = atom_orbital_first[displaced_atom];
969 int displaced_last = atom_orbital_last[displaced_atom];
970 auto belongs_to_displaced_atom = [&](int orbital) {
971 return orbital >= displaced_first && orbital < displaced_last;
972 };
973
974 for (int direction = 0; direction < 3; direction++) {
975 const auto& orbital_derivative =
976 orbital_derivative_overlaps[static_cast<size_t>(direction)];
977 const auto& band_derivative =
978 band_derivative_overlaps[static_cast<size_t>(direction)];
979
980 std::vector<Complex> overlap_derivative(
981 static_cast<size_t>(number_of_atomic_orbitals) * number_of_atomic_orbitals,
982 Complex(0.0, 0.0));
983 for (int column = 0; column < number_of_atomic_orbitals; column++)
984 for (int row = 0; row < number_of_atomic_orbitals; row++) {
985 bool row_moves = belongs_to_displaced_atom(row);
986 bool column_moves = belongs_to_displaced_atom(column);
987 if (row_moves == column_moves)
988 continue;
989 Complex value = orbital_derivative[element(row, column)];
990 overlap_derivative[element(row, column)] =
991 column_moves ? value : -value;
992 }
993 auto inverse_square_root_derivative =
994 lowdin_inverse_square_root_derivative(transform, overlap_derivative);
995
996 for (size_t manifold_index = 0; manifold_index < number_of_manifolds;
997 manifold_index++) {
998 const auto& manifold = manifolds[manifold_index];
999 int manifold_size = static_cast<int>(manifold.orbital_indices.size());
1000 int row_offset = manifold_row_offsets[manifold_index];
1001
1002 Matrix<Complex> manifold_overlaps(static_cast<size_t>(number_of_bands),
1003 static_cast<size_t>(manifold_size));
1004 Matrix<Complex> derivative_overlaps(static_cast<size_t>(number_of_bands),
1005 static_cast<size_t>(manifold_size));
1006 for (int m = 0; m < manifold_size; m++) {
1007 int orbital = manifold.orbital_indices[static_cast<size_t>(m)];
1008 for (int occupied_band_index = 0;
1009 occupied_band_index < number_of_occupied_bands;
1010 occupied_band_index++) {
1011 Complex transform_term(0.0, 0.0);
1012 for (int inner = 0; inner < number_of_atomic_orbitals; inner++)
1013 transform_term +=
1014 inverse_square_root_derivative[element(orbital, inner)] *
1015 atomic_band_overlaps[element(inner, occupied_band_index)];
1016 Complex own_orbital_term(0.0, 0.0);
1017 for (int inner = displaced_first; inner < displaced_last; inner++)
1018 own_orbital_term +=
1019 transform.inverse_square_root[element(orbital, inner)] *
1020 band_derivative[element(inner, occupied_band_index)];
1021 size_t band_row =
1022 static_cast<size_t>(occupied_band_indices[static_cast<size_t>(
1023 occupied_band_index)]);
1024 manifold_overlaps(band_row, static_cast<size_t>(m)) =
1025 projection[element(orbital, occupied_band_index)];
1026 derivative_overlaps(band_row, static_cast<size_t>(m)) =
1027 own_orbital_term - transform_term;
1028 }
1029 }
1030
1031 const auto& combined_d_matrix =
1032 scf_state_.hubbard_d_matrix_per_spin[spin_index];
1033 std::vector<Complex> manifold_d_matrix(static_cast<size_t>(manifold_size) *
1034 manifold_size);
1035 for (int m = 0; m < manifold_size; m++)
1036 for (int m_prime = 0; m_prime < manifold_size; m_prime++)
1037 manifold_d_matrix[static_cast<size_t>(m * manifold_size +
1038 m_prime)] =
1039 combined_d_matrix[static_cast<size_t>(row_offset + m) +
1040 static_cast<size_t>(row_offset + m_prime) *
1041 total_hubbard_orbitals];
1042
1043 hubbard_force[displaced_atom][direction] +=
1044 accumulate_hubbard_force_gradient(number_of_bands,
1045 manifold_size,
1046 manifold_overlaps,
1047 derivative_overlaps,
1048 band_occupations,
1049 kpoint_weight,
1050 manifold_d_matrix);
1051 }
1052 }
1053 }
1054 }
1055
1056 return hubbard_force;
1057 }
1058
1059 /**
1060 * The DFT+U force when the projector is S|phi> rather than |phi>.
1061 *
1062 * ── What is being differentiated ──────────────────────────────────────
1063 *
1064 * For ultrasoft and PAW the occupation is built from
1065 *
1066 * P_m = <phi_m|S|psi>, S = 1 + sum_ij q_ij |beta_i><beta_j|
1067 *
1068 * so displacing an atom moves three things: the Hubbard orbital itself
1069 * (only when the displaced atom *is* the Hubbard atom, since these
1070 * orbitals are localized), and both beta projectors inside S (whenever the
1071 * displaced atom carries projectors, Hubbard or not). That last part is
1072 * why this cannot reuse the norm-conserving loop: there, a Hubbard force
1073 * only ever appears on a Hubbard atom; here, displacing an ordinary atom
1074 * exerts one too.
1075 *
1076 * S|phi> is never built. Expanding it,
1077 *
1078 * P_m = <phi_m|psi> + sum_ij q_ij <phi_m|beta_i><beta_j|psi>
1079 *
1080 * needs only overlaps this file already forms plus <phi|beta>, which is
1081 * cheaper than carrying an augmented projector around and cannot drift
1082 * from the SCF's own definition.
1083 *
1084 * ── Sign convention ───────────────────────────────────────────────────
1085 *
1086 * `accumulate_hubbard_force_gradient` is shared with the norm-conserving
1087 * path, which hands it <phi_m|-i(k+G)|psi>. Writing T[X]_m for
1088 * <phi_m|-i(k+G)|X>, and using dphi/dtau = -i(k+G)phi (so
1089 * <dphi_m|X> = -T[X]_m) and the same for beta, the derivative of P is
1090 *
1091 * dP_m/dtau = -delta_a ( T[psi]_m + sum_ij q_ij T[beta_i]_m <beta_j|psi> )
1092 * + sum_{ij in alpha} q_ij ( T[beta_i]_m <beta_j|psi>
1093 * - <phi_m|beta_i> Bd_j )
1094 *
1095 * with delta_a = 1 when the displaced atom owns the manifold and
1096 * Bd_j = <beta_j|-i(k+G)|psi>. What the contraction wants is **-dP/dtau**:
1097 * that is exactly T[psi]_m when q_ij vanishes, which is what the
1098 * norm-conserving path passes today. Deriving the sign by reduction to an
1099 * already-validated path is safer than re-deriving it from the energy.
1100 *
1101 * The k in (k+G) cancels between the two halves of dn, so using (k+G)
1102 * throughout — as the norm-conserving path does — rather than QE's bare G
1103 * for the orbital term gives the same force.
1104 */
1105 [[nodiscard]] std::vector<Vector3D> compute_with_overlap_operator() const {
1106 int number_of_atoms = crystal_structure_.number_of_atoms();
1107 std::vector<Vector3D> hubbard_force(number_of_atoms, {0.0, 0.0, 0.0});
1108
1109 const auto& manifolds = scf_state_.hubbard_manifolds;
1110 const auto& basis = scf_state_.quantum_system().pseudopotential_basis();
1111 const auto& pseudopotential_set = basis.pseudopotential_set();
1112 const auto& ultrasoft_basis =
1113 static_cast<const UltrasoftPlanewaveKohnShamQuantumSystem&>(scf_state_.quantum_system())
1114 .ultrasoft_basis();
1115
1116 int total_projectors = basis.total_number_of_projectors();
1117 if (total_projectors == 0)
1118 return hubbard_force;
1119
1120 constexpr Real q_spacing = 0.01;
1121 constexpr int number_of_q_points = 900;
1122 auto atomic_wavefunction_table = tabulate_atomic_wavefunctions_in_q_space(
1123 pseudopotential_set, crystal_structure_.cell_volume(), q_spacing, number_of_q_points);
1124
1125 size_t number_of_manifolds = manifolds.size();
1126 std::vector<int> manifold_row_offsets(number_of_manifolds);
1127 int total_hubbard_orbitals = 0;
1128 for (size_t manifold_index = 0; manifold_index < number_of_manifolds; manifold_index++) {
1129 manifold_row_offsets[manifold_index] = total_hubbard_orbitals;
1130 total_hubbard_orbitals +=
1131 static_cast<int>(manifolds[manifold_index].orbital_indices.size());
1132 }
1133
1134 /**
1135 * Which atom each projector belongs to, and the q_ij block of each.
1136 * q_ij is diagonal in the combined lm index, the same rule
1137 * `augment_hubbard_projectors_with_overlap` applies when it builds
1138 * S|phi> during the SCF — the two must agree or the force differentiates
1139 * a different S than the energy used.
1140 */
1141 std::vector<int> projector_atom(total_projectors, 0);
1142 std::vector<int> projector_offset_by_atom(number_of_atoms, 0);
1143 std::vector<int> channels_by_atom(number_of_atoms, 0);
1144 {
1145 int offset = 0;
1146 for (int atom_index = 0; atom_index < number_of_atoms; atom_index++) {
1147 int species_index = crystal_structure_.species_type_index(atom_index);
1148 int number_of_channels = basis.number_of_projector_channels(species_index);
1149 projector_offset_by_atom[atom_index] = offset;
1150 channels_by_atom[atom_index] = number_of_channels;
1151 for (int channel = 0; channel < number_of_channels; channel++)
1152 projector_atom[offset + channel] = atom_index;
1153 offset += number_of_channels;
1154 }
1155 }
1156
1157 auto augmentation_charge = [&](int atom_index, int channel_i, int channel_j) -> Real {
1158 return hubbard_augmentation_overlap(ultrasoft_basis,
1159 basis,
1160 crystal_structure_.species_type_index(atom_index),
1161 channel_i,
1162 channel_j);
1163 };
1164
1165 int number_of_kpoints = kpoint_grid_.number_of_kpoints();
1166 int number_of_bands = scf_state_.number_of_bands();
1167 int number_of_kpoint_states = number_of_kpoints * scf_state_.number_of_spin_components();
1168 const Matrix3x3& reciprocal_lattice_vectors =
1169 crystal_structure_.reciprocal_lattice_vectors();
1170
1171 for (int kpoint_index_state = 0; kpoint_index_state < number_of_kpoint_states;
1172 kpoint_index_state++) {
1173 int spin_index = kpoint_index_state / number_of_kpoints;
1174 int kpoint_index = kpoint_index_state % number_of_kpoints;
1175 int number_of_plane_waves = scf_state_.number_of_plane_waves(kpoint_index);
1176 const auto& plane_wave_indices = scf_state_.g_vector_indices_at_kpoint(kpoint_index);
1177 Real kpoint_weight = kpoint_grid_[kpoint_index].weight;
1178 bool is_gamma_only = planewave_basis_.is_gamma_only();
1179 int g_zero_index = is_gamma_only ? find_g_zero_plane_wave_index(planewave_basis_,
1180 plane_wave_indices,
1181 number_of_plane_waves)
1182 : -1;
1183
1184 Vector3D kpoint_cartesian;
1185 {
1186 Vector3D kpoint_vector = kpoint_grid_[kpoint_index].coordinates;
1187 for (int direction = 0; direction < 3; direction++)
1188 kpoint_cartesian[direction] =
1189 kpoint_vector.x * reciprocal_lattice_vectors[0][direction] +
1190 kpoint_vector.y * reciprocal_lattice_vectors[1][direction] +
1191 kpoint_vector.z * reciprocal_lattice_vectors[2][direction];
1192 }
1193
1194 std::vector<std::vector<Complex>> all_atomic_orbitals;
1195 build_atomic_orbitals_g_space(kpoint_index,
1196 number_of_plane_waves,
1197 plane_wave_indices,
1198 crystal_structure_,
1199 kpoint_grid_,
1200 planewave_basis_,
1201 pseudopotential_set,
1202 atomic_wavefunction_table,
1203 q_spacing,
1204 number_of_q_points,
1205 all_atomic_orbitals);
1206
1207 Matrix<Complex> beta_projectors;
1208 compute_kleinman_bylander_projectors(kpoint_grid_[kpoint_index].coordinates,
1209 crystal_structure_,
1210 planewave_basis_,
1211 pseudopotential_set,
1212 basis.beta_function_table(),
1213 plane_wave_indices,
1214 number_of_plane_waves,
1215 total_projectors,
1216 beta_projectors);
1217
1218 std::vector<std::vector<Complex>> minus_i_k_plus_g(
1219 3, std::vector<Complex>(number_of_plane_waves));
1220 for (int g_index = 0; g_index < number_of_plane_waves; g_index++) {
1221 const auto& g_vector = planewave_basis_.g_vectors()[plane_wave_indices[g_index]];
1222 for (int direction = 0; direction < 3; direction++)
1223 minus_i_k_plus_g[direction][g_index] =
1224 Complex(0.0,
1225 -(g_vector.g_cart[direction] + kpoint_cartesian[direction]) *
1226 crystal_structure_.two_pi_over_lattice_parameter());
1227 }
1228
1229 std::vector<Real> band_occupations(number_of_bands);
1230 for (int band_index = 0; band_index < number_of_bands; band_index++)
1231 band_occupations[band_index] =
1232 scf_state_.occupation(kpoint_index_state, band_index);
1233
1234 /**
1235 * Every overlap this assembly needs, formed once per k-point.
1236 * Indexed [row * columns + column] throughout.
1237 *
1238 * Each is a partial sum under band (G-vector) parallelism and is
1239 * reduced before use — the same reason the occupation matrix and
1240 * S|phi> are reduced at their own G-sums rather than afterwards.
1241 *
1242 * Gamma-only sums half the G sphere, and every operand here is real
1243 * in real space, so each contraction is completed to its full-sphere
1244 * (and real) value on the spot. `modifier` decides the G = 0 rule:
1245 * a plain overlap has G = 0 in the half sphere and the doubling
1246 * counts it twice, while a modified one carries a factor of
1247 * (k+G)_alpha = G_alpha that vanishes there. Every consumer below
1248 * multiplies these together, so completing them individually is what
1249 * makes the products right.
1250 */
1251 auto contract = [&](const Complex* left_rows,
1252 int number_of_left,
1253 size_t left_stride,
1254 const Complex* right_columns,
1255 int number_of_right,
1256 size_t right_stride,
1257 const Complex* modifier) {
1258 std::vector<Complex> result(static_cast<size_t>(number_of_left) * number_of_right,
1259 Complex(0.0, 0.0));
1260 for (int left = 0; left < number_of_left; left++) {
1261 const Complex* left_row = left_rows + static_cast<size_t>(left) * left_stride;
1262 for (int right = 0; right < number_of_right; right++) {
1263 const Complex* right_column =
1264 right_columns + static_cast<size_t>(right) * right_stride;
1265 Complex accumulated(0.0, 0.0);
1266 for (int g_index = 0; g_index < number_of_plane_waves; g_index++) {
1267 Complex value = std::conj(left_row[g_index]) * right_column[g_index];
1268 accumulated += modifier ? value * modifier[g_index] : value;
1269 }
1270 if (is_gamma_only) {
1271 Real g_zero_product = (modifier == nullptr && g_zero_index >= 0)
1272 ? (std::conj(left_row[g_zero_index]) *
1273 right_column[g_zero_index])
1274 .real()
1275 : 0.0;
1276 accumulated =
1277 gamma_only_full_sphere_overlap(accumulated, g_zero_product);
1278 }
1279 result[static_cast<size_t>(left) * number_of_right + right] = accumulated;
1280 }
1281 }
1282 if (scf_state_.execution_policy != nullptr)
1283 scf_state_.execution_policy->reduce_within_group(
1284 result.data(), static_cast<int>(result.size()));
1285 return result;
1286 };
1287
1288 // Hubbard orbitals and wavefunctions packed contiguously per row/column.
1289 std::vector<Complex> orbitals(static_cast<size_t>(total_hubbard_orbitals) *
1290 number_of_plane_waves);
1291 for (size_t manifold_index = 0; manifold_index < number_of_manifolds;
1292 manifold_index++) {
1293 const auto& manifold = manifolds[manifold_index];
1294 int row_offset = manifold_row_offsets[manifold_index];
1295 for (size_t m = 0; m < manifold.orbital_indices.size(); m++) {
1296 const auto& orbital =
1297 all_atomic_orbitals[static_cast<size_t>(manifold.orbital_indices[m])];
1298 std::copy(orbital.begin(),
1299 orbital.begin() + number_of_plane_waves,
1300 orbitals.begin() +
1301 (static_cast<size_t>(row_offset) + m) * number_of_plane_waves);
1302 }
1303 }
1304 std::vector<Complex> wavefunctions(static_cast<size_t>(number_of_bands) *
1305 number_of_plane_waves);
1306 for (int band_index = 0; band_index < number_of_bands; band_index++)
1307 for (int g_index = 0; g_index < number_of_plane_waves; g_index++)
1308 wavefunctions[static_cast<size_t>(band_index) * number_of_plane_waves +
1309 g_index] = scf_state_.wavefunction_coefficient(kpoint_index_state,
1310 band_index,
1311 g_index);
1312
1313 auto phi_psi = contract(orbitals.data(),
1314 total_hubbard_orbitals,
1315 number_of_plane_waves,
1316 wavefunctions.data(),
1317 number_of_bands,
1318 number_of_plane_waves,
1319 nullptr);
1320 auto phi_beta = contract(orbitals.data(),
1321 total_hubbard_orbitals,
1322 number_of_plane_waves,
1323 beta_projectors.data(),
1324 total_projectors,
1325 number_of_plane_waves,
1326 nullptr);
1327 auto beta_psi = contract(beta_projectors.data(),
1328 total_projectors,
1329 number_of_plane_waves,
1330 wavefunctions.data(),
1331 number_of_bands,
1332 number_of_plane_waves,
1333 nullptr);
1334 std::vector<std::vector<Complex>> phi_derivative_psi(3), phi_derivative_beta(3),
1335 beta_derivative_psi(3);
1336 for (int direction = 0; direction < 3; direction++) {
1337 phi_derivative_psi[direction] = contract(orbitals.data(),
1338 total_hubbard_orbitals,
1339 number_of_plane_waves,
1340 wavefunctions.data(),
1341 number_of_bands,
1342 number_of_plane_waves,
1343 minus_i_k_plus_g[direction].data());
1344 phi_derivative_beta[direction] = contract(orbitals.data(),
1345 total_hubbard_orbitals,
1346 number_of_plane_waves,
1347 beta_projectors.data(),
1348 total_projectors,
1349 number_of_plane_waves,
1350 minus_i_k_plus_g[direction].data());
1351 beta_derivative_psi[direction] = contract(beta_projectors.data(),
1352 total_projectors,
1353 number_of_plane_waves,
1354 wavefunctions.data(),
1355 number_of_bands,
1356 number_of_plane_waves,
1357 minus_i_k_plus_g[direction].data());
1358 }
1359
1360 auto orbital_band = [&](int orbital, int band) {
1361 return static_cast<size_t>(orbital) * number_of_bands + band;
1362 };
1363 auto orbital_projector = [&](int orbital, int projector) {
1364 return static_cast<size_t>(orbital) * total_projectors + projector;
1365 };
1366 auto projector_band = [&](int projector, int band) {
1367 return static_cast<size_t>(projector) * number_of_bands + band;
1368 };
1369
1370 /**
1371 * P_m = <phi_m|S|psi>, and the delta_a bracket -- the part of
1372 * -dP/dtau that appears only on the manifold's own atom. Both sum
1373 * q_ij over *every* atom's projectors, because S does.
1374 */
1375 std::vector<Complex> projection = phi_psi;
1376 std::vector<std::vector<Complex>> own_atom_term(3, phi_derivative_psi[0]);
1377 for (int direction = 0; direction < 3; direction++)
1378 own_atom_term[direction] = phi_derivative_psi[direction];
1379
1380 for (int atom_index = 0; atom_index < number_of_atoms; atom_index++) {
1381 int offset = projector_offset_by_atom[atom_index];
1382 for (int channel_i = 0; channel_i < channels_by_atom[atom_index]; channel_i++) {
1383 for (int channel_j = 0; channel_j < channels_by_atom[atom_index]; channel_j++) {
1384 Real q_value = augmentation_charge(atom_index, channel_i, channel_j);
1385 if (std::abs(q_value) < math::tolerances::d_matrix_element_threshold)
1386 continue;
1387 int projector_i = offset + channel_i;
1388 int projector_j = offset + channel_j;
1389 for (int orbital = 0; orbital < total_hubbard_orbitals; orbital++) {
1390 Complex phi_beta_i = phi_beta[orbital_projector(orbital, projector_i)];
1391 for (int band = 0; band < number_of_bands; band++) {
1392 Complex beta_j_psi = beta_psi[projector_band(projector_j, band)];
1393 projection[orbital_band(orbital, band)] +=
1394 q_value * phi_beta_i * beta_j_psi;
1395 for (int direction = 0; direction < 3; direction++)
1396 own_atom_term[direction][orbital_band(orbital, band)] +=
1397 q_value *
1398 phi_derivative_beta[direction][orbital_projector(
1399 orbital, projector_i)] *
1400 beta_j_psi;
1401 }
1402 }
1403 }
1404 }
1405 }
1406
1407 Matrix<Complex> overlaps(static_cast<size_t>(number_of_bands),
1408 static_cast<size_t>(total_hubbard_orbitals));
1409 for (int orbital = 0; orbital < total_hubbard_orbitals; orbital++)
1410 for (int band = 0; band < number_of_bands; band++)
1411 overlaps(static_cast<size_t>(band), static_cast<size_t>(orbital)) =
1412 projection[orbital_band(orbital, band)];
1413
1414 /**
1415 * One displaced atom at a time. The dS/dtau part is non-zero for
1416 * every atom carrying projectors, which is exactly the term the
1417 * norm-conserving path has no analogue of.
1418 */
1419 for (int displaced_atom = 0; displaced_atom < number_of_atoms; displaced_atom++) {
1420 int offset = projector_offset_by_atom[displaced_atom];
1421 int number_of_channels = channels_by_atom[displaced_atom];
1422
1423 for (int direction = 0; direction < 3; direction++) {
1424 std::vector<Complex> derivative(static_cast<size_t>(total_hubbard_orbitals) *
1425 number_of_bands,
1426 Complex(0.0, 0.0));
1427
1428 for (int channel_i = 0; channel_i < number_of_channels; channel_i++) {
1429 for (int channel_j = 0; channel_j < number_of_channels; channel_j++) {
1430 Real q_value =
1431 augmentation_charge(displaced_atom, channel_i, channel_j);
1432 if (std::abs(q_value) < math::tolerances::d_matrix_element_threshold)
1433 continue;
1434 int projector_i = offset + channel_i;
1435 int projector_j = offset + channel_j;
1436 for (int orbital = 0; orbital < total_hubbard_orbitals; orbital++) {
1437 Complex phi_beta_i =
1438 phi_beta[orbital_projector(orbital, projector_i)];
1439 Complex phi_derivative_beta_i =
1440 phi_derivative_beta[direction]
1441 [orbital_projector(orbital, projector_i)];
1442 for (int band = 0; band < number_of_bands; band++) {
1443 Complex beta_j_psi =
1444 beta_psi[projector_band(projector_j, band)];
1445 Complex beta_j_derivative_psi =
1446 beta_derivative_psi[direction]
1447 [projector_band(projector_j, band)];
1448 derivative[orbital_band(orbital, band)] -=
1449 q_value * (phi_derivative_beta_i * beta_j_psi -
1450 phi_beta_i * beta_j_derivative_psi);
1451 }
1452 }
1453 }
1454 }
1455
1456 for (size_t manifold_index = 0; manifold_index < number_of_manifolds;
1457 manifold_index++) {
1458 const auto& manifold = manifolds[manifold_index];
1459 int manifold_size = static_cast<int>(manifold.orbital_indices.size());
1460 int row_offset = manifold_row_offsets[manifold_index];
1461 bool displaced_owns_manifold = manifold.atom_index == displaced_atom;
1462
1463 Matrix<Complex> derivative_overlaps(static_cast<size_t>(number_of_bands),
1464 static_cast<size_t>(manifold_size));
1465 for (int m = 0; m < manifold_size; m++) {
1466 int orbital = row_offset + m;
1467 for (int band = 0; band < number_of_bands; band++) {
1468 Complex value = derivative[orbital_band(orbital, band)];
1469 if (displaced_owns_manifold)
1470 value += own_atom_term[direction][orbital_band(orbital, band)];
1471 derivative_overlaps(static_cast<size_t>(band),
1472 static_cast<size_t>(m)) = value;
1473 }
1474 }
1475
1476 const auto& combined_d_matrix =
1477 scf_state_.hubbard_d_matrix_per_spin[spin_index];
1478 std::vector<Complex> manifold_d_matrix(static_cast<size_t>(manifold_size) *
1479 manifold_size);
1480 for (int m = 0; m < manifold_size; m++)
1481 for (int m_prime = 0; m_prime < manifold_size; m_prime++)
1482 manifold_d_matrix[static_cast<size_t>(m * manifold_size +
1483 m_prime)] =
1484 combined_d_matrix[static_cast<size_t>(row_offset + m) +
1485 static_cast<size_t>(row_offset + m_prime) *
1486 total_hubbard_orbitals];
1487
1488 Matrix<Complex> manifold_overlaps(static_cast<size_t>(number_of_bands),
1489 static_cast<size_t>(manifold_size));
1490 for (int m = 0; m < manifold_size; m++)
1491 for (int band = 0; band < number_of_bands; band++)
1492 manifold_overlaps(static_cast<size_t>(band),
1493 static_cast<size_t>(m)) =
1494 overlaps(static_cast<size_t>(band),
1495 static_cast<size_t>(row_offset + m));
1496
1497 hubbard_force[displaced_atom][direction] +=
1498 accumulate_hubbard_force_gradient(number_of_bands,
1499 manifold_size,
1500 manifold_overlaps,
1501 derivative_overlaps,
1502 band_occupations,
1503 kpoint_weight,
1504 manifold_d_matrix);
1505 }
1506 }
1507 }
1508 }
1509
1510 return hubbard_force;
1511 }
1512
1513 const Crystal& crystal_structure_;
1514 const PlanewaveBasis& planewave_basis_;
1515 const KPointsGrid& kpoint_grid_;
1516 const NormConservingPlanewaveKohnShamSCFState& scf_state_;
1517};
1518
1519} // namespace q3
atomic_orbital_builder.hpp
q3::Crystal::cell_volume
Real cell_volume() const
Definition crystal.hpp:193
q3::Crystal::two_pi_over_lattice_parameter
Real two_pi_over_lattice_parameter() const
Definition crystal.hpp:196
q3::Crystal::reciprocal_lattice_vectors
const Matrix3x3 & reciprocal_lattice_vectors() const
Definition crystal.hpp:205
q3::Crystal::species_type_index
int species_type_index(int atom_index) const
Definition crystal.hpp:273
q3::Crystal::number_of_atoms
int number_of_atoms() const
Definition crystal.hpp:209
q3::HubbardForceCalculator::HubbardForceCalculator
HubbardForceCalculator(const Crystal &crystal_structure, const PlanewaveBasis &planewave_basis, const KPointsGrid &kpoint_grid, const NormConservingPlanewaveKohnShamSCFState &scf_state)
Definition hubbard.hpp:45
q3::HubbardForceCalculator::compute
std::vector< Vector3D > compute() const
Definition hubbard.hpp:54
q3::KPointsGrid::number_of_kpoints
int number_of_kpoints() const
Definition kpoints_grid.hpp:71
q3::Matrix::data
T * data()
Returns the underlying flat data pointer.
Definition matrix.hpp:75
q3::NormConservingPlanewaveBasis::pseudopotential_set
const PseudopotentialSet & pseudopotential_set() const
Definition init.hpp:156
q3::NormConservingPlanewaveKohnShamQuantumSystem::pseudopotential_basis
const NormConservingPlanewaveBasis & pseudopotential_basis() const
Definition init.hpp:65
q3::PlanewaveBasis::g_vectors
const std::vector< GVector > & g_vectors() const
Definition init.hpp:148
q3::PlanewaveBasis::is_gamma_only
bool is_gamma_only() const
Definition init.hpp:157
q3::execution::ExecutionPolicy::reduce_within_group
virtual void reduce_within_group(ComplexValue *data, int number_of_elements) const =0
kpoints_grid.hpp
init.hpp
init.hpp
init.hpp
init.hpp
factory.hpp
hubbard.hpp
init.hpp
lowdin_orthogonalization.hpp
q3::math::tolerances::d_matrix_element_threshold
constexpr double d_matrix_element_threshold
Definition tolerances.hpp:99
q3::math::tolerances::occupation_threshold
constexpr double occupation_threshold
Definition tolerances.hpp:119
q3::hubbard_augmentation_overlap
Real hubbard_augmentation_overlap(const UltrasoftPlanewaveBasis &ultrasoft_basis, const NormConservingPlanewaveBasis &pseudopotential_basis, int species_index, int channel_row, int channel_column)
Definition hubbard.hpp:258
q3::build_atomic_orbitals_g_space
void build_atomic_orbitals_g_space(int kpoint_index, int number_of_plane_waves, const std::vector< int > &g_vector_indices, const Crystal &crystal_structure, const KPointsGrid &kpoint_grid, const PlanewaveBasis &planewave_basis, const PseudopotentialSet &pseudopotential_set, const std::vector< std::vector< std::vector< Real > > > &atomic_wavefunction_table, Real q_spacing, int number_of_q_points, std::vector< std::vector< Complex > > &orbitals)
Definition atomic_orbital_builder.hpp:125
q3::compute_hubbard_projector_overlaps
void compute_hubbard_projector_overlaps(int number_of_plane_waves, int total_hubbard_orbitals, int number_of_bands, const Complex *packed_orbitals, const Complex *band_functions, std::vector< Complex > &overlaps)
Definition hubbard.hpp:304
q3::Real
double Real
Definition type_aliases.hpp:20
q3::convert_gamma_only_projector_overlaps_to_full_sphere
void convert_gamma_only_projector_overlaps_to_full_sphere(int number_of_plane_waves, int total_hubbard_orbitals, int number_of_bands, const Complex *packed_orbitals, const Complex *band_functions, int g_zero_index, bool band_functions_carry_g_component, std::vector< Complex > &overlaps)
Definition hubbard.hpp:392
q3::accumulate_hubbard_force_gradient
Real accumulate_hubbard_force_gradient(int number_of_bands, int manifold_size, const Matrix< Complex > &overlaps, const Matrix< Complex > &direction_derivative_overlaps, const std::vector< Real > &band_occupations, Real kpoint_weight, const std::vector< Complex > &d_matrix)
Definition hubbard.hpp:489
q3::enumerate_atomic_orbitals
std::vector< AtomicOrbitalLabel > enumerate_atomic_orbitals(const Crystal &crystal_structure, const PseudopotentialSet &pseudopotential_set)
Definition atomic_orbital_builder.hpp:91
q3::create_overlap_operator
std::unique_ptr< OverlapOperator > create_overlap_operator(const Crystal &crystal_structure, const NormConservingPlanewaveKohnShamQuantumSystem &quantum_system, const Matrix< Complex > &kleinman_bylander_potential, const execution::ExecutionPolicy &execution_policy)
Definition factory.hpp:23
q3::compute_kleinman_bylander_projectors
void compute_kleinman_bylander_projectors(const Vector3D &k_crystal, const Crystal &crystal_structure, const PlanewaveBasis &planewave_basis, const PseudopotentialSet &pseudopotential_set, const BetaFunctionTable &beta_function_table, const std::vector< int > &g_vector_indices, int number_of_plane_waves, int total_number_of_projectors, Matrix< Complex > &kleinman_bylander_potential)
Definition projector.hpp:36
q3::tabulate_atomic_wavefunctions_in_q_space
std::vector< std::vector< std::vector< Real > > > tabulate_atomic_wavefunctions_in_q_space(const PseudopotentialSet &pseudopotential_set, Real cell_volume, Real q_spacing, int number_of_q_points)
Definition wavefunction_atomic.hpp:50
q3::find_g_zero_plane_wave_index
int find_g_zero_plane_wave_index(const PlanewaveBasis &planewave_basis, const std::vector< int > &plane_wave_indices, int number_of_plane_waves)
Definition hubbard.hpp:364
q3::lowdin_inverse_square_root_derivative
std::vector< Complex > lowdin_inverse_square_root_derivative(const LowdinTransform &transform, const std::vector< Complex > &overlap_derivative)
Definition lowdin_orthogonalization.hpp:221
q3::gamma_only_full_sphere_overlap
Complex gamma_only_full_sphere_overlap(Complex half_sphere_sum, Real g_zero_product)
Definition hubbard.hpp:352
q3::Complex
std::complex< double > Complex
Definition type_aliases.hpp:21
q3::compute_lowdin_transform
LowdinTransform compute_lowdin_transform(std::vector< Complex > overlap_column_major, int number_of_orbitals)
Definition lowdin_orthogonalization.hpp:68
wavefunction_atomic.hpp
q3::Matrix3x3
3×3 matrix stored row-major: row i = mat[i]
Definition matrix_3x3.hpp:15
q3::NormConservingPlanewaveKohnShamSCFState::quantum_system
const NormConservingPlanewaveKohnShamQuantumSystem & quantum_system() const override
Definition init.hpp:82
q3::NormConservingPlanewaveKohnShamSCFState::hubbard_manifolds
std::vector< HubbardManifold > hubbard_manifolds
DFT+U (Hubbard) resolved manifolds at convergence; empty means no correction is active.
Definition init.hpp:62
q3::NormConservingPlanewaveKohnShamSCFState::hubbard_d_matrix_per_spin
std::vector< std::vector< Complex > > hubbard_d_matrix_per_spin
Definition init.hpp:67
q3::PlanewaveKohnShamSCFState::g_vector_indices_at_kpoint
const std::vector< int > & g_vector_indices_at_kpoint(int kpoint_index) const
Definition init.hpp:80
q3::PlanewaveKohnShamSCFState::number_of_plane_waves
int number_of_plane_waves(int kpoint_index) const
Definition init.hpp:77
q3::PlanewaveKohnShamSCFState::occupation
Real occupation(int kpoint_index, int band_index) const
Definition init.hpp:90
q3::PlanewaveKohnShamSCFState::execution_policy
const execution::PlanewaveExecutionPolicy * execution_policy
Definition init.hpp:64
q3::PlanewaveKohnShamSCFState::number_of_spin_components
int number_of_spin_components() const
Definition init.hpp:106
q3::PlanewaveKohnShamSCFState::wavefunction_coefficient
Complex wavefunction_coefficient(int kpoint_index, int band_index, int g_vector_index) const
Definition init.hpp:84