pyxdh: UB3LYP grad

ajz34 · Jun 8, 2020 · 4e5c044 · 4e5c044
1 parent 2f4d8cb
commit 4e5c044
Show file tree

Hide file tree

Showing 6 changed files with 498 additions and 33 deletions.
diff --git a/.idea/workspace.xml b/.idea/workspace.xml
diff --git a/pyxdh/DerivOnce/deriv_once_uscf.py b/pyxdh/DerivOnce/deriv_once_uscf.py
@@ -33,7 +33,7 @@ def initialization_pyscf(self):
             self.scf_eng.kernel()
             if not self.scf_eng.converged:
                 warnings.warn("SCF not converged!")
-        if isinstance(self.scf_eng, dft.rks.RKS):
+        if isinstance(self.scf_eng, dft.rks.RKS) or isinstance(self.scf_eng, dft.uks.UKS):
             self.xc = self.scf_eng.xc
             self.grids = self.scf_eng.grids
             self.xc_type = dft.libxc.xc_type(self.xc)

diff --git a/pyxdh/DerivOnce/grad_uscf.py b/pyxdh/DerivOnce/grad_uscf.py
@@ -3,10 +3,9 @@
 import os
 
 from pyscf import grad
-from pyscf.scf import _vhf
 
 from pyxdh.DerivOnce import DerivOnceUSCF, DerivOnceNCDFT, GradSCF
-from pyxdh.Utilities import GridIterator, KernelHelper, timing
+from pyxdh.Utilities import GridIterator, KernelHelper
 
 MAXMEM = float(os.getenv("MAXMEM", 2))
 np.einsum = partial(np.einsum, optimize=["greedy", 1024 ** 3 * MAXMEM / 8])
@@ -28,16 +27,32 @@ def _get_E_1(self):
         F_0_mo = self.F_0_mo
         occ = self.occ
         D = self.D
+        grids = self.grids
+        grdit_memory = self.grdit_memory
 
         E_1 = (
                 + np.einsum("Auv, xuv -> A", H_1_ao, D)
                 + 0.5 * np.einsum("Auvkl, yuv, xkl -> A", eri1_ao, D, D)
                 - 0.5 * cx * np.einsum("Aukvl, xuv, xkl -> A", eri1_ao, D, D)
                 - np.einsum("xApq, xpq, xp, xq -> A", S_1_mo, F_0_mo, occ, occ)
                 + grad.rhf.grad_nuc(mol).reshape(-1)
-        ).reshape((natm, 3))
+        )
 
-        return E_1
+        # GGA part contiribution
+        if self.xc_type == "GGA":
+            grdit = zip(GridIterator(mol, grids, D[0], deriv=2), GridIterator(mol, grids, D[1], deriv=2))
+            for grdh in grdit:
+                kerh = KernelHelper(grdh, xc)
+                E_1 += (
+                    + np.einsum("g, Atg -> At", kerh.fr[0], grdh[0].A_rho_1)
+                    + np.einsum("g, Atg -> At", kerh.fr[1], grdh[1].A_rho_1)
+                    + 2 * np.einsum("g, rg, Atrg -> At", kerh.fg[0], grdh[0].rho_1, grdh[0].A_rho_2)
+                    + 2 * np.einsum("g, rg, Atrg -> At", kerh.fg[2], grdh[1].rho_1, grdh[1].A_rho_2)
+                    + 1 * np.einsum("g, rg, Atrg -> At", kerh.fg[1], grdh[1].rho_1, grdh[0].A_rho_2)
+                    + 1 * np.einsum("g, rg, Atrg -> At", kerh.fg[1], grdh[0].rho_1, grdh[1].A_rho_2)
+            ).reshape(-1)
+
+        return E_1.reshape((natm, 3))
 
 
 class Test_GradUSCF:
@@ -52,3 +67,14 @@ def test_UHF_grad(self):
         helper = GradUSCF({"scf_eng": CH3.hf_eng})
         formchk = FormchkInterface(resource_filename("pyxdh", "Validation/gaussian/CH3-HF-freq.fchk"))
         assert(np.allclose(helper.E_1, formchk.grad(), atol=1e-5, rtol=1e-4))
+
+    def test_UB3LYP_grad(self):
+
+        from pyxdh.Utilities.test_molecules import Mol_CH3
+        from pkg_resources import resource_filename
+        from pyxdh.Utilities import FormchkInterface
+
+        CH3 = Mol_CH3()
+        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))