Hybrid MC-PDFT Nuclear Gradients

.gitignore

@@ -131,3 +131,6 @@ dmypy.json
 # IDEs
 .idea
 .vscode
+
+# MacOS
+.DS_Store

doc/mcpdft/README.md

@@ -61,11 +61,13 @@ energy calculations for wave functions of various types.
 - Additional properties
   - Decomposition of total electronic energy into core, Coulomb, on-top
   components
-  - Analytical nuclear gradients (non-hybrid functionals only) for:
-    1. Single-state CASSCF wave function: [*JCTC* **2018**, *14*, 126]
-    1. State-averaged CASSCF wave functions: [*JCP* **2020**, *153*, 014106]
-    1. CMS-PDFT: [*Mol Phys* **2022**, 120]
-    1. L-PDFT (no meta-GGA): [*JCTC* **2024**, *20*, 3637]
+  - Analytical nuclear gradients for:
+    - non-hybrid functionals only
+        1. CMS-PDFT: [*Mol Phys* **2022**, 120]
+        1. L-PDFT (no meta-GGA): [*JCTC* **2024**, *20*, 3637]
+    - hybrid and non-hybrid functionals:
+        1. Single-state CASSCF wave function: [*JCTC* **2018**, *14*, 126]
+        1. State-averaged CASSCF wave functions: [*JCP* **2020**, *153*, 014106]
   - Permanent electric dipole moment (non-hybrid functionals only) for:
     1. Single-state CASSCF wave function: [*JCTC* **2021**, *17*, 7586]
     1. State-averaged CASSCF wave functions

pyscf/grad/mcpdft.py

@@ -64,7 +64,6 @@ def gfock_sym(mc, mo_coeff, casdm1, casdm2, h1e, eris):
 
     return gfock
 
-
 def xc_response(ot, vot, rho, Pi, weights, moval_occ, aoval, mo_occ, mo_occup, ncore, nocc, casdm2_pack, ndpi, mo_cas):
     vrho, vPi = vot
 
