pyxdh: UHF-UB3LYP grad

pyxdh/DerivOnce/__init__.py

@@ -1,17 +1,21 @@
 __all__ = [
     "DerivOnceSCF", "DerivOnceNCDFT",  # deriv_once_scf
-    "GradSCF", "GradNCDFT",  # grad_scf
     "DerivOnceMP2", "DerivOnceXDH",  # deriv_once_mp2
+    "GradSCF", "GradNCDFT",  # grad_scf
     "GradMP2", "GradXDH",  # grad_mp2
     "DipoleSCF", "DipoleNCDFT",  # dipole_scf
     "DipoleMP2", "DipoleXDH",  # dipole_mp2
-    "DerivOnceUSCF",  # deriv_once_uscf
+
+    "DerivOnceUSCF", "DerivOnceUNCDFT",  # deriv_once_uscf
+    "GradUSCF", "GradUNCDFT", # grad_uscf
 ]
 
 from pyxdh.DerivOnce.deriv_once_scf import DerivOnceSCF, DerivOnceNCDFT
-from pyxdh.DerivOnce.grad_scf import GradSCF, GradNCDFT
 from pyxdh.DerivOnce.deriv_once_mp2 import DerivOnceMP2, DerivOnceXDH
+from pyxdh.DerivOnce.grad_scf import GradSCF, GradNCDFT
 from pyxdh.DerivOnce.grad_mp2 import GradMP2, GradXDH
 from pyxdh.DerivOnce.dipole_scf import DipoleSCF, DipoleNCDFT
 from pyxdh.DerivOnce.dipole_mp2 import DipoleMP2, DipoleXDH
-from pyxdh.DerivOnce.deriv_once_uscf import DerivOnceUSCF
+
+from pyxdh.DerivOnce.deriv_once_uscf import DerivOnceUSCF, DerivOnceUNCDFT
+from pyxdh.DerivOnce.grad_uscf import GradUSCF, GradUNCDFT

pyxdh/DerivOnce/deriv_once_scf.py

@@ -387,7 +387,7 @@ def mol_slice(self, atm_id):
         _, _, p0, p1 = self.mol.aoslice_by_atom()[atm_id]
         return slice(p0, p1)
 
-    def Ax0_Core(self, si, sj, sk, sl, reshape=True, in_cphf=False):
+    def Ax0_Core(self, si, sj, sk, sl, reshape=True, in_cphf=False, C=None):
         """
 
         Parameters
@@ -408,7 +408,8 @@ def Ax0_Core(self, si, sj, sk, sl, reshape=True, in_cphf=False):
         -------
         fx : function which pass matrix, then return Ax @ X.
         """
-        C = self.C
+        if C is None:
+            C = self.C
         nao = self.nao
         resp = self.resp_cphf if in_cphf else self.resp
 

pyxdh/DerivOnce/deriv_once_uscf.py

@@ -1,16 +1,14 @@
 import numpy as np
-from abc import ABC, abstractmethod
+from abc import ABC
 from functools import partial
 import os
 import warnings
-import copy
-from pyscf.scf._response_functions import _gen_rhf_response
+from pyscf.scf._response_functions import _gen_uhf_response
 
-from pyscf import gto, dft, grad, hessian
-from pyscf.scf import cphf
+from pyscf import dft, grad, hessian
+from pyscf.scf import ucphf
 
-from pyxdh.Utilities import timing
-from pyxdh.DerivOnce.deriv_once_scf import DerivOnceSCF
+from pyxdh.DerivOnce.deriv_once_scf import DerivOnceSCF, DerivOnceNCDFT
 from pyxdh.Utilities import GridIterator, KernelHelper, timing
 
 MAXMEM = float(os.getenv("MAXMEM", 2))
@@ -45,6 +43,10 @@ def initialization_pyscf(self):
             self.scf_hess = hessian.UHF(self.scf_eng)
         return
 
+    @property
+    def nvir(self):
+        return self.nmo - self.nocc[0], self.nmo - self.nocc[1]
+
     @property
     def so(self):
         return slice(0, self.nocc[0]), slice(0, self.nocc[1])
@@ -100,6 +102,81 @@ def _get_F_1_mo(self):
             return 0
         return np.einsum("xAuv, xup, xvq -> xApq", self.F_1_ao, self.C, self.C)
 
+    @property
+    def resp(self):
+        if self._resp is NotImplemented:
+            self._resp = _gen_uhf_response(self.scf_eng, mo_coeff=self.C, mo_occ=self.mo_occ, hermi=1, max_memory=self.grdit_memory)
+        return self._resp
+
+    def _get_B_1(self):
+        sa = self.sa
+        so = self.so
+
+        B_1 = self.F_1_mo.copy()
+        if isinstance(self.S_1_mo, np.ndarray):
+            B_1 += (
+                - np.einsum("xApq, xq -> xApq", self.S_1_mo, self.e)
+                - 0.5 * np.array(self.Ax0_Core(sa, sa, so, so)((self.S_1_mo[0, :, so[0], so[0]], self.S_1_mo[1, :, so[1], so[1]])))
+            )
+        return B_1
 
