diff --git a/qstack/spahm/compute_spahm.py b/qstack/spahm/compute_spahm.py
index 541803b..25c7bdf 100644
--- a/qstack/spahm/compute_spahm.py
+++ b/qstack/spahm/compute_spahm.py
@@ -114,7 +114,7 @@ def get_guess_orbitals_grad(mol, guess, field=None):
     return e, eigenvalue_grad(mol, e, c, s1, h1), de_dfield
 
 
-def get_spahm_representation(mol, guess_in, xc="pbe", field=None):
+def get_spahm_representation(mol, guess_in, xc="pbe", field=None, verbose=1):
     """Compute the ε-SPAHM molecular representation.
 
     Reference:
@@ -128,6 +128,7 @@ def get_spahm_representation(mol, guess_in, xc="pbe", field=None):
         xc (str): Exchange-correlation functional name. Defaults to 'pbe'.
         field (numpy ndarray, optional): 3-component uniform electric field (∇φ) in atomic units.
             Defaults to None.
+        verbose (int): Level of verbosity
 
     Returns:
         numpy ndarray: SPAHM representation consisting of occupied orbital eigenvalues.
@@ -136,7 +137,7 @@ def get_spahm_representation(mol, guess_in, xc="pbe", field=None):
     """
     guess = get_guess(guess_in)
     e, _v = get_guess_orbitals(mol, guess, xc, field=field)
-    e1    = get_occ(e, mol.nelec, mol.spin)
+    e1    = get_occ(e, mol.nelec, mol.spin, verbose=verbose)
     return e1
 
 
diff --git a/qstack/spahm/guesses.py b/qstack/spahm/guesses.py
index 72ed240..a3bcaf6 100644
--- a/qstack/spahm/guesses.py
+++ b/qstack/spahm/guesses.py
@@ -197,13 +197,13 @@ def get_guess(arg):
     return guesses_dict[arg]
 
 
-def check_nelec(nelec, nao):
+def check_nelec(nelec, nao, full_shell_warning=1):
     """Validate that the number of electrons can be accommodated by available orbitals.
 
     Args:
         nelec (tuple or int): Number of electrons (alpha, beta) or total.
         nao (int): Number of atomic orbitals.
-
+        full_shell_warning (int): if < 1 supresses warnings about fully occupied shells
     Raises:
         RuntimeError: If there are more electrons than available orbitals.
 
@@ -212,12 +212,12 @@ def check_nelec(nelec, nao):
     """
     if np.any(np.array(nelec) > nao):
         raise RuntimeError(f'Too many electrons ({nelec}) for {nao} orbitals')
-    elif np.any(np.array(nelec) == nao):
+    elif np.any(np.array(nelec) == nao) and full_shell_warning>0:
         msg = f'{nelec} electrons for {nao} orbitals. Is the input intended to have a complete shell?'
         warnings.warn(msg, RuntimeWarning, stacklevel=2)
 
 
-def get_occ(e, nelec, spin):
+def get_occ(e, nelec, spin, verbose=1):
     """Extract occupied orbital eigenvalues/energies.
 
     Args:
@@ -225,13 +225,14 @@ def get_occ(e, nelec, spin):
             or possibly arrays of larger dimensionality.
         nelec (tuple): Number of (alpha, beta) electrons.
         spin (int or None): Spin multiplicity. If None, assumes closed-shell.
+        verbose (int): Level of verbosity (< 1 supresses full-shell warnings)
 
     Returns:
         numpy ndarray: Occupied eigenvalues. Shape depends on spin:
         - Closed-shell (spin=None): 1D array of occupied eigenvalues
         - Open-shell: 2D array (2, nocc) for alpha and beta separately
     """
-    check_nelec(nelec, e.shape[0])
+    check_nelec(nelec, e.shape[0], full_shell_warning=verbose)
     if spin is None:
         nocc = nelec[0]
         return e[:nocc,...]
@@ -243,20 +244,21 @@ def get_occ(e, nelec, spin):
         return e1
 
 
-def get_dm(v, nelec, spin):
+def get_dm(v, nelec, spin, verbose=1):
     """Construct density matrix from occupied molecular orbitals.
 
     Args:
         v (numpy ndarray): 2D array of MO coefficients (eigenvectors), columns are MOs.
         nelec (tuple): Number of (alpha, beta) electrons.
         spin (int or None): Spin multiplicity. If None, assumes closed-shell (RHF).
+        verbose (int): Level of verbosity (< 1 supresses full-shell warnings)
 
     Returns:
         numpy ndarray: Density matrix in AO basis.
         - Closed-shell: 2D array (nao, nao)
         - Open-shell: 3D array (2, nao, nao) for alpha and beta
     """
-    check_nelec(nelec, len(v))
+    check_nelec(nelec, len(v), full_shell_warning=verbose)
     if spin is None:
         nocc = nelec[0]
         dm = v[:,:nocc] @ v[:,:nocc].T
diff --git a/qstack/spahm/rho/utils.py b/qstack/spahm/rho/utils.py
index c87c2c5..93fec9b 100644
--- a/qstack/spahm/rho/utils.py
+++ b/qstack/spahm/rho/utils.py
@@ -106,7 +106,7 @@ def mols_guess(mols, xyzlist, guess, xc=defaults.xc, spin=None, readdm=None, pri
         xc (str): Exchange-correlation functional for guess. Defaults to defaults.xc.
         spin (list or None): List of spin multiplicities. Defaults to None.
         readdm (str, optional): Directory path to load pre-computed density matrices. Defaults to None.
-        printlevel (int): Verbosity level. Defaults to 0.
+        printlevel (int): Verbosity level. Defaults to 0. (if < 0 supresses warnings)
 
     Returns:
         list: List of density matrices (2D or 3D numpy arrays).
@@ -120,7 +120,7 @@ def mols_guess(mols, xyzlist, guess, xc=defaults.xc, spin=None, readdm=None, pri
             print(xyzfile, flush=True)
         if readdm is None:
             _e, v = spahm.get_guess_orbitals(mol, guess, xc=xc)
-            dm = guesses.get_dm(v, mol.nelec, mol.spin if sp is not None else None)
+            dm = guesses.get_dm(v, mol.nelec, mol.spin if sp is not None else None, verbose=printlevel+1)
         else:
             dm = np.load(f'{readdm}/{os.path.basename(xyzfile)}.npy')
             if spin and dm.ndim==2: