Fixed the DFTB+ workflow, added new Pythonic eigensolvers, small fixes to Bohmian module and step3 module

Author: alexvakimov

src/libra_py/__init__.py

@@ -23,6 +23,7 @@
            "data_stat",
            "data_visualize",
            "dynamics_plotting",
+           "eigensolvers",
            "fgr_py",
            "fit",
            "fix_motion",

src/libra_py/dynamics/bohmian/compute.py

@@ -68,7 +68,7 @@ def rho_lorentzian(q, Q, sigma):
     return y
 
 
-def quantum_potential_orginal(Q, sigma, mass, TBF):
+def quantum_potential_original(Q, sigma, mass, TBF):
     """
     Args:
     * Q (Tensor(ntraj, ndof)) - coordinates of all trajectories
@@ -92,6 +92,31 @@ def quantum_potential_orginal(Q, sigma, mass, TBF):
     return U
 
 
+def quantum_potential_original_gen(q, Q, sigma, mass, TBF):
+    """
+    Args:
+    * Q (Tensor(ntraj, ndof)) - coordinates of all trajectories
+    * sigma (Tensor(ndof)) - width parameters for each trajectory
+    * mass ( Tensor(1, ndof)) - masses of all DOFs, same for all trajectories
+    * TBF (object) - basis function reference (`rho_gaussian` or `rho_lorentzian`)
+
+    Returns:
+    Tensor(1) - quantum potential summed over all trajectory points
+    """
+
+    ntraj, ndof = Q.shape[0], Q.shape[1]
+    U = torch.zeros( (1,), requires_grad=True)
+    f = TBF(q, Q, sigma);
+    [deriv1] = torch.autograd.grad(f, [q], create_graph=True, retain_graph=True);
+    for i in range(ndof):
+        [deriv2] = torch.autograd.grad(deriv1[i], [q], create_graph=True, retain_graph=True);
+        u = -(0.25/mass[0,i])*( deriv2[i]/f  - 0.5 * (deriv1[i]/f)**2 );
+        U = U + u
+    return U
+
+
+
+
 def quantum_potential(Q, sigma, mass, TBF):
     """
     Compute quantum potential in a fully vectorized way.

src/libra_py/eigensolvers.py

