Skip to content

Commit

Permalink
lassi o1 soc debugged (so far)
Browse files Browse the repository at this point in the history
  • Loading branch information
@MatthewRHermes
MatthewRHermes committed Sep 10, 2024
1 parent c15699f commit 2798b10
Show file tree
Hide file tree
Showing 3 changed files with 57 additions and 32 deletions.
14 changes: 5 additions & 9 deletions debug/lassi/debug_soc.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -140,15 +140,17 @@ def case_soc_rdm12s_slow (self, opt=0):
self.assertAlmostEqual (lib.fp (rdm1s_test), lib.fp (rdm1s_ref), 10)
with self.subTest (sanity='dm2s'):
self.assertAlmostEqual (lib.fp (rdm2s_test), lib.fp (rdm2s_ref), 10)
# Stationary test has the issue that the 2nd and third states are degenerate.
# Stationary test has the issue of two doubly-degenerate manifolds: 1,2 and 4,5.
# Therefore their RDMs actually vary randomly. Average the second and third RDMs
# together to deal with this.
rdm1s_test[1:3] = rdm1s_test[1:3].sum (0) / 2
rdm2s_test[1:3] = rdm2s_test[1:3].sum (0) / 2
rdm1s_test[4:6] = rdm1s_test[4:6].sum (0) / 2
rdm2s_test[4:6] = rdm2s_test[4:6].sum (0) / 2
with self.subTest ('2-electron'):
self.assertAlmostEqual (linalg.norm (rdm2s_test), 13.767204017387787)
self.assertAlmostEqual (linalg.norm (rdm2s_test), 13.584509751113796)
with self.subTest ('1-electron'):
self.assertAlmostEqual (linalg.norm (rdm1s_test), 5.299117089455725)
self.assertAlmostEqual (linalg.norm (rdm1s_test), 5.298727485035966)
with lib.light_speed (5):
e0, h1, h2 = ham_2q (las2, las2.mo_coeff, soc=True)
rdm2_test = rdm2s_test.sum ((1,4))
Expand Down Expand Up @@ -246,8 +248,6 @@ def test_hso (hso_test, tag='kernel'):
test_hso (h0eff + h1eff + h2eff, 'make_stdm12s')
rdm1s, rdm2s = roots_make_rdm12s (las, las.ci, si, soc=True, break_symmetry=True,
opt=opt)
np.save ('rdm1s_o{}_1frag.npy'.format (opt), rdm1s)
np.save ('rdm2s_o{}_1frag.npy'.format (opt), rdm2s)
rdm2 = rdm2s.sum ((1,4))
e1eff = lib.einsum ('pq,iqp->i', h1, rdm1s)
e2eff = lib.einsum ('pqrs,ipqrs->i', h2, rdm2) * .5
Expand Down Expand Up @@ -289,11 +289,7 @@ def test_soc_rdm12s_slow_o0 (self):
def test_soc_rdm12s_slow_o1 (self):
#case_soc_rdm12s_slow (self, opt=1)
d_test = roots_make_rdm12s (las2, las2.ci, lsi2.si, opt=1)
np.save ('rdm1s_o1_2frag.npy', d_test[0])
np.save ('rdm2s_o1_2frag.npy', d_test[1])
d_ref = roots_make_rdm12s (las2, las2.ci, lsi2.si, opt=0)
np.save ('rdm1s_o0_2frag.npy', d_ref[0])
np.save ('rdm2s_o0_2frag.npy', d_ref[1])
for i in range (len (d_test[0])):
for r in range (2):
with self.subTest (state=i, rank=r+1):
Expand Down
2 changes: 1 addition & 1 deletion my_pyscf/lassi/op_o1/hams2ovlp.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -276,7 +276,7 @@ def _crunch_1c1d_(self, bra, ket, i, j, k, s1):
h_ = self.get_ham_2q (j,k,k,i).transpose (1,2,0,3) # BEWARE CONJ
h_ = np.tensordot (p_i, h_, axes=((-1),(-1)))
h_ = np.tensordot (h_j, h_, axes=((-1),(-1)))
ham -= np.tensordot (d1_k, h_, axes=((-1,-2),(-2,-1)))
ham -= np.tensordot (d1_k, h_, axes=((-2,-1),(-2,-1)))
ham *= fac
s2 = None
dt, dw = logger.process_clock () - t0, logger.perf_counter () - w0
Expand Down
73 changes: 51 additions & 22 deletions tests/lassi/test_soc.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -14,7 +14,11 @@
import itertools

