Update PDOS fallback, add EXX/RPA support

ChangeLog

@@ -2,6 +2,14 @@
 -Date
 -Name
 -changes
+--------------
+November 19, 2025
+Name: Qihao Cheng
+Changes: utils/pdos/, utils/atom/
+1. Added XC functional fallback in PDOS: unsupported SPARC functionals now default to GGA_PBE for atomic orbital generation.
+2. Extended atomic DFT module to support EXX and RPA (OEP-based implementation).
+3. Added thread-parallel RPA frequency evaluation and improved standalone atomic code performance.
+
 --------------
 November 07, 2025
 Name: Qihao Cheng

src/initialization.c

@@ -3722,7 +3722,7 @@ void write_output_init(SPARC_OBJ *pSPARC) {
     }
 
     fprintf(output_fp,"***************************************************************************\n");
-    fprintf(output_fp,"*                       SPARC (version November 07, 2025)                 *\n");
+    fprintf(output_fp,"*                       SPARC (version November 19, 2025)                 *\n");
     fprintf(output_fp,"*   Copyright (c) 2020 Material Physics & Mechanics Group, Georgia Tech   *\n");
     fprintf(output_fp,"*           Distributed under GNU General Public License 3 (GPL)          *\n");
     fprintf(output_fp,"*                   Start time: %s                  *\n",c_time_str);

utils/atom/__init__.py

@@ -15,6 +15,59 @@
 
 __version__ = "0.1.0"
 
+
+import sys
+import os
+
+# ---------------------------------------------------------------------------
+# Record whether NumPy was already imported BEFORE this package was imported.
+# This is the key information for determining BLAS thread safety later.
+# ---------------------------------------------------------------------------
+_NUMPY_IMPORTED_BEFORE_ATOMIC = ("numpy" in sys.modules)
+
+
+# ---------------------------------------------------------------------------
+# Record whether BLAS environment variables already indicate single-threaded
+# execution. This allows us to avoid warnings even when NumPy was imported
+# earlier, because the user may have correctly configured BLAS beforehand.
+# ---------------------------------------------------------------------------
+
+def _env_says_single_thread() -> bool:
+    for var in ("OMP_NUM_THREADS", "MKL_NUM_THREADS",
+                "OPENBLAS_NUM_THREADS", "NUMEXPR_NUM_THREADS"):
+        val = os.environ.get(var)
+        if val is None:
+            return False
+        if val.strip() not in ("1", "1.0"):
+            return False
+    return True
+
+_BLAS_ENV_SINGLE_THREADED = _env_says_single_thread()
+
+# ---------------------------------------------------------------------------
+# Record whether threadpoolctl is installed.
+# ---------------------------------------------------------------------------
+def _threadpoolctl_installed() -> bool:
+    try:
+        import threadpoolctl
+        return True
+    except ImportError:
+        return False
+
+_THREADPOOLCTL_INSTALLED = _threadpoolctl_installed()
+
+# Force MKL / OpenBLAS into single-thread mode without any extra dependency.
+if not _BLAS_ENV_SINGLE_THREADED:
+    os.environ["OMP_NUM_THREADS"] = "1"
+    os.environ["MKL_NUM_THREADS"] = "1"
+    os.environ["OPENBLAS_NUM_THREADS"] = "1"
+    os.environ["NUMEXPR_NUM_THREADS"] = "1"
+else:
+    # Do not override user settings when NumPy is already loaded.
+    pass
+
+
+
 # Main solver class
 from .solver import AtomicDFTSolver
 

utils/atom/mesh/builder.py

@@ -131,11 +131,11 @@ def gauss_legendre(n: int) -> Tuple[np.ndarray, np.ndarray]:
         assert n > 0, \
             QUADRATURE_POINT_NUMBER_NOT_GREATER_THAN_0_ERROR.format(n)
 
-        beta = 0.5/np.sqrt(1 - 1 / (2*(np.arange(1,n)))**(2))
-        T = np.diag(beta,k=1) + np.diag(beta,k=-1)
-        nodes, eigvecs = np.linalg.eigh(T,UPLO='L')
+        beta = 0.5 / np.sqrt(1.0 - 1.0 / (2.0 * (np.arange(1, n)))**2)
+        T = np.diag(beta, k=1) + np.diag(beta, k=-1)
+        nodes, eigvecs = np.linalg.eigh(T, UPLO='L')
         index = np.arange(nodes.size)
