Create my_pyscf/lassi/lassis.py

Uses ExcitationPSFCISolver for charge-separated rootspaces,
currently with equal-weights state-averaging. TODO: spin
fluctations and unittests.

my_pyscf/lassi/lassi.py

@@ -677,6 +677,7 @@ def __init__(self, las, mo_coeff=None, ci=None, veff_c=None, h2eff_sub=None, orb
         else: raise RuntimeError("LASSI requires las instance")
         self.__dict__.update(las.__dict__)
         self.ncore = las.ncore
+        self.nfrags = las.nfrags
         keys = set(('e_roots', 'si', 's2', 's2_mat', 'nelec', 'wfnsym', 'rootsym', 'break_symmetry', 'soc', 'opt'))
         self.e_roots = None
         self.si = None

my_pyscf/lassi/lassis.py

@@ -0,0 +1,74 @@
+import numpy as np
+from pyscf.lib import logger
+from mrh.my_pyscf.mcscf.lasci import get_space_info
+from mrh.my_pyscf.mcscf.productstate import ProductStateFCISolver
+from mrh.my_pyscf.lassi.excitations import ExcitationPSFCISolver
+from mrh.my_pyscf.lassi.states import spin_shuffle, spin_shuffle_ci
+from mrh.my_pyscf.lassi.states import all_single_excitations, SingleLASRootspace
+from mrh.my_pyscf.lassi.lassi import LASSI
+
+def prepare_states (lsi, nmax_charge=0):
+    las = lsi._las
+    las1 = spin_shuffle (las)
+    las1.ci = spin_shuffle_ci (las1, las1.ci)
+    # TODO: make states_energy_elec capable of handling lroots and address inconsistency
+    # between definition of e_states array for neutral and charge-separated rootspaces
+    las1.e_states = las1.energy_nuc () + np.array (las1.states_energy_elec ())
+    las2 = all_single_excitations (las1)
+    las2.ci, las2.e_states = single_excitations_ci (lsi, las2, las1, nmax_charge=nmax_charge)
+    return las2
+
+def single_excitations_ci (lsi, las2, las1, nmax_charge=0):
+    log = logger.new_logger (lsi, lsi.verbose)
+    mol = lsi.mol
+    nfrags = lsi.nfrags
+    e_roots = np.append (las1.e_states, np.zeros (las2.nroots-las1.nroots))
+    psrefs = []
+    ci = [[ci_ij for ci_ij in ci_i] for ci_i in las2.ci]
+    for j in range (las1.nroots):
+        solvers = [b.fcisolvers[j] for b in las1.fciboxes]
+        psrefs.append (ProductStateFCISolver (solvers, stdout=mol.stdout, verbose=mol.verbose))
+    spaces = [SingleLASRootspace (las2, m, s, c, 0, ci=[c[ix] for c in ci])
+              for ix, (c, m, s, w) in enumerate (zip (*get_space_info (las2)))]
+    ncsf = las2.get_ugg ().ncsf_sub
+    if isinstance (nmax_charge, np.ndarray): nmax_charge=nmax_charge[None,:]
+    lroots = np.minimum (1+nmax_charge, ncsf)
+    h0, h1, h2 = lsi.ham_2q ()
+    t0 = (logger.process_clock (), logger.perf_counter ())
+    for i in range (las1.nroots, las2.nroots):
+        psref = []
+        ciref = [[] for j in range (nfrags)]
+        excfrags = set ()
+        for j in range (las1.nroots):
+            if not spaces[i].is_single_excitation_of (spaces[j]): continue
+            excfrags = excfrags.union (spaces[i].list_excited_fragments (spaces[j]))
+            psref.append (psrefs[j])
+            for k in range (nfrags):
+                ciref[k].append (las1.ci[k][j])
+        psexc = ExcitationPSFCISolver (psref, ciref, las2.ncas_sub, las2.nelecas_sub,
+                                       stdout=mol.stdout, verbose=mol.verbose)
+        neleca = spaces[i].neleca
+        nelecb = spaces[i].nelecb
+        smults = spaces[i].smults
+        for k in excfrags:
+            weights = np.ones (lroots[k,i]) / lroots[k,i]
+            psexc.set_excited_fragment_(k, (neleca[k],nelecb[k]), smults[k], weights=weights)
+        conv, e_roots[i], ci1 = psexc.kernel (h1, h2, ecore=h0,
+                                              max_cycle_macro=las2.max_cycle_macro,
+                                              conv_tol_self=1)
+        if not conv: log.warn ("CI vectors for charge-separated rootspace %d not converged", i)
+        for k in range (nfrags):
+            ci[k][i] = ci1[k]
+        t0 = log.timer ("Space {} excitations".format (i), *t0)
+    return ci, e_roots
+
+class LASSIS (LASSI):
+    def __init__(self, las, nmax_charge=None, **kwargs):
+        self.nmax_charge = nmax_charge
+        LASSI.__init__(self, las, **kwargs)
+    def kernel (self, nmax_charge=None, **kwargs):
+        if nmax_charge is None: nmax_charge = self.nmax_charge
+        las = prepare_states (self, nmax_charge=nmax_charge)
+        self.__dict__.update(las.__dict__)
+        return LASSI.kernel (self, **kwargs)
+

my_pyscf/lassi/states.py

@@ -188,7 +188,35 @@ def get_ndet (self):
         return [(cistring.num_strings (self.nlas[i], self.neleca[i]),
                  cistring.num_strings (self.nlas[i], self.nelecb[i]))
                 for i in range (self.nfrag)]
-
+
+    def is_single_excitation_of (self, other):
+        # Same charge sector
+        if self.nelec.sum () != other.nelec.sum (): return False
+        # Same spinpol sector
+        if self.spins.sum () != other.spins.sum (): return False
+        # Only 2 fragments involved
+        idx_exc = self.excited_fragments (other)
+        if np.count_nonzero (idx_exc) != 2: return False
+        # Only 1 electron hops
+        dnelec = self.nelec[idx_exc] - other.nelec[idx_exc]
+        if tuple (np.sort (dnelec)) != (-1,1): return False
+        # No additional spin fluctuation between the two excited fragments
+        dspins = self.spins[idx_exc] - other.spins[idx_exc]
+        if tuple (np.sort (dspins)) != (-1,1): return False
+        dsmults = np.abs (self.smults[idx_exc] - other.smults[idx_exc])
+        if np.any (dsmults != 1): return False
+        return True
+
+    def excited_fragments (self, other):
+        dneleca = self.neleca - other.neleca
+        dnelecb = self.nelecb - other.nelecb
+        dsmults = self.smults - other.smults
+        idx_same = (dneleca==0) & (dnelecb==0) & (dsmults==0)
+        return ~idx_same
+
+    def list_excited_fragments (self, other):
+        return np.where (self.excited_fragments (other))[0]
+
 def all_single_excitations (las, verbose=None):
     '''Add states characterized by one electron hopping from one fragment to another fragment
     in all possible ways. Uses all states already present as reference states, so that calling
@@ -255,6 +283,13 @@ def spin_shuffle (las, verbose=None):
     return las.state_average (weights=weights, charges=charges, spins=spins, smults=smults)
 
 def spin_shuffle_ci (las, ci):
+    '''Fill out the CI vectors for rootspaces constructed by the spin_shuffle function.
+    Unallocated CI vectors (None elements in ci) for rootspaces which have the same
+    charge and spin-magnitude strings of rootspaces that do have allocated CI
+    vectors are set to the appropriate rotations of the latter. In the event that
+    more than one reference state for an unallocated rootspace is identified, the
+    rotated vectors are combined and orthogonalized. Unlike running las.lasci (),
+    doing this should ALWAYS guarantee good spin quantum number.'''
     from mrh.my_pyscf.mcscf.lasci import get_space_info
     spaces = [SingleLASRootspace (las, m, s, c, 0, ci=[c[ix] for c in ci])
               for ix, (c, m, s, w) in enumerate (zip (*get_space_info (las)))]

tests/lassi/test_c2h4n4.py

@@ -132,8 +132,10 @@ def test_singles_constructor (self):
     def test_spin_shuffle (self):
         from mrh.my_pyscf.lassi.states import spin_shuffle, spin_shuffle_ci
         mf = lsi._las._scf
-        las3 = LASSCF (mf, (4,2,4), (4,2,4), spin_sub=(5,3,5))
-        las3.lasci ()
+        las3 = LASSCF (mf, (4,2,4), ((4,0),(1,1),(0,4)), spin_sub=(5,3,5))
+        ci_quin = np.array ([[1.0,],])
+        ci_sing = np.diag ([1,-1])[::-1] / np.sqrt (2)
+        las3.ci = [[ci_quin], [ci_sing], [ci_quin.copy ()]]
         las3 = spin_shuffle (las3)
         las3.check_sanity ()
         # The number of states is the number of graphs connecting one number