Merge branch 'dev' into lassis_ncharge_rationalization

debug/lassi/debug_22.py

@@ -20,7 +20,7 @@
 from mrh.my_pyscf.mcscf.lasscf_o0 import LASSCF
 from mrh.my_pyscf.lassi import LASSI
 from mrh.my_pyscf.lassi.lassi import root_make_rdm12s, make_stdm12s
-from mrh.my_pyscf.lassi.states import all_single_excitations, SingleLASRootspace
+from mrh.my_pyscf.lassi.spaces import all_single_excitations, SingleLASRootspace
 from mrh.my_pyscf.mcscf.lasci import get_space_info
 from mrh.my_pyscf.lassi import op_o0, op_o1, lassis
 

debug/lassi/debug_c2h4n4.py

@@ -122,7 +122,7 @@ class KnownValues(unittest.TestCase):
     #        self.assertAlmostEqual (e1, e0, 8)
 
     #def test_singles_constructor (self):
-    #    from mrh.my_pyscf.lassi.states import all_single_excitations
+    #    from mrh.my_pyscf.lassi.spaces import all_single_excitations
     #    las2 = all_single_excitations (lsi._las)
     #    las2.check_sanity ()
     #    # Meaning of tuple: (na+nb,smult)
@@ -135,7 +135,7 @@ class KnownValues(unittest.TestCase):
     #    self.assertEqual (las2.nroots, 33)
 
     #def test_spin_shuffle (self):
-    #    from mrh.my_pyscf.lassi.states import spin_shuffle, spin_shuffle_ci
+    #    from mrh.my_pyscf.lassi.spaces import spin_shuffle, spin_shuffle_ci
     #    mf = lsi._las._scf
     #    las3 = spin_shuffle (las)
     #    las3.check_sanity ()

debug/lassi/debug_lassis_targets_slow.py

@@ -186,7 +186,7 @@ def test_kremer_cr2_model (self):
             mo_list = mc_avas.ncore + np.array ([5,6,7,8,9,10,15,16,17,18,19,20])
             mo_coeff = las.sort_mo (mo_list, mo_coeff)
             mo_coeff = las.localize_init_guess (([0],[1]), mo_coeff)
-        las = lassi.states.spin_shuffle (las) # generate direct-exchange states
+        las = lassi.spaces.spin_shuffle (las) # generate direct-exchange states
         las.weights = [1.0/las.nroots,]*las.nroots # set equal weights
         nroots_ref = las.nroots
         las.kernel (mo_coeff) # optimize orbitals

examples/laspdft/c2h4n4_si_laspdft.py

@@ -5,7 +5,7 @@
 from mrh.my_pyscf import mcpdft
 from mrh.my_pyscf.tools.molden import from_lasscf
 from c2h4n4_struct import structure as struct
-from mrh.my_pyscf.lassi.states import all_single_excitations
+from mrh.my_pyscf.lassi.spaces import all_single_excitations
 
 # Mean field calculation
 mol = struct(0, 0, '6-31g')
@@ -23,7 +23,7 @@
 mo0 = las.localize_init_guess((list(range (5)), list(range (5,10))), guess_mo)
 las.kernel(mo0)
 
-las = lassi.states.all_single_excitations (las)
+las = lassi.spaces.all_single_excitations (las)
 las.lasci ()
 
 lsi = lassi.LASSI(las)

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}")

examples/sa_si_lasscf/automatic_singles.py

@@ -44,7 +44,7 @@
 print ("LASSI(hand) energy =", e_roots[0])
 molden.from_lassi (las2, 'c2h4n4_las66si4_631g.molden', si=si_hand)
 
-from mrh.my_pyscf.lassi.states import all_single_excitations
+from mrh.my_pyscf.lassi.spaces import all_single_excitations
 las = all_single_excitations (las)
 las.lasci () # Optimize the CI vectors
 las.dump_spaces () # prints all state tables in the output file

examples/sa_si_lasscf/cr2o3n6h21_lassis66_smallbasis.py

@@ -46,12 +46,12 @@ def yamaguchi (e_roots, s2):
     mo_coeff = las.localize_init_guess (([0],[1]), mo_coeff)
 
 # Direct exchange only result
