From 109e4b4adbb067af8ddb208dff58771b93c1c0b9 Mon Sep 17 00:00:00 2001
From: Matthew R Hermes <mrhermes@uchicago.edu>
Date: Thu, 14 Nov 2024 16:37:55 -0600
Subject: [PATCH] fix grad

---
 my_pyscf/lassi/excitations.py | 15 ++++++++++-----
 1 file changed, 10 insertions(+), 5 deletions(-)

diff --git a/my_pyscf/lassi/excitations.py b/my_pyscf/lassi/excitations.py
index 46353ce1..796f16f0 100644
--- a/my_pyscf/lassi/excitations.py
+++ b/my_pyscf/lassi/excitations.py
@@ -534,7 +534,8 @@ def get_hpp_xp (self, h1, h2, ci0, si_p, norb_f, nelec_f, ecore=0, nroots=1, **k
 
     def _get_grad (self, ci0, si_p, hpq_xq, hpp_xp, nroots=None):
         # Compute the gradient of the target interacting energy
-        grad = []
+        grad_ext = []
+        grad_int = []
         for solver, c, hc1, hc2 in zip (self.fcisolvers, ci0, hpq_xq, hpp_xp):
             if nroots is None: nroots = solver.nroots
             hc = si_p[:,None,None] * (hc1 + hc2)
@@ -543,11 +544,15 @@ def _get_grad (self, ci0, si_p, hpq_xq, hpp_xp, nroots=None):
                 hc = solver.transformer.vec_det2csf (hc, normalize=False)
             chc = np.dot (c.conj (), hc.T)
             hc = hc - np.dot (chc.T, c)
-            grad.append (hc.flat)
+            grad_ext.append (hc.flat)
+            grad_int.append (chc)
+        grad_int[0] -= grad_int[1].T
+        grad_int[1] = -grad_int[0].T
+        grad = []
+        for i, e in zip (grad_int, grad_ext):
+            grad.append (e)
             if nroots>1:
-                # TODO: figure out how this gradient should actually work
-                chc -= chc.T
-                grad.append (np.zeros_like (chc[np.tril_indices (nroots, k=-1)]))
+                grad.append (i[np.tril_indices (nroots, k=-1)])
         return np.concatenate (grad)
 
     def _1shot (self, h0, h1, h2, ci0, hpq_xq, hpp_xp, nroots=1, ovlp_thresh=1e-3):