Merge branch 'dev' of https://github.com/MatthewRHermes/mrh into gpudev

examples/pbc/lasci_gamma_point.py

@@ -0,0 +1,78 @@
+import numpy
+from pyscf import gto, scf
+from pyscf import mcscf
+from mrh.my_pyscf.mcscf.lasscf_o0 import LASSCF
+from mrh.my_pyscf.tools.molden import from_lasscf
+
+'''
+Here I am doing first gamma point LASCI calculation
+For reference I did molecular calculation, periodic calculation
+with large unit cell (supercell) and at gamma point should be 
+equals to the molecular values
+Outcome of these calculations
+1. HF, CACI, and LASCI is converging to molecular value
+
+Molecular calculations are done without density fitting, but
+for periodic calculations FFTDF is by default on.
+Here i am using GDF
+Probably that's the reason of slow convergence!
+'''
+
+# Molecular Calculation
+mol = gto.M(atom = '''
+    H         -6.37665        2.20769        0.00000
+    H         -5.81119        2.63374       -0.00000
+    H         -6.37665        2.20769        2.00000
+    H         -5.81119        2.63374        2.00000
+    ''',
+    basis = '631g',
+    verbose = 4,max_memory=10000)
+
+mf = scf.RHF(mol)
+memf = mf.kernel()
+
+mc = mcscf.CASCI(mf, 4, 4)
+mecasci = mc.kernel()[0]
+
+las = LASSCF(mf, (2,2), (2,2))
+mo0 = las.localize_init_guess((list(range(2)), list(range(2, 4))))
+melasci = las.lasci(mo0)
+
+del mol, mf, mc, las
+
+# Periodic Calculation
+from pyscf.pbc import gto, scf
+
+def cellObject(x):
+    cell = gto.M(a = numpy.eye(3)*x,
+    atom = '''
+    H         -6.37665        2.20769        0.00000
+    H         -5.81119        2.63374       -0.00000
+    H         -6.37665        2.20769        2.00000
+    H         -5.81119        2.63374        2.00000
+    ''',
+    basis = '631g',
+    verbose = 1, max_memory=10000)
+    cell.build()
+
+    mf = scf.RHF(cell).density_fit() # GDF
+    mf.exxdiv = None
+    emf = mf.kernel()
+
+    mc = mcscf.CASCI(mf, 4, 4)
+    ecasci = mc.kernel()[0]
+
+    las = LASSCF(mf, (2,2), (2,2))
+    mo0 = las.localize_init_guess((list(range(2)),list(range(2, 4))))
+    elasci = las.lasci(mo0)
+
+    del cell, mc, mf, las
+    return x, emf, ecasci, elasci
+
+print(" Energy Comparision with cubic unit cell size ")
+print(f"{'LatticeVector(a)':<20} {'HF':<20} {'CASCI':<15} {'LASCI':<20}")
+print(f"{'Reference':<18} {memf:<18.9f} {mecasci:<18.9f} {melasci[1]:<18.9f}")
+
+for x in range(7,10,1):
+    x, emf, ecasci, elasci = cellObject(x)
+    print(f"{x:<18.1f} {emf:<18.9f} {ecasci:<18.9f} {elasci[1]:<18.9f}")

my_pyscf/lassi/chkfile.py

@@ -2,22 +2,23 @@
 
 KEYS_CONFIG_LASSI = las_chkfile.KEYS_CONFIG_LASSCF + ['nfrags', 'break_symmetry', 'soc', 'opt']
 KEYS_SACONSTR_LASSI = las_chkfile.KEYS_SACONSTR_LASSCF
-KEYS_RESULTS_LASSI = ['e_roots', 'si', 's2', 's2_mat', 'nelec', 'wfnsym', 'rootsym']
+KEYS_RESULTS_LASSI = ['e_states', 'e_roots', 'si', 's2', 's2_mat', 'nelec', 'wfnsym', 'rootsym']
 
 def load_lsi_(lsi, chkfile=None, method_key='lsi',
               keys_config=KEYS_CONFIG_LASSI,
               keys_saconstr=KEYS_SACONSTR_LASSI,
               keys_results=KEYS_RESULTS_LASSI):
