From 0c404252dddc48990ac50cbbe06c51ad323a3183 Mon Sep 17 00:00:00 2001 From: Matthew R Hermes Date: Fri, 22 Sep 2023 14:07:48 -0500 Subject: [PATCH] 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 | 1 + my_pyscf/lassi/lassis.py | 74 ++++++++++++++++++++++++++++++++++++++ my_pyscf/lassi/states.py | 37 ++++++++++++++++++- tests/lassi/test_c2h4n4.py | 6 ++-- 4 files changed, 115 insertions(+), 3 deletions(-) create mode 100644 my_pyscf/lassi/lassis.py diff --git a/my_pyscf/lassi/lassi.py b/my_pyscf/lassi/lassi.py index 9bd79dbc..b9c1de68 100644 --- a/my_pyscf/lassi/lassi.py +++ b/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 diff --git a/my_pyscf/lassi/lassis.py b/my_pyscf/lassi/lassis.py new file mode 100644 index 00000000..f4e53033 --- /dev/null +++ b/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) + diff --git a/my_pyscf/lassi/states.py b/my_pyscf/lassi/states.py index 2c22a974..5fad0398 100644 --- a/my_pyscf/lassi/states.py +++ b/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)))] diff --git a/tests/lassi/test_c2h4n4.py b/tests/lassi/test_c2h4n4.py index aea6a5a0..5b07c7b8 100644 --- a/tests/lassi/test_c2h4n4.py +++ b/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