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/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/molbe/be_parallel.py b/src/quemb/molbe/be_parallel.py index 5532fb8b73..0cc87f3204 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( @@ -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_: @@ -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_coeffs = mo_a + fobj_b.mo_coeffs = 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/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_) diff --git a/src/quemb/molbe/misc.py b/src/quemb/molbe/misc.py index b6be7838fd..9abab1b0be 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,10 @@ def be2puffin( # Run oneshot embedding and return system energy mybe.oneshot(solver=solver, nproc=nproc, ompnum=ompnum) - return mybe.ebe_tot - mybe.ebe_hf + if unrestricted: + return mybe.ebe_tot - mybe.hf_etot + else: + return mybe.ebe_tot - mybe.ebe_hf 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 09d826b9f0..76b23d47be 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,36 @@ 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 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)]) 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/solver.py b/src/quemb/molbe/solver.py index 55ec01a5b4..f41db8ed4a 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 @@ -635,10 +631,9 @@ def be_func_u( assert fobj_a._mf is not None and fobj_b._mf is not None 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 d5f8251da3..d70a0a2ba1 100644 --- a/src/quemb/molbe/ube.py +++ b/src/quemb/molbe/ube.py @@ -18,15 +18,16 @@ 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 import lib as pyscf_lib from pyscf.scf.uhf import UHF 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 @@ -45,6 +46,8 @@ def __init__( pop_method: str | None = None, compute_hf: bool = True, thr_bath: float = 1.0e-10, + equal_bath: bool = True, + int_transform: IntTransforms = "in-core", ) -> None: """Initialize Unrestricted BE Object (ube🍠) @@ -72,11 +75,36 @@ 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. + With equal_bath = False, alpha and beta bath sizes are only forced + to match by a custom-compiled PySCF (see + :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 + 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 + 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, ( + "int_transform='out-core-DF' requires a density-fitted mf: " + "construct as scf.UHF(mol).density_fit()" + ) self.fobj = fobj @@ -109,7 +137,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 @@ -117,8 +144,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 @@ -171,7 +196,9 @@ 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.int_transform == "out-core-DF" else mf._eri, compute_hf + ) def initialize(self, eri_, compute_hf): if compute_hf: @@ -200,52 +227,65 @@ 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_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 + 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)." ) - fobj_b.sd( - self.W, - self.lmo_coeff_b, - self.Nocc[1], - thr_bath=self.thr_bath, + 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( + (fobj_a.TA, fobj_a.TA, fobj_b.TA, fobj_b.TA), compact=True ) - - 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 + 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 + 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) @@ -255,7 +295,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) @@ -280,7 +319,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) @@ -341,8 +379,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") @@ -359,7 +396,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, @@ -368,19 +405,19 @@ 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: 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, nproc=nproc, ompnum=ompnum, - relax_density=False, + relax_density=relax_density, frozen=self.frozen_core, ) unused(E_comp) @@ -391,7 +428,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), @@ -403,6 +440,65 @@ def oneshot(self, solver="UCCSD", nproc=1, ompnum=4): ) ) + def urdm1_fullbasis(self): + """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. + """ + nao = self.S.shape[0] + rdm1a_AO = zeros((nao, nao)) + rdm1b_AO = zeros((nao, nao)) + + def get_mo(fobj): + 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 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 + 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) + + # 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 + + # 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 + + rdm1a_AO = (rdm1a_AO + rdm1a_AO.T) / 2.0 + rdm1b_AO = (rdm1b_AO + rdm1b_AO.T) / 2.0 + + 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_ = [] 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..9b11b7ff84 --- /dev/null +++ b/src/quemb/scripts/ube_hfcc.py @@ -0,0 +1,228 @@ +# 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} {'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("\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("\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() 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 diff --git a/tests/ube-oneshot_test.py b/tests/ube-oneshot_test.py index a86ab120bc..dacda4173a 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 """ @@ -12,15 +18,18 @@ 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(self): - # Linear Equidistant (r=1Å) H8 Chain, STO-3G + 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" @@ -28,21 +37,69 @@ 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.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 + 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 + 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(self): - # Linear Equidistant (r=1Å) H8 Chain, cc-pVDZ + 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" @@ -50,21 +107,69 @@ 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.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 + 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 + 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(self): - # Octane, STO-3G + 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" @@ -72,20 +177,240 @@ def test_hexene_anion_sto3g_unfrz_ben(self): mol.spin = 1 mol.build() self.molecular_unrestricted_oneshot_test( - mol, 1, "Hexene Anion Unfrz (BE1)", False, -0.38478279 + 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", + "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() + 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 + mol, 1, "Hexene Anion Frz (BE1)", True, -0.35753375 ) self.molecular_unrestricted_oneshot_test( - mol, 3, "Hexene Anion Unfrz (BE3)", False, -0.38960174 + mol, 2, "Hexene Anion Frz (BE2)", True, -0.34617685, delta=1e-4 + ) + + @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", "This test is known to fail.", ) + def test_hexene_cation_sto3g_frz_ben(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.40383505 + ) + self.molecular_unrestricted_oneshot_test( + mol, 2, "Hexene Cation Frz (BE2)", True, -0.36736494, delta=1e-4 + ) + + @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 + ) + + def test_hexene_anion_sto3g_unfrz_ben(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 + ) + self.molecular_unrestricted_oneshot_test( + mol, 2, "Hexene Anion Unfrz (BE2)", False, -0.39052993, delta=1e-4 + ) + + @unittest.skipUnless( + 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 + 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 + ) + 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" @@ -96,27 +421,48 @@ 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.39849056, delta=1e-4 ) + + @unittest.skipUnless( + 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 + 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 + 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, + frag_type="autogen", + equal_bath=True, ): mf = scf.UHF(mol) mf.kernel() - fobj = fragmentate(frag_type="autogen", n_BE=n_BE, mol=mol, frozen_core=frz) - mybe = UBE(mf, fobj) + 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.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, )