Pao (#180)

Correct errors in projected atomic orbitals

- do not othogonalise before norm calculation
- use absolute of eigenvalues for cutoff

Author: MIWdlB

docs/jupyter_execute/notebooks/localization.ipynb

@@ -9,15 +9,14 @@
 class LocalizedSystem:
     """Required data from localized system.
 
-    Args:
-        active_mo_inds (np.array): 1D array of active occupied MO indices.
-        enviro_mo_inds (np.array): 1D array of environment occupied MO indices.
-        c_active (np.array): C matrix of localized occupied active MOs (columns define MOs).
-        c_enviro (np.array): C matrix of localized occupied ennironment MOs.
-        c_loc_occ (np.array): C matrix of localized occupied MOs.
-        c_loc_virt (np.array | None): C matrix of localized virual MOs.
-        dm_active (np.array): active system density matrix.
-        dm_enviro (np.array): environment system density matrix.
+    active_mo_inds (np.array): 1D array of active occupied MO indices
+    enviro_mo_inds (np.array): 1D array of environment occupied MO indices
+    c_active (np.array): C matrix of localized occupied active MOs (columns define MOs)
+    c_enviro (np.array): C matrix of localized occupied ennironment MOs
+    c_loc_occ (np.array): C matrix of localized occupied MOs
+    c_loc_virt (np.array | None): C matrix of localized virual MOs.
+    dm_active (np.array): active system density matrix
+    dm_enviro (np.array): environment system density matrix
     """
 
     active_mo_inds: NDArray

nbed/localizers/system.py

@@ -5,7 +5,6 @@
 import numpy as np
 from numpy.typing import NDArray
 from pyscf.lib import StreamObject
-from scipy.linalg import fractional_matrix_power
 
 from nbed.localizers.virtual.base import VirtualLocalizer
 
@@ -106,9 +105,7 @@ def _localize_virtual_spin_pao(
     # Take the columns of C matrix (MOs)
     truncated_paos = pao_projector[:, np.abs(pao_norms) > norm_cutoff]
 
-    s_half = fractional_matrix_power(ao_overlap, 0.5)
-
-    renormalized_paos = s_half @ truncated_paos
+    renormalized_paos = truncated_paos
     renormalized_paos = renormalized_paos / np.sqrt(
         np.einsum("ij,ij->j", renormalized_paos, renormalized_paos)
     )
@@ -121,7 +118,7 @@ def _localize_virtual_spin_pao(
 
     logger.debug(f"{eigvecs.shape=}")
     # How to transform the truncated paos?
-    final_paos = renormalized_paos[:, eigvals > overlap_cutoff]
+    final_paos = renormalized_paos[:, np.abs(eigvals) > overlap_cutoff]
     logger.debug(f"{final_paos.shape=}")
 
     if (n_paos := final_paos.shape[-1]) == 0: