From 784ab7e4cb2aa4bde3df9d9a9d2d3ecfc892bef3 Mon Sep 17 00:00:00 2001 From: lweisburn Date: Wed, 14 May 2025 17:58:03 -0400 Subject: [PATCH 01/31] remove UBE df option --- src/quemb/molbe/ube.py | 27 ++++++++++----------------- 1 file changed, 10 insertions(+), 17 deletions(-) diff --git a/src/quemb/molbe/ube.py b/src/quemb/molbe/ube.py index d5f8251da3..80c93262e8 100644 --- a/src/quemb/molbe/ube.py +++ b/src/quemb/molbe/ube.py @@ -228,25 +228,18 @@ def initialize(self, eri_, compute_hf): thr_bath=self.thr_bath, ) - if eri_ is None and self.mf.with_df is not None: - # NOT IMPLEMENTED: should not be called, as no unrestricted DF tested - # for density-fitted integrals; if mf is provided, pyscf.ao2mo uses DF - # object in an outcore fashion - eri_a = ao2mo.kernel(self.mf.mol, fobj_a.TA, compact=True) - eri_b = ao2mo.kernel(self.mf.mol, fobj_b.TA, compact=True) - else: - eri_a = ao2mo.incore.full( - eri_, fobj_a.TA, compact=True - ) # otherwise, do an incore ao2mo - eri_b = ao2mo.incore.full(eri_, fobj_b.TA, compact=True) + assert eri_ is not None, "eri_ is None: set incore_anyway for UHF" - Csd_A = fobj_a.TA # may have to add in nibath here - Csd_B = fobj_b.TA + eri_a = ao2mo.incore.full(eri_, fobj_a.TA, compact=True) + eri_b = ao2mo.incore.full(eri_, fobj_b.TA, compact=True) - # cross-spin ERI term - eri_ab = ao2mo.incore.general( - eri_, (Csd_A, Csd_A, Csd_B, Csd_B), compact=True - ) + Csd_A = fobj_a.TA # may have to add in nibath here + Csd_B = fobj_b.TA + + # cross-spin ERI term + eri_ab = ao2mo.incore.general( + eri_, (Csd_A, Csd_A, Csd_B, Csd_B), compact=True + ) file_eri.create_dataset(fobj_a.dname[0], data=eri_a) file_eri.create_dataset(fobj_a.dname[1], data=eri_b) From d85a8a8efa1ec8dad918b25c19a6798fc0c6cd8e Mon Sep 17 00:00:00 2001 From: lweisburn Date: Wed, 14 May 2025 17:58:36 -0400 Subject: [PATCH 02/31] remove UBE df option --- src/quemb/molbe/ube.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/quemb/molbe/ube.py b/src/quemb/molbe/ube.py index 80c93262e8..5ea5d0d06f 100644 --- a/src/quemb/molbe/ube.py +++ b/src/quemb/molbe/ube.py @@ -228,6 +228,7 @@ def initialize(self, eri_, compute_hf): thr_bath=self.thr_bath, ) + assert fobj_a.TA is not None and fobj_b.TA is not None assert eri_ is not None, "eri_ is None: set incore_anyway for UHF" eri_a = ao2mo.incore.full(eri_, fobj_a.TA, compact=True) @@ -248,7 +249,6 @@ def initialize(self, eri_, compute_hf): # sab = self.C_a @ self.S @ self.C_b _ = fobj_a.get_nsocc(self.S, self.C_a, self.Nocc[0], ncore=self.ncore) - assert fobj_a.TA is not None fobj_a.h1 = multi_dot((fobj_a.TA.T, self.hcore, fobj_a.TA)) eri_a = ao2mo.restore(8, eri_a, fobj_a.nao) @@ -273,7 +273,6 @@ def initialize(self, eri_, compute_hf): _ = fobj_b.get_nsocc(self.S, self.C_b, self.Nocc[1], ncore=self.ncore) - assert fobj_b.TA is not None fobj_b.h1 = multi_dot((fobj_b.TA.T, self.hcore, fobj_b.TA)) eri_b = ao2mo.restore(8, eri_b, fobj_b.nao) fobj_b.cons_fock(self.hf_veff[1], self.S, self.hf_dm[1] * 2.0, eri_=eri_b) From c84bb888ee4e1199b3dcc169c3cfd0eebe5dee29 Mon Sep 17 00:00:00 2001 From: lweisburn Date: Fri, 16 May 2025 17:35:56 -0400 Subject: [PATCH 03/31] adding equal_bath option --- src/quemb/molbe/pfrag.py | 21 ++---- src/quemb/molbe/ube.py | 69 ++++++++++++------- .../shared/external/unrestricted_utils.py | 5 +- 3 files changed, 55 insertions(+), 40 deletions(-) diff --git a/src/quemb/molbe/pfrag.py b/src/quemb/molbe/pfrag.py index 09d826b9f0..a556796e0d 100644 --- a/src/quemb/molbe/pfrag.py +++ b/src/quemb/molbe/pfrag.py @@ -6,7 +6,6 @@ import numpy as np import scipy.linalg from numpy import ( - argsort, array, diag_indices, einsum, @@ -469,28 +468,22 @@ def schmidt_decomposition( # Identify significant environment orbitals based on eigenvalue threshold Bidx = [] - + for i in range(len(Eval)): + if thr_bath < np.abs(Eval[i]) < 1.0 - thr_bath: + Bidx.append(i) # Set the number of orbitals to be taken from the environment orbitals # Based on an eigenvalue threshold ordering if norb is not None: - n_frag_ind = len(Frag_sites1) - n_bath_ind = norb - n_frag_ind - ind_sort = argsort(np.abs(Eval)) - first_el = [x for x in ind_sort if x < 1.0 - thr_bath][-1 * n_bath_ind] - for i in range(len(Eval)): - if np.abs(Eval[i]) >= first_el: - Bidx.append(i) - else: - for i in range(len(Eval)): - if thr_bath < np.abs(Eval[i]) < 1.0 - thr_bath: - Bidx.append(i) + # add extra correlated orbital from environment + # this will likely have Eval = 1 + more_corr = Bidx[-1] + 1 + Bidx.append(more_corr) # Initialize the transformation matrix (TA) TA = zeros([Tot_sites, len(AO_in_frag) + len(Bidx)]) TA[AO_in_frag, : len(AO_in_frag)] = eye(len(AO_in_frag)) # Fragment part TA[Env_sites1, len(AO_in_frag) :] = Evec[:, Bidx] # Environment part - # return TA, norbs_frag, norbs_bath return TA, Frag_sites1.shape[0], len(Bidx) diff --git a/src/quemb/molbe/ube.py b/src/quemb/molbe/ube.py index 5ea5d0d06f..93f6699ce4 100644 --- a/src/quemb/molbe/ube.py +++ b/src/quemb/molbe/ube.py @@ -45,6 +45,7 @@ def __init__( pop_method: str | None = None, compute_hf: bool = True, thr_bath: float = 1.0e-10, + equal_bath: bool = True, ) -> None: """Initialize Unrestricted BE Object (ube🍠) @@ -72,11 +73,17 @@ def __init__( pop_method : Method for calculating orbital population, by default 'meta-lowdin' See pyscf.lo for more details and options - thr_bath : float, + thr_bath : Threshold for bath orbitals in Schmidt decomposition + equal_bath : + Whether to use a bath with the same number of alpha and beta orbitals. + Using equal_bath = False will require custom compiled functions in + PySCF to perform integral transformations. Default is True """ + self.unrestricted = True self.thr_bath = thr_bath + self.equal_bath = equal_bath self.fobj = fobj @@ -200,33 +207,45 @@ def initialize(self, eri_, compute_hf): if self.frozen_core: fobj_a.core_veff = self.core_veff[0] fobj_b.core_veff = self.core_veff[1] - fobj_a.sd( - self.W[0], - self.lmo_coeff_a, - self.Nocc[0], - thr_bath=self.thr_bath, - ) - fobj_b.sd( - self.W[1], - self.lmo_coeff_b, - self.Nocc[1], - thr_bath=self.thr_bath, - ) else: fobj_a.core_veff = None fobj_b.core_veff = None - fobj_a.sd( - self.W, - self.lmo_coeff_a, - self.Nocc[0], - thr_bath=self.thr_bath, - ) - fobj_b.sd( - self.W, - self.lmo_coeff_b, - self.Nocc[1], - thr_bath=self.thr_bath, - ) + + fobj_a.sd( + self.W[0] if self.frozen_core else self.W, + self.lmo_coeff_a, + self.Nocc[0], + thr_bath=self.thr_bath, + ) + fobj_b.sd( + self.W[1] if self.frozen_core else self.W, + self.lmo_coeff_b, + self.Nocc[1], + thr_bath=self.thr_bath, + ) + + if self.equal_bath: + # Enforce the same number of alpha and beta orbitals + # by augmenting the bath + tot_alpha = fobj_a.n_f + fobj_a.n_b + tot_beta = fobj_b.n_f + fobj_b.n_b + + if tot_alpha > tot_beta: + fobj_b.sd( + self.W[1] if self.frozen_core else self.W, + self.lmo_coeff_b, + self.Nocc[1], + thr_bath=self.thr_bath, + norb=fobj_a.n_b, + ) + elif tot_beta > tot_alpha: + fobj_a.sd( + self.W[0] if self.frozen_core else self.W, + self.lmo_coeff_a, + self.Nocc[0], + thr_bath=self.thr_bath, + norb=fobj_b.n_b, + ) assert fobj_a.TA is not None and fobj_b.TA is not None assert eri_ is not None, "eri_ is None: set incore_anyway for UHF" diff --git a/src/quemb/shared/external/unrestricted_utils.py b/src/quemb/shared/external/unrestricted_utils.py index 42033bd4dc..ba991c6fa9 100644 --- a/src/quemb/shared/external/unrestricted_utils.py +++ b/src/quemb/shared/external/unrestricted_utils.py @@ -65,7 +65,10 @@ def uccsd_restore_eris(symm, fobj_a, fobj_b): Vs = [None] * 3 Vs[0] = ao2mo.restore(symm[0], fVs[Vsname[0]], nf[0]) Vs[1] = ao2mo.restore(symm[1], fVs[Vsname[1]], nf[1]) - Vs[2] = restore_eri_gen(symm[2], fVs[Vsname[2]][()], nf[0], nf[1]) + if nf[0] == nf[1]: + Vs[2] = ao2mo.restore(symm[2], fVs[Vsname[2]][()], nf[0]) + else: + Vs[2] = restore_eri_gen(symm[2], fVs[Vsname[2]][()], nf[0], nf[1]) return Vs From 2ad2cc82c5385523c9243a184174fd1a7dcc0931 Mon Sep 17 00:00:00 2001 From: lweisburn Date: Fri, 16 May 2025 18:18:51 -0400 Subject: [PATCH 04/31] remove unused doc options --- src/quemb/molbe/solver.py | 4 ---- 1 file changed, 4 deletions(-) diff --git a/src/quemb/molbe/solver.py b/src/quemb/molbe/solver.py index 55ec01a5b4..cc96312648 100644 --- a/src/quemb/molbe/solver.py +++ b/src/quemb/molbe/solver.py @@ -596,10 +596,6 @@ def be_func_u( Whether to evaluate the energy. Defaults to False. relax_density : bool, optional Whether to relax the density. Defaults to False. - return_vec : bool, optional - Whether to return the error vector. Defaults to False. - ebe_hf : float, optional - Hartree-Fock energy. Defaults to 0. use_cumulant : bool, optional Whether to use the cumulant-based energy expression. Defaults to True. frozen : bool, optional From 882fd14e58bf8e130121cff917452763879e73b4 Mon Sep 17 00:00:00 2001 From: lweisburn Date: Fri, 16 May 2025 18:33:53 -0400 Subject: [PATCH 05/31] update UBE test --- tests/ube-oneshot_test.py | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/tests/ube-oneshot_test.py b/tests/ube-oneshot_test.py index a86ab120bc..713c67d667 100644 --- a/tests/ube-oneshot_test.py +++ b/tests/ube-oneshot_test.py @@ -20,7 +20,7 @@ class TestOneShot_Unrestricted(unittest.TestCase): "This test is known to fail.", ) def test_hexene_anion_sto3g_frz_ben(self): - # Linear Equidistant (r=1Å) H8 Chain, STO-3G + # Hexene anion with frozen core, STO-3G mol = gto.M() mol.atom = os.path.join(os.path.dirname(__file__), "xyz/hexene.xyz") mol.basis = "sto-3g" @@ -42,7 +42,7 @@ def test_hexene_anion_sto3g_frz_ben(self): "This test is known to fail.", ) def test_hexene_cation_sto3g_frz_ben(self): - # Linear Equidistant (r=1Å) H8 Chain, cc-pVDZ + # Hexene cation with frozen core, STO-3G mol = gto.M() mol.atom = os.path.join(os.path.dirname(__file__), "xyz/hexene.xyz") mol.basis = "sto-3g" @@ -64,7 +64,7 @@ def test_hexene_cation_sto3g_frz_ben(self): "This test is known to fail.", ) def test_hexene_anion_sto3g_unfrz_ben(self): - # Octane, STO-3G + # Hexene anion without frozen core, STO-3G mol = gto.M() mol.atom = os.path.join(os.path.dirname(__file__), "xyz/hexene.xyz") mol.basis = "sto-3g" @@ -86,6 +86,7 @@ def test_hexene_anion_sto3g_unfrz_ben(self): "This test is known to fail.", ) def test_hexene_cation_sto3g_unfrz_ben(self): + # Hexene cation without frozen core, STO-3G mol = gto.M() mol.atom = os.path.join(os.path.dirname(__file__), "xyz/hexene.xyz") mol.basis = "sto-3g" @@ -107,7 +108,7 @@ def molecular_unrestricted_oneshot_test( ): mf = scf.UHF(mol) mf.kernel() - fobj = fragmentate(frag_type="autogen", n_BE=n_BE, mol=mol, frozen_core=frz) + fobj = fragmentate(frag_type="chemgen", n_BE=n_BE, mol=mol, frozen_core=frz) mybe = UBE(mf, fobj) mybe.oneshot(solver="UCCSD", nproc=1) self.assertAlmostEqual( From 45903b5e8255c7b406f7fd817e069a0bfbd1beaf Mon Sep 17 00:00:00 2001 From: Leah Weisburn Date: Tue, 17 Jun 2025 13:17:15 -0400 Subject: [PATCH 06/31] update (but not fix) tests and fix equal_bath --- example/molbe_oneshot_ube_qmmm.py | 7 +-- src/quemb/molbe/misc.py | 79 +++++++++++++++++++++---------- src/quemb/molbe/pfrag.py | 4 +- src/quemb/molbe/ube.py | 15 +++++- tests/ube-oneshot_test.py | 38 ++++++++++----- 5 files changed, 100 insertions(+), 43 deletions(-) diff --git a/example/molbe_oneshot_ube_qmmm.py b/example/molbe_oneshot_ube_qmmm.py index 5d61018d36..acbc05e5b5 100644 --- a/example/molbe_oneshot_ube_qmmm.py +++ b/example/molbe_oneshot_ube_qmmm.py @@ -33,7 +33,8 @@ unrestricted=True, # specify unrestricted calculation from_chk=False, # can save the UHF as PySCF checkpoint. # Set to true if running from converged UHF chk - checkfile=None, -) # if not None, will save UHF calculation to a checkfile. + checkfile=None, # if not None, will save UHF calculation to checkfile + opt="SOSCF", # SOSCF and DAMP options: note that damping settings hard-coded +) # if rerunning from chk (from_chk=True), name the checkfile here -# ecp = ecp) # can add ECP for heavy atoms as: {'Ru': 'def2-SVP'} +# ecp = ecp) # can add ECP for heavy atoms as: {'Ru': 'def2-SVP'} diff --git a/src/quemb/molbe/misc.py b/src/quemb/molbe/misc.py index b6be7838fd..f89aec7240 100644 --- a/src/quemb/molbe/misc.py +++ b/src/quemb/molbe/misc.py @@ -265,6 +265,7 @@ def be2puffin( checkfile=None, ecp=None, frag_type="chemgen", + opt="soscf", ): """Front-facing API bridge tailored for SCINE Puffin @@ -369,36 +370,56 @@ def be2puffin( "Using QM/MM Point Charges: Assuming QM structure in Angstrom " "and MM Coordinates in Bohr !!!" ) - mf1 = scf.UHF(mol).set( - max_cycle=200 - ) # using SOSCF is more reliable - # mf1 = scf.UHF(mol).set(max_cycle = 200, level_shift = (0.3, 0.2)) - # using level shift helps, but not always. level_shift and - # scf.addons.dynamic_level_shift do not seem to work with QM/MM - # note: from the SCINE database, the structure is in Angstrom but - # the MM point charges are in Bohr !! - mf = qmmm.mm_charge( - mf1, pts_and_charges[0], pts_and_charges[1], unit="bohr" - ).newton() # mf object, coordinates, charges + mf1 = scf.UHF(mol).set(max_cycle=200) + if opt.upper() == "SOSCF": + # using SOSCF is more reliable + # note: from the SCINE database, the structure is in Angstrom + # but the MM point charges are in Bohr !! + mf = qmmm.mm_charge( + mf1, pts_and_charges[0], pts_and_charges[1], unit="bohr" + ).newton() # mf object, coordinates, charges + elif opt.upper() == "DAMP": + mf = qmmm.mm_charge( + mf1, pts_and_charges[0], pts_and_charges[1], unit="bohr" + ) + mf.damp = 0.5 + mf.diis_start_cycle = 5 else: - mf = scf.UHF(mol).set(max_cycle=200, level_shift=(0.3, 0.2)) + if opt.upper() == "SOSCF": + mf = scf.UHF(mol).set(max_cycle=200).newton() + elif opt.upper() == "DAMP": + mf = scf.UHF(mol).set(max_cycle=200) + mf.damp = 0.5 + mf.diis_start_cycle = 5 else: - mf = scf.UHF(mol).set(max_cycle=200).newton() + if opt.upper() == "SOSCF": + mf = scf.UHF(mol).set(max_cycle=200).newton() + elif opt.upper() == "DAMP": + mf = scf.UHF(mol).set(max_cycle=200) + mf.damp = 0.5 + mf.diis_start_cycle = 5 else: # restricted if pts_and_charges: # running QM/MM print( "Using QM/MM Point Charges: Assuming QM structure in Angstrom and " "MM Coordinates in Bohr !!!" ) - mf1 = scf.RHF(mol).set(max_cycle=200) - mf = qmmm.mm_charge( - mf1, pts_and_charges[0], pts_and_charges[1], unit="bohr" - ).newton() if use_df or jk is not None: raise ValueError( "Setting use_df to false and jk to none: have not tested DF " "and QM/MM from point charges at the same time" ) + mf1 = scf.RHF(mol).set(max_cycle=200) + if opt.upper() == "SOSCF": + mf = qmmm.mm_charge( + mf1, pts_and_charges[0], pts_and_charges[1], unit="bohr" + ).newton() + elif opt.upper() == "DAMP": + mf = qmmm.mm_charge( + mf1, pts_and_charges[0], pts_and_charges[1], unit="bohr" + ) + mf.damp = 0.5 + mf.diis_start_cycle = 5 elif use_df and jk is None: mf = scf.RHF(mol).density_fit(auxbasis=df_aux_basis) else: @@ -464,12 +485,22 @@ def be2puffin( "Using QM/MM Point Charges: Assuming QM structure in Angstrom and " "MM Coordinates in Bohr !!!" ) - mf = qmmm.mm_charge( - mf, - pts_and_charges[0], - pts_and_charges[1], - unit="bohr", - ).newton() + if opt.upper() == "SOSCF": + mf = qmmm.mm_charge( + mf, + pts_and_charges[0], + pts_and_charges[1], + unit="bohr", + ).newton() + elif opt.upper() == "DAMP": + mf = qmmm.mm_charge( + mf, + pts_and_charges[0], + pts_and_charges[1], + unit="bohr", + ) + mf.damp = 0.5 + mf.diis_start_cycle = 5 time_post_mf = time.time() print("Chkfile electronic energy:", mf.energy_elec(), flush=True) print("Chkfile e_tot:", mf.e_tot, flush=True) @@ -496,7 +527,7 @@ def be2puffin( # Run oneshot embedding and return system energy mybe.oneshot(solver=solver, nproc=nproc, ompnum=ompnum) - return mybe.ebe_tot - mybe.ebe_hf + return mybe.ebe_tot - mybe.uhf_full_e def print_energy_cumulant(ecorr, e_V_Kapprox, e_F_dg, e_hf): diff --git a/src/quemb/molbe/pfrag.py b/src/quemb/molbe/pfrag.py index a556796e0d..7346fc1f1b 100644 --- a/src/quemb/molbe/pfrag.py +++ b/src/quemb/molbe/pfrag.py @@ -476,8 +476,8 @@ def schmidt_decomposition( if norb is not None: # add extra correlated orbital from environment # this will likely have Eval = 1 - more_corr = Bidx[-1] + 1 - Bidx.append(more_corr) + while len(Bidx) < norb: + Bidx.append(Bidx[-1] + 1) # Initialize the transformation matrix (TA) TA = zeros([Tot_sites, len(AO_in_frag) + len(Bidx)]) diff --git a/src/quemb/molbe/ube.py b/src/quemb/molbe/ube.py index 93f6699ce4..c89d9ba77a 100644 --- a/src/quemb/molbe/ube.py +++ b/src/quemb/molbe/ube.py @@ -227,9 +227,13 @@ def initialize(self, eri_, compute_hf): if self.equal_bath: # Enforce the same number of alpha and beta orbitals # by augmenting the bath + print("self.equal_bath true") tot_alpha = fobj_a.n_f + fobj_a.n_b tot_beta = fobj_b.n_f + fobj_b.n_b - + print("tot_alpha", tot_alpha) + print("tot_beta", tot_beta) + print("fobj_a.n_b,", fobj_a.n_b) + print("fobj_b.n_b,", fobj_b.n_b) if tot_alpha > tot_beta: fobj_b.sd( self.W[1] if self.frozen_core else self.W, @@ -246,7 +250,14 @@ def initialize(self, eri_, compute_hf): thr_bath=self.thr_bath, norb=fobj_b.n_b, ) - + print("tot_alphab", tot_alpha) + print("tot_betab", tot_beta) + print("fobj_a.n_b,", fobj_a.n_b) + print("fobj_b.n_b,", fobj_b.n_b) + tot_alpha = fobj_a.n_f + fobj_a.n_b + tot_beta = fobj_b.n_f + fobj_b.n_b + print("tot_alphac", tot_alpha) + print("tot_betac", tot_beta) assert fobj_a.TA is not None and fobj_b.TA is not None assert eri_ is not None, "eri_ is None: set incore_anyway for UHF" diff --git a/tests/ube-oneshot_test.py b/tests/ube-oneshot_test.py index 713c67d667..3fff71f60e 100644 --- a/tests/ube-oneshot_test.py +++ b/tests/ube-oneshot_test.py @@ -3,6 +3,12 @@ This tests for hexene anion and cation in minimal basis with and without frozen core. +Note: we can now run these without custom PySCF and avoiding Numpy errors +by adding extra bath orbitals. The current prescription can +generate slightly different baths and isn't totally deterministic, +so some tests have a larger delta for now. This may be modified in the +future. + Author(s): Leah Weisburn """ @@ -28,14 +34,16 @@ def test_hexene_anion_sto3g_frz_ben(self): mol.spin = 1 mol.build() self.molecular_unrestricted_oneshot_test( - mol, 1, "Hexene Anion Frz (BE1)", True, -0.35753374 + mol, 1, "Hexene Anion Frz (BE1)", True, -0.35753375 ) self.molecular_unrestricted_oneshot_test( - mol, 2, "Hexene Anion Frz (BE2)", True, -0.34725961 + mol, 2, "Hexene Anion Frz (BE2)", True, -0.34617257, delta=1e-4 ) + """ Cut for expense self.molecular_unrestricted_oneshot_test( - mol, 3, "Hexene Anion Frz (BE3)", True, -0.34300834 + mol, 3, "Hexene Anion Frz (BE3)", True, -0.34300832 ) + """ @unittest.skipUnless( os.getenv("QUEMB_DO_KNOWN_TO_FAIL_TESTS") == "true", @@ -50,14 +58,16 @@ def test_hexene_cation_sto3g_frz_ben(self): mol.spin = 1 mol.build() self.molecular_unrestricted_oneshot_test( - mol, 1, "Hexene Cation Frz (BE1)", True, -0.40383508 + mol, 1, "Hexene Cation Frz (BE1)", True, -0.40383505 ) self.molecular_unrestricted_oneshot_test( - mol, 2, "Hexene Cation Frz (BE2)", True, -0.36496690 + mol, 2, "Hexene Cation Frz (BE2)", True, -0.36736494, delta=1e-4 ) + """ Cut for expense self.molecular_unrestricted_oneshot_test( - mol, 3, "Hexene Cation Frz (BE3)", True, -0.36996484 + mol, 3, "Hexene Cation Frz (BE3)", True, -0.36996482 ) + """ @unittest.skipUnless( os.getenv("QUEMB_DO_KNOWN_TO_FAIL_TESTS") == "true", @@ -75,11 +85,13 @@ def test_hexene_anion_sto3g_unfrz_ben(self): mol, 1, "Hexene Anion Unfrz (BE1)", False, -0.38478279 ) self.molecular_unrestricted_oneshot_test( - mol, 2, "Hexene Anion Unfrz (BE2)", False, -0.39053689 + mol, 2, "Hexene Anion Unfrz (BE2)", False, -0.39053331, delta=1e-4 ) + """ Cut for expense self.molecular_unrestricted_oneshot_test( - mol, 3, "Hexene Anion Unfrz (BE3)", False, -0.38960174 + mol, 3, "Hexene Anion Unfrz (BE3)", False, -0.3895924 ) + """ @unittest.skipUnless( os.getenv("QUEMB_DO_KNOWN_TO_FAIL_TESTS") == "true", @@ -97,19 +109,21 @@ def test_hexene_cation_sto3g_unfrz_ben(self): mol, 1, "Hexene Cation Frz (BE1)", False, -0.39471433 ) self.molecular_unrestricted_oneshot_test( - mol, 2, "Hexene Cation Frz (BE2)", False, -0.39846777 + mol, 2, "Hexene Cation Frz (BE2)", False, -0.39846793, delta=1e-4 ) + """ Cut for expense self.molecular_unrestricted_oneshot_test( - mol, 3, "Hexene Cation Frz (BE3)", False, -0.39729184 + mol, 3, "Hexene Cation Frz (BE3)", False, -0.39729215 ) + """ def molecular_unrestricted_oneshot_test( - self, mol, n_BE, test_name, frz, exp_result, delta=1e-4 + self, mol, n_BE, test_name, frz, exp_result, delta=1e-5 ): mf = scf.UHF(mol) mf.kernel() fobj = fragmentate(frag_type="chemgen", n_BE=n_BE, mol=mol, frozen_core=frz) - mybe = UBE(mf, fobj) + mybe = UBE(mf, fobj, equal_bath=True) mybe.oneshot(solver="UCCSD", nproc=1) self.assertAlmostEqual( mybe.ebe_tot - mybe.uhf_full_e, From 7ce2e31b7a27ea31aed54a50fc24e27b3dda87d8 Mon Sep 17 00:00:00 2001 From: lweisburn Date: Tue, 17 Jun 2025 13:33:49 -0400 Subject: [PATCH 07/31] fix uhf energy change in misc --- src/quemb/molbe/misc.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/src/quemb/molbe/misc.py b/src/quemb/molbe/misc.py index f89aec7240..b076ff3cff 100644 --- a/src/quemb/molbe/misc.py +++ b/src/quemb/molbe/misc.py @@ -527,7 +527,10 @@ def be2puffin( # Run oneshot embedding and return system energy mybe.oneshot(solver=solver, nproc=nproc, ompnum=ompnum) - return mybe.ebe_tot - mybe.uhf_full_e + if unrestricted: + return mybe.ebe_tot - mybe.uhf_full_e + else: + return mybe.ebe_tot - mybe.ebe_hf def print_energy_cumulant(ecorr, e_V_Kapprox, e_F_dg, e_hf): From 9bb980aef6181ae76aa2efbbb182b6e3b67d542c Mon Sep 17 00:00:00 2001 From: lweisburn Date: Tue, 17 Jun 2025 14:02:54 -0400 Subject: [PATCH 08/31] get rid of redundant uhf_full_e and remove prints --- src/quemb/molbe/ube.py | 18 ++---------------- 1 file changed, 2 insertions(+), 16 deletions(-) diff --git a/src/quemb/molbe/ube.py b/src/quemb/molbe/ube.py index c89d9ba77a..ed5422760c 100644 --- a/src/quemb/molbe/ube.py +++ b/src/quemb/molbe/ube.py @@ -124,8 +124,6 @@ def __init__( self.P_core = None self.core_veff = None - self.uhf_full_e = mf.e_tot - if self.frozen_core: assert not ( fobj.ncore is None or fobj.no_core_idx is None or fobj.core_list is None @@ -227,13 +225,8 @@ def initialize(self, eri_, compute_hf): if self.equal_bath: # Enforce the same number of alpha and beta orbitals # by augmenting the bath - print("self.equal_bath true") tot_alpha = fobj_a.n_f + fobj_a.n_b tot_beta = fobj_b.n_f + fobj_b.n_b - print("tot_alpha", tot_alpha) - print("tot_beta", tot_beta) - print("fobj_a.n_b,", fobj_a.n_b) - print("fobj_b.n_b,", fobj_b.n_b) if tot_alpha > tot_beta: fobj_b.sd( self.W[1] if self.frozen_core else self.W, @@ -250,14 +243,7 @@ def initialize(self, eri_, compute_hf): thr_bath=self.thr_bath, norb=fobj_b.n_b, ) - print("tot_alphab", tot_alpha) - print("tot_betab", tot_beta) - print("fobj_a.n_b,", fobj_a.n_b) - print("fobj_b.n_b,", fobj_b.n_b) - tot_alpha = fobj_a.n_f + fobj_a.n_b - tot_beta = fobj_b.n_f + fobj_b.n_b - print("tot_alphac", tot_alpha) - print("tot_betac", tot_beta) + assert fobj_a.TA is not None and fobj_b.TA is not None assert eri_ is not None, "eri_ is None: set incore_anyway for UHF" @@ -413,7 +399,7 @@ def oneshot(self, solver="UCCSD", nproc=1, ompnum=4): print("-----------------------------------------------------", flush=True) print(flush=True) - self.ebe_tot = E + self.uhf_full_e + self.ebe_tot = E + self.hf_etot print( "Total Energy : {:>12.8f} Ha".format( (self.ebe_tot), From 560e681a4c5cf68a4aaf2a71bd7c80485c478ddd Mon Sep 17 00:00:00 2001 From: lweisburn Date: Tue, 17 Jun 2025 16:11:53 -0400 Subject: [PATCH 09/31] more standardized hf_etot --- src/quemb/molbe/misc.py | 2 +- tests/ube-oneshot_test.py | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/src/quemb/molbe/misc.py b/src/quemb/molbe/misc.py index b076ff3cff..9abab1b0be 100644 --- a/src/quemb/molbe/misc.py +++ b/src/quemb/molbe/misc.py @@ -528,7 +528,7 @@ def be2puffin( mybe.oneshot(solver=solver, nproc=nproc, ompnum=ompnum) if unrestricted: - return mybe.ebe_tot - mybe.uhf_full_e + return mybe.ebe_tot - mybe.hf_etot else: return mybe.ebe_tot - mybe.ebe_hf diff --git a/tests/ube-oneshot_test.py b/tests/ube-oneshot_test.py index 3fff71f60e..5e19838c21 100644 --- a/tests/ube-oneshot_test.py +++ b/tests/ube-oneshot_test.py @@ -126,12 +126,12 @@ def molecular_unrestricted_oneshot_test( mybe = UBE(mf, fobj, equal_bath=True) mybe.oneshot(solver="UCCSD", nproc=1) self.assertAlmostEqual( - mybe.ebe_tot - mybe.uhf_full_e, + mybe.ebe_tot - mybe.hf_etot, exp_result, msg="Unrestricted One-Shot Energy for " + test_name + " is incorrect by" - + str(mybe.ebe_tot - mybe.uhf_full_e - exp_result), + + str(mybe.ebe_tot - mybe.hf_etot - exp_result), delta=delta, ) From eb168deb12b0e2775939dcd13965108a65822f95 Mon Sep 17 00:00:00 2001 From: lweisburn Date: Wed, 25 Jun 2025 13:05:34 -0400 Subject: [PATCH 10/31] modify test targets --- src/quemb/molbe/ube.py | 4 +--- tests/ube-oneshot_test.py | 6 +++--- 2 files changed, 4 insertions(+), 6 deletions(-) diff --git a/src/quemb/molbe/ube.py b/src/quemb/molbe/ube.py index ed5422760c..fdbc80b292 100644 --- a/src/quemb/molbe/ube.py +++ b/src/quemb/molbe/ube.py @@ -116,7 +116,6 @@ def __init__( self.pot = initialize_pot(self.fobj.n_frag, self.fobj.relAO_per_edge_per_frag) self.eri_file = Path(eri_file) - self.ek = 0.0 self.frozen_core = fobj.frozen_core self.ncore = 0 self.E_core = 0 @@ -349,8 +348,7 @@ def initialize(self, eri_, compute_hf): ) if compute_hf: hf_err = self.hf_etot - (E_hf + self.enuc + self.E_core) - - self.ebe_hf = E_hf + self.enuc + self.E_core - self.ek + self.ebe_hf = E_hf + self.enuc + self.E_core print(f"HF-in-HF error : {hf_err:>.4e} Ha") if abs(hf_err) > 1.0e-5: warn("Large HF-in-HF energy error") diff --git a/tests/ube-oneshot_test.py b/tests/ube-oneshot_test.py index 5e19838c21..507c9592ba 100644 --- a/tests/ube-oneshot_test.py +++ b/tests/ube-oneshot_test.py @@ -37,7 +37,7 @@ def test_hexene_anion_sto3g_frz_ben(self): mol, 1, "Hexene Anion Frz (BE1)", True, -0.35753375 ) self.molecular_unrestricted_oneshot_test( - mol, 2, "Hexene Anion Frz (BE2)", True, -0.34617257, delta=1e-4 + mol, 2, "Hexene Anion Frz (BE2)", True, -0.34617685, delta=1e-4 ) """ Cut for expense self.molecular_unrestricted_oneshot_test( @@ -85,7 +85,7 @@ def test_hexene_anion_sto3g_unfrz_ben(self): mol, 1, "Hexene Anion Unfrz (BE1)", False, -0.38478279 ) self.molecular_unrestricted_oneshot_test( - mol, 2, "Hexene Anion Unfrz (BE2)", False, -0.39053331, delta=1e-4 + mol, 2, "Hexene Anion Unfrz (BE2)", False, -0.39052993, delta=1e-4 ) """ Cut for expense self.molecular_unrestricted_oneshot_test( @@ -109,7 +109,7 @@ def test_hexene_cation_sto3g_unfrz_ben(self): mol, 1, "Hexene Cation Frz (BE1)", False, -0.39471433 ) self.molecular_unrestricted_oneshot_test( - mol, 2, "Hexene Cation Frz (BE2)", False, -0.39846793, delta=1e-4 + mol, 2, "Hexene Cation Frz (BE2)", False, -0.39849056, delta=1e-4 ) """ Cut for expense self.molecular_unrestricted_oneshot_test( From ad440fcf29c4d06179ce7b83c573c6fb50735aa7 Mon Sep 17 00:00:00 2001 From: lweisburn Date: Wed, 23 Jul 2025 11:59:06 -0400 Subject: [PATCH 11/31] add schmidt_decomposition description --- src/quemb/molbe/pfrag.py | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/src/quemb/molbe/pfrag.py b/src/quemb/molbe/pfrag.py index 7346fc1f1b..0253a47d91 100644 --- a/src/quemb/molbe/pfrag.py +++ b/src/quemb/molbe/pfrag.py @@ -474,9 +474,16 @@ def schmidt_decomposition( # Set the number of orbitals to be taken from the environment orbitals # Based on an eigenvalue threshold ordering if norb is not None: - # add extra correlated orbital from environment + # add extra orbital from environment # this will likely have Eval = 1 + # note: there are normally very few orbitals with a Eval[i] <= thr_bath, + # so adding Bidx from the "front of the list" doesn't work. Instead, we add + # Bidx corresponding to a high eigenvalue from the environment + # (this is analagous to tightening up the threshold of the bath for the alpha + # or beta orbitals until they are the same size) while len(Bidx) < norb: + # Bidx corresponds to sorted Eval and Evec, so this simply adds indices + # corresponding to larger eigenvectors until the bath size reaches norb Bidx.append(Bidx[-1] + 1) # Initialize the transformation matrix (TA) From efa0394ef9cbc17cef2718b43866dd27e909278b Mon Sep 17 00:00:00 2001 From: Lea Northcote Date: Tue, 26 May 2026 12:05:50 -0400 Subject: [PATCH 12/31] fix bath augmentation out-of-bounds for open-shell systems --- src/quemb/molbe/pfrag.py | 31 ++++++++++++++++++++----------- 1 file changed, 20 insertions(+), 11 deletions(-) diff --git a/src/quemb/molbe/pfrag.py b/src/quemb/molbe/pfrag.py index 0253a47d91..7d87fca3c5 100644 --- a/src/quemb/molbe/pfrag.py +++ b/src/quemb/molbe/pfrag.py @@ -474,17 +474,26 @@ def schmidt_decomposition( # Set the number of orbitals to be taken from the environment orbitals # Based on an eigenvalue threshold ordering if norb is not None: - # add extra orbital from environment - # this will likely have Eval = 1 - # note: there are normally very few orbitals with a Eval[i] <= thr_bath, - # so adding Bidx from the "front of the list" doesn't work. Instead, we add - # Bidx corresponding to a high eigenvalue from the environment - # (this is analagous to tightening up the threshold of the bath for the alpha - # or beta orbitals until they are the same size) - while len(Bidx) < norb: - # Bidx corresponds to sorted Eval and Evec, so this simply adds indices - # corresponding to larger eigenvectors until the bath size reaches norb - Bidx.append(Bidx[-1] + 1) + # add extra orbital(s) from the environment; these will likely have + # Eval close to 1. note: there are normally very few orbitals with a + # Eval[i] <= thr_bath, so adding Bidx from the "front of the list" + # doesn't work. Instead, we add the excluded orbitals closest to the + # thr_bath/1-thr_bath boundary (this is analagous to tightening up + # the bath threshold for the alpha or beta orbitals until they are + # the same size) + excluded = [i for i in range(len(Eval)) if i not in set(Bidx)] + excluded_sorted = sorted( + excluded, + key=lambda i: min( + abs(Eval[i] - (1.0 - thr_bath)), abs(Eval[i] - thr_bath) + ), + ) + # Bidx corresponds to sorted Eval and Evec, so this adds indices + # closest to the bath threshold until the bath size reaches norb + for idx in excluded_sorted: + if len(Bidx) >= norb: + break + Bidx.append(idx) # Initialize the transformation matrix (TA) TA = zeros([Tot_sites, len(AO_in_frag) + len(Bidx)]) From 6888abb6e6affbcbed5bf044f1c9966cf45eb7d8 Mon Sep 17 00:00:00 2001 From: Lea Northcote Date: Wed, 27 May 2026 12:20:10 -0400 Subject: [PATCH 13/31] add urdm1_fullbasis for spin density assembly from BE-UCCSD fragments --- src/quemb/molbe/ube.py | 40 ++++++++++++++++++++++++++++++++++++++++ 1 file changed, 40 insertions(+) diff --git a/src/quemb/molbe/ube.py b/src/quemb/molbe/ube.py index fdbc80b292..2e546f4dc3 100644 --- a/src/quemb/molbe/ube.py +++ b/src/quemb/molbe/ube.py @@ -409,6 +409,46 @@ def oneshot(self, solver="UCCSD", nproc=1, ompnum=4): ) ) + def urdm1_fullbasis(self, return_ao=True): + """Assemble full-system alpha and beta 1-RDMs via democratic partitioning. + + Returns + ------- + rdm1a_AO, rdm1b_AO : numpy.ndarray + Alpha and beta 1-RDMs in the AO basis. + Spin density = rdm1a_AO - rdm1b_AO. + """ + from numpy import zeros + + nao = self.S.shape[0] + rdm1a_AO = zeros((nao, nao)) + rdm1b_AO = zeros((nao, nao)) + + W = self.W # localization matrix, unfrozen core + + for fobj_a, fobj_b in zip(self.Fobjs_a, self.Fobjs_b): + # Fragment AO centers - same for alpha and beta + cind = [fobj_a.AO_in_frag[i] + for i in fobj_a.weight_and_relAO_per_center[1]] + + # Democratic partitioning projection in full AO space + Proj = self.S @ self.W[:, cind] @ self.W[:, cind].T @ self.S + + # Transform fragment RDM to full AO space first + rdm1a_full = fobj_a.TA @ fobj_a.mo_coeffs @ fobj_a.rdm1__ \ + @ fobj_a.mo_coeffs.T @ fobj_a.TA.T + rdm1b_full = fobj_b.TA @ fobj_b.mo_coeffs @ fobj_b.rdm1__ \ + @ fobj_b.mo_coeffs.T @ fobj_b.TA.T + + # Apply democratic weight in full AO space + rdm1a_AO += Proj @ rdm1a_full + rdm1b_AO += Proj @ rdm1b_full + + # Symmetrize + rdm1a_AO = (rdm1a_AO + rdm1a_AO.T) / 2.0 + rdm1b_AO = (rdm1b_AO + rdm1b_AO.T) / 2.0 + + return rdm1a_AO, rdm1b_AO def initialize_pot(n_frag, relAO_per_edge): pot_ = [] From fabaf6da16dcb3c743d733d20523be41bb865b01 Mon Sep 17 00:00:00 2001 From: Lea Northcote Date: Wed, 3 Jun 2026 22:00:05 -0400 Subject: [PATCH 14/31] save mo_coeff_uccsd on fragment after UCCSD for RDM assembly --- src/quemb/molbe/solver.py | 2 ++ src/quemb/molbe/ube.py | 52 +++++++++++++++++++++++++++++---------- 2 files changed, 41 insertions(+), 13 deletions(-) diff --git a/src/quemb/molbe/solver.py b/src/quemb/molbe/solver.py index cc96312648..9c049cb8bf 100644 --- a/src/quemb/molbe/solver.py +++ b/src/quemb/molbe/solver.py @@ -630,6 +630,8 @@ def be_func_u( raise ValueError("Solver not implemented") assert fobj_a._mf is not None and fobj_b._mf is not None + fobj_a.mo_coeff_uccsd = fobj_a._mf.mo_coeff.copy() # save before _mf can be cleared + fobj_b.mo_coeff_uccsd = fobj_b._mf.mo_coeff.copy() fobj_a.rdm1__ = rdm1_tmp[0].copy() fobj_b._rdm1 = ( multi_dot((fobj_a._mf.mo_coeff, rdm1_tmp[0], fobj_a._mf.mo_coeff.T)) * 0.5 diff --git a/src/quemb/molbe/ube.py b/src/quemb/molbe/ube.py index 2e546f4dc3..eb0db919b0 100644 --- a/src/quemb/molbe/ube.py +++ b/src/quemb/molbe/ube.py @@ -365,7 +365,7 @@ def initialize(self, eri_, compute_hf): fobj.udim = couti couti = fobj.set_udim(couti) - def oneshot(self, solver="UCCSD", nproc=1, ompnum=4): + def oneshot(self, solver="UCCSD", nproc=1, ompnum=4, relax_density=False): if nproc == 1: E, E_comp = be_func_u( None, @@ -374,7 +374,7 @@ def oneshot(self, solver="UCCSD", nproc=1, ompnum=4): self.enuc, hf_veff=self.hf_veff, eeval=True, - relax_density=False, + relax_density=relax_density, frozen=self.frozen_core, ) else: @@ -386,7 +386,7 @@ def oneshot(self, solver="UCCSD", nproc=1, ompnum=4): hf_veff=self.hf_veff, nproc=nproc, ompnum=ompnum, - relax_density=False, + relax_density=relax_density, frozen=self.frozen_core, ) unused(E_comp) @@ -409,7 +409,7 @@ def oneshot(self, solver="UCCSD", nproc=1, ompnum=4): ) ) - def urdm1_fullbasis(self, return_ao=True): + def urdm1_fullbasis(self): """Assemble full-system alpha and beta 1-RDMs via democratic partitioning. Returns @@ -425,24 +425,50 @@ def urdm1_fullbasis(self, return_ao=True): rdm1b_AO = zeros((nao, nao)) W = self.W # localization matrix, unfrozen core - + + def get_mo(fobj): + if hasattr(fobj, 'mo_coeff_uccsd'): + return fobj.mo_coeff_uccsd + if fobj._mf is not None: + return fobj._mf.mo_coeff + return fobj._mo_coeffs + for fobj_a, fobj_b in zip(self.Fobjs_a, self.Fobjs_b): # Fragment AO centers - same for alpha and beta cind = [fobj_a.AO_in_frag[i] for i in fobj_a.weight_and_relAO_per_center[1]] - # Democratic partitioning projection in full AO space - Proj = self.S @ self.W[:, cind] @ self.W[:, cind].T @ self.S + # Democratic partitioning projection in full AO space, built + # per-spin since self.W is split into [Wa, Wb] under frozen_core + Wa = self.W[0] if self.frozen_core else self.W + Wb = self.W[1] if self.frozen_core else self.W + Proj_a = self.S @ Wa[:, cind] @ Wa[:, cind].T @ self.S + Proj_b = self.S @ Wb[:, cind] @ Wb[:, cind].T @ self.S + + # Get MO coefficients — _mf.mo_coeff is most reliable when available + if fobj_a._mf is not None: + mca = fobj_a._mf.mo_coeff + elif fobj_a.mo_coeffs is not None: + mca = fobj_a.mo_coeffs + else: + mca = fobj_a._mo_coeffs + + if fobj_b._mf is not None: + mcb = fobj_b._mf.mo_coeff + elif fobj_b.mo_coeffs is not None: + mcb = fobj_b.mo_coeffs + else: + mcb = fobj_b._mo_coeffs # Transform fragment RDM to full AO space first - rdm1a_full = fobj_a.TA @ fobj_a.mo_coeffs @ fobj_a.rdm1__ \ - @ fobj_a.mo_coeffs.T @ fobj_a.TA.T - rdm1b_full = fobj_b.TA @ fobj_b.mo_coeffs @ fobj_b.rdm1__ \ - @ fobj_b.mo_coeffs.T @ fobj_b.TA.T + rdm1a_full = fobj_a.TA @ fobj_a._mf.mo_coeff @ fobj_a.rdm1__ \ + @ fobj_a._mf.mo_coeff.T @ fobj_a.TA.T + rdm1b_full = fobj_b.TA @ fobj_b._mf.mo_coeff @ fobj_b.rdm1__ \ + @ fobj_b._mf.mo_coeff.T @ fobj_b.TA.T # Apply democratic weight in full AO space - rdm1a_AO += Proj @ rdm1a_full - rdm1b_AO += Proj @ rdm1b_full + rdm1a_AO += Proj_a @ rdm1a_full + rdm1b_AO += Proj_b @ rdm1b_full # Symmetrize rdm1a_AO = (rdm1a_AO + rdm1a_AO.T) / 2.0 From 62410d7b871ff93004d8f69656c4b543d343cc9b Mon Sep 17 00:00:00 2001 From: Lea Northcote Date: Tue, 16 Jun 2026 16:44:05 -0400 Subject: [PATCH 15/31] UBE: Parallel RDM propagation, DF support, bug fixes --- src/quemb/molbe/be_parallel.py | 21 +++++++---- src/quemb/molbe/solver.py | 5 ++- src/quemb/molbe/ube.py | 65 +++++++++++++++------------------- 3 files changed, 44 insertions(+), 47 deletions(-) diff --git a/src/quemb/molbe/be_parallel.py b/src/quemb/molbe/be_parallel.py index 5532fb8b73..b011070252 100644 --- a/src/quemb/molbe/be_parallel.py +++ b/src/quemb/molbe/be_parallel.py @@ -407,7 +407,7 @@ def run_solver_u( gcores=full_uhf.full_gcore, frozen=frozen, ) - return e_f + return e_f, rdm1_tmp, fobj_a._mf.mo_coeff.copy(), fobj_b._mf.mo_coeff.copy() def be_func_parallel( @@ -623,14 +623,21 @@ def be_func_parallel_u( ) results.append(result) - energy_list = [result.get() for result in results] - + results_list = [result.get() for result in results] + # Store RDMs back into fragment objects + for i, (fobj_a, fobj_b) in enumerate(Fobjs): + e_f, rdm1_tmp, mo_a, mo_b = results_list[i] + fobj_a.rdm1__ = rdm1_tmp[0].copy() + fobj_b.rdm1__ = rdm1_tmp[1].copy() + fobj_a.mo_coeff_uccsd = mo_a + fobj_b.mo_coeff_uccsd = mo_b # Compute and return fragment energy e_1 = 0.0 e_2 = 0.0 e_c = 0.0 - for i in range(len(energy_list)): - e_1 += energy_list[i][0] - e_2 += energy_list[i][1] - e_c += energy_list[i][2] + for i in range(len(results_list)): + e_f = results_list[i][0] + e_1 += e_f[0] + e_2 += e_f[1] + e_c += e_f[2] return (e_1 + e_2 + e_c, (e_1, e_2, e_c)) diff --git a/src/quemb/molbe/solver.py b/src/quemb/molbe/solver.py index 9c049cb8bf..332f797f87 100644 --- a/src/quemb/molbe/solver.py +++ b/src/quemb/molbe/solver.py @@ -630,13 +630,12 @@ def be_func_u( raise ValueError("Solver not implemented") assert fobj_a._mf is not None and fobj_b._mf is not None - fobj_a.mo_coeff_uccsd = fobj_a._mf.mo_coeff.copy() # save before _mf can be cleared + fobj_a.mo_coeff_uccsd = fobj_a._mf.mo_coeff.copy() fobj_b.mo_coeff_uccsd = fobj_b._mf.mo_coeff.copy() fobj_a.rdm1__ = rdm1_tmp[0].copy() - fobj_b._rdm1 = ( + fobj_a._rdm1 = ( multi_dot((fobj_a._mf.mo_coeff, rdm1_tmp[0], fobj_a._mf.mo_coeff.T)) * 0.5 ) - fobj_b.rdm1__ = rdm1_tmp[1].copy() fobj_b._rdm1 = ( multi_dot((fobj_b._mf.mo_coeff, rdm1_tmp[1], fobj_b._mf.mo_coeff.T)) * 0.5 diff --git a/src/quemb/molbe/ube.py b/src/quemb/molbe/ube.py index eb0db919b0..fb08780ebc 100644 --- a/src/quemb/molbe/ube.py +++ b/src/quemb/molbe/ube.py @@ -46,6 +46,7 @@ def __init__( compute_hf: bool = True, thr_bath: float = 1.0e-10, equal_bath: bool = True, + use_df: bool = False, ) -> None: """Initialize Unrestricted BE Object (ube🍠) @@ -84,6 +85,9 @@ def __init__( self.unrestricted = True self.thr_bath = thr_bath self.equal_bath = equal_bath + self.use_df = use_df + if use_df: + assert hasattr(mf, "with_df") and mf.with_df is not None, "use_df=True requires a density-fitted mf: construct as scf.UHF(mol).density_fit()" self.fobj = fobj @@ -175,7 +179,7 @@ def __init__( self.scratch_dir = scratch_dir self.eri_file = self.scratch_dir / eri_file - self.initialize(mf._eri, compute_hf) + self.initialize(None if self.use_df else mf._eri, compute_hf) def initialize(self, eri_, compute_hf): if compute_hf: @@ -244,18 +248,20 @@ def initialize(self, eri_, compute_hf): ) assert fobj_a.TA is not None and fobj_b.TA is not None - assert eri_ is not None, "eri_ is None: set incore_anyway for UHF" - - eri_a = ao2mo.incore.full(eri_, fobj_a.TA, compact=True) - eri_b = ao2mo.incore.full(eri_, fobj_b.TA, compact=True) - - Csd_A = fobj_a.TA # may have to add in nibath here - Csd_B = fobj_b.TA - - # cross-spin ERI term - eri_ab = ao2mo.incore.general( - eri_, (Csd_A, Csd_A, Csd_B, Csd_B), compact=True - ) + if self.use_df: + eri_a = self.mf.with_df.ao2mo(fobj_a.TA, compact=True) + eri_b = self.mf.with_df.ao2mo(fobj_b.TA, compact=True) + eri_ab = self.mf.with_df.ao2mo( + (fobj_a.TA, fobj_a.TA, fobj_b.TA, fobj_b.TA), compact=True + ) + else: + assert eri_ is not None, "eri_ is None: set incore_anyway for UHF" + eri_a = ao2mo.incore.full(eri_, fobj_a.TA, compact=True) + eri_b = ao2mo.incore.full(eri_, fobj_b.TA, compact=True) + # cross-spin ERI term + eri_ab = ao2mo.incore.general( + eri_, (fobj_a.TA, fobj_a.TA, fobj_b.TA, fobj_b.TA), compact=True + ) file_eri.create_dataset(fobj_a.dname[0], data=eri_a) file_eri.create_dataset(fobj_a.dname[1], data=eri_b) @@ -411,7 +417,7 @@ def oneshot(self, solver="UCCSD", nproc=1, ompnum=4, relax_density=False): def urdm1_fullbasis(self): """Assemble full-system alpha and beta 1-RDMs via democratic partitioning. - + Returns ------- rdm1a_AO, rdm1b_AO : numpy.ndarray @@ -419,15 +425,13 @@ def urdm1_fullbasis(self): Spin density = rdm1a_AO - rdm1b_AO. """ from numpy import zeros - + nao = self.S.shape[0] rdm1a_AO = zeros((nao, nao)) rdm1b_AO = zeros((nao, nao)) - W = self.W # localization matrix, unfrozen core - def get_mo(fobj): - if hasattr(fobj, 'mo_coeff_uccsd'): + if hasattr(fobj, "mo_coeff_uccsd"): return fobj.mo_coeff_uccsd if fobj._mf is not None: return fobj._mf.mo_coeff @@ -445,27 +449,13 @@ def get_mo(fobj): Proj_a = self.S @ Wa[:, cind] @ Wa[:, cind].T @ self.S Proj_b = self.S @ Wb[:, cind] @ Wb[:, cind].T @ self.S - # Get MO coefficients — _mf.mo_coeff is most reliable when available - if fobj_a._mf is not None: - mca = fobj_a._mf.mo_coeff - elif fobj_a.mo_coeffs is not None: - mca = fobj_a.mo_coeffs - else: - mca = fobj_a._mo_coeffs - - if fobj_b._mf is not None: - mcb = fobj_b._mf.mo_coeff - elif fobj_b.mo_coeffs is not None: - mcb = fobj_b.mo_coeffs - else: - mcb = fobj_b._mo_coeffs + mca = get_mo(fobj_a) + mcb = get_mo(fobj_b) # Transform fragment RDM to full AO space first - rdm1a_full = fobj_a.TA @ fobj_a._mf.mo_coeff @ fobj_a.rdm1__ \ - @ fobj_a._mf.mo_coeff.T @ fobj_a.TA.T - rdm1b_full = fobj_b.TA @ fobj_b._mf.mo_coeff @ fobj_b.rdm1__ \ - @ fobj_b._mf.mo_coeff.T @ fobj_b.TA.T - + rdm1a_full = fobj_a.TA @ mca @ fobj_a.rdm1__ @ mca.T @ fobj_a.TA.T + rdm1b_full = fobj_b.TA @ mcb @ fobj_b.rdm1__ @ mcb.T @ fobj_b.TA.T + # Apply democratic weight in full AO space rdm1a_AO += Proj_a @ rdm1a_full rdm1b_AO += Proj_b @ rdm1b_full @@ -476,6 +466,7 @@ def get_mo(fobj): return rdm1a_AO, rdm1b_AO + def initialize_pot(n_frag, relAO_per_edge): pot_ = [] From 5e155e209d2178a72b0c63c9efc5e7fecf964462 Mon Sep 17 00:00:00 2001 From: Lea Northcote Date: Tue, 16 Jun 2026 16:46:21 -0400 Subject: [PATCH 16/31] scripts: Add ube_hfcc as console entry point --- pyproject.toml | 3 + src/quemb/scripts/__init__.py | 0 src/quemb/scripts/ube_hfcc.py | 198 ++++++++++++++++++++++++++++++++++ 3 files changed, 201 insertions(+) create mode 100644 src/quemb/scripts/__init__.py create mode 100644 src/quemb/scripts/ube_hfcc.py diff --git a/pyproject.toml b/pyproject.toml index 520599236f..c0ca7b8394 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -46,3 +46,6 @@ where = ["./src"] include = ["quemb*"] exclude = [] namespaces = false + +[project.scripts] +ube-hfcc = "quemb.scripts.ube_hfcc:main" diff --git a/src/quemb/scripts/__init__.py b/src/quemb/scripts/__init__.py new file mode 100644 index 0000000000..e69de29bb2 diff --git a/src/quemb/scripts/ube_hfcc.py b/src/quemb/scripts/ube_hfcc.py new file mode 100644 index 0000000000..b449171984 --- /dev/null +++ b/src/quemb/scripts/ube_hfcc.py @@ -0,0 +1,198 @@ +# Author(s): Lea Kjaergaard + +""" +Compute the full (isotropic + anisotropic) hyperfine coupling tensor from a +BE-UCCSD spin density. + +Built from PySCF's stable, core one-electron integrals (mol.intor). +The tensor comes from a single set of second-derivative-of-1/r integrals, +int1e_ipiprinv + int1e_iprinvip, which together give . +Via the Poisson equation, nabla^2(1/r) = -4*pi*delta(r), so this object's +trace channel already is the isotropic Fermi-contact delta-function operator, +and its traceless part is the spin-dipolar operator. This comes directly from +the integral construction. + +Usage: + python ube_hfcc.py --xyz molecule.xyz --rdm1a rdm1a.npy --rdm1b rdm1b.npy + --charge -2 --spin 5 --basis def2-svp [--atoms Fe C] + +Example: + python ube_hfcc.py --xyz model_II.xyz --rdm1a model_II_rdm1a_BE2.npy + --rdm1b model_II_rdm1b_BE2.npy --charge -2 --spin 5 +""" + +import argparse +import numpy as np +from pyscf import gto +from pyscf.data import nist +from pyscf.data.gyro import get_nuc_g_factor + + +def _fcdip_integral(mol, atom_id): + nao = mol.nao + with mol.with_rinv_origin(mol.atom_coord(atom_id)): + ipipv = mol.intor("int1e_ipiprinv", 9).reshape(3, 3, nao, nao) + ipvip = mol.intor("int1e_iprinvip", 9).reshape(3, 3, nao, nao) + h1ao = ipipv + ipvip + h1ao = h1ao + h1ao.transpose(0, 1, 3, 2) + trace = h1ao[0, 0] + h1ao[1, 1] + h1ao[2, 2] + idx = np.arange(3) + h1ao[idx, idx] -= trace + return h1ao + +def hfc_tensor(mol, atom_id, dma, dmb): + """Full (Fermi-contact + spin-dipolar) hyperfine coupling tensor for + one nucleus, in MHz. + + Parameters + ---------- + mol : pyscf.gto.Mole + atom_id : int + 0-indexed atom for which to compute the tensor. + dma, dmb : numpy.ndarray, shape (nao, nao) + Alpha and beta density matrices in the AO basis. + + Returns + ------- + numpy.ndarray, shape (3, 3) + The total HFC tensor in MHz. + """ + spindm = dma - dmb + effspin = mol.spin * 0.5 + + e_gyro = 0.5 * nist.G_ELECTRON + nuc_mag = 0.5 * (nist.E_MASS / nist.PROTON_MASS) + au2MHz = nist.HARTREE2J / nist.PLANCK * 1e-6 + fac = nist.ALPHA**2 / 2 / effspin * e_gyro * au2MHz + nuc_gyro = get_nuc_g_factor(mol.atom_symbol(atom_id)) * nuc_mag + + h1_fcdip = _fcdip_integral(mol, atom_id) + fcsd = np.einsum("xyij,ji->xy", h1_fcdip, spindm) + + return fac * nuc_gyro * fcsd + +def hfc_principal_values(tensor): + """Diagonalize a (3,3) HFC tensor. Returns (A_principal sorted + ascending, A_iso = trace/3).""" + evals = np.linalg.eigvalsh(0.5 * (tensor + tensor.T)) + a_iso = float(np.trace(tensor)) / 3.0 + return evals, a_iso + +def compute_hfcc(mol, dma, dmb, atoms=None): + """Compute the full (isotropic + anisotropic) HFC tensor for a set of + atoms from alpha/beta density matrices in the AO basis. + + Parameters + ---------- + mol : pyscf.gto.Mole + dma, dmb : numpy.ndarray, shape (nao, nao) + atoms : list of int, optional. Default is all atoms. + + Returns + ------- + list of dict with keys : index, symbol, tensor, evals, a_iso + """ + if atoms is None: + atoms = range(mol.natm) + + results = [] + for i in atoms: + symbol = mol.atom_symbol(i) + if get_nuc_g_factor(symbol) == 0: + continue + + tensor = hfc_tensor(mol, i, dma, dmb) + evals, a_iso = hfc_principal_values(tensor) + results.append( + { + "index": i, + "symbol": symbol, + "tensor": tensor, + "evals": evals, + "a_iso": a_iso, + } + ) + + return results + +def print_hfcc_table(results, title="Hyperfine Coupling Constants"): + print(f"\n{title}") + print("=" * 56) + print( + f"{'Atom':>4} {'Symbol':>6} {'A_1':>10} {'A_2':>10} " + f"{'A_3':>10} {'A_iso':>10}" + ) + print("-" * 56) + for r in results: + e = r["evals"] + print( + f"{r['index']:>4} {r['symbol']:>6} " + f"{e[0]:>10.4f} {e[1]:>10.4f} {e[2]:>10.4f} {r['a_iso']:>10.4f}" + ) + print("=" * 56) + +def main(): + parser = argparse.ArgumentParser( + description="Compute the full HFC tensor from a BE-UCCSD spin density" + ) + parser.add_argument("--xyz", required=True, help="XYZ geometry file") + parser.add_argument("--name", default=None, help="Shortcut: sets --rdm1a to NAME_rdm1a.npy and --rdm1b to NAME_rdm1b.npy") + parser.add_argument("--rdm1a", required=False, default=None, help="Alpha 1-RDM .npy file") + parser.add_argument("--rdm1b", required=False, default=None, help="Beta 1-RDM .npy file") + parser.add_argument("--charge", required=True, type=int, help="Molecular charge") + parser.add_argument("--spin", required=True, type=int, help="2S (number of unpaired electrons)") + parser.add_argument("--basis", default="def2-svp", help="Basis set (default: def2-svp)") + parser.add_argument("--atoms", nargs="+", help="Atom symbols to print e.g. Fe C H (default: all)") + parser.add_argument("--unit", default="angstrom", help="Coordinate unit (default: angstrom)") + args = parser.parse_args() + if args.name is not None: + if args.rdm1a is None: + args.rdm1a = f"{args.name}_rdm1a.npy" + if args.rdm1b is None: + args.rdm1b = f"{args.name}_rdm1b.npy" + + # Build molecule + mol = gto.M() + with open(args.xyz) as f: + lines = f.readlines() + # Handle both raw xyz (no header) and standard xyz (2-line header) + try: + int(lines[0].strip()) + mol.atom = ''.join(lines[2:]) # standard xyz with natom + comment lines + except ValueError: + mol.atom = ''.join(lines) # raw xyz with no header + mol.basis = args.basis + mol.charge = args.charge + mol.spin = args.spin + mol.unit = args.unit + mol.build() + + print(f"Molecule: {args.xyz}", flush=True) + print(f"Basis: {args.basis}, charge={args.charge}, spin={args.spin}", flush=True) + print(f"nao: {mol.nao_nr()}, nelec: {mol.nelec}", flush=True) + + # Load RDMs + rdm1a = np.load(args.rdm1a) + rdm1b = np.load(args.rdm1b) + spin_density = rdm1a - rdm1b + + print(f"\nRDM validation:", flush=True) + S = mol.intor('int1e_ovlp') + print(f" Trace rdm1a: {np.trace(rdm1a @ S):.4f} (expected {mol.nelec[0]})", flush=True) + print(f" Trace rdm1b: {np.trace(rdm1b @ S):.4f} (expected {mol.nelec[1]})", flush=True) + print(f" Net spin: {np.trace(spin_density @ S):.4f} (expected {mol.spin})", flush=True) + + # Filter atoms by symbol if requested + if args.atoms: + atom_indices = [i for i in range(mol.natm) + if mol.atom_symbol(i) in args.atoms] + print(f"\nComputing HFCCs for atoms: {args.atoms}", flush=True) + else: + atom_indices = None + print(f"\nComputing HFCCs for all atoms", flush=True) + + results = compute_hfcc(mol, rdm1a, rdm1b, atoms=atom_indices) + print_hfcc_table(results, title=f"BE-UCCSD Hyperfine Coupling ({args.basis})") + +if __name__ == "__main__": + main() From d4407f14240476e17450f11b77a48523f60ce30b Mon Sep 17 00:00:00 2001 From: Lea Northcote Date: Mon, 22 Jun 2026 13:58:22 -0400 Subject: [PATCH 17/31] fix: convert zip to list in oneshot parallel path so rdm1__ is stored correctly --- src/quemb/molbe/ube.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/quemb/molbe/ube.py b/src/quemb/molbe/ube.py index fb08780ebc..a53103b9b8 100644 --- a/src/quemb/molbe/ube.py +++ b/src/quemb/molbe/ube.py @@ -386,7 +386,7 @@ def oneshot(self, solver="UCCSD", nproc=1, ompnum=4, relax_density=False): else: E, E_comp = be_func_parallel_u( pot=None, - Fobjs=zip(self.Fobjs_a, self.Fobjs_b), + Fobjs=list(zip(self.Fobjs_a, self.Fobjs_b)), solver=solver, enuc=self.enuc, hf_veff=self.hf_veff, From 77ca4fb502f74a38f7df1f55f66428761d28b6d1 Mon Sep 17 00:00:00 2001 From: Lea Northcote Date: Fri, 14 Aug 2026 11:53:50 -0400 Subject: [PATCH 18/31] Fix no-op os.system() call in be_func_parallel and be_func_parallel_u: OMP_NUM_THREADS must be set via os.environ in the parent process before Pool creation for forked workers to inherit it --- src/quemb/molbe/be_parallel.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/quemb/molbe/be_parallel.py b/src/quemb/molbe/be_parallel.py index b011070252..d2c911c99f 100644 --- a/src/quemb/molbe/be_parallel.py +++ b/src/quemb/molbe/be_parallel.py @@ -473,7 +473,7 @@ def be_func_parallel( the error norm, error vector, and the computed energy. """ # Set the number of OpenMP threads - os.system("export OMP_NUM_THREADS=" + str(ompnum)) + os.environ["OMP_NUM_THREADS"] = str(ompnum) nprocs = nproc // ompnum # Update the effective Hamiltonian with potentials @@ -603,7 +603,7 @@ def be_func_parallel_u( Returns the computed energy """ # Set the number of OpenMP threads - os.system("export OMP_NUM_THREADS=" + str(ompnum)) + os.environ["OMP_NUM_THREADS"] = str(ompnum) nprocs = nproc // ompnum with ProcessPool(nprocs) as pool_: From 62107b4e46305fa5a2a83aa049b321c6bc76d07e Mon Sep 17 00:00:00 2001 From: Lea Northcote Date: Fri, 14 Aug 2026 16:27:08 -0400 Subject: [PATCH 19/31] Apply ruff format to pass formatting check --- src/quemb/molbe/pfrag.py | 4 +-- src/quemb/molbe/ube.py | 7 ++-- src/quemb/scripts/ube_hfcc.py | 67 +++++++++++++++++++++++++---------- 3 files changed, 53 insertions(+), 25 deletions(-) diff --git a/src/quemb/molbe/pfrag.py b/src/quemb/molbe/pfrag.py index 7d87fca3c5..76b23d47be 100644 --- a/src/quemb/molbe/pfrag.py +++ b/src/quemb/molbe/pfrag.py @@ -484,9 +484,7 @@ def schmidt_decomposition( excluded = [i for i in range(len(Eval)) if i not in set(Bidx)] excluded_sorted = sorted( excluded, - key=lambda i: min( - abs(Eval[i] - (1.0 - thr_bath)), abs(Eval[i] - thr_bath) - ), + key=lambda i: min(abs(Eval[i] - (1.0 - thr_bath)), abs(Eval[i] - thr_bath)), ) # Bidx corresponds to sorted Eval and Evec, so this adds indices # closest to the bath threshold until the bath size reaches norb diff --git a/src/quemb/molbe/ube.py b/src/quemb/molbe/ube.py index a53103b9b8..2ce9acc2ca 100644 --- a/src/quemb/molbe/ube.py +++ b/src/quemb/molbe/ube.py @@ -87,7 +87,9 @@ def __init__( self.equal_bath = equal_bath self.use_df = use_df if use_df: - assert hasattr(mf, "with_df") and mf.with_df is not None, "use_df=True requires a density-fitted mf: construct as scf.UHF(mol).density_fit()" + assert hasattr(mf, "with_df") and mf.with_df is not None, ( + "use_df=True requires a density-fitted mf: construct as scf.UHF(mol).density_fit()" + ) self.fobj = fobj @@ -439,8 +441,7 @@ def get_mo(fobj): for fobj_a, fobj_b in zip(self.Fobjs_a, self.Fobjs_b): # Fragment AO centers - same for alpha and beta - cind = [fobj_a.AO_in_frag[i] - for i in fobj_a.weight_and_relAO_per_center[1]] + cind = [fobj_a.AO_in_frag[i] for i in fobj_a.weight_and_relAO_per_center[1]] # Democratic partitioning projection in full AO space, built # per-spin since self.W is split into [Wa, Wb] under frozen_core diff --git a/src/quemb/scripts/ube_hfcc.py b/src/quemb/scripts/ube_hfcc.py index b449171984..069d25aafd 100644 --- a/src/quemb/scripts/ube_hfcc.py +++ b/src/quemb/scripts/ube_hfcc.py @@ -40,6 +40,7 @@ def _fcdip_integral(mol, atom_id): h1ao[idx, idx] -= trace return h1ao + def hfc_tensor(mol, atom_id, dma, dmb): """Full (Fermi-contact + spin-dipolar) hyperfine coupling tensor for one nucleus, in MHz. @@ -71,6 +72,7 @@ def hfc_tensor(mol, atom_id, dma, dmb): return fac * nuc_gyro * fcsd + def hfc_principal_values(tensor): """Diagonalize a (3,3) HFC tensor. Returns (A_principal sorted ascending, A_iso = trace/3).""" @@ -78,6 +80,7 @@ def hfc_principal_values(tensor): a_iso = float(np.trace(tensor)) / 3.0 return evals, a_iso + def compute_hfcc(mol, dma, dmb, atoms=None): """Compute the full (isotropic + anisotropic) HFC tensor for a set of atoms from alpha/beta density matrices in the AO basis. @@ -115,12 +118,12 @@ def compute_hfcc(mol, dma, dmb, atoms=None): return results + def print_hfcc_table(results, title="Hyperfine Coupling Constants"): print(f"\n{title}") print("=" * 56) print( - f"{'Atom':>4} {'Symbol':>6} {'A_1':>10} {'A_2':>10} " - f"{'A_3':>10} {'A_iso':>10}" + f"{'Atom':>4} {'Symbol':>6} {'A_1':>10} {'A_2':>10} {'A_3':>10} {'A_iso':>10}" ) print("-" * 56) for r in results: @@ -131,19 +134,36 @@ def print_hfcc_table(results, title="Hyperfine Coupling Constants"): ) print("=" * 56) + def main(): parser = argparse.ArgumentParser( description="Compute the full HFC tensor from a BE-UCCSD spin density" ) - parser.add_argument("--xyz", required=True, help="XYZ geometry file") - parser.add_argument("--name", default=None, help="Shortcut: sets --rdm1a to NAME_rdm1a.npy and --rdm1b to NAME_rdm1b.npy") - parser.add_argument("--rdm1a", required=False, default=None, help="Alpha 1-RDM .npy file") - parser.add_argument("--rdm1b", required=False, default=None, help="Beta 1-RDM .npy file") - parser.add_argument("--charge", required=True, type=int, help="Molecular charge") - parser.add_argument("--spin", required=True, type=int, help="2S (number of unpaired electrons)") - parser.add_argument("--basis", default="def2-svp", help="Basis set (default: def2-svp)") - parser.add_argument("--atoms", nargs="+", help="Atom symbols to print e.g. Fe C H (default: all)") - parser.add_argument("--unit", default="angstrom", help="Coordinate unit (default: angstrom)") + parser.add_argument("--xyz", required=True, help="XYZ geometry file") + parser.add_argument( + "--name", + default=None, + help="Shortcut: sets --rdm1a to NAME_rdm1a.npy and --rdm1b to NAME_rdm1b.npy", + ) + parser.add_argument( + "--rdm1a", required=False, default=None, help="Alpha 1-RDM .npy file" + ) + parser.add_argument( + "--rdm1b", required=False, default=None, help="Beta 1-RDM .npy file" + ) + parser.add_argument("--charge", required=True, type=int, help="Molecular charge") + parser.add_argument( + "--spin", required=True, type=int, help="2S (number of unpaired electrons)" + ) + parser.add_argument( + "--basis", default="def2-svp", help="Basis set (default: def2-svp)" + ) + parser.add_argument( + "--atoms", nargs="+", help="Atom symbols to print e.g. Fe C H (default: all)" + ) + parser.add_argument( + "--unit", default="angstrom", help="Coordinate unit (default: angstrom)" + ) args = parser.parse_args() if args.name is not None: if args.rdm1a is None: @@ -158,9 +178,9 @@ def main(): # Handle both raw xyz (no header) and standard xyz (2-line header) try: int(lines[0].strip()) - mol.atom = ''.join(lines[2:]) # standard xyz with natom + comment lines + mol.atom = "".join(lines[2:]) # standard xyz with natom + comment lines except ValueError: - mol.atom = ''.join(lines) # raw xyz with no header + mol.atom = "".join(lines) # raw xyz with no header mol.basis = args.basis mol.charge = args.charge mol.spin = args.spin @@ -177,15 +197,23 @@ def main(): spin_density = rdm1a - rdm1b print(f"\nRDM validation:", flush=True) - S = mol.intor('int1e_ovlp') - print(f" Trace rdm1a: {np.trace(rdm1a @ S):.4f} (expected {mol.nelec[0]})", flush=True) - print(f" Trace rdm1b: {np.trace(rdm1b @ S):.4f} (expected {mol.nelec[1]})", flush=True) - print(f" Net spin: {np.trace(spin_density @ S):.4f} (expected {mol.spin})", flush=True) + S = mol.intor("int1e_ovlp") + print( + f" Trace rdm1a: {np.trace(rdm1a @ S):.4f} (expected {mol.nelec[0]})", + flush=True, + ) + print( + f" Trace rdm1b: {np.trace(rdm1b @ S):.4f} (expected {mol.nelec[1]})", + flush=True, + ) + print( + f" Net spin: {np.trace(spin_density @ S):.4f} (expected {mol.spin})", + flush=True, + ) # Filter atoms by symbol if requested if args.atoms: - atom_indices = [i for i in range(mol.natm) - if mol.atom_symbol(i) in args.atoms] + atom_indices = [i for i in range(mol.natm) if mol.atom_symbol(i) in args.atoms] print(f"\nComputing HFCCs for atoms: {args.atoms}", flush=True) else: atom_indices = None @@ -194,5 +222,6 @@ def main(): results = compute_hfcc(mol, rdm1a, rdm1b, atoms=atom_indices) print_hfcc_table(results, title=f"BE-UCCSD Hyperfine Coupling ({args.basis})") + if __name__ == "__main__": main() From 6677d2b38c54ef5091235008d3458cec93307bcb Mon Sep 17 00:00:00 2001 From: Lea Northcote Date: Fri, 14 Aug 2026 16:35:22 -0400 Subject: [PATCH 20/31] Fix ruff check errors: line length, local import, import sort, f-strings --- src/quemb/molbe/ube.py | 7 +++---- src/quemb/scripts/ube_hfcc.py | 5 +++-- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/src/quemb/molbe/ube.py b/src/quemb/molbe/ube.py index 2ce9acc2ca..c9dfe2e19a 100644 --- a/src/quemb/molbe/ube.py +++ b/src/quemb/molbe/ube.py @@ -18,7 +18,7 @@ from warnings import warn import h5py -from numpy import array, einsum, zeros_like +from numpy import array, einsum, zeros, zeros_like from numpy.linalg import multi_dot from pyscf import ao2mo from pyscf.scf.uhf import UHF @@ -88,7 +88,8 @@ def __init__( self.use_df = use_df if use_df: assert hasattr(mf, "with_df") and mf.with_df is not None, ( - "use_df=True requires a density-fitted mf: construct as scf.UHF(mol).density_fit()" + "use_df=True requires a density-fitted mf: " + "construct as scf.UHF(mol).density_fit()" ) self.fobj = fobj @@ -426,8 +427,6 @@ def urdm1_fullbasis(self): Alpha and beta 1-RDMs in the AO basis. Spin density = rdm1a_AO - rdm1b_AO. """ - from numpy import zeros - nao = self.S.shape[0] rdm1a_AO = zeros((nao, nao)) rdm1b_AO = zeros((nao, nao)) diff --git a/src/quemb/scripts/ube_hfcc.py b/src/quemb/scripts/ube_hfcc.py index 069d25aafd..9b11b7ff84 100644 --- a/src/quemb/scripts/ube_hfcc.py +++ b/src/quemb/scripts/ube_hfcc.py @@ -22,6 +22,7 @@ """ import argparse + import numpy as np from pyscf import gto from pyscf.data import nist @@ -196,7 +197,7 @@ def main(): rdm1b = np.load(args.rdm1b) spin_density = rdm1a - rdm1b - print(f"\nRDM validation:", flush=True) + print("\nRDM validation:", flush=True) S = mol.intor("int1e_ovlp") print( f" Trace rdm1a: {np.trace(rdm1a @ S):.4f} (expected {mol.nelec[0]})", @@ -217,7 +218,7 @@ def main(): print(f"\nComputing HFCCs for atoms: {args.atoms}", flush=True) else: atom_indices = None - print(f"\nComputing HFCCs for all atoms", flush=True) + print("\nComputing HFCCs for all atoms", flush=True) results = compute_hfcc(mol, rdm1a, rdm1b, atoms=atom_indices) print_hfcc_table(results, title=f"BE-UCCSD Hyperfine Coupling ({args.basis})") From 977967ea1a90e4a189d77b54fc53d5185c6e23dd Mon Sep 17 00:00:00 2001 From: Lea Northcote Date: Fri, 14 Aug 2026 19:33:20 -0400 Subject: [PATCH 21/31] Fix mypy S_ type assignment errors inherited from main (#191) --- src/quemb/molbe/mbe.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/quemb/molbe/mbe.py b/src/quemb/molbe/mbe.py index b4e8330acb..224fb62a01 100644 --- a/src/quemb/molbe/mbe.py +++ b/src/quemb/molbe/mbe.py @@ -1422,7 +1422,7 @@ def localize( Cpop = multi_dot((C_.T, self.S, C_)) no_core_idx = where(diag(Cpop) > 0.7)[0] C_ = C_[:, no_core_idx] - S_ = multi_dot((C_.T, self.S, C_)) + S_ = multi_dot((C_.T, self.S, C_)) # type: ignore[assignment] es_, vs_ = eigh(S_) s_ = sqrt(es_) s_ = diag(1.0 / s_) @@ -1459,7 +1459,7 @@ def localize( Cpop = diag(Cpop) no_core_idx = where(Cpop > 0.55)[0] C_ = C_[:, no_core_idx] - S_ = multi_dot((C_.T, self.S, C_)) + S_ = multi_dot((C_.T, self.S, C_)) # type: ignore[assignment] es_, vs_ = eigh(S_) s_ = sqrt(es_) s_ = diag(1.0 / s_) From 5aca29b65f72ec02e271e0738b47f0a382c5b9b8 Mon Sep 17 00:00:00 2001 From: Lea Northcote Date: Thu, 20 Aug 2026 10:35:40 -0400 Subject: [PATCH 22/31] Use fobj.mo_coeffs consistently with BE instead of redundant mo_coeff_uccsd --- src/quemb/molbe/be_parallel.py | 4 ++-- src/quemb/molbe/solver.py | 2 -- src/quemb/molbe/ube.py | 6 +----- 3 files changed, 3 insertions(+), 9 deletions(-) diff --git a/src/quemb/molbe/be_parallel.py b/src/quemb/molbe/be_parallel.py index d2c911c99f..0cc87f3204 100644 --- a/src/quemb/molbe/be_parallel.py +++ b/src/quemb/molbe/be_parallel.py @@ -629,8 +629,8 @@ def be_func_parallel_u( e_f, rdm1_tmp, mo_a, mo_b = results_list[i] fobj_a.rdm1__ = rdm1_tmp[0].copy() fobj_b.rdm1__ = rdm1_tmp[1].copy() - fobj_a.mo_coeff_uccsd = mo_a - fobj_b.mo_coeff_uccsd = mo_b + fobj_a.mo_coeffs = mo_a + fobj_b.mo_coeffs = mo_b # Compute and return fragment energy e_1 = 0.0 e_2 = 0.0 diff --git a/src/quemb/molbe/solver.py b/src/quemb/molbe/solver.py index 332f797f87..f41db8ed4a 100644 --- a/src/quemb/molbe/solver.py +++ b/src/quemb/molbe/solver.py @@ -630,8 +630,6 @@ def be_func_u( raise ValueError("Solver not implemented") assert fobj_a._mf is not None and fobj_b._mf is not None - fobj_a.mo_coeff_uccsd = fobj_a._mf.mo_coeff.copy() - fobj_b.mo_coeff_uccsd = fobj_b._mf.mo_coeff.copy() fobj_a.rdm1__ = rdm1_tmp[0].copy() fobj_a._rdm1 = ( multi_dot((fobj_a._mf.mo_coeff, rdm1_tmp[0], fobj_a._mf.mo_coeff.T)) * 0.5 diff --git a/src/quemb/molbe/ube.py b/src/quemb/molbe/ube.py index c9dfe2e19a..0020de5328 100644 --- a/src/quemb/molbe/ube.py +++ b/src/quemb/molbe/ube.py @@ -432,11 +432,7 @@ def urdm1_fullbasis(self): rdm1b_AO = zeros((nao, nao)) def get_mo(fobj): - if hasattr(fobj, "mo_coeff_uccsd"): - return fobj.mo_coeff_uccsd - if fobj._mf is not None: - return fobj._mf.mo_coeff - return fobj._mo_coeffs + return fobj.mo_coeffs for fobj_a, fobj_b in zip(self.Fobjs_a, self.Fobjs_b): # Fragment AO centers - same for alpha and beta From fe78296fbd8a0490b6b00a433580763477b4ecb6 Mon Sep 17 00:00:00 2001 From: Lea Northcote Date: Thu, 20 Aug 2026 11:36:54 -0400 Subject: [PATCH 23/31] Assert custom PySCF ERI support before opposite-spin transform can fail late --- src/quemb/molbe/ube.py | 17 +++++++++++++++++ 1 file changed, 17 insertions(+) diff --git a/src/quemb/molbe/ube.py b/src/quemb/molbe/ube.py index 0020de5328..738bbd58be 100644 --- a/src/quemb/molbe/ube.py +++ b/src/quemb/molbe/ube.py @@ -21,6 +21,7 @@ from numpy import array, einsum, zeros, zeros_like from numpy.linalg import multi_dot from pyscf import ao2mo +from pyscf import lib as pyscf_lib from pyscf.scf.uhf import UHF from quemb.molbe.be_parallel import be_func_parallel_u @@ -91,6 +92,11 @@ def __init__( "use_df=True requires a density-fitted mf: " "construct as scf.UHF(mol).density_fit()" ) + if not equal_bath: + assert _opposite_spin_eri_supported(), ( + "equal_bath=False requires a custom-compiled PySCF " + "(see unrestricted_utils._convert_eri_gen for the patch)." + ) self.fobj = fobj @@ -251,6 +257,12 @@ def initialize(self, eri_, compute_hf): ) assert fobj_a.TA is not None and fobj_b.TA is not None + if fobj_a.TA.shape[1] != fobj_b.TA.shape[1]: + assert _opposite_spin_eri_supported(), ( + "alpha/beta bath sizes differ despite equal_bath=" + f"{self.equal_bath} " + "(see unrestricted_utils._convert_eri_gen for the patch)." + ) if self.use_df: eri_a = self.mf.with_df.ao2mo(fobj_a.TA, compact=True) eri_b = self.mf.with_df.ao2mo(fobj_b.TA, compact=True) @@ -463,6 +475,11 @@ def get_mo(fobj): return rdm1a_AO, rdm1b_AO +def _opposite_spin_eri_supported() -> bool: + libao2mo = pyscf_lib.load_library("libao2mo") + return hasattr(libao2mo, "AO2MOrestore_nr4to1_gen") + + def initialize_pot(n_frag, relAO_per_edge): pot_ = [] From 917b12954983d986a1f1a3748dd906c1da57eaec Mon Sep 17 00:00:00 2001 From: Lea Northcote Date: Thu, 20 Aug 2026 11:44:54 -0400 Subject: [PATCH 24/31] Replace use_df bool with int_transform from UBE --- src/quemb/molbe/ube.py | 31 ++++++++++++++++++------------- 1 file changed, 18 insertions(+), 13 deletions(-) diff --git a/src/quemb/molbe/ube.py b/src/quemb/molbe/ube.py index 738bbd58be..6aaa3b5070 100644 --- a/src/quemb/molbe/ube.py +++ b/src/quemb/molbe/ube.py @@ -27,7 +27,7 @@ from quemb.molbe.be_parallel import be_func_parallel_u from quemb.molbe.fragment import FragPart from quemb.molbe.lo import LocMethods -from quemb.molbe.mbe import BE +from quemb.molbe.mbe import BE, IntTransforms from quemb.molbe.pfrag import Frags from quemb.molbe.solver import be_func_u from quemb.shared.helper import unused @@ -47,7 +47,7 @@ def __init__( compute_hf: bool = True, thr_bath: float = 1.0e-10, equal_bath: bool = True, - use_df: bool = False, + int_transform: IntTransforms = "in-core", ) -> None: """Initialize Unrestricted BE Object (ube🍠) @@ -81,15 +81,24 @@ def __init__( Whether to use a bath with the same number of alpha and beta orbitals. Using equal_bath = False will require custom compiled functions in PySCF to perform integral transformations. Default is True + int_transform : + The integral transformation strategy. UBE currently supports + "in-core" and "out-core-DF" (see :class:`quemb.molbe.mbe.IntTransforms` + for the full set of options restricted BE supports). """ self.unrestricted = True self.thr_bath = thr_bath self.equal_bath = equal_bath - self.use_df = use_df - if use_df: + if int_transform not in ("in-core", "out-core-DF"): + raise NotImplementedError( + f"UBE currently only supports int_transform in " + f"('in-core', 'out-core-DF'), got {int_transform!r}" + ) + self.int_transform = int_transform + if int_transform == "out-core-DF": assert hasattr(mf, "with_df") and mf.with_df is not None, ( - "use_df=True requires a density-fitted mf: " + "int_transform='out-core-DF' requires a density-fitted mf: " "construct as scf.UHF(mol).density_fit()" ) if not equal_bath: @@ -188,7 +197,9 @@ def __init__( self.scratch_dir = scratch_dir self.eri_file = self.scratch_dir / eri_file - self.initialize(None if self.use_df else mf._eri, compute_hf) + self.initialize( + None if self.int_transform == "out-core-DF" else mf._eri, compute_hf + ) def initialize(self, eri_, compute_hf): if compute_hf: @@ -257,13 +268,7 @@ def initialize(self, eri_, compute_hf): ) assert fobj_a.TA is not None and fobj_b.TA is not None - if fobj_a.TA.shape[1] != fobj_b.TA.shape[1]: - assert _opposite_spin_eri_supported(), ( - "alpha/beta bath sizes differ despite equal_bath=" - f"{self.equal_bath} " - "(see unrestricted_utils._convert_eri_gen for the patch)." - ) - if self.use_df: + if self.int_transform == "out-core-DF": eri_a = self.mf.with_df.ao2mo(fobj_a.TA, compact=True) eri_b = self.mf.with_df.ao2mo(fobj_b.TA, compact=True) eri_ab = self.mf.with_df.ao2mo( From 45be38c52553e3dec3fdffc0538a2041c22f0878 Mon Sep 17 00:00:00 2001 From: Lea Northcote Date: Fri, 21 Aug 2026 16:19:57 -0400 Subject: [PATCH 25/31] Fix equal_bath assert checks, activate new test set (Test B) and add the old baseline tests (Test A) --- src/quemb/molbe/ube.py | 19 +- tests/ube-oneshot_test.py | 366 ++++++++++++++++++++++++++++++++++++-- 2 files changed, 367 insertions(+), 18 deletions(-) diff --git a/src/quemb/molbe/ube.py b/src/quemb/molbe/ube.py index 6aaa3b5070..8ef988f9fc 100644 --- a/src/quemb/molbe/ube.py +++ b/src/quemb/molbe/ube.py @@ -79,8 +79,12 @@ def __init__( Threshold for bath orbitals in Schmidt decomposition equal_bath : Whether to use a bath with the same number of alpha and beta orbitals. - Using equal_bath = False will require custom compiled functions in - PySCF to perform integral transformations. Default is True + With equal_bath = False, alpha and beta bath sizes are only forced + to match by a custom-compiled PySCF (see + :func:`quemb.shared.external.unrestricted_utils._convert_eri_gen`) + if they naturally come out unequal; this is checked in + :meth:`initialize` once the real sizes are known, not here. + Default is True int_transform : The integral transformation strategy. UBE currently supports "in-core" and "out-core-DF" (see :class:`quemb.molbe.mbe.IntTransforms` @@ -101,11 +105,6 @@ def __init__( "int_transform='out-core-DF' requires a density-fitted mf: " "construct as scf.UHF(mol).density_fit()" ) - if not equal_bath: - assert _opposite_spin_eri_supported(), ( - "equal_bath=False requires a custom-compiled PySCF " - "(see unrestricted_utils._convert_eri_gen for the patch)." - ) self.fobj = fobj @@ -268,6 +267,12 @@ def initialize(self, eri_, compute_hf): ) assert fobj_a.TA is not None and fobj_b.TA is not None + if fobj_a.TA.shape[1] != fobj_b.TA.shape[1]: + assert _opposite_spin_eri_supported(), ( + "alpha/beta bath sizes differ despite equal_bath=" + f"{self.equal_bath} " + "(see unrestricted_utils._convert_eri_gen for the patch)." + ) if self.int_transform == "out-core-DF": eri_a = self.mf.with_df.ao2mo(fobj_a.TA, compact=True) eri_b = self.mf.with_df.ao2mo(fobj_b.TA, compact=True) diff --git a/tests/ube-oneshot_test.py b/tests/ube-oneshot_test.py index 507c9592ba..b455ac1cf7 100644 --- a/tests/ube-oneshot_test.py +++ b/tests/ube-oneshot_test.py @@ -18,9 +18,293 @@ from pyscf import gto, scf from quemb.molbe import UBE, fragmentate +from quemb.molbe.ube import _opposite_spin_eri_supported class TestOneShot_Unrestricted(unittest.TestCase): + # Test A + + @unittest.skipUnless( + os.getenv("QUEMB_DO_KNOWN_TO_FAIL_TESTS") == "true", + "This test is known to fail.", + ) + def test_hexene_anion_sto3g_frz_ben_autogen(self): + # Hexene anion with frozen core, STO-3G + mol = gto.M() + mol.atom = os.path.join(os.path.dirname(__file__), "xyz/hexene.xyz") + mol.basis = "sto-3g" + mol.charge = -1 + mol.spin = 1 + mol.build() + self.molecular_unrestricted_oneshot_test( + mol, + 1, + "Hexene Anion Frz (BE1)", + True, + -0.35753374, + delta=1e-4, + frag_type="autogen", + equal_bath=False, + ) + + @unittest.skipUnless( + os.getenv("QUEMB_DO_KNOWN_TO_FAIL_TESTS") == "true" + and _opposite_spin_eri_supported(), + "This test is known to fail and/or requires custom-compiled PySCF.", + ) + def test_hexene_anion_sto3g_frz_ben_autogen_be2(self): + mol = gto.M() + mol.atom = os.path.join(os.path.dirname(__file__), "xyz/hexene.xyz") + mol.basis = "sto-3g" + mol.charge = -1 + mol.spin = 1 + mol.build() + self.molecular_unrestricted_oneshot_test( + mol, + 2, + "Hexene Anion Frz (BE2)", + True, + -0.34725961, + delta=1e-4, + frag_type="autogen", + equal_bath=False, + ) + + @unittest.skipUnless( + os.getenv("QUEMB_DO_KNOWN_TO_FAIL_TESTS") == "true" + and os.getenv("QUEMB_DO_EXPENSIVE_TESTS") == "true" + and _opposite_spin_eri_supported(), + "This test is known to fail, expensive, and/or requires custom-compiled PySCF.", + ) + def test_hexene_anion_sto3g_frz_ben_autogen_be3(self): + mol = gto.M() + mol.atom = os.path.join(os.path.dirname(__file__), "xyz/hexene.xyz") + mol.basis = "sto-3g" + mol.charge = -1 + mol.spin = 1 + mol.build() + self.molecular_unrestricted_oneshot_test( + mol, + 3, + "Hexene Anion Frz (BE3)", + True, + -0.34300834, + delta=1e-4, + frag_type="autogen", + equal_bath=False, + ) + + @unittest.skipUnless( + os.getenv("QUEMB_DO_KNOWN_TO_FAIL_TESTS") == "true", + "This test is known to fail.", + ) + def test_hexene_cation_sto3g_frz_ben_autogen(self): + # Hexene cation with frozen core, STO-3G + mol = gto.M() + mol.atom = os.path.join(os.path.dirname(__file__), "xyz/hexene.xyz") + mol.basis = "sto-3g" + mol.charge = 1 + mol.spin = 1 + mol.build() + self.molecular_unrestricted_oneshot_test( + mol, + 1, + "Hexene Cation Frz (BE1)", + True, + -0.40383508, + delta=1e-4, + frag_type="autogen", + equal_bath=False, + ) + + @unittest.skipUnless( + os.getenv("QUEMB_DO_KNOWN_TO_FAIL_TESTS") == "true" + and _opposite_spin_eri_supported(), + "This test is known to fail and/or requires custom-compiled PySCF.", + ) + def test_hexene_cation_sto3g_frz_ben_autogen_be2(self): + mol = gto.M() + mol.atom = os.path.join(os.path.dirname(__file__), "xyz/hexene.xyz") + mol.basis = "sto-3g" + mol.charge = 1 + mol.spin = 1 + mol.build() + self.molecular_unrestricted_oneshot_test( + mol, + 2, + "Hexene Cation Frz (BE2)", + True, + -0.36496690, + delta=1e-4, + frag_type="autogen", + equal_bath=False, + ) + + @unittest.skipUnless( + os.getenv("QUEMB_DO_KNOWN_TO_FAIL_TESTS") == "true" + and os.getenv("QUEMB_DO_EXPENSIVE_TESTS") == "true" + and _opposite_spin_eri_supported(), + "This test is known to fail, expensive, and/or requires custom-compiled PySCF.", + ) + def test_hexene_cation_sto3g_frz_ben_autogen_be3(self): + mol = gto.M() + mol.atom = os.path.join(os.path.dirname(__file__), "xyz/hexene.xyz") + mol.basis = "sto-3g" + mol.charge = 1 + mol.spin = 1 + mol.build() + self.molecular_unrestricted_oneshot_test( + mol, + 3, + "Hexene Cation Frz (BE3)", + True, + -0.36996484, + delta=1e-4, + frag_type="autogen", + equal_bath=False, + ) + + @unittest.skipUnless( + os.getenv("QUEMB_DO_KNOWN_TO_FAIL_TESTS") == "true", + "This test is known to fail.", + ) + def test_hexene_anion_sto3g_unfrz_ben_autogen(self): + # Hexene anion without frozen core, STO-3G + mol = gto.M() + mol.atom = os.path.join(os.path.dirname(__file__), "xyz/hexene.xyz") + mol.basis = "sto-3g" + mol.charge = -1 + mol.spin = 1 + mol.build() + self.molecular_unrestricted_oneshot_test( + mol, + 1, + "Hexene Anion Unfrz (BE1)", + False, + -0.38478279, + delta=1e-4, + frag_type="autogen", + equal_bath=False, + ) + + @unittest.skipUnless( + os.getenv("QUEMB_DO_KNOWN_TO_FAIL_TESTS") == "true" + and _opposite_spin_eri_supported(), + "This test is known to fail and/or requires custom-compiled PySCF.", + ) + def test_hexene_anion_sto3g_unfrz_ben_autogen_be2(self): + mol = gto.M() + mol.atom = os.path.join(os.path.dirname(__file__), "xyz/hexene.xyz") + mol.basis = "sto-3g" + mol.charge = -1 + mol.spin = 1 + mol.build() + self.molecular_unrestricted_oneshot_test( + mol, + 2, + "Hexene Anion Unfrz (BE2)", + False, + -0.39053689, + delta=1e-4, + frag_type="autogen", + equal_bath=False, + ) + + @unittest.skipUnless( + os.getenv("QUEMB_DO_KNOWN_TO_FAIL_TESTS") == "true" + and os.getenv("QUEMB_DO_EXPENSIVE_TESTS") == "true" + and _opposite_spin_eri_supported(), + "This test is known to fail, expensive, and/or requires custom-compiled PySCF.", + ) + def test_hexene_anion_sto3g_unfrz_ben_autogen_be3(self): + mol = gto.M() + mol.atom = os.path.join(os.path.dirname(__file__), "xyz/hexene.xyz") + mol.basis = "sto-3g" + mol.charge = -1 + mol.spin = 1 + mol.build() + self.molecular_unrestricted_oneshot_test( + mol, + 3, + "Hexene Anion Unfrz (BE3)", + False, + -0.38960174, + delta=1e-4, + frag_type="autogen", + equal_bath=False, + ) + + @unittest.skipUnless( + os.getenv("QUEMB_DO_KNOWN_TO_FAIL_TESTS") == "true", + "This test is known to fail.", + ) + def test_hexene_cation_sto3g_unfrz_ben_autogen(self): + mol = gto.M() + mol.atom = os.path.join(os.path.dirname(__file__), "xyz/hexene.xyz") + mol.basis = "sto-3g" + mol.charge = 1 + mol.spin = 1 + mol.build() + self.molecular_unrestricted_oneshot_test( + mol, + 1, + "Hexene Cation Frz (BE1)", + False, + -0.39471433, + delta=1e-4, + frag_type="autogen", + equal_bath=False, + ) + + @unittest.skipUnless( + os.getenv("QUEMB_DO_KNOWN_TO_FAIL_TESTS") == "true" + and _opposite_spin_eri_supported(), + "This test is known to fail and/or requires custom-compiled PySCF.", + ) + def test_hexene_cation_sto3g_unfrz_ben_autogen_be2(self): + mol = gto.M() + mol.atom = os.path.join(os.path.dirname(__file__), "xyz/hexene.xyz") + mol.basis = "sto-3g" + mol.charge = 1 + mol.spin = 1 + mol.build() + self.molecular_unrestricted_oneshot_test( + mol, + 2, + "Hexene Cation Frz (BE2)", + False, + -0.39846777, + delta=1e-4, + frag_type="autogen", + equal_bath=False, + ) + + @unittest.skipUnless( + os.getenv("QUEMB_DO_KNOWN_TO_FAIL_TESTS") == "true" + and os.getenv("QUEMB_DO_EXPENSIVE_TESTS") == "true" + and _opposite_spin_eri_supported(), + "This test is known to fail, expensive, and/or requires custom-compiled PySCF.", + ) + def test_hexene_cation_sto3g_unfrz_ben_autogen_be3(self): + mol = gto.M() + mol.atom = os.path.join(os.path.dirname(__file__), "xyz/hexene.xyz") + mol.basis = "sto-3g" + mol.charge = 1 + mol.spin = 1 + mol.build() + self.molecular_unrestricted_oneshot_test( + mol, + 3, + "Hexene Cation Frz (BE3)", + False, + -0.39729184, + delta=1e-4, + frag_type="autogen", + equal_bath=False, + ) + + # Test B + @unittest.skipUnless( os.getenv("QUEMB_DO_KNOWN_TO_FAIL_TESTS") == "true", "This test is known to fail.", @@ -39,11 +323,23 @@ def test_hexene_anion_sto3g_frz_ben(self): self.molecular_unrestricted_oneshot_test( mol, 2, "Hexene Anion Frz (BE2)", True, -0.34617685, delta=1e-4 ) - """ Cut for expense + + @unittest.skipUnless( + os.getenv("QUEMB_DO_KNOWN_TO_FAIL_TESTS") == "true" + and os.getenv("QUEMB_DO_EXPENSIVE_TESTS") == "true", + "This test is known to fail and/or expensive.", + ) + def test_hexene_anion_sto3g_frz_ben_be3(self): + # Hexene anion with frozen core, STO-3G + mol = gto.M() + mol.atom = os.path.join(os.path.dirname(__file__), "xyz/hexene.xyz") + mol.basis = "sto-3g" + mol.charge = -1 + mol.spin = 1 + mol.build() self.molecular_unrestricted_oneshot_test( mol, 3, "Hexene Anion Frz (BE3)", True, -0.34300832 ) - """ @unittest.skipUnless( os.getenv("QUEMB_DO_KNOWN_TO_FAIL_TESTS") == "true", @@ -63,11 +359,23 @@ def test_hexene_cation_sto3g_frz_ben(self): self.molecular_unrestricted_oneshot_test( mol, 2, "Hexene Cation Frz (BE2)", True, -0.36736494, delta=1e-4 ) - """ Cut for expense + + @unittest.skipUnless( + os.getenv("QUEMB_DO_KNOWN_TO_FAIL_TESTS") == "true" + and os.getenv("QUEMB_DO_EXPENSIVE_TESTS") == "true", + "This test is known to fail and/or expensive.", + ) + def test_hexene_cation_sto3g_frz_ben_be3(self): + # Hexene cation with frozen core, STO-3G + mol = gto.M() + mol.atom = os.path.join(os.path.dirname(__file__), "xyz/hexene.xyz") + mol.basis = "sto-3g" + mol.charge = 1 + mol.spin = 1 + mol.build() self.molecular_unrestricted_oneshot_test( mol, 3, "Hexene Cation Frz (BE3)", True, -0.36996482 ) - """ @unittest.skipUnless( os.getenv("QUEMB_DO_KNOWN_TO_FAIL_TESTS") == "true", @@ -87,11 +395,23 @@ def test_hexene_anion_sto3g_unfrz_ben(self): self.molecular_unrestricted_oneshot_test( mol, 2, "Hexene Anion Unfrz (BE2)", False, -0.39052993, delta=1e-4 ) - """ Cut for expense + + @unittest.skipUnless( + os.getenv("QUEMB_DO_KNOWN_TO_FAIL_TESTS") == "true" + and os.getenv("QUEMB_DO_EXPENSIVE_TESTS") == "true", + "This test is known to fail and/or expensive.", + ) + def test_hexene_anion_sto3g_unfrz_ben_be3(self): + # Hexene anion without frozen core, STO-3G + mol = gto.M() + mol.atom = os.path.join(os.path.dirname(__file__), "xyz/hexene.xyz") + mol.basis = "sto-3g" + mol.charge = -1 + mol.spin = 1 + mol.build() self.molecular_unrestricted_oneshot_test( mol, 3, "Hexene Anion Unfrz (BE3)", False, -0.3895924 ) - """ @unittest.skipUnless( os.getenv("QUEMB_DO_KNOWN_TO_FAIL_TESTS") == "true", @@ -111,19 +431,43 @@ def test_hexene_cation_sto3g_unfrz_ben(self): self.molecular_unrestricted_oneshot_test( mol, 2, "Hexene Cation Frz (BE2)", False, -0.39849056, delta=1e-4 ) - """ Cut for expense + + @unittest.skipUnless( + os.getenv("QUEMB_DO_KNOWN_TO_FAIL_TESTS") == "true" + and os.getenv("QUEMB_DO_EXPENSIVE_TESTS") == "true", + "This test is known to fail and/or expensive.", + ) + def test_hexene_cation_sto3g_unfrz_ben_be3(self): + # Hexene cation without frozen core, STO-3G + mol = gto.M() + mol.atom = os.path.join(os.path.dirname(__file__), "xyz/hexene.xyz") + mol.basis = "sto-3g" + mol.charge = 1 + mol.spin = 1 + mol.build() self.molecular_unrestricted_oneshot_test( mol, 3, "Hexene Cation Frz (BE3)", False, -0.39729215 ) - """ + + # --- Test A: pre-#166 baseline (autogen fragmentation, equal_bath=False) --- + # Kept alongside Test B (chemgen, equal_bath=True) above so both code paths + # retain coverage. Values are the pre-existing ones from `main`, unchanged. def molecular_unrestricted_oneshot_test( - self, mol, n_BE, test_name, frz, exp_result, delta=1e-5 + self, + mol, + n_BE, + test_name, + frz, + exp_result, + delta=1e-5, + frag_type="autogen", + equal_bath=True, ): mf = scf.UHF(mol) mf.kernel() - fobj = fragmentate(frag_type="chemgen", n_BE=n_BE, mol=mol, frozen_core=frz) - mybe = UBE(mf, fobj, equal_bath=True) + fobj = fragmentate(frag_type=frag_type, n_BE=n_BE, mol=mol, frozen_core=frz) + mybe = UBE(mf, fobj, equal_bath=equal_bath) mybe.oneshot(solver="UCCSD", nproc=1) self.assertAlmostEqual( mybe.ebe_tot - mybe.hf_etot, From 6489d08d87eb023dc9c0472141feb6542db413a6 Mon Sep 17 00:00:00 2001 From: Lea Northcote Date: Fri, 21 Aug 2026 16:36:18 -0400 Subject: [PATCH 26/31] Fix Sphinx doc build: use :python: role instead of broken :func: reference to private _convert_eri_gen --- src/quemb/molbe/ube.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/quemb/molbe/ube.py b/src/quemb/molbe/ube.py index 8ef988f9fc..6c0153e7d2 100644 --- a/src/quemb/molbe/ube.py +++ b/src/quemb/molbe/ube.py @@ -81,7 +81,7 @@ def __init__( Whether to use a bath with the same number of alpha and beta orbitals. With equal_bath = False, alpha and beta bath sizes are only forced to match by a custom-compiled PySCF (see - :func:`quemb.shared.external.unrestricted_utils._convert_eri_gen`) + :python:`quemb.shared.external.unrestricted_utils._convert_eri_gen`) if they naturally come out unequal; this is checked in :meth:`initialize` once the real sizes are known, not here. Default is True From 13717d9ddb2e262ccf3c492b56b184c1858daf06 Mon Sep 17 00:00:00 2001 From: Lea Northcote Date: Fri, 21 Aug 2026 22:36:47 -0400 Subject: [PATCH 27/31] Activate Test B for the three passing hexene tests, cation frozen is known to fail --- tests/ube-oneshot_test.py | 16 ---------------- 1 file changed, 16 deletions(-) diff --git a/tests/ube-oneshot_test.py b/tests/ube-oneshot_test.py index b455ac1cf7..cfc83a84a9 100644 --- a/tests/ube-oneshot_test.py +++ b/tests/ube-oneshot_test.py @@ -305,10 +305,6 @@ def test_hexene_cation_sto3g_unfrz_ben_autogen_be3(self): # Test B - @unittest.skipUnless( - os.getenv("QUEMB_DO_KNOWN_TO_FAIL_TESTS") == "true", - "This test is known to fail.", - ) def test_hexene_anion_sto3g_frz_ben(self): # Hexene anion with frozen core, STO-3G mol = gto.M() @@ -377,10 +373,6 @@ def test_hexene_cation_sto3g_frz_ben_be3(self): mol, 3, "Hexene Cation Frz (BE3)", True, -0.36996482 ) - @unittest.skipUnless( - os.getenv("QUEMB_DO_KNOWN_TO_FAIL_TESTS") == "true", - "This test is known to fail.", - ) def test_hexene_anion_sto3g_unfrz_ben(self): # Hexene anion without frozen core, STO-3G mol = gto.M() @@ -413,10 +405,6 @@ def test_hexene_anion_sto3g_unfrz_ben_be3(self): mol, 3, "Hexene Anion Unfrz (BE3)", False, -0.3895924 ) - @unittest.skipUnless( - os.getenv("QUEMB_DO_KNOWN_TO_FAIL_TESTS") == "true", - "This test is known to fail.", - ) def test_hexene_cation_sto3g_unfrz_ben(self): # Hexene cation without frozen core, STO-3G mol = gto.M() @@ -449,10 +437,6 @@ def test_hexene_cation_sto3g_unfrz_ben_be3(self): mol, 3, "Hexene Cation Frz (BE3)", False, -0.39729215 ) - # --- Test A: pre-#166 baseline (autogen fragmentation, equal_bath=False) --- - # Kept alongside Test B (chemgen, equal_bath=True) above so both code paths - # retain coverage. Values are the pre-existing ones from `main`, unchanged. - def molecular_unrestricted_oneshot_test( self, mol, From d2b6447eeba1dd596fe2aeed5f54c4ac12c7cff9 Mon Sep 17 00:00:00 2001 From: Lea Northcote Date: Tue, 25 Aug 2026 08:29:03 -0400 Subject: [PATCH 28/31] Gate anion_frz BE2 as known to fail --- tests/ube-oneshot_test.py | 14 ++++++++------ 1 file changed, 8 insertions(+), 6 deletions(-) diff --git a/tests/ube-oneshot_test.py b/tests/ube-oneshot_test.py index cfc83a84a9..bfcfd6866d 100644 --- a/tests/ube-oneshot_test.py +++ b/tests/ube-oneshot_test.py @@ -305,6 +305,10 @@ def test_hexene_cation_sto3g_unfrz_ben_autogen_be3(self): # Test B + @unittest.skipUnless( + os.getenv("QUEMB_DO_KNOWN_TO_FAIL_TESTS") == "true", + "BE2 bath selection is not deterministic across Python/numpy versions with frozen core", + ) def test_hexene_anion_sto3g_frz_ben(self): # Hexene anion with frozen core, STO-3G mol = gto.M() @@ -389,9 +393,8 @@ def test_hexene_anion_sto3g_unfrz_ben(self): ) @unittest.skipUnless( - os.getenv("QUEMB_DO_KNOWN_TO_FAIL_TESTS") == "true" - and os.getenv("QUEMB_DO_EXPENSIVE_TESTS") == "true", - "This test is known to fail and/or expensive.", + os.getenv("QUEMB_DO_EXPENSIVE_TESTS") == "true", + "Skipped expensive tests for QuEmb.", ) def test_hexene_anion_sto3g_unfrz_ben_be3(self): # Hexene anion without frozen core, STO-3G @@ -421,9 +424,8 @@ def test_hexene_cation_sto3g_unfrz_ben(self): ) @unittest.skipUnless( - os.getenv("QUEMB_DO_KNOWN_TO_FAIL_TESTS") == "true" - and os.getenv("QUEMB_DO_EXPENSIVE_TESTS") == "true", - "This test is known to fail and/or expensive.", + os.getenv("QUEMB_DO_EXPENSIVE_TESTS") == "true", + "Skipped expensive tests for QuEmb.", ) def test_hexene_cation_sto3g_unfrz_ben_be3(self): # Hexene cation without frozen core, STO-3G From 1f6c6ede2c2996dddaf93eec22652958a90b6a22 Mon Sep 17 00:00:00 2001 From: Lea Northcote Date: Tue, 25 Aug 2026 08:40:22 -0400 Subject: [PATCH 29/31] ruff formatting --- tests/ube-oneshot_test.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/tests/ube-oneshot_test.py b/tests/ube-oneshot_test.py index bfcfd6866d..dacda4173a 100644 --- a/tests/ube-oneshot_test.py +++ b/tests/ube-oneshot_test.py @@ -307,7 +307,8 @@ def test_hexene_cation_sto3g_unfrz_ben_autogen_be3(self): @unittest.skipUnless( os.getenv("QUEMB_DO_KNOWN_TO_FAIL_TESTS") == "true", - "BE2 bath selection is not deterministic across Python/numpy versions with frozen core", + "BE2 bath selection is not deterministic across Python/numpy versions " + "with frozen core", ) def test_hexene_anion_sto3g_frz_ben(self): # Hexene anion with frozen core, STO-3G From c6a54eaf5cb95517fe18c3fd5904c4bc71a2371f Mon Sep 17 00:00:00 2001 From: Lea Northcote Date: Fri, 28 Aug 2026 14:07:17 -0400 Subject: [PATCH 30/31] Fix urdm1_fullbasis(): Apply fragment projector in local embedding space before expanding to full AO basis, not after --- src/quemb/molbe/ube.py | 28 +++++++++++++++++----------- 1 file changed, 17 insertions(+), 11 deletions(-) diff --git a/src/quemb/molbe/ube.py b/src/quemb/molbe/ube.py index 6c0153e7d2..092bb75934 100644 --- a/src/quemb/molbe/ube.py +++ b/src/quemb/molbe/ube.py @@ -460,25 +460,31 @@ def get_mo(fobj): # Fragment AO centers - same for alpha and beta cind = [fobj_a.AO_in_frag[i] for i in fobj_a.weight_and_relAO_per_center[1]] - # Democratic partitioning projection in full AO space, built - # per-spin since self.W is split into [Wa, Wb] under frozen_core + # Democratic partitioning projector, built per-spin since self.W + # is split into [Wa, Wb] under frozen_core. Built in the LOCAL + # embedding-space basis (matching the restricted rdm1_fullbasis + # pattern) rather than full AO space -- Proj is not a proper + # idempotent projector once conjugated the other way across the + # non-orthogonal TA embedding transformation. Wa = self.W[0] if self.frozen_core else self.W Wb = self.W[1] if self.frozen_core else self.W - Proj_a = self.S @ Wa[:, cind] @ Wa[:, cind].T @ self.S - Proj_b = self.S @ Wb[:, cind] @ Wb[:, cind].T @ self.S + Pc_a = fobj_a.TA.T @ self.S @ Wa[:, cind] @ Wa[:, cind].T @ self.S @ fobj_a.TA + Pc_b = fobj_b.TA.T @ self.S @ Wb[:, cind] @ Wb[:, cind].T @ self.S @ fobj_b.TA mca = get_mo(fobj_a) mcb = get_mo(fobj_b) - # Transform fragment RDM to full AO space first - rdm1a_full = fobj_a.TA @ mca @ fobj_a.rdm1__ @ mca.T @ fobj_a.TA.T - rdm1b_full = fobj_b.TA @ mcb @ fobj_b.rdm1__ @ mcb.T @ fobj_b.TA.T + # Local density in the embedding-orbital AO-equivalent space + rdm1a_eo = mca @ fobj_a.rdm1__ @ mca.T + rdm1b_eo = mcb @ fobj_b.rdm1__ @ mcb.T - # Apply democratic weight in full AO space - rdm1a_AO += Proj_a @ rdm1a_full - rdm1b_AO += Proj_b @ rdm1b_full + # Project in the SMALL local space, THEN expand to full AO space + rdm1a_center = Pc_a @ rdm1a_eo + rdm1b_center = Pc_b @ rdm1b_eo + + rdm1a_AO += fobj_a.TA @ rdm1a_center @ fobj_a.TA.T + rdm1b_AO += fobj_b.TA @ rdm1b_center @ fobj_b.TA.T - # Symmetrize rdm1a_AO = (rdm1a_AO + rdm1a_AO.T) / 2.0 rdm1b_AO = (rdm1b_AO + rdm1b_AO.T) / 2.0 From 95961df3934ffd5dc30b95c8fb09ce56a854bd98 Mon Sep 17 00:00:00 2001 From: Lea Northcote Date: Fri, 28 Aug 2026 14:41:34 -0400 Subject: [PATCH 31/31] Fix formatting --- src/quemb/molbe/ube.py | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/src/quemb/molbe/ube.py b/src/quemb/molbe/ube.py index 092bb75934..d70a0a2ba1 100644 --- a/src/quemb/molbe/ube.py +++ b/src/quemb/molbe/ube.py @@ -468,8 +468,12 @@ def get_mo(fobj): # non-orthogonal TA embedding transformation. Wa = self.W[0] if self.frozen_core else self.W Wb = self.W[1] if self.frozen_core else self.W - Pc_a = fobj_a.TA.T @ self.S @ Wa[:, cind] @ Wa[:, cind].T @ self.S @ fobj_a.TA - Pc_b = fobj_b.TA.T @ self.S @ Wb[:, cind] @ Wb[:, cind].T @ self.S @ fobj_b.TA + Pc_a = ( + fobj_a.TA.T @ self.S @ Wa[:, cind] @ Wa[:, cind].T @ self.S @ fobj_a.TA + ) + Pc_b = ( + fobj_b.TA.T @ self.S @ Wb[:, cind] @ Wb[:, cind].T @ self.S @ fobj_b.TA + ) mca = get_mo(fobj_a) mcb = get_mo(fobj_b)