Notes: projected orbital selection in QO calculated from LCAO may be highly nontrivial

[kirk0830]

Describe the Code Quality Issue

Related issues

Connected with reported problem of Quaisatomic Orbital (QO) implementation: #3632 #3484

Theoretical background

QO method is for finding subspace from a space spanned by states calculated from DFT which is relatively large. It is expected that QO can extract subspace composed of bands of interest. However, the projection from the large space to the subspace might be a non-trivial task.

Notation

General

• $H{\mu\nu}(\mathbf{k})$: element of Hamiltonian matrix in representation of ABACUS numerical atomic orbitals (NAOs) at $\mathbf{k}$-point: $H{\mu\nu}(\mathbf{k}) =\langle\phi{\mu}(\mathbf{k}) |\hat{H} |\phi{\nu}(\mathbf{k})\rangle$
• $\phi_{\mu}(\mathbf{k})$: numerical atomic orbital (NAO) at $\mathbf{k}$-point
• $\mu,\nu$: indices of NAOs
• $\mathbf{k}$: wave vector, $\mathbf{k} = (k_x, k_y, k_z)$, representing one certain kind of translational symmetry
• $\hat{H}$: Hamiltonian operator
• $\hat{P}$: projection operator
• $\hat{P} =\sum{\mu\nu} |\phi{\mu}\rangle\langle\phi{\nu} |$
  Similarly the overlap matrix $S{\mu\nu}(\mathbf{k}) =\langle\phi{\mu}(\mathbf{k}) |\phi{\nu}(\mathbf{k})\rangle$ is defined.
• $|\psi_n(\mathbf{k})\rangle$: Bloch state at $\mathbf{k}$-point
• $n$: band index
• $E_n(\mathbf{k})$: energy of Bloch state $|\psi_n(\mathbf{k})\rangle$

QO-specific

Also should define things related to the subspace.

• $|A\mu\rangle$: Atom-centered orbital that should be constructed manually. More specifically, the atom-centered orbital in real space in cell $\mathbf{R}$ is $|A{\mu}(\mathbf{R})\rangle$, and the corresponding reciprocal space representation is $|A{\mu}(\mathbf{k})\rangle =\frac{1}{\sqrt{N\text{cell}}}\sum{\mathbf{R}} e^{i\mathbf{k}\cdot\mathbf{R}} |A{\mu}(\mathbf{R})\rangle$.
• $|Q{\mu}\rangle$: Quasiatomic orbital (QO) constructed from atom-centered orbitals. Correspondingly, $|Q{\mu}(\mathbf{k})\rangle =\frac{1}{\sqrt{N\text{cell}}}\sum{\mathbf{R}} e^{i\mathbf{k}\cdot\mathbf{R}} |Q_{\mu}(\mathbf{R})\rangle$
• $\hat{P}_{\text{QO}}$: projection operator to the subspace spanned by QOs Can also define Hamiltonian and overlap matrix in QO representation, in k-space and R counterparts.
• $H{\mu\nu}^{\text{QO}}(\mathbf{k}) =\langle Q{\mu}(\mathbf{k}) |\hat{H} | Q_{\nu}(\mathbf{k})\rangle$
• $S{\mu\nu}^{\text{QO}}(\mathbf{k}) =\langle Q{\mu}(\mathbf{k}) | Q_{\nu}(\mathbf{k})\rangle$
  Additionally the overlap between NAO and atom-centered orbital is,
• $S{\mu\nu}^{\phi A}(\mathbf{k}) =\langle\phi{\mu}(\mathbf{k}) | A_{\nu}(\mathbf{k})\rangle$

In ABACUS, the Hamiltonian and overlap matrix in NAO representation is calculated and output to file, along with the NAO-atom-centered orbital overlap matrix. In Python code, QOs are constructed and Hamilontian, overlap matrices are calculated.

Derivation

