diff --git a/debug/lassi/debug_4frag.py b/debug/lassi/debug_4frag.py index 50d5ba9e..f02c82d7 100644 --- a/debug/lassi/debug_4frag.py +++ b/debug/lassi/debug_4frag.py @@ -18,7 +18,7 @@ import numpy as np from scipy import linalg from copy import deepcopy -from itertools import product, combinations_with_replacement +from itertools import product from pyscf import lib, gto, scf, dft, fci, mcscf, df from pyscf.tools import molden from pyscf.fci import cistring @@ -29,13 +29,14 @@ from mrh.my_pyscf.lassi.lassi import make_stdm12s, ham_2q, las_symm_tuple from mrh.my_pyscf.lassi import op_o0 from mrh.my_pyscf.lassi import op_o1 -from mrh.my_pyscf.lassi import op_o2 from mrh.my_pyscf.lassi import LASSIS -from mrh.my_pyscf.lassi.op_o2 import get_fdm1_maker +from mrh.my_pyscf.lassi.op_o1 import get_fdm1_maker from mrh.my_pyscf.lassi.sitools import make_sdm1 +from mrh.tests.lassi.addons import case_contract_hlas_ci, case_lassis_fbf_2_model_state +from mrh.tests.lassi.addons import case_lassis_fbfdm def setUpModule (): - global mol, mf, las, nroots, nelec_frs, si, lroots + global mol, mf, las, nroots, nelec_frs, si # State list contains a couple of different 4-frag interactions states = {'charges': [[0,0,0,0],[1,1,-1,-1],[2,-1,-1,0],[1,0,0,-1],[1,1,-1,-1],[2,-1,-1,0],[1,0,0,-1]], 'spins': [[0,0,0,0],[1,1,-1,-1],[0,1,-1,0], [1,0,0,-1],[-1,-1,1,1],[0,-1,1,0], [-1,0,0,1]], @@ -63,8 +64,8 @@ def setUpModule (): 1 -3.206320000 -3.233120000 0.000000000''' mol = gto.M (atom = xyz, basis='STO-3G', symmetry=False, - #verbose=5, output='test_4frag.log') - verbose=0, output='/dev/null') + verbose=5, output='test_4frag.log') + #verbose=0, output='/dev/null') mf = scf.RHF (mol).run () las = LASSCF (mf, (2,2,2,2),((1,1),(1,1),(1,1),(1,1))) las.state_average_(weights=weights, **states) @@ -73,7 +74,6 @@ def setUpModule (): las.mo_coeff = las.localize_init_guess (frags, mf.mo_coeff) las.ci = las.get_init_guess_ci (las.mo_coeff, las.get_h2eff (las.mo_coeff)) lroots = np.minimum (2, las.get_ugg ().ncsf_sub) - lroots[:] = 1 nelec_frs = np.array ( [[_unpack_nelec (fcibox._get_nelec (solver, nelecas)) for solver in fcibox.fcisolvers] for fcibox, nelecas in zip (las.fciboxes, las.nelecas_sub)] @@ -107,111 +107,95 @@ def setUpModule (): si = lib.tag_array (si, rootsym=las_symm_tuple (las)[0]) def tearDownModule(): - global mol, mf, las, nroots, nelec_frs, si, lroots + global mol, mf, las, nroots, nelec_frs, si mol.stdout.close () - del mol, mf, las, nroots, nelec_frs, si, lroots + del mol, mf, las, nroots, nelec_frs, si class KnownValues(unittest.TestCase): - #def test_stdm12s (self): - # d12_o0 = make_stdm12s (las, opt=0) - # d12_o1 = make_stdm12s (las, opt=1) - # for r in range (2): - # for i, j in product (range (si.shape[0]), repeat=2): - # with self.subTest (rank=r+1, bra=i, ket=j): - # self.assertAlmostEqual (lib.fp (d12_o0[r][i,...,j]), - # lib.fp (d12_o1[r][i,...,j]), 9) + def test_stdm12s (self): + d12_o0 = make_stdm12s (las, opt=0) + d12_o1 = make_stdm12s (las, opt=1) + for r in range (2): + for i, j in product (range (si.shape[0]), repeat=2): + with self.subTest (rank=r+1, bra=i, ket=j): + self.assertAlmostEqual (lib.fp (d12_o0[r][i,...,j]), + lib.fp (d12_o1[r][i,...,j]), 9) def test_ham_s2_ovlp (self): - nprods = np.prod (lroots, axis=0) - loff1 = np.cumsum (nprods) - loff0 = loff1 - nprods h1, h2 = ham_2q (las, las.mo_coeff, veff_c=None, h2eff_sub=None)[1:] - h1[:] = 0 lbls = ('ham','s2','ovlp') mats_o0 = op_o0.ham (las, h1, h2, las.ci, nelec_frs)#, orbsym=orbsym, wfnsym=wfnsym) - np.save ('ham_test.npy', mats_o0[0]) fps_o0 = [lib.fp (mat) for mat in mats_o0] mats_o1 = op_o1.ham (las, h1, h2, las.ci, nelec_frs)#, orbsym=orbsym, wfnsym=wfnsym) for lbl, mat, fp in zip (lbls, mats_o1, fps_o0): with self.subTest(opt=1, matrix=lbl): self.assertAlmostEqual (lib.fp (mat), fp, 9) - mats_o2 = op_o2.ham (las, h1, h2, las.ci, nelec_frs)#, orbsym=orbsym, wfnsym=wfnsym) - np.save ('ham_ref.npy', mats_o2[0]) - for lbl, mt, mr in zip (lbls, mats_o2, mats_o0): - for i, j in combinations_with_replacement (range (nroots), 2): - p, q = loff0[i], loff1[i] - r, s = loff0[j], loff1[j] - mtfp = lib.fp (mt[p:q,r:s]) - mrfp = lib.fp (mr[p:q,r:s]) - with self.subTest((i,j), opt=2, matrix=lbl): - self.assertAlmostEqual (mtfp, mrfp, 9) - #def test_rdm12s (self): - # d12_o0 = op_o0.roots_make_rdm12s (las, las.ci, nelec_frs, si)#, orbsym=orbsym, wfnsym=wfnsym) - # d12_o1 = op_o1.roots_make_rdm12s (las, las.ci, nelec_frs, si)#, orbsym=orbsym, wfnsym=wfnsym) - # d12_o2 = op_o2.roots_make_rdm12s (las, las.ci, nelec_frs, si)#, orbsym=orbsym, wfnsym=wfnsym) - # for r in range (2): - # for i in range (nroots): - # with self.subTest (rank=r+1, root=i, opt=1): - # self.assertAlmostEqual (lib.fp (d12_o0[r][i]), - # lib.fp (d12_o1[r][i]), 9) - # with self.subTest (rank=r+1, root=i, opt=2): - # self.assertAlmostEqual (lib.fp (d12_o0[r][i]), - # lib.fp (d12_o2[r][i]), 9) - # with self.subTest ('single matrix constructor', opt=0, rank=r+1, root=i): - # d12_o0_test = root_make_rdm12s (las, las.ci, si, state=i, soc=False, - # break_symmetry=False, opt=0)[r] - # self.assertAlmostEqual (lib.fp (d12_o0_test), lib.fp (d12_o0[r][i]), 9) - # with self.subTest ('single matrix constructor', opt=1, rank=r+1, root=i): - # d12_o1_test = root_make_rdm12s (las, las.ci, si, state=i, soc=False, - # break_symmetry=False, opt=1)[r] - # self.assertAlmostEqual (lib.fp (d12_o1_test), lib.fp (d12_o0[r][i]), 9) - # # I don't need to test o2 with this separately because the index-down is - # # entirely within the o1 parent function - # #with self.subTest ('single matrix constructor', opt=2, rank=r+1, root=i): - # # d12_o2_test = root_make_rdm12s (las, las.ci, si, state=i, soc=False, - # # break_symmetry=False, opt=1)[r] - # # self.assertAlmostEqual (lib.fp (d12_o2_test), lib.fp (d12_o0[r][i]), 9) + def test_rdm12s (self): + d12_o0 = op_o0.roots_make_rdm12s (las, las.ci, nelec_frs, si)#, orbsym=orbsym, wfnsym=wfnsym) + d12_o1 = op_o1.roots_make_rdm12s (las, las.ci, nelec_frs, si)#, orbsym=orbsym, wfnsym=wfnsym) + for r in range (2): + for i in range (nroots): + with self.subTest (rank=r+1, root=i, opt=1): + self.assertAlmostEqual (lib.fp (d12_o0[r][i]), + lib.fp (d12_o1[r][i]), 9) + with self.subTest ('single matrix constructor', opt=0, rank=r+1, root=i): + d12_o0_test = root_make_rdm12s (las, las.ci, si, state=i, soc=False, + break_symmetry=False, opt=0)[r] + self.assertAlmostEqual (lib.fp (d12_o0_test), lib.fp (d12_o0[r][i]), 9) + with self.subTest ('single matrix constructor', opt=1, rank=r+1, root=i): + d12_o1_test = root_make_rdm12s (las, las.ci, si, state=i, soc=False, + break_symmetry=False, opt=1)[r] + self.assertAlmostEqual (lib.fp (d12_o1_test), lib.fp (d12_o0[r][i]), 9) - #def test_lassis (self): - # las0 = las.get_single_state_las (state=0) - # for ifrag in range (len (las0.ci)): - # las0.ci[ifrag][0] = las0.ci[ifrag][0][0] - # lsi = LASSIS (las0) - # lsi.prepare_states_() - # self.assertTrue (lsi.converged) + def test_lassis_fast (self): + las0 = las.get_single_state_las (state=0) + for ifrag in range (len (las0.ci)): + las0.ci[ifrag][0] = las0.ci[ifrag][0][0] + lsi = LASSIS (las0) + lsi.prepare_states_() + self.assertTrue (lsi.converged) + case_lassis_fbf_2_model_state (self, lsi) - #def test_lassis_1111 (self): - # xyz='''H 0 0 0 - # H 3 0 0 - # H 6 0 0 - # H 9 0 0''' - # mol1 = gto.M (atom=xyz, basis='sto3g', symmetry=False, verbose=0, output='/dev/null') - # mf1 = scf.RHF (mol1).run () + def test_lassis_1111 (self): + xyz='''H 0 0 0 + H 3 0 0 + H 6 0 0 + H 9 0 0''' + mol1 = gto.M (atom=xyz, basis='sto3g', symmetry=False, verbose=0, output='/dev/null') + mf1 = scf.RHF (mol1).run () - # las1 = LASSCF (mf1, (1,1,1,1), ((0,1),(1,0),(0,1),(1,0))) - # mo_coeff = las1.localize_init_guess ([[0,],[1,],[2,],[3,]]) - # las1.lasci_(mo_coeff) - # lsi = LASSIS (las1).run () - # self.assertTrue (lsi.converged) - # self.assertAlmostEqual (lsi.e_roots[0], -1.867291372401379, 6) + las1 = LASSCF (mf1, (1,1,1,1), ((0,1),(1,0),(0,1),(1,0))) + mo_coeff = las1.localize_init_guess ([[0,],[1,],[2,],[3,]]) + las1.lasci_(mo_coeff) + lsi = LASSIS (las1).run () + self.assertTrue (lsi.converged) + self.assertAlmostEqual (lsi.e_roots[0], -1.867291372401379, 6) + case_lassis_fbf_2_model_state (self, lsi) + case_lassis_fbfdm (self, lsi) - #def test_lassis_slow (self): - # las0 = las.get_single_state_las (state=0) - # for ifrag in range (len (las0.ci)): - # las0.ci[ifrag][0] = las0.ci[ifrag][0][0] - # lsi = LASSIS (las0).run () - # self.assertTrue (lsi.converged) - # self.assertAlmostEqual (lsi.e_roots[0], -304.5372586630968, 3) + def test_lassis_slow (self): + las0 = las.get_single_state_las (state=0) + for ifrag in range (len (las0.ci)): + las0.ci[ifrag][0] = las0.ci[ifrag][0][0] + lsi = LASSIS (las0).run () + self.assertTrue (lsi.converged) + self.assertAlmostEqual (lsi.e_roots[0], -304.5372586630968, 3) + case_lassis_fbf_2_model_state (self, lsi) + #case_lassis_fbfdm (self, lsi) - #def test_fdm1 (self): - # make_fdm1 = get_fdm1_maker (las, las.ci, nelec_frs, si) - # for iroot in range (nroots): - # for ifrag in range (4): - # with self.subTest (iroot=iroot, ifrag=ifrag): - # fdm1 = make_fdm1 (iroot, ifrag) - # sdm1 = make_sdm1 (las, iroot, ifrag, si=si) - # self.assertAlmostEqual (lib.fp (fdm1), lib.fp (sdm1), 7) + def test_fdm1 (self): + make_fdm1 = get_fdm1_maker (las, las.ci, nelec_frs, si) + for iroot in range (nroots): + for ifrag in range (4): + with self.subTest (iroot=iroot, ifrag=ifrag): + fdm1 = make_fdm1 (iroot, ifrag) + sdm1 = make_sdm1 (las, iroot, ifrag, si=si) + self.assertAlmostEqual (lib.fp (fdm1), lib.fp (sdm1), 7) + + def test_contract_hlas_ci (self): + h0, h1, h2 = ham_2q (las, las.mo_coeff) + case_contract_hlas_ci (self, las, h0, h1, h2, las.ci, nelec_frs) if __name__ == "__main__": diff --git a/my_pyscf/lassi/excitations.py b/my_pyscf/lassi/excitations.py index ec88b843..46353ce1 100644 --- a/my_pyscf/lassi/excitations.py +++ b/my_pyscf/lassi/excitations.py @@ -433,7 +433,7 @@ def fixedpoint (self, h1, h2, norb_f, nelec_f, ecore=0, ci0=None, orbsym=None, e, si_p, si_q, ci0 = self._eig (h0, h1, h2, ci1, nroots=nroots)[:4] hpq_xq = self.get_hpq_xq (h1, h2, ci0, si_q) hpp_xp = self.get_hpp_xp (h1, h2, ci0, si_p, norb_f, nelec_f, ecore=h0, nroots=nroots) - grad = self._get_grad (si_p, hpq_xq, hpp_xp) + grad = self._get_grad (ci0, si_p, hpq_xq, hpp_xp, nroots=nroots) grad_max = np.amax (np.abs (grad)) log.info ('Cycle %d: max grad = %e ; e = %e, |delta| = %e', it, grad_max, e, e - e_last) @@ -447,7 +447,7 @@ def fixedpoint (self, h1, h2, norb_f, nelec_f, ecore=0, ci0=None, orbsym=None, if not converged: ci1 = self.get_init_guess (ci1, norb_f, nelec_f, h1, h2, nroots=nroots) # Issue #86: see above, same problem - self._debug_csfs (log, ci0, ci1, norb_f, nelec_f, grad) + self._debug_csfs (log, ci0, ci1, norb_f, nelec_f, grad, nroots=nroots) return converged, e, ci1 def _eig (self, h0, h1, h2, ci0, ovlp_thresh=1e-3, nroots=1): @@ -478,7 +478,7 @@ def get_hpq_xq (self, h1, h2, ci0, si_q): hci_f_pab = [] for ci, hci_pabq in zip (ci0, hci_f_pabq): hci_pab = np.dot (hci_pabq, si_q) - hci_pr = np.tensordot (hci_pab, ci, axes=((1,2),(1,2))) + hci_pr = np.tensordot (hci_pab, ci.conj (), axes=((1,2),(1,2))) hci_pab -= np.tensordot (hci_pr, ci, axes=1) hci_f_pab.append (hci_pab) return hci_f_pab @@ -527,19 +527,27 @@ def get_hpp_xp (self, h1, h2, ci0, si_p, norb_f, nelec_f, ecore=0, nroots=1, **k hc[iroot] += sj * contract_1e_nosym (veff, jcol, no, nelec) c, hc = np.asarray (c), np.asarray (hc) chc = np.dot (np.asarray (c).reshape (nroots,-1).conj (), - np.asarray (hc).reshape (nroots,-1).T) + np.asarray (hc).reshape (nroots,-1).T).T hc = hc - np.tensordot (chc, c, axes=1) hci_f_pab.append (hc) return hci_f_pab - def _get_grad (self, si_p, hpq_xq, hpp_xp): + def _get_grad (self, ci0, si_p, hpq_xq, hpp_xp, nroots=None): # Compute the gradient of the target interacting energy grad = [] - for solver, hc1, hc2 in zip (self.fcisolvers, hpq_xq, hpp_xp): + for solver, c, hc1, hc2 in zip (self.fcisolvers, ci0, hpq_xq, hpp_xp): + if nroots is None: nroots = solver.nroots hc = si_p[:,None,None] * (hc1 + hc2) if isinstance (solver, CSFFCISolver): + c = solver.transformer.vec_det2csf (c, normalize=True) hc = solver.transformer.vec_det2csf (hc, normalize=False) + chc = np.dot (c.conj (), hc.T) + hc = hc - np.dot (chc.T, c) grad.append (hc.flat) + if nroots>1: + # TODO: figure out how this gradient should actually work + chc -= chc.T + grad.append (np.zeros_like (chc[np.tril_indices (nroots, k=-1)])) return np.concatenate (grad) def _1shot (self, h0, h1, h2, ci0, hpq_xq, hpp_xp, nroots=1, ovlp_thresh=1e-3): diff --git a/my_pyscf/mcscf/productstate.py b/my_pyscf/mcscf/productstate.py index f18ff7d1..3bf8d441 100644 --- a/my_pyscf/mcscf/productstate.py +++ b/my_pyscf/mcscf/productstate.py @@ -146,15 +146,16 @@ def _check_init_guess (self, ci0, norb_f, nelec_f, nroots=None): snroots, solver.transformer.ncsf, ix, solver.nelec, solver.smult)) return ci1 - def _debug_csfs (self, log, ci0, ci1, norb_f, nelec_f, grad): + def _debug_csfs (self, log, ci0, ci1, norb_f, nelec_f, grad, nroots=None): if not all ([isinstance (s, CSFFCISolver) for s in self.fcisolvers]): return if log.verbose < lib.logger.INFO: return transformers = [s.transformer for s in self.fcisolvers] grad_f = [] for s,t in zip (self.fcisolvers, transformers): - grad_f.append (grad[:t.ncsf*s.nroots].reshape (s.nroots, t.ncsf)) - offs = (t.ncsf*s.nroots) + (s.nroots*(s.nroots-1)//2) + snroots = nroots if nroots is not None else s.nroots + grad_f.append (grad[:t.ncsf*snroots].reshape (snroots, t.ncsf)) + offs = (t.ncsf*snroots) + (snroots*(snroots-1)//2) grad = grad[offs:] assert (len (grad) == 0) log.info ('Debugging CI and gradient vectors...')