@@ -0,0 +1,169 @@
+
+import numpy as np
+from scipy.linalg import cholesky, solve_triangular
+from scipy.linalg import eigh, eig
+
+def generalized_eigensolve_scipy(A, B, hermitian=True, sort=True):
+    """
+    Solve the generalized eigenvalue problem A v = λ B v using SciPy.
+
+    Args:
+        A (np.ndarray): Matrix A (N x N)
+        B (np.ndarray): Matrix B (N x N)
+        hermitian (bool): Set True if A and B are Hermitian/symmetric.
+        sort (bool): If True, sort eigenpairs by ascending eigenvalue.
+
+    Returns:
+        eigvals (np.ndarray): Eigenvalues (N,)
+        eigvecs (np.ndarray): Corresponding eigenvectors (N, N)
+
+    Example:
+
+    import numpy as np
+
+    # Symmetric A and B
+    A = np.array([[6., 2.], [2., 3.]])
+    B = np.array([[1., 0.], [0., 2.]])
+    eigvals, eigvecs = generalized_eigensolve_scipy(A, B)
+    print("Eigenvalues:", eigvals)
+    print("Eigenvectors (columns):\n", eigvecs)
+
+    """
+
+    if hermitian:
+        eigvals, eigvecs = eigh(A, B)
+    else:
+        eigvals, eigvecs = eig(A, B)
+
+    if sort:
+        idx = np.argsort(eigvals.real if not hermitian else eigvals)
+        eigvals = eigvals[idx]
+        eigvecs = eigvecs[:, idx]
+
+    # Normalize eigenvectors with respect to B: v† B v = 1
+    for i in range(eigvecs.shape[1]):
+        vec = eigvecs[:, i]
+        #norm = np.sqrt(np.conj(vec).T @ B @ vec)
+        norm = np.sqrt(np.vdot(vec, B @ vec))  # vdot does conjugation
+        eigvecs[:, i] = vec / norm
+
+
+    return eigvals, eigvecs
+
+
+
+def generalized_eigensolve_scipy_cholesky(A, B, hermitian=True, sort=True):
+    """
+    Robust solver for A v = λ B v, with B-orthonormal eigenvectors:
+        v_i^† B v_j = δ_ij (up to numerical error)
+    
+    For Hermitian A, B > 0.
+    """
+    if not hermitian:
+        raise NotImplementedError("Only Hermitian positive-definite B is supported in stable form")
+
+    # Cholesky decomposition: B = L L^H
+    L = cholesky(B, lower=True)
+
+    # Transform to standard eigenvalue problem: C = L^{-1} A L^{-H}
+    Linv_A = solve_triangular(L, A, lower=True, trans='N')
+    C = solve_triangular(L, Linv_A.T, lower=True, trans='C').T
+
+    # Solve C y = λ y
+    eigvals, y = eigh(C)
+
+    if sort:
+        idx = np.argsort(eigvals)
+        eigvals = eigvals[idx]
+        y = y[:, idx]
+
+    # Recover v = L^{-H} y
+    eigvecs = solve_triangular(L, y, lower=True, trans='C')
+
+    # At this point, eigvecs should satisfy v.T.conj() @ B @ v ≈ I
+    return eigvals, eigvecs
+
+
+"""
+
+Tests and conventions for this module:
+
+from liblibra_core import *
+from libra_py import eigensolvers, data_conv
+import numpy as np
+
+# ===== Do the solution ==========
+A = np.array([[0, 1.0], 
+              [1.0, 1.0]])
+B = np.array([[1.0, 0.0], 
+              [0.0, 1.0]])
+
+eigvals, eigvecs = eigensolvers.generalized_eigensolve_scipy(A, B, hermitian=False)
+
+
+# ======= Numpy multiplications ==========
+E = np.diag(eigvals)
+C = eigvecs
+print( E )
+print( C )
+print("LHS = ", A @ C)
+print("RHS = ", B @ C @ E)
+
+Selection deleted
+E = np.diag(eigvals)
+C = eigvecs
+print( E )
+print( C )
+print("LHS = ", A @ C)
+
+print("RHS = ", B @ C @ E)
+
+>>[[-0.61803399+0.j  0.        +0.j]
+>> [ 0.        +0.j  1.61803399+0.j]]
+>>[[ 0.85065081  0.52573111]
+>> [-0.52573111  0.85065081]]
+>>LHS =  [[-0.52573111  0.85065081]
+>> [ 0.3249197   1.37638192]]
+>>RHS =  [[-0.52573111+0.j  0.85065081+0.j]
+>> [ 0.3249197 +0.j  1.37638192+0.j]]
+
+
+# ======= Conversion to MATRIX ===========
+
+H = data_conv.nparray2MATRIX(A)
+U = data_conv.nparray2MATRIX(C)
+S = data_conv.nparray2MATRIX(B)
+e = data_conv.nparray2MATRIX(E)
+
+# ======== Libra multiplications ========
+
+(H * U).show_matrix()
+
+>>-0.52573111  0.85065081  
+>>0.32491970  1.3763819   
+
+(S * U * e).show_matrix()
+
+>>-0.52573111  0.85065081
+>>0.32491970  1.3763819
+
+# ========= Eigenvectors ==============
+print("First eigenvector")
+print(U.get(0, 0))
+print(U.get(1, 0))
+
+>>First eigenvector
+>>0.8506508083520399
+>>-0.5257311121191337
+
+print("First eigenvector")
+print(C[0, 0])
+print(C[1, 0])
+
+>>First eigenvector
+>>0.8506508083520399
+>>-0.5257311121191337
+
+
+"""
+

src/libra_py/packages/dftbplus/methods.py

