From 84ee0e57324690fa20461254965272f6dd4addae Mon Sep 17 00:00:00 2001
From: YAY-C <yannick.calvinoalonso@epfl.ch>
Date: Mon, 22 Apr 2024 20:19:01 +0200
Subject: [PATCH 01/11] Checkout test from yannick/fix_sad

---
 qstack/spahm/guesses.py                 | 104 ++----------------------
 tests/data/645688c7c139c6414b5c0b00.xyz |  20 +++++
 tests/test_SAD.py                       |  35 ++++++++
 3 files changed, 60 insertions(+), 99 deletions(-)
 create mode 100644 tests/data/645688c7c139c6414b5c0b00.xyz
 create mode 100644 tests/test_SAD.py

diff --git a/qstack/spahm/guesses.py b/qstack/spahm/guesses.py
index 8fb5671..bf183e3 100644
--- a/qstack/spahm/guesses.py
+++ b/qstack/spahm/guesses.py
@@ -1,5 +1,4 @@
 import sys
-import warnings
 import numpy
 import scipy
 import pyscf
@@ -7,27 +6,11 @@
 from qstack.spahm.LB2020guess import LB2020guess as LB20
 
 def hcore(mol, *_):
-  """Uses the diagonalization of the core to compute the guess Hamiltonian.
-
-  Args:
-    mol (pyscf Mole): pyscf Mole object.
-
-  Returns:
-    A numpy ndarray containing the computed approximate Hamiltonian.
-  """
   h  = mol.intor_symmetric('int1e_kin')
   h += mol.intor_symmetric('int1e_nuc')
   return h
 
 def GWH(mol, *_):
-  """Uses the generalized Wolfsberg-Helmholtz to compute the guess Hamiltonian.
-
-  Args:
-    mol (pyscf Mole): pyscf Mole object.
-
-  Returns:
-    A numpy ndarray containing the computed approximate Hamiltonian.
-  """
   h = hcore(mol)
   S = mol.intor_symmetric('int1e_ovlp')
   K = 1.75 # See J. Chem. Phys. 1952, 20, 837
@@ -41,32 +24,19 @@ def GWH(mol, *_):
   return h_gwh
 
 def SAD(mol, func):
-  """Uses the superposition of atomic densities to compute the guess Hamiltonian.
-
-  Args:
-    mol (pyscf Mole): pyscf Mole object.
-    func (str): Exchange-correlation functional.
-
-  Returns:
-    A numpy ndarray containing the computed approximate Hamiltonian.
-  """
   hc = hcore(mol)
+  print(f"SAD-Hcore = {hc.shape}")
   dm =  pyscf.scf.hf.init_guess_by_atom(mol)
+  print(f"SAD-DM = {dm.shape}")
   mf = pyscf.dft.RKS(mol)
   mf.xc = func
   vhf = mf.get_veff(dm=dm)
+  print(f"XC-dunc = {func} (converged = {mf.converged})")
+  print(f"SAD-vhf = {vhf.shape}")
   fock = hc + vhf
   return fock
 
 def SAP(mol, *_):
-  """Uses the superposition of atomic potentials to compute the guess Hamiltonian.
-
-  Args:
-    mol (pyscf Mole): pyscf Mole object.
-
-  Returns:
-    A numpy ndarray containing the computed approximate Hamiltonian.
-  """
   mf = pyscf.dft.RKS(mol)
   vsap = mf.get_vsap()
   t = mol.intor_symmetric('int1e_kin')
@@ -74,46 +44,17 @@ def SAP(mol, *_):
   return fock
 
 def LB(mol, *_):
-  """Uses the Laikov-Briling model with HF-based parameters to compute the guess Hamiltonian.
-
-  Args:
-    mol (pyscf Mole): pyscf Mole object.
-
-  Returns:
-    A numpy ndarray containing the computed approximate Hamiltonian.
-  """
   return LB20(parameters='HF').Heff(mol)
 
 def LB_HFS(mol, *_):
-  """ Laikov-Briling using HFS-based parameters
-
-  Args:
-    mol (pyscf Mole): pyscf Mole object.
-
-  Returns:
-    A numpy ndarray containing the computed approximate Hamiltonian.
-  """
   return LB20(parameters='HFS').Heff(mol)
 
 def solveF(mol, fock):
-  """Computes the eigenvalues and eigenvectors corresponding to the given Hamiltonian.
-
-  Args:
-    mol (pyscf Mole): pyscf Mole object.
-    fock (numpy ndarray): Approximate Hamiltonian.
-  """
   s1e = mol.intor_symmetric('int1e_ovlp')
+  print(f"1e-overlap = {s1e.shape}, Fock = {fock.shape}")
   return scipy.linalg.eigh(fock, s1e)
 
 def get_guess(arg):
-  """Returns the function of the method selected to compute the approximate hamiltoninan
-
-  Args:
-    arg (str): Approximate Hamiltonian
-
-  Returns:
-    The function of the selected method.
-  """
   arg = arg.lower()
   guesses = {'core':hcore, 'sad':SAD, 'sap':SAP, 'gwh':GWH, 'lb':LB, 'huckel':'huckel', 'lb-hfs':LB_HFS}
   if arg not in guesses.keys():
@@ -121,29 +62,7 @@ def get_guess(arg):
     exit(1)
   return guesses[arg]
 
-
-def check_nelec(nelec, nao):
-    """ Checks if there is enough orbitals
-    for the electrons"""
-    if numpy.any(numpy.array(nelec) > nao):
-        raise RuntimeError(f'Too many electrons ({nelec}) for {nao} orbitals')
-    elif numpy.any(numpy.array(nelec) == nao):
-        msg = f'{nelec} electrons for {nao} orbitals. Is the input intended to have a complete shell?'
-        warnings.warn(msg)
-
-
 def get_occ(e, nelec, spin):
-    """Returns the occupied subset of e
-
-    Args:
-      e (numpy ndarray): Energy eigenvalues.
-      nelec(tuple): Number of alpha and beta electrons.
-      spin(int): Spin.
-
-    Returns:
-      A numpy ndarray containing the occupied eigenvalues.
-    """
-    check_nelec(nelec, e.shape[0])
     if spin==None:
         nocc = nelec[0]
         return e[:nocc,...]
@@ -154,20 +73,7 @@ def get_occ(e, nelec, spin):
         e1[1,:nocc[1],...] = e[:nocc[1],...]
         return e1
 
-
 def get_dm(v, nelec, spin):
-  """Computes the density matrix.
-
-  Args:
-    v (numpy ndarray): Eigenvectors of a previously solve Hamiltoinan.
-    nelec(tuple): Number of alpha and beta electrons.
-    spin(int): Spin.
-
-  Return:
-    A numpy ndarray containing the density matrix computed using the guess Hamiltonian.
-  """
-
-  check_nelec(nelec, len(v))
   if spin==None:
     nocc = nelec[0]
     dm = v[:,:nocc] @ v[:,:nocc].T
diff --git a/tests/data/645688c7c139c6414b5c0b00.xyz b/tests/data/645688c7c139c6414b5c0b00.xyz
new file mode 100644
index 0000000..0f8d329
--- /dev/null
+++ b/tests/data/645688c7c139c6414b5c0b00.xyz
@@ -0,0 +1,20 @@
+18
+## STRUCTURE-ID: 645688c7c139c6414b5c0b00 ; CHARGE = -1; SPIN = 3
+I -1.4436121126949015 0.42029126592817107 -1.6632891671244727
+C -3.33580743801813 -0.15353604789979267 -0.7339542845183026
+H -3.1287423587721643 -0.9434022824024976 -0.01933696409311693
+H -3.9964207871651225 -0.4901283664418908 -1.5229488133526843
+H -3.7274422375299143 0.7217041183959343 -0.2282996045106841
+O 2.7223903943392993 0.8261760398398512 1.3050834155297588
+C 0.9888087204630703 -2.379094167649743 -0.5398992858009954
+H 1.6940986505047444 -3.152743685891341 -0.26448601467465105
+H 1.4948333510508294 -1.5155711156764566 -0.9595290502133175
+H 0.23509595271835407 -2.753118014964276 -1.2220792965839438
+H 3.6453537356643806 0.6064576399990429 1.069461740053683
+Rh 0.014459630959218359 0.9679146681456635 0.7456381302221328
+I -1.880589614365632 1.6750428531997783 2.5349367788411588
+I -0.030250731343837742 -1.749134256908 1.27855077077407
+O 2.4116878114136866 0.7236029720756493 -0.9162960476741423
+C 1.9293704614065443 0.9262914436358692 0.19859748723285423
+O 1.6398052199498723 3.5271918982221226 0.580225781514007
+C 0.766961351419702 2.742055038391913 0.35762442437864966
diff --git a/tests/test_SAD.py b/tests/test_SAD.py
new file mode 100644
index 0000000..7279dab
--- /dev/null
+++ b/tests/test_SAD.py
@@ -0,0 +1,35 @@
+#!/usr/bin/env python3
+
+import os
+import numpy as np
+from qstack import compound
+from qstack.spahm.rho import atom
+
+
+def test_water():
+    print("Running water-test")
+    path = os.path.dirname(os.path.realpath(__file__))
+    mol = compound.xyz_to_mol(path+'/data/H2O.xyz', 'sto3g', charge=0, spin=0)
+
+    Xsad = atom.get_repr(mol, ["H", "O"], 0, 0, dm=None,
+                      xc = 'hf', guess='sad', model='lowdin-long-x', auxbasis='ccpvdzjkfit')
+
+    print(np.array([*Xsad[:,1]]).shape)
+    return 0
+
+
+def test_other():
+    print("Running weird-test")
+    path = os.path.dirname(os.path.realpath(__file__))
+    mol = compound.xyz_to_mol(path+'/data/645688c7c139c6414b5c0b00.xyz', 'sto3g', charge=-1, spin=2)
+
+    Xsad = atom.get_repr(mol, ["H", "O", "C", "I", "Rh"], -1, 2, dm=None,
+                      xc = 'hf', guess='sad', model='lowdin-long-x', auxbasis='def2tzvpjkfit')
+
+    print(np.array([*Xsad[:,1]]).shape)
+    return 0
+
+
+if __name__ == '__main__':
+    test_water()
+    test_other()

