Added the module implementing Lowdin orthogonalization, added this correction to the cp2k workflow

Author: alexvakimov

src/libra_py/__init__.py

@@ -48,6 +48,7 @@
            "namd",
            "normal_modes",
            "nve_md",
+           "orthogonalizations",
            "packages",
            "parse_gamess",
            "pdos",

src/libra_py/orthogonalizations.py

@@ -0,0 +1,97 @@
+# *********************************************************************************
+# * Copyright (C) 2026 Alexey V. Akimov
+# *
+# * This file is distributed under the terms of the GNU General Public License
+# * as published by the Free Software Foundation, either version 3 of
+# * the License, or (at your option) any later version.
+# * See the file LICENSE in the root directory of this distribution
+# * or <http://www.gnu.org/licenses/>.
+# *
+# *********************************************************************************/
+
+"""
+.. module:: orthogonalizations
+   :platform: Unix, Windows
+   :synopsis: this module implements various general-purpose orthogonalization procedures
+
+.. moduleauthor:: Alexey V. Akimov, ChatGPT
+
+"""
+
+import numpy as np
+
+def lowdin_inverse_sqrt(S, thresh=1e-10):
+    """
+    Compute the inverse square root of an overlap matrix using Löwdin symmetric orthogonalization.
+
+    This function evaluates S^{-1/2} via eigenvalue decomposition:
+        S = U diag(λ) U†  →  S^{-1/2} = U diag(λ^{-1/2}) U†
+
+    Small eigenvalues below a specified threshold are treated as zero to ensure
+    numerical stability and avoid divergences. The input matrix is explicitly
+    symmetrized to enforce Hermiticity before diagonalization.
+
+    Parameters
+    ----------
+    S : array_like (N, N)
+        Real or complex overlap matrix. It is assumed to be Hermitian
+        (or close to Hermitian within numerical precision).
+    thresh : float, optional
+        Eigenvalue cutoff below which eigenvalues are treated as zero.
+        This prevents instabilities due to near-linear dependencies.
+        Default is 1e-10.
+
+    Returns
+    -------
+    S_inv_sqrt : ndarray (N, N)
+        The inverse square root of the overlap matrix, S^{-1/2}, computed
+        in the Löwdin (symmetric) orthogonalization scheme.
+
+    Notes
+    -----
+    - The function enforces Hermiticity via (S + S†)/2 before diagonalization.
+    - Input arrays are cast to float64 or complex128 to ensure compatibility
+      with NumPy linear algebra routines.
+    - Eigenvalues below `thresh` are set to zero, effectively projecting out
+      linearly dependent components of the basis.
+    - The resulting matrix can be used to orthonormalize a non-orthogonal basis:
+          C_orth = S^{-1/2} C
+
+    Raises
+    ------
+    LinAlgError
+        If the eigenvalue decomposition fails.
+
+    Examples
+    --------
+    >>> S = np.array([[1.0, 0.2], [0.2, 1.0]])
+    >>> S_inv_sqrt = lowdin_inverse_sqrt(S)
+    >>> np.allclose(S_inv_sqrt @ S @ S_inv_sqrt, np.eye(2))
+    True
+    """
+
+    # Ensure numpy array
+    S = np.array(S)
+
+    # Force supported dtype
+    if np.iscomplexobj(S):
+        S = S.astype(np.complex128)
+    else:
+        S = S.astype(np.float64)
+
+    # Enforce Hermiticity 
+    S = 0.5 * (S + S.conj().T)
+
+    eigvals, eigvecs = np.linalg.eigh(S)
+
+    eigvals_inv_sqrt = np.zeros_like(eigvals)
+
+    for i, val in enumerate(eigvals):
+        if val > thresh:
+            eigvals_inv_sqrt[i] = 1.0 / np.sqrt(val)
+        else:
+            eigvals_inv_sqrt[i] = 0.0
+
+    S_inv_sqrt = eigvecs @ np.diag(eigvals_inv_sqrt) @ eigvecs.conj().T
+
+    return S_inv_sqrt

src/libra_py/packages/cp2k/methods.py

@@ -33,6 +33,7 @@
 
 from libra_py import data_outs, data_stat, data_conv
 from libra_py import units, molden_methods
+import libra_py.orthogonalizations as ortho
 
 
 def ndigits(integer_number: int):
@@ -2867,6 +2868,7 @@ def cp2k_compute_adi(q, params, full_id):
     params.setdefault("MO_prev", {})
     params.setdefault("data_prev", {})
     params.setdefault("coordinates_prev", {})
+    params.setdefault("s_ci_inv_prev", {})
     params.setdefault("is_first_time", {})
     params.setdefault("act_state", {})
 
@@ -2918,6 +2920,7 @@ def cp2k_compute_adi(q, params, full_id):
     obj.hvib_adi = CMATRIX(nstates, nstates)
     obj.basis_transform = CMATRIX(nstates, nstates)
     obj.time_overlap_adi = CMATRIX(nstates, nstates)
+    obj.overlap_adi = CMATRIX(nstates, nstates)
 
     # ================= Compute overlaps =================
     st_mo = sp.load_npz(params["time_overlap_filename"])
@@ -2927,6 +2930,9 @@ def cp2k_compute_adi(q, params, full_id):
     else:
         st_mo = np.asarray(st_mo)
     st_mo = st_mo.astype(np.float64, copy=False)
+
+    s_mo = sp.load_npz(params["overlap_filename"]).todense()
+
 
     ndim = st_mo.shape[0] // 2
 
@@ -2964,6 +2970,20 @@ def cp2k_compute_adi(q, params, full_id):
 
     st_ci = ci.overlap(st_mo, data_prev, data_curr, ovlp_params)
 
+    s_ci = ci.overlap(s_mo, data_curr, data_curr, ovlp_params)
+    
+    s_ci_inv_curr = ortho.lowdin_inverse_sqrt(s_ci)
+    s_ci_inv_prev = None
+    if is_first_time:
+        s_ci_inv_prev = copy.deepcopy(s_ci_inv_curr)
+    else:
+        s_ci_inv_prev = params["s_ci_inv_prev"].get(itraj, s_ci_inv_curr)
+    
+    s_ci = s_ci_inv_curr @ s_ci @ s_ci_inv_curr
+    st_ci = s_ci_inv_prev @ st_ci @ s_ci_inv_curr
+
+
+
     # ================= Populate Hamiltonian =================
     for i in range(nstates):
 
@@ -2977,6 +2997,7 @@ def cp2k_compute_adi(q, params, full_id):
 
         for j in range(nstates):
             obj.time_overlap_adi.set(i, j, float(st_ci[i, j]))
+            obj.overlap_adi.set(i,  j, float(s_ci[i,j]) )
 
     # ================= Forces =================
     obj.d1ham_adi = CMATRIXList()
@@ -3018,6 +3039,7 @@ def cp2k_compute_adi(q, params, full_id):
     # params["MO_prev"][itraj] = MO_curr.copy()
     params["data_prev"][itraj] = copy.deepcopy(data_curr)
     params["coordinates_prev"][itraj] = coordinates.copy()
+    params["s_ci_inv_prev"][itraj] = copy.deepcopy(s_ci_inv_curr)
     params["is_first_time"][itraj] = False
 
     return obj