-        legendre_weights = 2*(eigvecs[0,index])**2
+        legendre_weights = 2*(eigvecs[0, index])**2
         return nodes, legendre_weights
 
 

utils/atom/pseudo/local.py

@@ -185,7 +185,7 @@ def get_rho_core_correction(self, r_nodes: np.ndarray) -> np.ndarray:
                 self.density_nlcc, 
                 k=3
             )
-            rho_nlcc[indices_below_cutoff[0]] = rho_nlcc_interpolator(r_below_cutoff)
+            rho_nlcc[indices_below_cutoff[0]] = rho_nlcc_interpolator(r_below_cutoff).reshape(-1)
 
             return rho_nlcc
 
@@ -318,8 +318,8 @@ def print_info(self):
 
 
         # Basic atomic information
-        print(f"Valence Charge               (z_valence)             : {self.z_valence}")
-        print(f"Nuclear Charge               (z_nuclear)             : {self.z_nuclear}")
+        print(f"  Valence Charge             (z_valence)             : {self.z_valence}")
+        print(f"  Nuclear Charge             (z_nuclear)             : {self.z_nuclear}")
         print()
 
         if self.load_psp:

utils/atom/scf/convergence.py

@@ -10,6 +10,8 @@
 # Convergence Checker Error Messages
 LOOP_TYPE_ERROR_MESSAGE = \
     "parameter loop_type must be 'Inner' or 'Outer', get {} instead"
+FOOTER_BLANK_LINE_TYPE_ERROR_MESSAGE = \
+    "parameter footer_blank_line must be a boolean, get type {} instead"
 
 @dataclass
 class ConvergenceHistory:
@@ -36,11 +38,12 @@ class ConvergenceChecker:
 
     def __init__(
         self,
-        tolerance: float = 1e-6,
-        n_consecutive: int = 1,
-        norm_type: int = 2,
-        loop_type: str = "Inner"
-    ):
+        tolerance        : float = 1e-6,
+        n_consecutive    : int   = 1,
+        norm_type        : int   = 2,
+        loop_type        : str   = "Inner",
+        footer_blank_line: bool  = False,
+        ):
         """
         Parameters
         ----------
@@ -52,17 +55,28 @@ def __init__(
             Norm type for residual calculation (1, 2, or np.inf)
         loop_type : str
             Type of loop for display purposes ("Inner" or "Outer")
+        footer_blank_line: bool
+            If True, print a blank line after the footer
         """
         assert loop_type in ["Inner", "Outer"], \
             LOOP_TYPE_ERROR_MESSAGE.format(loop_type)
 
-        self.tolerance     = tolerance
-        self.n_consecutive = n_consecutive
-        self.norm_type     = norm_type
-        self.loop_type     = loop_type
-
+        self.tolerance          = tolerance
+        self.n_consecutive      = n_consecutive
+        self.norm_type          = norm_type
+        self.loop_type          = loop_type
+        self._footer_blank_line = footer_blank_line
+
         self._consecutive_count = 0
         self._history = ConvergenceHistory(iterations=[], residuals=[])
+
+    def set_footer_blank_line(self, footer_blank_line: bool):
+        """
+        Set the footer blank line flag
+        """
+        assert isinstance(footer_blank_line, bool), \
+            FOOTER_BLANK_LINE_TYPE_ERROR_MESSAGE.format(type(footer_blank_line))
+        self._footer_blank_line = footer_blank_line
 
 
     def check(
@@ -170,8 +184,10 @@ def print_footer(self, converged: bool, n_iterations: int, prefix: str = ""):
         """Print table footer with final status"""
         status_msg = "converged" if converged else "not converged"
         if self.loop_type == "Outer":
-            print(f"[Outer] Converged after {n_iterations} iteration(s)")
+            print(f"{prefix} Outer SCF iteration converged after {n_iterations} iteration(s)")
         else:
             print(f"{prefix}\t SCF Iteration {status_msg} after {n_iterations} iteration(s)")
+
+        if self._footer_blank_line:
             print()