-
+    @property
+    def resp_cphf(self):
+        if self._resp_cphf is NotImplemented:
+            if self.xc_type == "HF":
+                self._resp_cphf = self.resp
+            else:
+                mf = dft.RKS(self.mol)
+                mf.xc = self.scf_eng.xc
+                mf.grids = self.cphf_grids
+                self._resp_cphf = _gen_uhf_response(mf, mo_coeff=self.C, mo_occ=self.mo_occ, hermi=1, max_memory=self.grdit_memory)
+        return self._resp_cphf
+
+    def Ax0_Core(self, si, sj, sk, sl, reshape=True, in_cphf=False, C=None):
+        if C is None:
+            C = self.C
+        nao = self.nao
+        resp = self.resp_cphf if in_cphf else self.resp
+
+        @timing
+        def fx(X_):
+            # For simplicity, shape of X should be (2, dim_prop, sk, sl)
+            have_first_dim = len(X_[0].shape) == 3
+            prop_dim = X_[0].shape[0] if have_first_dim else 1
+
+            dm = np.zeros((2, prop_dim, nao, nao))
+            dm[0] = C[0][:, sk[0]] @ X_[0] @ C[0][:, sl[0]].T
+            dm[1] = C[1][:, sk[1]] @ X_[1] @ C[1][:, sl[1]].T
+            dm += dm.swapaxes(-1, -2)
+            if (not have_first_dim):
+                dm.shape = (2, nao, nao)
+
+            ax_ao = resp(dm)
+            Ax = (
+                C[0][:, si[0]].T @ ax_ao[0] @ C[0][:, sj[0]],
+                C[1][:, si[1]].T @ ax_ao[1] @ C[1][:, sj[1]]
+            )
+            return Ax
+
+        return fx
+
+
+class DerivOnceUNCDFT(DerivOnceUSCF, DerivOnceNCDFT):
+
+    def _get_Z(self):
+        so, sv = self.so, self.sv
+        nocc, nvir = self.nocc, self.nvir
+        Ax0_Core = self.Ax0_Core
+        e, mo_occ = self.e, self.mo_occ
+        F_0_mo = self.nc_deriv.F_0_mo
+        def fx(X):
+            X_ = (
+                X[:, :nocc[0] * nvir[0]].reshape((nvir[0], nocc[0])),
+                X[:, nocc[0] * nvir[0]:].reshape((nvir[1], nocc[1]))
+            )
+            return np.concatenate([v.ravel() for v in Ax0_Core(sv, so, sv, so, in_cphf=True)(X_)])
+        Z = ucphf.solve(fx, e, mo_occ, (F_0_mo[0, None, sv[0], so[0]], F_0_mo[1, None, sv[1], so[1]]), max_cycle=100, tol=self.cphf_tol)[0]
+        # output Z shape is (1, nvir, nocc), we remove the first dimension
+        Z = (Z[0][0], Z[1][0])
+        return Z
 

pyxdh/DerivOnce/grad_uscf.py

@@ -4,7 +4,7 @@
 
 from pyscf import grad
 
-from pyxdh.DerivOnce import DerivOnceUSCF, DerivOnceNCDFT, GradSCF
+from pyxdh.DerivOnce import DerivOnceUSCF, GradSCF, DerivOnceUNCDFT
 from pyxdh.Utilities import GridIterator, KernelHelper
 
 MAXMEM = float(os.getenv("MAXMEM", 2))
@@ -55,6 +55,25 @@ def _get_E_1(self):
         return E_1.reshape((natm, 3))
 
 
+class GradUNCDFT(DerivOnceUNCDFT, GradUSCF):
+
+    @property
+    def DerivOnceMethod(self):
+        return GradUSCF
+
+    def _get_E_1(self):
+        natm = self.natm
+        so, sv = self.so, self.sv
+        B_1 = self.B_1
+        Z = self.Z
+        E_1 = (
+            + self.nc_deriv.E_1
+            + 2 * np.einsum("ai, Aai -> A", Z[0], B_1[0, :, sv[0], so[0]]).reshape((natm, 3))
+            + 2 * np.einsum("ai, Aai -> A", Z[1], B_1[1, :, sv[1], so[1]]).reshape((natm, 3))
+        )
+        return E_1
+
+
 class Test_GradUSCF:
 
     def test_UHF_grad(self):
@@ -78,3 +97,15 @@ def test_UB3LYP_grad(self):
         helper = GradUSCF({"scf_eng": CH3.gga_eng})
         formchk = FormchkInterface(resource_filename("pyxdh", "Validation/gaussian/CH3-B3LYP-freq.fchk"))
         assert(np.allclose(helper.E_1, formchk.grad(), atol=1e-6, rtol=1e-4))