QO is defined to be calculated in the following way:
$|Q{\mu}(\mathbf{R})\rangle=\sum{\mathbf{k}}{\sum{n}{|\varphi{n}(\mathbf{k})\rangle\langle\varphi{n}(\mathbf{k})|A{\mu}(\mathbf{R})\rangle}}$ where $|A{\mu}(\mathbf{R})\rangle$ is atom-centered orbital in real space. The overlap between QO in realspace is obviously to be: $S{\mu\nu}^{\mathrm{QO}}(\mathbf{R})=\sum_{\mathbf{k}}{\sumn{\langle A{\mu}(\mathbf{0})|\varphi_n(\mathbf{k})\rangle\langle\varphin(\mathbf{k})|A{\nu}(\mathbf{R})\rangle}}$ , of which the physical meaning is clear: the overlap between QOs in real space is the that of atom-centered orbital represented in space of ${\varphi_n(\mathbf{k})}$. It is the same for Hamiltonian matrix in QO representation but this time write in matrix form:
$\mathbf{H}^{QO}\left(\mathbf{R}\right)=\left[\mathbf{D}\left(\mathbf{k}\right)\mathbf{S}^{\phi A}\left(\mathbf{k}\right)\right]^{\dagger}\mathbf{H}\left(\mathbf{k}\right)\left[\mathbf{D}\left(\mathbf{k}\right)\mathbf{S}^{\phi A}\left(\mathbf{k}\right)\right]$ where $\mathbf{D}\left(\mathbf{k}\right)$ is density matrix in ${\varphi_n(\mathbf{k})}$ space, and $\mathbf{S}^{\phi A}\left(\mathbf{k}\right)$ is overlap matrix between NAO and atom-centered orbital.
However ${\varphi_n(\mathbf{k})}$ and ${\psi_n(\mathbf{k})}$ are not the same. The latter is subspace of the former, the rest of the former is constructed from atom-centered orbitals, which more specifically speaking are those cannot be represented by ${\psin(\mathbf{k})}$, and treated as empty states introduced. Keep in mind there must be orthogonality between solved occupied states and introduced empty states, thus to get the correct set of empty states, from atom-centered orbitals, one should remove the occupied states from the set of atom-centered orbitals: $|A{\mu}^{\text{empty}}\rangle = |A{\mu}\rangle-\sum{n}{|\varphi{n}\rangle\langle\varphi{n} | A_{\mu}\rangle}$ Orthogonality should also be kept between introduced empty states, this is guranateed by explicit diagonalization of overlap matrix between empty states.
Remember the dimension of QO (or say the subspace) is the same as the number of atom-centered orbitals, therefore the number of empty states ($Nc$) plus dimension of manually selected subspace spanned by ${\psi_n(\mathbf{k})}$ ($Nb$) should be the same as the number of atom-centered orbitals ($Nq$), thus only $Nq-Nb=Nc$ empty states are needed to be introduced, which in actual implementation only the $Nc$ eigenvectors with largest eigenvalues are kept.
Combining ${\psi_n(\mathbf{k})}$ with empty states, the set of ${\varphi_n(\mathbf{k})}$ is complete, and the density matrix $\mathbf{D}\left(\mathbf{k}\right)$ can be constructed.

Problem

Basic

In ABACUS, two different ways to construct atomic orbitals are implemented: pswfc and hydrogen. The former use pseudo-wavefunction from pseudopotential file, while the latter directly construct hydrogen-like (single-electron) raidal function with Laguerre Polynomials. The band reproduction result is poor (blue circles are eigenvalues solved from eigenequation in rep. of QO, lines are result of LCAO):
(pswfc, use all pswfcs for projection)
At all other kpoints, diagonalization fails due to non-positive-definite of QO overlap.
Change kpoint mesh from 999 to 555 or 777:
Change pswfc to hydrogen, the strategy for generating orbitals set to minimal-nodeless, for Si will produce 1s, 2p and 3d orbitals (nodeless radial functions, generate to nmax the principle quantum number):
Change to minimal-valence, for Si generate 3s, 3p and 3d orbitals:
Change to energy-valence, for Si generate 3s and 3p orbitals:

Update

• change the method for finding eigenvals and eigenvecs of overlap between empty states to SVD
• reduce the dimension of matrices in QO rep. according to one manually set threshold (canonical diagonalization)
  For pswfc:
  For hydrogen with minimal-nodeless:
  For hydrogen with minimal-valence:
  For hydrogen with energy-valence:
  

Additional Context

No response

Task list for Issue attackers (only for developers)

• [X] Identify the specific code file or section with the code quality issue.
• [X] Investigate the issue and determine the root cause.
• [ ] Research best practices and potential solutions for the identified issue.
• [ ] Refactor the code to improve code quality, following the suggested solution.
• [ ] Ensure the refactored code adheres to the project's coding standards.
• [ ] Test the refactored code to ensure it functions as expected.
• [ ] Update any relevant documentation, if necessary.
• [ ] Submit a pull request with the refactored code and a description of the changes made.

[kirk0830]

With single zeta orbitals, I get this result: Seems the space spanned by numerical atomic orbital is small and any other constructed orbital would have significant part of space spanned outside the space of numerical atomic orbitals.

[kirk0830]

I will submit a pull request to update ABACUS QO module and close this issue temporarily. cc. @QG-phys