Merge pull request #120 from matthew-hennefarth/mcpdft-mgga-nucgrad

PDFT meta-GGA nuclear gradients and MS-PDFT meta-GGA tests

Author: MatthewRHermes

doc/mcpdft/README.md

@@ -65,7 +65,7 @@ energy calculations for wave functions of various types.
     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: [*JCTC* **2024**, *20*, 3637]
+    1. L-PDFT (no meta-GGA): [*JCTC* **2024**, *20*, 3637]
   - 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

@@ -68,11 +68,19 @@ def gfock_sym(mc, mo_coeff, casdm1, casdm2, h1e, eris):
 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
 
+
     # Vpq + Vpqrs * Drs ; I'm not sure why the list comprehension down
     # there doesn't break ao's stride order but I'm not complaining
     vrho = _contract_eff_rho(vPi, rho.sum(0), add_eff_rho=vrho)
-    tmp_dv = np.stack([ot.get_veff_1body(rho, Pi, [ao_i, moval_occ], weights, kern=vrho) for ao_i in aoval], axis=0)
-    tmp_dv = (tmp_dv * mo_occ[None,:,:] * mo_occup[None, None,:nocc]).sum(2)
+
+    tmp_dv = np.stack(
+        [
+            ot.get_veff_1body(rho, Pi, [ao_i, moval_occ], weights, kern=vrho)
+            for ao_i in aoval
+        ],
+        axis=0,
+    )
+    tmp_dv = (tmp_dv * mo_occ[None, :, :] * mo_occup[None, None, :nocc]).sum(2)
 
     # Vpuvx * Lpuvx ; remember the stupid slowest->fastest->medium
     # stride order of the ao grid arrays
@@ -202,7 +210,7 @@ def mcpdft_HellmanFeynman_grad (mc, ot, veff1, veff2, mo_coeff=None, ci=None,
     casdm1, casdm2 = mc.fcisolver.make_rdm12(ci, ncas, nelecas)
     twoCDM = _dms.dm2_cumulant (casdm2, casdm1)
     dm1 = tag_array (dm1, mo_coeff=mo_occ, mo_occ=mo_occup[:nocc])
-    make_rho = ot._numint._gen_rho_evaluator (mol, dm1, 1)[0]
+    make_rho = ot._numint._gen_rho_evaluator (mol, dm1, hermi=1, with_lapl=False)[0]
     dvxc = np.zeros ((3,nao))
     idx = np.array ([[1,4,5,6],[2,5,7,8],[3,6,8,9]], dtype=np.int_)
     # For addressing particular ao derivatives
@@ -215,9 +223,13 @@ def mcpdft_HellmanFeynman_grad (mc, ot, veff1, veff2, mo_coeff=None, ci=None,
     for k, ia in enumerate (atmlst):
         full_atmlst[ia] = k
 
-    ndao = (1, 4)[ot.dens_deriv]
+    # for LDA we need 1 deriv, GGA: 2 deriv
+    # for mGGA with tau, we only need 2 deriv
+    ao_deriv = (1, 2, 2)[ot.dens_deriv]
+    ndao = (1, 4, 10)[ao_deriv]
+    ndrho = (1, 4, 5)[ot.dens_deriv]
     ndpi = (1, 4)[ot.Pi_deriv]
-    ncols = 1.05 * 3 * (ndao * (nao + nocc) + max(ndao * nao, ndpi * ncas * ncas))
+    ncols = 1.05 * 3 * (ndao * (nao + nocc) + max(ndrho * nao, ndpi * ncas * ncas))
 
     for ia, (coords, w0, w1) in enumerate (rks_grad.grids_response_cc (
             ot.grids)):
@@ -235,27 +247,39 @@ def mcpdft_HellmanFeynman_grad (mc, ot, veff1, veff2, mo_coeff=None, ci=None,
         blksize = max (BLKSIZE, min (blksize, ngrids, BLKSIZE*1200))
         t1 = logger.timer (mc, 'PDFT HlFn quadrature atom {} mask and memory '
             'setup'.format (ia), *t1)
-        for ip0 in range (0, ngrids, blksize):
-            ip1 = min (ngrids, ip0+blksize)
-            mask = gen_grid.make_mask (mol, coords[ip0:ip1])
-            logger.info (mc, ('PDFT gradient atom {} slice {}-{} of {} '
-                'total').format (ia, ip0, ip1, ngrids))
-            ao = ot._numint.eval_ao (mol, coords[ip0:ip1],
-                deriv=ot.dens_deriv+1, non0tab=mask)
+
+        for ip0 in range(0, ngrids, blksize):
+            ip1 = min(ngrids, ip0 + blksize)
+            mask = gen_grid.make_mask(mol, coords[ip0:ip1])
+            logger.info(
+                mc,
+                ("PDFT gradient atom {} slice {}-{} of {} total").format(
+                    ia, ip0, ip1, ngrids
+                ),
+            )
+
+            ao = ot._numint.eval_ao(mol, coords[ip0:ip1], deriv=ao_deriv, non0tab=mask)
+
             # Need 1st derivs for LDA, 2nd for GGA, etc.
-            t1 = logger.timer (mc, ('PDFT HlFn quadrature atom {} ao '
-                'grids').format (ia), *t1)
+            t1 = logger.timer(
+                mc, ("PDFT HlFn quadrature atom {} ao " "grids").format(ia), *t1
+            )
             # Slice down ao so as not to confuse the rho and Pi generators
-            if ot.xctype == 'LDA':
+            if ot.xctype == "LDA":
                 aoval = ao[0]
-            if ot.xctype == 'GGA':
+            elif ot.xctype == "GGA":
                 aoval = ao[:4]
+            elif ot.xctype == "MGGA":
+                aoval = ao[:4]
+            else:
+                raise ValueError("Unknown xctype: {}".format(ot.xctype))
+
             rho = make_rho (0, aoval, mask, ot.xctype) / 2.0
             rho = np.stack ((rho,)*2, axis=0)
             t1 = logger.timer (mc, ('PDFT HlFn quadrature atom {} rho '
                 'calc').format (ia), *t1)
             Pi = get_ontop_pair_density (ot, rho, aoval, twoCDM, mo_cas,
-                ot.dens_deriv, mask)
+                ot.Pi_deriv, mask)
             t1 = logger.timer (mc, ('PDFT HlFn quadrature atom {} Pi '
                 'calc').format (ia), *t1)
 

pyscf/grad/test/test_diatomic_gradients.py

@@ -13,24 +13,24 @@
 # See the License for the specific language governing permissions and
 # limitations under the License.
 #
-import numpy as np
-from scipy import linalg
-from pyscf import gto, scf, df, dft, mcscf, lib
+from pyscf import gto, scf, df, dft
 from pyscf.data.nist import BOHR
 from pyscf import mcpdft
-#from pyscf.fci import csf_solver
-from pyscf.grad.cmspdft import diab_response, diab_grad, diab_response_o0, diab_grad_o0
-from pyscf.grad import mspdft as mspdft_grad
-import unittest, math
+import unittest
 
 def diatomic (atom1, atom2, r, fnal, basis, ncas, nelecas, nstates,
               charge=None, spin=None, symmetry=False, cas_irrep=None,
-              density_fit=False):
+              density_fit=False, grids_level=9):
+    global mols
     xyz = '{:s} 0.0 0.0 0.0; {:s} {:.3f} 0.0 0.0'.format (atom1, atom2, r)
     mol = gto.M (atom=xyz, basis=basis, charge=charge, spin=spin, symmetry=symmetry, verbose=0, output='/dev/null')
+    mols.append(mol)
     mf = scf.RHF (mol)
-    if density_fit: mf = mf.density_fit (auxbasis = df.aug_etb (mol))
-    mc = mcpdft.CASSCF (mf.run (), fnal, ncas, nelecas, grids_level=9)
+
+    if density_fit: 
+        mf = mf.density_fit (auxbasis = df.aug_etb (mol))
+
+    mc = mcpdft.CASSCF (mf.run (), fnal, ncas, nelecas, grids_level=grids_level)
     #if spin is not None: smult = spin+1
     #else: smult = (mol.nelectron % 2) + 1
     #mc.fcisolver = csf_solver (mol, smult=smult)
@@ -46,14 +46,16 @@ def diatomic (atom1, atom2, r, fnal, basis, ncas, nelecas, nstates,
     return mc.nuc_grad_method ()
 
 def setUpModule():
-    global diatomic, original_grids
+    global mols, diatomic, original_grids
+    mols = []
     original_grids = dft.radi.ATOM_SPECIFIC_TREUTLER_GRIDS
     dft.radi.ATOM_SPECIFIC_TREUTLER_GRIDS = False
 
 def tearDownModule():
-    global diatomic, original_grids
+    global mols, diatomic, original_grids
     dft.radi.ATOM_SPECIFIC_TREUTLER_GRIDS = original_grids
-    del diatomic, original_grids
+    [m.stdout.close() for m in mols]
+    del diatomic, original_grids, mols
 
 # The purpose of these separate test functions is to narrow down an error to specific degrees of
 # freedom. But if the lih_cms2ftpbe and _df cases pass, then almost certainly, they all pass.
@@ -155,6 +157,20 @@ def test_grad_lih_cms2ftpbe22_sto3g (self):
          with self.subTest (state=i):
             de = mc_grad.kernel (state=i) [1,0] / BOHR
             self.assertAlmostEqual (de, de_ref[i], 5)
+    
+    def test_grad_lih_cms2tm06l22_sto3g (self):
+        # z_orb:    yes
+        # z_ci:     yes
+        # z_is:     yes
+        mc_grad = diatomic ('Li', 'H', 0.8, 'tM06L', 'STO-3G', 2, 2, 2, grids_level=1)
+        de_ref = [-1.03428938, -0.88278628]
+
+        # Numerical from this software
+        for i in range (2):
+         with self.subTest (state=i):
+            de = mc_grad.kernel (state=i) [1,0]/BOHR
+        
+            self.assertAlmostEqual (de, de_ref[i], 5)
 
     # MRH 05/05/2023: currently, the only other test which uses DF-MC-PDFT features is
     # test_grad_h2co, which is slower than this. Therefore I am restoring it.

pyscf/grad/test/test_grad_metagga_mcpdft.py

@@ -0,0 +1,127 @@
+#!/usr/bin/env python
+# Copyright 2014-2025 The PySCF Developers. All Rights Reserved.
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#     http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+#
+# Author: Matthew Hennefarth <[email protected]>
+
+import unittest
+
+from pyscf import scf, gto, df, dft
+from pyscf.data.nist import BOHR
+from pyscf import mcpdft
+
+
+def diatomic(
+    atom1,
+    atom2,
+    r,
+    fnal,
+    basis,
+    ncas,
+    nelecas,
+    nstates,
+    charge=None,
+    spin=None,
+    symmetry=False,
+    cas_irrep=None,
+    density_fit=False,
+    grids_level=9,
+):
+    """Used for checking diatomic systems to see if the Lagrange Multipliers are working properly."""
+    global mols
+    xyz = "{:s} 0.0 0.0 0.0; {:s} {:.3f} 0.0 0.0".format(atom1, atom2, r)
+    mol = gto.M(
+        atom=xyz,
+        basis=basis,
+        charge=charge,
+        spin=spin,
+        symmetry=symmetry,
+        verbose=0,
+        output="/dev/null",
+    )
+    mols.append(mol)
+    mf = scf.RHF(mol)
+    if density_fit:
+        mf = mf.density_fit(auxbasis=df.aug_etb(mol))
+
+    mc = mcpdft.CASSCF(mf.run(), fnal, ncas, nelecas, grids_level=grids_level)
+    if spin is None:
+        spin = mol.nelectron % 2
+
+    ss = spin * (spin + 2) * 0.25
+    mc.fix_spin_(ss=ss, shift=2)
+
+    if nstates > 1:
+        mc = mc.state_average(
+            [
+                1.0 / float(nstates),
+            ]
+            * nstates,
+        )
+
+    mc.conv_tol = 1e-12
+    mc.conv_grad_tol = 1e-6
+    mo = None
+    if symmetry and (cas_irrep is not None):
+        mo = mc.sort_mo_by_irrep(cas_irrep)
+
+    mc_grad = mc.run(mo).nuc_grad_method()
+    mc_grad.conv_rtol = 1e-12
+    return mc_grad
+
+
+def setUpModule():
+    global mols, original_grids
+    mols = []
+    original_grids = dft.radi.ATOM_SPECIFIC_TREUTLER_GRIDS
+    dft.radi.ATOM_SPECIFIC_TREUTLER_GRIDS = False
+
+
+def tearDownModule():
+    global mols, diatomic, original_grids
+    [m.stdout.close() for m in mols]
+    dft.radi.ATOM_SPECIFIC_TREUTLER_GRIDS = original_grids
+    del mols, diatomic, original_grids
+
+
+class KnownValues(unittest.TestCase):
+
+    def test_grad_lih_sstm06l22_sto3g(self):
+        mc = diatomic("Li", "H", 0.8, "tM06L", "STO-3G", 2, 2, 1, grids_level=1)
+        de = mc.kernel()[1, 0] / BOHR
+
+        # Numerical from this software
+        # PySCF commit:         f2c2d3f963916fb64ae77241f1b44f24fa484d96
+        # PySCF-forge commit:   4015363355dc691a80bc94d4b2b094318b213e36
+        DE_REF = -1.0546009263404388
+
+        self.assertAlmostEqual(de, DE_REF, 5)
+
+    def test_grad_lih_sa2tm06l22_sto3g(self):
+        mc = diatomic("Li", "H", 0.8, "tM06L", "STO-3G", 2, 2, 2, grids_level=1)
+
+        # Numerical from this software
+        # PySCF commit:         f2c2d3f963916fb64ae77241f1b44f24fa484d96
+        # PySCF-forge commit:   4015363355dc691a80bc94d4b2b094318b213e36
+        DE_REF = [-1.0351271000, -0.8919881992]
+
+        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")
+    unittest.main()