Allow Gaussian MM charges (pyscf#1677)

* finite size gaussian radii for MM charges in QM/MM

* minor fixes

* minor fix 2

* enable QM/MM to compute nuclear gradients w.r.t MM particles (pyscf#4)

* finite size gaussian radii for MM charges in QM/MM

* add mm gradient in qmmm

* add CP-CI for FCI

* remove cpci for a different PR

---------

Co-authored-by: MoleOrbitalHybridAnalyst <[email protected]>

examples/qmmm/31-gaussian_charge.py

@@ -0,0 +1,37 @@
+"""
+QMMM calculation with the MM charges being Gaussian distributions
+"""
+
+from pyscf import gto, scf, qmmm, grad
+
+mol = gto.M(
+    verbose = 3,
+    atom = '''O       -1.464   0.099   0.300
+              H       -1.956   0.624  -0.340
+              H       -1.797  -0.799   0.206''',
+    basis = '631G')
+
+mm_coords = [(1.369, 0.146,-0.395),  # O
+              (1.894, 0.486, 0.335), # H
+              (0.451, 0.165,-0.083)] # H
+mm_charges = [-1.040, 0.520, 0.520]
+mm_radii = [0.63, 0.32, 0.32] # radii of Gaussians
+
+
+mf = qmmm.mm_charge(scf.RHF(mol), mm_coords, mm_charges, mm_radii)
+e_hf = mf.kernel() # -76.00028338468692
+
+# QM grad
+g_hf = qmmm.mm_charge_grad(grad.RHF(mf), mm_coords, mm_charges, mm_radii)
+g_hf_qm = g_hf.kernel()
+print('Nuclear gradients of QM atoms:')
+print(g_hf_qm)
+
+# MM grad
+# NOTE For post-HF methods, the density matrix should
+# include the orbital response. See more details
+# in examples/qmmm/30-force_on_mm_particles.py
+g_hf_mm_h1 = g_hf.grad_hcore_mm(mf.make_rdm1())
+g_hf_mm_nuc = g_hf.grad_nuc_mm()
+print('Nuclear gradients of MM atoms:')
+print(g_hf_mm_h1 + g_hf_mm_nuc)

pyscf/gto/mole.py

@@ -64,7 +64,8 @@
 CHARGE_OF  = 0
 PTR_COORD  = 1
 NUC_MOD_OF = 2
-PTR_ZETA   = PTR_FRAC_CHARGE = 3
+PTR_ZETA   = 3
+PTR_FRAC_CHARGE = 4
 ATM_SLOTS  = 6
 ATOM_OF    = 0
 ANG_OF     = 1
@@ -3951,6 +3952,11 @@ def fakemol_for_charges(coords, expnt=1e16):
     distribution with the Gaussian exponent (expnt).
     '''
     nbas = coords.shape[0]
+    expnt = numpy.asarray(expnt).ravel()
+    if expnt.size == 1:
+        expnt = numpy.repeat(expnt, nbas)
+    assert expnt.size == nbas
+
     fakeatm = numpy.zeros((nbas,ATM_SLOTS), dtype=numpy.int32)
     fakebas = numpy.zeros((nbas,BAS_SLOTS), dtype=numpy.int32)
     fakeenv = [0] * PTR_ENV_START
@@ -3961,11 +3967,12 @@ def fakemol_for_charges(coords, expnt=1e16):
     fakebas[:,ATOM_OF] = numpy.arange(nbas)
     fakebas[:,NPRIM_OF] = 1
     fakebas[:,NCTR_OF] = 1
-# approximate point charge with gaussian distribution exp(-1e16*r^2)
-    fakebas[:,PTR_EXP] = ptr
-    fakebas[:,PTR_COEFF] = ptr+1
-    fakeenv.append([expnt, 1/(2*numpy.sqrt(numpy.pi)*gaussian_int(2,expnt))])
-    ptr += 2
+# approximate point charge with gaussian distribution exp(-expnt*r^2)
+    fakebas[:,PTR_EXP] = ptr + numpy.arange(nbas) * 2
+    fakebas[:,PTR_COEFF] = ptr + numpy.arange(nbas) * 2 + 1
+    coeff = 1 / (2 * numpy.sqrt(numpy.pi) * gaussian_int(2, expnt))
+    fakeenv.append(numpy.vstack((expnt, coeff)).T.ravel())
+
     fakemol = Mole()
     fakemol._atm = fakeatm
     fakemol._bas = fakebas

pyscf/qmmm/itrf.py

@@ -32,7 +32,7 @@
 from pyscf.qmmm import mm_mole
 
 
-def add_mm_charges(scf_method, atoms_or_coords, charges, unit=None):
+def add_mm_charges(scf_method, atoms_or_coords, charges, radii=None, unit=None):
     '''Embedding the one-electron (non-relativistic) potential generated by MM
     point charges into QM Hamiltonian.
 
@@ -50,7 +50,10 @@ def add_mm_charges(scf_method, atoms_or_coords, charges, unit=None):
             MM particle coordinates
         charges : 1D array
             MM particle charges
+
     Kwargs:
+        radii : 1D array
+            The Gaussian charge distribution radii of MM atoms.
         unit : str
             Bohr, AU, Ang (case insensitive). Default is the same to mol.unit
 
@@ -74,7 +77,8 @@ def add_mm_charges(scf_method, atoms_or_coords, charges, unit=None):
     mol = scf_method.mol
     if unit is None:
         unit = mol.unit
-    mm_mol = mm_mole.create_mm_mol(atoms_or_coords, charges, unit)
+    mm_mol = mm_mole.create_mm_mol(atoms_or_coords, charges,
+                                   radii=radii, unit=unit)
     return qmmm_for_scf(scf_method, mm_mol)
 
 # Define method mm_charge for backward compatibility
@@ -123,20 +127,42 @@ def dump_flags(self, verbose=None):
             return self
 
         def get_hcore(self, mol=None):
-            if mol is None: mol = self.mol
+            if mol is None:
+                mol = self.mol
+            mm_mol = self.mm_mol
+
             if getattr(method_class, 'get_hcore', None):
                 h1e = method_class.get_hcore(self, mol)
             else:  # DO NOT modify post-HF objects to avoid the MM charges applied twice
                 raise RuntimeError('mm_charge function cannot be applied on post-HF methods')
 
-            coords = self.mm_mol.atom_coords()
-            charges = self.mm_mol.atom_charges()
+            coords = mm_mol.atom_coords()
+            charges = mm_mol.atom_charges()
             nao = mol.nao
             max_memory = self.max_memory - lib.current_memory()[0]
             blksize = int(min(max_memory*1e6/8/nao**2, 200))
-            for i0, i1 in lib.prange(0, charges.size, blksize):
-                j3c = mol.intor('int1e_grids', hermi=1, grids=coords[i0:i1])
-                h1e += numpy.einsum('kpq,k->pq', j3c, -charges[i0:i1])
+            blksize = max(blksize, 1)
+            if mm_mol.charge_model == 'gaussian':
+                expnts = mm_mol.get_zetas()
+
+                if mol.cart:
+                    intor = 'int3c2e_cart'
+                else:
+                    intor = 'int3c2e_sph'
+                cintopt = gto.moleintor.make_cintopt(mol._atm, mol._bas,
+                                                     mol._env, intor)
+                v = 0
+                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=intor,
+                                           aosym='s2ij', cintopt=cintopt)
+                    v += numpy.einsum('xk,k->x', j3c, -charges[i0:i1])
+                v = lib.unpack_tril(v)
+                h1e += v
+            else:
+                for i0, i1 in lib.prange(0, charges.size, blksize):
+                    j3c = mol.intor('int1e_grids', hermi=1, grids=coords[i0:i1])
+                    h1e += numpy.einsum('kpq,k->pq', j3c, -charges[i0:i1])
             return h1e
 
         def energy_nuc(self):
@@ -163,7 +189,7 @@ def nuc_grad_method(self):
         scf_method.mo_energy = scf_method._scf.mo_energy
         return scf_method
 
-def add_mm_charges_grad(scf_grad, atoms_or_coords, charges, unit=None):
+def add_mm_charges_grad(scf_grad, atoms_or_coords, charges, radii=None, unit=None):
     '''Apply the MM charges in the QM gradients' method.  It affects both the
     electronic and nuclear parts of the QM fragment.
 
@@ -176,6 +202,8 @@ def add_mm_charges_grad(scf_grad, atoms_or_coords, charges, unit=None):
         charges : 1D array
             MM particle charges
     Kwargs:
+        radii : 1D array
+            The Gaussian charge distribution radii of MM atoms.
         unit : str
             Bohr, AU, Ang (case insensitive). Default is the same to mol.unit
 
@@ -198,7 +226,8 @@ def add_mm_charges_grad(scf_grad, atoms_or_coords, charges, unit=None):
     mol = scf_grad.mol
     if unit is None:
         unit = mol.unit
-    mm_mol = mm_mole.create_mm_mol(atoms_or_coords, charges, unit)
+    mm_mol = mm_mole.create_mm_mol(atoms_or_coords, charges,
+                                   radii=radii, unit=unit)
     mm_grad = qmmm_grad_for_scf(scf_grad)
     mm_grad.base.mm_mol = mm_mol
     return mm_grad
@@ -238,19 +267,77 @@ def dump_flags(self, verbose=None):
 
         def get_hcore(self, mol=None):
             ''' (QM 1e grad) + <-d/dX i|q_mm/r_mm|j>'''
-            if mol is None: mol = self.mol
-            coords = self.base.mm_mol.atom_coords()
-            charges = self.base.mm_mol.atom_charges()
+            if mol is None:
+                mol = self.mol
+            mm_mol = self.base.mm_mol
+            coords = mm_mol.atom_coords()
+            charges = mm_mol.atom_charges()
 
             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)
             g_qm = grad_class.get_hcore(self, mol)
-            for i0, i1 in lib.prange(0, charges.size, blksize):
-                j3c = mol.intor('int1e_grids_ip', grids=coords[i0:i1])
-                g_qm += numpy.einsum('ikpq,k->ipq', j3c, charges[i0:i1])
+            if mm_mol.charge_model == 'gaussian':
+                expnts = mm_mol.get_zetas()
+                if mol.cart:
+                    intor = 'int3c2e_ip1_cart'
+                else:
+                    intor = 'int3c2e_ip1_sph'
+                cintopt = gto.moleintor.make_cintopt(mol._atm, mol._bas,
+                                                     mol._env, intor)
+                v = 0
+                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)
+                    v += numpy.einsum('ipqk,k->ipq', j3c, charges[i0:i1])
+                g_qm += v
+            else:
+                for i0, i1 in lib.prange(0, charges.size, blksize):
+                    j3c = mol.intor('int1e_grids_ip', grids=coords[i0:i1])
+                    g_qm += numpy.einsum('ikpq,k->ipq', j3c, charges[i0:i1])
             return g_qm
 
+        def grad_hcore_mm(self, dm, mol=None):
+            r'''Nuclear gradients of the electronic energy
+            with respect to MM atoms:
+
+            ... math::
+                g = \sum_{ij} \frac{\partial hcore_{ij}}{\partial R_{I}} P_{ji},
+
+            where I represents MM atoms.
+
+            Args:
+                dm : array
+                    The QM density matrix.
+            '''
+            if mol is None:
+                mol = self.mol
+            mm_mol = self.base.mm_mol
+
+            coords = mm_mol.atom_coords()
+            charges = mm_mol.atom_charges()
+            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)
+
+            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
+            return g
+
+        contract_hcore_mm = grad_hcore_mm # for backward compatibility
+
         def grad_nuc(self, mol=None, atmlst=None):
             if mol is None: mol = self.mol
             coords = self.base.mm_mol.atom_coords()
@@ -267,6 +354,24 @@ def grad_nuc(self, mol=None, atmlst=None):
             if atmlst is not None:
                 g_mm = g_mm[atmlst]
             return g_qm + g_mm
+
+        def grad_nuc_mm(self, mol=None):
+            '''Nuclear gradients of the QM-MM nuclear energy
+            (in the form of point charge Coulomb interactions)
+            with respect to MM atoms.
+            '''
+            if mol is None:
+                mol = self.mol
+            mm_mol = self.base.mm_mol
+            coords = mm_mol.atom_coords()
+            charges = mm_mol.atom_charges()
+            g_mm = numpy.zeros_like(coords)
+            for i in range(mol.natm):
+                q1 = mol.atom_charge(i)
+                r1 = mol.atom_coord(i)
+                r = lib.norm(r1-coords, axis=1)
+                g_mm += q1 * numpy.einsum('i,ix,i->ix', charges, r1-coords, 1/r**3)
+            return g_mm
     return QMMM(scf_grad)
 
 # A tag to label the derived class

pyscf/qmmm/mm_mole.py

@@ -22,7 +22,9 @@
 
 import numpy
 from pyscf import gto
+from pyscf.gto.mole import is_au
 from pyscf.data.elements import charge
+from pyscf.lib import param
 
 class Mole(gto.Mole):
     '''Mole class for MM particles.
@@ -36,13 +38,18 @@ class Mole(gto.Mole):
             |  [atomN, (x, y, z)]]
 
     Kwargs:
-        charges : fractional charges of MM particles
+        charges : 1D array
+            fractional charges of MM particles
+        zeta : 1D array
+            Gaussian charge distribution parameter.
+            rho(r) = charge * Norm * exp(-zeta * r^2)
 
     '''
-    def __init__(self, atoms, charges=None):
+    def __init__(self, atoms, charges=None, zeta=None):
         gto.Mole.__init__(self)
         self.atom = self._atom = atoms
         self.unit = 'Bohr'
+        self.charge_model = 'point'
 
         # Initialize ._atm and ._env to save the coordinates and charges and
         # other info of MM particles
@@ -61,13 +68,38 @@ def __init__(self, atoms, charges=None):
                                  numpy.hstack((coords, charges)).ravel())
         _atm[:,gto.PTR_COORD] = gto.PTR_ENV_START + numpy.arange(natm) * 4
         _atm[:,gto.PTR_FRAC_CHARGE] = gto.PTR_ENV_START + numpy.arange(natm) * 4 + 3
-        self._atm = _atm
 
+        if zeta is not None:
+            self.charge_model = 'gaussian'
+            zeta = numpy.asarray(zeta, dtype=float).ravel()
+            self._env = numpy.append(self._env, zeta)
+            _atm[:,gto.PTR_ZETA] = gto.PTR_ENV_START + natm*4 + numpy.arange(natm)
+
+        self._atm = _atm
         self._built = True
 
-def create_mm_mol(atoms_or_coords, charges=None, unit='Angstrom'):
+    def get_zetas(self):
+        if self.charge_model == 'gaussian':
+            return self._env[self._atm[:,gto.PTR_ZETA]]
+        else:
+            return 1e16
+
+def create_mm_mol(atoms_or_coords, charges=None, radii=None, unit='Angstrom'):
     '''Create an MM object based on the given coordinates and charges of MM
     particles.
+
+    Args:
+        atoms_or_coords : array-like
+            Cartesian coordinates of MM atoms, in the form of a 2D array:
+            [(x1, y1, z1), (x2, y2, z2), ...]
+
+    Kwargs:
+        charges : 1D array
+            The charges of MM atoms.
+        radii : 1D array
+            The Gaussian charge distribuction radii of MM atoms.
+        unit : string
+            The unit of the input. Default is 'Angstrom'.
     '''
     if isinstance(atoms_or_coords, numpy.ndarray):
         # atoms_or_coords == np.array([(xx, xx, xx)])
@@ -82,5 +114,12 @@ def create_mm_mol(atoms_or_coords, charges=None, unit='Angstrom'):
     else:
         atoms = atoms_or_coords
     atoms = gto.format_atom(atoms, unit=unit)
-    return Mole(atoms, charges)
 
+    if radii is None:
+        zeta = None
+    else:
+        radii = numpy.asarray(radii, dtype=float).ravel()
+        if not is_au(unit):
+            radii = radii / param.BOHR
+        zeta = 1 / radii**2
+    return Mole(atoms, charges=charges, zeta=zeta)