- Generated by
1.9.8
|
Quantum 3xplorer (Q3)
An AI-native quantum simulation engine
|
#include <initializer.hpp>
Public Member Functions | |
| void | compute_atomic_wavefunctions (int kpoint_index, int number_of_plane_waves, int number_of_bands, const std::vector< int > &g_vector_indices, const Crystal &crystal_structure, const KPointsGrid &kpoint_grid, const PlanewaveBasis &planewave_basis, const PseudopotentialSet &pseudopotential_set, std::vector< std::vector< Complex > > &atomic_wavefunctions) |
| void | initialize_interpolation_table (const PseudopotentialSet &pseudopotential_set, Real cell_volume) |
| void | rotate_wavefunctions (int kpoint_index, int kpoint_index_state, int number_of_plane_waves, int number_of_bands, const std::vector< int > &g_vector_indices, const std::vector< Real > &kinetic_energy_squared, const std::vector< Real > &smooth_potential, const Crystal &crystal_structure, const KPointsGrid &kpoint_grid, const PlanewaveBasis &planewave_basis, const execution::PlanewaveExecutionPolicy &execution_policy, const NormConservingPlanewaveKohnShamQuantumSystem &pseudopotential_quantum_system, const std::vector< std::vector< Complex > > &atomic_wavefunctions, PlanewaveKohnShamQuantumSystem &quantum_system, const Matrix< Complex > &precomputed_kleinman_bylander_potential) |
Computes starting wavefunctions from atomic orbitals (chi) in the pseudopotential.
The procedure is the standard LCAO starting guess for a plane-wave SCF (Payne et al., Rev. Mod. Phys. 64, 1045 (1992), sec. V): superpose the pseudo-atomic orbitals of every atom, pad the set to the requested band count, then take the Rayleigh-Ritz solution of H c = e S c on that subspace. Rayleigh-Ritz is used because it is the variationally optimal selection of number_of_bands vectors out of a larger non-orthogonal subspace, and because it delivers them S-orthonormal, which is what the Davidson solver requires on entry.
Quantum ESPRESSO's pw.x reaches the same three stages. Its routine names are not reproduced here and its stage boundaries are not the ones below: its driver also performs wavefunction file I/O and Hubbard-projector orthogonalisation, neither of which this class does, and the per-k-point driver loop lives in the SCF iterator rather than here.
Definition at line 71 of file initializer.hpp.
|
inline |
Compute the starting vector set in G-space for one k-point. Returns atomic_wavefunctions[starting_vector_count][number_of_plane_waves], where starting_vector_count = max(atomic orbital count, band count).
Two steps: the pure atomic orbitals come from the shared builder (see atomic_orbital_builder.hpp for their derivation), and this method then brings the set up to the band count – padding the tail with random vectors when the atomic orbitals do not reach it.
Build the pure atomic orbital set (no padding) via the shared builder – identical construction used by PDOS/Löwdin projection (layer_3_representation/pseudopotential/planewave/internal/atomic_orbital_builder.hpp).
Symmetry-breaking hook, currently INERT: the coefficient below is 0.0, so the loop that follows multiplies every coefficient by exactly 1 and no perturbation is applied. It is kept as the place a perturbation would go, not as one that happens. See the audit note (R-D5-1) – either the branch goes, or a magnitude criterion gets written down.
The branch exists at all because a perturbation is only useful when the atomic orbitals already cover the band count: when they do not, the random padding below breaks the same degeneracy for free.
Fill the missing vectors with random coefficients. Each is drawn from a hash of the k-point, the vector and the plane wave's Miller indices, not from a generator stepped along the plane-wave list: ranks sharing a pool each hold a different part of that list, and k-point pools number their k-points locally, so a stepped generator gave every topology its own starting guess. The ultrasoft silicon dataset carries no atomic orbitals and starts from these vectors alone, and took a different SCF path on every rank count to the same fixed point.
Definition at line 110 of file initializer.hpp.
References q3::build_atomic_orbitals_g_space(), q3::PlanewaveBasis::g_vectors(), and q3::constants::two_pi.
Referenced by q3::PlanewavePseudopotentialKohnShamSCFIterator::initialize_single_particle_states().
|
inline |
Precompute the radial Fourier transform interpolation table for all species:
table[species][shell][q] = 4π/√Ω ∫ chi(r) j_l(q·r) r dr
The 4π and the j_l come from the Rayleigh expansion of a plane wave, the 1/√Ω from the plane-wave normalisation, and the single power of r from UPF storing r·R_nl(r) rather than R_nl(r).
The q axis is innermost because the only consumer is a 4-point Lagrange stencil that reads four consecutive q entries, and lagrange_interpolation_4point takes a one-dimensional q table. Species is outermost because the shell count varies per species, so the jagged axis has to be the inner of the two; that leaves the shell axis with no remaining freedom.
Definition at line 88 of file initializer.hpp.
References q3::math::tolerances::default_number_of_q_points_for_atomic_wavefunctions, q3::log::info(), and q3::tabulate_atomic_wavefunctions_in_q_space().
Referenced by q3::PlanewavePseudopotentialKohnShamSCFIterator::initialize_single_particle_states().
|
inline |
Rayleigh-Ritz on the starting subspace: form H_ij and S_ij in the basis of the starting vectors, solve the generalised eigenproblem H c = e S c, and rotate. Taking the lowest number_of_bands Ritz pairs is the variationally optimal selection out of that subspace, and the Ritz vectors come out S-orthonormal, which is the eigensolver's entry condition.
Gamma-only takes the real path throughout: under psi(-G) = conj(psi(G)) the subspace matrices are real symmetric, so a real solver is both cheaper and the one that keeps the rotated basis gamma-symmetric.
This overload accepts pre-computed KB projectors to avoid redundant recomputation.
smooth_potential is the local potential on this rank's part of the wavefunction grid, already multiplied by 1/N of that grid – what the Hamiltonian application holds, not the density-grid potential the SCF stores.
Pack atomic wavefunctions into contiguous column-major matrix: contiguous_wavefunctions[g + wfc * npw] = atomic_wavefunctions[wfc][g]
Apply S to every starting vector. For norm-conserving pseudopotentials S = I, so the overlap-applied set IS the input set and nothing is allocated: everything downstream reads it, so pointing at the input is enough. Ultrasoft needs the q_ij augmentation and so needs its own buffer. This block is the tallest point of a run's setup, and memory released here is not returned to the system, so the array not taken is worth more than its lifetime suggests.
Gamma-only: use DGEMM on interleaved real/imaginary parts. Reinterpret complex arrays as double arrays of length 2*npw. DGEMM("T","N", n, n, 2*npw, alpha=2) computes 2*Re(Σ conj(a)*b) which is the gamma inner product before G=0 correction.
Non-gamma: use ZGEMM H_ij = χ†·Hχ: ZGEMM('C', 'N', n_wfc, n_wfc, npw)
Reduce H_ij and S_ij across band ranks within this pool. Each band rank holds the same G-vectors (replicated), but in the future when G-space is distributed, each rank holds a slab and we need reduce_within_pool_vector to sum the partial dot products.
Diagonalize: H*c = e*S*c For gamma-only, use dsygv (real symmetric) to guarantee real eigenvectors that preserve gamma symmetry after ZGEMM rotation.
Rotate wavefunctions: evc = χ · eigenvectors via ZGEMM evc(npw × n_bands) = χ(npw × n_wfc) · c(n_wfc × n_bands) The eigenvectors are stored column-major in hamiltonian_matrix after zhegv.
Gamma-only: enforce real G=0 for the rotated wavefunctions, which the two-band FFT trick requires.
The atomic orbitals are already real at G=0 here – at gamma the k-phase and the structure factor are both unity, and chi_l(0) = 0 for l >= 1 because j_l(0) = 0, so only the l = 0 shells contribute and they carry i^0 = 1. The complex G=0 component this removes comes from the random padding vectors, which are generated with an arbitrary phase. Because the subspace matrices above were formed from the uncleaned vectors, the G=0 correction there uses the general Re[conj(a0)*b0] form rather than the real-parts-only form that would need this cleanup to run first.
Definition at line 225 of file initializer.hpp.
References q3::Matrix< T >::data(), q3::blas::dgemm_(), q3::lapack::dsygv_(), q3::log::error(), q3::PlanewaveBasis::gamma_origin_is_local(), q3::NormConservingPlanewaveKohnShamQuantumSystem::has_ultrasoft(), q3::PlanewaveBasis::is_gamma_only(), q3::execution::PlanewaveExecutionPolicy::reduce_within_pool_vector(), q3::PlanewaveKohnShamQuantumSystem::set_eigenvalue(), q3::PlanewaveKohnShamQuantumSystem::wavefunction_band_data(), q3::blas::zgemm_(), and q3::lapack::zhegv_().
Referenced by q3::PlanewavePseudopotentialKohnShamSCFIterator::initialize_single_particle_states().
1.9.8