diff --git a/examples/qmmm/31-gaussian_charge.py b/examples/qmmm/31-gaussian_charge.py new file mode 100644 index 0000000000..dc60a19574 --- /dev/null +++ b/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) diff --git a/pyscf/gto/mole.py b/pyscf/gto/mole.py index b95da64a79..97e3c650ba 100644 --- a/pyscf/gto/mole.py +++ b/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 diff --git a/pyscf/qmmm/itrf.py b/pyscf/qmmm/itrf.py index 1833150078..afd9ff05a8 100644 --- a/pyscf/qmmm/itrf.py +++ b/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 diff --git a/pyscf/qmmm/mm_mole.py b/pyscf/qmmm/mm_mole.py index c4712cb19d..b71beb13d0 100644 --- a/pyscf/qmmm/mm_mole.py +++ b/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) diff --git a/pyscf/qmmm/test/test_grad_mm.py b/pyscf/qmmm/test/test_grad_mm.py new file mode 100644 index 0000000000..58f9835f74 --- /dev/null +++ b/pyscf/qmmm/test/test_grad_mm.py @@ -0,0 +1,79 @@ +#!/usr/bin/env python +# Copyright 2014-2023 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. + +import unittest +import numpy +import pyscf +from pyscf import lib +from pyscf import gto +from pyscf import scf +from pyscf import grad +from pyscf.qmmm import itrf + +def setUpModule(): + global mol, mm_coords, mm_charges, mm_radii + mol = gto.M( + verbose = 5, + output = '/dev/null', + 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), + (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] + +def tearDownModule(): + global mol, mm_coords, mm_charges, mm_radii + + +class KnowValues(unittest.TestCase): + def test_grad_mm(self): + mf = itrf.mm_charge(scf.RHF(mol), mm_coords, mm_charges, mm_radii) + e_hf = mf.kernel() + self.assertAlmostEqual(e_hf, -76.00028338468692, 8) + + # qm grad + g_hf = itrf.mm_charge_grad(grad.RHF(mf), mm_coords, mm_charges, mm_radii) + g_hf_qm = g_hf.kernel() + self.assertAlmostEqual(numpy.linalg.norm(g_hf_qm), 0.031003880070795818, 6) + + # mm grad + g_hf_mm_h1 = g_hf.grad_hcore_mm(mf.make_rdm1()) + g_hf_mm_nuc = g_hf.grad_nuc_mm() + self.assertAlmostEqual(numpy.linalg.norm(g_hf_mm_h1), 0.5108777013806877, 6) + self.assertAlmostEqual(numpy.linalg.norm(g_hf_mm_nuc), 0.4915404602273757, 6) + + # finite difference for MM atoms + mm_coords1 = [(1.369, 0.147,-0.395), + (1.894, 0.486, 0.335), + (0.451, 0.165,-0.083)] + mf1 = itrf.mm_charge(scf.RHF(mol), mm_coords1, mm_charges, mm_radii) + e1 = mf1.kernel() + + mm_coords2 = [(1.369, 0.145,-0.395), + (1.894, 0.486, 0.335), + (0.451, 0.165,-0.083)] + mf2 = itrf.mm_charge(scf.RHF(mol), mm_coords2, mm_charges, mm_radii) + e2 = mf2.kernel() + self.assertAlmostEqual((e1 - e2) / 0.002*lib.param.BOHR, + (g_hf_mm_h1+g_hf_mm_nuc)[0,1], 6) + +if __name__ == "__main__": + print("Full Tests for qmmm MM force.") + unittest.main()