From 645d4ea6259d19231c75d0652398deb172f3d64f Mon Sep 17 00:00:00 2001 From: Matthew Hennefarth Date: Thu, 31 Oct 2024 12:23:52 -0500 Subject: [PATCH 1/3] allow hybrid energy decomp and check for SS --- pyscf/mcpdft/mcpdft.py | 34 ++++++++++++++++++++++++++------ pyscf/mcpdft/test/test_mcpdft.py | 22 +++++++++++++++++++++ 2 files changed, 50 insertions(+), 6 deletions(-) diff --git a/pyscf/mcpdft/mcpdft.py b/pyscf/mcpdft/mcpdft.py index ed7327a..0c80695 100644 --- a/pyscf/mcpdft/mcpdft.py +++ b/pyscf/mcpdft/mcpdft.py @@ -206,9 +206,13 @@ def get_energy_decomposition(mc, mo_coeff=None, ci=None, ot=None, otxc=None, E(MC-SCF) = E0 + E1 + E2c + Enc E(MC-PDFT) = E0 + E1 + E2c + EOTx + EOTc - Only compatible with pure translated or fully-translated fnals. If - mc.fcisolver.nroots > 1, lists are returned for everything except - the nuclear potential energy. + For hybrid functionals, + + E(MC-PDFT) = E0 + E1 + E2c + EOTx + EOTc + Enc + + Where the Enc and EOTx/c terms are premultiplied by the hybrid factor. If + mc.fcisolver.nroots > 1, lists are returned for everything except the + nuclear potential energy. Args: mc : an instance of CASSCF or CASCI class @@ -246,6 +250,8 @@ def get_energy_decomposition(mc, mo_coeff=None, ci=None, ot=None, otxc=None, EOTc = correlation part of translated functional e_ncwfn : float or list of length nroots E2ncc = - E0 - E1 - E2c + If hybrid functional, this term is weighted appropriately. For pure + functionals, it is the full NC component ''' if verbose is None: verbose = mc.verbose log = logger.new_logger(mc, verbose) @@ -271,8 +277,8 @@ def get_energy_decomposition(mc, mo_coeff=None, ci=None, ot=None, otxc=None, ) hyb_x, hyb_c = ot._numint.hybrid_coeff(ot.otxc) - if hyb_x > 1e-10 or hyb_c > 1e-10: - raise NotImplementedError("Decomp for hybrid PDFT fnals") + if abs(hyb_x - hyb_c) > 1e-11: + raise NotImplementedError("hybrid functionals with different exchange, correlations components") if not isinstance(ot, transfnal): raise NotImplementedError("Decomp for non-translated PDFT fnals") @@ -309,11 +315,23 @@ def get_energy_decomposition(mc, mo_coeff=None, ci=None, ot=None, otxc=None, def _get_e_decomp(mc, mo_coeff=None, ci=None, ot=None, state=0, verbose=None): ncore = mc.ncore ncas = mc.ncas - if ot is None: ot = mc.otfnal + if ot is None: ot = [mc.otfnal,] if mo_coeff is None: mo_coeff = mc.mo_coeff if ci is None: ci = mc.ci if verbose is None: verbose = mc.verbose + if len(ot) == 1: + hyb_x, hyb_c = ot[0]._numint.hybrid_coeff(ot[0].otxc) + + elif len(ot) == 2: + hyb_x, hyb_c = [fnal._numint.hybrid_coeff(fnal.otxc)[idx] for idx, fnal in enumerate(ot)] + + else: + raise ValueError("ot must be length of 1 or 2") + + if abs(hyb_x - hyb_c) > 1e-11: + raise NotImplementedError("hybrid functionals with different exchange, correlations components") + casdm1s = mc.make_one_casdm1s(ci, state=state) casdm1 = casdm1s[0] + casdm1s[1] casdm2 = mc.make_one_casdm2(ci, state=state) @@ -334,6 +352,10 @@ def _get_e_decomp(mc, mo_coeff=None, ci=None, ot=None, state=0, verbose=None): max_memory=mc.max_memory) for fnal in ot] e_ncwfn = e_mcscf - e_nuc - e_1e - e_coul + + if abs(hyb_x) > 1e-10: + e_ncwfn = hyb_x*e_ncwfn + return e_1e, e_coul, e_otxc, e_ncwfn class _mcscf_env: diff --git a/pyscf/mcpdft/test/test_mcpdft.py b/pyscf/mcpdft/test/test_mcpdft.py index b6b3e4e..40bdad5 100644 --- a/pyscf/mcpdft/test/test_mcpdft.py +++ b/pyscf/mcpdft/test/test_mcpdft.py @@ -256,6 +256,28 @@ def test_decomposition_ss(self): # TODO with self.subTest(objtype=objtype, symmetry=bool(ix), term="sanity"): self.assertAlmostEqual(np.sum(test[:-1]), obj.e_tot, 9) + def test_decomposition_hybrid(self): + ref = [ + 1.0583544218, + -12.5375911135, + 5.8093938665, + -1.6287258807, + -0.0619586538, + -0.5530763650, + ] + terms = ["nuc", "core", "Coulomb", "OT(X)", "OT(C)", "WFN(XC)"] + for ix, mc in enumerate(mcp[0]): + mc_scf = mcpdft.CASSCF(mc, "tPBE0", 5, 2).run() + mc_ci = mcpdft.CASCI(mc, "tPBE0", 5, 2).run() + for obj, objtype in zip((mc_scf, mc_ci), ("CASSCF", "CASCI")): + test = obj.get_energy_decomposition() + for t, r, term in zip(test, ref, terms): + with self.subTest(objtype=objtype, symmetry=bool(ix), term=term): + self.assertAlmostEqual(t, r, delta=1e-5) + with self.subTest(objtype=objtype, symmetry=bool(ix), term="sanity"): + self.assertAlmostEqual(np.sum(test), obj.e_tot, 9) + + def test_decomposition_sa(self): ref_nuc = 1.0583544218 ref_states = np.array( From 7611319d4fcd256cb2ad3ccf8371bdcd77dbb2f7 Mon Sep 17 00:00:00 2001 From: Matthew Hennefarth Date: Thu, 31 Oct 2024 12:53:32 -0500 Subject: [PATCH 2/3] add test for sa hybrid decomp --- pyscf/mcpdft/test/test_mcpdft.py | 98 ++++++++++++++++++++++++++++++++ 1 file changed, 98 insertions(+) diff --git a/pyscf/mcpdft/test/test_mcpdft.py b/pyscf/mcpdft/test/test_mcpdft.py index 40bdad5..c535b54 100644 --- a/pyscf/mcpdft/test/test_mcpdft.py +++ b/pyscf/mcpdft/test/test_mcpdft.py @@ -369,6 +369,104 @@ def test_decomposition_sa(self): self.assertAlmostEqual( np.sum(test[:-1]) + test_nuc, e_ref[state], 9 ) + + def test_decomposition_hybrid_sa(self): + ref_nuc = 1.0583544218 + ref_states = np.array( + [ + [ + -12.5385413915, + 5.8109724796, + -1.6290249990, + -0.0619850920, + -0.5531991067, + ], + [ + -12.1706553996, + 5.5463231972, + -1.6201152933, + -0.0469559736, + -0.5485771470, + ], + [ + -12.1768195314, + 5.5632261670, + -1.6164436229, + -0.0492571730, + -0.5471763843, + ], + [ + -12.1874226655, + 5.5856701424, + -1.6111471613, + -0.0474456546, + -0.5422714244, + ], + [ + -12.1874226655, + 5.5856701424, + -1.6111480360, + -0.0474456745, + -0.5422714244, + ], + ] + ) + terms = ["core", "Coulomb", "OT(X)", "OT(C)", "WFN(XC)"] + for ix, (mc, ms) in enumerate(zip(mcp[1], [0, 1, "mixed"])): + s = bool(ix // 2) + mc_scf = mcpdft.CASSCF(mc, "tPBE0", 5, 2) + if ix == 0: + mc_scf = mc_scf.state_average(mc.weights) + else: + mc_scf = mc_scf.state_average_mix(mc.fcisolver.fcisolvers, mc.weights) + mc_scf.run(ci=mc.ci, mo_coeff=mc.mo_coeff) + objs = [ + mc_scf, + ] + objtypes = [ + "CASSCF", + ] + if ix != 1: # There is no CASCI equivalent to mcp[1][1] + mc_ci = mcpdft.CASCI(mc, "tPBE0", 5, 2) + mc_ci.fcisolver.nroots = ( + 5 - ix + ) # Just check the A1 roots when symmetry is enabled + mc_ci.ci = mc.ci[: 5 - ix] + mc_ci.kernel() + objs.append(mc_ci) + objtypes.append("CASCI") + for obj, objtype in zip(objs, objtypes): + test = obj.get_energy_decomposition() + test_nuc, test_states = test[0], np.array(test[1:]).T + # Arrange states in ascending energy order + e_states = getattr(obj, "e_states", obj.e_tot) + idx = np.argsort(e_states) + test_states = test_states[idx, :] + e_ref = np.array(e_states)[idx] + with self.subTest( + objtype=objtype, symmetry=s, triplet_ms=ms, term="nuc" + ): + self.assertAlmostEqual(test_nuc, ref_nuc, 9) + for state, (test, ref) in enumerate(zip(test_states, ref_states)): + for t, r, term in zip(test, ref, terms): + with self.subTest( + objtype=objtype, + symmetry=s, + triplet_ms=ms, + term=term, + state=state, + ): + self.assertAlmostEqual(t, r, delta=1e-5) + with self.subTest( + objtype=objtype, + symmetry=s, + triplet_ms=ms, + term="sanity", + state=state, + ): + self.assertAlmostEqual( + np.sum(test) + test_nuc, e_ref[state], 9 + ) def test_energy_tot(self): # Test both correctness and energy_tot function purity From 6bd55913d3259b1b24e2cde47e4abd9d532a9a61 Mon Sep 17 00:00:00 2001 From: Matthew Hennefarth Date: Thu, 31 Oct 2024 13:00:26 -0500 Subject: [PATCH 3/3] fix flake --- pyscf/mcpdft/mcpdft.py | 18 +++++++++++------- 1 file changed, 11 insertions(+), 7 deletions(-) diff --git a/pyscf/mcpdft/mcpdft.py b/pyscf/mcpdft/mcpdft.py index 0c80695..c92ccb0 100644 --- a/pyscf/mcpdft/mcpdft.py +++ b/pyscf/mcpdft/mcpdft.py @@ -206,7 +206,7 @@ def get_energy_decomposition(mc, mo_coeff=None, ci=None, ot=None, otxc=None, E(MC-SCF) = E0 + E1 + E2c + Enc E(MC-PDFT) = E0 + E1 + E2c + EOTx + EOTc - For hybrid functionals, + For hybrid functionals, E(MC-PDFT) = E0 + E1 + E2c + EOTx + EOTc + Enc @@ -322,16 +322,20 @@ def _get_e_decomp(mc, mo_coeff=None, ci=None, ot=None, state=0, verbose=None): if len(ot) == 1: hyb_x, hyb_c = ot[0]._numint.hybrid_coeff(ot[0].otxc) - + elif len(ot) == 2: - hyb_x, hyb_c = [fnal._numint.hybrid_coeff(fnal.otxc)[idx] for idx, fnal in enumerate(ot)] + hyb_x, hyb_c = [ + fnal._numint.hybrid_coeff(fnal.otxc)[idx] for idx, fnal in enumerate(ot) + ] else: raise ValueError("ot must be length of 1 or 2") if abs(hyb_x - hyb_c) > 1e-11: - raise NotImplementedError("hybrid functionals with different exchange, correlations components") - + raise NotImplementedError( + "hybrid functionals with different exchange, correlations components" + ) + casdm1s = mc.make_one_casdm1s(ci, state=state) casdm1 = casdm1s[0] + casdm1s[1] casdm2 = mc.make_one_casdm2(ci, state=state) @@ -354,8 +358,8 @@ def _get_e_decomp(mc, mo_coeff=None, ci=None, ot=None, state=0, verbose=None): e_ncwfn = e_mcscf - e_nuc - e_1e - e_coul if abs(hyb_x) > 1e-10: - e_ncwfn = hyb_x*e_ncwfn - + e_ncwfn = hyb_x * e_ncwfn + return e_1e, e_coul, e_otxc, e_ncwfn class _mcscf_env: