CNEO-QM/MM: support pure point-charge model MM gradients

Though the computational efficiency of current implementation is
questionable because the parallelism might be limited.

pyscf/neo/grad.py

@@ -1,7 +1,7 @@
 #!/usr/bin/env python
 
 '''
-Analytical nuclear gradient for constrained nuclear-electronic orbital
+Analytic nuclear gradient for constrained nuclear-electronic orbital
 '''
 import numpy
 from pyscf import df, gto, lib, neo
@@ -219,22 +219,33 @@ def hcore_deriv(atm_id):
 def grad_hcore_mm(mm_mol, mol_n, dm_n):
     coords = mm_mol.atom_coords()
     charges = mm_mol.atom_charges()
-    expnts = mm_mol.get_zetas()
-
-    intor = 'int3c2e_ip2'
-    nao = mol_n.nao
-    max_memory = mol_n.super_mol.max_memory - lib.current_memory()[0]
-    blksize = int(min(max_memory*1e6/8/nao**2/3, 200))
-    blksize = max(blksize, 1)
-    cintopt = gto.moleintor.make_cintopt(mol_n._atm, mol_n._bas,
-                                         mol_n._env, intor)
-
     g = numpy.empty_like(coords)
-    for i0, i1 in lib.prange(0, charges.size, blksize):
-        fakemol = gto.fakemol_for_charges(coords[i0:i1], expnts[i0:i1])
-        j3c = df.incore.aux_e2(mol_n, fakemol, intor, aosym='s1',
-                               comp=3, cintopt=cintopt)
-        g[i0:i1] = numpy.einsum('ipqk,qp->ik', j3c * charges[i0:i1], dm_n).T
+    if mm_mol.charge_model == 'gaussian':
+        expnts = mm_mol.get_zetas()
+
+        intor = 'int3c2e_ip2'
+        nao = mol_n.nao
+        max_memory = mol_n.super_mol.max_memory - lib.current_memory()[0]
+        blksize = int(min(max_memory*1e6/8/nao**2/3, 200))
+        blksize = max(blksize, 1)
+        cintopt = gto.moleintor.make_cintopt(mol_n._atm, mol_n._bas,
+                                             mol_n._env, intor)
+
+        for i0, i1 in lib.prange(0, charges.size, blksize):
+            fakemol = gto.fakemol_for_charges(coords[i0:i1], expnts[i0:i1])
+            j3c = df.incore.aux_e2(mol_n, fakemol, intor, aosym='s1',
+                                   comp=3, cintopt=cintopt)
+            g[i0:i1] = numpy.einsum('ipqk,qp->ik', j3c * charges[i0:i1], dm_n).T
+    else:
+        # From examples/qmmm/30-force_on_mm_particles.py
+        # The interaction between electron density and MM particles
+        # d/dR <i| (1/|r-R|) |j> = <i| d/dR (1/|r-R|) |j> = <i| -d/dr (1/|r-R|) |j>
+        #   = <d/dr i| (1/|r-R|) |j> + <i| (1/|r-R|) |d/dr j>
+        for i, q in enumerate(charges):
+            with mol_n.with_rinv_origin(coords[i]):
+                v = mol_n.intor('int1e_iprinv')
+            g[i] = (numpy.einsum('ij,xji->x', dm_n, v) +
+                    numpy.einsum('ij,xij->x', dm_n, v.conj())) * -q
     ia = mol_n.atom_index
     charge = mol_n.super_mol.atom_charge(ia)
     return -charge * g

pyscf/neo/test/test_qmmm.py

@@ -25,13 +25,16 @@
 
 
 def setUpModule():
-    global mol, mol_1e6, mol_dft, mm_coords, mm_charges, mm_radii, mm_mol
+    global mol, mol_1e6, mol_dft, mm_coords, mm_charges, mm_radii, mm_mol, \
+           mm_mol_point, mm_mol_small_radius, mol_point, mol_small_radius
     mm_coords = [(1.369, 0.146,-0.395),
                  (1.894, 0.486, 0.335),
                  (0.451, 0.165,-0.083)]
     mm_charges = [-1.040, 0.520, 0.520]
     mm_radii = [0.63, 0.32, 0.32]
     mm_mol = create_mm_mol(mm_coords, mm_charges, mm_radii)
+    mm_mol_point = create_mm_mol(mm_coords, mm_charges)
+    mm_mol_small_radius = create_mm_mol(mm_coords, mm_charges, [1e-8]*3)
     atom='''O       -1.464   0.099   0.300
             H       -1.956   0.624  -0.340
             H       -1.797  -0.799   0.206'''
@@ -40,9 +43,14 @@ def setUpModule():
     mol_1e6 = neo.M(atom=atom, basis='631G', nuc_basis='1e6',
                     quantum_nuc=['H'], mm_mol=mm_mol)
     mol_dft = gto.M(atom=atom, basis='631G')
+    mol_point = neo.M(atom=atom, basis='631G', nuc_basis='pb4d',
+                      quantum_nuc=['H'], mm_mol=mm_mol_point)
+    mol_small_radius = neo.M(atom=atom, basis='631G', nuc_basis='pb4d',
+                             quantum_nuc=['H'], mm_mol=mm_mol_small_radius)
 
 def tearDownModule():
-    global mol, mol_1e6, mol_dft, mm_coords, mm_charges, mm_radii, mm_mol
+    global mol, mol_1e6, mol_dft, mm_coords, mm_charges, mm_radii, mm_mol, \
+           mm_mol_point, mm_mol_small_radius, mol_point, mol_small_radius
 
 class KnowValues(unittest.TestCase):
     def test(self):
@@ -155,6 +163,19 @@ def test_no_mm(self):
         numpy.testing.assert_array_almost_equal(g_qm0, g_qm2)
         self.assertAlmostEqual(numpy.linalg.norm(g_mm2), 0.0, 6)
 
+    def test_point_and_small_radius_gaussian(self):
+        mf_point = neo.CDFT(mol_point)
+        mf_point.mf_elec.xc = 'PBE0'
+        mf_small_radius = neo.CDFT(mol_small_radius)
+        mf_small_radius.mf_elec.xc = 'PBE0'
+        self.assertAlmostEqual(mf_point.kernel(), mf_small_radius.kernel())
+
+        g_point = mf_point.Gradients()
+        g_small_radius = mf_small_radius.Gradients()
+        numpy.testing.assert_array_almost_equal(g_point.grad(), g_small_radius.grad())
+        numpy.testing.assert_array_almost_equal(g_point.grad_mm(), g_small_radius.grad_mm())
+
+
 if __name__ == "__main__":
     print("Full Tests for CNEO-QMMM.")
     unittest.main()

pyscf/qmmm/itrf.py

@@ -341,22 +341,33 @@ def grad_hcore_mm(self, dm, mol=None):
 
         coords = mm_mol.atom_coords()
         charges = mm_mol.atom_charges()
-        expnts = mm_mol.get_zetas()
+        g = numpy.empty_like(coords)
+        if mm_mol.charge_model == 'gaussian':
+            expnts = mm_mol.get_zetas()
 
-        intor = 'int3c2e_ip2'
-        nao = mol.nao
-        max_memory = self.max_memory - lib.current_memory()[0]
-        blksize = int(min(max_memory*1e6/8/nao**2/3, 200))
-        blksize = max(blksize, 1)
-        cintopt = gto.moleintor.make_cintopt(mol._atm, mol._bas,
-                                             mol._env, intor)
+            intor = 'int3c2e_ip2'
+            nao = mol.nao
+            max_memory = self.max_memory - lib.current_memory()[0]
+            blksize = int(min(max_memory*1e6/8/nao**2/3, 200))
+            blksize = max(blksize, 1)
+            cintopt = gto.moleintor.make_cintopt(mol._atm, mol._bas,
+                                                 mol._env, intor)
 
-        g = numpy.empty_like(coords)
-        for i0, i1 in lib.prange(0, charges.size, blksize):
-            fakemol = gto.fakemol_for_charges(coords[i0:i1], expnts[i0:i1])
-            j3c = df.incore.aux_e2(mol, fakemol, intor, aosym='s1',
-                                   comp=3, cintopt=cintopt)
-            g[i0:i1] = numpy.einsum('ipqk,qp->ik', j3c * charges[i0:i1], dm).T
+            for i0, i1 in lib.prange(0, charges.size, blksize):
+                fakemol = gto.fakemol_for_charges(coords[i0:i1], expnts[i0:i1])
+                j3c = df.incore.aux_e2(mol, fakemol, intor, aosym='s1',
+                                       comp=3, cintopt=cintopt)
+                g[i0:i1] = numpy.einsum('ipqk,qp->ik', j3c * charges[i0:i1], dm).T
+        else:
+            # From examples/qmmm/30-force_on_mm_particles.py
+            # The interaction between electron density and MM particles
+            # d/dR <i| (1/|r-R|) |j> = <i| d/dR (1/|r-R|) |j> = <i| -d/dr (1/|r-R|) |j>
+            #   = <d/dr i| (1/|r-R|) |j> + <i| (1/|r-R|) |d/dr j>
+            for i, q in enumerate(charges):
+                with mol.with_rinv_origin(coords[i]):
+                    v = mol.intor('int1e_iprinv')
+                g[i] = (numpy.einsum('ij,xji->x', dm, v) +
+                        numpy.einsum('ij,xij->x', dm, v.conj())) * -q
         return g
 
     contract_hcore_mm = grad_hcore_mm # for backward compatibility