From 0f4ba3988a96ea681bdd7134f3a27146c9a94916 Mon Sep 17 00:00:00 2001
From: calvinoY <yannick.calvinoalonso@epfl.ch>
Date: Tue, 23 Apr 2024 17:47:38 +0200
Subject: [PATCH 02/11] refined failed test (water example)

---
 tests/test_SAD.py | 31 ++++++++++++++++++++++---------
 1 file changed, 22 insertions(+), 9 deletions(-)

diff --git a/tests/test_SAD.py b/tests/test_SAD.py
index 7279dab..edd6622 100644
--- a/tests/test_SAD.py
+++ b/tests/test_SAD.py
@@ -6,10 +6,10 @@
 from qstack.spahm.rho import atom
 
 
-def test_water():
+def test_water_open_shell():
     print("Running water-test")
     path = os.path.dirname(os.path.realpath(__file__))
-    mol = compound.xyz_to_mol(path+'/data/H2O.xyz', 'sto3g', charge=0, spin=0)
+    mol = compound.xyz_to_mol(path+'/data/H2O.xyz', 'sto3g', charge=1, spin=1) ## test breaks when effective open-shell caluclation is needed
 
     Xsad = atom.get_repr(mol, ["H", "O"], 0, 0, dm=None,
                       xc = 'hf', guess='sad', model='lowdin-long-x', auxbasis='ccpvdzjkfit')
@@ -18,18 +18,31 @@ def test_water():
     return 0
 
 
-def test_other():
-    print("Running weird-test")
+def test_water_close_shell():
+    print("Running water-test")
     path = os.path.dirname(os.path.realpath(__file__))
-    mol = compound.xyz_to_mol(path+'/data/645688c7c139c6414b5c0b00.xyz', 'sto3g', charge=-1, spin=2)
+    mol = compound.xyz_to_mol(path+'/data/H2O.xyz', 'sto3g', charge=0, spin=0) ## test breaks when effective open-shell caluclation is needed
 
-    Xsad = atom.get_repr(mol, ["H", "O", "C", "I", "Rh"], -1, 2, dm=None,
-                      xc = 'hf', guess='sad', model='lowdin-long-x', auxbasis='def2tzvpjkfit')
+    Xsad = atom.get_repr(mol, ["H", "O"], 0, 0, dm=None,
+                      xc = 'hf', guess='sad', model='lowdin-long-x', auxbasis='ccpvdzjkfit')
 
     print(np.array([*Xsad[:,1]]).shape)
     return 0
 
+## This is the original failure (spotted) for a xyz-structutre from Reiher-group Monstanto DB
+#def test_other():
+#    print("Running weird-test")
+#    path = os.path.dirname(os.path.realpath(__file__))
+#    mol = compound.xyz_to_mol(path+'/data/645688c7c139c6414b5c0b00.xyz', 'sto3g', charge=-1, spin=2)
+#
+#    Xsad = atom.get_repr(mol, ["H", "O", "C", "I", "Rh"], -1, 2, dm=None,
+#                      xc = 'hf', guess='sad', model='lowdin-long-x', auxbasis='def2tzvpjkfit')
+#
+#    print(np.array([*Xsad[:,1]]).shape)
+#    return 0
+
 
 if __name__ == '__main__':
-    test_water()
-    test_other()
+    test_water_open_shell()
+    test_water_close_shell()
+    #test_other()

From 7486329b3958a9ff3324fb3abbd3a7b82d51ff55 Mon Sep 17 00:00:00 2001
From: calvinoY <yannick.calvinoalonso@epfl.ch>
Date: Wed, 24 Apr 2024 18:40:47 +0200
Subject: [PATCH 03/11] Added print-details

---
 qstack/spahm/compute_spahm.py | 2 ++
 qstack/spahm/guesses.py       | 2 ++
 tests/test_SAD.py             | 2 +-
 3 files changed, 5 insertions(+), 1 deletion(-)

diff --git a/qstack/spahm/compute_spahm.py b/qstack/spahm/compute_spahm.py
index 62e14f0..372d83a 100644
--- a/qstack/spahm/compute_spahm.py
+++ b/qstack/spahm/compute_spahm.py
@@ -39,6 +39,8 @@ def get_guess_dm(mol, guess, xc="pbe", openshell=None):
         A numpy ndarray containing the density matrix computed using the guess Hamiltonian.
     """
     e,v = get_guess_orbitals(mol, guess, xc)
+    print(f"Orbitals =  {v.shape}")
+    exit()
     return get_dm(v, mol.nelec, mol.spin if mol.spin>0 or not openshell is None else None)
 
 def get_spahm_representation(mol, guess_in, xc="pbe"):
diff --git a/qstack/spahm/guesses.py b/qstack/spahm/guesses.py
index bf183e3..100bb74 100644
--- a/qstack/spahm/guesses.py
+++ b/qstack/spahm/guesses.py
@@ -33,7 +33,9 @@ def SAD(mol, func):
   vhf = mf.get_veff(dm=dm)
   print(f"XC-dunc = {func} (converged = {mf.converged})")
   print(f"SAD-vhf = {vhf.shape}")
+  if vhf.ndim >2 : print(f"alpha and beta components are equal {(vhf[0] == vhf[1]).all()}")
   fock = hc + vhf
+  print(f"SAD-Fock = {fock.shape}")
   return fock
 
 def SAP(mol, *_):
diff --git a/tests/test_SAD.py b/tests/test_SAD.py
index edd6622..0facecb 100644
--- a/tests/test_SAD.py
+++ b/tests/test_SAD.py
@@ -11,7 +11,7 @@ def test_water_open_shell():
     path = os.path.dirname(os.path.realpath(__file__))
     mol = compound.xyz_to_mol(path+'/data/H2O.xyz', 'sto3g', charge=1, spin=1) ## test breaks when effective open-shell caluclation is needed
 
-    Xsad = atom.get_repr(mol, ["H", "O"], 0, 0, dm=None,
+    Xsad = atom.get_repr(mol, ["H", "O"], 1, 1, dm=None,
                       xc = 'hf', guess='sad', model='lowdin-long-x', auxbasis='ccpvdzjkfit')
 
     print(np.array([*Xsad[:,1]]).shape)

From ef83cdae1a302dfef82bdfd9b314d2ebc452c085 Mon Sep 17 00:00:00 2001
From: calvinoY <yannick.calvinoalonso@epfl.ch>
Date: Thu, 25 Apr 2024 13:03:07 +0200
Subject: [PATCH 04/11] Fixed SAS issue and cleaned test_SAD

---
 qstack/spahm/compute_spahm.py |  1 -
 qstack/spahm/guesses.py       |  2 +-
 tests/test_SAD.py             | 32 +++++++-------------------------
 3 files changed, 8 insertions(+), 27 deletions(-)

diff --git a/qstack/spahm/compute_spahm.py b/qstack/spahm/compute_spahm.py
index 372d83a..c3d1523 100644
--- a/qstack/spahm/compute_spahm.py
+++ b/qstack/spahm/compute_spahm.py
@@ -40,7 +40,6 @@ def get_guess_dm(mol, guess, xc="pbe", openshell=None):
     """
     e,v = get_guess_orbitals(mol, guess, xc)
     print(f"Orbitals =  {v.shape}")
-    exit()
     return get_dm(v, mol.nelec, mol.spin if mol.spin>0 or not openshell is None else None)
 
 def get_spahm_representation(mol, guess_in, xc="pbe"):
diff --git a/qstack/spahm/guesses.py b/qstack/spahm/guesses.py
index 100bb74..d76be56 100644
--- a/qstack/spahm/guesses.py
+++ b/qstack/spahm/guesses.py
@@ -34,7 +34,7 @@ def SAD(mol, func):
   print(f"XC-dunc = {func} (converged = {mf.converged})")
   print(f"SAD-vhf = {vhf.shape}")
   if vhf.ndim >2 : print(f"alpha and beta components are equal {(vhf[0] == vhf[1]).all()}")
