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()