Becke weight Hessian (#413)

* Becke Hessian complete offdiagonal term working

* Becke Hessian same derivative direction diagonal term working

* Becke Hessian main diagonal working

* Improvement according to suggestions

Author: henryw7

gpu4pyscf/hessian/rks.py

@@ -1383,6 +1383,40 @@ def get_dweight_dA(mol, grids):
 
     return dweight_dA
 
+def get_d2weight_dAdB(mol, grids):
+    ngrids = grids.coords.shape[0]
+    assert grids.atm_idx.shape[0] == ngrids
+    assert grids.quadrature_weights.shape[0] == ngrids
+    atm_coords = cupy.asarray(mol.atom_coords(), order = "C")
+
+    from gpu4pyscf.dft import radi
+    a_factor = radi.get_treutler_fac(mol, grids.atomic_radii)
+
+    d2weight_dAdB = cupy.zeros([mol.natm, mol.natm, 3, 3, ngrids], order = "C")
+    libgdft.GDFTbecke_partition_weight_second_derivative(
+        ctypes.cast(d2weight_dAdB.data.ptr, ctypes.c_void_p),
+        ctypes.cast(grids.coords.data.ptr, ctypes.c_void_p),
+        ctypes.cast(grids.quadrature_weights.data.ptr, ctypes.c_void_p),
+        ctypes.cast(atm_coords.data.ptr, ctypes.c_void_p),
+        ctypes.cast(a_factor.data.ptr, ctypes.c_void_p),
+        ctypes.cast(grids.atm_idx.data.ptr, ctypes.c_void_p),
+        ctypes.c_int(ngrids),
+        ctypes.c_int(mol.natm),
+    )
+
+    range_ngrids = cupy.arange(ngrids)
+    for i_atom in range(mol.natm):
+        for i_xyz in range(3):
+            for j_xyz in range(3):
+                d2weight_dAdB[i_atom, grids.atm_idx, i_xyz, j_xyz, range_ngrids] = -cupy.sum(d2weight_dAdB[i_atom, :, i_xyz, j_xyz, :], axis=[0])
+
+    for i_atom in range(mol.natm):
+        for i_xyz in range(3):
+            for j_xyz in range(3):
+                d2weight_dAdB[grids.atm_idx, i_atom, i_xyz, j_xyz, range_ngrids] = -cupy.sum(d2weight_dAdB[:, i_atom, i_xyz, j_xyz, :], axis=[0])
+
+    return d2weight_dAdB
+
 def _get_vnlc_deriv1(hessobj, mo_coeff, mo_occ, max_memory):
     """
         Equation notation follows:

gpu4pyscf/hessian/tests/test_vv10_hessian.py

@@ -14,10 +14,12 @@
 
 import unittest
 import numpy as np
+import cupy as cp
 import pyscf
 from gpu4pyscf.dft import rks
 from gpu4pyscf.hessian.rks import _get_vnlc_deriv1, _get_vnlc_deriv1_numerical, \
-                                  _get_enlc_deriv2, _get_enlc_deriv2_numerical
+                                  _get_enlc_deriv2, _get_enlc_deriv2_numerical, \
+                                  get_dweight_dA, get_d2weight_dAdB
 
 def setUpModule():
     global mol
@@ -425,6 +427,41 @@ def test_vv10_fock_first_derivative(self):
 
         assert np.linalg.norm(test_dF - reference_dF) < 1e-8
 
+    def test_becke_second_derivative(self):
+        mf = rks.RKS(mol, xc = "PBE")
+        mf.grids.atom_grid = (50,194)
+        mf.grids.build()
+        grids = mf.grids
+
+        test_d2w = get_d2weight_dAdB(mol, grids)
+
+        reference_d2w = cp.empty([mol.natm, mol.natm, 3, 3, grids.coords.shape[0]])
+        dx = 1e-5
+        mol_copy = mol.copy()
+        for i_atom in range(mol.natm):
+            for i_xyz in range(3):
+                xyz_p = mol.atom_coords()
+                xyz_p[i_atom, i_xyz] += dx
+                mol_copy.set_geom_(xyz_p, unit='Bohr')
+                mol_copy.build()
+                grids.reset(mol_copy)
+                grids.build()
+                w_p = get_dweight_dA(mol_copy, grids)
+
+                xyz_m = mol.atom_coords()
+                xyz_m[i_atom, i_xyz] -= dx
+                mol_copy.set_geom_(xyz_m, unit='Bohr')
+                mol_copy.build()
+                grids.reset(mol_copy)
+                grids.build()
+                w_m = get_dweight_dA(mol_copy, grids)
+
+                reference_d2w[i_atom, :, i_xyz, :, :] = (w_p - w_m) / (2 * dx)
+        grids.build(mol)
+
+        assert cp.max(cp.abs(test_d2w - reference_d2w)) < 1e-7
+
+
 if __name__ == "__main__":
     print("Full Tests for RKS Hessian with VV10")
     unittest.main()