Merge pull request #76 from matthew-hennefarth/hybrid-decomp

MC-PDFT Hybrid Energy Decomposition

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 = <H> - 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,27 @@ 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 +356,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:

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(
@@ -347,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