Loading...
Searching...
No Matches
hubbard.hpp
Go to the documentation of this file.
1#pragma once
2
3/**
4 * DFT+U (Hubbard) stress tensor contribution, in the rotationally invariant simplified
5 * form of Dudarev et al., Phys. Rev. B 57, 1505 (1998).
6 *
7 * The strain analogue of HubbardForceCalculator
8 * (forces/planewave/pseudopotential/internal/hubbard.hpp): the same D-matrix-weighted
9 * contraction of dn/d(parameter), with the parameter being the strain rather than an
10 * atomic position. E_U = (U/2)Tr[n(1-n)] depends on strain only through the projection
11 *
12 * P_m = <phi_m|S|psi_v>, S = 1 + sum_ij q_ij |beta_i><beta_j|
13 *
14 * with psi and the occupations held fixed, so
15 *
16 * dP_m/deps_ab = <dphi_m/deps_ab|S|psi> + <phi_m|dS/deps_ab|psi>
17 * dS/deps_ab = sum_ij q_ij ( |dbeta_i/deps_ab><beta_j| + |beta_i><dbeta_j/deps_ab| )
18 *
19 * q_ij = int Q_ij(r)dr is a pure number and carries no strain derivative. Both phi and
20 * beta have the form A(k+G) = f(|k+G|) Y_lm(k+G) S_I(k+G) i^l / sqrt(Omega), whose strain
21 * derivative is
22 *
23 * dA/deps_ab = -[ f'(q) (k+G)_a (k+G)_b / q Y_lm + f(q) dY_lm/dq_a (k+G)_b ]
24 * - delta_ab A / 2
25 *
26 * the last term from 1/sqrt(Omega) (dOmega/deps_ab = Omega delta_ab), and the structure
27 * factor contributing nothing because (k+G).tau is strain-invariant. Quantum ESPRESSO's
28 * Hubbard stress builds the same object from the same three pieces: an angular term
29 * weighted by -(k+G)_b, a radial term weighted by -(k+G)_a (k+G)_b/|k+G|, and -1/2 times
30 * the orbital itself on the diagonal (a == b).
31 *
32 * ── Sign and normalization ────────────────────────────────────────────────
33 *
34 * Pinned by reduction to the already-validated nonlocal stress rather than re-derived
35 * from the energy. NonlocalPseudopotentialStressTensorCalculator (internal/nonlocal.hpp)
36 * accumulates, for E_nl = sum_kv w sum_ij D_ij <psi|beta_i><beta_j|psi>,
37 *
38 * sigma_ab = (1/Omega) sum D_ij [ 2Re(conj(becp_i)(bessel_j + ylm_j))
39 * + delta_ab Re(conj(becp_i) becp_j) ]
40 *
41 * which is identically -(1/Omega) dE_nl/deps_ab once dbeta/deps is substituted in the
42 * form above -- the factor 2 on the derivative overlaps and the factor 1 on the diagonal
43 * are what the -[bessel + ylm] - delta_ab/2 split produces. So Q3's convention is
44 * sigma = -(1/Omega) dE/deps, the same sign and normalization Quantum ESPRESSO reports
45 * for its own Hubbard stress tensor, and this calculator writes the Hubbard term the
46 * same way: form dP/deps, hand it to the
47 * shared `accumulate_hubbard_force_gradient` -- which returns
48 * sum_v w f 2Re[sum_mm' D_mm' conj(X_m) P_m'], i.e. dE_U/dlambda when X = dP/dlambda --
49 * and multiply by -1/Omega. Setting q_ij = 0 collapses dP to <dphi/deps|psi> and the
50 * whole expression to the norm-conserving limit, where dS/deps vanishes identically.
51 *
52 * ── Time reversal ─────────────────────────────────────────────────────────
53 *
54 * `accumulate_hubbard_force_gradient` contracts 2*Re[...], which *is* the time-reversal
55 * average of dn, matching how the occupation matrix itself is built
56 * (symmetrize_hubbard_occupation_matrices_under_time_reversal). Using it rather than a
57 * private contraction is what keeps dn/deps consistent with n, the inconsistency that
58 * produced a 0.3% force error on the symmetry-reduced wedge (KNOWN-ISSUES.md,
59 * 2026-08-09).
60 *
61 * The full 3x3 strain tensor is accumulated, not just its upper triangle. Of the three
62 * pieces above only the dY_lm/dq_a (k+G)_b term is asymmetric in (a,b), and summed over
63 * a symmetry-reduced k-set that asymmetry survives, so sigma_ab and sigma_ba are
64 * genuinely different numbers whose average is the stress. `symmetrize_stress_tensor`
65 * therefore averages rather than mirrors, cancelling the antisymmetric half instead of
66 * substituting one triangle for it -- unlike the nonlocal calculator's same-named
67 * helper, which fills the upper triangle only and mirrors it down.
68 */
69
70#include <algorithm>
71#include <array>
72#include <cmath>
73#include <vector>
74
75#include "layer_0_core/math/interpolation/lagrange.hpp"
76#include "layer_0_core/math/phase_factor.hpp"
77#include "layer_0_core/math/spherical/harmonics.hpp"
78#include "layer_0_core/math/tolerances.hpp"
79#include "layer_0_core/types/init.hpp"
80#include "layer_0_core/types/matrix.hpp"
81#include "layer_1_structure/crystal/crystal.hpp"
82#include "layer_1_structure/reciprocal_space/kpoints_grid.hpp"
83#include "layer_3_representation/basis/planewave/init.hpp"
84#include "layer_3_representation/pseudopotential/planewave/internal/atomic_orbital_builder.hpp"
85#include "layer_3_representation/pseudopotential/planewave/internal/lowdin_orthogonalization.hpp"
86#include "layer_3_representation/pseudopotential/planewave/internal/projector.hpp"
87#include "layer_3_representation/pseudopotential/planewave/internal/wavefunction_atomic.hpp"
88#include "layer_3_representation/quantum_system/kohn_sham/planewave/ultrasoft/init.hpp"
89#include "layer_4_computation/diagonalizer/planewave/internal/overlap/pseudopotential/factory.hpp"
90#include "layer_4_computation/effective_potential/internal/hubbard/hubbard.hpp"
91#include "layer_4_computation/state/kohn_sham/planewave/norm_conserving/init.hpp"
92#include "layer_5_properties/stress_tensor/planewave/contributions.hpp"
93
94namespace q3 {
95
96class HubbardStressTensorCalculator {
97 public:
98 HubbardStressTensorCalculator(const Crystal& crystal_structure,
99 const PlanewaveBasis& planewave_basis,
100 const KPointsGrid& kpoint_grid,
101 const NormConservingPlanewaveKohnShamSCFState& scf_state,
102 StressTensorContributions& result)
103 : crystal_structure_(crystal_structure)
104 , planewave_basis_(planewave_basis)
105 , kpoint_grid_(kpoint_grid)
106 , scf_state_(scf_state)
107 , result_(result) {}
108
109 /**
110 * Whether a Hubbard stress term exists for this run's projector and dataset choice.
111 *
112 * Ortho-atomic projectors are a joint, system-wide Lowdin transform of every atomic
113 * orbital in the cell, so their strain derivative carries a d(O^-1/2)/deps term over
114 * orbitals belonging to no Hubbard manifold at all. The norm-conserving path evaluates
115 * that term with `lowdin_inverse_square_root_derivative`. Ultrasoft and PAW use the same
116 * matrix-function derivative with the generalized <phi|S|phi> metric and its complete
117 * orbital-plus-overlap strain derivative.
118 */
119 [[nodiscard]] static bool supports(const NormConservingPlanewaveKohnShamSCFState& scf_state) {
120 return scf_state.hubbard_manifolds.empty() ||
121 scf_state.hubbard_projector_type == "atomic" ||
122 scf_state.hubbard_projector_type == "ortho-atomic";
123 }
124
125 void compute() {
126 const auto& manifolds = scf_state_.hubbard_manifolds;
127 if (manifolds.empty() || !supports(scf_state_))
128 return;
129 if (scf_state_.hubbard_projector_type == "ortho-atomic") {
130 compute_with_orthogonal_projectors();
131 return;
132 }
133
134 const auto& basis = scf_state_.quantum_system().pseudopotential_basis();
135 const auto& pseudopotential_set = basis.pseudopotential_set();
136
137 /**
138 * dS/deps is identically zero without augmentation charges, so the whole
139 * beta-projector half of the assembly is switched off by giving it zero
140 * projectors rather than by branching through it. `ultrasoft_basis()` is only
141 * reachable once has_ultrasoft() is true.
142 */
143 bool has_overlap_operator = scf_state_.quantum_system().has_ultrasoft();
144 int total_projectors = has_overlap_operator ? basis.total_number_of_projectors() : 0;
145
146 constexpr Real q_spacing = 0.01;
147 constexpr int number_of_q_points = 900;
148 auto atomic_wavefunction_table = tabulate_atomic_wavefunctions_in_q_space(
149 pseudopotential_set, crystal_structure_.cell_volume(), q_spacing, number_of_q_points);
150
151 auto orbital_descriptors = enumerate_atomic_orbital_descriptors();
152
153 size_t number_of_manifolds = manifolds.size();
154 std::vector<int> manifold_row_offsets(number_of_manifolds);
155 int total_hubbard_orbitals = 0;
156 for (size_t manifold_index = 0; manifold_index < number_of_manifolds; manifold_index++) {
157 manifold_row_offsets[manifold_index] = total_hubbard_orbitals;
158 total_hubbard_orbitals +=
159 static_cast<int>(manifolds[manifold_index].orbital_indices.size());
160 }
161 if (total_hubbard_orbitals == 0)
162 return;
163
164 // The descriptor of each packed Hubbard orbital, in the packed order.
165 std::vector<AtomicOrbitalDescriptor> hubbard_descriptors(
166 static_cast<size_t>(total_hubbard_orbitals));
167 for (size_t manifold_index = 0; manifold_index < number_of_manifolds; manifold_index++) {
168 const auto& manifold = manifolds[manifold_index];
169 int row_offset = manifold_row_offsets[manifold_index];
170 for (size_t m = 0; m < manifold.orbital_indices.size(); m++)
171 hubbard_descriptors[static_cast<size_t>(row_offset) + m] =
172 orbital_descriptors[static_cast<size_t>(manifold.orbital_indices[m])];
173 }
174
175 int number_of_atoms = crystal_structure_.number_of_atoms();
176 std::vector<int> projector_offset_by_atom(number_of_atoms, 0);
177 std::vector<int> channels_by_atom(number_of_atoms, 0);
178 // Left all-zero for norm-conserving, which makes every q_ij loop below empty.
179 if (total_projectors > 0) {
180 int offset = 0;
181 for (int atom_index = 0; atom_index < number_of_atoms; atom_index++) {
182 int species_index = crystal_structure_.species_type_index(atom_index);
183 int number_of_channels = basis.number_of_projector_channels(species_index);
184 projector_offset_by_atom[atom_index] = offset;
185 channels_by_atom[atom_index] = number_of_channels;
186 offset += number_of_channels;
187 }
188 }
189
190 auto augmentation_charge = [&](int atom_index, int channel_i, int channel_j) -> Real {
191 return hubbard_augmentation_overlap(
192 static_cast<const UltrasoftPlanewaveKohnShamQuantumSystem&>(
193 scf_state_.quantum_system())
194 .ultrasoft_basis(),
195 basis,
196 crystal_structure_.species_type_index(atom_index),
197 channel_i,
198 channel_j);
199 };
200
201 int number_of_kpoints = kpoint_grid_.number_of_kpoints();
202 int number_of_bands = scf_state_.number_of_bands();
203 int number_of_kpoint_states = number_of_kpoints * scf_state_.number_of_spin_components();
204 const Matrix3x3& reciprocal_lattice_vectors =
205 crystal_structure_.reciprocal_lattice_vectors();
206 Real inverse_cell_volume = 1.0 / crystal_structure_.cell_volume();
207
208 for (int kpoint_index_state = 0; kpoint_index_state < number_of_kpoint_states;
209 kpoint_index_state++) {
210 int spin_index = kpoint_index_state / number_of_kpoints;
211 int kpoint_index = kpoint_index_state % number_of_kpoints;
212 int number_of_plane_waves = scf_state_.number_of_plane_waves(kpoint_index);
213 const auto& plane_wave_indices = scf_state_.g_vector_indices_at_kpoint(kpoint_index);
214 Real kpoint_weight = kpoint_grid_[kpoint_index].weight;
215 bool is_gamma_only = planewave_basis_.is_gamma_only();
216 int g_zero_index = is_gamma_only ? find_g_zero_plane_wave_index(planewave_basis_,
217 plane_wave_indices,
218 number_of_plane_waves)
219 : -1;
220
221 Vector3D kpoint_cartesian;
222 {
223 Vector3D kpoint_vector = kpoint_grid_[kpoint_index].coordinates;
224 for (int direction = 0; direction < 3; direction++)
225 kpoint_cartesian[direction] =
226 kpoint_vector.x * reciprocal_lattice_vectors[0][direction] +
227 kpoint_vector.y * reciprocal_lattice_vectors[1][direction] +
228 kpoint_vector.z * reciprocal_lattice_vectors[2][direction];
229 }
230
231 PlaneWaveGeometry geometry = build_plane_wave_geometry(
232 kpoint_cartesian, plane_wave_indices, number_of_plane_waves);
233
234 std::vector<Real> band_occupations(number_of_bands);
235 std::vector<int> occupied_band_indices;
236 for (int band_index = 0; band_index < number_of_bands; band_index++) {
237 band_occupations[band_index] =
238 scf_state_.occupation(kpoint_index_state, band_index);
239 if (band_occupations[band_index] >= math::tolerances::occupation_threshold)
240 occupied_band_indices.push_back(band_index);
241 }
242 int number_of_occupied_bands = static_cast<int>(occupied_band_indices.size());
243 if (number_of_occupied_bands == 0)
244 continue;
245
246 // Column-major (number_of_plane_waves x N) operands, the layout the ZGEMM wants.
247 std::vector<Complex> wavefunctions(static_cast<size_t>(number_of_plane_waves) *
248 number_of_occupied_bands);
249 for (int occupied_index = 0; occupied_index < number_of_occupied_bands;
250 occupied_index++)
251 for (int g_index = 0; g_index < number_of_plane_waves; g_index++)
252 wavefunctions[static_cast<size_t>(occupied_index) * number_of_plane_waves +
253 g_index] =
254 scf_state_.wavefunction_coefficient(
255 kpoint_index_state,
256 occupied_band_indices[static_cast<size_t>(occupied_index)],
257 g_index);
258
259 std::vector<Complex> orbitals;
260 std::vector<Complex> orbital_radial_derivative;
261 std::array<std::vector<Complex>, 3> orbital_angular_derivative;
262 build_atomic_orbital_strain_parts(kpoint_index,
263 plane_wave_indices,
264 number_of_plane_waves,
265 geometry,
266 hubbard_descriptors,
267 atomic_wavefunction_table,
268 q_spacing,
269 number_of_q_points,
270 orbitals,
271 orbital_radial_derivative,
272 orbital_angular_derivative);
273
274 std::vector<Complex> beta_projectors;
275 std::vector<Complex> beta_radial_derivative;
276 std::array<std::vector<Complex>, 3> beta_angular_derivative;
277 build_beta_strain_parts(kpoint_index,
278 plane_wave_indices,
279 number_of_plane_waves,
280 total_projectors,
281 geometry,
282 beta_projectors,
283 beta_radial_derivative,
284 beta_angular_derivative);
285
286 auto contract = [&](const Complex* left,
287 int number_of_left,
288 const Complex* right,
289 int number_of_right) {
290 std::vector<Complex> result;
291 compute_hubbard_projector_overlaps(
292 number_of_plane_waves, number_of_left, number_of_right, left, right, result);
293 if (is_gamma_only)
294 complete_gamma_only_overlaps(number_of_plane_waves,
295 number_of_left,
296 number_of_right,
297 left,
298 right,
299 g_zero_index,
300 result);
301 if (scf_state_.execution_policy != nullptr)
302 scf_state_.execution_policy->reduce_within_group(
303 result.data(), static_cast<int>(result.size()));
304 return result;
305 };
306
307 // Column-major throughout: element (left, right) lives at left + right*number_of_left.
308 auto orbital_band = [&](int orbital, int band) {
309 return static_cast<size_t>(orbital) +
310 static_cast<size_t>(band) * total_hubbard_orbitals;
311 };
312 auto orbital_projector = [&](int orbital, int projector) {
313 return static_cast<size_t>(orbital) +
314 static_cast<size_t>(projector) * total_hubbard_orbitals;
315 };
316 auto projector_band = [&](int projector, int band) {
317 return static_cast<size_t>(projector) +
318 static_cast<size_t>(band) * total_projectors;
319 };
320
321 auto phi_psi = contract(orbitals.data(),
322 total_hubbard_orbitals,
323 wavefunctions.data(),
324 number_of_occupied_bands);
325 auto phi_beta = contract(
326 orbitals.data(), total_hubbard_orbitals, beta_projectors.data(), total_projectors);
327 auto beta_psi = contract(beta_projectors.data(),
328 total_projectors,
329 wavefunctions.data(),
330 number_of_occupied_bands);
331
332 // P_m = <phi_m|S|psi_v>, exactly the projection the occupation matrix is built from.
333 std::vector<Complex> projection = phi_psi;
334 accumulate_augmentation(number_of_atoms,
335 channels_by_atom,
336 projector_offset_by_atom,
337 augmentation_charge,
338 total_hubbard_orbitals,
339 number_of_occupied_bands,
340 phi_beta,
341 beta_psi,
342 orbital_projector,
343 projector_band,
344 orbital_band,
345 projection);
346
347 std::vector<Complex> strained_orbitals(orbitals.size());
348 std::vector<Complex> strained_betas(beta_projectors.size());
349
350 for (int direction_a = 0; direction_a < 3; direction_a++) {
351 for (int direction_b = 0; direction_b < 3; direction_b++) {
352 apply_strain_derivative(number_of_plane_waves,
353 total_hubbard_orbitals,
354 geometry,
355 direction_a,
356 direction_b,
357 orbitals,
358 orbital_radial_derivative,
359 orbital_angular_derivative[direction_a],
360 strained_orbitals);
361 apply_strain_derivative(number_of_plane_waves,
362 total_projectors,
363 geometry,
364 direction_a,
365 direction_b,
366 beta_projectors,
367 beta_radial_derivative,
368 beta_angular_derivative[direction_a],
369 strained_betas);
370
371 /**
372 * dP_m/deps = <dphi_m|S|psi> + <phi_m|dS/deps|psi>, the same two
373 * pieces Quantum ESPRESSO assembles separately and adds.
374 * <dphi_m|S|psi> is expanded rather than built through an explicit
375 * S|psi>, exactly as the force does, so it cannot drift from the
376 * SCF's own definition of S.
377 */
378 auto derivative_projection = contract(strained_orbitals.data(),
379 total_hubbard_orbitals,
380 wavefunctions.data(),
381 number_of_occupied_bands);
382 auto strained_orbital_beta = contract(strained_orbitals.data(),
383 total_hubbard_orbitals,
384 beta_projectors.data(),
385 total_projectors);
386 auto orbital_strained_beta = contract(orbitals.data(),
387 total_hubbard_orbitals,
388 strained_betas.data(),
389 total_projectors);
390 auto strained_beta_psi = contract(strained_betas.data(),
391 total_projectors,
392 wavefunctions.data(),
393 number_of_occupied_bands);
394
395 accumulate_augmentation(number_of_atoms,
396 channels_by_atom,
397 projector_offset_by_atom,
398 augmentation_charge,
399 total_hubbard_orbitals,
400 number_of_occupied_bands,
401 strained_orbital_beta,
402 beta_psi,
403 orbital_projector,
404 projector_band,
405 orbital_band,
406 derivative_projection);
407 accumulate_augmentation(number_of_atoms,
408 channels_by_atom,
409 projector_offset_by_atom,
410 augmentation_charge,
411 total_hubbard_orbitals,
412 number_of_occupied_bands,
413 orbital_strained_beta,
414 beta_psi,
415 orbital_projector,
416 projector_band,
417 orbital_band,
418 derivative_projection);
419 accumulate_augmentation(number_of_atoms,
420 channels_by_atom,
421 projector_offset_by_atom,
422 augmentation_charge,
423 total_hubbard_orbitals,
424 number_of_occupied_bands,
425 phi_beta,
426 strained_beta_psi,
427 orbital_projector,
428 projector_band,
429 orbital_band,
430 derivative_projection);
431
432 for (size_t manifold_index = 0; manifold_index < number_of_manifolds;
433 manifold_index++) {
434 const auto& manifold = manifolds[manifold_index];
435 int manifold_size = static_cast<int>(manifold.orbital_indices.size());
436 int row_offset = manifold_row_offsets[manifold_index];
437
438 Matrix<Complex> manifold_overlaps(static_cast<size_t>(number_of_bands),
439 static_cast<size_t>(manifold_size));
440 Matrix<Complex> derivative_overlaps(static_cast<size_t>(number_of_bands),
441 static_cast<size_t>(manifold_size));
442 for (int m = 0; m < manifold_size; m++)
443 for (int occupied_index = 0; occupied_index < number_of_occupied_bands;
444 occupied_index++) {
445 size_t band_row = static_cast<size_t>(
446 occupied_band_indices[static_cast<size_t>(occupied_index)]);
447 manifold_overlaps(band_row, static_cast<size_t>(m)) =
448 projection[orbital_band(row_offset + m, occupied_index)];
449 derivative_overlaps(band_row, static_cast<size_t>(m)) =
450 derivative_projection[orbital_band(row_offset + m,
451 occupied_index)];
452 }
453
454 const auto& combined_d_matrix =
455 scf_state_.hubbard_d_matrix_per_spin[spin_index];
456 std::vector<Complex> manifold_d_matrix(static_cast<size_t>(manifold_size) *
457 manifold_size);
458 for (int m = 0; m < manifold_size; m++)
459 for (int m_prime = 0; m_prime < manifold_size; m_prime++)
460 manifold_d_matrix[static_cast<size_t>(m * manifold_size +
461 m_prime)] =
462 combined_d_matrix[static_cast<size_t>(row_offset + m) +
463 static_cast<size_t>(row_offset + m_prime) *
464 total_hubbard_orbitals];
465
466 /**
467 * The accumulator returns dE_U/dlambda for a projection derivative
468 * handed to it; sigma = -(1/Omega) dE/deps.
469 */
470 result_.hubbard[direction_a][direction_b] -=
471 inverse_cell_volume *
472 accumulate_hubbard_force_gradient(number_of_bands,
473 manifold_size,
474 manifold_overlaps,
475 derivative_overlaps,
476 band_occupations,
477 kpoint_weight,
478 manifold_d_matrix);
479 }
480 }
481 }
482 }
483
484 symmetrize_stress_tensor(result_.hubbard);
485 }
486
487 private:
488 /// One atomic orbital's place in the pseudopotential tables and in the cell.
489 struct AtomicOrbitalDescriptor {
490 int global_orbital_index{}; // row of build_atomic_orbitals_g_space' output
491 int atom_index{};
492 int species_type_index{};
493 int wavefunction_index{}; // row of tabulate_atomic_wavefunctions_in_q_space
494 int angular_momentum{};
495 int combined_lm_index{}; // l*l + m, the Y_lm storage order
496 };
497
498 /// Per-G geometry shared by every radial and angular derivative below.
499 struct PlaneWaveGeometry {
500 std::vector<std::array<Real, 3>> reciprocal_units; ///< k+G in units of 2pi/a
501 std::vector<Real> magnitude_squared_reciprocal_units; ///< |k+G|^2 in units of (2pi/a)^2
502 std::vector<std::array<Real, 3>> cartesian; ///< k+G in bohr^-1
503 std::vector<Real> magnitude_bohr; ///< |k+G| in bohr^-1
504 std::vector<Real> inverse_magnitude_bohr; ///< 1/|k+G|, zero at k+G = 0
505 std::vector<std::vector<Real>> spherical_harmonics; ///< [g][lm]
506 std::array<std::vector<std::vector<Real>>, 3>
507 spherical_harmonic_derivatives; ///< [a][g][lm]
508 std::vector<std::vector<Complex>> structure_factors; ///< [atom][g], exp(i(k+G).tau)
509 };
510
511 /**
512 * Compute ortho-atomic Hubbard stress.
513 *
514 * The SCF projector is the system-wide Lowdin transform
515 *
516 * |orthogonal_I> = sum_J |atomic_J> (O^-1/2)_JI,
517 * O_IJ = <atomic_I|S|atomic_J>.
518 *
519 * A homogeneous strain moves every atomic orbital, so both terms in
520 *
521 * dP_I = sum_J d(O^-1/2)_IJ <atomic_J|psi>
522 * + sum_J (O^-1/2)_IJ <d atomic_J|psi>
523 *
524 * run over the complete cell-wide orbital set. The overlap derivative is assembled from
525 * the same analytic orbital derivatives used by the atomic stress path, reduced before
526 * diagonalization, and handed to the shared Lowdin matrix-function derivative. Only after
527 * the complete derivative is formed are the Hubbard-manifold rows selected.
528 */
529 void compute_with_orthogonal_projectors() {
530 const auto& manifolds = scf_state_.hubbard_manifolds;
531 const auto& pseudopotential_basis = scf_state_.quantum_system().pseudopotential_basis();
532 const auto& pseudopotential_set = pseudopotential_basis.pseudopotential_set();
533 bool has_augmented_overlap = scf_state_.quantum_system().has_ultrasoft();
534
535 constexpr Real q_spacing = 0.01;
536 constexpr int number_of_q_points = 900;
537 auto atomic_wavefunction_table = tabulate_atomic_wavefunctions_in_q_space(
538 pseudopotential_set, crystal_structure_.cell_volume(), q_spacing, number_of_q_points);
539
540 auto atomic_orbital_descriptors = enumerate_atomic_orbital_descriptors();
541 int number_of_atomic_orbitals = static_cast<int>(atomic_orbital_descriptors.size());
542
543 size_t number_of_manifolds = manifolds.size();
544 std::vector<int> manifold_row_offsets(number_of_manifolds);
545 int total_hubbard_orbitals = 0;
546 for (size_t manifold_index = 0; manifold_index < number_of_manifolds; manifold_index++) {
547 manifold_row_offsets[manifold_index] = total_hubbard_orbitals;
548 total_hubbard_orbitals +=
549 static_cast<int>(manifolds[manifold_index].orbital_indices.size());
550 }
551 if (number_of_atomic_orbitals == 0 || total_hubbard_orbitals == 0)
552 return;
553
554 int number_of_atoms = crystal_structure_.number_of_atoms();
555 int total_projectors = pseudopotential_basis.total_number_of_projectors();
556 std::vector<int> projector_offset_by_atom(number_of_atoms, 0);
557 std::vector<int> projector_channels_by_atom(number_of_atoms, 0);
558 {
559 int projector_offset = 0;
560 for (int atom_index = 0; atom_index < number_of_atoms; atom_index++) {
561 int species_index = crystal_structure_.species_type_index(atom_index);
562 int number_of_projector_channels =
563 pseudopotential_basis.number_of_projector_channels(species_index);
564 projector_offset_by_atom[atom_index] = projector_offset;
565 projector_channels_by_atom[atom_index] = number_of_projector_channels;
566 projector_offset += number_of_projector_channels;
567 }
568 }
569 auto augmentation_overlap = [&](int atom_index, int channel_row, int channel_column) {
570 return hubbard_augmentation_overlap(
571 static_cast<const UltrasoftPlanewaveKohnShamQuantumSystem&>(
572 scf_state_.quantum_system())
573 .ultrasoft_basis(),
574 pseudopotential_basis,
575 crystal_structure_.species_type_index(atom_index),
576 channel_row,
577 channel_column);
578 };
579
580 int number_of_kpoints = kpoint_grid_.number_of_kpoints();
581 int number_of_bands = scf_state_.number_of_bands();
582 int number_of_kpoint_states = number_of_kpoints * scf_state_.number_of_spin_components();
583 const Matrix3x3& reciprocal_lattice_vectors =
584 crystal_structure_.reciprocal_lattice_vectors();
585 Real inverse_cell_volume = 1.0 / crystal_structure_.cell_volume();
586
587 auto orbital_matrix_element = [number_of_atomic_orbitals](int row, int column) {
588 return static_cast<size_t>(row + column * number_of_atomic_orbitals);
589 };
590
591 for (int kpoint_index_state = 0; kpoint_index_state < number_of_kpoint_states;
592 kpoint_index_state++) {
593 int spin_index = kpoint_index_state / number_of_kpoints;
594 int kpoint_index = kpoint_index_state % number_of_kpoints;
595 int number_of_plane_waves = scf_state_.number_of_plane_waves(kpoint_index);
596 const auto& plane_wave_indices = scf_state_.g_vector_indices_at_kpoint(kpoint_index);
597 Real kpoint_weight = kpoint_grid_[kpoint_index].weight;
598 bool is_gamma_only = planewave_basis_.is_gamma_only();
599 int g_zero_index = is_gamma_only ? find_g_zero_plane_wave_index(planewave_basis_,
600 plane_wave_indices,
601 number_of_plane_waves)
602 : -1;
603
604 Vector3D kpoint_cartesian;
605 Vector3D kpoint_vector = kpoint_grid_[kpoint_index].coordinates;
606 for (int direction = 0; direction < 3; direction++)
607 kpoint_cartesian[direction] =
608 kpoint_vector.x * reciprocal_lattice_vectors[0][direction] +
609 kpoint_vector.y * reciprocal_lattice_vectors[1][direction] +
610 kpoint_vector.z * reciprocal_lattice_vectors[2][direction];
611
612 PlaneWaveGeometry geometry = build_plane_wave_geometry(
613 kpoint_cartesian, plane_wave_indices, number_of_plane_waves);
614
615 std::vector<Real> band_occupations(number_of_bands);
616 std::vector<int> occupied_band_indices;
617 for (int band_index = 0; band_index < number_of_bands; band_index++) {
618 band_occupations[band_index] =
619 scf_state_.occupation(kpoint_index_state, band_index);
620 if (band_occupations[band_index] >= math::tolerances::occupation_threshold)
621 occupied_band_indices.push_back(band_index);
622 }
623 int number_of_occupied_bands = static_cast<int>(occupied_band_indices.size());
624 if (number_of_occupied_bands == 0)
625 continue;
626
627 std::vector<Complex> wavefunctions(static_cast<size_t>(number_of_plane_waves) *
628 number_of_occupied_bands);
629 for (int occupied_band_index = 0; occupied_band_index < number_of_occupied_bands;
630 occupied_band_index++)
631 for (int g_index = 0; g_index < number_of_plane_waves; g_index++)
632 wavefunctions[static_cast<size_t>(occupied_band_index) * number_of_plane_waves +
633 g_index] =
634 scf_state_.wavefunction_coefficient(
635 kpoint_index_state,
636 occupied_band_indices[static_cast<size_t>(occupied_band_index)],
637 g_index);
638
639 std::vector<Complex> atomic_orbitals;
640 std::vector<Complex> atomic_orbital_radial_derivative;
641 std::array<std::vector<Complex>, 3> atomic_orbital_angular_derivative;
642 build_atomic_orbital_strain_parts(kpoint_index,
643 plane_wave_indices,
644 number_of_plane_waves,
645 geometry,
646 atomic_orbital_descriptors,
647 atomic_wavefunction_table,
648 q_spacing,
649 number_of_q_points,
650 atomic_orbitals,
651 atomic_orbital_radial_derivative,
652 atomic_orbital_angular_derivative);
653
654 std::vector<Complex> beta_projectors;
655 std::vector<Complex> beta_radial_derivative;
656 std::array<std::vector<Complex>, 3> beta_angular_derivative;
657 if (has_augmented_overlap) {
658 build_beta_strain_parts(kpoint_index,
659 plane_wave_indices,
660 number_of_plane_waves,
661 total_projectors,
662 geometry,
663 beta_projectors,
664 beta_radial_derivative,
665 beta_angular_derivative);
666 }
667
668 auto contract = [&](const Complex* left,
669 int number_of_left,
670 const Complex* right,
671 int number_of_right) {
672 std::vector<Complex> contraction;
673 compute_hubbard_projector_overlaps(number_of_plane_waves,
674 number_of_left,
675 number_of_right,
676 left,
677 right,
678 contraction);
679 if (is_gamma_only)
680 complete_gamma_only_overlaps(number_of_plane_waves,
681 number_of_left,
682 number_of_right,
683 left,
684 right,
685 g_zero_index,
686 contraction);
687 if (scf_state_.execution_policy != nullptr)
688 scf_state_.execution_policy->reduce_within_group(
689 contraction.data(), static_cast<int>(contraction.size()));
690 return contraction;
691 };
692
693 if (has_augmented_overlap) {
694 if (scf_state_.execution_policy == nullptr) {
695 throw std::runtime_error(
696 "Ortho-atomic ultrasoft Hubbard stress requires an execution policy");
697 }
698
699 Matrix<Complex> canonical_beta_projectors(
700 static_cast<size_t>(total_projectors),
701 static_cast<size_t>(number_of_plane_waves));
702 std::copy(beta_projectors.begin(),
703 beta_projectors.end(),
704 canonical_beta_projectors.data());
705 auto overlap_operator = create_overlap_operator(crystal_structure_,
706 scf_state_.quantum_system(),
707 canonical_beta_projectors,
708 *scf_state_.execution_policy);
709
710 std::vector<Complex> overlap_applied_orbitals(atomic_orbitals.size());
711 overlap_operator->apply_batch(number_of_plane_waves,
712 number_of_atomic_orbitals,
713 atomic_orbitals.data(),
714 overlap_applied_orbitals.data());
715 std::vector<Complex> overlap_applied_wavefunctions(wavefunctions.size());
716 overlap_operator->apply_batch(number_of_plane_waves,
717 number_of_occupied_bands,
718 wavefunctions.data(),
719 overlap_applied_wavefunctions.data());
720
721 auto overlap_matrix = contract(atomic_orbitals.data(),
722 number_of_atomic_orbitals,
723 overlap_applied_orbitals.data(),
724 number_of_atomic_orbitals);
725 auto lowdin_transform =
726 compute_lowdin_transform(overlap_matrix, number_of_atomic_orbitals);
727 auto atomic_band_overlaps = contract(atomic_orbitals.data(),
728 number_of_atomic_orbitals,
729 overlap_applied_wavefunctions.data(),
730 number_of_occupied_bands);
731
732 std::vector<Complex> projection(static_cast<size_t>(number_of_atomic_orbitals) *
733 number_of_occupied_bands,
734 Complex(0.0, 0.0));
735 for (int occupied_band_index = 0; occupied_band_index < number_of_occupied_bands;
736 occupied_band_index++) {
737 for (int orbital_index = 0; orbital_index < number_of_atomic_orbitals;
738 orbital_index++) {
739 for (int inner_orbital_index = 0;
740 inner_orbital_index < number_of_atomic_orbitals;
741 inner_orbital_index++) {
742 projection[orbital_matrix_element(orbital_index,
743 occupied_band_index)] +=
744 lowdin_transform.inverse_square_root[orbital_matrix_element(
745 orbital_index, inner_orbital_index)] *
746 atomic_band_overlaps[orbital_matrix_element(inner_orbital_index,
747 occupied_band_index)];
748 }
749 }
750 }
751
752 auto apply_overlap_derivative = [&](const Complex* input_vectors,
753 int number_of_vectors,
754 const std::vector<Complex>&
755 strained_beta_projectors) {
756 auto projector_overlaps = contract(
757 beta_projectors.data(), total_projectors, input_vectors, number_of_vectors);
758 auto strained_projector_overlaps = contract(strained_beta_projectors.data(),
759 total_projectors,
760 input_vectors,
761 number_of_vectors);
762 std::vector<Complex> result(static_cast<size_t>(number_of_plane_waves) *
763 number_of_vectors,
764 Complex(0.0, 0.0));
765 for (int atom_index = 0; atom_index < number_of_atoms; atom_index++) {
766 int projector_offset = projector_offset_by_atom[atom_index];
767 int number_of_projector_channels = projector_channels_by_atom[atom_index];
768 for (int channel_row = 0; channel_row < number_of_projector_channels;
769 channel_row++) {
770 int projector_row = projector_offset + channel_row;
771 for (int channel_column = 0;
772 channel_column < number_of_projector_channels;
773 channel_column++) {
774 Real overlap_value =
775 augmentation_overlap(atom_index, channel_row, channel_column);
776 if (std::abs(overlap_value) <
777 math::tolerances::d_matrix_element_threshold) {
778 continue;
779 }
780 int projector_column = projector_offset + channel_column;
781 for (int vector_index = 0; vector_index < number_of_vectors;
782 vector_index++) {
783 size_t projector_overlap_index =
784 static_cast<size_t>(projector_column) +
785 static_cast<size_t>(vector_index) * total_projectors;
786 Complex strained_projector_coefficient =
787 overlap_value * projector_overlaps[projector_overlap_index];
788 Complex projector_coefficient =
789 overlap_value *
790 strained_projector_overlaps[projector_overlap_index];
791 for (int g_index = 0; g_index < number_of_plane_waves;
792 g_index++) {
793 size_t result_index = static_cast<size_t>(vector_index) *
794 number_of_plane_waves +
795 g_index;
796 result[result_index] +=
797 strained_beta_projectors[static_cast<size_t>(
798 projector_row) *
799 number_of_plane_waves +
800 g_index] *
801 strained_projector_coefficient +
802 beta_projectors[static_cast<size_t>(projector_row) *
803 number_of_plane_waves +
804 g_index] *
805 projector_coefficient;
806 }
807 }
808 }
809 }
810 }
811 return result;
812 };
813
814 std::vector<Complex> strained_atomic_orbitals(atomic_orbitals.size());
815 std::vector<Complex> strained_beta_projectors(beta_projectors.size());
816 for (int direction_a = 0; direction_a < 3; direction_a++) {
817 for (int direction_b = 0; direction_b < 3; direction_b++) {
818 apply_strain_derivative(number_of_plane_waves,
819 number_of_atomic_orbitals,
820 geometry,
821 direction_a,
822 direction_b,
823 atomic_orbitals,
824 atomic_orbital_radial_derivative,
825 atomic_orbital_angular_derivative[direction_a],
826 strained_atomic_orbitals);
827 apply_strain_derivative(number_of_plane_waves,
828 total_projectors,
829 geometry,
830 direction_a,
831 direction_b,
832 beta_projectors,
833 beta_radial_derivative,
834 beta_angular_derivative[direction_a],
835 strained_beta_projectors);
836
837 std::vector<Complex> overlap_applied_strained_orbitals(
838 strained_atomic_orbitals.size());
839 overlap_operator->apply_batch(number_of_plane_waves,
840 number_of_atomic_orbitals,
841 strained_atomic_orbitals.data(),
842 overlap_applied_strained_orbitals.data());
843 auto overlap_derivative_on_orbitals =
844 apply_overlap_derivative(atomic_orbitals.data(),
845 number_of_atomic_orbitals,
846 strained_beta_projectors);
847 auto overlap_derivative_on_wavefunctions =
848 apply_overlap_derivative(wavefunctions.data(),
849 number_of_occupied_bands,
850 strained_beta_projectors);
851
852 auto left_orbital_metric_derivative =
853 contract(strained_atomic_orbitals.data(),
854 number_of_atomic_orbitals,
855 overlap_applied_orbitals.data(),
856 number_of_atomic_orbitals);
857 auto overlap_metric_derivative =
858 contract(atomic_orbitals.data(),
859 number_of_atomic_orbitals,
860 overlap_derivative_on_orbitals.data(),
861 number_of_atomic_orbitals);
862 auto right_orbital_metric_derivative =
863 contract(atomic_orbitals.data(),
864 number_of_atomic_orbitals,
865 overlap_applied_strained_orbitals.data(),
866 number_of_atomic_orbitals);
867 std::vector<Complex> metric_derivative(overlap_matrix.size());
868 for (size_t matrix_index = 0; matrix_index < metric_derivative.size();
869 matrix_index++) {
870 metric_derivative[matrix_index] =
871 left_orbital_metric_derivative[matrix_index] +
872 overlap_metric_derivative[matrix_index] +
873 right_orbital_metric_derivative[matrix_index];
874 }
875 for (int column = 0; column < number_of_atomic_orbitals; column++) {
876 for (int row = 0; row <= column; row++) {
877 Complex hermitian_value =
878 0.5 *
879 (metric_derivative[orbital_matrix_element(row, column)] +
880 std::conj(
881 metric_derivative[orbital_matrix_element(column, row)]));
882 metric_derivative[orbital_matrix_element(row, column)] =
883 hermitian_value;
884 metric_derivative[orbital_matrix_element(column, row)] =
885 std::conj(hermitian_value);
886 }
887 }
888 auto inverse_square_root_derivative = lowdin_inverse_square_root_derivative(
889 lowdin_transform, metric_derivative);
890
891 auto left_orbital_projection_derivative =
892 contract(strained_atomic_orbitals.data(),
893 number_of_atomic_orbitals,
894 overlap_applied_wavefunctions.data(),
895 number_of_occupied_bands);
896 auto overlap_projection_derivative =
897 contract(atomic_orbitals.data(),
898 number_of_atomic_orbitals,
899 overlap_derivative_on_wavefunctions.data(),
900 number_of_occupied_bands);
901 std::vector<Complex> raw_projection_derivative(atomic_band_overlaps.size());
902 for (size_t overlap_index = 0;
903 overlap_index < raw_projection_derivative.size();
904 overlap_index++) {
905 raw_projection_derivative[overlap_index] =
906 left_orbital_projection_derivative[overlap_index] +
907 overlap_projection_derivative[overlap_index];
908 }
909
910 std::vector<Complex> projection_derivative(projection.size(),
911 Complex(0.0, 0.0));
912 for (int occupied_band_index = 0;
913 occupied_band_index < number_of_occupied_bands;
914 occupied_band_index++) {
915 for (int orbital_index = 0; orbital_index < number_of_atomic_orbitals;
916 orbital_index++) {
917 for (int inner_orbital_index = 0;
918 inner_orbital_index < number_of_atomic_orbitals;
919 inner_orbital_index++) {
920 projection_derivative[orbital_matrix_element(
921 orbital_index, occupied_band_index)] +=
922 inverse_square_root_derivative[orbital_matrix_element(
923 orbital_index, inner_orbital_index)] *
924 atomic_band_overlaps[orbital_matrix_element(
925 inner_orbital_index, occupied_band_index)] +
926 lowdin_transform.inverse_square_root[orbital_matrix_element(
927 orbital_index, inner_orbital_index)] *
928 raw_projection_derivative[orbital_matrix_element(
929 inner_orbital_index, occupied_band_index)];
930 }
931 }
932 }
933
934 for (size_t manifold_index = 0; manifold_index < number_of_manifolds;
935 manifold_index++) {
936 const auto& manifold = manifolds[manifold_index];
937 int manifold_size = static_cast<int>(manifold.orbital_indices.size());
938 int row_offset = manifold_row_offsets[manifold_index];
939 Matrix<Complex> manifold_overlaps(static_cast<size_t>(number_of_bands),
940 static_cast<size_t>(manifold_size));
941 Matrix<Complex> derivative_overlaps(
942 static_cast<size_t>(number_of_bands),
943 static_cast<size_t>(manifold_size));
944 for (int orbital_within_manifold = 0;
945 orbital_within_manifold < manifold_size;
946 orbital_within_manifold++) {
947 int orbital_index = manifold.orbital_indices[static_cast<size_t>(
948 orbital_within_manifold)];
949 for (int occupied_band_index = 0;
950 occupied_band_index < number_of_occupied_bands;
951 occupied_band_index++) {
952 size_t band_row = static_cast<size_t>(
953 occupied_band_indices[static_cast<size_t>(
954 occupied_band_index)]);
955 manifold_overlaps(
956 band_row, static_cast<size_t>(orbital_within_manifold)) =
957 projection[orbital_matrix_element(orbital_index,
958 occupied_band_index)];
959 derivative_overlaps(
960 band_row, static_cast<size_t>(orbital_within_manifold)) =
961 projection_derivative[orbital_matrix_element(
962 orbital_index, occupied_band_index)];
963 }
964 }
965
966 const auto& combined_d_matrix =
967 scf_state_.hubbard_d_matrix_per_spin[spin_index];
968 std::vector<Complex> manifold_d_matrix(
969 static_cast<size_t>(manifold_size) * manifold_size);
970 for (int row = 0; row < manifold_size; row++) {
971 for (int column = 0; column < manifold_size; column++) {
972 manifold_d_matrix[static_cast<size_t>(row * manifold_size +
973 column)] =
974 combined_d_matrix[static_cast<size_t>(row_offset + row) +
975 static_cast<size_t>(row_offset + column) *
976 total_hubbard_orbitals];
977 }
978 }
979
980 result_.hubbard[direction_a][direction_b] -=
981 inverse_cell_volume *
982 accumulate_hubbard_force_gradient(number_of_bands,
983 manifold_size,
984 manifold_overlaps,
985 derivative_overlaps,
986 band_occupations,
987 kpoint_weight,
988 manifold_d_matrix);
989 }
990 }
991 }
992 continue;
993 }
994
995 auto overlap_matrix = contract(atomic_orbitals.data(),
996 number_of_atomic_orbitals,
997 atomic_orbitals.data(),
998 number_of_atomic_orbitals);
999 auto lowdin_transform =
1000 compute_lowdin_transform(overlap_matrix, number_of_atomic_orbitals);
1001 auto atomic_band_overlaps = contract(atomic_orbitals.data(),
1002 number_of_atomic_orbitals,
1003 wavefunctions.data(),
1004 number_of_occupied_bands);
1005
1006 std::vector<Complex> projection(static_cast<size_t>(number_of_atomic_orbitals) *
1007 number_of_occupied_bands,
1008 Complex(0.0, 0.0));
1009 for (int occupied_band_index = 0; occupied_band_index < number_of_occupied_bands;
1010 occupied_band_index++)
1011 for (int orbital_index = 0; orbital_index < number_of_atomic_orbitals;
1012 orbital_index++)
1013 for (int inner_orbital_index = 0;
1014 inner_orbital_index < number_of_atomic_orbitals;
1015 inner_orbital_index++)
1016 projection[orbital_matrix_element(orbital_index, occupied_band_index)] +=
1017 lowdin_transform.inverse_square_root[orbital_matrix_element(
1018 orbital_index, inner_orbital_index)] *
1019 atomic_band_overlaps[orbital_matrix_element(inner_orbital_index,
1020 occupied_band_index)];
1021
1022 std::vector<Complex> strained_atomic_orbitals(atomic_orbitals.size());
1023 for (int direction_a = 0; direction_a < 3; direction_a++) {
1024 for (int direction_b = 0; direction_b < 3; direction_b++) {
1025 apply_strain_derivative(number_of_plane_waves,
1026 number_of_atomic_orbitals,
1027 geometry,
1028 direction_a,
1029 direction_b,
1030 atomic_orbitals,
1031 atomic_orbital_radial_derivative,
1032 atomic_orbital_angular_derivative[direction_a],
1033 strained_atomic_orbitals);
1034
1035 auto atomic_band_derivative = contract(strained_atomic_orbitals.data(),
1036 number_of_atomic_orbitals,
1037 wavefunctions.data(),
1038 number_of_occupied_bands);
1039 auto left_overlap_derivative = contract(strained_atomic_orbitals.data(),
1040 number_of_atomic_orbitals,
1041 atomic_orbitals.data(),
1042 number_of_atomic_orbitals);
1043 std::vector<Complex> overlap_derivative(
1044 static_cast<size_t>(number_of_atomic_orbitals) * number_of_atomic_orbitals);
1045 for (int column = 0; column < number_of_atomic_orbitals; column++)
1046 for (int row = 0; row < number_of_atomic_orbitals; row++)
1047 overlap_derivative[orbital_matrix_element(row, column)] =
1048 left_overlap_derivative[orbital_matrix_element(row, column)] +
1049 std::conj(
1050 left_overlap_derivative[orbital_matrix_element(column, row)]);
1051
1052 auto inverse_square_root_derivative =
1053 lowdin_inverse_square_root_derivative(lowdin_transform, overlap_derivative);
1054 std::vector<Complex> projection_derivative(
1055 static_cast<size_t>(number_of_atomic_orbitals) * number_of_occupied_bands,
1056 Complex(0.0, 0.0));
1057 for (int occupied_band_index = 0;
1058 occupied_band_index < number_of_occupied_bands;
1059 occupied_band_index++)
1060 for (int orbital_index = 0; orbital_index < number_of_atomic_orbitals;
1061 orbital_index++)
1062 for (int inner_orbital_index = 0;
1063 inner_orbital_index < number_of_atomic_orbitals;
1064 inner_orbital_index++)
1065 projection_derivative[orbital_matrix_element(
1066 orbital_index, occupied_band_index)] +=
1067 inverse_square_root_derivative[orbital_matrix_element(
1068 orbital_index, inner_orbital_index)] *
1069 atomic_band_overlaps[orbital_matrix_element(
1070 inner_orbital_index, occupied_band_index)] +
1071 lowdin_transform.inverse_square_root[orbital_matrix_element(
1072 orbital_index, inner_orbital_index)] *
1073 atomic_band_derivative[orbital_matrix_element(
1074 inner_orbital_index, occupied_band_index)];
1075
1076 for (size_t manifold_index = 0; manifold_index < number_of_manifolds;
1077 manifold_index++) {
1078 const auto& manifold = manifolds[manifold_index];
1079 int manifold_size = static_cast<int>(manifold.orbital_indices.size());
1080 int row_offset = manifold_row_offsets[manifold_index];
1081 Matrix<Complex> manifold_overlaps(static_cast<size_t>(number_of_bands),
1082 static_cast<size_t>(manifold_size));
1083 Matrix<Complex> derivative_overlaps(static_cast<size_t>(number_of_bands),
1084 static_cast<size_t>(manifold_size));
1085 for (int orbital_within_manifold = 0;
1086 orbital_within_manifold < manifold_size;
1087 orbital_within_manifold++) {
1088 int orbital_index =
1089 manifold
1090 .orbital_indices[static_cast<size_t>(orbital_within_manifold)];
1091 for (int occupied_band_index = 0;
1092 occupied_band_index < number_of_occupied_bands;
1093 occupied_band_index++) {
1094 size_t band_row =
1095 static_cast<size_t>(occupied_band_indices[static_cast<size_t>(
1096 occupied_band_index)]);
1097 manifold_overlaps(band_row,
1098 static_cast<size_t>(orbital_within_manifold)) =
1099 projection[orbital_matrix_element(orbital_index,
1100 occupied_band_index)];
1101 derivative_overlaps(band_row,
1102 static_cast<size_t>(orbital_within_manifold)) =
1103 projection_derivative[orbital_matrix_element(
1104 orbital_index, occupied_band_index)];
1105 }
1106 }
1107
1108 const auto& combined_d_matrix =
1109 scf_state_.hubbard_d_matrix_per_spin[spin_index];
1110 std::vector<Complex> manifold_d_matrix(static_cast<size_t>(manifold_size) *
1111 manifold_size);
1112 for (int row = 0; row < manifold_size; row++)
1113 for (int column = 0; column < manifold_size; column++)
1114 manifold_d_matrix[static_cast<size_t>(row * manifold_size +
1115 column)] =
1116 combined_d_matrix[static_cast<size_t>(row_offset + row) +
1117 static_cast<size_t>(row_offset + column) *
1118 total_hubbard_orbitals];
1119
1120 result_.hubbard[direction_a][direction_b] -=
1121 inverse_cell_volume *
1122 accumulate_hubbard_force_gradient(number_of_bands,
1123 manifold_size,
1124 manifold_overlaps,
1125 derivative_overlaps,
1126 band_occupations,
1127 kpoint_weight,
1128 manifold_d_matrix);
1129 }
1130 }
1131 }
1132 }
1133
1134 symmetrize_stress_tensor(result_.hubbard);
1135 }
1136
1137 /**
1138 * Every atomic orbital's descriptor, in `enumerate_atomic_orbitals`' order -- the order
1139 * `build_atomic_orbitals_g_space` fills, and the order `HubbardManifold::orbital_indices`
1140 * indexes into. Walked here rather than taken from `AtomicOrbitalLabel` because the
1141 * radial table is addressed by the shell's index within the dataset, which the label
1142 * does not carry.
1143 */
1144 [[nodiscard]] std::vector<AtomicOrbitalDescriptor>
1145 enumerate_atomic_orbital_descriptors() const {
1146 const auto& pseudopotential_set =
1147 scf_state_.quantum_system().pseudopotential_basis().pseudopotential_set();
1148 std::vector<AtomicOrbitalDescriptor> descriptors;
1149 for (int atom_index = 0; atom_index < crystal_structure_.number_of_atoms(); atom_index++) {
1150 int species_type_index = crystal_structure_.species_type_index(atom_index);
1151 const auto& pseudopotential = pseudopotential_set[species_type_index];
1152 for (int wavefunction_index = 0;
1153 wavefunction_index < static_cast<int>(pseudopotential.wavefunctions_atomic.size());
1154 wavefunction_index++) {
1155 const auto& atomic_wavefunction =
1156 pseudopotential.wavefunctions_atomic[static_cast<size_t>(wavefunction_index)];
1157 if (atomic_wavefunction.occupations < 0.0)
1158 continue;
1159 int angular_momentum = atomic_wavefunction.angular_momentum;
1160 for (int magnetic_index = 0; magnetic_index < 2 * angular_momentum + 1;
1161 magnetic_index++)
1162 descriptors.push_back({static_cast<int>(descriptors.size()),
1163 atom_index,
1164 species_type_index,
1165 wavefunction_index,
1166 angular_momentum,
1167 angular_momentum * angular_momentum + magnetic_index});
1168 }
1169 }
1170 return descriptors;
1171 }
1172
1173 /**
1174 * k+G geometry, Y_lm, dY_lm/dq and the per-atom structure factors at one k-point.
1175 *
1176 * Y_lm and its derivative are evaluated on k+G in units of 2pi/a, which is what
1177 * `build_atomic_orbitals_g_space` and `compute_kleinman_bylander_projectors` both do;
1178 * the radial tables are addressed in bohr^-1. `lmax` covers both the atomic shells and
1179 * the beta channels so one batch serves both.
1180 */
1181 [[nodiscard]] PlaneWaveGeometry
1182 build_plane_wave_geometry(const Vector3D& kpoint_cartesian,
1183 const std::vector<int>& plane_wave_indices,
1184 int number_of_plane_waves) const {
1185 const auto& basis = scf_state_.quantum_system().pseudopotential_basis();
1186 const auto& pseudopotential_set = basis.pseudopotential_set();
1187 Real two_pi_over_lattice_parameter = crystal_structure_.two_pi_over_lattice_parameter();
1188
1189 PlaneWaveGeometry geometry;
1190 geometry.reciprocal_units.resize(number_of_plane_waves);
1191 geometry.magnitude_squared_reciprocal_units.resize(number_of_plane_waves);
1192 geometry.cartesian.resize(number_of_plane_waves);
1193 geometry.magnitude_bohr.resize(number_of_plane_waves);
1194 geometry.inverse_magnitude_bohr.resize(number_of_plane_waves);
1195
1196 for (int g_index = 0; g_index < number_of_plane_waves; g_index++) {
1197 const auto& g_vector =
1198 planewave_basis_.g_vectors()[plane_wave_indices[static_cast<size_t>(g_index)]];
1199 std::array<Real, 3> k_plus_g = {kpoint_cartesian.x + g_vector.g_cart.x,
1200 kpoint_cartesian.y + g_vector.g_cart.y,
1201 kpoint_cartesian.z + g_vector.g_cart.z};
1202 geometry.reciprocal_units[g_index] = k_plus_g;
1203 Real magnitude_squared =
1204 k_plus_g[0] * k_plus_g[0] + k_plus_g[1] * k_plus_g[1] + k_plus_g[2] * k_plus_g[2];
1205 geometry.magnitude_squared_reciprocal_units[g_index] = magnitude_squared;
1206 for (int direction = 0; direction < 3; direction++)
1207 geometry.cartesian[g_index][direction] =
1208 k_plus_g[direction] * two_pi_over_lattice_parameter;
1209 Real magnitude_bohr = std::sqrt(magnitude_squared) * two_pi_over_lattice_parameter;
1210 geometry.magnitude_bohr[g_index] = magnitude_bohr;
1211 geometry.inverse_magnitude_bohr[g_index] =
1212 (magnitude_bohr > math::tolerances::bessel_argument_zero_threshold)
1213 ? 1.0 / magnitude_bohr
1214 : 0.0;
1215 }
1216
1217 int maximum_angular_momentum = 0;
1218 for (int species_index = 0; species_index < static_cast<int>(pseudopotential_set.size());
1219 species_index++) {
1220 for (const auto& atomic_wavefunction :
1221 pseudopotential_set[species_index].wavefunctions_atomic)
1222 if (atomic_wavefunction.occupations >= 0.0)
1223 maximum_angular_momentum =
1224 std::max(maximum_angular_momentum, atomic_wavefunction.angular_momentum);
1225 for (const auto& projector : pseudopotential_set[species_index].projectors)
1226 maximum_angular_momentum =
1227 std::max(maximum_angular_momentum, projector.angular_momentum);
1228 }
1229 int number_of_lm_components =
1230 (maximum_angular_momentum + 1) * (maximum_angular_momentum + 1);
1231
1232 compute_real_spherical_harmonics_batch(number_of_lm_components,
1233 number_of_plane_waves,
1234 geometry.reciprocal_units,
1235 geometry.magnitude_squared_reciprocal_units,
1236 geometry.spherical_harmonics);
1237 for (int direction = 0; direction < 3; direction++)
1238 compute_spherical_harmonic_derivatives_batch(
1239 number_of_lm_components,
1240 number_of_plane_waves,
1241 geometry.reciprocal_units,
1242 geometry.magnitude_squared_reciprocal_units,
1243 direction,
1244 geometry.spherical_harmonic_derivatives[static_cast<size_t>(direction)]);
1245
1246 geometry.structure_factors.resize(
1247 static_cast<size_t>(crystal_structure_.number_of_atoms()));
1248 for (int atom_index = 0; atom_index < crystal_structure_.number_of_atoms(); atom_index++) {
1249 Vector3D atomic_position =
1250 crystal_structure_.convert_from_crystal_to_cartesian(atom_index);
1251 auto& factors = geometry.structure_factors[static_cast<size_t>(atom_index)];
1252 factors.resize(number_of_plane_waves);
1253 for (int g_index = 0; g_index < number_of_plane_waves; g_index++) {
1254 const auto& k_plus_g = geometry.reciprocal_units[g_index];
1255 Real phase = constants::two_pi *
1256 (k_plus_g[0] * atomic_position.x + k_plus_g[1] * atomic_position.y +
1257 k_plus_g[2] * atomic_position.z);
1258 factors[g_index] = compute_phase_factor(phase);
1259 }
1260 }
1261 return geometry;
1262 }
1263
1264 /**
1265 * The Hubbard orbitals and the two halves of their strain derivative, packed
1266 * column-major (number_of_plane_waves x total_hubbard_orbitals).
1267 *
1268 * `orbitals` comes from the shared `build_atomic_orbitals_g_space` -- the same call the
1269 * SCF and the force make -- so what is differentiated here is exactly the projector the
1270 * energy used. The radial half is dchi/dq from the same tabulated chi (0.01 spacing,
1271 * 900 points, smooth), the angular half is chi times dY/dq_a with the 2pi/a unit
1272 * conversion `NonlocalPseudopotentialStressTensorCalculator` applies to its own.
1273 */
1274 void build_atomic_orbital_strain_parts(
1275 int kpoint_index,
1276 const std::vector<int>& plane_wave_indices,
1277 int number_of_plane_waves,
1278 const PlaneWaveGeometry& geometry,
1279 const std::vector<AtomicOrbitalDescriptor>& hubbard_descriptors,
1280 const std::vector<std::vector<std::vector<Real>>>& atomic_wavefunction_table,
1281 Real q_spacing,
1282 int number_of_q_points,
1283 std::vector<Complex>& orbitals,
1284 std::vector<Complex>& radial_derivative,
1285 std::array<std::vector<Complex>, 3>& angular_derivative) const {
1286 const auto& pseudopotential_set =
1287 scf_state_.quantum_system().pseudopotential_basis().pseudopotential_set();
1288 int total_hubbard_orbitals = static_cast<int>(hubbard_descriptors.size());
1289 Real two_pi_over_lattice_parameter = crystal_structure_.two_pi_over_lattice_parameter();
1290
1291 std::vector<std::vector<Complex>> all_atomic_orbitals;
1292 build_atomic_orbitals_g_space(kpoint_index,
1293 number_of_plane_waves,
1294 plane_wave_indices,
1295 crystal_structure_,
1296 kpoint_grid_,
1297 planewave_basis_,
1298 pseudopotential_set,
1299 atomic_wavefunction_table,
1300 q_spacing,
1301 number_of_q_points,
1302 all_atomic_orbitals);
1303
1304 size_t block = static_cast<size_t>(number_of_plane_waves) * total_hubbard_orbitals;
1305 orbitals.assign(block, Complex(0.0, 0.0));
1306 radial_derivative.assign(block, Complex(0.0, 0.0));
1307 for (int direction = 0; direction < 3; direction++)
1308 angular_derivative[static_cast<size_t>(direction)].assign(block, Complex(0.0, 0.0));
1309
1310 for (int packed = 0; packed < total_hubbard_orbitals; packed++) {
1311 const auto& descriptor = hubbard_descriptors[static_cast<size_t>(packed)];
1312 int global_orbital = descriptor.global_orbital_index;
1313
1314 /**
1315 * i^l, NOT `compute_angular_phase_factor`, which is (-i)^l.
1316 * `build_atomic_orbitals_g_space` carries i^l on the atomic orbitals while
1317 * `compute_kleinman_bylander_projectors` carries (-i)^l on the betas; the two
1318 * conventions differ, and taking the beta one here would multiply every l = 1
1319 * derivative part by -1 relative to the orbital it differentiates -- worth a
1320 * factor of 16 in the silicon 3p stress, since it flips the radial and angular
1321 * halves against the 1/sqrt(Omega) half rather than dropping them.
1322 */
1323 Complex angular_phase =
1324 std::conj(compute_angular_phase_factor(descriptor.angular_momentum));
1325 const auto& structure_factor =
1326 geometry.structure_factors[static_cast<size_t>(descriptor.atom_index)];
1327 const auto& radial_table =
1328 atomic_wavefunction_table[static_cast<size_t>(descriptor.species_type_index)]
1329 [static_cast<size_t>(descriptor.wavefunction_index)];
1330 size_t column = static_cast<size_t>(packed) * number_of_plane_waves;
1331
1332 for (int g_index = 0; g_index < number_of_plane_waves; g_index++) {
1333 orbitals[column + g_index] =
1334 all_atomic_orbitals[static_cast<size_t>(global_orbital)]
1335 [static_cast<size_t>(g_index)];
1336
1337 Real magnitude = geometry.magnitude_bohr[g_index];
1338 Real radial_value = lagrange_interpolation_4point(
1339 radial_table, magnitude, q_spacing, number_of_q_points);
1340 Real radial_slope = lagrange_interpolation_4point_derivative(
1341 radial_table, magnitude, q_spacing, number_of_q_points);
1342 Real harmonic =
1343 geometry.spherical_harmonics[static_cast<size_t>(g_index)]
1344 [static_cast<size_t>(descriptor.combined_lm_index)];
1345 Complex common = structure_factor[static_cast<size_t>(g_index)] * angular_phase;
1346
1347 radial_derivative[column + g_index] = radial_slope * harmonic * common;
1348 for (int direction = 0; direction < 3; direction++) {
1349 Real harmonic_slope =
1350 geometry.spherical_harmonic_derivatives[static_cast<size_t>(direction)]
1351 [static_cast<size_t>(g_index)]
1352 [static_cast<size_t>(
1353 descriptor.combined_lm_index)];
1354 angular_derivative[static_cast<size_t>(direction)][column + g_index] =
1355 (radial_value / two_pi_over_lattice_parameter) * harmonic_slope * common;
1356 }
1357 }
1358 }
1359 }
1360
1361 /**
1362 * The beta projectors and the two halves of their strain derivative, packed
1363 * column-major (number_of_plane_waves x total_projectors).
1364 *
1365 * The same construction `NonlocalPseudopotentialStressTensorCalculator` uses for its
1366 * own dbeta/deps, re-expressed against the shared geometry so that beta and phi are
1367 * differentiated identically. The undifferentiated beta comes from the shared
1368 * `compute_kleinman_bylander_projectors` rather than being rebuilt here, for the same
1369 * reason phi comes from `build_atomic_orbitals_g_space`: it is what the SCF's own S is
1370 * made of. `total_projectors == 0` (norm-conserving) leaves every array empty and
1371 * switches the whole dS/deps assembly off.
1372 */
1373 void build_beta_strain_parts(int kpoint_index,
1374 const std::vector<int>& plane_wave_indices,
1375 int number_of_plane_waves,
1376 int total_projectors,
1377 const PlaneWaveGeometry& geometry,
1378 std::vector<Complex>& beta_projectors,
1379 std::vector<Complex>& radial_derivative,
1380 std::array<std::vector<Complex>, 3>& angular_derivative) const {
1381 size_t block = static_cast<size_t>(number_of_plane_waves) * total_projectors;
1382 beta_projectors.assign(block, Complex(0.0, 0.0));
1383 radial_derivative.assign(block, Complex(0.0, 0.0));
1384 for (int direction = 0; direction < 3; direction++)
1385 angular_derivative[static_cast<size_t>(direction)].assign(block, Complex(0.0, 0.0));
1386 if (total_projectors == 0)
1387 return;
1388
1389 const auto& basis = scf_state_.quantum_system().pseudopotential_basis();
1390 Matrix<Complex> canonical_projectors;
1391 compute_kleinman_bylander_projectors(kpoint_grid_[kpoint_index].coordinates,
1392 crystal_structure_,
1393 planewave_basis_,
1394 basis.pseudopotential_set(),
1395 basis.beta_function_table(),
1396 plane_wave_indices,
1397 number_of_plane_waves,
1398 total_projectors,
1399 canonical_projectors);
1400 for (int projector = 0; projector < total_projectors; projector++)
1401 for (int g_index = 0; g_index < number_of_plane_waves; g_index++)
1402 beta_projectors[static_cast<size_t>(projector) * number_of_plane_waves + g_index] =
1403 canonical_projectors[projector][g_index];
1404
1405 Real inverse_sqrt_cell_volume = 1.0 / std::sqrt(crystal_structure_.cell_volume());
1406 Real two_pi_over_lattice_parameter = crystal_structure_.two_pi_over_lattice_parameter();
1407
1408 int projector_index = 0;
1409 for (int atom_index = 0; atom_index < crystal_structure_.number_of_atoms(); atom_index++) {
1410 int species_type_index = crystal_structure_.species_type_index(atom_index);
1411 const auto& structure_factor =
1412 geometry.structure_factors[static_cast<size_t>(atom_index)];
1413 for (int beta_index = 0; beta_index < basis.number_of_projectors(species_type_index);
1414 beta_index++) {
1415 int angular_momentum = basis.beta_angular_momentum(species_type_index, beta_index);
1416 Complex angular_phase = compute_angular_phase_factor(angular_momentum);
1417 Complex prefactor = inverse_sqrt_cell_volume * angular_phase;
1418
1419 for (int magnetic_projection : harmonic_order_sequence(angular_momentum)) {
1420 int lm_index =
1421 harmonic_combined_degree_order_index(angular_momentum, magnetic_projection);
1422 size_t column = static_cast<size_t>(projector_index) * number_of_plane_waves;
1423 for (int g_index = 0; g_index < number_of_plane_waves; g_index++) {
1424 Real magnitude = geometry.magnitude_bohr[g_index];
1425 Real beta_value =
1426 basis.interpolate_beta(species_type_index, beta_index, magnitude);
1427 Real beta_slope =
1428 basis.interpolate_dbeta(species_type_index, beta_index, magnitude);
1429 Real harmonic = geometry.spherical_harmonics[static_cast<size_t>(g_index)]
1430 [static_cast<size_t>(lm_index)];
1431 Complex common = structure_factor[static_cast<size_t>(g_index)] * prefactor;
1432
1433 radial_derivative[column + g_index] = beta_slope * harmonic * common;
1434 for (int direction = 0; direction < 3; direction++) {
1435 Real harmonic_slope =
1436 geometry
1437 .spherical_harmonic_derivatives[static_cast<size_t>(direction)]
1438 [static_cast<size_t>(g_index)]
1439 [static_cast<size_t>(lm_index)];
1440 angular_derivative[static_cast<size_t>(direction)][column + g_index] =
1441 (beta_value / two_pi_over_lattice_parameter) * harmonic_slope *
1442 common;
1443 }
1444 }
1445 projector_index++;
1446 }
1447 }
1448 }
1449 }
1450
1451 /**
1452 * dA/deps_ab = -[ radial_a_b + angular_a * (k+G)_b ] - delta_ab A/2, in place.
1453 *
1454 * `radial` is f'(q) Y_lm S i^l /sqrt(Omega) and is weighted here by
1455 * (k+G)_a (k+G)_b/|k+G|; `angular` is f(q) (dY_lm/dq_a) S i^l /sqrt(Omega) and is
1456 * weighted by (k+G)_b. Quantum ESPRESSO carries those same two weights, with the
1457 * overall minus already folded into them.
1458 */
1459 static void apply_strain_derivative(int number_of_plane_waves,
1460 int number_of_functions,
1461 const PlaneWaveGeometry& geometry,
1462 int direction_a,
1463 int direction_b,
1464 const std::vector<Complex>& base,
1465 const std::vector<Complex>& radial,
1466 const std::vector<Complex>& angular,
1467 std::vector<Complex>& strained) {
1468 strained.assign(base.size(), Complex(0.0, 0.0));
1469 Real diagonal = (direction_a == direction_b) ? 0.5 : 0.0;
1470 for (int function_index = 0; function_index < number_of_functions; function_index++) {
1471 size_t column = static_cast<size_t>(function_index) * number_of_plane_waves;
1472 for (int g_index = 0; g_index < number_of_plane_waves; g_index++) {
1473 Real component_a = geometry.cartesian[g_index][static_cast<size_t>(direction_a)];
1474 Real component_b = geometry.cartesian[g_index][static_cast<size_t>(direction_b)];
1475 Real radial_weight =
1476 component_a * component_b * geometry.inverse_magnitude_bohr[g_index];
1477 strained[column + g_index] = -(radial[column + g_index] * radial_weight +
1478 angular[column + g_index] * component_b) -
1479 diagonal * base[column + g_index];
1480 }
1481 }
1482 }
1483
1484 /**
1485 * Add sum_ij q_ij left(m,i) right(j,v) into `target`, over every atom's own projector
1486 * block. S runs over every atom that carries projectors, Hubbard or not, so this sum
1487 * does too -- the same reason a Hubbard *force* appears on non-Hubbard atoms.
1488 */
1489 template <typename AugmentationCharge,
1490 typename OrbitalProjectorIndex,
1491 typename ProjectorBandIndex,
1492 typename OrbitalBandIndex>
1493 static void accumulate_augmentation(int number_of_atoms,
1494 const std::vector<int>& channels_by_atom,
1495 const std::vector<int>& projector_offset_by_atom,
1496 const AugmentationCharge& augmentation_charge,
1497 int total_hubbard_orbitals,
1498 int number_of_bands,
1499 const std::vector<Complex>& left,
1500 const std::vector<Complex>& right,
1501 const OrbitalProjectorIndex& orbital_projector,
1502 const ProjectorBandIndex& projector_band,
1503 const OrbitalBandIndex& orbital_band,
1504 std::vector<Complex>& target) {
1505 for (int atom_index = 0; atom_index < number_of_atoms; atom_index++) {
1506 int offset = projector_offset_by_atom[static_cast<size_t>(atom_index)];
1507 int number_of_channels = channels_by_atom[static_cast<size_t>(atom_index)];
1508 for (int channel_i = 0; channel_i < number_of_channels; channel_i++) {
1509 for (int channel_j = 0; channel_j < number_of_channels; channel_j++) {
1510 Real q_value = augmentation_charge(atom_index, channel_i, channel_j);
1511 if (std::abs(q_value) < math::tolerances::d_matrix_element_threshold)
1512 continue;
1513 int projector_i = offset + channel_i;
1514 int projector_j = offset + channel_j;
1515 for (int orbital = 0; orbital < total_hubbard_orbitals; orbital++) {
1516 Complex left_value = left[orbital_projector(orbital, projector_i)];
1517 if (left_value == Complex(0.0, 0.0))
1518 continue;
1519 for (int band = 0; band < number_of_bands; band++)
1520 target[orbital_band(orbital, band)] +=
1521 q_value * left_value * right[projector_band(projector_j, band)];
1522 }
1523 }
1524 }
1525 }
1526 }
1527
1528 /**
1529 * Rewrite a Gamma-only half-sphere ZGEMM result as its full-sphere value.
1530 *
1531 * Generalizes `convert_gamma_only_projector_overlaps_to_full_sphere`'s flag: the G = 0
1532 * term is read from the operands themselves rather than assumed present or absent. It
1533 * is neither here -- a strain derivative is -[(k+G)-weighted terms] - delta_ab A/2, so
1534 * it vanishes at G = 0 off the diagonal and does not on it, and hard-coding either rule
1535 * would be wrong for half the tensor.
1536 */
1537 static void complete_gamma_only_overlaps(int number_of_plane_waves,
1538 int number_of_left,
1539 int number_of_right,
1540 const Complex* left,
1541 const Complex* right,
1542 int g_zero_index,
1543 std::vector<Complex>& overlaps) {
1544 for (int left_index = 0; left_index < number_of_left; left_index++) {
1545 for (int right_index = 0; right_index < number_of_right; right_index++) {
1546 Real g_zero_product =
1547 g_zero_index >= 0
1548 ? (std::conj(left[static_cast<size_t>(left_index) * number_of_plane_waves +
1549 g_zero_index]) *
1550 right[static_cast<size_t>(right_index) * number_of_plane_waves +
1551 g_zero_index])
1552 .real()
1553 : 0.0;
1554 size_t flat_index = static_cast<size_t>(left_index) +
1555 static_cast<size_t>(right_index) * number_of_left;
1556 overlaps[flat_index] =
1557 gamma_only_full_sphere_overlap(overlaps[flat_index], g_zero_product);
1558 }
1559 }
1560 }
1561
1562 /**
1563 * Replace the tensor by its symmetric part.
1564 *
1565 * Averaged, not mirrored: the dY_lm/dq_a (k+G)_b term is genuinely asymmetric in
1566 * (a,b) before symmetrization, and the full 3x3 is accumulated above so that the
1567 * antisymmetric half is cancelled rather than discarded.
1568 */
1569 static void symmetrize_stress_tensor(Real stress[3][3]) {
1570 for (int direction_a = 0; direction_a < 3; direction_a++) {
1571 for (int direction_b = direction_a + 1; direction_b < 3; direction_b++) {
1572 Real average =
1573 0.5 * (stress[direction_a][direction_b] + stress[direction_b][direction_a]);
1574 stress[direction_a][direction_b] = average;
1575 stress[direction_b][direction_a] = average;
1576 }
1577 }
1578 }
1579
1580 const Crystal& crystal_structure_;
1581 const PlanewaveBasis& planewave_basis_;
1582 const KPointsGrid& kpoint_grid_;
1583 const NormConservingPlanewaveKohnShamSCFState& scf_state_;
1584 StressTensorContributions& result_;
1585};
1586
1587} // namespace q3
atomic_orbital_builder.hpp
q3::Crystal::cell_volume
Real cell_volume() const
Definition crystal.hpp:193
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::HubbardStressTensorCalculator::supports
static bool supports(const NormConservingPlanewaveKohnShamSCFState &scf_state)
Definition hubbard.hpp:119
q3::HubbardStressTensorCalculator::HubbardStressTensorCalculator
HubbardStressTensorCalculator(const Crystal &crystal_structure, const PlanewaveBasis &planewave_basis, const KPointsGrid &kpoint_grid, const NormConservingPlanewaveKohnShamSCFState &scf_state, StressTensorContributions &result)
Definition hubbard.hpp:98
q3::KPointsGrid::number_of_kpoints
int number_of_kpoints() const
Definition kpoints_grid.hpp:71
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::is_gamma_only
bool is_gamma_only() const
Definition init.hpp:157
q3::UltrasoftPlanewaveKohnShamQuantumSystem::ultrasoft_basis
const UltrasoftPlanewaveBasis & ultrasoft_basis() const
Definition init.hpp:59
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
factory.hpp
hubbard.hpp
init.hpp
lowdin_orthogonalization.hpp
q3::constants::two_pi
constexpr Real two_pi
Definition constants.hpp:25
q3::lattice_fourier_transform::internal::phase
Complex phase(const Vector3D &kpoint_crystal, const IntegerVector3D &translation, Real sign)
e^{sign * 2 pi i k.R} for crystal-coordinate k and integer R.
Definition lattice_fourier_transform.hpp:48
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::math::tolerances::bessel_argument_zero_threshold
constexpr double bessel_argument_zero_threshold
Definition tolerances.hpp:158
q3::paw_one_center::radial_derivative
std::vector< Real > radial_derivative(const std::vector< Real > &values, const std::vector< Real > &measures, int mesh_size)
Definition gradient.hpp:45
q3::compute_real_spherical_harmonics_batch
void compute_real_spherical_harmonics_batch(int number_of_lm_components, int number_of_g_vectors, const std::vector< std::array< Real, 3 > > &g_vectors, const std::vector< Real > &g_squared, std::vector< std::vector< Real > > &spherical_harmonics_values)
Definition harmonics.hpp:218
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::harmonic_combined_degree_order_index
int harmonic_combined_degree_order_index(int degree, int order)
Definition harmonics.hpp:61
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_spherical_harmonic_derivatives_batch
void compute_spherical_harmonic_derivatives_batch(int number_of_lm_components, int number_of_g_vectors, const std::vector< std::array< Real, 3 > > &g_vectors, const std::vector< Real > &g_squared, int polarization_index, std::vector< std::vector< Real > > &spherical_harmonic_derivatives)
Definition harmonics.hpp:353
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::compute_phase_factor
Complex compute_phase_factor(Real phase)
Compute e^{-iφ} = cos(φ) - i·sin(φ) (standard structure factor phase).
Definition phase_factor.hpp:17
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::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::compute_angular_phase_factor
Complex compute_angular_phase_factor(int angular_momentum)
Definition phase_factor.hpp:35
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::lagrange_interpolation_4point_derivative
Real lagrange_interpolation_4point_derivative(const std::vector< Real > &table, Real q_value, Real q_spacing, int table_size)
Definition lagrange.hpp:60
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
q3::lagrange_interpolation_4point
Real lagrange_interpolation_4point(const std::vector< Real > &table, Real q_value, Real q_spacing, int table_size)
Definition lagrange.hpp:31
q3::harmonic_order_sequence
std::vector< int > harmonic_order_sequence(int degree)
Definition harmonics.hpp:40
phase_factor.hpp
wavefunction_atomic.hpp
contributions.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
q3::StressTensorContributions
Per-contribution stress tensor breakdown (kbar after conversion).
Definition contributions.hpp:33