-  fock = hc + vhf
+  fock = hc + (vhf if vhf.ndim == 2 else vhf[0])
   print(f"SAD-Fock = {fock.shape}")
   return fock
 
diff --git a/tests/test_SAD.py b/tests/test_SAD.py
index 0facecb..847c438 100644
--- a/tests/test_SAD.py
+++ b/tests/test_SAD.py
@@ -5,44 +5,26 @@
 from qstack import compound
 from qstack.spahm.rho import atom
 
-
 def test_water_open_shell():
-    print("Running water-test")
     path = os.path.dirname(os.path.realpath(__file__))
     mol = compound.xyz_to_mol(path+'/data/H2O.xyz', 'sto3g', charge=1, spin=1) ## test breaks when effective open-shell caluclation is needed
 
     Xsad = atom.get_repr(mol, ["H", "O"], 1, 1, dm=None,
                       xc = 'hf', guess='sad', model='lowdin-long-x', auxbasis='ccpvdzjkfit')
-
-    print(np.array([*Xsad[:,1]]).shape)
-    return 0
-
+    return Xsad
 
 def test_water_close_shell():
-    print("Running water-test")
     path = os.path.dirname(os.path.realpath(__file__))
     mol = compound.xyz_to_mol(path+'/data/H2O.xyz', 'sto3g', charge=0, spin=0) ## test breaks when effective open-shell caluclation is needed
 
     Xsad = atom.get_repr(mol, ["H", "O"], 0, 0, dm=None,
                       xc = 'hf', guess='sad', model='lowdin-long-x', auxbasis='ccpvdzjkfit')
+    return Xsad
 
-    print(np.array([*Xsad[:,1]]).shape)
-    return 0
-
-## This is the original failure (spotted) for a xyz-structutre from Reiher-group Monstanto DB
-#def test_other():
-#    print("Running weird-test")
-#    path = os.path.dirname(os.path.realpath(__file__))
-#    mol = compound.xyz_to_mol(path+'/data/645688c7c139c6414b5c0b00.xyz', 'sto3g', charge=-1, spin=2)
-#
-#    Xsad = atom.get_repr(mol, ["H", "O", "C", "I", "Rh"], -1, 2, dm=None,
-#                      xc = 'hf', guess='sad', model='lowdin-long-x', auxbasis='def2tzvpjkfit')
-#
-#    print(np.array([*Xsad[:,1]]).shape)
-#    return 0
-
+def test_equivalence():
+    Xos = test_water_open_shell()
+    Xcs = test_water_close_shell()
+    assert(not np.array_equal(Xos, Xcs))
 
 if __name__ == '__main__':
-    test_water_open_shell()
-    test_water_close_shell()
-    #test_other()
+    test_equivalence()

From cd3f1b12fef5336c5c3742c470c9059efe26124b Mon Sep 17 00:00:00 2001
From: calvinoY <yannick.calvinoalonso@epfl.ch>
Date: Thu, 25 Apr 2024 13:14:55 +0200
Subject: [PATCH 05/11] Recovered doc-strings from my-fork/yannick

---
 qstack/spahm/guesses.py | 105 +++++++++++++++++++++++++++++++++++++---
 1 file changed, 99 insertions(+), 6 deletions(-)

diff --git a/qstack/spahm/guesses.py b/qstack/spahm/guesses.py
index d76be56..0fd16df 100644
--- a/qstack/spahm/guesses.py
+++ b/qstack/spahm/guesses.py
@@ -1,4 +1,5 @@
 import sys
+import warnings
 import numpy
 import scipy
 import pyscf
@@ -6,11 +7,27 @@
 from qstack.spahm.LB2020guess import LB2020guess as LB20
 
 def hcore(mol, *_):
+  """Uses the diagonalization of the core to compute the guess Hamiltonian.
+
+  Args:
+    mol (pyscf Mole): pyscf Mole object.
+
+  Returns:
+    A numpy ndarray containing the computed approximate Hamiltonian.
+  """
   h  = mol.intor_symmetric('int1e_kin')
   h += mol.intor_symmetric('int1e_nuc')
   return h
 
 def GWH(mol, *_):
+  """Uses the generalized Wolfsberg-Helmholtz to compute the guess Hamiltonian.
+
+  Args:
+    mol (pyscf Mole): pyscf Mole object.
+
+  Returns:
+    A numpy ndarray containing the computed approximate Hamiltonian.
+  """
   h = hcore(mol)
   S = mol.intor_symmetric('int1e_ovlp')
   K = 1.75 # See J. Chem. Phys. 1952, 20, 837
@@ -24,21 +41,33 @@ def GWH(mol, *_):
   return h_gwh
 
 def SAD(mol, func):
+  """Uses the superposition of atomic densities to compute the guess Hamiltonian.
+
+  Args:
+    mol (pyscf Mole): pyscf Mole object.
+    func (str): Exchange-correlation functional.
+
+  Returns:
+    A numpy ndarray containing the computed approximate Hamiltonian.
+  """
   hc = hcore(mol)
   print(f"SAD-Hcore = {hc.shape}")
   dm =  pyscf.scf.hf.init_guess_by_atom(mol)
-  print(f"SAD-DM = {dm.shape}")
   mf = pyscf.dft.RKS(mol)
   mf.xc = func
   vhf = mf.get_veff(dm=dm)
-  print(f"XC-dunc = {func} (converged = {mf.converged})")
-  print(f"SAD-vhf = {vhf.shape}")
-  if vhf.ndim >2 : print(f"alpha and beta components are equal {(vhf[0] == vhf[1]).all()}")
   fock = hc + (vhf if vhf.ndim == 2 else vhf[0])
-  print(f"SAD-Fock = {fock.shape}")
   return fock
 
 def SAP(mol, *_):
+  """Uses the superposition of atomic potentials to compute the guess Hamiltonian.
+
+  Args:
+    mol (pyscf Mole): pyscf Mole object.
+
+  Returns:
+    A numpy ndarray containing the computed approximate Hamiltonian.
+  """
   mf = pyscf.dft.RKS(mol)
   vsap = mf.get_vsap()
   t = mol.intor_symmetric('int1e_kin')
@@ -46,17 +75,46 @@ def SAP(mol, *_):
   return fock
 
 def LB(mol, *_):
+  """Uses the Laikov-Briling model with HF-based parameters to compute the guess Hamiltonian.
+
+  Args:
+    mol (pyscf Mole): pyscf Mole object.
+
+  Returns:
+    A numpy ndarray containing the computed approximate Hamiltonian.
+  """
   return LB20(parameters='HF').Heff(mol)
 
 def LB_HFS(mol, *_):
+  """ Laikov-Briling using HFS-based parameters
+
+  Args:
+    mol (pyscf Mole): pyscf Mole object.
+
+  Returns:
+    A numpy ndarray containing the computed approximate Hamiltonian.
+  """
   return LB20(parameters='HFS').Heff(mol)
 
 def solveF(mol, fock):
+  """Computes the eigenvalues and eigenvectors corresponding to the given Hamiltonian.
+
+  Args:
+    mol (pyscf Mole): pyscf Mole object.
+    fock (numpy ndarray): Approximate Hamiltonian.
+  """
   s1e = mol.intor_symmetric('int1e_ovlp')
-  print(f"1e-overlap = {s1e.shape}, Fock = {fock.shape}")
   return scipy.linalg.eigh(fock, s1e)
 
 def get_guess(arg):
+  """Returns the function of the method selected to compute the approximate hamiltoninan
+
+  Args:
+    arg (str): Approximate Hamiltonian
+
+  Returns:
+    The function of the selected method.
+  """
   arg = arg.lower()
   guesses = {'core':hcore, 'sad':SAD, 'sap':SAP, 'gwh':GWH, 'lb':LB, 'huckel':'huckel', 'lb-hfs':LB_HFS}
   if arg not in guesses.keys():
@@ -64,7 +122,29 @@ def get_guess(arg):
     exit(1)
   return guesses[arg]
 
+
+def check_nelec(nelec, nao):
+    """ Checks if there is enough orbitals
+    for the electrons"""
+    if numpy.any(numpy.array(nelec) > nao):
+        raise RuntimeError(f'Too many electrons ({nelec}) for {nao} orbitals')
+    elif numpy.any(numpy.array(nelec) == nao):
+        msg = f'{nelec} electrons for {nao} orbitals. Is the input intended to have a complete shell?'
+        warnings.warn(msg)
+
+
 def get_occ(e, nelec, spin):
+    """Returns the occupied subset of e
+
+    Args:
+      e (numpy ndarray): Energy eigenvalues.
+      nelec(tuple): Number of alpha and beta electrons.
+      spin(int): Spin.
+
+    Returns:
+      A numpy ndarray containing the occupied eigenvalues.
+    """
+    check_nelec(nelec, e.shape[0])
     if spin==None:
         nocc = nelec[0]
         return e[:nocc,...]
@@ -75,7 +155,20 @@ def get_occ(e, nelec, spin):
         e1[1,:nocc[1],...] = e[:nocc[1],...]
         return e1
 