@@ -164,11 +163,43 @@ def mcpdft_HellmanFeynman_grad (mc, ot, veff1, veff2, mo_coeff=None, ci=None,
 
     casdm1, casdm2 = mc.fcisolver.make_rdm12(ci, ncas, nelecas)
 
+    spin = abs(nelecas[0] - nelecas[1])
+    omega, _, hyb = ot._numint.rsh_and_hybrid_coeff(ot.otxc, spin=spin)
+    if abs(omega) > 1e-11:
+        raise NotImplementedError("range-separated on-top functionals")
+    if abs(hyb[0] - hyb[1]) > 1e-11:
+        raise NotImplementedError(
+            "hybrid on-top functionals with different exchange,correlation components"
+        )
+
+    cas_hyb = hyb[0]
+    ot_hyb = 1.0 - cas_hyb
+
+    if cas_hyb > 1e-11:
+        # TODO: actually implement this in a more efficient manner
+        # That is, lets not call grad_elec, but lets actually do this our self maybe?
+        # Can then use get_pdft_veff with drop_mcwfn = False to automatically get
+        # Generalized fock matrix terms, but then have to deal explicitly with the
+        # derivative of eri terms and auxbasis stuff. Also, get_pdft_veff with
+        # drop_mcwfn=False means the get_wfn_response doesn't have to add the
+        # eris to the veff2 since it should just include it already
+        if auxbasis_response:
+            from pyscf.df.grad import casscf as casscf_grad
+        else:
+            from pyscf.grad import casscf as casscf_grad
+
+        cas_grad = casscf_grad.Gradients(mc)
+
+        de_cas = cas_hyb * cas_grad.grad_elec(
+            mo_coeff=mo_coeff, ci=ci, atmlst=atmlst, verbose=verbose
+        )
+
     # gfock = Generalized Fock, Adv. Chem. Phys., 69, 63
     dm_core = 2 * mo_core @ mo_core.T
     dm_cas = mo_cas @ casdm1 @ mo_cas.T
 
-    gfock = gfock_sym(mc, mo_coeff, casdm1, casdm2, mc.get_hcore() + veff1, veff2)
+    gfock = gfock_sym(mc, mo_coeff, casdm1, casdm2, ot_hyb*mc.get_hcore() + veff1, veff2)
+
     dme0 = mo_coeff @ (0.5*(gfock+gfock.T)) @ mo_coeff.T
     del gfock
 
@@ -190,7 +221,7 @@ def mcpdft_HellmanFeynman_grad (mc, ot, veff1, veff2, mo_coeff=None, ci=None,
     # MRH: vhf1c and vhf1a should be the TRUE vj_c and vj_a (no vk!)
     vj = mf_grad.get_jk (dm=dm1)[0]
     if auxbasis_response:
-        de_aux += np.squeeze (vj.aux[:,:,atmlst,:])
+        de_aux += ot_hyb*np.squeeze (vj.aux[:,:,atmlst,:])
 
     # MRH: Now I have to compute the gradient of the on-top energy
     # This involves derivatives of the orbitals that construct rho and Pi and
@@ -330,6 +361,9 @@ def coul_term(p0, p1):
 
     de_hcore, de_coul, de_xc, de_nuc, de_renorm = sum_terms(mf_grad, mol, atmlst, dm1, dme0, coul_term,
                                                                         dvxc)
+    # Deal with hybridization
+    de_hcore *= ot_hyb
+    de_coul *= ot_hyb
 
     logger.debug (mc, "MC-PDFT Hellmann-Feynman nuclear:\n{}".format (de_nuc))
     logger.debug (mc, "MC-PDFT Hellmann-Feynman hcore component:\n{}".format (
@@ -353,6 +387,10 @@ def coul_term(p0, p1):
         logger.debug (mc, "MC-PDFT Hellmann-Feynman aux component:\n{}".format
             (de_aux))
 
+    if cas_hyb > 1e-11:
+        de += de_cas
+        logger.debug(mc, "MC-PDFT Hellmann-Feynman CAS component:\n{}".format(de_cas))
+
     t1 = logger.timer (mc, 'PDFT HlFn total', *t0)
 
     return de
@@ -382,9 +420,9 @@ def _not_implemented_check (self):
         omega, alpha, hyb = ot._numint.rsh_and_hybrid_coeff (
             otxc, spin=spin)
         hyb_x, hyb_c = hyb
-        if hyb_x or hyb_c:
+        if abs(hyb_x - hyb_c) >1e-11:
             raise NotImplementedError (
-                "{} for hybrid MC-PDFT functionals".format (name)
+                "{} for hybrid MC-PDFT functionals with different exchange, correlation".format (name)
             )
         if omega or alpha:
             raise NotImplementedError (
@@ -400,7 +438,7 @@ def get_wfn_response (self, state=None, verbose=None, mo=None,
         if nlag is None: nlag = self.nlag
         if (veff1 is None) or (veff2 is None):
             veff1, veff2 = self.base.get_pdft_veff (mo, ci[state],
-                incl_coul=True, paaa_only=True)
+                incl_coul=True, paaa_only=True, drop_mcwfn=True)
 
         log = logger.new_logger(self, verbose)
 
@@ -412,6 +450,26 @@ def my_hcore ():
             return self.base.get_hcore () + veff1
         fcasscf.get_hcore = my_hcore
 
+        nelecas = self.base.nelecas
+        spin = abs (nelecas[0]-nelecas[1])
+        omega, alpha, hyb = self.base.otfnal._numint.rsh_and_hybrid_coeff(self.base.otxc, spin=spin)
+        if omega or alpha:
+            raise NotImplementedError("range-separated MC-PDFT functionals")
+        if abs(hyb[0] - hyb[1]) > 1e-11:
+            raise NotImplementedError("hybrid on-top functional with different exchange,correlation components")
+        cas_hyb = hyb[0]
+
+        if cas_hyb > 1e-11:# and len(ci) > 1:
+            # For SS-CASSCF, there are no Lagrange multipliers
+            # This is only needed in SA case
+            eris = self.base.ao2mo(mo_coeff=mo)
+            terms = ["vhf_c", "papa", "ppaa", "j_pc", "k_pc"]
+            for term in terms:
+                setattr(eris, term, getattr(veff2, term) + cas_hyb*getattr(eris, term)[:])
+            veff2.vhf_c
+            veff2 = eris
+
+
         g_all_state = newton_casscf.gen_g_hop (fcasscf, mo, ci[state], veff2, verbose)[0]
 
         g_all = np.zeros (nlag)
@@ -455,9 +513,18 @@ def get_ham_response (self, state=None, atmlst=None, verbose=None, mo=None,
         fcasscf = self.make_fcasscf (state)
         fcasscf.mo_coeff = mo
         fcasscf.ci = ci[state]
-        return mcpdft_HellmanFeynman_grad (fcasscf, self.base.otfnal, veff1,
-            veff2, mo_coeff=mo, ci=ci[state], atmlst=atmlst, mf_grad=mf_grad,
-            verbose=verbose)
+
+        return mcpdft_HellmanFeynman_grad(
+            fcasscf,
+            self.base.otfnal,
+            veff1,
+            veff2,
+            mo_coeff=mo,
+            ci=ci[state],
+            atmlst=atmlst,
+            mf_grad=mf_grad,
+            verbose=verbose,
+        )
 
     def get_init_guess (self, bvec, Adiag, Aop, precond):
         '''Initial guess should solve the problem for SA-SA rotations'''
@@ -527,19 +594,26 @@ def kernel (self, **kwargs):
         '''Cache the effective Hamiltonian terms so you don't have to
         calculate them twice
         '''
-        state = kwargs['state'] if 'state' in kwargs else self.state
+
+        state = kwargs["state"] if "state" in kwargs else self.state
         if state is None:
-            raise NotImplementedError ('Gradient of PDFT state-average energy')
-        self.state = state # Not the best code hygiene maybe
-        mo = kwargs['mo'] if 'mo' in kwargs else self.base.mo_coeff
-        ci = kwargs['ci'] if 'ci' in kwargs else self.base.ci
-        if isinstance (ci, np.ndarray): ci = [ci] # hack hack hack...
-        kwargs['ci'] = ci
-        if ('veff1' not in kwargs) or ('veff2' not in kwargs):
-            kwargs['veff1'], kwargs['veff2'] = self.base.get_pdft_veff (mo,
-                ci, incl_coul=True, paaa_only=True, state=state)
-        if 'mf_grad' not in kwargs:
-            kwargs['mf_grad'] = self.base.get_rhf_base ().nuc_grad_method ()
+            raise NotImplementedError("Gradient of PDFT state-average energy")
+
+        self.state = state  # Not the best code hygiene maybe
+        mo = kwargs["mo"] if "mo" in kwargs else self.base.mo_coeff
+        ci = kwargs["ci"] if "ci" in kwargs else self.base.ci
+        if isinstance(ci, np.ndarray):
+            ci = [ci]  # hack hack hack...
+
+        kwargs["ci"] = ci
+        if ("veff1" not in kwargs) or ("veff2" not in kwargs):
+            kwargs["veff1"], kwargs["veff2"] = self.base.get_pdft_veff(
+                mo, ci, incl_coul=True, paaa_only=True, state=state, drop_mcwfn=True
+            )
+
+        if "mf_grad" not in kwargs:
+            kwargs["mf_grad"] = self.base.get_rhf_base().nuc_grad_method()
+
         return super().kernel (**kwargs)
 
     def project_Aop (self, Aop, ci, state):

pyscf/grad/test/test_grad_metagga_mcpdft.py

@@ -21,7 +21,6 @@
 from pyscf.data.nist import BOHR
 from pyscf import mcpdft
 
-
 def diatomic(
     atom1,
     atom2,
@@ -121,6 +120,45 @@ def test_grad_lih_sa2tm06l22_sto3g(self):
                 de = mc.kernel(state=state)[1, 0] / BOHR
                 self.assertAlmostEqual(de, DE_REF[state], 5)
 
+    def test_grad_lih_ssmc2322_sto3g(self):
+        mc = diatomic("Li", "H", 0.8, "MC23", "STO-3G", 2, 2, 1, grids_level=1)
+        de = mc.kernel()[1, 0] / BOHR
+
+        # Numerical from this software
+        # PySCF commit:         f2c2d3f963916fb64ae77241f1b44f24fa484d96
+        # PySCF-forge commit:   e82ba940654cd0b91f799e889136a316fda34b10
+        DE_REF = -1.0641645070
+
+        self.assertAlmostEqual(de, DE_REF, 5)
+
+    def test_grad_lih_sa2mc2322_sto3g(self):
+        mc = diatomic("Li", "H", 0.8, "MC23", "STO-3G", 2, 2, 2, grids_level=1)
+
+        # Numerical from this software
+        # PySCF commit:         f2c2d3f963916fb64ae77241f1b44f24fa484d96
+        # PySCF-forge commit:   e82ba940654cd0b91f799e889136a316fda34b10
+        DE_REF = [-1.0510225010, -0.8963063432]
+
+        for state in range(2):
+            with self.subTest(state=state):
+                de = mc.kernel(state=state)[1, 0] / BOHR
+                self.assertAlmostEqual(de, DE_REF[state], 5)
+
+    def test_grad_lih_sa2mc2322_sto3g_df(self):
+        mc = diatomic(
+            "Li", "H", 0.8, "MC23", "STO-3G", 2, 2, 2, grids_level=1, density_fit=df
+        )
+
+        # Numerical from this software
+        # PySCF commit:         f2c2d3f963916fb64ae77241f1b44f24fa484d96
+        # PySCF-forge commit:   ee6ac742fbc79d170bc4b63ef2b2c4b49478c53a
+        DE_REF = [-1.0510303416, -0.8963992331]
+
+        for state in range(2):
+            with self.subTest(state=state):
+                de = mc.kernel(state=state)[1, 0] / BOHR
+                self.assertAlmostEqual(de, DE_REF[state], 5)
+
 
 if __name__ == "__main__":
     print("Full Tests for MC-PDFT gradients with meta-GGA functionals")

pyscf/mcpdft/mcpdft.py

@@ -136,15 +136,10 @@ def energy_mcwfn(mc, mo_coeff=None, ci=None, ot=None, state=0, casdm1s=None,
     if casdm1s is None: casdm1s = mc.make_one_casdm1s(ci=ci, state=state)
     if casdm2 is None: casdm2 = mc.make_one_casdm2(ci=ci, state=state)
     log = logger.new_logger(mc, verbose=verbose)
-    nelecas = mc.nelecas
     dm1s = _dms.casdm1s_to_dm1s(mc, casdm1s, mo_coeff=mo_coeff)
     cascm2 = _dms.dm2_cumulant(casdm2, casdm1s)
 
-    spin = abs(nelecas[0] - nelecas[1])
-    omega, alpha, hyb = ot._numint.rsh_and_hybrid_coeff(ot.otxc, spin=spin)
-    if abs(omega) > 1e-11:
-        raise NotImplementedError("range-separated on-top functionals")
-    hyb_x, hyb_c = hyb
+    hyb_x, hyb_c = ot._numint.rsh_and_hybrid_coeff(ot.otxc, mc.mol.spin)[2]
 
     Vnn = mc._scf.energy_nuc()
     h = mc._scf.get_hcore()
@@ -552,7 +547,7 @@ def dump_flags(self, verbose=None):
 
     def get_pdft_veff(self, mo=None, ci=None, state=0, casdm1s=None,
                       casdm2=None, incl_coul=False, paaa_only=False, aaaa_only=False,
-                      jk_pc=False, drop_mcwfn=False, incl_energy=False):
+                      jk_pc=False, drop_mcwfn=False, incl_energy=False, ot=None):
         '''Get the 1- and 2-body MC-PDFT effective potentials for a set
         of mos and ci vectors
 
@@ -588,6 +583,7 @@ def get_pdft_veff(self, mo=None, ci=None, state=0, casdm1s=None,
                 (ie the ``Hartree exchange-correlation'') from the response
             incl_energy : logical
                 If true, includes the on-top potential energy as a 3rd return argument
+            ot : an instance of otfnal class
 
         Returns:
             veff1 : ndarray of shape (nao, nao)
@@ -600,22 +596,51 @@ def get_pdft_veff(self, mo=None, ci=None, state=0, casdm1s=None,
                 On-top energy. Only included if incl_energy is true.
         '''
         t0 = (logger.process_clock(), logger.perf_counter())
-        if mo is None: mo = self.mo_coeff
-        if ci is None: ci = self.ci
-        if casdm1s is None: casdm1s = self.make_one_casdm1s(ci, state=state)
-        if casdm2 is None: casdm2 = self.make_one_casdm2(ci, state=state)
+        if mo is None:
+            mo = self.mo_coeff
+        if ci is None:
+            ci = self.ci
+        if casdm1s is None:
+            casdm1s = self.make_one_casdm1s(ci, state=state)
+        if casdm2 is None:
+            casdm2 = self.make_one_casdm2(ci, state=state)
+        if ot is None:
+            ot = self.otfnal
+
         ncore, ncas = self.ncore, self.ncas
         dm1s = _dms.casdm1s_to_dm1s(self, casdm1s, mo_coeff=mo,
                                     ncore=ncore, ncas=ncas)
         cascm2 = _dms.dm2_cumulant(casdm2, casdm1s)
 
-        E_ot, pdft_veff1, pdft_veff2 = pdft_veff.kernel(self.otfnal, dm1s,
-                                                        cascm2, mo, ncore, ncas, max_memory=self.max_memory,
-                                                        paaa_only=paaa_only, aaaa_only=aaaa_only, jk_pc=jk_pc,
-                                                        drop_mcwfn=drop_mcwfn)
+        spin = abs(self.nelecas[0] - self.nelecas[1])
+        omega, _, hyb = ot._numint.rsh_and_hybrid_coeff(ot.otxc, spin=spin)
+
+        if abs(omega) > 1e-11:
+            raise NotImplementedError("range-separated on-top 1e potentials")
+        if abs(hyb[0] > hyb[1]) > 1e-11:
+            raise NotImplementedError("hybrid functionals with different exchange, correlation components")
+
+        cas_hyb = hyb[0]
+        ot_hyb = 1.0-cas_hyb
+
+        E_ot, pdft_veff1, pdft_veff2 = pdft_veff.kernel(
+            ot,
+            dm1s,
+            cascm2,
+            mo,
+            ncore,
+            ncas,
+            max_memory=self.max_memory,
+            paaa_only=paaa_only,
+            aaaa_only=aaaa_only,
+            jk_pc=jk_pc,
+        )
 
         if incl_coul:
-            pdft_veff1 += self._scf.get_j(self.mol, dm1s[0] + dm1s[1])
+            pdft_veff1 += ot_hyb*self._scf.get_j(self.mol, dm1s[0] + dm1s[1])
+
+        if not drop_mcwfn and cas_hyb > 1e-11:
+            raise NotImplementedError("adding mcwfn response to pdft_veff2")
 
         logger.timer(self, 'get_pdft_veff', *t0)