@@ -144,6 +144,7 @@ def get_dftb_matrices(filename, act_sp1=None, act_sp2=None, correction=0, tol=0.
     """
 
     # Copy the content of the file to a temporary location
+    print(F"Reading file {filename}")
     f = open(filename, "r")
     A = f.readlines()
     f.close()
@@ -167,6 +168,8 @@ def get_dftb_matrices(filename, act_sp1=None, act_sp2=None, correction=0, tol=0.
     X = []
     norbs2 = int(norbs / 2)
 
+    print(F"norbs = {norbs}, norbs2 = {norbs2}, act_sp1 = {act_sp1}, act_sp2 = {act_sp2}")
+
     # For each k-point
     for ikpt in range(0, nkpts):
 
@@ -176,6 +179,8 @@ def get_dftb_matrices(filename, act_sp1=None, act_sp2=None, correction=0, tol=0.
         B = A[5 + ikpt * norbs: 5 + (1 + ikpt) * norbs]
 
         f1 = open(filename + "_tmp", "w")
+        x = MATRIX(norbs, norbs)
+
         for i in range(0, norbs):
 
             tmp = B[i].split()
@@ -186,18 +191,23 @@ def get_dftb_matrices(filename, act_sp1=None, act_sp2=None, correction=0, tol=0.
                 if tmp[j] == "NaN" or tmp[j] == "NAN":
                     z = 1000.0
                     # z = float(tmp[int(j%norbs2)]) #-1000.0
+                    sys.exit(0)
                 else:
                     try:
                         z = float(tmp[j])
+                        x.set(i,j, z)
                     except ValueError:
-                        z = 1000.0
+                        sys.exit(0)
+                        #z = 1000.0
                         # z = float(tmp[int(j%norbs2)]) #-1000.0
-                        pass
+                        #pass
                     except TypeError:
-                        z = 1000.0
+                        sys.exit(0)
+                        #z = 1000.0
                         # z = float(tmp[int(j%norbs2)]) #-1000.0
-                        pass
+                        #pass
 
+                """
                 if correction == 1:  # if we want to enforce correction of the matrix elements in blocks 01 and 10
                     #           norbs2    norbs
                     #       ______________
@@ -225,16 +235,18 @@ def get_dftb_matrices(filename, act_sp1=None, act_sp2=None, correction=0, tol=0.
                                 z = 0.0
                         else:
                             pass
-
-                line = line + "%10.8f  " % (z)
+                """
+                #line = line + "%10.8f  " % (z)
+                line = F"{line} {z}" 
             line = line + "\n"
 
             f1.write(line)
         f1.close()
 
         # Read in the temporary file - get the entire matrix
-        x = MATRIX(norbs, norbs)
-        x.Load_Matrix_From_File(filename + "_tmp")
+        #x = MATRIX(norbs, norbs)
+        #x.Load_Matrix_From_File(filename + "_tmp")
+      
 
         # Extract the sub-matrix of interest
         x_sub = MATRIX(nstates1, nstates2)

src/libra_py/workflows/nbra/mapping.py

@@ -215,17 +215,20 @@ def energy_arb(SD, e):
             This is the approximation of the SD energy
 
     """
+    if isinstance(e, np.ndarray)==True:
+        nbasis = e.shape[0]
+        sd = sd2indx(SD, nbasis)
+        res = 0.0
+        for i in sd:
+            res = res + e[i, i]
 
-    # nbasis = e.num_of_rows
-    nbasis = e.shape[0]
 
-    sd = sd2indx(SD, nbasis)
-    # res = 0.0+0.0j
-    res = 0.0
-
-    for i in sd:
-        # res = res + e.get(i,i)
-        res = res + e[i, i]
+    elif isinstance(e, CMATRIX)==True or isinstance(e, MATRIX)==True:
+        nbasis = e.num_of_rows
+        sd = sd2indx(SD, nbasis)
+        res = 0.0+0.0j  
+        for i in sd:
+            res = res + e.get(i,i)
 
     return res
 
@@ -246,14 +249,21 @@ def energy_mat_arb(SD, e, dE):
         CMATRIX(N,N): the matrix of energies in the SD basis. Here, N = len(SD) - the number of SDs.
 
     """
-    n = len(SD)
-    # E = CMATRIX(n,n)
-    E = np.zeros((n, n))
-
-    E0 = energy_arb(SD[0], e) + dE[0]  # *(1.0+0.0j)
-    for i in range(0, n):
-        # E.set(i,i, energy_arb(SD[i], e) + dE[i]*(1.0+0.0j) - E0 )
-        E[i, i] = energy_arb(SD[i], e) + dE[i] - E0
+    E = None
+
+    if isinstance(e, np.ndarray)==True:
+        n = len(SD)
+        E = np.zeros((n, n))
+        E0 = energy_arb(SD[0], e) + dE[0]
+        for i in range(0, n):
+            E[i, i] = energy_arb(SD[i], e) + dE[i] - E0
+
+    elif isinstance(e, CMATRIX)==True or isinstance(e, MATRIX)==True:
+        n = len(SD)
+        E = CMATRIX(n,n)
+        E0 = energy_arb(SD[0], e) + dE[0]*(1.0+0.0j)
+        for i in range(0, n):
+            E.set(i,i, energy_arb(SD[i], e) + dE[i]*(1.0+0.0j) - E0 )
 
     return E