Merge pull request #121 from obackhouse/adc2_bug_fix

Fixes bug in initialisation of ADC guess vectors

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)

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)

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