+
+    def test_UHF_UB3LYP_grad(self):
+
+        from pyxdh.Utilities.test_molecules import Mol_CH3
+        from pkg_resources import resource_filename
+        import pickle
+
+        CH3 = Mol_CH3()
+        helper = GradUNCDFT({"scf_eng": CH3.hf_eng, "nc_eng": CH3.gga_eng, "cphf_tol": 1e-10})
+        with open(resource_filename("pyxdh", "Validation/numerical_deriv/ncdft_derivonce_uhf_ub3lyp.dat"), "rb") as f:
+            ref_grad = pickle.load(f)["grad"].reshape(-1, 3)
+        assert (np.allclose(helper.E_1, ref_grad, atol=1e-5, rtol=1e-4))

pyxdh/Utilities/deriv_numerical.py

@@ -50,7 +50,7 @@ def move_mol(self, movelist):
         for tup in movelist:
             mol_ret_coords[tup[0], tup[1]] += tup[2] * self.interval * lib.param.BOHR
         mol_ret.set_geom_(mol_ret_coords)
-        return mol_ret
+        return mol_ret.build()
 
     def perform_mf(self):
         natm = self.mol.natm

pyxdh/Utilities/test_molecules.py

@@ -3,7 +3,7 @@
 
 class Mol_H2O2:
 
-    def __init__(self, xc="B3LYPg", mol=None):
+    def __init__(self, xc="B3LYPg", mol=None, atom_grid=(99, 590)):
 
         if mol is None:
             mol = gto.Mole()
@@ -19,6 +19,7 @@ def __init__(self, xc="B3LYPg", mol=None):
 
         self.mol = mol
         self.xc = xc
+        self.atom_grid = atom_grid
 
         self._hf_eng = NotImplemented
         self._hf_grad = NotImplemented
@@ -48,7 +49,7 @@ def gga_eng(self):
         if self._gga_eng is not NotImplemented:
             return self._gga_eng
 
-        grids = self.gen_grids()
+        grids = self.gen_grids(atom_grid=self.atom_grid)
 
         gga_eng = scf.RKS(self.mol)
         gga_eng.grids = grids
@@ -70,9 +71,12 @@ def gga_grad(self):
         self._gga_grad = gga_grad
         return self._gga_grad
 
-    def gen_grids(self, rad_points=99, sph_points=590):
+    def gen_grids(self, rad_points=99, sph_points=590, atom_grid=None):
         grids = dft.Grids(self.mol)
-        grids.atom_grid = (rad_points, sph_points)
+        grids.atom_grid = atom_grid
+        if atom_grid is None:
+            grids.atom_grid = (rad_points, sph_points)
+        grids.radi_method = dft.gauss_chebyshev
         grids.becke_scheme = dft.gen_grid.stratmann
         grids.build()
         return grids

pyxdh/Validation/numerical_deriv/ncdft_derivonce_uhf_ub3lyp.py

@@ -0,0 +1,54 @@
+import pickle
+from pyscf import scf
+import numpy as np
+
+from pyxdh.Utilities.test_molecules import Mol_CH3
+from pyxdh.Utilities import NumericDiff, NucCoordDerivGenerator, DipoleDerivGenerator
+from pyxdh.DerivOnce import GradUNCDFT
+
+
+def mol_to_grad_helper(mol):
+    CH3 = Mol_CH3(mol=mol)
+    CH3.hf_eng.kernel()
+    config = {
+        "scf_eng": CH3.hf_eng,
+        "nc_eng": CH3.gga_eng
+    }
+    helper = GradUNCDFT(config)
+    return helper
+
+
+def dipole_generator(component, interval):
+    CH3 = Mol_CH3()
+    mol = CH3.mol
+
+    def get_hcore(mol=mol):
+        return scf.rhf.get_hcore(mol) - interval * mol.intor("int1e_r")[component]
+
+    CH3.hf_eng.get_hcore = get_hcore
+    CH3.gga_eng.get_hcore = get_hcore
+
+    config = {
+        "scf_eng": CH3.hf_eng,
+        "nc_eng": CH3.gga_eng
+    }
+    helper = GradUNCDFT(config)
+    return helper
+
+
+if __name__ == '__main__':
+    result_dict = {}
+    CH3 = Mol_CH3()
+    mol = CH3.mol
+
+    num_obj = NucCoordDerivGenerator(CH3.mol, mol_to_grad_helper)
+    num_dif = NumericDiff(num_obj, lambda helper: helper.eng)
+    result_dict["grad"] = num_dif.derivative
+
+    num_obj = DipoleDerivGenerator(dipole_generator)
+    num_dif = NumericDiff(num_obj, lambda helper: helper.eng)
+    dip_nuc = np.einsum("A, At -> t", mol.atom_charges(), mol.atom_coords())
+    result_dict["dipole"] = num_dif.derivative + dip_nuc
+
+    with open("ncdft_derivonce_uhf_ub3lyp.dat", "wb") as f:
+        pickle.dump(result_dict, f, pickle.HIGHEST_PROTOCOL)