+
 def get_dm(v, nelec, spin):
+  """Computes the density matrix.
+
+  Args:
+    v (numpy ndarray): Eigenvectors of a previously solve Hamiltoinan.
+    nelec(tuple): Number of alpha and beta electrons.
+    spin(int): Spin.
+
+  Return:
+    A numpy ndarray containing the density matrix computed using the guess Hamiltonian.
+  """
+
+  check_nelec(nelec, len(v))
   if spin==None:
     nocc = nelec[0]
     dm = v[:,:nocc] @ v[:,:nocc].T

From 2f04c5169386f86908a433ffb4d1e575b1edae9c Mon Sep 17 00:00:00 2001
From: calvinoY <yannick.calvinoalonso@epfl.ch>
Date: Thu, 25 Apr 2024 13:44:00 +0200
Subject: [PATCH 06/11] Cleanup and added warning for SAD guess

---
 qstack/spahm/compute_spahm.py | 1 -
 qstack/spahm/guesses.py       | 9 +++++++--
 2 files changed, 7 insertions(+), 3 deletions(-)

diff --git a/qstack/spahm/compute_spahm.py b/qstack/spahm/compute_spahm.py
index c3d1523..62e14f0 100644
--- a/qstack/spahm/compute_spahm.py
+++ b/qstack/spahm/compute_spahm.py
@@ -39,7 +39,6 @@ def get_guess_dm(mol, guess, xc="pbe", openshell=None):
         A numpy ndarray containing the density matrix computed using the guess Hamiltonian.
     """
     e,v = get_guess_orbitals(mol, guess, xc)
-    print(f"Orbitals =  {v.shape}")
     return get_dm(v, mol.nelec, mol.spin if mol.spin>0 or not openshell is None else None)
 
 def get_spahm_representation(mol, guess_in, xc="pbe"):
diff --git a/qstack/spahm/guesses.py b/qstack/spahm/guesses.py
index 0fd16df..f37d3b8 100644
--- a/qstack/spahm/guesses.py
+++ b/qstack/spahm/guesses.py
@@ -51,12 +51,17 @@ def SAD(mol, func):
     A numpy ndarray containing the computed approximate Hamiltonian.
   """
   hc = hcore(mol)
-  print(f"SAD-Hcore = {hc.shape}")
   dm =  pyscf.scf.hf.init_guess_by_atom(mol)
   mf = pyscf.dft.RKS(mol)
   mf.xc = func
   vhf = mf.get_veff(dm=dm)
-  fock = hc + (vhf if vhf.ndim == 2 else vhf[0])
+  if vhf.ndim == 2:
+      fock = hc + vhf
+  elif vhf.ndim == 3 and numpy.array_equal(vhf[0], vhf[1]):
+      fock = hc + vhf[0]
+  else:
+      msg = f'The effective potential ({func}) return different alpha and beta matrix components from atomicHF DM'
+      warnings.warn(msg)
   return fock
 
 def SAP(mol, *_):

From aea95445d5c91f71a81048e34d8605d41e924970 Mon Sep 17 00:00:00 2001
From: Yannick Calvino Alonso <yannick.calvinoalonso@epfl.ch>
Date: Fri, 26 Apr 2024 14:00:38 +0200
Subject: [PATCH 07/11] Imporved warning

---
 qstack/spahm/guesses.py | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/qstack/spahm/guesses.py b/qstack/spahm/guesses.py
index f37d3b8..a5446ed 100644
--- a/qstack/spahm/guesses.py
+++ b/qstack/spahm/guesses.py
@@ -57,11 +57,11 @@ def SAD(mol, func):
   vhf = mf.get_veff(dm=dm)
   if vhf.ndim == 2:
       fock = hc + vhf
-  elif vhf.ndim == 3 and numpy.array_equal(vhf[0], vhf[1]):
-      fock = hc + vhf[0]
   else:
-      msg = f'The effective potential ({func}) return different alpha and beta matrix components from atomicHF DM'
-      warnings.warn(msg)
+      fock = hc + vhf[0]
+      if numpy.array_equal(vhf[0], vhf[1]):
+        msg = f'The effective potential ({func}) return different alpha and beta matrix components from atomicHF DM'
+        warnings.warn(msg)
   return fock
 
 def SAP(mol, *_):

From 494c9dddf252e45325dd04329a4c608490ad5efa Mon Sep 17 00:00:00 2001
From: Yannick Calvino Alonso <yannick.calvinoalonso@epfl.ch>
Date: Mon, 29 Apr 2024 10:27:36 +0200
Subject: [PATCH 08/11] Removed unecessary file

---
 tests/data/645688c7c139c6414b5c0b00.xyz | 20 --------------------
 1 file changed, 20 deletions(-)
 delete mode 100644 tests/data/645688c7c139c6414b5c0b00.xyz

diff --git a/tests/data/645688c7c139c6414b5c0b00.xyz b/tests/data/645688c7c139c6414b5c0b00.xyz
deleted file mode 100644
index 0f8d329..0000000
--- a/tests/data/645688c7c139c6414b5c0b00.xyz
+++ /dev/null
@@ -1,20 +0,0 @@
-18
-## STRUCTURE-ID: 645688c7c139c6414b5c0b00 ; CHARGE = -1; SPIN = 3
-I -1.4436121126949015 0.42029126592817107 -1.6632891671244727
-C -3.33580743801813 -0.15353604789979267 -0.7339542845183026
-H -3.1287423587721643 -0.9434022824024976 -0.01933696409311693
-H -3.9964207871651225 -0.4901283664418908 -1.5229488133526843
-H -3.7274422375299143 0.7217041183959343 -0.2282996045106841
-O 2.7223903943392993 0.8261760398398512 1.3050834155297588
-C 0.9888087204630703 -2.379094167649743 -0.5398992858009954
-H 1.6940986505047444 -3.152743685891341 -0.26448601467465105
-H 1.4948333510508294 -1.5155711156764566 -0.9595290502133175
-H 0.23509595271835407 -2.753118014964276 -1.2220792965839438
-H 3.6453537356643806 0.6064576399990429 1.069461740053683
-Rh 0.014459630959218359 0.9679146681456635 0.7456381302221328
-I -1.880589614365632 1.6750428531997783 2.5349367788411588
-I -0.030250731343837742 -1.749134256908 1.27855077077407
-O 2.4116878114136866 0.7236029720756493 -0.9162960476741423
-C 1.9293704614065443 0.9262914436358692 0.19859748723285423
-O 1.6398052199498723 3.5271918982221226 0.580225781514007
-C 0.766961351419702 2.742055038391913 0.35762442437864966

From 7ffc81322abfc79ec630c155397ec49e0ebbe55d Mon Sep 17 00:00:00 2001
From: Yannick Calvino Alonso <yannick.calvinoalonso@epfl.ch>
Date: Mon, 29 Apr 2024 10:32:34 +0200
Subject: [PATCH 09/11] Fixed warning-msg

---
 qstack/spahm/guesses.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/qstack/spahm/guesses.py b/qstack/spahm/guesses.py
index a5446ed..156f900 100644
--- a/qstack/spahm/guesses.py
+++ b/qstack/spahm/guesses.py
@@ -59,7 +59,7 @@ def SAD(mol, func):
       fock = hc + vhf
   else:
       fock = hc + vhf[0]
-      if numpy.array_equal(vhf[0], vhf[1]):
+      if not numpy.array_equal(vhf[0], vhf[1]):
         msg = f'The effective potential ({func}) return different alpha and beta matrix components from atomicHF DM'
         warnings.warn(msg)
   return fock

From caa08dbfe3f2dfdff8c62473c4fd1fce31f48df3 Mon Sep 17 00:00:00 2001
From: Yannick Calvino Alonso <yannick.calvinoalonso@epfl.ch>
Date: Wed, 1 May 2024 21:15:36 +0200
Subject: [PATCH 10/11] Move SAD-test to test_saphm_a.py (renamed) and add test
 data

---
 tests/data/SPAHM_a_H2O/X_H2O-RC_SAD.npy | Bin 0 -> 10405 bytes
 tests/data/SPAHM_a_H2O/X_H2O_SAD.npy    | Bin 0 -> 5434 bytes
 tests/test_SAD.py                       |  30 --------------
 tests/test_atomic_spahm.py              |  24 -----------
 tests/test_spahm_a.py                   |  51 ++++++++++++++++++++++++
 5 files changed, 51 insertions(+), 54 deletions(-)
 create mode 100644 tests/data/SPAHM_a_H2O/X_H2O-RC_SAD.npy
 create mode 100644 tests/data/SPAHM_a_H2O/X_H2O_SAD.npy
 delete mode 100644 tests/test_SAD.py
 delete mode 100755 tests/test_atomic_spahm.py
 create mode 100755 tests/test_spahm_a.py