+    lsi._las.load_chk (chkfile=chkfile)
     return las_chkfile.load_las_(lsi, chkfile=chkfile, method_key=method_key,
                                  keys_config=keys_config,
                                  keys_saconstr=keys_saconstr, keys_results=keys_results)
 
-
-def dump_lsi (lsi, chkfile='None', method_key='lsi', mo_coeff=None, ci=None,
+def dump_lsi (lsi, chkfile=None, method_key='lsi', mo_coeff=None, ci=None,
               overwrite_mol=True, keys_config=KEYS_CONFIG_LASSI,
               keys_saconstr=KEYS_SACONSTR_LASSI,
               keys_results=KEYS_RESULTS_LASSI,
               **kwargs):
+    lsi._las.dump_chk (chkfile=chkfile)
     return las_chkfile.dump_las (lsi, chkfile=chkfile, method_key=method_key, mo_coeff=mo_coeff,
                                  ci=ci, overwrite_mol=overwrite_mol, keys_config=keys_config,
                                  keys_saconstr=keys_saconstr, keys_results=keys_results, **kwargs)

my_pyscf/lassi/lassi.py

@@ -3,8 +3,9 @@
 from scipy import linalg
 from mrh.my_pyscf.lassi import op_o0
 from mrh.my_pyscf.lassi import op_o1
+from mrh.my_pyscf.lassi import chkfile
 from mrh.my_pyscf.lassi.citools import get_lroots
-from pyscf import lib, symm
+from pyscf import lib, symm, ao2mo
 from pyscf.scf.addons import canonical_orth_
 from pyscf.lib.numpy_helper import tag_array
 from pyscf.fci.direct_spin1 import _unpack_nelec
@@ -80,8 +81,6 @@ def ham_2q (las, mo_coeff, veff_c=None, h2eff_sub=None, soc=0):
     if veff_c is None: 
         dm_core = 2 * mo_core @ mo_core.conj ().T
         veff_c = las.get_veff (dm1s=dm_core)
-    if h2eff_sub is None:
-        h2eff_sub = las.ao2mo (mo_coeff)
 
     h0 = las._scf.energy_nuc () + 2 * (((hcore + veff_c/2) @ mo_core) * mo_core).sum ()
 
@@ -97,8 +96,11 @@ def ham_2q (las, mo_coeff, veff_c=None, h2eff_sub=None, soc=0):
         h1[0:ncas,0:ncas] += hso[2] # a'a
         h1[ncas:2*ncas,ncas:2*ncas] -= hso[2] # b'b
 
-    h2 = h2eff_sub[ncore:nocc].reshape (ncas*ncas, ncas * (ncas+1) // 2)
-    h2 = lib.numpy_helper.unpack_tril (h2).reshape (ncas, ncas, ncas, ncas)
+    if h2eff_sub is None:
+        h2 = las.get_h2cas (mo_coeff)
+    else:
+        h2 = h2eff_sub[ncore:nocc].reshape (ncas*ncas, ncas * (ncas+1) // 2)
+        h2 = lib.numpy_helper.unpack_tril (h2).reshape (ncas, ncas, ncas, ncas)
 
     return h0, h1, h2
 
@@ -248,7 +250,10 @@ def lassi (las, mo_coeff=None, ci=None, veff_c=None, h2eff_sub=None, orbsym=None
         raise RuntimeError ('Insufficient memory to use o0 LASSI algorithm')
 
     # Construct second-quantization Hamiltonian
-    e0, h1, h2 = ham_2q (las, mo_coeff, veff_c=veff_c, h2eff_sub=h2eff_sub, soc=soc)
+    if callable (getattr (las, 'ham_2q', None)):
+        e0, h1, h2 = las.ham_2q (mo_coeff, veff_c=veff_c, h2eff_sub=h2eff_sub, soc=soc)
+    else:
+        e0, h1, h2 = ham_2q (las, mo_coeff, veff_c=veff_c, h2eff_sub=h2eff_sub, soc=soc)
 
     # Symmetry tuple: neleca, nelecb, irrep
     statesym, s2_states = las_symm_tuple (las, break_spin=soc, break_symmetry=break_symmetry)
@@ -582,8 +587,10 @@ def roots_make_rdm12s (las, ci, si, orbsym=None, soc=None, break_symmetry=None,
             orbsym = las.label_symmetry_(las.mo_coeff).orbsym
         if orbsym is not None:
             orbsym = orbsym[las.ncore:las.ncore+las.ncas]
-    if soc is None: soc = si.soc
-    if break_symmetry is None: break_symmetry = si.break_symmetry
+    if soc is None:
+        soc = getattr (si, 'soc', getattr (las, 'soc', False))
+    if break_symmetry is None:
+        break_symmetry = getattr (si, 'break_symmetry', getattr (las, 'break_symmetry', False))
     o0_memcheck = op_o0.memcheck (las, ci, soc=soc)
     if opt == 0 and o0_memcheck == False:
         raise RuntimeError ('Insufficient memory to use o0 LASSI algorithm')
@@ -681,8 +688,10 @@ def root_make_rdm12s (las, ci, si, state=0, orbsym=None, soc=None, break_symmetr
     '''
     states = np.atleast_1d (state)
     si_column = si[:,states]
-    if soc is None: soc = si.soc
-    if break_symmetry is None: break_symmetry = si.break_symmetry
+    if soc is None:
+        soc = getattr (si, 'soc', getattr (las, 'soc', False))
+    if break_symmetry is None:
+        break_symmetry = getattr (si, 'break_symmetry', getattr (las, 'break_symmetry', False))
     if rootsym is None:
         rootsym = getattr (si, 'rootsym', getattr (las, 'rootsym', None))
     rootsym = [rootsym[s] for s in states]
@@ -710,6 +719,7 @@ def __init__(self, las, mo_coeff=None, ci=None, soc=False, break_symmetry=False,
         self.ncore, self.ncas = las.ncore, las.ncas
         self.nfrags, self.nroots = las.nfrags, las.nroots
         self.ncas_sub, self.nelecas_sub, self.fciboxes = las.ncas_sub, las.nelecas_sub, las.fciboxes
+        self.nelecas = sum (self.nelecas_sub)
         self.weights, self.e_states, self.e_lexc = las.weights, las.e_states, las.e_lexc
         self.converged = las.converged
         # I/O data from las parent
@@ -752,13 +762,13 @@ def ham_2q (self, mo_coeff=None, veff_c=None, h2eff_sub=None, soc=0):
 
     def get_nelec_frs (self, las=None):
         if las is None: las = self
-        nelec_frs = []
-        for fcibox, nelec in zip (las.fciboxes, las.nelecas_sub):
-            nelec_rs = []
-            for solver in fcibox.fcisolvers:
-                nelec_rs.append (_unpack_nelec (fcibox._get_nelec (solver, nelec)))
-            nelec_frs.append (nelec_rs)
-        return np.asarray (nelec_frs)
+        from mrh.my_pyscf.mcscf.lasci import get_nelec_frs
+        return get_nelec_frs (las)
+
+    def get_smult_fr (self, las=None):
+        if las is None: las = self
+        from mrh.my_pyscf.mcscf.lasci import get_space_info
+        return get_space_info (las)[2].T
 
     def get_lroots (self, ci=None):
         if ci is None: ci = self.ci
@@ -800,12 +810,15 @@ def get_sivec_vacuum_shuffle (self, state=None, nelec_vac=None, si=None, ci=None
         nelec_frs = self.get_nelec_frs ()
         return sivec_vacuum_shuffle (si, nelec_frs, lroots, nelec_vac=nelec_vac, state=state)
 
-    def analyze (self, state=None, **kwargs):
+    def analyze (self, state=0, **kwargs):
         from mrh.my_pyscf.lassi.sitools import analyze
-        analyze (self, self.si, state=state, **kwargs)
+        return analyze (self, self.si, state=state, **kwargs)
 
     def reset (self, mol=None):
         if mol is not None:
             self.mol = mol
         self._las.reset (mol)
 
+    dump_chk = chkfile.dump_lsi
+    load_chk = load_chk_ = chkfile.load_lsi_
+

my_pyscf/lassi/lassis.py

@@ -31,6 +31,7 @@ def prepare_states (lsi, ncharge=1, nspin=0, sa_heff=True, deactivate_vrv=False,
     # TODO: make states_energy_elec capable of handling lroots and address inconsistency
     # between definition of e_states array for neutral and charge-separated rootspaces
     t0 = (logger.process_clock (), logger.perf_counter ())
+    ham_2q = lsi.ham_2q ()
     log = logger.new_logger (lsi, lsi.verbose)
     las = lsi._las.get_single_state_las (state=0)
     # 1. Spin shuffle step
@@ -39,6 +40,7 @@ def prepare_states (lsi, ncharge=1, nspin=0, sa_heff=True, deactivate_vrv=False,
         las1 = las
     else:
         las1 = spin_shuffle (las, equal_weights=True)
+        # TODO: memory efficiency; the line below makes copies
         las1.ci = spin_shuffle_ci (las1, las1.ci)
         las1.converged = las.converged
     nroots_ref = las1.nroots
@@ -50,7 +52,7 @@ def prepare_states (lsi, ncharge=1, nspin=0, sa_heff=True, deactivate_vrv=False,
         log.info ("Reference space %d:", ix)
         SingleLASRootspace (las1, m, s, c, 0, ci=[c[ix] for c in las1.ci]).table_printlog ()
     # 2. Spin excitations part 1
-    spin_flips = all_spin_flips (lsi, las1, nspin=nspin) if nspin else None
+    spin_flips = all_spin_flips (lsi, las1, nspin=nspin, ham_2q=ham_2q) if nspin else None
     las1.e_states = las1.energy_nuc () + np.array (las1.states_energy_elec ())
     # 3. Charge excitations
     # TODO: Store the irreducible degrees of freedom of the charge excitations more transparently,
@@ -59,7 +61,7 @@ def prepare_states (lsi, ncharge=1, nspin=0, sa_heff=True, deactivate_vrv=False,
         las2 = all_single_excitations (las1)
         converged, spaces2 = single_excitations_ci (
             lsi, las2, las1, ncharge=ncharge, sa_heff=sa_heff, deactivate_vrv=deactivate_vrv,
-            spin_flips=spin_flips, crash_locmin=crash_locmin
+            spin_flips=spin_flips, crash_locmin=crash_locmin, ham_2q=ham_2q
         )
     else:
         converged = las1.converged
@@ -84,7 +86,7 @@ def prepare_states (lsi, ncharge=1, nspin=0, sa_heff=True, deactivate_vrv=False,
     return converged, las3
 
 def single_excitations_ci (lsi, las2, las1, ncharge=1, sa_heff=True, deactivate_vrv=False,
-                           spin_flips=None, crash_locmin=False):
+                           spin_flips=None, crash_locmin=False, ham_2q=None):
     log = logger.new_logger (lsi, lsi.verbose)
     mol = lsi.mol
     nfrags = lsi.nfrags
@@ -102,7 +104,10 @@ def single_excitations_ci (lsi, las2, las1, ncharge=1, sa_heff=True, deactivate_
         else:
             raise RuntimeError ("Valid ncharge values are integers or 's'")
     lroots = np.minimum (ncharge, ncsf)
-    h0, h1, h2 = lsi.ham_2q ()
+    if ham_2q is None:
+        h0, h1, h2 = lsi.ham_2q ()
+    else:
+        h0, h1, h2 = ham_2q
     t0 = (logger.process_clock (), logger.perf_counter ())
     converged = True
     # Prefilter spin-shuffles
@@ -184,7 +189,7 @@ def __init__(self, ci, spins, smults):
         self.spins = spins
         self.smults = smults
 
-def all_spin_flips (lsi, las, nspin=1):
+def all_spin_flips (lsi, las, nspin=1, ham_2q=None):
     # NOTE: this actually only uses the -first- rootspace in las, so it can be done before
     # the initial spin shuffle
     log = logger.new_logger (lsi, lsi.verbose)
@@ -202,7 +207,10 @@ def all_spin_flips (lsi, las, nspin=1):
     smults1 = []
     spins1 = []
     ci1 = []
-    h0, h1, h2 = lsi.ham_2q ()
+    if ham_2q is None:
+        h0, h1, h2 = lsi.ham_2q ()
+    else:
+        h0, h1, h2 = ham_2q
     casdm1s = las.make_casdm1s ()
     f1 = h1 + np.tensordot (h2, casdm1s.sum (0), axes=2)
     f1 = f1[None,:,:] - np.tensordot (casdm1s, h2, axes=((1,2),(2,1)))
@@ -328,6 +336,8 @@ def _spin_shuffle_ci_(spaces, spin_flips, nroots_ref, nroots_refc):
         ci_szrot_1c.append (space.get_ci_szrot (ifrags=(ifrag,jfrag)))
     charges0 = spaces[0].charges
     smults0 = spaces[0].smults
+    # Prepare reference szrots
+    ci_szrot_ref = spaces[0].get_ci_szrot ()
     for ix in new_idx:
         idx = spaces[ix].excited_fragments (spaces[0])
         space = spaces[ix]
@@ -343,9 +353,7 @@ def _spin_shuffle_ci_(spaces, spin_flips, nroots_ref, nroots_refc):
                 iflp = np.where (iflp)[0][0]
                 space.ci[ifrag] = sf.ci[iflp]
             else: # Reference-state spin-shuffles
-                spaces0 = [sp for sp in spaces[:nroots_ref] if sp.spins[ifrag]==space.spins[ifrag]]
-                assert (len (spaces0))
-                space.ci[ifrag] = spaces0[0].ci[ifrag]
+                space.ci[ifrag] = ci_szrot_ref[ifrag][space.spins[ifrag]]
         for (ci_i, ci_j), sp_1c in zip (ci_szrot_1c, spaces_1c):
             ijfrag = sp_1c.entmap[0]
             if ijfrag not in spaces[ix].entmap: continue
@@ -484,10 +492,15 @@ def kernel (self, ncharge=None, nspin=None, sa_heff=None, deactivate_vrv=None,
                 crash_locmin=None, **kwargs):
         t0 = (logger.process_clock (), logger.perf_counter ())
         log = logger.new_logger (self, self.verbose)
-        self.converged = self.prepare_states_(ncharge=ncharge, nspin=nspin,
-                                              sa_heff=sa_heff, deactivate_vrv=deactivate_vrv,
-                                              crash_locmin=crash_locmin)
-        self.e_roots, self.si = self.eig (**kwargs)
+        h0, h1, h2 = self.ham_2q ()
+        t1 = log.timer ("LASSIS integral transformation", *t0)
+        with lib.temporary_env (self, ham_2q=lambda *args, **kwargs: (h0, h1, h2)):
+            self.converged = self.prepare_states_(ncharge=ncharge, nspin=nspin,
+                                                  sa_heff=sa_heff, deactivate_vrv=deactivate_vrv,
+                                                  crash_locmin=crash_locmin)
+            t1 = log.timer ("LASSIS state preparation", *t1)
+            self.e_roots, self.si = self.eig (**kwargs)
+            t1 = log.timer ("LASSIS diagonalization", *t1)
         log.timer ("LASSIS", *t0)
         return self.e_roots, self.si
 

my_pyscf/lassi/op_o0.py

@@ -144,6 +144,25 @@ def ci1_r_gen (buf=None):
     return ci1_r_gen, ss2spinless, spinless2ss
 
 def civec_spinless_repr (ci0_r, norb, nelec_r):
+    '''Put CI vectors in the spinless representation; i.e., map
+        norb -> 2 * norb
+        (neleca, nelecb) -> (neleca+nelecb, 0)
+    This permits linear combinations of CI vectors with different
+    M == neleca-nelecb at the price of higher memory cost. This function
+    does NOT change the datatype.
+
+    Args:
+        ci0_r: sequence or generator of ndarray of length nprods
+            CAS-CI vectors in the spin-pure representation
+        norb: integer
+            Number of orbitals
+        nelec_r: sequence of tuple of length (2)
+            (neleca, nelecb) for each element of ci0_r
+
+    Returns:
+        ci1_r: ndarray of shape (nprods, ndet_spinless)
+            spinless CAS-CI vectors
+    '''
     ci1_r_gen, _, _ = civec_spinless_repr_generator (ci0_r, norb, nelec_r)
     ci1_r = np.stack ([x.copy () for x in ci1_r_gen ()], axis=0)
     return ci1_r