diff --git a/gpu4pyscf/gto/tests/test_int1e_grids_ip.py b/gpu4pyscf/gto/tests/test_int1e_grids_ip.py index 56f87e4b..de68266b 100644 --- a/gpu4pyscf/gto/tests/test_int1e_grids_ip.py +++ b/gpu4pyscf/gto/tests/test_int1e_grids_ip.py @@ -364,5 +364,5 @@ def test_int1e_grids_ip1_density_contracted(self): cp.testing.assert_allclose(ref_int1e_dA, test_int1e_dA, atol = integral_threshold) if __name__ == "__main__": - print("Full Tests for One Electron Coulomb Integrals") + print("Full Tests for One Electron Coulomb Integrals 1st Derivative") unittest.main() diff --git a/gpu4pyscf/gto/tests/test_int1e_grids_ipip.py b/gpu4pyscf/gto/tests/test_int1e_grids_ipip.py new file mode 100644 index 00000000..18f9fed9 --- /dev/null +++ b/gpu4pyscf/gto/tests/test_int1e_grids_ipip.py @@ -0,0 +1,480 @@ +# Copyright 2024 The GPU4PySCF Authors. All Rights Reserved. +# +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . + +import unittest +import numpy as np +import cupy as cp +import pyscf +from pyscf import lib, gto, df +from gpu4pyscf.gto.int3c1e_ipip import int1e_grids_ipip1, int1e_grids_ipvip1, int1e_grids_ip1ip2, int1e_grids_ipip2 + +def setUpModule(): + global mol_sph, mol_cart, grid_points, integral_threshold, density_contraction_threshold, charge_contraction_threshold + atom = ''' +O 0.0000 0.7375 -0.0528 +O 0.0000 -0.7375 -0.1528 +H 0.8190 0.8170 0.4220 +H -0.8190 -0.8170 0.4220 +''' + bas='def2-qzvpp' + + mol_sph = pyscf.M(atom=atom, basis=bas, max_memory=32000) + mol_sph.output = '/dev/null' + mol_sph.verbose = 0 + mol_sph.build() + + mol_cart = pyscf.M(atom=atom, basis=bas, max_memory=32000, cart=True) + mol_cart.output = '/dev/null' + mol_cart.verbose = 0 + mol_cart.build() + + xs = np.arange(-2.01, 2.0, 0.5) + ys = np.arange(-2.02, 2.0, 0.5) + zs = np.arange(-2.03, 2.0, 0.5) + grid_points = lib.cartesian_prod([xs, ys, zs]) + + # All of the following thresholds bound the max value of the corresponding matrix / tensor. + integral_threshold = 1e-12 + density_contraction_threshold = 1e-10 + charge_contraction_threshold = 1e-12 + +def tearDownModule(): + global mol_sph, mol_cart, grid_points + mol_sph.stdout.close() + mol_cart.stdout.close() + del mol_sph, mol_cart, grid_points + +class KnownValues(unittest.TestCase): + ''' + Values are compared to PySCF CPU intor() function + ''' + def test_int1e_grids_ipip1_charge_contracted_cart(self): + np.random.seed(12345) + charges = np.random.uniform(-2.0, 2.0, grid_points.shape[0]) + + mol = mol_cart + fakemol = gto.fakemol_for_charges(grid_points) + + int3c2e_ipip1 = mol._add_suffix('int3c2e_ipip1') + cintopt = gto.moleintor.make_cintopt(mol._atm, mol._bas, mol._env, int3c2e_ipip1) + v_nj = df.incore.aux_e2(mol, fakemol, intor=int3c2e_ipip1, aosym='s1', cintopt=cintopt) + ref_int1e_dAdA = np.einsum('dijq,q->dij', v_nj, charges) + ref_int1e_dAdA = ref_int1e_dAdA.reshape(3, 3, mol.nao, mol.nao) + + test_int1e_dAdA = int1e_grids_ipip1(mol, grid_points, charges = charges) + + assert isinstance(test_int1e_dAdA, cp.ndarray) + cp.testing.assert_allclose(ref_int1e_dAdA, test_int1e_dAdA, atol = integral_threshold) + + def test_int1e_grids_ipip1_charge_contracted_sph(self): + np.random.seed(12345) + charges = np.random.uniform(-2.0, 2.0, grid_points.shape[0]) + + mol = mol_sph + fakemol = gto.fakemol_for_charges(grid_points) + + int3c2e_ipip1 = mol._add_suffix('int3c2e_ipip1') + cintopt = gto.moleintor.make_cintopt(mol._atm, mol._bas, mol._env, int3c2e_ipip1) + v_nj = df.incore.aux_e2(mol, fakemol, intor=int3c2e_ipip1, aosym='s1', cintopt=cintopt) + ref_int1e_dAdA = np.einsum('dijq,q->dij', v_nj, charges) + ref_int1e_dAdA = ref_int1e_dAdA.reshape(3, 3, mol.nao, mol.nao) + + test_int1e_dAdA = int1e_grids_ipip1(mol, grid_points, charges = charges) + + assert isinstance(test_int1e_dAdA, cp.ndarray) + cp.testing.assert_allclose(ref_int1e_dAdA, test_int1e_dAdA, atol = integral_threshold) + + def test_int1e_grids_ipip1_charge_contracted_gaussian_charge(self): + np.random.seed(12345) + charges = np.random.uniform(-2.0, 2.0, grid_points.shape[0]) + charge_exponents = np.random.uniform(0.5, 1.0, grid_points.shape[0]) + + mol = mol_sph + fakemol = gto.fakemol_for_charges(grid_points, expnt=charge_exponents) + + int3c2e_ipip1 = mol._add_suffix('int3c2e_ipip1') + cintopt = gto.moleintor.make_cintopt(mol._atm, mol._bas, mol._env, int3c2e_ipip1) + v_nj = df.incore.aux_e2(mol, fakemol, intor=int3c2e_ipip1, aosym='s1', cintopt=cintopt) + ref_int1e_dAdA = np.einsum('dijq,q->dij', v_nj, charges) + ref_int1e_dAdA = ref_int1e_dAdA.reshape(3, 3, mol.nao, mol.nao) + + test_int1e_dAdA = int1e_grids_ipip1(mol, grid_points, charges = charges, charge_exponents = charge_exponents) + + assert isinstance(test_int1e_dAdA, cp.ndarray) + cp.testing.assert_allclose(ref_int1e_dAdA, test_int1e_dAdA, atol = integral_threshold) + + def test_int1e_grids_ipip1_charge_contracted_omega(self): + np.random.seed(12345) + charges = np.random.uniform(-2.0, 2.0, grid_points.shape[0]) + + omega = 0.8 + mol_sph_omega = mol_sph.copy() + mol_sph_omega.set_range_coulomb(omega) + + mol = mol_sph_omega + fakemol = gto.fakemol_for_charges(grid_points) + + int3c2e_ipip1 = mol._add_suffix('int3c2e_ipip1') + cintopt = gto.moleintor.make_cintopt(mol._atm, mol._bas, mol._env, int3c2e_ipip1) + v_nj = df.incore.aux_e2(mol, fakemol, intor=int3c2e_ipip1, aosym='s1', cintopt=cintopt) + ref_int1e_dAdA = np.einsum('dijq,q->dij', v_nj, charges) + ref_int1e_dAdA = ref_int1e_dAdA.reshape(3, 3, mol.nao, mol.nao) + + test_int1e_dAdA = int1e_grids_ipip1(mol, grid_points, charges = charges) + + assert isinstance(test_int1e_dAdA, cp.ndarray) + cp.testing.assert_allclose(ref_int1e_dAdA, test_int1e_dAdA, atol = integral_threshold) + + def test_int1e_grids_ipip1_charge_contracted_gaussian_charge_omega(self): + np.random.seed(12345) + charges = np.random.uniform(-2.0, 2.0, grid_points.shape[0]) + charge_exponents = np.random.uniform(0.5, 1.0, grid_points.shape[0]) + + omega = 0.8 + mol_sph_omega = mol_sph.copy() + mol_sph_omega.set_range_coulomb(omega) + + mol = mol_sph_omega + fakemol = gto.fakemol_for_charges(grid_points, expnt=charge_exponents) + + int3c2e_ipip1 = mol._add_suffix('int3c2e_ipip1') + cintopt = gto.moleintor.make_cintopt(mol._atm, mol._bas, mol._env, int3c2e_ipip1) + v_nj = df.incore.aux_e2(mol, fakemol, intor=int3c2e_ipip1, aosym='s1', cintopt=cintopt) + ref_int1e_dAdA = np.einsum('dijq,q->dij', v_nj, charges) + ref_int1e_dAdA = ref_int1e_dAdA.reshape(3, 3, mol.nao, mol.nao) + + test_int1e_dAdA = int1e_grids_ipip1(mol, grid_points, charges = charges, charge_exponents = charge_exponents) + + assert isinstance(test_int1e_dAdA, cp.ndarray) + cp.testing.assert_allclose(ref_int1e_dAdA, test_int1e_dAdA, atol = integral_threshold) + + # ^ ipip1 v ipvip1 + + def test_int1e_grids_ipvip1_charge_contracted_cart(self): + np.random.seed(12345) + charges = np.random.uniform(-2.0, 2.0, grid_points.shape[0]) + + mol = mol_cart + fakemol = gto.fakemol_for_charges(grid_points) + + int3c2e_ipvip1 = mol._add_suffix('int3c2e_ipvip1') + cintopt = gto.moleintor.make_cintopt(mol._atm, mol._bas, mol._env, int3c2e_ipvip1) + v_nj = df.incore.aux_e2(mol, fakemol, intor=int3c2e_ipvip1, aosym='s1', cintopt=cintopt) + ref_int1e_dAdB = np.einsum('dijq,q->dij', v_nj, charges) + ref_int1e_dAdB = ref_int1e_dAdB.reshape(3, 3, mol.nao, mol.nao) + + test_int1e_dAdB = int1e_grids_ipvip1(mol, grid_points, charges = charges) + + assert isinstance(test_int1e_dAdB, cp.ndarray) + cp.testing.assert_allclose(ref_int1e_dAdB, test_int1e_dAdB, atol = integral_threshold) + + def test_int1e_grids_ipvip1_charge_contracted_sph(self): + np.random.seed(12345) + charges = np.random.uniform(-2.0, 2.0, grid_points.shape[0]) + + mol = mol_sph + fakemol = gto.fakemol_for_charges(grid_points) + + int3c2e_ipvip1 = mol._add_suffix('int3c2e_ipvip1') + cintopt = gto.moleintor.make_cintopt(mol._atm, mol._bas, mol._env, int3c2e_ipvip1) + v_nj = df.incore.aux_e2(mol, fakemol, intor=int3c2e_ipvip1, aosym='s1', cintopt=cintopt) + ref_int1e_dAdB = np.einsum('dijq,q->dij', v_nj, charges) + ref_int1e_dAdB = ref_int1e_dAdB.reshape(3, 3, mol.nao, mol.nao) + + test_int1e_dAdB = int1e_grids_ipvip1(mol, grid_points, charges = charges) + + assert isinstance(test_int1e_dAdB, cp.ndarray) + cp.testing.assert_allclose(ref_int1e_dAdB, test_int1e_dAdB, atol = integral_threshold) + + def test_int1e_grids_ipvip1_charge_contracted_gaussian_charge(self): + np.random.seed(12345) + charges = np.random.uniform(-2.0, 2.0, grid_points.shape[0]) + charge_exponents = np.random.uniform(0.5, 1.0, grid_points.shape[0]) + + mol = mol_sph + fakemol = gto.fakemol_for_charges(grid_points, expnt=charge_exponents) + + int3c2e_ipvip1 = mol._add_suffix('int3c2e_ipvip1') + cintopt = gto.moleintor.make_cintopt(mol._atm, mol._bas, mol._env, int3c2e_ipvip1) + v_nj = df.incore.aux_e2(mol, fakemol, intor=int3c2e_ipvip1, aosym='s1', cintopt=cintopt) + ref_int1e_dAdB = np.einsum('dijq,q->dij', v_nj, charges) + ref_int1e_dAdB = ref_int1e_dAdB.reshape(3, 3, mol.nao, mol.nao) + + test_int1e_dAdB = int1e_grids_ipvip1(mol, grid_points, charges = charges, charge_exponents = charge_exponents) + + assert isinstance(test_int1e_dAdB, cp.ndarray) + cp.testing.assert_allclose(ref_int1e_dAdB, test_int1e_dAdB, atol = integral_threshold) + + def test_int1e_grids_ipvip1_charge_contracted_omega(self): + np.random.seed(12345) + charges = np.random.uniform(-2.0, 2.0, grid_points.shape[0]) + + omega = 0.8 + mol_sph_omega = mol_sph.copy() + mol_sph_omega.set_range_coulomb(omega) + + mol = mol_sph_omega + fakemol = gto.fakemol_for_charges(grid_points) + + int3c2e_ipvip1 = mol._add_suffix('int3c2e_ipvip1') + cintopt = gto.moleintor.make_cintopt(mol._atm, mol._bas, mol._env, int3c2e_ipvip1) + v_nj = df.incore.aux_e2(mol, fakemol, intor=int3c2e_ipvip1, aosym='s1', cintopt=cintopt) + ref_int1e_dAdB = np.einsum('dijq,q->dij', v_nj, charges) + ref_int1e_dAdB = ref_int1e_dAdB.reshape(3, 3, mol.nao, mol.nao) + + test_int1e_dAdB = int1e_grids_ipvip1(mol, grid_points, charges = charges) + + assert isinstance(test_int1e_dAdB, cp.ndarray) + cp.testing.assert_allclose(ref_int1e_dAdB, test_int1e_dAdB, atol = integral_threshold) + + def test_int1e_grids_ipvip1_charge_contracted_gaussian_charge_omega(self): + np.random.seed(12345) + charges = np.random.uniform(-2.0, 2.0, grid_points.shape[0]) + charge_exponents = np.random.uniform(0.5, 1.0, grid_points.shape[0]) + + omega = 0.8 + mol_sph_omega = mol_sph.copy() + mol_sph_omega.set_range_coulomb(omega) + + mol = mol_sph_omega + fakemol = gto.fakemol_for_charges(grid_points, expnt=charge_exponents) + + int3c2e_ipvip1 = mol._add_suffix('int3c2e_ipvip1') + cintopt = gto.moleintor.make_cintopt(mol._atm, mol._bas, mol._env, int3c2e_ipvip1) + v_nj = df.incore.aux_e2(mol, fakemol, intor=int3c2e_ipvip1, aosym='s1', cintopt=cintopt) + ref_int1e_dAdB = np.einsum('dijq,q->dij', v_nj, charges) + ref_int1e_dAdB = ref_int1e_dAdB.reshape(3, 3, mol.nao, mol.nao) + + test_int1e_dAdB = int1e_grids_ipvip1(mol, grid_points, charges = charges, charge_exponents = charge_exponents) + + assert isinstance(test_int1e_dAdB, cp.ndarray) + cp.testing.assert_allclose(ref_int1e_dAdB, test_int1e_dAdB, atol = integral_threshold) + + # ^ ipvip1 v ip1ip2 + + def test_int1e_grids_ip1ip2_charge_contracted_cart(self): + np.random.seed(12345) + charges = np.random.uniform(-2.0, 2.0, grid_points.shape[0]) + + mol = mol_cart + fakemol = gto.fakemol_for_charges(grid_points) + + int3c2e_ip1ip2 = mol._add_suffix('int3c2e_ip1ip2') + cintopt = gto.moleintor.make_cintopt(mol._atm, mol._bas, mol._env, int3c2e_ip1ip2) + v_nj = df.incore.aux_e2(mol, fakemol, intor=int3c2e_ip1ip2, aosym='s1', cintopt=cintopt) + ref_int1e_dAdC = np.einsum('dijq,q->dij', v_nj, charges) + ref_int1e_dAdC = ref_int1e_dAdC.reshape(3, 3, mol.nao, mol.nao) + + test_int1e_dAdC = int1e_grids_ip1ip2(mol, grid_points, charges = charges) + + assert isinstance(test_int1e_dAdC, cp.ndarray) + cp.testing.assert_allclose(ref_int1e_dAdC, test_int1e_dAdC, atol = integral_threshold) + + def test_int1e_grids_ip1ip2_charge_contracted_sph(self): + np.random.seed(12345) + charges = np.random.uniform(-2.0, 2.0, grid_points.shape[0]) + + mol = mol_sph + fakemol = gto.fakemol_for_charges(grid_points) + + int3c2e_ip1ip2 = mol._add_suffix('int3c2e_ip1ip2') + cintopt = gto.moleintor.make_cintopt(mol._atm, mol._bas, mol._env, int3c2e_ip1ip2) + v_nj = df.incore.aux_e2(mol, fakemol, intor=int3c2e_ip1ip2, aosym='s1', cintopt=cintopt) + ref_int1e_dAdC = np.einsum('dijq,q->dij', v_nj, charges) + ref_int1e_dAdC = ref_int1e_dAdC.reshape(3, 3, mol.nao, mol.nao) + + test_int1e_dAdC = int1e_grids_ip1ip2(mol, grid_points, charges = charges) + + assert isinstance(test_int1e_dAdC, cp.ndarray) + cp.testing.assert_allclose(ref_int1e_dAdC, test_int1e_dAdC, atol = integral_threshold) + + def test_int1e_grids_ip1ip2_charge_contracted_gaussian_charge(self): + np.random.seed(12345) + charges = np.random.uniform(-2.0, 2.0, grid_points.shape[0]) + charge_exponents = np.random.uniform(0.5, 1.0, grid_points.shape[0]) + + mol = mol_sph + fakemol = gto.fakemol_for_charges(grid_points, expnt=charge_exponents) + + int3c2e_ip1ip2 = mol._add_suffix('int3c2e_ip1ip2') + cintopt = gto.moleintor.make_cintopt(mol._atm, mol._bas, mol._env, int3c2e_ip1ip2) + v_nj = df.incore.aux_e2(mol, fakemol, intor=int3c2e_ip1ip2, aosym='s1', cintopt=cintopt) + ref_int1e_dAdC = np.einsum('dijq,q->dij', v_nj, charges) + ref_int1e_dAdC = ref_int1e_dAdC.reshape(3, 3, mol.nao, mol.nao) + + test_int1e_dAdC = int1e_grids_ip1ip2(mol, grid_points, charges = charges, charge_exponents = charge_exponents) + + assert isinstance(test_int1e_dAdC, cp.ndarray) + cp.testing.assert_allclose(ref_int1e_dAdC, test_int1e_dAdC, atol = integral_threshold) + + def test_int1e_grids_ip1ip2_charge_contracted_omega(self): + np.random.seed(12345) + charges = np.random.uniform(-2.0, 2.0, grid_points.shape[0]) + + omega = 0.8 + mol_sph_omega = mol_sph.copy() + mol_sph_omega.set_range_coulomb(omega) + + mol = mol_sph_omega + fakemol = gto.fakemol_for_charges(grid_points) + + int3c2e_ip1ip2 = mol._add_suffix('int3c2e_ip1ip2') + cintopt = gto.moleintor.make_cintopt(mol._atm, mol._bas, mol._env, int3c2e_ip1ip2) + v_nj = df.incore.aux_e2(mol, fakemol, intor=int3c2e_ip1ip2, aosym='s1', cintopt=cintopt) + ref_int1e_dAdC = np.einsum('dijq,q->dij', v_nj, charges) + ref_int1e_dAdC = ref_int1e_dAdC.reshape(3, 3, mol.nao, mol.nao) + + test_int1e_dAdC = int1e_grids_ip1ip2(mol, grid_points, charges = charges) + + assert isinstance(test_int1e_dAdC, cp.ndarray) + cp.testing.assert_allclose(ref_int1e_dAdC, test_int1e_dAdC, atol = integral_threshold) + + def test_int1e_grids_ip1ip2_charge_contracted_gaussian_charge_omega(self): + np.random.seed(12345) + charges = np.random.uniform(-2.0, 2.0, grid_points.shape[0]) + charge_exponents = np.random.uniform(0.5, 1.0, grid_points.shape[0]) + + omega = 0.8 + mol_sph_omega = mol_sph.copy() + mol_sph_omega.set_range_coulomb(omega) + + mol = mol_sph_omega + fakemol = gto.fakemol_for_charges(grid_points, expnt=charge_exponents) + + int3c2e_ip1ip2 = mol._add_suffix('int3c2e_ip1ip2') + cintopt = gto.moleintor.make_cintopt(mol._atm, mol._bas, mol._env, int3c2e_ip1ip2) + v_nj = df.incore.aux_e2(mol, fakemol, intor=int3c2e_ip1ip2, aosym='s1', cintopt=cintopt) + ref_int1e_dAdC = np.einsum('dijq,q->dij', v_nj, charges) + ref_int1e_dAdC = ref_int1e_dAdC.reshape(3, 3, mol.nao, mol.nao) + + test_int1e_dAdC = int1e_grids_ip1ip2(mol, grid_points, charges = charges, charge_exponents = charge_exponents) + + assert isinstance(test_int1e_dAdC, cp.ndarray) + cp.testing.assert_allclose(ref_int1e_dAdC, test_int1e_dAdC, atol = integral_threshold) + + # ^ ip1ip2 v ipip2 + + def test_int1e_grids_ipip2_charge_contracted_cart(self): + np.random.seed(12345) + dm = np.random.uniform(-2.0, 2.0, (mol_cart.nao, mol_cart.nao)) + + mol = mol_cart + fakemol = gto.fakemol_for_charges(grid_points) + + # Note: we cannot compute ipip2 (dCdC) directly due to numerical problems, + # pyscf treat a point charge as a sharp Gaussian, and we cannot take 2nd derivative of it. + int3c2e_ip1ip2 = mol._add_suffix('int3c2e_ip1ip2') + cintopt = gto.moleintor.make_cintopt(mol._atm, mol._bas, mol._env, int3c2e_ip1ip2) + v_nj = df.incore.aux_e2(mol, fakemol, intor=int3c2e_ip1ip2, aosym='s1', cintopt=cintopt) + v_nj = -v_nj - v_nj.transpose(0, 2, 1, 3) # dCdC = -dAdC - dBdC + ref_int1e_dCdC = np.einsum('dijq,ij->dq', v_nj, dm) + ref_int1e_dCdC = ref_int1e_dCdC.reshape(3, 3, grid_points.shape[0]) + + test_int1e_dCdC = int1e_grids_ipip2(mol, grid_points, dm = dm) + + assert isinstance(test_int1e_dCdC, cp.ndarray) + cp.testing.assert_allclose(ref_int1e_dCdC, test_int1e_dCdC, atol = integral_threshold) + + def test_int1e_grids_ipip2_charge_contracted_sph(self): + np.random.seed(12345) + dm = np.random.uniform(-2.0, 2.0, (mol_sph.nao, mol_sph.nao)) + + mol = mol_sph + fakemol = gto.fakemol_for_charges(grid_points) + + # Note: we cannot compute ipip2 (dCdC) directly due to numerical problems, + # pyscf treat a point charge as a sharp Gaussian, and we cannot take 2nd derivative of it. + int3c2e_ip1ip2 = mol._add_suffix('int3c2e_ip1ip2') + cintopt = gto.moleintor.make_cintopt(mol._atm, mol._bas, mol._env, int3c2e_ip1ip2) + v_nj = df.incore.aux_e2(mol, fakemol, intor=int3c2e_ip1ip2, aosym='s1', cintopt=cintopt) + v_nj = -v_nj - v_nj.transpose(0, 2, 1, 3) # dCdC = -dAdC - dBdC + ref_int1e_dCdC = np.einsum('dijq,ij->dq', v_nj, dm) + ref_int1e_dCdC = ref_int1e_dCdC.reshape(3, 3, grid_points.shape[0]) + + test_int1e_dCdC = int1e_grids_ipip2(mol, grid_points, dm = dm) + + assert isinstance(test_int1e_dCdC, cp.ndarray) + cp.testing.assert_allclose(ref_int1e_dCdC, test_int1e_dCdC, atol = integral_threshold) + + def test_int1e_grids_ipip2_charge_contracted_gaussian_charge(self): + np.random.seed(12345) + dm = np.random.uniform(-2.0, 2.0, (mol_sph.nao, mol_sph.nao)) + charge_exponents = np.random.uniform(0.5, 1.0, grid_points.shape[0]) + + mol = mol_sph + fakemol = gto.fakemol_for_charges(grid_points, expnt=charge_exponents) + + int3c2e_ipip2 = mol._add_suffix('int3c2e_ipip2') + cintopt = gto.moleintor.make_cintopt(mol._atm, mol._bas, mol._env, int3c2e_ipip2) + v_nj = df.incore.aux_e2(mol, fakemol, intor=int3c2e_ipip2, aosym='s1', cintopt=cintopt) + ref_int1e_dCdC = np.einsum('dijq,ij->dq', v_nj, dm) + ref_int1e_dCdC = ref_int1e_dCdC.reshape(3, 3, grid_points.shape[0]) + + test_int1e_dCdC = int1e_grids_ipip2(mol, grid_points, dm = dm, charge_exponents = charge_exponents) + + assert isinstance(test_int1e_dCdC, cp.ndarray) + cp.testing.assert_allclose(ref_int1e_dCdC, test_int1e_dCdC, atol = integral_threshold) + + def test_int1e_grids_ipip2_charge_contracted_omega(self): + np.random.seed(12345) + dm = np.random.uniform(-2.0, 2.0, (mol_sph.nao, mol_sph.nao)) + + omega = 0.8 + mol_sph_omega = mol_sph.copy() + mol_sph_omega.set_range_coulomb(omega) + + mol = mol_sph_omega + fakemol = gto.fakemol_for_charges(grid_points) + + # Note: we cannot compute ipip2 (dCdC) directly due to numerical problems, + # pyscf treat a point charge as a sharp Gaussian, and we cannot take 2nd derivative of it. + int3c2e_ip1ip2 = mol._add_suffix('int3c2e_ip1ip2') + cintopt = gto.moleintor.make_cintopt(mol._atm, mol._bas, mol._env, int3c2e_ip1ip2) + v_nj = df.incore.aux_e2(mol, fakemol, intor=int3c2e_ip1ip2, aosym='s1', cintopt=cintopt) + v_nj = -v_nj - v_nj.transpose(0, 2, 1, 3) # dCdC = -dAdC - dBdC + ref_int1e_dCdC = np.einsum('dijq,ij->dq', v_nj, dm) + ref_int1e_dCdC = ref_int1e_dCdC.reshape(3, 3, grid_points.shape[0]) + + test_int1e_dCdC = int1e_grids_ipip2(mol, grid_points, dm = dm) + + assert isinstance(test_int1e_dCdC, cp.ndarray) + cp.testing.assert_allclose(ref_int1e_dCdC, test_int1e_dCdC, atol = integral_threshold) + + def test_int1e_grids_ipip2_charge_contracted_gaussian_charge_omega(self): + np.random.seed(12345) + dm = np.random.uniform(-2.0, 2.0, (mol_sph.nao, mol_sph.nao)) + charge_exponents = np.random.uniform(0.5, 1.0, grid_points.shape[0]) + + omega = 0.8 + mol_sph_omega = mol_sph.copy() + mol_sph_omega.set_range_coulomb(omega) + + mol = mol_sph_omega + fakemol = gto.fakemol_for_charges(grid_points, expnt=charge_exponents) + + int3c2e_ipip2 = mol._add_suffix('int3c2e_ipip2') + cintopt = gto.moleintor.make_cintopt(mol._atm, mol._bas, mol._env, int3c2e_ipip2) + v_nj = df.incore.aux_e2(mol, fakemol, intor=int3c2e_ipip2, aosym='s1', cintopt=cintopt) + ref_int1e_dCdC = np.einsum('dijq,ij->dq', v_nj, dm) + ref_int1e_dCdC = ref_int1e_dCdC.reshape(3, 3, grid_points.shape[0]) + + test_int1e_dCdC = int1e_grids_ipip2(mol, grid_points, dm = dm, charge_exponents = charge_exponents) + + assert isinstance(test_int1e_dCdC, cp.ndarray) + cp.testing.assert_allclose(ref_int1e_dCdC, test_int1e_dCdC, atol = integral_threshold) + +if __name__ == "__main__": + print("Full Tests for One Electron Coulomb Integrals 2nd Derivative") + unittest.main() diff --git a/gpu4pyscf/solvent/hessian/pcm.py b/gpu4pyscf/solvent/hessian/pcm.py index 719f0e02..a6e73b38 100644 --- a/gpu4pyscf/solvent/hessian/pcm.py +++ b/gpu4pyscf/solvent/hessian/pcm.py @@ -151,7 +151,7 @@ def get_d2D_d2S(surface, with_S=True, with_D=False, stream=None): raise RuntimeError('Failed in generating PCM d2D and d2S matrices.') return d2D, d2S -def hess_nuc(pcmobj, dm, verbose=None): +def analytical_hess_nuc(pcmobj, dm, verbose=None): if not pcmobj._intermediates: pcmobj.build() dm_cache = pcmobj._intermediates.get('dm', None) @@ -221,7 +221,7 @@ def hess_nuc(pcmobj, dm, verbose=None): t1 = log.timer_debug1('solvent hessian d(dVnuc/dx * q)/dx contribution', *t1) return d2e -def hess_qv(pcmobj, dm, verbose=None): +def analytical_hess_qv(pcmobj, dm, verbose=None): if not pcmobj._intermediates: pcmobj.build() dm_cache = pcmobj._intermediates.get('dm', None) @@ -238,9 +238,6 @@ def hess_qv(pcmobj, dm, verbose=None): grid_coords = pcmobj.surface['grid_coords'] q_sym = pcmobj._intermediates['q_sym'] - ngrids = q_sym.shape[0] - nao = mol.nao - aoslice = mol.aoslice_by_atom() aoslice = numpy.array(aoslice) @@ -368,7 +365,7 @@ def get_v_dot_d2D_dot_q(d2D, v_left, q_right, natom, gridslice): def get_v_dot_d2DT_dot_q(d2D, v_left, q_right, natom, gridslice): return get_v_dot_d2D_dot_q(d2D.transpose(0,1,3,2), v_left, q_right, natom, gridslice) -def hess_solver(pcmobj, dm): +def analytical_hess_solver(pcmobj, dm): if not pcmobj._intermediates: pcmobj.build() dm_cache = pcmobj._intermediates.get('dm', None) @@ -625,44 +622,6 @@ def hess_solver(pcmobj, dm): t1 = log.timer_debug1('solvent hessian d(V * dK-1R/dx * V)/dx contribution', *t1) return d2e -def hess_elec(pcmobj, dm, verbose=None): - ''' - slow version with finite difference - TODO: use analytical hess_nuc - ''' - log = logger.new_logger(pcmobj, verbose) - t1 = log.init_timer() - pmol = pcmobj.mol.copy() - mol = pmol.copy() - coords = mol.atom_coords(unit='Bohr') - - def pcm_grad_scanner(mol): - # TODO: use more analytical forms - pcmobj.reset(mol) - e, v = pcmobj._get_vind(dm) - #return grad_elec(pcmobj, dm) - pcm_grad = grad_nuc(pcmobj, dm) - pcm_grad+= grad_solver(pcmobj, dm) - pcm_grad+= grad_qv(pcmobj, dm) - return pcm_grad - - mol.verbose = 0 - de = numpy.zeros([mol.natm, mol.natm, 3, 3]) - eps = 1e-3 - for ia in range(mol.natm): - for ix in range(3): - dv = numpy.zeros_like(coords) - dv[ia,ix] = eps - mol.set_geom_(coords + dv, unit='Bohr') - g0 = pcm_grad_scanner(mol) - - mol.set_geom_(coords - dv, unit='Bohr') - g1 = pcm_grad_scanner(mol) - de[ia,:,ix] = (g0 - g1)/2.0/eps - t1 = log.timer_debug1('solvent energy', *t1) - pcmobj.reset(pmol) - return de - def get_dqsym_dx_fix_vgrids(pcmobj, atmlst, inverse_K): assert pcmobj._intermediates is not None @@ -852,7 +811,7 @@ def get_dqsym_dx(pcmobj, dm, atmlst, intopt_derivative): inverse_K = cupy.linalg.inv(K) return get_dqsym_dx_fix_vgrids(pcmobj, atmlst, inverse_K) + get_dqsym_dx_fix_K_R(pcmobj, dm, atmlst, inverse_K, intopt_derivative) -def analytic_grad_vmat(pcmobj, dm, mo_coeff, mo_occ, atmlst=None, verbose=None): +def analytical_grad_vmat(pcmobj, dm, mo_coeff, mo_occ, atmlst=None, verbose=None): ''' dv_solv / da ''' @@ -961,8 +920,9 @@ def kernel(self, *args, dm=None, atmlst=None, **kwargs): dm = dm[0] + dm[1] is_equilibrium = self.base.with_solvent.equilibrium_solvation self.base.with_solvent.equilibrium_solvation = True - self.de_solvent = hess_elec(self.base.with_solvent, dm, verbose=self.verbose) - #self.de_solvent+= hess_nuc(self.base.with_solvent) + self.de_solvent = analytical_hess_nuc(self.base.with_solvent, dm, verbose=self.verbose) + self.de_solvent += analytical_hess_qv(self.base.with_solvent, dm, verbose=self.verbose) + self.de_solvent += analytical_hess_solver(self.base.with_solvent, dm, verbose=self.verbose) self.de_solute = super().kernel(*args, **kwargs) self.de = self.de_solute + self.de_solvent self.base.with_solvent.equilibrium_solvation = is_equilibrium @@ -974,7 +934,7 @@ def make_h1(self, mo_coeff, mo_occ, chkfile=None, atmlst=None, verbose=None): h1ao = super().make_h1(mo_coeff, mo_occ, atmlst=atmlst, verbose=verbose) if isinstance(self.base, scf.hf.RHF): dm = self.base.make_rdm1(ao_repr=True) - dv = analytic_grad_vmat(self.base.with_solvent, dm, mo_coeff, mo_occ, atmlst=atmlst, verbose=verbose) + dv = analytical_grad_vmat(self.base.with_solvent, dm, mo_coeff, mo_occ, atmlst=atmlst, verbose=verbose) for i0, ia in enumerate(atmlst): h1ao[i0] += dv[i0] return h1ao @@ -983,8 +943,8 @@ def make_h1(self, mo_coeff, mo_occ, chkfile=None, atmlst=None, verbose=None): solvent = self.base.with_solvent dm = self.base.make_rdm1(ao_repr=True) dm = dm[0] + dm[1] - dva = analytic_grad_vmat(solvent, dm, mo_coeff[0], mo_occ[0], atmlst=atmlst, verbose=verbose) - dvb = analytic_grad_vmat(solvent, dm, mo_coeff[1], mo_occ[1], atmlst=atmlst, verbose=verbose) + dva = analytical_grad_vmat(solvent, dm, mo_coeff[0], mo_occ[0], atmlst=atmlst, verbose=verbose) + dvb = analytical_grad_vmat(solvent, dm, mo_coeff[1], mo_occ[1], atmlst=atmlst, verbose=verbose) for i0, ia in enumerate(atmlst): h1aoa[i0] += dva[i0] h1aob[i0] += dvb[i0] diff --git a/gpu4pyscf/solvent/hessian/smd.py b/gpu4pyscf/solvent/hessian/smd.py index 49897d74..6039a730 100644 --- a/gpu4pyscf/solvent/hessian/smd.py +++ b/gpu4pyscf/solvent/hessian/smd.py @@ -154,7 +154,7 @@ def make_h1(self, mo_coeff, mo_occ, chkfile=None, atmlst=None, verbose=None): h1ao = super().make_h1(mo_coeff, mo_occ, atmlst=atmlst, verbose=verbose) if isinstance(self.base, scf.hf.RHF): dm = self.base.make_rdm1(ao_repr=True) - dv = pcm_hess.analytic_grad_vmat(self.base.with_solvent, dm, mo_coeff, mo_occ, atmlst=atmlst, verbose=verbose) + dv = pcm_hess.analytical_grad_vmat(self.base.with_solvent, dm, mo_coeff, mo_occ, atmlst=atmlst, verbose=verbose) for i0, ia in enumerate(atmlst): h1ao[i0] += dv[i0] return h1ao @@ -163,8 +163,8 @@ def make_h1(self, mo_coeff, mo_occ, chkfile=None, atmlst=None, verbose=None): solvent = self.base.with_solvent dm = self.base.make_rdm1(ao_repr=True) dm = dm[0] + dm[1] - dva = pcm_hess.analytic_grad_vmat(solvent, dm, mo_coeff[0], mo_occ[0], atmlst=atmlst, verbose=verbose) - dvb = pcm_hess.analytic_grad_vmat(solvent, dm, mo_coeff[1], mo_occ[1], atmlst=atmlst, verbose=verbose) + dva = pcm_hess.analytical_grad_vmat(solvent, dm, mo_coeff[0], mo_occ[0], atmlst=atmlst, verbose=verbose) + dvb = pcm_hess.analytical_grad_vmat(solvent, dm, mo_coeff[1], mo_occ[1], atmlst=atmlst, verbose=verbose) for i0, ia in enumerate(atmlst): h1aoa[i0] += dva[i0] h1aob[i0] += dvb[i0] diff --git a/gpu4pyscf/solvent/tests/test_pcm_hessian.py b/gpu4pyscf/solvent/tests/test_pcm_hessian.py index c7076f29..6e19ec96 100644 --- a/gpu4pyscf/solvent/tests/test_pcm_hessian.py +++ b/gpu4pyscf/solvent/tests/test_pcm_hessian.py @@ -21,7 +21,7 @@ from gpu4pyscf.solvent import pcm from gpu4pyscf import scf, dft from packaging import version -from gpu4pyscf.solvent.hessian.pcm import analytic_grad_vmat +from gpu4pyscf.solvent.hessian.pcm import analytical_grad_vmat, analytical_hess_nuc, analytical_hess_solver, analytical_hess_qv from gpu4pyscf.lib.cupy_helper import contract pyscf_25 = version.parse(pyscf.__version__) <= version.parse('2.5.0') @@ -130,6 +130,37 @@ def pcm_vmat_scanner(mol): pcmobj.reset(pmol) return vmat +def _fd_hess_contribution(pcmobj, dm, gradient_function): + pmol = pcmobj.mol.copy() + mol = pmol.copy() + coords = mol.atom_coords(unit='Bohr') + + def pcm_grad_scanner(mol): + pcmobj.reset(mol) + e, v = pcmobj._get_vind(dm) + pcm_grad = gradient_function(pcmobj, dm) + # pcm_grad = grad_nuc(pcmobj, dm) + # pcm_grad+= grad_solver(pcmobj, dm) + # pcm_grad+= grad_qv(pcmobj, dm) + return pcm_grad + + mol.verbose = 0 + de = np.zeros([mol.natm, mol.natm, 3, 3]) + eps = 1e-5 + for ia in range(mol.natm): + for ix in range(3): + dv = np.zeros_like(coords) + dv[ia,ix] = eps + mol.set_geom_(coords + dv, unit='Bohr') + g0 = pcm_grad_scanner(mol) + + mol.set_geom_(coords - dv, unit='Bohr') + g1 = pcm_grad_scanner(mol) + + de[ia,:,ix,:] = (g0 - g1)/2.0/eps + pcmobj.reset(pmol) + return de + @unittest.skipIf(pcm.libsolvent is None, "solvent extension not compiled") class KnownValues(unittest.TestCase): def test_df_hess_cpcm(self): @@ -192,7 +223,7 @@ def test_grad_vmat_cpcm(self): mo_coeff = mf.mo_coeff mo_occ = mf.mo_occ - test_grad_vmat = analytic_grad_vmat(hobj.base.with_solvent, dm, mo_coeff, mo_occ) + test_grad_vmat = analytical_grad_vmat(hobj.base.with_solvent, dm, mo_coeff, mo_occ) ref_grad_vmat = _fd_grad_vmat(hobj.base.with_solvent, dm, mo_coeff, mo_occ) cp.testing.assert_allclose(ref_grad_vmat, test_grad_vmat, atol = 1e-10) @@ -206,7 +237,7 @@ def test_grad_vmat_iefpcm(self): mo_coeff = mf.mo_coeff mo_occ = mf.mo_occ - test_grad_vmat = analytic_grad_vmat(hobj.base.with_solvent, dm, mo_coeff, mo_occ) + test_grad_vmat = analytical_grad_vmat(hobj.base.with_solvent, dm, mo_coeff, mo_occ) ref_grad_vmat = _fd_grad_vmat(hobj.base.with_solvent, dm, mo_coeff, mo_occ) cp.testing.assert_allclose(ref_grad_vmat, test_grad_vmat, atol = 1e-10) @@ -220,11 +251,71 @@ def test_grad_vmat_ssvpe(self): mo_coeff = mf.mo_coeff mo_occ = mf.mo_occ - test_grad_vmat = analytic_grad_vmat(hobj.base.with_solvent, dm, mo_coeff, mo_occ) + test_grad_vmat = analytical_grad_vmat(hobj.base.with_solvent, dm, mo_coeff, mo_occ) ref_grad_vmat = _fd_grad_vmat(hobj.base.with_solvent, dm, mo_coeff, mo_occ) cp.testing.assert_allclose(ref_grad_vmat, test_grad_vmat, atol = 1e-10) + def test_hess_nuc_iefpcm(self): + print("testing IEF-PCM d2E_nuc/dx2") + mf = _make_mf(method='IEF-PCM') + hobj = mf.Hessian() + dm = mf.make_rdm1() + + test_grad_vmat = analytical_hess_nuc(hobj.base.with_solvent, dm) + from gpu4pyscf.solvent.grad.pcm import grad_nuc + ref_grad_vmat = _fd_hess_contribution(hobj.base.with_solvent, dm, grad_nuc) + + cp.testing.assert_allclose(ref_grad_vmat, test_grad_vmat, atol = 1e-10) + + def test_hess_qv_iefpcm(self): + print("testing IEF-PCM d2E_elec/dx2") + mf = _make_mf(method='IEF-PCM') + hobj = mf.Hessian() + dm = mf.make_rdm1() + + test_grad_vmat = analytical_hess_qv(hobj.base.with_solvent, dm) + from gpu4pyscf.solvent.grad.pcm import grad_qv + ref_grad_vmat = _fd_hess_contribution(hobj.base.with_solvent, dm, grad_qv) + + cp.testing.assert_allclose(ref_grad_vmat, test_grad_vmat, atol = 1e-10) + + def test_hess_solver_cpcm(self): + print("testing C-PCM d2E_KR/dx2") + mf = _make_mf(method='C-PCM') + hobj = mf.Hessian() + dm = mf.make_rdm1() + + test_grad_vmat = analytical_hess_solver(hobj.base.with_solvent, dm) + from gpu4pyscf.solvent.grad.pcm import grad_solver + ref_grad_vmat = _fd_hess_contribution(hobj.base.with_solvent, dm, grad_solver) + + cp.testing.assert_allclose(ref_grad_vmat, test_grad_vmat, atol = 1e-10) + + def test_hess_solver_iefpcm(self): + print("testing IEF-PCM d2E_KR/dx2") + mf = _make_mf(method='IEF-PCM') + hobj = mf.Hessian() + dm = mf.make_rdm1() + + test_grad_vmat = analytical_hess_solver(hobj.base.with_solvent, dm) + from gpu4pyscf.solvent.grad.pcm import grad_solver + ref_grad_vmat = _fd_hess_contribution(hobj.base.with_solvent, dm, grad_solver) + + cp.testing.assert_allclose(ref_grad_vmat, test_grad_vmat, atol = 1e-10) + + def test_hess_solver_ssvpe(self): + print("testing SS(V)PE d2E_KR/dx2") + mf = _make_mf(method='SS(V)PE') + hobj = mf.Hessian() + dm = mf.make_rdm1() + + test_grad_vmat = analytical_hess_solver(hobj.base.with_solvent, dm) + from gpu4pyscf.solvent.grad.pcm import grad_solver + ref_grad_vmat = _fd_hess_contribution(hobj.base.with_solvent, dm, grad_solver) + + cp.testing.assert_allclose(ref_grad_vmat, test_grad_vmat, atol = 1e-10) + @pytest.mark.skipif(pyscf_25, reason='requires pyscf 2.6 or higher') def test_to_gpu(self): import pyscf