diff --git a/tests/data/SPAHM_a_H2O/X_H2O-RC_SAD.npy b/tests/data/SPAHM_a_H2O/X_H2O-RC_SAD.npy
new file mode 100644
index 0000000000000000000000000000000000000000..311909ec2534763a7fe77df0636f9f0569ef38e3
GIT binary patch
literal 10405
zcmeI0do)#D)c+51SCP9YB=l55(y6rOK}m!{QX&aalA}XJrBDi872VvDR6N}zlyel-
z=^-H%5|RqJ-^%THYp<m-e*gZ)c;CPBhduV1bA9)B&bh}}`*Y4(ueVyek&hNg3(<1-
z*t5-FYpJSMsEyV<RV~lm{w)8kzFT(tyL<dzHrwj6$K%)f9<Qx_9>4Cj7tT}FpEFN&
zzv}=0pIDf{Hs63<e!)81cKdti><aK<d2jXi-x@5o#ouGwZr?pD|A1{QF(wUXiTS$Y
zD3i~U7Rlr{<#S^S9A!?jWeRFr@>vR4@*iOekyLk9u%Cw*Q`nK8Mx)u#nIgxcjxt4U
znPS?O0*-vY)<T)$*4EY*<G=o>f0n>vO1K?!W=__o;x^1FUNkT2K~uecJw?)%DfR0q
zR*8I(%xPX&*OSha{;exx%aonaHDb#Bzq&LQQ{L^ENw)-z_U_B~^4L~gDBs;?x71-7
z(92pRjA|IrZ(;c4aOPs*l67~vuA6Q0X21=2?q9Ra63kRQ4oy`z1qR7~N!Cr4=!*<<
z*Mhlh(!-nQR>4oAV@mSZko7gM4dV(=GoUeK={FmnJTz?RvifRdJ^Os2OWxMhFKpYz
z_<l1}FT3-l1pKHAd34@Z2GU7<klefMRNn4?R~htQ@W#5;PNKnrMe?(94znA{RpAv1
zO-lKJ@}GsBZ#w)GDjgOIMH|(@Bk`0};n%G&IvVCSOS%i}Z=PxyUE2ylr^EL9_}2k?
zmo!=)SqYxL{c1YGZ$R$g&9cub6_9>a{rFQ}5rh$)#DU{O^2Y{X(mrQ9>DvN~V7{rR
zjPDCGT*LeMd#wc+?bj~|>)soHfA(C=5HN0stYG>58-zN*_eO0XZ+aW7x|`QNcCa4e
zh#y`s)-el_TU|aT#mEY;hWHz@jBJK@bMG}NhW(BUb*WVX45RJ?wJL9hAWZ*+&Ay|<
z;FEoQmet1}(A}eGnbG+DPhCO8&n==;-Rh<~Q&x!~miq6CQjlfPvOy}OQJB#<sw+KJ
zVI1h%-YN8+ahUxvdr9r1QQ)^u3fjT#h9#sw4y3O8#E&@9dAPpi!c!Y5#{14H<^T2w
zFs4ci&K<WOf<<YSC2fMkz(}e5=ZV)3uuv`(3U~YtnxxM|q(5s(T_VKKYohbK^~b^U
z9|k~0Tt1^CupM~zmxS_{cYtHE#+!1hHaP#{puApwJ(!VsiI6^Pk^X!nb@3BFazsZ)
z^-K*XstM?0L;Li^K0}U{YSsnSPvDs5A>YAQ1+?9G?rCIR>7>sS{n00N@vp~z5Yd_9
zc2{fe$<LtZbd5$&_yoHD-R;eYtpa`W{;<g0cawQ(kUp;>{n<?Fx=s8j5S@hj4X4ju
zZ-Dew21T>?SHguf{2nHd_otuCU5U)gjPyB)^k<?j_Dk$1k?0s!1@g&Vcm~IOw|C@K
zR-u0rLf3k}??vWyEl&?fzd^oQT16`rCvoN~Jf~}%li@gT2-@z`$j_<M_!9AmRm;vK
z`8JC;jc<~W0^V4l)gzN0RCu;!Z^8m~PQ(?SNLG?QN7K~QC{jj^BQw=GaBlYxcJ^DW
z?>~(5b>H!czYMNGbzIkb7rgm743{j9^oMDjP-*jX;}(3JuujpaNTX`@Y2wF#=u4F1
ze2Y<jvz49;L5eZyc-2Q!c-jA5*P5;m$=_A3D<2buP<iaXmiS>fVf{+L`D~lO>nj}X
zp;EYZcS42<WIQgJvR_jPS{-n`Gl>5vS?s5W=)WfUb*oL>?knqo$Xl;HsUyOWtBU(m
zJq6b*@*VqM(}w+=DZu)YRXCq3{7u!XtP!m$^S-{PM-D7fJ35-X7eEZm=)SvQAy`~6
zZCK;20*=R*+8tihjh2)6j?lQd-;B(l{#JZ%PKO1!rNim^{RR-(&~&!snGl>H`x{97
z&W>Ze1nadhd$ar@S}7g&$t<;*<tPpkMBj_-pNJc8WOp`nsYBTjt=zcW0c1;_pFs34
zll|2sF0L3I_#&+wbq33x<*aBztr6aVw}rm_xhmQStB$mGpoHCh#lgb5(6{ox_H)ze
z5E3VHHM@HmDCM-;M~^RrF0#%JkTapFnZVC0LK{o@SHc6~2x+!7f=H4-MApw9=zP|+
z(Hw?f9JBhm%NUf2PBzKsk#!;U&?yFebU6HN!@bRSYS59?peJd;jqF`UAD{BxcC(Yo
zeQSl^%D&TW2sZlBIhjMqkHn40eJK6Ls{vI}=#ee%@m!XT(loqIv|W3~zD}+ir!##I
z#5W?o5I)yT?o-eT6LEjKxem;8)r8a}+u-;dlU)|;y5Ro4wKe}`wSw!>JN&|?bs!@*
zm6!Iu9EQcq?3DG&;O5JXx%rJRK}WYiW9@V{6d%ModGnrbO7tE8SF3e%mt+YrCY1;E
zK1&m3RR64v8F3R}&=jOp^c)7@_}9vsRlcpT^r4FUqLQz`eKA?z@7QP9C*TQ%`R_oE
z_;DdRbqAhw)yhm`Oslh>ES)dQu+@2?lv*yuP?H^+zU7htqd)al!+w_^V4JukdUkI=
zbl8`q{p{<3eMRv`LZ6%A1gVSn8vBtXI#$g`xBE?3VlZ6GG!?yM8J(A6B5h9!Gq?o;
z&#uUfLiVQZ%MIoZfpDzCiIZmrps@WPrifxITqgZ_N9vmJgNRP_9D3G{SSiNCm_ya>
z>jW4*b3F6b3H^Ws(Y3=WSNq|JjD6FxJ3YXvG|pcf*9=Oe&lCNTBXx}sKbb_QrPuA?
zTay9!cmC{+i#D_ZYyE}EM{~Zyle%dabM}4)+Ytc`we)wOOy=cB`rJ<X!zOjL5kD7+
z&f^_(T!$CdLuSc8!b_%ofb4NQ?%2ZjVAC)yF^l&fNKVFc*C+F0l0M%d{W(wSiXwh)
z5uJDI#qV4+t%j^O5!WV{Rsnyhq@tWY7kbG1Q$gn5Mdl?*`W#34lS%5jL;SoVI*(44
zY1zK10Ef>3??nrWp;-^Vhg|ah93gZ6O6KK4`n-?y$AHvDFTs91iOvk=H(H6RxiG~+
zAg<5qJ$hsvdiLV}K9v2;y<(5!Th#J+hEh0NkVAX7e7VAAInGoQTbYwHML7HND$X2i
z`^*j?`E#uFzKU8&0yjhm&9Ci2yusoJQNPZe<xOIr?e6Py0+PxjTjptS+9etvDQD0)
zb%j{JJ{adGuS^qEj(vkxSX~y2Zxi5LHd&_7u}_A>Q~q$rPeF*&Yxk<TZFv*hoA}8k
z`U`Dve)=SL!HU%eP`G+fx2A3d$mIU$kUeY&H+d4qv(8O}i*K=iU*hKk1M8nB`Lacl
z-^Ez#!S&0ym%;2+Kx<AfDw{qF*v7bC7V#f71N+e=`V;wIGN0ER{Yi&7_F!-F6&a{~
zi2EZq8P|K&7W+3Oe%26ubCR!7-sUk+um|;8aASl-m&5j$Ca0{l^^nOtl6>vuW=ONu
zHt5^24k{lcm_$!9fUP9{$z0`3r@9~f^t4M6-WveC><>OlJa2Hz-;%@lX7y+PMBF9X
z)wVV^6dHsx$_*@;a4ob@Kkl0yoFe)Y`x{C8K$>(_+_*dFa4Yi18Rk$#p6^NY)5-pM
zByRMfUGc2;JXm+GRUpy%EpmDp!$qvRKi96&k80(HpAl1m#Z+{oaLSyiG=Wpof8&(0
zEmLIzr)GG`Vm_kZ={w$c*cBE<^|M#fcY&_5(Tc|f+rg)=wO~Q+CYVpw9R%$P%4>sx
zuKSK%!xsQw?6b;a(ze4QlE0L!=gt#RN|a(kL5pH>)1qyEbec*2#QImwdPlyb^)TO2
zy5zy=6nH<q(WgUA3SLA^zWK9~4{FFgr(ve#Xx}stzaLOC{8#|EBtCJ!YoxdO^yPju
zXJ`JJ7}HM(RKFf5w|S0?$n{^(^DRy1ve1r8v6Vk=?Llc^(xfqd4$XCpxi#`635B^;
z&*igDL2E4!RTT3Fpu=J^lXPFFAWhNk0$$@Ylr6JY#S&GaCFa?j49ha)k-9;A`RZKc
zCxvxpj+GeRNXSJ;a;&a-KY4<LzLuOjwXq0k3)?aS4#cD0+oxLu+KNyin11UGY(cMM
zURxEV4x-{ni72)1ZuG`ytU*xw19IhJKSS5BPPR_S4p!<*bbV(IH?OM%nGf!cy*H4I
zDzr_@%nD1<W+(k*=^H(0n%S(IN*ic!E2JdXiuDr})2qtEE`COpYPhaxJnZMNFV@kB
zc5ZKBmLb*eib*FzQqh)BpP>coU!xPf%Cp6x56QpvQ^+gh13fW+KW6$cGERI}t?v8<
zHR>Iv`olT8j;d>YD)tlQq(td7_M`=W^9w|Gu3KJ-61<9<zS^}LySzZM$~Q(PC%r{n
z?fn-@=9i;W&Uah<yYtZHZfUB|-KDrcW(#m#qkFNR?&DY|%}lK(dS*HDbMpwvE$c(E
z&bl=jcKndL%*)Ds*9dY@O6E@*{eqqg;CXc}!+nk+{aMwG>xvr0e#G3d&OXDq(jDs+
z;I1Cy+v9LexRzq()K;wuG9K=0U-65Bk362c`w2WRot3!HneT9aa;M|E`kk<!VxrTK
zcmGX&l>z8{H=e>je*p|VRu=m_rT}y!{QjgL!E<jT^AZ)oeRfpE{Rv{=x=s;4KiaSk
zEyTx9(q0QJy*u~#-Jb>-x%fS7bjR;c#8N!>S~4&52;AqrGPpnfq%O$<>_>^{IFvU8
zB#exq?A~p6gy&?TAkj|2do2Y>GeAkZP(BSAPUU{!oTy<dzG58ExcrH|HS%M1^Kdb{
zd2dRlWb0KlgXEXLco>(LQ-&rf=$2{3WTVR$?D!WfZ)3l=?TS4i-^;#WvE+=;xwq`F
z``hjXZ_7p!EUYgfkMkLJ4m0DMFCndEdaq(nSFn{+qz)~ZTg^^Z8Si(!_n58gtSf32
za|!7YKT$;An&ej=6Vj{LAq2`!sG_oE3e>j3ySu)<h;6#U=G3kV_T^ISzm)j7nS}Kx
zk$lg5O9{_)vLKPQ-)oFHghZSL=SJ|}vn@!yF_zfB?o8}Qt_|x~W#arOotOf_o%JZ@
zVv$o~)MNIQVB8;+gX{I%jQu}ff&D}~Vg1oWoX?n>zBdUjA%CzqdDx{6=}5jd@NDfx
zr=E`;9H{I>jz{Vk?Gh!Z*f(v~&ioX%7KtzF(i<&UEDFMxel|`%ISC#aRP@f2X+j10
z>B$n|57}{KzX*xD>b0k6zMz5s(vwn~Zhu6dl2ZmpYVNT|BeDKtvVSj$Ka;dxoNN9F
zRhfqe-QF0>{*OFAp6GMQ{<u{*o+DOMx<GOly6r96&#Nv*%GoQj$B#VxbNzIW+Zz#o
z6a77R{+>I3&z-;L&fjzA@455)T>HN|cbGF-OjS3pKV+o#n~c<LnX@LyNMnMG>K<b<
zQooJK=;&5VM%z8xC^AwY>$?px866JBWTec+WOSV5hm&<yI3^>WyAnl4UPMQU<d2i}
zm*tp@&>R_xjG|gG8O@(Fl_I0#<o>&+8bwB<-}h2vw1C9N$US}OAVo$AhL0#RdOp3C
zBBLU59n8B-k&&NM3`IsQ4~;1@%Gn%Ekx_Q>dy0%ij^QtVqaggn-W1qHk<qT$K#GiF
zmGGCXV~-<6M!CnZ&JE{Mii|WrVxnpLiAltKgFQtM)^G8bzu{T@#V%EgzqAeHCsSl}
zbYPq!qv}%ZC!rMU#4o{Q#OHvCW<>;(BBS7Wm>_b3Fd3O_#b4}o>ZKGJaoXNeWW)={
zb#+RXQlwVrh;?W-n2cPHOs9yZWD_Qlog<haHhsij{+`MBi~XeY3yO@^?hdEOsQtYi
zRo60I?5AoF>(B*nP-NuiYDJNeg0=%iMk4yRC^C|lI7pF^dCvoijOZ@7&#QZJf7<kM
zUG5IpPnapzX)wiP#M+OE=6E0`k;z?{Ai{@JC^9<1#Pf1JiTm74`lCLC>&jS+{Rkh#
zI)kB@jE4LP8O>Qmk&$nmCq+iur|{g<Kj3+}xZ*x{))O+)!*va@u%9V=u+I8HOh%nw
zO(-&2AUKmEqp&IX{gG3^b61*$=fz0GeNGm@{VD6ib)6)B9uOUgJG&?{(mmKskx{bn
zM(RC`>%i~N!FD`%RTn%jnnDiM=j;8rKd(q#=EM*00M?P~Yo^F3+TNWaqcB-aM*UYY
zxd=I7GICQlrpTxup^+k^HA?3wGMYv5Uj$(?iV?@8bBmD6M;0a{Nkb!wjH2#3P-Ik|
ziuIGlaDG@2CL_8tCL``BCYMt>0t3IvC^>v4MMgB@CxPg1Ao)j%Fd4mT#$<Fl-h?8f
zOot^D8GTj5{>zD<H6vJ`P4fRWUPF<Q-04FU8HMf~{O30rWsrIW7i0gby4cSvqTg@|
z=c^XRQDl^``T<2o{K~jL%d>I4G~)l>2JC0l8|zD0<9u!B2#Smti!d4WN@FtmsKlel
zD4y>VMMjez#8YI{PU5Th^(iunsK;cKyQGRDqb(oKQDjt4_D{so8B9h@&p3*VmgvM$
zWQ2(RW3oSh#1AcZpvY)mhCf9{66E<yiT+HozcU@joB1D5Wc1~d2t`IA5hfHF?IYK-
z&x$EB`g`vDJ$HVS(cg3D@455)7rDRZ&fjzA_b+<?+vg5*HjAn0=D}iWx#|25gKV=;

literal 0
HcmV?d00001

diff --git a/tests/data/SPAHM_a_H2O/X_H2O_SAD.npy b/tests/data/SPAHM_a_H2O/X_H2O_SAD.npy
new file mode 100644
index 0000000000000000000000000000000000000000..39b132def9b2ad0954866609a2a5b15e77bc3871
GIT binary patch
literal 5434
zcmbW2dpwm{7sn5B2^F~<q;XtAql7m)NzZ;z5}`pTY2+3qJjWrTQYeaIqL^uN8A-S4
zCMh|mrkjQ&qoRl@m0LnmavX2%wd?bF|9#IN=lMK)t^Ip`Ypwm+--E8JTvo1@V8yZ$
zjDmur19?W)9HT@xqd6R-kVqb%w?5+MNM2Cz*JX$G;nBfg)}uq$M+JX*H~!BYj>YUb
zoNb){|0ka$6&Mk-IV#>XFp?K+x;ZACAGV&yTOTj~GcPzWG9sGKiwWe*Gc3rGj|hS&
zBjL<SWhCt-0vM@1%p`Y4+SpmbS;|>*HzPw*gZS}L!Sal(k0gu5a^o^`>1lhI$?lB2
zv9pwq#Fw>1M#0t9)oJ3(5C1tMKBE|r?#E0q#&I`hYA7ode`s3hm#-+fGs<7T@*tI&
z9tyG{Tt?-qtg1Vs_D$A?QU70A7N21Uq}#>Guvm{@^)_U7uS5-zJs#FxzKC1jrD)UQ
zkA|JBZtObkjcUkxkovomJ?s!v6LW{xvS<@>m>IltnqDBHN&Y>uZnw}v?sH8Pa#Sn0
zd@^?f`b2c>N&X75zT$yZ&h>16)RAEQ#x1;s`fRms`5&ob!J+HxuB<PX7P>Eh_qB|9
zdH=mR=woX_)gkx!sFcJ<$oslY+8cQ1ZIH!X@%Vu5S!$%(lC57+Ea)Iljr&w$K@l5W
zy`6OE#jY2q(d$1l`)ow0N};G&_Q7OYBuWa<SD8k8p3LeJt(;89WhZS5=ZO&YxIom9
zI)p+ZhILJ4-y`)MC+lC%97Lsm=<WSW{1PP*oueQ2j<VC$=t*1sJXFr=)A171vTQFJ
z(f-N9l0&Zgbl>rOS+fgj^pDZMmPy$v(r4n?+q`9_&=Dt|#)@ai)633Z=^Ni6MmfZf
zI39Ey64aMBjGNPElAF=t6P9#2t<X7i%$)YznBQ95q)*!nY=1iQA5A*RV!zwgJz8{l
z`Eh-hHa0ylsOenRq42$~IO69N(P{1u;OMAX)AGfCE~!J7G^-pbCv+Im9U?Q8X&O2-
z*Ep<*JE%j?YAav(v`U+nbS;S6SffNQB>nLsbzLNW6o}5{RbBaixtY^X`kS=>j@GBA
zsYshlcxuv?C5^Q`(pofK)cE6#P&VzPeO)HmM}amZeI6zKSxM@WBYqwbojcthcN}`A
zM$c4Wm%WWuq{W`U$y{ALh4v{l_@}`|o<4MU2iyFr7&(x6$&o%Ak^Z!ix+IAob)uun
z`K3jeCP{P0N4J{GOVJfZoHO|x3EHP5nEh5_7_lN@?j>YirKHc_`eQ-rl3WFTD55hp
z;JlH^0V!J3_ZW+tCqbM2eSWAca~N5W^TQ`|k0kRlAbnm&`m>hQb(Z+iAUb)^y|WJ<
zpFpL{=HHyPZ3yL;z&Xq#=VzGAU5m`if%Lh6^ygb$f_vcSDABQMij`2$e}&Q`HoU#k
zI86PSm$)+I$#mq{+V$6VmG@MHk<rbing&7>jXPX}-1$O3@3;-&9eTo6gI6h4{Lg~Z
zB;Repnu#@v=14pq>r!Q>jIJN-+LAZVPMGq$Snf=LyU@_y-X>MmR;W77FV<v$Ei8Wo
z`n=DOZ}wP1;a>b8)mr0!AwNt{NUu91RC#70OjL2qop90>CiPEFOSSnV$R>W~6Me-x
z$afM+IxG$GK+5!_y-jU_=-zOznJd>F752_NuAM#`C9=W)Q{spA1^q?@`2x4t<4b(_
zs8RN5WL{YSD!W!Yb(^6z>h^+qe<1$T)WFXm(SJblTbJ7fT-5#v$vp~<E*>*Nl^p0#
z^Hiu;t{42T=m9^!RD-@!6Xe&(zTt${v(U2ou;bB#R>-ON?c2@)H<W>X7&!0kft>Q~
z+gF4wMm~G3J$6|tp~WP=F)_#FjZFx8e(La0#oLW2pcG|4-{y%@+g}{4y={p0ll|=^
zer{vBIL|c-&05Re$*S`~TUD*y^nFZ_BGC^e`@h9)R}1>v-`S!1-A0u;m1@YHd_Rxq
zA0hjjNnAm5U+mqI0jfV<?VxapB<fBHi#;nN`~9iuCd(P?o`M(+KBF0cspiONeWTip
zuT;}^XJ!(r{Sc}KnBty$zxVyFBxIR3ELh6TL}uDHORiP#LE-PatLIg2MRUozm$XN9
z!_$1kHG3>*k;p}_JkK;_l<Yy4B;T5>SI&{sI;xzFs=G9Ay|7IFuG2~KzpcMEeC{Jr
z5QFCWsMKB(nW85%I>O)TnxnfZQ%-(r)I}}iUD&RpBzkXw6fVZpe!ixUYDoOs`{uEs
z=IkS?X!fS7D>CdQ5aPVv-r)9{vLVkuLk@Ml$h}W({4KNb<EevG39@@(FmZ`8@yR$f
zc0))d1vHyTxZa>vI`4dVOEQ<*C9gWk?7<DnaPmO4xZ@R7uDWHWGc`;tbSxK^Irme+
z#oh{wmp4;U%Alh&UTbwCubJ9i;c_hO#=n%z>)PC`)h{VyS$8I8`(<kAY<8Da&r9k$
zvVSuaD~BFrJaD;Ltd4G_DyHcUD4~DC$J?ci$0+|A@H2W0bjnQ=Hu8(RspFd}YOcKN
zr5r~hGcSC&PCYcXuXniKN3Hd>C{#J2jHWy2pVac^pi>F8l`eb@bc@^6ko2n*YSe|g
zri;PPt_aXE*yq>R#q?90Ud@9234c*PCx(yCTlJ3GKcqcN0jVJNqbQ9l^}5Jho)^WG
zY9ZUBx106+q)~_YF5DmCo}X}CtBS!-ny(ht=@=}De-m|xI(OXp_cZC7)Qi_1eYWe~
zP-@yI#-<c}plXb_<=4&~ptAhVcku>Vs3QX^xX%N1&>x3+P?u;6_!-y>IwcOeE&Fr^
zsHlM8gvxpql<8;IQs$wDDi?*i1Z|#yytE1>3q;cBjugzRe-ZRKgY;+F0MwN>0)FI!
zKxeB}PTj^;Ty)-?ese9^0UaxH@a<`)k!o<z%KMTgD4Y#*584m&GF=LNKK&T_Q#k|b
z8ukT0w}?*rm5cv8Z}LQ@y|z;&=enWMYufTJ$GM1W1Lvo7H_W|<%xkh7^x204{fVQY
zt}No`V-M)C62haDJeMKou>R<%ix#M?63$^q5S*VBYnb~}GB3vz=<^m;=ns$7rBn@m
zw1|#ZLwii#*i2MD6nIW{_I)aDa=-M2t~SarM$7m*`xa$2t>&3<|ENIoKE2)Gh=g!`
zYFl&X=T1TAmZH;2-8ZQpNPffJ%Q;sn`l(48X7vUcP1KQm56O9p<%Lh&-(~J+PZ#Dp
zE&L@s_k$qmV&H}Nz$QwO5BiF1$frHLbPjh`QAUf*?`LKY3bc!qcg{2UBq*FYG3<Y#
zRbb|4HrXYkiZUmD(ulq*$!|=TF@Lzx5NZ2T4;#Bo(bI19_<Y24B(PuNmbH0MaHJ0W
z*AYJ_3qXGo$q%{etQhi>B~m=IEp(jGL~?%8CMn{N0w+>$hBNp#(*Zx~J)qxo8uHUj
zGpeOGiK&cVZ~7ihYZd$+5B;Gkpx&sp;Q!7N@Uzbs^hHM@pEfDoQh=%`9&$RcYn_NP
zReCT#q<cEbx--7xL*rEBv%8h<Q|zT~MU?1ox_U!kMB<k3%th4;W+U0(K6OkvV1TOT
zKOE9gl|<E7OA8f~p9pfuemN5NH}5Mlyvsp6>jTPb&W=<678Z?+wKNJusi1$2?B7D-
zx0PHMR60JRnjDkk&aOTq_?vwHFww6e`*W5-yh6UEZk|#mbvA7Bu(-LC(k@?8KC%1B
z_vgPCYKBq{|4lJEe1;Pc`W;ZZUje1(&gc_B8GHj&>oowB-dO<Dp7j8z4Iw=kP#R=C
z(h2~zD;@x)T?2sHOY)P+IzJfzB@WWUfC?o#S|ophtlw(@KvA<*F`&}A0Z?;iPs4!P
zOWu18bupksy<0G#=8^a~dFNV>U_j+rRbfEgnbC~_b(1_tt{lODit^3Cfa<zzivd-!
zb{_^*`K>1yP;%*T0o%mE#rY!k9R}3q%vcPlOf9%zd<K0mpeoZr=Y(G!29#kNfTiaX
z0Kw7Q6JO-6kKh8fItUl1@>96*=Ch|@K<)W3fdSQA2Y&MEK<Ds60F;CmfMrPvg8>yk
z2QKo8H~^I0dbl`S_3AL7gguWipv1{gSHDsn24<@d=&;-XQ2x7TV6fD#0U&G|gNuAk
z8(hF4g>Z2m@VkovwK6gp1FG+d2d-<88Te@$0UfUN2@I$xe-{iW4P!41C^?H$7*K4*
z9T-rKgO@O%xa*+L%ZH#pJr+<`kQew#vIm`ZdjJ%F8-QhREC6B3JGjV`KNn#@?Pp+K
z{s*AXouof{qfl4b0`Mcd19V0b0Z^kn0;t)`FrXq@LolF>vtaI}&tP8b{GreN&k3N+
zp{`Lr_?a3FI;%ziQ2np%Frem1>tH}7O@;HLt^sq`(uaA`N1@M!QqZ6J_fXdX;^z|4
zQ9QR91Ild200vZ{>}q@tbKb)F+0h4c=d6Qyu{0`hpN|hif9{jI9El(CcF<9O--!XW
z&oc-EDoG6hHT(wvM8+2Y6`*H}0acyXfdREbD;EPwpXA?-13+ac0B}wbfVA-eP)b%d
z7*J{Fy)d8}ib20n9`cjo08m^%091_#0Fq@Y_2Da^3X^p(pjgCD9?|zE`MYldpdNPu
zpt29!VL+YsT8IJlS{M8`5I-x%Kwm)e|Fm6!0i~Y369X!7)5wos0aZrom0ke;IcDJJ
zKGAQ_f_%>P91N(u<(Dv^B(<SGi_4*27V&?<8~ljEKwr@n@{RpcFra8l0Mw8Q0IE$(
zi~)66;$IA?NtX^|K=qOMGD!;zsFddbsLF*+7*Iby%f*0tPWFF`Q@;S9n2;O{sD-9E
z7*G_^zee`QkoeBUUKmhw%6J%1isbv&L|=#O?=OXTr{pCJs8_$qVL&CM*kM3zCC`Jm
SZvmiY@fpK_U_N6MVESKhq0Alt

literal 0
HcmV?d00001

diff --git a/tests/test_SAD.py b/tests/test_SAD.py
deleted file mode 100644
index 847c438..0000000
--- a/tests/test_SAD.py
+++ /dev/null
@@ -1,30 +0,0 @@
-#!/usr/bin/env python3
-
-import os
-import numpy as np
-from qstack import compound
-from qstack.spahm.rho import atom
-
-def test_water_open_shell():
-    path = os.path.dirname(os.path.realpath(__file__))
-    mol = compound.xyz_to_mol(path+'/data/H2O.xyz', 'sto3g', charge=1, spin=1) ## test breaks when effective open-shell caluclation is needed
-
-    Xsad = atom.get_repr(mol, ["H", "O"], 1, 1, dm=None,
-                      xc = 'hf', guess='sad', model='lowdin-long-x', auxbasis='ccpvdzjkfit')
-    return Xsad
-
-def test_water_close_shell():
-    path = os.path.dirname(os.path.realpath(__file__))
-    mol = compound.xyz_to_mol(path+'/data/H2O.xyz', 'sto3g', charge=0, spin=0) ## test breaks when effective open-shell caluclation is needed
-
-    Xsad = atom.get_repr(mol, ["H", "O"], 0, 0, dm=None,
-                      xc = 'hf', guess='sad', model='lowdin-long-x', auxbasis='ccpvdzjkfit')
-    return Xsad
-
-def test_equivalence():
-    Xos = test_water_open_shell()
-    Xcs = test_water_close_shell()
-    assert(not np.array_equal(Xos, Xcs))
-
-if __name__ == '__main__':
-    test_equivalence()
diff --git a/tests/test_atomic_spahm.py b/tests/test_atomic_spahm.py
deleted file mode 100755
index 79f8351..0000000
--- a/tests/test_atomic_spahm.py
+++ /dev/null
@@ -1,24 +0,0 @@
-#!/usr/bin/env python3
-
-import os
-import numpy as np
-from qstack import compound
-from qstack.spahm.rho import atom
-
-
-def test_water():
-    path = os.path.dirname(os.path.realpath(__file__))
-    mol = compound.xyz_to_mol(path+'/data/H2O.xyz', 'minao', charge=0, spin=None)
-
-    X = atom.get_repr(mol, ["H", "O"], 0, None, dm=None,
-                      guess='LB', model='lowdin-long-x', auxbasis='ccpvdzjkfit')
-
-    X_true = np.load(path+'/data/SPAHM_a_H2O/X_H2O.npy', allow_pickle=True)
-    assert(X.shape == X_true.shape)
-    for a, a_true in zip(X, X_true):
-        assert(a[0] == a_true[0])                        # atom type
-        assert(np.linalg.norm(a[1]-a_true[1]) < 1e-08)   # atom representations
-
-
-if __name__ == '__main__':
-    test_water()
diff --git a/tests/test_spahm_a.py b/tests/test_spahm_a.py
new file mode 100755
index 0000000..d208adc
--- /dev/null
+++ b/tests/test_spahm_a.py
@@ -0,0 +1,51 @@
+#!/usr/bin/env python3
+
+import os
+import numpy as np
+from qstack import compound
+from qstack.spahm.rho import atom
+
+
+def test_water():
+    path = os.path.dirname(os.path.realpath(__file__))
+    mol = compound.xyz_to_mol(path+'/data/H2O.xyz', 'minao', charge=0, spin=None)
+
+    X = atom.get_repr(mol, ["H", "O"], 0, None, dm=None,
+                      guess='LB', model='lowdin-long-x', auxbasis='ccpvdzjkfit')
+
+    X_true = np.load(path+'/data/SPAHM_a_H2O/X_H2O.npy', allow_pickle=True)
+    assert(X.shape == X_true.shape)
+    for a, a_true in zip(X, X_true):
+        assert(a[0] == a_true[0])                        # atom type
+        assert(np.linalg.norm(a[1]-a_true[1]) < 1e-08)   # atom representations
+
+def test_water_SAD_guess_open_shell():
+    path = os.path.dirname(os.path.realpath(__file__))
+    mol = compound.xyz_to_mol(path+'/data/H2O.xyz', 'sto3g', charge=1, spin=1) ## test breaks when effective open-shell caluclation is needed
+
+    Xsad = atom.get_repr(mol, ["H", "O"], 1, 1, dm=None,
+                      xc = 'hf', guess='sad', model='lowdin-long-x', auxbasis='ccpvdzjkfit')
+    Xtrue = np.load(path+'/data/SPAHM_a_H2O/X_H2O-RC_SAD.npy', allow_pickle=True)
+    assert(Xsad.shape == Xtrue.shape)
+    for a, a_true in zip(Xsad, Xtrue):
+        assert(a[0] == a_true[0])                        # atom type
+        assert(np.linalg.norm(a[1]-a_true[1]) < 1e-08)   # atom representations
+
+def test_water_SAD_guess_close_shell():
+    path = os.path.dirname(os.path.realpath(__file__))
+    mol = compound.xyz_to_mol(path+'/data/H2O.xyz', 'sto3g', charge=0, spin=0) ## test breaks when effective open-shell caluclation is needed
+
+    Xsad = atom.get_repr(mol, ["H", "O"], 0, None, dm=None,
+                      xc = 'hf', guess='sad', model='lowdin-long-x', auxbasis='ccpvdzjkfit')
+    Xtrue = np.load(path+'/data/SPAHM_a_H2O/X_H2O_SAD.npy', allow_pickle=True)
+    assert(Xsad.shape == Xtrue.shape)
+    for a, a_true in zip(Xsad, Xtrue):
+        assert(a[0] == a_true[0])                        # atom type
+        assert(np.linalg.norm(a[1]-a_true[1]) < 1e-08)   # atom representations
+
+
+
+if __name__ == '__main__':
+    test_water()
+    test_water_SAD_guess_close_shell()
+    test_water_SAD_guess_open_shell()

From 4fcbecc47a78be20b9ff36bd16b755480c486561 Mon Sep 17 00:00:00 2001
From: Yannick Calvino Alonso <yannick.calvinoalonso@epfl.ch>
Date: Wed, 1 May 2024 21:32:46 +0200
Subject: [PATCH 11/11] Fix test_global (broken by caa08db)

---
 tests/test_global.py | 6 ++++--
 1 file changed, 4 insertions(+), 2 deletions(-)

diff --git a/tests/test_global.py b/tests/test_global.py
index ed84f09..ce4bf86 100644
--- a/tests/test_global.py
+++ b/tests/test_global.py
@@ -8,7 +8,8 @@ def test_avg_kernel():
 
     path = os.path.dirname(os.path.realpath(__file__))
     X_dir = os.path.join(path, 'data/SPAHM_a_H2O/')
-    mols = [np.load(os.path.join(X_dir, f), allow_pickle=True) for f in os.listdir(X_dir) if os.path.isfile(os.path.join(X_dir,f))]
+    mollist = [os.path.join(X_dir, s) for s in ['X_H2O.npy', 'X_rotated_H2O.npy', 'X_H2O_dist.npy']]
+    mols = [np.load(f, allow_pickle=True) for f in mollist]
     K = kernel.kernel(mols, akernel='L', gkernel='avg', sigma=1.0)
 
     true_K = np.array(  [[1.       ,  1.       ,  0.79179528], \
@@ -22,7 +23,8 @@ def test_avg_kernel():
 def test_rem_kernel():
     path = os.path.dirname(os.path.realpath(__file__))
     X_dir = os.path.join(path, 'data/SPAHM_a_H2O/')
-    mols = [np.load(os.path.join(X_dir, f), allow_pickle=True) for f in os.listdir(X_dir) if os.path.isfile(os.path.join(X_dir,f))]
+    mollist = [os.path.join(X_dir, s) for s in ['X_H2O.npy', 'X_rotated_H2O.npy', 'X_H2O_dist.npy']]
+    mols = [np.load(f, allow_pickle=True) for f in mollist]
     K = kernel.kernel(mols, akernel='L', gkernel='rem', sigma=1.0, gdict={'alpha':1.0, 'normalize':1})
 
     true_K = np.array(  [[1.      , 0.6528238, 1.      ], \