def setUpModule():
global mol1, mf1, mol2, mf2, las2, lsi2
global mol1, mf1, mol2, mf2, las2, lsi2, oldvars
oldvars = {}
from mrh.my_pyscf.lassi.op_o1 import frag
oldvars['SCREEN_THRESH'] = frag.SCREEN_THRESH
frag.SCREEN_THRESH = 1e-32
mol1 = gto.M (atom="""
O 0.000000 0.000000 0.000000
H 0.758602 0.000000 0.504284
Expand All @@ -38,9 +42,10 @@ def setUpModule():
las2.mo_coeff = las2.localize_init_guess ((list (range (3)), list (range (9,12))), mf2.mo_coeff)
# NOTE: for 2-fragment and above, you will ALWAYS need to remember to do the line above.
# If you skip it, you can expect your orbitals not to converge.
las2.state_average_(weights=[1,0,0,0,0],
spins=[[0,0],[2,0],[-2,0],[0,2],[0,-2]],
smults=[[1,1],[3,1],[3,1],[1,3],[1,3]])
las2.state_average_(weights=[1,0,0,0,0,0,0],
spins=[[0,0],[2,0],[-2,0],[0,2],[0,-2],[0,0],[0,0]],
smults=[[1,1],[3,1],[3,1],[1,3],[1,3],[3,1],[1,3]])

# NOTE: Be careful about state selection. You have to select states that can actually be coupled
# by a 1-body SOC operator. For instance, spins=[0,0] and spins=[2,2] would need at least a 2-body
# operator to couple.
Expand All @@ -52,18 +57,19 @@ def setUpModule():
lsi2.kernel (opt=0)

def tearDownModule():
global mol1, mf1, mol2, mf2, las2, lsi2
global mol1, mf1, mol2, mf2, las2, lsi2, oldvars
from mrh.my_pyscf.lassi.op_o1 import frag
for key, val in oldvars.items (): setattr (frag, key, val)
mol1.stdout.close()
mol2.stdout.close()
del mol1, mf1, mol2, mf2, las2, lsi2
del mol1, mf1, mol2, mf2, las2, lsi2, oldvars

def case_soc_stdm12s_slow (self, opt=0):
stdm1s_test, stdm2s_test = make_stdm12s (las2, soc=True, opt=opt)
np.save ('stdm1s.{}.npy'.format (opt), stdm1s_test)
with self.subTest ('2-electron'):
self.assertAlmostEqual (linalg.norm (stdm2s_test), 12.835690990485933, 6)
self.assertAlmostEqual (linalg.norm (stdm2s_test), 16.901692823561433, 6)
with self.subTest ('1-electron'):
self.assertAlmostEqual (linalg.norm (stdm1s_test), 5.719302720036956, 6)
self.assertAlmostEqual (linalg.norm (stdm1s_test), 7.075259874940101, 6)
dm1s_test = lib.einsum ('ipqi->ipq', stdm1s_test)
with self.subTest (oneelectron_sanity='diag'):
# LAS states are spin-pure: there should be nothing in the spin-breaking sector
Expand Down Expand Up @@ -100,7 +106,7 @@ def case_soc_stdm12s_slow (self, opt=0):
self.assertAlmostEqual (np.amax (np.abs (d1[1,:,0,:])), 0, 16)
d1=d1[0,:,1,:] # [alpha,beta] -> beta' alpha
with self.subTest (oneelectron_sanity='nonzero S.O.C.', ket=i):
self.assertAlmostEqual (linalg.norm (d1), 1.1539612706310192, 8)
self.assertAlmostEqual (linalg.norm (d1), 1.1539612627187337, 8)
# lassi_dms.make_trans and total electron count
ncas = las2.ncas
nelec_fr = [[_unpack_nelec (fcibox._get_nelec (solver, nelecas))
Expand Down Expand Up @@ -134,15 +140,17 @@ def case_soc_rdm12s_slow (self, opt=0):
self.assertAlmostEqual (lib.fp (rdm1s_test), lib.fp (rdm1s_ref), 10)
with self.subTest (sanity='dm2s'):
self.assertAlmostEqual (lib.fp (rdm2s_test), lib.fp (rdm2s_ref), 10)
# Stationary test has the issue that the 2nd and third states are degenerate.
# Stationary test has the issue of two doubly-degenerate manifolds: 1,2 and 4,5.
# Therefore their RDMs actually vary randomly. Average the second and third RDMs
# together to deal with this.
rdm1s_test[1:3] = rdm1s_test[1:3].sum (0) / 2
rdm2s_test[1:3] = rdm2s_test[1:3].sum (0) / 2
rdm1s_test[4:6] = rdm1s_test[4:6].sum (0) / 2
rdm2s_test[4:6] = rdm2s_test[4:6].sum (0) / 2
with self.subTest ('2-electron'):
self.assertAlmostEqual (linalg.norm (rdm2s_test), 11.39986880775559)
self.assertAlmostEqual (linalg.norm (rdm2s_test), 13.584509751113796)
with self.subTest ('1-electron'):
self.assertAlmostEqual (linalg.norm (rdm1s_test), 4.47832557674054)
self.assertAlmostEqual (linalg.norm (rdm1s_test), 5.298727485035966)
with lib.light_speed (5):
e0, h1, h2 = ham_2q (las2, las2.mo_coeff, soc=True)
rdm2_test = rdm2s_test.sum ((1,4))
Expand Down Expand Up @@ -197,10 +205,12 @@ def test_soc_1frag (self):
with self.subTest (deltaE='SF'):
esf_test = las.e_states - e0
self.assertAlmostEqual (lib.fp (esf_test), lib.fp (esf_ref), 6)
ham = [None, None]
for opt in (0,1):
with lib.light_speed (10):
e_roots, si = las.lassi (opt=opt, soc=True, break_symmetry=True)
h0, h1, h2 = ham_2q (las, las.mo_coeff, soc=True)
ham[opt] = (si * e_roots[None,:]) @ si.conj ().T
eso_test = e_roots - e0
with self.subTest (opt=opt, deltaE='SO'):
self.assertAlmostEqual (lib.fp (eso_test), lib.fp (eso_ref), 6)
Expand Down Expand Up @@ -236,37 +246,56 @@ def test_hso (hso_test, tag='kernel'):
h1eff = lib.einsum ('pq,iqpj->ij', h1, stdm1s)
h2eff = lib.einsum ('pqrs,ipqrsj->ij', h2, stdm2) * .5
test_hso (h0eff + h1eff + h2eff, 'make_stdm12s')
if opt==1: continue # TODO: fix bug and remove this line
rdm1s, rdm2s = roots_make_rdm12s (las, las.ci, si, soc=True, break_symmetry=True,
opt=opt)
rdm2 = rdm2s.sum ((1,4))
e1eff = lib.einsum ('pq,iqp->i', h1, rdm1s)
e2eff = lib.einsum ('pqrs,ipqrs->i', h2, rdm2) * .5
test_hso ((si * (e0eff+e1eff+e2eff)[None,:]) @ si.conj ().T, 'roots_make_rdm12s')
with self.subTest ('o0-o1 ham agreement'):
self.assertAlmostEqual (lib.fp (ham[0]), lib.fp (ham[1]), 8)

def test_soc_2frag (self):
## stationary test for >1 frag calc
with self.subTest (deltaE='SF'):
self.assertAlmostEqual (lib.fp (lsi2._las.e_states), 154.09610921621356, 8)
self.assertAlmostEqual (lib.fp (lsi2._las.e_states), -214.8686632658775, 8)
with self.subTest (opt=0, deltaE='SO'):
self.assertAlmostEqual (lib.fp (lsi2.e_roots), 154.09559506105586, 8)
self.assertAlmostEqual (lib.fp (lsi2.e_roots), -214.8684319949548, 8)
lsi = lassi.LASSI (lsi2._las, soc=True, break_symmetry=True, opt=1)
with lib.light_speed (5): lsi.kernel (opt=1)
with self.subTest (opt=1, deltaE='SO'):
lsi = lassi.LASSI (lsi2._las, soc=True, break_symmetry=True, opt=1)
with lib.light_speed (5): lsi.kernel (opt=1)
self.assertAlmostEqual (lib.fp (lsi.e_roots), 154.09559506105586, 8)
self.assertAlmostEqual (lib.fp (lsi.e_roots), -214.8684319949548, 8)
with self.subTest ('hamiltonian', opt=1):
ham_o0 = (lsi2.si * lsi2.e_roots[None,:]) @ lsi2.si.conj ().T
ham_o1 = (lsi.si * lsi.e_roots[None,:]) @ lsi.si.conj ().T
self.assertAlmostEqual (lib.fp (ham_o1), lib.fp (ham_o0), 8)

def test_soc_stdm12s_slow_o0 (self):
case_soc_stdm12s_slow (self, opt=0)

def test_soc_stdm12s_slow_o1 (self):
case_soc_stdm12s_slow (self, opt=1)
#case_soc_stdm12s_slow (self, opt=1)
d_test = make_stdm12s (las2, soc=True, opt=1)
d_ref = make_stdm12s (las2, soc=True, opt=0)
for i,j in itertools.product (range (len (d_test[0])), repeat=2):
for r in range (2):
with self.subTest (rank=r+1, element=(i,j)):
self.assertAlmostEqual (lib.fp (d_test[r][i,...,j]),
lib.fp (d_ref[r][i,...,j]), 8)

def test_soc_rdm12s_slow_o0 (self):
case_soc_rdm12s_slow (self, opt=0)

#def test_soc_rdm12s_slow_o1 (self):
# case_soc_rdm12s_slow (self, opt=1)
def test_soc_rdm12s_slow_o1 (self):
#case_soc_rdm12s_slow (self, opt=1)
d_test = roots_make_rdm12s (las2, las2.ci, lsi2.si, opt=1)
d_ref = roots_make_rdm12s (las2, las2.ci, lsi2.si, opt=0)
for i in range (len (d_test[0])):
for r in range (2):
with self.subTest (state=i, rank=r+1):
self.assertAlmostEqual (lib.fp (d_test[r][i]), lib.fp (d_ref[r][i]), 8)

if __name__ == "__main__":
print("Full Tests for SOC")
unittest.main()

0 comments on commit 2798b10

Please sign in to comment.