diff --git a/pyscf/neo/grad.py b/pyscf/neo/grad.py index 31de81ece4..8d97bd9012 100644 --- a/pyscf/neo/grad.py +++ b/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 = = + # = + + 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 diff --git a/pyscf/neo/test/test_qmmm.py b/pyscf/neo/test/test_qmmm.py index f80bc417b0..3f5274a66d 100644 --- a/pyscf/neo/test/test_qmmm.py +++ b/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() diff --git a/pyscf/qmmm/itrf.py b/pyscf/qmmm/itrf.py index 3bb79202fb..f50e891321 100644 --- a/pyscf/qmmm/itrf.py +++ b/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 = = + # = + + 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