-las = lassi.states.spin_shuffle (las) # generate direct-exchange states
+las = lassi.spaces.spin_shuffle (las) # generate direct-exchange states
 las.weights = [1.0/las.nroots,]*las.nroots # set equal weights
 las.kernel (mo_coeff) # optimize orbitals
 lsi0 = lassi.LASSI (las).run ()
 print (("Direct exchange only is modeled by {} states constructed with "
-        "lassi.states.spin_shuffle").format (las.nroots))
+        "lassi.spaces.spin_shuffle").format (las.nroots))
 print ("J(LASSI, direct) = %.2f cm^-1" % yamaguchi (lsi0.e_roots, lsi0.s2))
 print ("{} rootspaces and {} product states total\n".format (lsi0.nroots, lsi0.si.shape[1]))
 
@@ -70,8 +70,8 @@ def yamaguchi (e_roots, s2):
     mc.ci.shape[0], mc.ci.shape[1], mc.ci.size))
 
 # Direct exchange & kinetic exchange result
-las = lassi.states.all_single_excitations (las) # generate kinetic-exchange states
-print (("Use of lassi.states.all_single_excitations generates "
+las = lassi.spaces.all_single_excitations (las) # generate kinetic-exchange states
+print (("Use of lassi.spaces.all_single_excitations generates "
         "{} additional kinetic-exchange (i.e., charge-transfer) "
         "states").format (las.nroots-4))
 las.lasci () # do not reoptimize orbitals at this step - not likely to converge

exploratory/citools/fockspace.py

@@ -78,11 +78,11 @@ def fermion_spin_shuffle (norb, norb_f):
     nelec_f = np.zeros ((ndet, nfrag), dtype=np.int8)
     for ifrag, n in enumerate (ndet_f):
         addrs, fragaddrs = np.divmod (addrs, n)
-        nelec_f[:,ifrag] = ADDRS_NELEC[fragaddrs]
+        nelec_f[:,ifrag] = ADDRS_NELEC[fragaddrs.astype (int)]
     ne_f, ne_f_idx = np.unique (nelec_f, return_inverse=True, axis=0)
     for (ia, na_f), (ib, nb_f) in product (enumerate (ne_f), repeat=2):
-        idx_a = ne_f_idx == ia
-        idx_b = ne_f_idx == ib
+        idx_a = np.squeeze (ne_f_idx == ia)
+        idx_b = np.squeeze (ne_f_idx == ib)
         idx = np.ix_(idx_a,idx_b)
         sgn[idx] = _fss (na_f, nb_f)
     assert (np.count_nonzero (sgn==0) == 0)
@@ -127,7 +127,7 @@ def fermion_frag_shuffle (norb, i, j):
     # Lower orbital indices change faster than higher ones
     addrs_p = np.arange (ndet, dtype=np.uint64) // (2**i)
     addrs_q, addrs_p = np.divmod (addrs_p, 2**(j-i))
-    nperms = ADDRS_NELEC[addrs_p] * ADDRS_NELEC[addrs_q]
+    nperms = ADDRS_NELEC[addrs_p.astype (int)] * ADDRS_NELEC[addrs_q.astype (int)]
     sgn = np.array ([1,-1], dtype=np.int8)[nperms%2]
     assert (sgn.size==ndet)
     return sgn

exploratory/unitary_cc/uccsd_sym1.py

@@ -46,7 +46,7 @@ def spincases (p_idxs, norb):
     m = np.array ([0])
     for ielec in range (nelec):
         q_idxs = p_idxs.copy ()
-        q_idxs[:,ielec] += norb
+        q_idxs[:,ielec] += np.array ([norb], dtype=q_idxs.dtype)
         p_idxs = np.append (p_idxs, q_idxs, axis=0)
         m = np.append (m, m+1)
     p_sorted = np.stack ([np.sort (prow) for prow in p_idxs], axis=0)

my_pyscf/fci/csfstring.py

@@ -563,7 +563,7 @@ def _transform_det2csf (inparr, norb, neleca, nelecb, smult, reverse=False, csd_
     ncol_out = (ncsf_all, ndet_all)[reverse or project]
     ncol_in = (ncsf_all, ndet_all)[~reverse or project]
     if not project:
-        outarr = np.ascontiguousarray (np.zeros ((nrow, ncol_out), dtype=np.float_))
+        outarr = np.ascontiguousarray (np.zeros ((nrow, ncol_out), dtype=np.float64))
         csf_addrs = np.zeros (ncsf_all, dtype=np.bool_)
     # Initialization is necessary because not all determinants have a csf for all spin states
 
@@ -874,7 +874,7 @@ def get_spin_evecs (nspin, neleca, nelecb, smult):
     scstrs = addrs2str (nspin, smult, list (range (ncsf)))
     assert (len (scstrs) == ncsf), "should have {} coupling strings; have {} (nspin={}, s={})".format (ncsf, len (scstrs), nspin, s)
 
-    umat = np.ones ((ndet, ncsf), dtype=np.float_)
+    umat = np.ones ((ndet, ncsf), dtype=np.float64)
     twoS = smult-1
     twoMS = neleca - nelecb
 
@@ -904,7 +904,7 @@ def test_spin_evecs (nspin, neleca, nelecb, smult, S2mat=None):
     spinstrs = cistring.addrs2str (nspin, na, list (range (ndet)))
 
     if S2mat is None:
-        S2mat = np.zeros ((ndet, ndet), dtype=np.float_)
+        S2mat = np.zeros ((ndet, ndet), dtype=np.float64)
         twoS = smult - 1
         twoMS = int (round (2 * ms))
 

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/citools.py

@@ -48,20 +48,64 @@ def get_rootaddr_fragaddr (lroots):
             The ordinal designation local to each fragment of each LAS state.
     '''
     nfrags, nroots = lroots.shape
-    nprods = np.product (lroots, axis=0)
+    nprods = np.prod (lroots, axis=0)
     fragaddr = np.zeros ((nfrags, sum(nprods)), dtype=int)
     rootaddr = np.zeros (sum(nprods), dtype=int)
     offs = np.cumsum (nprods)
     for iroot in range (nroots):
         j = offs[iroot]
         i = j - nprods[iroot]
         for ifrag in range (nfrags):
-            prods_before = np.product (lroots[:ifrag,iroot], axis=0)
-            prods_after = np.product (lroots[ifrag+1:,iroot], axis=0)
+            prods_before = np.prod (lroots[:ifrag,iroot], axis=0)
+            prods_after = np.prod (lroots[ifrag+1:,iroot], axis=0)
             addrs = np.repeat (np.arange (lroots[ifrag,iroot]), prods_before)
             addrs = np.tile (addrs, prods_after)
             fragaddr[ifrag,i:j] = addrs
         rootaddr[i:j] = iroot
     return rootaddr, fragaddr
 
+def umat_dot_1frag_(target, umat, lroots, ifrag, iroot, axis=0):
+    '''Apply a unitary transformation for 1 fragment in 1 rootspace to a tensor
+    whose target axis spans all model states.
 
+    Args:
+        target: ndarray whose length on axis 'axis' is nstates
+            The object to which the unitary transformation is to be applied.
+            Modified in-place.
+        umat: ndarray of shape (lroots[ifrag,iroot],lroots[ifrag,iroot])
+            A unitary matrix; the row axis is contracted
+        lroots: ndarray of shape (nfrags, nroots)
+            Number of basis states in each fragment in each rootspace
+        ifrag: integer
+            Fragment index for the targeted block
+        iroot: integer
+            Rootspace index for the targeted block
+
+    Kwargs:
+        axis: integer
+            The axis of target to which umat is to be applied
+
+    Returns:
+        target: same as input target
+            After application of unitary transformation'''
+    nprods = np.prod (lroots, axis=0)
+    offs = [0,] + list (np.cumsum (nprods))
+    i, j = offs[iroot], offs[iroot+1]
+    newaxes = [axis,] + list (range (axis)) + list (range (axis+1, target.ndim))
+    oldaxes = list (np.argsort (newaxes))
+    target = target.transpose (*newaxes)
+    target[i:j] = _umat_dot_1frag (target[i:j], umat, lroots[:,iroot], ifrag)
+    target = target.transpose (*oldaxes)
+    return target
+
+def _umat_dot_1frag (target, umat, lroots, ifrag):
+    # Remember: COLUMN-MAJOR ORDER!!
+    old_shape = target.shape
+    new_shape = tuple (lroots[::-1]) + old_shape[1:]
+    target = target.reshape (*new_shape)
+    iifrag = len (lroots) - ifrag - 1
+    newaxes = [iifrag,] + list (range (iifrag)) + list (range (iifrag+1, target.ndim))
+    oldaxes = list (np.argsort (newaxes))
+    target = target.transpose (*newaxes)
+    target = np.tensordot (umat.T, target, axes=1).transpose (*oldaxes)
+    return target.reshape (*old_shape)