diff --git a/.idea/workspace.xml b/.idea/workspace.xml
index 87948b8..b6a47f7 100644
--- a/.idea/workspace.xml
+++ b/.idea/workspace.xml
@@ -19,10 +19,15 @@
-
+
+
+
+
+
+
@@ -50,7 +55,7 @@
-
+
@@ -76,7 +81,12 @@
-
+
+
+
+
+
+
@@ -99,20 +109,21 @@
-
+
+
-
-
+
+
-
+
@@ -121,7 +132,7 @@
-
+
@@ -129,12 +140,12 @@
-
+
-
+
@@ -183,6 +194,26 @@
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
@@ -238,22 +269,23 @@
-
+
@@ -287,7 +319,9 @@
-
+
+
+
1591017061437
@@ -310,7 +344,14 @@
1591613805618
-
+
+ 1591619557688
+
+
+
+ 1591619557688
+
+
@@ -332,84 +373,98 @@
-
+
+
-
+
+
-
+
-
+
+
+
+
+
+
+
-
-
+
+
+
+
-
+
-
-
-
-
+
+
+
+
-
-
-
-
+
+
+
+
-
-
-
-
+
+
+
+
-
-
-
-
-
+
+
+
+
+
-
-
-
-
+
+
+
+
+
-
-
-
-
+
+
+
+
+
-
-
-
-
+
+
+
+
+
-
-
+
+
+
@@ -419,14 +474,15 @@
-
-
+
+
-
+
+
-
+
@@ -438,20 +494,48 @@
17
+
+ file://$USER_HOME$/AppData/Local/JetBrains/PyCharm2020.1/remote_sources/-1053058360/1179071206/pyscf/hessian/rhf.py
+ 58
+
+
+
+ file://$USER_HOME$/AppData/Local/JetBrains/PyCharm2020.1/remote_sources/-1053058360/1179071206/pyscf/hessian/rhf.py
+ 576
+
+
+
+ file://$USER_HOME$/AppData/Local/JetBrains/PyCharm2020.1/remote_sources/-1053058360/1179071206/pyscf/hessian/uhf.py
+ 306
+
+
+
+ file://$USER_HOME$/AppData/Local/JetBrains/PyCharm2020.1/remote_sources/-1053058360/1179071206/pyscf/hessian/uhf.py
+ 341
+
+
+
+ file://$PROJECT_DIR$/pyxdh/DerivOnce/grad_scf.py
+ 317
+
+
-
-
-
+
+
+
+
+
+
+
-
\ No newline at end of file
diff --git a/pyxdh/DerivOnce/__init__.py b/pyxdh/DerivOnce/__init__.py
index 223825f..3f810cc 100644
--- a/pyxdh/DerivOnce/__init__.py
+++ b/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
diff --git a/pyxdh/DerivOnce/deriv_once_scf.py b/pyxdh/DerivOnce/deriv_once_scf.py
index 8d6f4cc..0b176b5 100644
--- a/pyxdh/DerivOnce/deriv_once_scf.py
+++ b/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
diff --git a/pyxdh/DerivOnce/deriv_once_uscf.py b/pyxdh/DerivOnce/deriv_once_uscf.py
index e3961c4..6b2fe1a 100644
--- a/pyxdh/DerivOnce/deriv_once_uscf.py
+++ b/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
diff --git a/pyxdh/DerivOnce/grad_uscf.py b/pyxdh/DerivOnce/grad_uscf.py
index 8b9e41e..419de96 100644
--- a/pyxdh/DerivOnce/grad_uscf.py
+++ b/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))
diff --git a/pyxdh/Utilities/deriv_numerical.py b/pyxdh/Utilities/deriv_numerical.py
index 0f718bb..16fc99c 100644
--- a/pyxdh/Utilities/deriv_numerical.py
+++ b/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
diff --git a/pyxdh/Utilities/test_molecules.py b/pyxdh/Utilities/test_molecules.py
index 76b8494..2800d79 100644
--- a/pyxdh/Utilities/test_molecules.py
+++ b/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
diff --git a/pyxdh/Validation/numerical_deriv/ncdft_derivonce_uhf_ub3lyp.dat b/pyxdh/Validation/numerical_deriv/ncdft_derivonce_uhf_ub3lyp.dat
new file mode 100644
index 0000000..a0148e6
Binary files /dev/null and b/pyxdh/Validation/numerical_deriv/ncdft_derivonce_uhf_ub3lyp.dat differ
diff --git a/pyxdh/Validation/numerical_deriv/ncdft_derivonce_uhf_ub3lyp.py b/pyxdh/Validation/numerical_deriv/ncdft_derivonce_uhf_ub3lyp.py
new file mode 100644
index 0000000..26654de
--- /dev/null
+++ b/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)