From 0b6cb1f7ea18d0172985efc944043c5396bb1ca7 Mon Sep 17 00:00:00 2001
From: Matthew R Hermes <mrhermes@uchicago.edu>
Date: Fri, 27 Oct 2023 15:36:11 -0500
Subject: [PATCH] lassis "fullauto" bugfix counting spin states

count permitted spin states based on m=s coordinate frame, not
whatever "spaces[0]" happens to be
---
 my_pyscf/lassi/lassis.py | 10 +++++++---
 my_pyscf/lassi/states.py |  6 ++++++
 2 files changed, 13 insertions(+), 3 deletions(-)

diff --git a/my_pyscf/lassi/lassis.py b/my_pyscf/lassi/lassis.py
index 85f1d289..75358fcc 100644
--- a/my_pyscf/lassi/lassis.py
+++ b/my_pyscf/lassi/lassis.py
@@ -179,8 +179,8 @@ def all_spin_halfexcitations (lsi, las, nspin=1):
     f1 = f1[None,:,:] - np.tensordot (casdm1s, h2, axes=((1,2),(2,1)))
     i = 0
     auto_singles = isinstance (nspin, str) and 's' in nspin.lower ()
-    nup0 = np.minimum (spaces[0].nelecb, spaces[0].nholea)
-    ndn0 = np.minimum (spaces[0].neleca, spaces[0].nholeb)
+    nup0 = np.minimum (spaces[0].nelecd, spaces[0].nholeu)
+    ndn0 = np.minimum (spaces[0].nelecu, spaces[0].nholed)
     if not auto_singles: # integer supplied by caller
         nup0[:] = nspin
         ndn0[:] = nspin
@@ -214,7 +214,11 @@ def cisolve (sm, nroots):
                       ifrag, nelec, norb, smult-2)
             smults1_i.extend ([smult-2,]*(smult-2))
             spins1_i.extend (list (range (smult-3, -(smult-3)-1, -2)))
-            ci1_i.extend (cisolve (smult-2, ndn0[ifrag]))
+            try:
+                ci1_i.extend (cisolve (smult-2, ndn0[ifrag]))
+            except ValueError as err:
+                print (ndn0[ifrag], nelec, norb, smult)
+                raise (err)
         min_npair = max (0, nelec-norb)
         max_smult = (nelec - 2*min_npair) + 1
         if smult < max_smult: # spin-raised
diff --git a/my_pyscf/lassi/states.py b/my_pyscf/lassi/states.py
index dd4e60a4..c4be84dd 100644
--- a/my_pyscf/lassi/states.py
+++ b/my_pyscf/lassi/states.py
@@ -35,6 +35,12 @@ def __init__(self, las, spins, smults, charges, weight, nlas=None, nelelas=None,
         self.nholea = self.nlas - self.neleca
         self.nholeb = self.nlas - self.nelecb
 
+        # "u", "d": like "a", "b", but presuming spins+1==smults everywhere
+        self.nelecu = (self.nelec + (self.smults-1)) // 2
+        self.nelecd = (self.nelec + (self.smults-1)) // 2
+        self.nholeu = self.nlas - self.nelecu
+        self.nholed = self.nlas - self.nelecd
+
     def __eq__(self, other):
         if self.nfrag != other.nfrag: return False
         return (np.all (self.spins==other.spins) and