diff --git a/Algebraic-Diagrammatic-Construction/EE_ADC2.py b/Algebraic-Diagrammatic-Construction/EE_ADC2.py
index 7d00229..50d0469 100644
--- a/Algebraic-Diagrammatic-Construction/EE_ADC2.py
+++ b/Algebraic-Diagrammatic-Construction/EE_ADC2.py
@@ -139,7 +139,9 @@ def matvec(y):
 d_ia -= einsum('klac,klac->a', eri[o, o, v, v], t2)[None, :] * 0.25
 
 arg = np.argsort(np.absolute(diag))
-guess = np.eye(diag.size)[:, arg[:n_states]]
+guess = np.zeros((diag.size, n_states))
+for i in range(n_states):
+    guess[arg[i], i] = 1.0
 
 # Computes the EEs
 e_ee, v_ee = davidson(matvec, guess, diag, tol=tol)
diff --git a/Algebraic-Diagrammatic-Construction/IP_EA_ADC2.py b/Algebraic-Diagrammatic-Construction/IP_EA_ADC2.py
index b854d82..8964c58 100644
--- a/Algebraic-Diagrammatic-Construction/IP_EA_ADC2.py
+++ b/Algebraic-Diagrammatic-Construction/IP_EA_ADC2.py
@@ -142,7 +142,9 @@ def ea_matvec(y):
 # for the Davidson algorithm, and to generate the guess vectors
 diag = np.concatenate([np.diag(h_hh), e_ija.ravel()])
 arg = np.argsort(np.absolute(diag))
-guess = np.eye(diag.size)[:, arg[:n_states]]
+guess = np.zeros((diag.size, n_states))
+for i in range(n_states):
+    guess[arg[i], i] = 1.0
 
 # Compute the IPs
 e_ip, v_ip = davidson(ip_matvec, guess, diag, tol=tol)
@@ -159,7 +161,9 @@ def ea_matvec(y):
 # for the Davidson algorithm, and to generate the guess vectors
 diag = np.concatenate([np.diag(h_pp), -e_iab.ravel()])
 arg = np.argsort(np.absolute(diag))
-guess = np.eye(diag.size)[:, arg[:n_states]]
+guess = np.zeros((diag.size, n_states))
+for i in range(n_states):
+    guess[arg[i], i] = 1.0
 
 # Compute the EAs
 e_ea, v_ea = davidson(ea_matvec, guess, diag, tol=tol)
diff --git a/Algebraic-Diagrammatic-Construction/adc_helper.py b/Algebraic-Diagrammatic-Construction/adc_helper.py
index 1d80fc2..6631862 100644
--- a/Algebraic-Diagrammatic-Construction/adc_helper.py
+++ b/Algebraic-Diagrammatic-Construction/adc_helper.py
@@ -57,6 +57,7 @@ def davidson(matrix, guesses, diag, maxiter=None, sort_via_abs=True, tol=1e-10,
         de = np.linalg.norm(theta[:k] - theta_old)
         if de < tol:
             conv = True
+            b = np.dot(b, alpha)
             break
         else:
             if b.shape[1] >= (k * nvec_per_root):