Skip to content

Commit

Permalink
Fix grad
Browse files Browse the repository at this point in the history
  • Loading branch information
@MatthewRHermes
MatthewRHermes committed Nov 14, 2024
1 parent f699832 commit d01f76f
Show file tree
Hide file tree
Showing 3 changed files with 95 additions and 102 deletions.
170 changes: 77 additions & 93 deletions debug/lassi/debug_4frag.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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]],
Expand Down Expand Up @@ -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)
Expand All @@ -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)]
Expand Down Expand Up @@ -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__":
Expand Down
20 changes: 14 additions & 6 deletions my_pyscf/lassi/excitations.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand All @@ -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):
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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):
Expand Down
7 changes: 4 additions & 3 deletions my_pyscf/mcscf/productstate.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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...')
Expand Down

0 comments on commit d01f76f

Please sign in to comment.