From 98552a31cf963b9893052ae457897c3d48304abd Mon Sep 17 00:00:00 2001 From: ajz34 Date: Sat, 6 May 2023 14:56:25 +0800 Subject: [PATCH 01/20] feat(dh resp): rmp2ri gamma --- pyscf/dh/response/mp2/rmp2ri.py | 277 ++++++++++++++++++++++++++++++++ pyscf/dh/util/df_addon.py | 11 +- 2 files changed, 287 insertions(+), 1 deletion(-) create mode 100644 pyscf/dh/response/mp2/rmp2ri.py diff --git a/pyscf/dh/response/mp2/rmp2ri.py b/pyscf/dh/response/mp2/rmp2ri.py new file mode 100644 index 0000000..7ebfbf9 --- /dev/null +++ b/pyscf/dh/response/mp2/rmp2ri.py @@ -0,0 +1,277 @@ +""" RI-MP2 Response-Related Utilities. """ + +from pyscf.dh import RMP2RI +from pyscf.dh import util +from pyscf import lib, __config__ +import h5py +import numpy as np + +CONFIG_incore_t_oovv_mp2 = getattr(__config__, "incore_t_oovv_mp2", "auto") +CONFIG_incore_cderi_uaa_mp2 = getattr(__config__, "incore_cderi_uaa_mp2", "auto") +CONFIG_max_cycle_cpks = getattr(__config__, "max_cycle_cpks", 20) +CONFIG_tol_cpks = getattr(__config__, "tol_cpks", 1e-9) + + +def get_mp2_integrals( + cderi_uov, mo_occ, mo_energy, mask_act, + c_os, c_ss, incore_t_oovv, + verbose=lib.logger.NOTE, max_memory=2000, h5file=NotImplemented): + r""" Get 3-index transformed MP2 amplitude. + + .. math:: + \Gamma_{ia, P} = T_{ij}^{ab} Y_{jb, P} + + Parameters + ---------- + cderi_uov : np.ndarray + Cholesky decomposed ERI in MO basis (occ-vir block). + dim: (naux, nocc_act, nvir_act). + mo_occ : np.ndarray + Occupation number. + dim: (nmo, ). + mo_energy : np.ndarray + Molecular orbital energy. + dim: (nmo, ). + mask_act : np.ndarray + Active molecular orbital mask. + dim: (nmo, ); dtype: bool. + + c_os : float + Oppo-spin coefficient. + c_ss : float + Same-spin coefficient. + incore_t_oovv : bool or str or float or int + Whether save :math:`t_{ij}^{ab}` into memory (True) or disk (False). + + verbose : int + Print verbosity. + max_memory : int or float + Maximum memory available for molecule object. + h5file : h5py.File + HDF5 file, used when ``t_oovv`` is stored in disk. + + Returns + ------- + dict[str, np.ndarray], dict[str, float] + + Tensors and results are returned. Tensors included: + + - `t_oovv`: MP2 amplitude :math:`t_{ij}^{ab}`. dim: (nocc_act, nocc_act, nvir_act, nvir_act). + - `G_uov`: MP2 amplitude contracted by ERI :math:`\Gamma_{ia, P}`. dim: (naux, nocc_act, nvir_act). + - `rdm1_corr`: Response density matrix of MP2 correlation contribution in MO (non-response) + :math:`D_{pq}^{(2)}`. dim: (nmo, nmo). + + Results are the MP2 energy contributions (including oppo-spin and same-spin). + """ + log = lib.logger.new_logger(verbose=verbose) + time0 = lib.logger.process_clock(), lib.logger.perf_counter() + + # dimension definition and check sanity + mask_occ_act = (mo_occ > 0) & mask_act + mask_vir_act = (mo_occ == 0) & mask_act + nocc_act = mask_occ_act.sum() + nvir_act = mask_vir_act.sum() + nmo = len(mask_act) + naux = cderi_uov.shape[0] + assert len(mo_occ) == nmo + assert (naux, nocc_act, nvir_act) == cderi_uov.shape + assert incore_t_oovv is not None, "t_oovv must be stored for evaluation of MP2 response" + + # preparation + eo = mo_energy[mask_occ_act] + ev = mo_energy[mask_vir_act] + D_ovv = lib.direct_sum("j - a - b -> jab", eo, ev, ev) + + block_oo_act = np.ix_(mask_occ_act, mask_occ_act) + block_vv_act = np.ix_(mask_vir_act, mask_vir_act) + + # allocate results + rdm1_corr = np.zeros((nmo, nmo)) + G_uov = np.zeros((naux, nocc_act, nvir_act)) + t_oovv = util.allocate_array( + incore_t_oovv, (nocc_act, nocc_act, nvir_act, nvir_act), max_memory, + h5file=h5file, name="t_oovv", zero_init=False, chunk=(1, 1, nvir_act, nvir_act)) + eng_bi1 = eng_bi2 = 0 + + # for async write t_oovv + def write_t_oovv(slc, buf): + t_oovv[slc] = buf + + # prepare batch + mem_avail = max_memory - lib.current_memory()[0] + nbatch = util.calc_batch_size(4 * nocc_act * nvir_act**2 + naux * nvir_act, mem_avail) + + with lib.call_in_background(write_t_oovv) as async_write_t_oovv: + for sI in util.gen_batch(0, nocc_act, nbatch): + log.debug(f"[DEBUG] Loop in get_mp2_integrals, slice {sI} of {nocc_act} orbitals.") + D_Oovv = lib.direct_sum("I + jab -> Ijab", eo[sI], D_ovv) + g_Oovv = lib.einsum("PIa, Pjb -> Ijab", cderi_uov[:, sI], cderi_uov) + t_Oovv = g_Oovv / D_Oovv + async_write_t_oovv(sI, t_Oovv) + T_Oovv = util.restricted_biorthogonalize(t_Oovv, 1, c_os, c_ss) + eng_bi1 += lib.einsum("Ijab, Ijab ->", t_Oovv, g_Oovv) + eng_bi2 += lib.einsum("Ijab, Ijba ->", t_Oovv, g_Oovv) + rdm1_corr[block_vv_act] += 2 * lib.einsum("ijac, ijbc -> ab", T_Oovv, t_Oovv) + rdm1_corr[block_oo_act] -= 2 * lib.einsum("ijab, ikab -> jk", T_Oovv, t_Oovv) + G_uov[:, sI] += lib.einsum("ijab, Pjb -> Pia", T_Oovv, cderi_uov) + + # results + eng_os = eng_bi1 + eng_ss = eng_bi1 - eng_bi2 + eng_mp2 = c_os * eng_os + c_ss * eng_ss + results = dict() + results["eng_corr_MP2_bi1"] = eng_bi1 + results["eng_corr_MP2_bi2"] = eng_bi2 + results["eng_corr_MP2_OS"] = eng_os + results["eng_corr_MP2_SS"] = eng_ss + results["eng_corr_MP2"] = eng_mp2 + log.note(f"[RESULT] Energy corr MP2 of same-spin: {eng_ss :18.10f}") + log.note(f"[RESULT] Energy corr MP2 of oppo-spin: {eng_os :18.10f}") + log.note(f"[RESULT] Energy corr MP2 of total: {eng_mp2:18.10f}") + + # tensors + tensors = { + "t_oovv": t_oovv, + "G_uov": G_uov, + "rdm1_corr": rdm1_corr} + + log.timer("get_mp2_integrals", *time0) + return tensors, results + + +class RMP2RespRI(RMP2RI): + + def __init__(self, *args, **kwargs): + super().__init__(*args, **kwargs) + self.incore_cderi_uaa = CONFIG_incore_cderi_uaa_mp2 + self.c_os = kwargs.get("c_os", 1) + self.c_ss = kwargs.get("c_ss", 1) + self.incore_t_oovv = kwargs.get("incore_t_oovv", CONFIG_incore_t_oovv_mp2) + + def make_cderi_uaa(self): + """ Generate cholesky decomposed ERI (all block, active only, s1 symm, in memory/disk). """ + if "cderi_uaa" in self.tensors: + return self.tensors["cderi_uaa"] + + log = lib.logger.new_logger(verbose=self.verbose) + + # dimension and mask + mask = self.get_frozen_mask() + nocc_act = (mask & (self.mo_occ != 0)).sum() + nvir_act = (mask & (self.mo_occ == 0)).sum() + nact = nocc_act + nvir_act + mo_coeff_act = self.mo_coeff[:, mask] + + # density fitting preparation + omega = self.omega + with_df = util.get_with_df_omega(self.with_df, omega) + naux = with_df.get_naoaux() + + # array allocation + mem_avail = self.max_memory - lib.current_memory()[0] + incore_cderi_uaa = self.incore_cderi_uaa + cderi_uaa = util.allocate_array( + incore_cderi_uaa, (naux, nact, nact), mem_avail, + h5file=self._tmpfile, name="cderi_uaa", zero_init=False) + log.info(f"[INFO] Store type of cderi_uaa: {type(cderi_uaa)}") + + # generate array + util.get_cderi_mo(with_df, mo_coeff_act, cderi_uaa, max_memory=self.max_memory) + + tensors = {"cderi_uaa": cderi_uaa} + self.tensors.update(tensors) + return cderi_uaa + + def make_cderi_uov(self): + """ Generate cholesky decomposed ERI (occ-vir block, active only, in memory). """ + if "cderi_uov" in self.tensors: + return self.tensors["cderi_uov"] + + # dimension and mask + mask = self.get_frozen_mask() + nocc_act = (mask & (self.mo_occ != 0)).sum() + + # In response evaluation, ERI of all orbitals are generally required. + # Thus, cderi_uaa is a must-evaluated tensor. + cderi_uaa = self.tensors.get("cderi_uaa", self.make_cderi_uaa()) + cderi_uov = cderi_uaa[:, :nocc_act, nocc_act:] + + tensors = {"cderi_uov": cderi_uov} + self.tensors.update(tensors) + return cderi_uov + + def _make_mp2_integrals(self): + cderi_uov = self.tensors.get("cderi_uov", self.make_cderi_uov()) + mo_occ = self.mo_occ + mo_energy = self.mo_energy + mask_act = self.get_frozen_mask() + c_os = self.c_os + c_ss = self.c_ss + incore_t_oovv = self.incore_t_oovv + verbose = self.verbose + max_memory = self.max_memory + h5file = self._tmpfile + tensors, results = get_mp2_integrals( + cderi_uov, mo_occ, mo_energy, mask_act, + c_os, c_ss, incore_t_oovv, + verbose=verbose, max_memory=max_memory, h5file=h5file) + self.tensors.update(tensors) + util.update_results(self.results, results) + + def make_t_oovv(self): + r""" MP2 amplitude :math:`t_{ij}^{ab}`. + + dim: (nocc_act, nocc_act, nvir_act, nvir_act) + """ + if "t_oovv" in self.tensors: + return self.tensors["t_oovv"] + + self._make_mp2_integrals() + return self.tensors["t_oovv"] + + def make_G_uov(self): + r""" MP2 amplitude contracted by ERI :math:`\Gamma_{ia, P}`. + + dim: (naux, nocc_act, nvir_act) + """ + if "G_uov" in self.tensors: + return self.tensors["G_uov"] + + self._make_mp2_integrals() + return self.tensors["G_uov"] + + def make_rdm1_corr(self): + r""" Response density matrix of MP2 correlation contribution in MO (non-response) :math:`D_{pq}^{(2)}`. + + dim: (nmo, nmo) + """ + if "rdm1_corr" in self.tensors: + return self.tensors["rdm1_corr"] + + self._make_mp2_integrals() + return self.tensors["rdm1_corr"] + + def driver_eng_mp2(self, **kwargs): + """ Driver of MP2 energy. + + For response computation of MP2, energy evaluation is accompanied with evaluation of MP2 integrals. + """ + if "eng_corr_MP2" not in self.results: + self._make_mp2_integrals() + + return self.results + + get_mp2_integrals = staticmethod(get_mp2_integrals) + + +if __name__ == '__main__': + def main_1(): + from pyscf import gto, scf + mol = gto.Mole(atom="O; H 1 0.94; H 1 0.94 2 104.5", basis="6-31G").build() + mf = scf.RHF(mol).run() + mf_mp2 = RMP2RespRI(mf) + print(mf_mp2.make_rdm1_corr()) + + main_1() + + diff --git a/pyscf/dh/util/df_addon.py b/pyscf/dh/util/df_addon.py index b05f610..f1e7879 100644 --- a/pyscf/dh/util/df_addon.py +++ b/pyscf/dh/util/df_addon.py @@ -5,7 +5,9 @@ from pyscf.ao2mo import _ao2mo from pyscf.dh.util import calc_batch_size -def get_cderi_mo(with_df, mo_coeff, Y_mo=None, pqslice=None, max_memory=2000): +def get_cderi_mo( + with_df, mo_coeff, Y_mo=None, pqslice=None, + max_memory=2000, verbose=lib.logger.NOTE): """ Get Cholesky decomposed ERI in MO basis. Parameters @@ -20,12 +22,17 @@ def get_cderi_mo(with_df, mo_coeff, Y_mo=None, pqslice=None, max_memory=2000): Slice of orbitals to be transformed. max_memory : float Maximum memory available. + verbose : int + Print verbosity. Returns ------- np.ndarray Cholesky decomposed ERI in MO basis. """ + log = lib.logger.new_logger(verbose=verbose) + time0 = lib.logger.process_clock(), lib.logger.perf_counter() + naux = with_df.get_naoaux() nmo = mo_coeff.shape[-1] if pqslice is None: @@ -55,6 +62,8 @@ def get_cderi_mo(with_df, mo_coeff, Y_mo=None, pqslice=None, max_memory=2000): Y_ao.reshape(naux, nao, nao), mo_coeff[:, pqslice[0]:pqslice[1]].conj(), mo_coeff[:, pqslice[2]:pqslice[3]]) p0 = p1 + + log.timer("get_mp2_integrals", *time0) return Y_mo -- Gitee From 85fa91bccca8415c6bd019ac439d4bd02a933043 Mon Sep 17 00:00:00 2001 From: ajz34 Date: Mon, 8 May 2023 14:37:19 +0800 Subject: [PATCH 02/20] feat(dh resp): rmp2 dipole --- pyscf/dh/energy/rhdft.py | 10 +- pyscf/dh/response/hdft/__init__.py | 1 + pyscf/dh/response/hdft/rhdft.py | 99 ++++++++ pyscf/dh/response/mp2/rmp2ri.py | 381 +++++++++++++++++++++++++++-- 4 files changed, 472 insertions(+), 19 deletions(-) create mode 100644 pyscf/dh/response/hdft/__init__.py create mode 100644 pyscf/dh/response/hdft/rhdft.py diff --git a/pyscf/dh/energy/rhdft.py b/pyscf/dh/energy/rhdft.py index 4a971e4..63c10a2 100644 --- a/pyscf/dh/energy/rhdft.py +++ b/pyscf/dh/energy/rhdft.py @@ -394,9 +394,14 @@ class RHDFT(EngBase): It's better to initialize this object first, before actually running SCF iterations. """ - def __init__(self, mf, xc): + def __init__(self, mf, xc=None): super().__init__(mf) - self.xc = xc # type: XCList + if xc is not None: + self.xc = xc + elif not hasattr(mf, "xc"): + self.xc = "HF" + else: + self.xc = self.scf.xc if isinstance(self.xc, str): xc_scf = XCList(self.xc, code_scf=True) xc_eng = XCList(self.xc, code_scf=False) @@ -405,6 +410,7 @@ class RHDFT(EngBase): self.xc = xc_scf else: xc_scf = self.xc + self._scf = custom_mf(mf, xc_scf) @property diff --git a/pyscf/dh/response/hdft/__init__.py b/pyscf/dh/response/hdft/__init__.py new file mode 100644 index 0000000..253f044 --- /dev/null +++ b/pyscf/dh/response/hdft/__init__.py @@ -0,0 +1 @@ +from pyscf.dh.response.hdft.rhdft import RHDFTResp diff --git a/pyscf/dh/response/hdft/rhdft.py b/pyscf/dh/response/hdft/rhdft.py new file mode 100644 index 0000000..9bdc9f7 --- /dev/null +++ b/pyscf/dh/response/hdft/rhdft.py @@ -0,0 +1,99 @@ +""" Hybrid-DFT Response-Related Utilities. """ + +from pyscf.dh import RHDFT +from pyscf.scf import _response_functions # this import is not unnecessary +from functools import cached_property + + +def Ax0_Core_resp(sp, sq, sr, ss, vresp, mo_coeff): + r""" Convenient function for evaluation of Fock response in MO basis + :math:`\sum_{rs} A_{pq, rs} X_{rs}^\mathbb{A}` by PySCF's response function. + + Parameters + ---------- + sp, sq, sr, ss : slice or list + Slice of molecular orbital indices. + vresp : callable + Fock response function in AO basis (generated by ``mf.scf.gen_response``). + mo_coeff : np.ndarray + Molecular orbital coefficients. + + Returns + ------- + callable + A function where input is :math:`X_{rs}^\mathbb{A}`, and output is + :math:`\sum_{rs} A_{pq, rs} X_{rs}^\mathbb{A}`. + """ + C = mo_coeff + + def Ax0_Core_resp_inner(X): + X_shape = X.shape + X = X.reshape((-1, X_shape[-2], X_shape[-1])) + dmX = C[:, sr] @ X @ C[:, ss].T + dmX += dmX.swapaxes(-1, -2) + ax_ao = vresp(dmX) + res = 2 * C[:, sp].T @ ax_ao @ C[:, sq] + res.shape = list(X_shape[:-2]) + [res.shape[-2], res.shape[-1]] + return res + + return Ax0_Core_resp_inner + + +class RHDFTResp(RHDFT): + + @cached_property + def vresp(self): + """ Fock response function (derivative w.r.t. molecular coefficient in AO basis). """ + return self.scf.gen_response() + + def Ax0_Core(self, sp, sq, sr, ss): + r""" Convenient function for evaluation of Fock response in MO basis + :math:`\sum_{rs} A_{pq, rs} X_{rs}^\mathbb{A}`. + + Parameters + ---------- + sp, sq, sr, ss : slice or list + Slice of molecular orbital indices. + + Returns + ------- + callable + A function where input is :math:`X_{rs}^\mathbb{A}`, and output is + :math:`\sum_{rs} A_{pq, rs} X_{rs}^\mathbb{A}`. + + Notes + ----- + This function acts as a wrapper of various possible Fock response algorithms. + + TODO: Functionality of this method is to be implemented. + """ + # currently we just implemented response by PySCF + return self.Ax0_Core_resp(sp, sq, sr, ss) + + def Ax0_Core_resp(self, sp, sq, sr, ss, vresp=None, mo_coeff=None): + r""" Convenient function for evaluation of Fock response in MO basis + :math:`\sum_{rs} A_{pq, rs} X_{rs}^\mathbb{A}` by PySCF's response function. + + Parameters + ---------- + sp, sq, sr, ss : slice or list + Slice of molecular orbital indices. + vresp : callable + Fock response function in AO basis (generated by ``mf.scf.gen_response``). + mo_coeff : np.ndarray + Molecular orbital coefficients. + + Returns + ------- + callable + A function where input is :math:`X_{rs}^\mathbb{A}`, and output is + :math:`\sum_{rs} A_{pq, rs} X_{rs}^\mathbb{A}`. + + Notes + ----- + This function calls PySCF's gen_response function. For cases that requires large virtual contraction + (such as :math:`A_{ai, pq} X_{pq}`), this function should be somehow quicker. + """ + vresp = vresp if vresp is not None else self.vresp + mo_coeff = mo_coeff if mo_coeff is not None else self.mo_coeff + return Ax0_Core_resp(sp, sq, sr, ss, vresp, mo_coeff) diff --git a/pyscf/dh/response/mp2/rmp2ri.py b/pyscf/dh/response/mp2/rmp2ri.py index 7ebfbf9..65a3491 100644 --- a/pyscf/dh/response/mp2/rmp2ri.py +++ b/pyscf/dh/response/mp2/rmp2ri.py @@ -2,7 +2,8 @@ from pyscf.dh import RMP2RI from pyscf.dh import util -from pyscf import lib, __config__ +from pyscf import scf, lib, __config__ +from pyscf.scf import cphf import h5py import numpy as np @@ -19,12 +20,14 @@ def get_mp2_integrals( r""" Get 3-index transformed MP2 amplitude. .. math:: - \Gamma_{ia, P} = T_{ij}^{ab} Y_{jb, P} + t_{ij}^{ab} &= Y_{ia, P} Y_{jb, P} / D_{ij}^{ab} \\ + T_{ij}^{ab} &= 2 c_\mathrm{os} t_{ij}^{ab} - c_\mathrm{ss} t_{ij}^{ba} \\ + \Gamma_{ia, P} &= T_{ij}^{ab} Y_{jb, P} \\ Parameters ---------- cderi_uov : np.ndarray - Cholesky decomposed ERI in MO basis (occ-vir block). + Cholesky decomposed ERI in MO basis (occ-vir block) :math:`Y_{ia, P}`. dim: (naux, nocc_act, nvir_act). mo_occ : np.ndarray Occupation number. @@ -139,28 +142,205 @@ def get_mp2_integrals( return tensors, results +def get_W_I(cderi_uov, cderi_uoo, G_uov, verbose=lib.logger.NOTE): + r""" Part I of MO-energy-weighted density matrix. + + .. math:: + W_{ij} [\mathrm{I}] &= -2 \Gamma_{ia, P} Y_{ja, P} \\ + W_{ab} [\mathrm{I}] &= -2 \Gamma_{ia, P} Y_{ib, P} \\ + W_{ai} [\mathrm{I}] &= -4 \Gamma_{ja, P} Y_{ij, P} + + This energy-weighted density matrix is not symmetric by current implementation. + + Parameters + ---------- + cderi_uov : np.ndarray + cderi_uoo : np.ndarray + G_uov : np.ndarray + verbose : int + + Returns + ------- + dict[str, np.ndarray] + """ + + log = lib.logger.new_logger(verbose=verbose) + time0 = lib.logger.process_clock(), lib.logger.perf_counter() + + # dimension definition and check sanity + naux, nocc, nvir = cderi_uov.shape + nmo = nocc + nvir + assert cderi_uoo.shape == (naux, nocc, nocc) + assert G_uov.shape == (naux, nocc, nvir) + + # prepare essential matrices and slices + so, sv = slice(0, nocc), slice(nocc, nmo) + + W_I = np.zeros((nmo, nmo)) + W_I[so, so] = - 2 * lib.einsum("Pia, Pja -> ij", G_uov, cderi_uov) + W_I[sv, sv] = - 2 * lib.einsum("Pia, Pib -> ab", G_uov, cderi_uov) + W_I[sv, so] = - 4 * lib.einsum("Pja, Pij -> ai", G_uov, cderi_uoo) + + log.timer("get_W_I", *time0) + tensors = {"W_I": W_I} + return tensors + + +def get_lag_vo( + G_uov, cderi_uaa, W_I, rdm1_corr, Ax0_Core, + max_memory=2000, verbose=lib.logger.NOTE): + r""" MP2 contribution to Lagrangian vir-occ block. + + .. math:: + L_{ai} = - 4 \Gamma_{ja, P} Y_{ij, P} + 4 \Gamma_{ib, P} Y_{ab, P} + A_{ai, pq} D_{pq}^\mathrm{(2)} + + where :math:`- 4 \Gamma_{ja, P} Y_{ij, P}` is evaluated by :math:`W_{ai} [\mathrm{I}]`. + + Some pratical meaning of lagrangian, is that for vir-occ block, it can be defined by MP2 energy derivative wrt + MO coefficients: + + .. math:: + \mathscr{L}_{pq} &= C_{\mu p} \frac{\partial E^\mathrm{(2)}}{\partial C_{\mu q}} \\ + L_{ai} &= \mathscr{L}_{ai} - \mathscr{L}_{ia} + + Parameters + ---------- + G_uov : np.ndarray + cderi_uaa : np.ndarray + W_I : np.ndarray + rdm1_corr : np.ndarray + Ax0_Core : callable + max_memory : float + verbose : int + + Returns + ------- + dict[str, np.ndarray] + """ + + log = lib.logger.new_logger(verbose=verbose) + time0 = lib.logger.process_clock(), lib.logger.perf_counter() + + # dimension definition and check sanity + naux, nocc, nvir = G_uov.shape + nmo = nocc + nvir + assert cderi_uaa.shape == (naux, nmo, nmo) + if W_I is not None: + assert W_I.shape == (nmo, nmo) + + # prepare essential matrices and slices + so, sv, sa = slice(0, nocc), slice(nocc, nmo), slice(0, nmo) + + # generate lagrangian occ-vir block + lag_vo = np.zeros((nvir, nocc)) + lag_vo += W_I[sv, so] + lag_vo += Ax0_Core(sv, so, sa, sa)(rdm1_corr) + mem_avail = max_memory - lib.current_memory()[0] + nbatch = util.calc_batch_size(nvir ** 2 + nocc * nvir, mem_avail) + for saux in util.gen_batch(0, naux, nbatch): + lag_vo += 4 * lib.einsum("Pib, Pab -> ai", G_uov[saux], cderi_uaa[saux, sv, sv]) + + log.timer("get_lag_vo", *time0) + tensors = {"lag_vo": lag_vo} + return tensors + + +def get_rdm1_corr_resp( + rdm1_corr, lag_vo, + mo_energy, mo_occ, Ax0_Core, + max_cycle=20, tol=1e-9, verbose=lib.logger.NOTE): + r""" 1-RDM of response MP2 correlation by MO basis. + + For other parts, response density matrix is the same to the usual 1-RDM. + + .. math:: + A'_{ai, bj} D_{bj}^\mathrm{(2)} &= L_{ai} + + Parameters + ---------- + rdm1_corr : np.ndarray + lag_vo : np.ndarray + mo_energy : np.ndarray + mo_occ : np.ndarray + Ax0_Core : callable + max_cycle : int + tol : float + verbose : int + + Returns + ------- + dict[str, np.ndarray] + + See Also + -------- + get_rdm1_corr + """ + + log = lib.logger.new_logger(verbose=verbose) + time0 = lib.logger.process_clock(), lib.logger.perf_counter() + + # dimension definition and check sanity + nocc = (mo_occ > 0).sum() + nvir = (mo_occ == 0).sum() + nmo = nocc + nvir + assert lag_vo.shape == (nvir, nocc) + assert rdm1_corr.shape == (nmo, nmo) + assert len(mo_energy) == nmo + assert len(mo_occ) == nmo + + # prepare essential matrices and slices + so, sv = slice(0, nocc), slice(nocc, nmo) + + # cp-ks evaluation + rdm1_corr_resp = rdm1_corr.copy() + rdm1_corr_resp_vo = cphf.solve( + Ax0_Core(sv, so, sv, so), mo_energy, mo_occ, lag_vo, + max_cycle=max_cycle, tol=tol)[0] + rdm1_corr_resp[sv, so] = rdm1_corr_resp_vo + + log.timer("get_rdm1_corr_resp", *time0) + tensors = {"rdm1_corr_resp": rdm1_corr_resp} + return tensors + + class RMP2RespRI(RMP2RI): def __init__(self, *args, **kwargs): super().__init__(*args, **kwargs) - self.incore_cderi_uaa = CONFIG_incore_cderi_uaa_mp2 self.c_os = kwargs.get("c_os", 1) self.c_ss = kwargs.get("c_ss", 1) + self.incore_cderi_uaa = CONFIG_incore_cderi_uaa_mp2 self.incore_t_oovv = kwargs.get("incore_t_oovv", CONFIG_incore_t_oovv_mp2) + self.max_cycle_cpks = CONFIG_max_cycle_cpks + self.tol_cpks = CONFIG_tol_cpks + self._Ax0_Core = NotImplemented + + @property + def mask_occ(self): + return self.mo_occ != 0 + + @property + def mask_occ_act(self): + return self.get_frozen_mask() & (self.mo_occ != 0) + + @property + def mask_vir(self): + return self.mo_occ == 0 + + @property + def mask_vir_act(self): + return self.get_frozen_mask() & (self.mo_occ == 0) def make_cderi_uaa(self): - """ Generate cholesky decomposed ERI (all block, active only, s1 symm, in memory/disk). """ + """ Generate cholesky decomposed ERI (all block, full orbitals, s1 symm, in memory/disk). """ if "cderi_uaa" in self.tensors: return self.tensors["cderi_uaa"] log = lib.logger.new_logger(verbose=self.verbose) # dimension and mask - mask = self.get_frozen_mask() - nocc_act = (mask & (self.mo_occ != 0)).sum() - nvir_act = (mask & (self.mo_occ == 0)).sum() - nact = nocc_act + nvir_act - mo_coeff_act = self.mo_coeff[:, mask] + mo_coeff = self.mo_coeff + nmo = self.nmo # density fitting preparation omega = self.omega @@ -171,12 +351,12 @@ class RMP2RespRI(RMP2RI): mem_avail = self.max_memory - lib.current_memory()[0] incore_cderi_uaa = self.incore_cderi_uaa cderi_uaa = util.allocate_array( - incore_cderi_uaa, (naux, nact, nact), mem_avail, + incore_cderi_uaa, (naux, nmo, nmo), mem_avail, h5file=self._tmpfile, name="cderi_uaa", zero_init=False) log.info(f"[INFO] Store type of cderi_uaa: {type(cderi_uaa)}") # generate array - util.get_cderi_mo(with_df, mo_coeff_act, cderi_uaa, max_memory=self.max_memory) + util.get_cderi_mo(with_df, mo_coeff, cderi_uaa, max_memory=self.max_memory) tensors = {"cderi_uaa": cderi_uaa} self.tensors.update(tensors) @@ -188,18 +368,33 @@ class RMP2RespRI(RMP2RI): return self.tensors["cderi_uov"] # dimension and mask - mask = self.get_frozen_mask() - nocc_act = (mask & (self.mo_occ != 0)).sum() + mask_occ_act = self.mask_occ_act + mask_vir_act = self.mask_vir_act # In response evaluation, ERI of all orbitals are generally required. # Thus, cderi_uaa is a must-evaluated tensor. + # note that cderi_uaa may be stored by h5py, so indexing by list should be done in following way cderi_uaa = self.tensors.get("cderi_uaa", self.make_cderi_uaa()) - cderi_uov = cderi_uaa[:, :nocc_act, nocc_act:] + cderi_uov = cderi_uaa[:, mask_occ_act, :][:, :, mask_vir_act] tensors = {"cderi_uov": cderi_uov} self.tensors.update(tensors) return cderi_uov + def make_cderi_uoo(self): + """ Generate cholesky decomposed ERI (occ-occ block, in memory). """ + # todo: frozen core not applied + if "cderi_uoo" in self.tensors: + return self.tensors["cderi_uoo"] + + mo_coeff_occ = self.mo_coeff[:, self.mask_occ] + with_df = self.with_df + cderi_uoo = util.get_cderi_mo(with_df, mo_coeff_occ, max_memory=self.max_memory) + + tensors = {"cderi_uoo": cderi_uoo} + self.tensors.update(tensors) + return cderi_uoo + def _make_mp2_integrals(self): cderi_uov = self.tensors.get("cderi_uov", self.make_cderi_uov()) mo_occ = self.mo_occ @@ -261,6 +456,124 @@ class RMP2RespRI(RMP2RI): return self.results + def make_W_I(self): + r""" Part I of MO-energy-weighted density matrix W_{pq} [\mathrm{I}]. + + dim: (nmo, nmo) + """ + if "W_I" in self.tensors: + return self.tensors["W_I"] + + # prepare tensors + cderi_uov = self.tensors.get("cderi_uov", self.make_cderi_uov()) + cderi_uoo = self.tensors.get("cderi_uoo", self.make_cderi_uoo()) + G_uov = self.tensors.get("G_uov", self.make_G_uov()) + + # main function + tensors = get_W_I(cderi_uov, cderi_uoo, G_uov, verbose=self.verbose) + self.tensors.update(tensors) + return tensors["W_I"] + + def make_lag_vo(self): + r""" Generate MP2 contribution to Lagrangian vir-occ block :math:`L_{ai}`. """ + if "lag_vo" in self.tensors: + return self.tensors["lag_vo"] + + # prepare tensors + W_I = self.tensors.get("W_I", self.make_W_I()) + cderi_uaa = self.tensors.get("cderi_uaa", self.make_cderi_uaa()) + G_uov = self.tensors["G_uov"] + rdm1_corr = self.tensors["rdm1_corr"] + Ax0_Core = self.Ax0_Core + + # main function + tensors = get_lag_vo( + G_uov, cderi_uaa, W_I, rdm1_corr, Ax0_Core, + max_memory=self.max_memory, verbose=self.verbose) + + self.tensors.update(tensors) + return tensors["lag_vo"] + + def make_rdm1_corr_resp(self): + """ Generate 1-RDM (response) of MP2 correlation :math:`D_{pq}^{(2)}`. """ + if "rdm1_corr_resp" in self.tensors: + return self.tensors["rdm1_corr_resp"] + + # prepare input + lag_vo = self.tensors.get("lag_vo", self.make_lag_vo()) + rdm1_corr = self.tensors["rdm1_corr"] + mask = self.get_frozen_mask() + mo_energy = self.mo_energy[mask] + mo_occ = self.mo_occ[mask] + Ax0_Core = self.Ax0_Core + max_cycle = self.max_cycle_cpks + tol = self.tol_cpks + verbose = self.verbose + + # main function + tensors = get_rdm1_corr_resp( + rdm1_corr, lag_vo, mo_energy, mo_occ, Ax0_Core, + max_cycle=max_cycle, tol=tol, verbose=verbose) + + self.tensors.update(tensors) + return tensors["rdm1_corr_resp"] + + def make_rdm1(self, ao=False): + r""" Generate 1-RDM (non-response) of MP2 :math:`D_{pq}^{MP2}` in MO or :math:`D_{\mu \nu}^{MP2}` in AO. """ + # prepare input + rdm1_corr = self.tensors.get("rdm1_corr", self.make_rdm1_corr()) + + rdm1 = np.diag(self.mo_occ) + mask = self.get_frozen_mask() + ix_act = np.ix_(mask, mask) + rdm1[ix_act] += rdm1_corr + self.tensors["rdm1"] = rdm1 + if ao: + rdm1 = self.mo_coeff @ rdm1 @ self.mo_coeff.T + return rdm1 + + def make_rdm1_resp(self, ao=False): + r""" Generate 1-RDM (response) of MP2 :math:`D_{pq}^{MP2}` in MO or :math:`D_{\mu \nu}^{MP2}` in AO. """ + # prepare input + rdm1_corr_resp = self.tensors.get("rdm1_corr_resp", self.make_rdm1_corr_resp()) + + rdm1 = np.diag(self.mo_occ) + mask = self.get_frozen_mask() + ix_act = np.ix_(mask, mask) + rdm1[ix_act] += rdm1_corr_resp + self.tensors["rdm1_resp"] = rdm1 + if ao: + rdm1 = self.mo_coeff @ rdm1 @ self.mo_coeff.T + return rdm1 + + def make_dipole(self): + # prepare input + mol = self.mol + rdm1_ao = self.make_rdm1_resp(ao=True) + int1e_r = mol.intor("int1e_r") + + dip_elec = - np.einsum("uv, tuv -> t", rdm1_ao, int1e_r) + dip_nuc = np.einsum("A, At -> t", mol.atom_charges(), mol.atom_coords()) + dip = dip_elec + dip_nuc + self.tensors["dipole"] = dip + return dip + + @property + def Ax0_Core(self): + """ Fock response of underlying SCF object in MO basis. """ + if self._Ax0_Core is NotImplemented: + restricted = isinstance(self.scf, scf.rhf.RHF) + from pyscf.dh.response.hdft import RHDFTResp + UHDFTResp = NotImplemented + HDFTResp = RHDFTResp if restricted else UHDFTResp + mf_scf = HDFTResp(self.scf) + self._Ax0_Core = mf_scf.Ax0_Core + return self._Ax0_Core + + @Ax0_Core.setter + def Ax0_Core(self, Ax0_Core): + self._Ax0_Core = Ax0_Core + get_mp2_integrals = staticmethod(get_mp2_integrals) @@ -270,8 +583,42 @@ if __name__ == '__main__': mol = gto.Mole(atom="O; H 1 0.94; H 1 0.94 2 104.5", basis="6-31G").build() mf = scf.RHF(mol).run() mf_mp2 = RMP2RespRI(mf) - print(mf_mp2.make_rdm1_corr()) + print(mf_mp2.make_dipole()) - main_1() + def main_2(): + # test numerical dipole + from pyscf import gto, scf, mp, df + np.set_printoptions(5, suppress=True, linewidth=150) + mol = gto.Mole(atom="O; H 1 0.94; H 1 0.94 2 104.5", basis="6-31G").build() + + def eng_with_dipole_field(t, h): + mf_scf = scf.RHF(mol).density_fit("cc-pVDZ-jkfit") + mf_scf.get_hcore = lambda *args, **kwargs: scf.hf.get_hcore(mol) - h * mol.intor("int1e_r")[t] + mf_scf.run(conv_tol=1e-12) + mf_mp = RMP2RI(mf_scf).run() + # mf_mp = mp.dfmp2.DFMP2(mf_scf) + # mf_mp.with_df = df.DF(mol, df.make_auxbasis(mol, mp2fit=True)) + # mf_mp.run() + return mf_mp.e_tot + + mf_scf = scf.RHF(mol).density_fit("cc-pVDZ-jkfit").run() + mf_mp = RMP2RespRI(mf_scf).run() + dip_anal = mf_mp.make_dipole() + + eng_array = np.zeros((2, 3)) + h = 1e-4 + for idx, h in [(0, h), [1, -h]]: + for t in (0, 1, 2): + eng_array[idx, t] = eng_with_dipole_field(t, h) + dip_elec_num = - (eng_array[0] - eng_array[1]) / (2 * h) + dip_nuc = np.einsum("A, At -> t", mol.atom_charges(), mol.atom_coords()) + dip_num = dip_elec_num + dip_nuc + print(mf_scf.dip_moment(unit="AU")) + print(dip_num) + print(dip_anal) + main_2() + + + main_1() -- Gitee From c1a148201b49987cc7328bb6e24a12248862ceb8 Mon Sep 17 00:00:00 2001 From: ajz34 Date: Mon, 8 May 2023 15:54:19 +0800 Subject: [PATCH 03/20] feat(dh resp): add get_eri_cpks --- pyscf/dh/response/hdft/rhdft.py | 72 +++++++++++++++++++++++++++++++++ 1 file changed, 72 insertions(+) diff --git a/pyscf/dh/response/hdft/rhdft.py b/pyscf/dh/response/hdft/rhdft.py index 9bdc9f7..883d093 100644 --- a/pyscf/dh/response/hdft/rhdft.py +++ b/pyscf/dh/response/hdft/rhdft.py @@ -1,10 +1,82 @@ """ Hybrid-DFT Response-Related Utilities. """ from pyscf.dh import RHDFT +from pyscf.dh import util +from pyscf import lib, __config__ from pyscf.scf import _response_functions # this import is not unnecessary +import numpy as np from functools import cached_property +CONFIG_dh_verbose = getattr(__config__, "dh_verbose", lib.logger.NOTE) + + +def get_eri_cpks( + cderi_uaa, mo_occ, cx, eri_cpks=None, + max_memory=2000, verbose=CONFIG_dh_verbose): + r""" Obtain ERI term :math:`A_{ai, bj}^\mathrm{HF}`. + + For given + + Parameters + ---------- + cderi_uaa : np.ndarray + mo_occ : np.ndarray + cx : float + eri_cpks : np.ndarray or None + max_memory : float + verbose : int + + Returns + ------- + np.ndarray + """ + log = lib.logger.new_logger(verbose=verbose) + time0 = lib.logger.process_clock(), lib.logger.perf_counter() + + # dimension and sanity check + naux = cderi_uaa.shape[0] + nmo = len(mo_occ) + nocc = (mo_occ != 0).sum() + nvir = (mo_occ == 0).sum() + so, sv = slice(0, nocc), slice(nocc, nmo) + assert cderi_uaa.shape == (naux, nmo, nmo) + + # allocate eri_cpks if necessary + if eri_cpks is None: + eri_cpks = np.zeros((nvir, nocc, nvir, nocc)) + else: + assert eri_cpks.shape == (nvir, nocc, nvir, nocc) + + # copy some tensors to memory + cderi_uvo = np.asarray(cderi_uaa[:, sv, so]) + cderi_uoo = np.asarray(cderi_uaa[:, so, so]) + + # prepare for async and batch + mem_avail = max_memory - lib.current_memory()[0] + nbatch = util.calc_batch_size(nvir*naux + 2*nocc**2*nvir, mem_avail) + + def load_cderi_uaa(slc): + return np.asarray(cderi_uaa[:, slc, :]) + + def write_eri_cpks(slc, blk): + eri_cpks[slc] += blk + + batches = util.gen_batch(nocc, nmo, nbatch) + with lib.call_in_background(write_eri_cpks) as async_write_eri_cpks: + for sA, cderi_uAa in zip(batches, lib.map_with_prefetch(load_cderi_uaa, batches)): + log.debug(f"[DEBUG] in get_eri_cpks, slice {sA} of virtual orbitals {nocc, nmo}") + sAvir = slice(sA.start - nocc, sA.stop - nocc) + blk = ( + + 4 * lib.einsum("Pai, Pbj -> aibj", cderi_uvo[:, sAvir], cderi_uvo) + - cx * lib.einsum("Paj, Pbi -> aibj", cderi_uvo[:, sAvir], cderi_uvo) + - cx * lib.einsum("Pij, Pab -> aibj", cderi_uoo, cderi_uAa)) + async_write_eri_cpks(sAvir, blk) + + log.timer("get_eri_cpks", *time0) + return eri_cpks + + def Ax0_Core_resp(sp, sq, sr, ss, vresp, mo_coeff): r""" Convenient function for evaluation of Fock response in MO basis :math:`\sum_{rs} A_{pq, rs} X_{rs}^\mathbb{A}` by PySCF's response function. -- Gitee From 68ee4f30530373cf695cb1f18135381a43ceef5e Mon Sep 17 00:00:00 2001 From: ajz34 Date: Mon, 8 May 2023 15:56:17 +0800 Subject: [PATCH 04/20] style(dh resp): add verbose to Ax0_Core_resp --- pyscf/dh/response/hdft/rhdft.py | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/pyscf/dh/response/hdft/rhdft.py b/pyscf/dh/response/hdft/rhdft.py index 883d093..514b89d 100644 --- a/pyscf/dh/response/hdft/rhdft.py +++ b/pyscf/dh/response/hdft/rhdft.py @@ -77,7 +77,7 @@ def get_eri_cpks( return eri_cpks -def Ax0_Core_resp(sp, sq, sr, ss, vresp, mo_coeff): +def Ax0_Core_resp(sp, sq, sr, ss, vresp, mo_coeff, verbose=CONFIG_dh_verbose): r""" Convenient function for evaluation of Fock response in MO basis :math:`\sum_{rs} A_{pq, rs} X_{rs}^\mathbb{A}` by PySCF's response function. @@ -89,6 +89,8 @@ def Ax0_Core_resp(sp, sq, sr, ss, vresp, mo_coeff): Fock response function in AO basis (generated by ``mf.scf.gen_response``). mo_coeff : np.ndarray Molecular orbital coefficients. + verbose : int + Print verbosity. Returns ------- @@ -99,6 +101,9 @@ def Ax0_Core_resp(sp, sq, sr, ss, vresp, mo_coeff): C = mo_coeff def Ax0_Core_resp_inner(X): + log = lib.logger.new_logger(verbose=verbose) + time0 = lib.logger.process_clock(), lib.logger.perf_counter() + X_shape = X.shape X = X.reshape((-1, X_shape[-2], X_shape[-1])) dmX = C[:, sr] @ X @ C[:, ss].T @@ -106,6 +111,8 @@ def Ax0_Core_resp(sp, sq, sr, ss, vresp, mo_coeff): ax_ao = vresp(dmX) res = 2 * C[:, sp].T @ ax_ao @ C[:, sq] res.shape = list(X_shape[:-2]) + [res.shape[-2], res.shape[-1]] + + log.timer("Ax0_Core_resp_inner", *time0) return res return Ax0_Core_resp_inner -- Gitee From e4b0d87c1fad247eb0b730282892196a9727d09d Mon Sep 17 00:00:00 2001 From: ajz34 Date: Mon, 8 May 2023 23:05:04 +0800 Subject: [PATCH 05/20] style(dh resp): valid that rsh and meta-GGA is supported for Ax0_cpks --- pyscf/dh/energy/rresp.py | 98 --------- pyscf/dh/response/hdft/rhdft.py | 349 ++++++++++++++++++++++++++++++-- 2 files changed, 327 insertions(+), 120 deletions(-) delete mode 100644 pyscf/dh/energy/rresp.py diff --git a/pyscf/dh/energy/rresp.py b/pyscf/dh/energy/rresp.py deleted file mode 100644 index d6ba04e..0000000 --- a/pyscf/dh/energy/rresp.py +++ /dev/null @@ -1,98 +0,0 @@ -import numpy as np -from pyscf import lib - -from pyscf.dh import util - - -def kernel_resp_eri(params, Y_mo, nocc, hyb, - Y_mo_lr=None, alpha=None, - max_memory=2000, verbose=lib.logger.NOTE): - """ Prepare eri evluated in response function - - Cholesky decomposed ERIs in this function should be jk form - (density-fitting object used in SCF object). - - Parameters - ---------- - params : util.Params - Y_mo : np.ndarray - nocc : int - hyb : float - Y_mo_lr : np.ndarray - alpha : float - max_memory : float - verbose : int - - Returns - ------- - np.ndarray - """ - log = lib.logger.new_logger(verbose=verbose) - naux, nmo, _ = Y_mo.shape - nvir = nmo - nocc - is_rsh = Y_mo_lr is not None - einsum = lib.einsum - so, sv = slice(0, nocc), slice(nocc, nmo) - # sanity check - assert Y_mo.dtype is not np.complex - assert Y_mo.shape == (naux, nmo, nmo) - if is_rsh: - assert Y_mo_lr.shape == (naux, nmo, nmo) - assert alpha is not None - # create space bulk - incore_resp_eri = params.flags["incore_resp_eri"] - incore_resp_eri = util.parse_incore_flag( - incore_resp_eri, nocc**2 * nvir**2, max_memory) - log.debug("Actual incore strategy of incore_resp_eri: {:}".format(incore_resp_eri)) - log.debug("Creating `resp_eri`, shape {:}.".format((nvir, nocc, nvir, nocc))) - resp_eri = params.tensors.create("resp_eri", shape=(nvir, nocc, nvir, nocc), chunk=(1, 1, nvir, nocc)) - # take some parts of cholesky eri into memory - Y_vo = np.asarray(Y_mo[:, sv, so]) - Y_oo = np.asarray(Y_mo[:, so, so]) - Y_vo_lr = Y_oo_lr = None - if is_rsh: - Y_vo_lr = np.asarray(Y_mo_lr[:, sv, so]) - Y_oo_lr = np.asarray(Y_mo_lr[:, so, so]) - - nbatch = util.calc_batch_size( - nvir * naux + 3 * nocc ** 2 * nvir, max_memory, 2 * Y_mo.size) - - def save_eri_cpks(sA, buf): - resp_eri[sA] = buf - - with lib.call_in_background(save_eri_cpks) as async_write: - for sA in util.gen_batch(nocc, nmo, nbatch): - sAvir = slice(sA.start - nocc, sA.stop - nocc) - print(sAvir) - resp_eri_buf = ( - + 4 * einsum("Pai, Pbj -> aibj", Y_vo[:, sAvir], Y_vo) - - hyb * einsum("Paj, Pbi -> aibj", Y_vo[:, sAvir], Y_vo) - - hyb * einsum("Pij, Pab -> aibj", Y_oo, Y_mo[:, sA, sv])) - if is_rsh: - resp_eri_buf += ( - - (alpha - hyb) * einsum("Paj, Pbi -> aibj", Y_vo_lr[:, sAvir], Y_vo_lr) - - (alpha - hyb) * einsum("Pij, Pab -> aibj", Y_oo_lr, Y_mo_lr[:, sA, sv])) - async_write(sAvir, resp_eri_buf) - - return dict() - - -def resp_Ax0_HF(si, sa, sj, sb, hyb, Y_mo, - alpha=None, Y_mo_lr=None, - max_memory=2000, verbose=lib.logger.NOTE): - naux, nmo, _ = Y_mo.shape - ni, na = si.stop - si.start, sa.stop - sa.start - # sanity check - is_rsh = Y_mo_lr is not None - - def resp_Ax0_HF_inner(X): - X_shape = X.shape - X = X.reshape((-1, X_shape[-2], X_shape[-1])) - res = np.zeros((X.shape[0], ni, na)) - nbatch = util.calc_batch_size(nmo**2, max_memory, X.size + res.size) - for saux in util.gen_batch(0, naux, nbatch): - pass - - - - diff --git a/pyscf/dh/response/hdft/rhdft.py b/pyscf/dh/response/hdft/rhdft.py index 514b89d..95f7d12 100644 --- a/pyscf/dh/response/hdft/rhdft.py +++ b/pyscf/dh/response/hdft/rhdft.py @@ -2,17 +2,21 @@ from pyscf.dh import RHDFT from pyscf.dh import util -from pyscf import lib, __config__ +from pyscf import gto, dft, lib, __config__ from pyscf.scf import _response_functions # this import is not unnecessary +from pyscf.dh.energy.rhdft import get_rho import numpy as np from functools import cached_property CONFIG_dh_verbose = getattr(__config__, "dh_verbose", lib.logger.NOTE) +CONFIG_incore_cderi_uaa_hdft = getattr(__config__, "incore_cderi_uaa_hdft", "auto") +CONFIG_incore_eri_cpks_vovo = getattr(__config__, "incore_eri_cpks_vovo", "auto") +CONFIG_dft_gen_grid_Grids_level = getattr(__config__, 'dft_gen_grid_Grids_level', 3) -def get_eri_cpks( - cderi_uaa, mo_occ, cx, eri_cpks=None, +def get_eri_cpks_vovo( + cderi_uaa, mo_occ, cj, ck, eri_cpks_vovo=None, max_memory=2000, verbose=CONFIG_dh_verbose): r""" Obtain ERI term :math:`A_{ai, bj}^\mathrm{HF}`. @@ -22,8 +26,9 @@ def get_eri_cpks( ---------- cderi_uaa : np.ndarray mo_occ : np.ndarray - cx : float - eri_cpks : np.ndarray or None + cj : float + ck : float + eri_cpks_vovo : np.ndarray or None max_memory : float verbose : int @@ -42,11 +47,11 @@ def get_eri_cpks( so, sv = slice(0, nocc), slice(nocc, nmo) assert cderi_uaa.shape == (naux, nmo, nmo) - # allocate eri_cpks if necessary - if eri_cpks is None: - eri_cpks = np.zeros((nvir, nocc, nvir, nocc)) + # allocate eri_cpks_vovo if necessary + if eri_cpks_vovo is None: + eri_cpks_vovo = np.zeros((nvir, nocc, nvir, nocc)) else: - assert eri_cpks.shape == (nvir, nocc, nvir, nocc) + assert eri_cpks_vovo.shape == (nvir, nocc, nvir, nocc) # copy some tensors to memory cderi_uvo = np.asarray(cderi_uaa[:, sv, so]) @@ -57,27 +62,117 @@ def get_eri_cpks( nbatch = util.calc_batch_size(nvir*naux + 2*nocc**2*nvir, mem_avail) def load_cderi_uaa(slc): - return np.asarray(cderi_uaa[:, slc, :]) + return np.asarray(cderi_uaa[:, slc, sv]) - def write_eri_cpks(slc, blk): - eri_cpks[slc] += blk + def write_eri_cpks_vovo(slc, blk): + eri_cpks_vovo[slc] += blk batches = util.gen_batch(nocc, nmo, nbatch) - with lib.call_in_background(write_eri_cpks) as async_write_eri_cpks: + with lib.call_in_background(write_eri_cpks_vovo) as async_write_eri_cpks_vovo: for sA, cderi_uAa in zip(batches, lib.map_with_prefetch(load_cderi_uaa, batches)): - log.debug(f"[DEBUG] in get_eri_cpks, slice {sA} of virtual orbitals {nocc, nmo}") + log.debug(f"[DEBUG] in get_eri_cpks_vovo, slice {sA} of virtual orbitals {nocc, nmo}") sAvir = slice(sA.start - nocc, sA.stop - nocc) - blk = ( - + 4 * lib.einsum("Pai, Pbj -> aibj", cderi_uvo[:, sAvir], cderi_uvo) - - cx * lib.einsum("Paj, Pbi -> aibj", cderi_uvo[:, sAvir], cderi_uvo) - - cx * lib.einsum("Pij, Pab -> aibj", cderi_uoo, cderi_uAa)) - async_write_eri_cpks(sAvir, blk) + blk = np.zeros((sA.stop - sA.start, nocc, nvir, nocc)) + if abs(cj) > 1e-10: + blk += 4 * cj * lib.einsum("Pai, Pbj -> aibj", cderi_uvo[:, sAvir], cderi_uvo) + if abs(ck) > 1e-10: + blk -= ck * ( + + lib.einsum("Paj, Pbi -> aibj", cderi_uvo[:, sAvir], cderi_uvo) + + lib.einsum("Pij, Pab -> aibj", cderi_uoo, cderi_uAa)) + async_write_eri_cpks_vovo(sAvir, blk) + + log.timer("get_eri_cpks_vovo", *time0) + return eri_cpks_vovo + + +def Ax0_cpks_HF(eri_cpks_vovo, max_memory=2000, verbose=CONFIG_dh_verbose): + r""" Convenient function for evaluation of HF contribution of Fock response in MO basis + :math:`\sum_{rs} A_{ai, bj} X_{bj}^\mathbb{A}` by explicitly contraction to MO ERI :math:`(ai|bj)`. + + Parameters + ---------- + eri_cpks_vovo : np.ndarray + max_memory : float + verbose : int + + Returns + ------- + callable + """ + nvir, nocc = eri_cpks_vovo.shape[:2] + + def Ax0_cpks_HF_inner(X): + log = lib.logger.new_logger(verbose=verbose) + time0 = lib.logger.process_clock(), lib.logger.perf_counter() + + X_shape = X.shape + X = X.reshape((-1, X_shape[-2], X_shape[-1])) + res = np.zeros_like(X) + + def load_eri_cpks_vovo(slc): + return eri_cpks_vovo[slc] + + mem_avail = max_memory - lib.current_memory()[0] + nbatch = util.calc_batch_size(nocc**2 * nvir, mem_avail) + batches = util.gen_batch(0, nvir, nbatch) + + for sA, eri_cpks_Vovo in zip(batches, lib.map_with_prefetch(load_eri_cpks_vovo, batches)): + res[:, sA] = lib.einsum("aibj, Abj -> Aai", eri_cpks_Vovo, X) + + res.shape = list(X_shape[:-2]) + [res.shape[-2], res.shape[-1]] + + log.timer("Ax0_cpks_HF_inner", *time0) + return res + + return Ax0_cpks_HF_inner - log.timer("get_eri_cpks", *time0) - return eri_cpks +def Ax0_Core_KS( + sp, sq, sr, ss, + mo_coeff, xc_setting, xc_kernel, + max_memory=2000, verbose=CONFIG_dh_verbose): + r""" Convenient function for evaluation of pure DFT contribution of Fock response in MO basis + :math:`\sum_{rs} A_{pq, rs} X_{rs}^\mathbb{A}`. + + Parameters + ---------- + sp, sq, sr, ss : slice or list + mo_coeff : np.ndarray + xc_setting : tuple[dft.numint.NumInt, gto.Mole, dft.Grids, str, np.ndarray] + xc_kernel : tuple[np.ndarray, np.ndarray, np.ndarray] + max_memory : float + verbose : int -def Ax0_Core_resp(sp, sq, sr, ss, vresp, mo_coeff, verbose=CONFIG_dh_verbose): + Returns + ------- + callable + """ + C = mo_coeff + ni, mol, grids, xc, dm = xc_setting + rho, _, fxc = xc_kernel + + def Ax0_Core_KS_inner(X): + log = lib.logger.new_logger(verbose=verbose) + time0 = lib.logger.process_clock(), lib.logger.perf_counter() + mem_avail = max_memory - lib.current_memory()[0] + + X_shape = X.shape + X = X.reshape((-1, X_shape[-2], X_shape[-1])) + dmX = C[:, sr] @ X @ C[:, ss].T + dmX += dmX.swapaxes(-1, -2) + ax_ao = ni.nr_rks_fxc(mol, grids, xc, dm, dmX, hermi=1, rho0=rho, fxc=fxc, max_memory=mem_avail) + res = 2 * C[:, sp].T @ ax_ao @ C[:, sq] + res.shape = list(X_shape[:-2]) + [res.shape[-2], res.shape[-1]] + + log.timer("Ax0_Core_KS_inner", *time0) + return res + + return Ax0_Core_KS_inner + + +def Ax0_Core_resp( + sp, sq, sr, ss, vresp, mo_coeff, + verbose=CONFIG_dh_verbose): r""" Convenient function for evaluation of Fock response in MO basis :math:`\sum_{rs} A_{pq, rs} X_{rs}^\mathbb{A}` by PySCF's response function. @@ -118,13 +213,75 @@ def Ax0_Core_resp(sp, sq, sr, ss, vresp, mo_coeff, verbose=CONFIG_dh_verbose): return Ax0_Core_resp_inner +def get_xc_integral(ni, mol, grids, xc, dm): + rho = get_rho(mol, grids, dm) + try: + exc, vxc, fxc, kxc = ni.eval_xc_eff(xc, rho, deriv=3) + tensors = { + "rho": rho, + f"exc_{xc}": exc, + f"vxc_{xc}": vxc, + f"fxc_{xc}": fxc, + f"kxc_{xc}": kxc, + } + except NotImplementedError: + exc, vxc, fxc, _ = ni.eval_xc_eff(xc, rho, deriv=2) + tensors = { + "rho": rho, + f"exc_{xc}": exc, + f"vxc_{xc}": vxc, + f"fxc_{xc}": fxc, + } + return tensors + + class RHDFTResp(RHDFT): + def __init__(self, *args, **kwargs): + super().__init__(*args, **kwargs) + self.incore_cderi_uaa = CONFIG_incore_cderi_uaa_hdft + self.incore_eri_cpks_vovo = CONFIG_incore_eri_cpks_vovo + + grid_level_cpks = max(CONFIG_dft_gen_grid_Grids_level - 1, 1) + self.grids_cpks = dft.Grids(self.mol) + self.grids_cpks.level = grid_level_cpks + self.grids_cpks.build() + @cached_property def vresp(self): """ Fock response function (derivative w.r.t. molecular coefficient in AO basis). """ return self.scf.gen_response() + def make_cderi_uaa(self, omega=0): + """ Generate cholesky decomposed ERI (all block, full orbitals, s1 symm, in memory/disk). """ + if util.pad_omega("cderi_uaa", omega) in self.tensors: + return self.tensors[util.pad_omega("cderi_uaa", omega)] + + log = lib.logger.new_logger(verbose=self.verbose) + + # dimension and mask + mo_coeff = self.mo_coeff + nmo = self.nmo + + # density fitting preparation + with_df = util.get_with_df_omega(self.scf.with_df, omega) + naux = with_df.get_naoaux() + + # array allocation + mem_avail = self.max_memory - lib.current_memory()[0] + incore_cderi_uaa = self.incore_cderi_uaa + cderi_uaa = util.allocate_array( + incore_cderi_uaa, (naux, nmo, nmo), mem_avail, + h5file=self._tmpfile, name=util.pad_omega("cderi_uaa", omega), zero_init=False) + log.info(f"[INFO] Store type of cderi_uaa: {type(cderi_uaa)}") + + # generate array + util.get_cderi_mo(with_df, mo_coeff, cderi_uaa, max_memory=self.max_memory) + + tensors = {util.pad_omega("cderi_uaa", omega): cderi_uaa} + self.tensors.update(tensors) + return cderi_uaa + def Ax0_Core(self, sp, sq, sr, ss): r""" Convenient function for evaluation of Fock response in MO basis :math:`\sum_{rs} A_{pq, rs} X_{rs}^\mathbb{A}`. @@ -176,3 +333,151 @@ class RHDFTResp(RHDFT): vresp = vresp if vresp is not None else self.vresp mo_coeff = mo_coeff if mo_coeff is not None else self.mo_coeff return Ax0_Core_resp(sp, sq, sr, ss, vresp, mo_coeff) + + def make_eri_cpks_vovo(self): + r""" Generate ERI for CP-KS evaluation :math:`(ai, bj)` for current exch-corr setting. """ + if "eri_cpks_vovo" in self.tensors: + return self.tensors["eri_cpks_vovo"] + + # This function will utilize mf.scf._numint to specify hybrid and omega coefficients + + ni = self.scf._numint # type: dft.numint.NumInt + omega, alpha, hyb = ni.rsh_and_hybrid_coeff(self.scf.xc) + nvir, nocc = self.nvir, self.nocc + cderi_uaa = self.tensors.get("cderi_uaa", self.make_cderi_uaa()) + + eri_cpks_vovo = util.allocate_array( + incore=self.incore_eri_cpks_vovo, + shape=(nvir, nocc, nvir, nocc), + max_memory=self.max_memory, + h5file=self._tmpfile, + name="eri_cpks_vovo", + chunk=(1, nocc, nvir, nocc)) + + get_eri_cpks_vovo( + cderi_uaa=cderi_uaa, + mo_occ=self.mo_occ, + cj=1, ck=hyb, + eri_cpks_vovo=eri_cpks_vovo, + max_memory=self.max_memory, + verbose=self.verbose + ) + + if abs(omega) > 1e-10: + cderi_uaa = self.make_cderi_uaa(omega=omega) + get_eri_cpks_vovo( + cderi_uaa=cderi_uaa, + mo_occ=self.mo_occ, + cj=0, ck=alpha-hyb, + eri_cpks_vovo=eri_cpks_vovo, + max_memory=self.max_memory, + verbose=self.verbose + ) + + self.tensors["eri_cpks_vovo"] = eri_cpks_vovo + return eri_cpks_vovo + + def Ax0_cpks_HF(self): + eri_cpks_vovo = self.tensors.get("eri_cpks_vovo", self.make_eri_cpks_vovo()) + ax0_cpks_hf = Ax0_cpks_HF(eri_cpks_vovo, self.max_memory, self.verbose) + return ax0_cpks_hf + + def make_xc_integral(self): + ni = self.scf._numint # type: dft.numint.NumInt + mol = self.mol + grids = self.scf.grids + xc_token = self.xc.token + dm = self.scf.make_rdm1() + tensors = get_xc_integral(ni, mol, grids, xc_token, dm) + self.tensors.update(tensors) + return tensors + + def make_xc_integral_cpks(self): + ni = self.scf._numint # type: dft.numint.NumInt + mol = self.mol + grids = self.grids_cpks + xc_token = self.xc.token + dm = self.scf.make_rdm1() + tensors = get_xc_integral(ni, mol, grids, xc_token, dm) + tensors = {key + "/cpks": val for key, val in tensors.items()} + self.tensors.update(tensors) + return tensors + + def Ax0_Core_KS(self, sp, sq, sr, ss): + ni = self.scf._numint + mol = self.mol + grids = self.grids_cpks + xc_token = self.xc.token + dm = self.scf.make_rdm1() + if f"fxc_{xc_token}" not in self.tensors: + self.make_xc_integral_cpks() + rho = self.tensors[f"rho/cpks"] + vxc = self.tensors[f"vxc_{xc_token}/cpks"] + fxc = self.tensors[f"fxc_{xc_token}/cpks"] + + xc_setting = ni, mol, grids, xc_token, dm + xc_kernel = rho, vxc, fxc + mo_coeff = self.mo_coeff + + ax0_core_ks = Ax0_Core_KS( + sp, sq, sr, ss, + mo_coeff, xc_setting, xc_kernel, + max_memory=self.max_memory, + verbose=self.verbose) + return ax0_core_ks + + def Ax0_cpks(self): + nocc, nmo = self.nocc, self.nmo + so, sv = slice(0, nocc), slice(nocc, nmo) + ax0_core_ks = self.Ax0_Core_KS(sv, so, sv, so) + ax0_cpks_hf = self.Ax0_cpks_HF() + + def Ax0_cpks_inner(X): + res = ax0_cpks_hf(X) + ax0_core_ks(X) + return res + return Ax0_cpks_inner + + +if __name__ == '__main__': + def main_1(): + # test that RSH functional is correct for Ax0_cpks + from pyscf import gto, scf + np.set_printoptions(5, suppress=True) + np.random.seed(0) + + mol = gto.Mole(atom="O; H 1 0.94; H 1 0.94 2 104.5", basis="6-31G").build() + mf = dft.RKS(mol, xc="CAM-B3LYP").density_fit().run() + mf_scf = RHDFTResp(mf) + mf_scf.grids_cpks = mf_scf.scf.grids + + nocc, nvir, nmo = mf_scf.nocc, mf_scf.nvir, mf_scf.nmo + so, sv = slice(0, nocc), slice(nocc, nmo) + X = np.random.randn(3, nvir, nocc) + ax_cpks = mf_scf.Ax0_cpks()(X) + ax_core = mf_scf.Ax0_Core(sv, so, sv, so)(X) + print(np.allclose(ax_cpks, ax_core)) + # print(ax_core[0]) + # print(ax_cpks[0]) + + def main_2(): + # test that meta-GGA functional is correct for Ax0_cpks + from pyscf import gto, scf + np.set_printoptions(5, suppress=True) + np.random.seed(0) + + mol = gto.Mole(atom="O; H 1 0.94; H 1 0.94 2 104.5", basis="6-31G").build() + mf = dft.RKS(mol, xc="TPSS0").density_fit().run() + mf_scf = RHDFTResp(mf) + mf_scf.grids_cpks = mf_scf.scf.grids + + nocc, nvir, nmo = mf_scf.nocc, mf_scf.nvir, mf_scf.nmo + so, sv = slice(0, nocc), slice(nocc, nmo) + X = np.random.randn(3, nvir, nocc) + ax_cpks = mf_scf.Ax0_cpks()(X) + ax_core = mf_scf.Ax0_Core(sv, so, sv, so)(X) + print(np.allclose(ax_cpks, ax_core)) + # print(ax_core[0]) + # print(ax_cpks[0]) + + main_1() + main_2() -- Gitee From b3d32c4fb4b8df1dae85a885f57929331ea39a1d Mon Sep 17 00:00:00 2001 From: ajz34 Date: Tue, 9 May 2023 08:19:29 +0800 Subject: [PATCH 06/20] test(dh resp): add test for ax0_cpks --- pyscf/dh/response/hdft/rhdft.py | 29 +++++++++--- pyscf/dh/response/hdft/test/test_rhdft.py | 58 +++++++++++++++++++++++ 2 files changed, 81 insertions(+), 6 deletions(-) create mode 100644 pyscf/dh/response/hdft/test/test_rhdft.py diff --git a/pyscf/dh/response/hdft/rhdft.py b/pyscf/dh/response/hdft/rhdft.py index 95f7d12..60c00a9 100644 --- a/pyscf/dh/response/hdft/rhdft.py +++ b/pyscf/dh/response/hdft/rhdft.py @@ -13,6 +13,7 @@ CONFIG_dh_verbose = getattr(__config__, "dh_verbose", lib.logger.NOTE) CONFIG_incore_cderi_uaa_hdft = getattr(__config__, "incore_cderi_uaa_hdft", "auto") CONFIG_incore_eri_cpks_vovo = getattr(__config__, "incore_eri_cpks_vovo", "auto") CONFIG_dft_gen_grid_Grids_level = getattr(__config__, 'dft_gen_grid_Grids_level', 3) +CONFIG_use_eri_cpks = getattr(__config__, "use_eri_cpks", True) def get_eri_cpks_vovo( @@ -241,6 +242,7 @@ class RHDFTResp(RHDFT): super().__init__(*args, **kwargs) self.incore_cderi_uaa = CONFIG_incore_cderi_uaa_hdft self.incore_eri_cpks_vovo = CONFIG_incore_eri_cpks_vovo + self.use_eri_cpks = CONFIG_use_eri_cpks grid_level_cpks = max(CONFIG_dft_gen_grid_Grids_level - 1, 1) self.grids_cpks = dft.Grids(self.mol) @@ -300,11 +302,26 @@ class RHDFTResp(RHDFT): Notes ----- This function acts as a wrapper of various possible Fock response algorithms. - - TODO: Functionality of this method is to be implemented. """ - # currently we just implemented response by PySCF - return self.Ax0_Core_resp(sp, sq, sr, ss) + # if not RI, then use general Ax0_Core_resp + if not hasattr(self.scf, "with_df") or not self.use_eri_cpks: + return self.Ax0_Core_resp(sp, sq, sr, ss) + + # try if satisfies CPKS evaluation (vovo) + lst_nmo = np.arange(self.nmo) + lst_occ = lst_nmo[:self.nocc] + lst_vir = lst_nmo[self.nocc:] + try: + if ( + np.all(lst_nmo[sp] == lst_vir) and np.all(lst_nmo[sq] == lst_occ) and + np.all(lst_nmo[sr] == lst_vir) and np.all(lst_nmo[ss] == lst_occ)): + return self.Ax0_cpks() + else: + # otherwise, use response by PySCF + return self.Ax0_Core_resp(sp, sq, sr, ss) + except ValueError: + # dimension not match + return self.Ax0_Core_resp(sp, sq, sr, ss) def Ax0_Core_resp(self, sp, sq, sr, ss, vresp=None, mo_coeff=None): r""" Convenient function for evaluation of Fock response in MO basis @@ -454,7 +471,7 @@ if __name__ == '__main__': so, sv = slice(0, nocc), slice(nocc, nmo) X = np.random.randn(3, nvir, nocc) ax_cpks = mf_scf.Ax0_cpks()(X) - ax_core = mf_scf.Ax0_Core(sv, so, sv, so)(X) + ax_core = mf_scf.Ax0_Core_resp(sv, so, sv, so)(X) print(np.allclose(ax_cpks, ax_core)) # print(ax_core[0]) # print(ax_cpks[0]) @@ -474,7 +491,7 @@ if __name__ == '__main__': so, sv = slice(0, nocc), slice(nocc, nmo) X = np.random.randn(3, nvir, nocc) ax_cpks = mf_scf.Ax0_cpks()(X) - ax_core = mf_scf.Ax0_Core(sv, so, sv, so)(X) + ax_core = mf_scf.Ax0_Core_resp(sv, so, sv, so)(X) print(np.allclose(ax_cpks, ax_core)) # print(ax_core[0]) # print(ax_cpks[0]) diff --git a/pyscf/dh/response/hdft/test/test_rhdft.py b/pyscf/dh/response/hdft/test/test_rhdft.py new file mode 100644 index 0000000..ad48566 --- /dev/null +++ b/pyscf/dh/response/hdft/test/test_rhdft.py @@ -0,0 +1,58 @@ +import unittest +from pyscf.dh.response.hdft.rhdft import RHDFTResp +from pyscf import gto, dft +import numpy as np + + +class TestRHDFTResp(unittest.TestCase): + + def test_rsh_ax0_cpks(self): + # test that RSH functional is correct for Ax0_cpks + np.random.seed(0) + + mol = gto.Mole(atom="O; H 1 0.94; H 1 0.94 2 104.5", basis="6-31G").build() + mf = dft.RKS(mol, xc="CAM-B3LYP").density_fit().run() + mf_scf = RHDFTResp(mf) + mf_scf.grids_cpks = mf_scf.scf.grids + + nocc, nvir, nmo = mf_scf.nocc, mf_scf.nvir, mf_scf.nmo + so, sv = slice(0, nocc), slice(nocc, nmo) + X = np.random.randn(3, nvir, nocc) + ax_cpks = mf_scf.Ax0_cpks()(X) + ax_core = mf_scf.Ax0_Core_resp(sv, so, sv, so)(X) + self.assertTrue("eri_cpks_vovo" in mf_scf.tensors) + self.assertTrue(np.allclose(ax_cpks, ax_core)) + + def test_meta_ax0_cpks(self): + # test that meta-GGA functional is correct for Ax0_cpks + np.random.seed(0) + + mol = gto.Mole(atom="O; H 1 0.94; H 1 0.94 2 104.5", basis="6-31G").build() + mf = dft.RKS(mol, xc="TPSS0").density_fit().run() + mf_scf = RHDFTResp(mf) + mf_scf.grids_cpks = mf_scf.scf.grids + + nocc, nvir, nmo = mf_scf.nocc, mf_scf.nvir, mf_scf.nmo + so, sv = slice(0, nocc), slice(nocc, nmo) + X = np.random.randn(3, nvir, nocc) + ax_cpks = mf_scf.Ax0_Core(sv, so, sv, so)(X) + ax_core = mf_scf.Ax0_Core_resp(sv, so, sv, so)(X) + self.assertTrue("eri_cpks_vovo" in mf_scf.tensors) + self.assertTrue(np.allclose(ax_cpks, ax_core)) + + def test_no_density_fit(self): + # test that RSH functional is correct for Ax0_cpks if no density fitting + np.random.seed(0) + + mol = gto.Mole(atom="O; H 1 0.94; H 1 0.94 2 104.5", basis="6-31G").build() + mf = dft.RKS(mol, xc="CAM-B3LYP").run() + mf_scf = RHDFTResp(mf) + mf_scf.grids_cpks = mf_scf.scf.grids + + nocc, nvir, nmo = mf_scf.nocc, mf_scf.nvir, mf_scf.nmo + so, sv = slice(0, nocc), slice(nocc, nmo) + X = np.random.randn(3, nvir, nocc) + ax_cpks = mf_scf.Ax0_Core(sv, so, sv, so)(X) + ax_core = mf_scf.Ax0_Core_resp(sv, so, sv, so)(X) + self.assertTrue(np.allclose(ax_cpks, ax_core)) + -- Gitee From 5b2265fe2d1857c1869a40ca848bb80ac8505fca Mon Sep 17 00:00:00 2001 From: ajz34 Date: Tue, 9 May 2023 08:41:10 +0800 Subject: [PATCH 07/20] test(dh resp): add test for rmp2ri dipole --- pyscf/dh/response/mp2/rmp2ri.py | 35 +++++++++++-- pyscf/dh/response/mp2/test/test_rmp2ri.py | 62 +++++++++++++++++++++++ 2 files changed, 94 insertions(+), 3 deletions(-) create mode 100644 pyscf/dh/response/mp2/test/test_rmp2ri.py diff --git a/pyscf/dh/response/mp2/rmp2ri.py b/pyscf/dh/response/mp2/rmp2ri.py index 65a3491..90abbdd 100644 --- a/pyscf/dh/response/mp2/rmp2ri.py +++ b/pyscf/dh/response/mp2/rmp2ri.py @@ -588,7 +588,7 @@ if __name__ == '__main__': def main_2(): # test numerical dipole from pyscf import gto, scf, mp, df - np.set_printoptions(5, suppress=True, linewidth=150) + np.set_printoptions(8, suppress=True, linewidth=150) mol = gto.Mole(atom="O; H 1 0.94; H 1 0.94 2 104.5", basis="6-31G").build() def eng_with_dipole_field(t, h): @@ -618,7 +618,36 @@ if __name__ == '__main__': print(dip_num) print(dip_anal) - main_2() + def main_3(): + # test numerical dipole for DFT PT2 + from pyscf import gto, scf, dft + np.set_printoptions(8, suppress=True, linewidth=150) + mol = gto.Mole(atom="O; H 1 0.94; H 1 0.94 2 104.5", basis="6-31G").build() + + c_os, c_ss = 1.3, 0.6 + def eng_with_dipole_field(t, h): + mf_scf = dft.RKS(mol, xc="B3LYPg").density_fit("cc-pVDZ-jkfit") + mf_scf.get_hcore = lambda *args, **kwargs: scf.hf.get_hcore(mol) - h * mol.intor("int1e_r")[t] + mf_scf.run(conv_tol=1e-12) + mf_mp = RMP2RI(mf_scf).run() + return mf_mp.scf.e_tot + c_os * mf_mp.results["eng_corr_MP2_OS"] + c_ss * mf_mp.results["eng_corr_MP2_SS"] + + mf_scf = dft.RKS(mol, xc="B3LYPg").density_fit("cc-pVDZ-jkfit").run() + mf_mp = RMP2RespRI(mf_scf).run(c_os=c_os, c_ss=c_ss) + dip_anal = mf_mp.make_dipole() + + eng_array = np.zeros((2, 3)) + h = 1e-4 + for idx, h in [(0, h), [1, -h]]: + for t in (0, 1, 2): + eng_array[idx, t] = eng_with_dipole_field(t, h) + dip_elec_num = - (eng_array[0] - eng_array[1]) / (2 * h) + dip_nuc = np.einsum("A, At -> t", mol.atom_charges(), mol.atom_coords()) + dip_num = dip_elec_num + dip_nuc + + print(mf_scf.dip_moment(unit="AU")) + print(dip_num) + print(dip_anal) - main_1() + main_3() diff --git a/pyscf/dh/response/mp2/test/test_rmp2ri.py b/pyscf/dh/response/mp2/test/test_rmp2ri.py new file mode 100644 index 0000000..f9fbce0 --- /dev/null +++ b/pyscf/dh/response/mp2/test/test_rmp2ri.py @@ -0,0 +1,62 @@ +import unittest +from pyscf.dh.response.mp2.rmp2ri import RMP2RespRI +from pyscf.dh.energy.mp2.rmp2 import RMP2RI +import numpy as np + + +class TestRMP2RI(unittest.TestCase): + def test_num_dipole_ref_hf(self): + from pyscf import gto, scf + mol = gto.Mole(atom="O; H 1 0.94; H 1 0.94 2 104.5", basis="6-31G").build() + + def eng_with_dipole_field(t, h): + mf_scf = scf.RHF(mol).density_fit("cc-pVDZ-jkfit") + mf_scf.get_hcore = lambda *args, **kwargs: scf.hf.get_hcore(mol) - h * mol.intor("int1e_r")[t] + mf_scf.run(conv_tol=1e-12) + mf_mp = RMP2RI(mf_scf).run() + return mf_mp.e_tot + + mf_scf = scf.RHF(mol).density_fit("cc-pVDZ-jkfit").run() + mf_mp = RMP2RespRI(mf_scf).run() + dip_anal = mf_mp.make_dipole() + + eng_array = np.zeros((2, 3)) + h = 1e-4 + for idx, h in [(0, h), [1, -h]]: + for t in (0, 1, 2): + eng_array[idx, t] = eng_with_dipole_field(t, h) + dip_elec_num = - (eng_array[0] - eng_array[1]) / (2 * h) + dip_nuc = np.einsum("A, At -> t", mol.atom_charges(), mol.atom_coords()) + dip_num = dip_elec_num + dip_nuc + + self.assertTrue(np.allclose(dip_num, dip_anal, atol=1e-5, rtol=1e-7)) + + def test_num_dipole_ref_b3lyp(self): + from pyscf import gto, scf, dft + mol = gto.Mole(atom="O; H 1 0.94; H 1 0.94 2 104.5", basis="6-31G").build() + + c_os, c_ss = 1.3, 0.6 + + def eng_with_dipole_field(t, h): + mf_scf = dft.RKS(mol, xc="B3LYPg").density_fit("cc-pVDZ-jkfit") + mf_scf.get_hcore = lambda *args, **kwargs: scf.hf.get_hcore(mol) - h * mol.intor("int1e_r")[t] + mf_scf.run(conv_tol=1e-12) + mf_mp = RMP2RI(mf_scf).run() + return mf_mp.scf.e_tot + c_os * mf_mp.results["eng_corr_MP2_OS"] + c_ss * mf_mp.results["eng_corr_MP2_SS"] + + mf_scf = dft.RKS(mol, xc="B3LYPg").density_fit("cc-pVDZ-jkfit").run() + mf_mp = RMP2RespRI(mf_scf).run(c_os=c_os, c_ss=c_ss) + dip_anal = mf_mp.make_dipole() + + eng_array = np.zeros((2, 3)) + h = 1e-4 + for idx, h in [(0, h), [1, -h]]: + for t in (0, 1, 2): + eng_array[idx, t] = eng_with_dipole_field(t, h) + dip_elec_num = - (eng_array[0] - eng_array[1]) / (2 * h) + dip_nuc = np.einsum("A, At -> t", mol.atom_charges(), mol.atom_coords()) + dip_num = dip_elec_num + dip_nuc + + self.assertTrue(np.allclose(dip_num, dip_anal, atol=1e-5, rtol=1e-7)) + + -- Gitee From 1d819b144a44e87595e456798ffcd02473990159 Mon Sep 17 00:00:00 2001 From: ajz34 Date: Tue, 9 May 2023 09:02:47 +0800 Subject: [PATCH 08/20] refactor(dh eng): move hdft to new folder --- pyscf/dh/energy/__init__.py | 6 ++---- pyscf/dh/energy/dh.py | 4 ++-- pyscf/dh/energy/hdft/__init__.py | 4 ++++ pyscf/dh/energy/{ => hdft}/rhdft.py | 0 pyscf/dh/energy/{ => hdft}/uhdft.py | 3 +-- 5 files changed, 9 insertions(+), 8 deletions(-) create mode 100644 pyscf/dh/energy/hdft/__init__.py rename pyscf/dh/energy/{ => hdft}/rhdft.py (100%) rename pyscf/dh/energy/{ => hdft}/uhdft.py (97%) diff --git a/pyscf/dh/energy/__init__.py b/pyscf/dh/energy/__init__.py index 4edca96..c7d8654 100644 --- a/pyscf/dh/energy/__init__.py +++ b/pyscf/dh/energy/__init__.py @@ -1,13 +1,11 @@ from .dhbase import EngBase -from .rhdft import RHDFT -from .uhdft import UHDFT - -from . import mp2, iepa, ring_ccd +from . import mp2, iepa, ring_ccd, hdft from .mp2 import ( RMP2RIPySCF, RMP2ConvPySCF, UMP2ConvPySCF, RMP2Conv, RMP2RI, UMP2Conv, UMP2RI) from .iepa import (RIEPAConv, RIEPARI, UIEPAConv, UIEPARI) from .ring_ccd import (RRingCCDConv, URingCCDConv) +from .hdft import RHDFT, UHDFT from .dh import DH diff --git a/pyscf/dh/energy/dh.py b/pyscf/dh/energy/dh.py index 6025e65..8e39c35 100644 --- a/pyscf/dh/energy/dh.py +++ b/pyscf/dh/energy/dh.py @@ -2,9 +2,9 @@ from pyscf.dh.energy import EngBase from typing import Tuple, List from pyscf.dh import util from pyscf.dh.util import XCType, XCList, XCDH, update_results, pad_omega -from pyscf.dh.energy.rhdft import get_rho, numint_customized +from pyscf.dh.energy.hdft.rhdft import get_rho, numint_customized from pyscf import lib, scf, gto, dft, df, __config__ -from pyscf.dh.energy.rhdft import custom_mf +from pyscf.dh.energy.hdft.rhdft import custom_mf CONFIG_etb_first = getattr(__config__, "etb_first", False) diff --git a/pyscf/dh/energy/hdft/__init__.py b/pyscf/dh/energy/hdft/__init__.py new file mode 100644 index 0000000..93bffda --- /dev/null +++ b/pyscf/dh/energy/hdft/__init__.py @@ -0,0 +1,4 @@ +from . import rhdft, uhdft + +from .rhdft import RHDFT +from .uhdft import UHDFT diff --git a/pyscf/dh/energy/rhdft.py b/pyscf/dh/energy/hdft/rhdft.py similarity index 100% rename from pyscf/dh/energy/rhdft.py rename to pyscf/dh/energy/hdft/rhdft.py diff --git a/pyscf/dh/energy/uhdft.py b/pyscf/dh/energy/hdft/uhdft.py similarity index 97% rename from pyscf/dh/energy/uhdft.py rename to pyscf/dh/energy/hdft/uhdft.py index ecab719..039a5c2 100644 --- a/pyscf/dh/energy/uhdft.py +++ b/pyscf/dh/energy/hdft/uhdft.py @@ -1,6 +1,5 @@ -from pyscf.dh.energy.rhdft import RHDFT +from pyscf.dh.energy.hdft.rhdft import RHDFT from pyscf.dh import util -from pyscf import dft import numpy as np -- Gitee From bd234a77a018abf5caec03338976e7aa76325660 Mon Sep 17 00:00:00 2001 From: ajz34 Date: Tue, 9 May 2023 09:20:29 +0800 Subject: [PATCH 09/20] fix(dh resp): update kernel function --- pyscf/dh/response/mp2/rmp2ri.py | 1 + 1 file changed, 1 insertion(+) diff --git a/pyscf/dh/response/mp2/rmp2ri.py b/pyscf/dh/response/mp2/rmp2ri.py index 90abbdd..5e58876 100644 --- a/pyscf/dh/response/mp2/rmp2ri.py +++ b/pyscf/dh/response/mp2/rmp2ri.py @@ -574,6 +574,7 @@ class RMP2RespRI(RMP2RI): def Ax0_Core(self, Ax0_Core): self._Ax0_Core = Ax0_Core + kernel = driver_eng_mp2 get_mp2_integrals = staticmethod(get_mp2_integrals) -- Gitee From b4501c680820c2368ed58e452594f7cde73bd890 Mon Sep 17 00:00:00 2001 From: ajz34 Date: Tue, 9 May 2023 10:54:47 +0800 Subject: [PATCH 10/20] fix(dh resp): add option of xc to DH.to_scf --- pyscf/dh/energy/dh.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/pyscf/dh/energy/dh.py b/pyscf/dh/energy/dh.py index 8e39c35..715dc8e 100644 --- a/pyscf/dh/energy/dh.py +++ b/pyscf/dh/energy/dh.py @@ -614,15 +614,17 @@ class DH(EngBase): def e_tot(self): return self.results[f"eng_dh_{self.xc.xc_eng.token}"] - def to_scf(self, **kwargs): + def to_scf(self, xc=None, **kwargs): # import if self.restricted: from pyscf.dh import RHDFT as HDFT else: from pyscf.dh import UHDFT as HDFT + xc = xc if xc is not None else self.xc.xc_scf + # generate instance - mf = HDFT.from_rdh(self, self.scf, self.xc.xc_scf, **kwargs) + mf = HDFT.from_rdh(self, self.scf, xc, **kwargs) return mf -- Gitee From cd5a79e080dcbcc32b4ef664839c8e03a9d6cee1 Mon Sep 17 00:00:00 2001 From: ajz34 Date: Tue, 9 May 2023 10:58:03 +0800 Subject: [PATCH 11/20] fix(dh resp): fix refactor of energy.hdft --- pyscf/dh/response/hdft/rhdft.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyscf/dh/response/hdft/rhdft.py b/pyscf/dh/response/hdft/rhdft.py index 60c00a9..799b664 100644 --- a/pyscf/dh/response/hdft/rhdft.py +++ b/pyscf/dh/response/hdft/rhdft.py @@ -4,7 +4,7 @@ from pyscf.dh import RHDFT from pyscf.dh import util from pyscf import gto, dft, lib, __config__ from pyscf.scf import _response_functions # this import is not unnecessary -from pyscf.dh.energy.rhdft import get_rho +from pyscf.dh.energy.hdft.rhdft import get_rho import numpy as np from functools import cached_property -- Gitee From f6e69e17c4b917bf93a48bf3c6f7d6c07d1cbdbc Mon Sep 17 00:00:00 2001 From: ajz34 Date: Tue, 9 May 2023 11:29:17 +0800 Subject: [PATCH 12/20] fix(dh eng): copy method instance instead of use mf directly in customize_mf --- pyscf/dh/energy/hdft/rhdft.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/pyscf/dh/energy/hdft/rhdft.py b/pyscf/dh/energy/hdft/rhdft.py index 4a971e4..d3f1b94 100644 --- a/pyscf/dh/energy/hdft/rhdft.py +++ b/pyscf/dh/energy/hdft/rhdft.py @@ -3,6 +3,7 @@ from pyscf.dh.util import XCList, XCType, XCInfo from pyscf.dh.energy import EngBase from pyscf import dft, lib, scf, df, __config__ import numpy as np +import copy from types import MethodType CONFIG_ssr_x_fr = getattr(__config__, "ssr_x_fr", "LDA_X") @@ -333,6 +334,7 @@ def custom_mf(mf, xc, auxbasis_or_with_df=None): Note that if no with_df object passed in, the density-fitting setting of an SCF object is left as is. So leaving option ``auxbasis_or_with_df=None`` will not convert density-fitting SCF to conventional SCF. """ + mf = copy.copy(mf) verbose = mf.verbose log = lib.logger.new_logger(verbose=verbose) restricted = isinstance(mf, scf.hf.RHF) -- Gitee From 088c5e642e39d40a1c059ca73f127bfcabce3970 Mon Sep 17 00:00:00 2001 From: ajz34 Date: Tue, 9 May 2023 13:14:59 +0800 Subject: [PATCH 13/20] fix(dh eng): change name and logic of from_cls --- pyscf/dh/energy/dh.py | 12 ++++++------ pyscf/dh/energy/dhbase.py | 12 ++++++------ 2 files changed, 12 insertions(+), 12 deletions(-) diff --git a/pyscf/dh/energy/dh.py b/pyscf/dh/energy/dh.py index 715dc8e..dbe6761 100644 --- a/pyscf/dh/energy/dh.py +++ b/pyscf/dh/energy/dh.py @@ -624,7 +624,7 @@ class DH(EngBase): xc = xc if xc is not None else self.xc.xc_scf # generate instance - mf = HDFT.from_rdh(self, self.scf, xc, **kwargs) + mf = HDFT.from_cls(self, self.scf, xc, **kwargs) return mf @@ -643,9 +643,9 @@ class DH(EngBase): # generate instance if route_mp2.lower().startswith("ri"): - mf = MP2RI.from_rdh(self, self.scf, **kwargs) + mf = MP2RI.from_cls(self, self.scf, **kwargs) elif route_mp2.lower().startswith("conv"): - mf = MP2Conv.from_rdh(self, self.scf, **kwargs) + mf = MP2Conv.from_cls(self, self.scf, **kwargs) else: assert False, "Not recognized route_mp2." @@ -672,9 +672,9 @@ class DH(EngBase): # generate instance if route_iepa.lower().startswith("ri"): - mf = IEPARI.from_rdh(self, self.scf, **kwargs) + mf = IEPARI.from_cls(self, self.scf, **kwargs) elif route_iepa.lower().startswith("conv"): - mf = IEPAConv.from_rdh(self, self.scf, **kwargs) + mf = IEPAConv.from_cls(self, self.scf, **kwargs) else: assert False, "Not recognized route_iepa." @@ -701,7 +701,7 @@ class DH(EngBase): max_cycle_ring_ccd = self.flags.get("max_cycle_ring_ccd", NotImplemented) # generated instance - mf = RingCCDConv.from_rdh(self, self.scf, **kwargs) + mf = RingCCDConv.from_cls(self, self.scf, **kwargs) # fill configurations if tol_eng_ring_ccd is not NotImplemented: diff --git a/pyscf/dh/energy/dhbase.py b/pyscf/dh/energy/dhbase.py index db81883..7b7d082 100644 --- a/pyscf/dh/energy/dhbase.py +++ b/pyscf/dh/energy/dhbase.py @@ -149,7 +149,7 @@ class EngBase(lib.StreamObject, ABC): return frozen_core.mask @classmethod - def from_rdh(cls, mf_dh: "EngBase", *args, **kwargs): + def from_cls(cls, mf_dh: "EngBase", *args, **kwargs): """ Build a new (inherited) instance by a base EngBase instance. Only some attributes are copied. Excluded attributes including @@ -159,9 +159,9 @@ class EngBase(lib.StreamObject, ABC): - ``e_corr`` """ mf_new = cls(*args, **kwargs) - mf_new.with_df = mf_dh.with_df - mf_new.frozen = mf_dh.frozen - mf_new.verbose = mf_dh.verbose - mf_new.max_memory = mf_dh.max_memory - mf_new.base = mf_dh + mf_new.__dict__.update(mf_dh.__dict__) + mf_new.results = dict() + mf_new.tensors = dict() + mf_new._tmpfile = lib.H5TmpFile() + mf_new.e_corr = NotImplemented return mf_new -- Gitee From 9c2797b0bc8a92e2e7bf7e9b0ca5c5042fda380a Mon Sep 17 00:00:00 2001 From: ajz34 Date: Tue, 9 May 2023 16:34:38 +0800 Subject: [PATCH 14/20] fix(dh eng): accept default xc for HDFT class --- pyscf/dh/energy/hdft/rhdft.py | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/pyscf/dh/energy/hdft/rhdft.py b/pyscf/dh/energy/hdft/rhdft.py index d3f1b94..726f4cc 100644 --- a/pyscf/dh/energy/hdft/rhdft.py +++ b/pyscf/dh/energy/hdft/rhdft.py @@ -396,9 +396,14 @@ class RHDFT(EngBase): It's better to initialize this object first, before actually running SCF iterations. """ - def __init__(self, mf, xc): + def __init__(self, mf, xc=None): super().__init__(mf) - self.xc = xc # type: XCList + if xc is not None: + self.xc = xc + elif not hasattr(mf, "xc"): + self.xc = "HF" + else: + self.xc = self.scf.xc if isinstance(self.xc, str): xc_scf = XCList(self.xc, code_scf=True) xc_eng = XCList(self.xc, code_scf=False) @@ -407,6 +412,7 @@ class RHDFT(EngBase): self.xc = xc_scf else: xc_scf = self.xc + self._scf = custom_mf(mf, xc_scf) @property -- Gitee From 8339513cd242951a3edcfebf8046c22c12cdedc4 Mon Sep 17 00:00:00 2001 From: ajz34 Date: Tue, 9 May 2023 16:39:10 +0800 Subject: [PATCH 15/20] fix(dh eng): new attribute `hdft` make a new attribute `hdft` for storing energy evaluation of hybrid DFT instead of using attribute `_scf` --- pyscf/dh/energy/hdft/rhdft.py | 20 ++++++++++++-------- 1 file changed, 12 insertions(+), 8 deletions(-) diff --git a/pyscf/dh/energy/hdft/rhdft.py b/pyscf/dh/energy/hdft/rhdft.py index 726f4cc..2cda438 100644 --- a/pyscf/dh/energy/hdft/rhdft.py +++ b/pyscf/dh/energy/hdft/rhdft.py @@ -1,3 +1,5 @@ +from functools import cached_property + from pyscf.dh import util from pyscf.dh.util import XCList, XCType, XCInfo from pyscf.dh.energy import EngBase @@ -413,15 +415,15 @@ class RHDFT(EngBase): else: xc_scf = self.xc - self._scf = custom_mf(mf, xc_scf) + self.hdft = custom_mf(mf, xc_scf) - @property + @cached_property def restricted(self): - return isinstance(self.scf, scf.rhf.RHF) + return isinstance(self.hdft, scf.rhf.RHF) - @property + @cached_property def e_tot(self) -> float: - return self.scf.e_tot + return self.hdft.energy_tot() def make_energy_purexc(self, xc_lists, numint=None, dm=None): """ Evaluate energy contributions of pure (DFT) exchange-correlation effects. @@ -439,9 +441,9 @@ class RHDFT(EngBase): -------- get_energy_purexc """ - grids = self.scf.grids + grids = self.hdft.grids if dm is None: - dm = self.scf.make_rdm1() + dm = self.hdft.make_rdm1() dm = np.asarray(dm) if (self.restricted and dm.ndim != 2) or (not self.restricted and (dm.ndim != 3 or dm.shape[0] != 2)): raise ValueError("Dimension of input density matrix is not correct.") @@ -450,7 +452,9 @@ class RHDFT(EngBase): xc_lists, rho, grids.weights, self.restricted, numint=numint) def kernel(self, *args, **kwargs): - return self.scf.kernel(*args, **kwargs) + if not self.scf.converged: + self.scf.kernel(*args, **kwargs) + return self.e_tot get_energy_exactx = staticmethod(get_energy_restricted_exactx) get_energy_noxc = staticmethod(get_energy_restricted_noxc) -- Gitee From 2121d6336e11ecbbc1035cf89bd1ed1853cda5eb Mon Sep 17 00:00:00 2001 From: ajz34 Date: Wed, 10 May 2023 13:29:41 +0800 Subject: [PATCH 16/20] feat(dh resp): rdhresp (wip) - add RespBase as base class - add rdhresp and support of dipole (XYG3 validated) - changed to_xxx, from_cls, and added to_resp --- pyscf/dh/energy/dhbase.py | 16 +- pyscf/dh/energy/hdft/rhdft.py | 4 + pyscf/dh/energy/mp2/rmp2.py | 4 + pyscf/dh/response/__init__.py | 4 + pyscf/dh/response/dh/rdh.py | 197 ++++++++++++++++++++++ pyscf/dh/response/hdft/rhdft.py | 58 +++++-- pyscf/dh/response/hdft/test/test_rhdft.py | 12 +- pyscf/dh/response/mp2/rmp2ri.py | 20 +-- pyscf/dh/response/respbase.py | 33 ++++ 9 files changed, 302 insertions(+), 46 deletions(-) create mode 100644 pyscf/dh/response/__init__.py create mode 100644 pyscf/dh/response/dh/rdh.py create mode 100644 pyscf/dh/response/respbase.py diff --git a/pyscf/dh/energy/dhbase.py b/pyscf/dh/energy/dhbase.py index 7b7d082..4ef5b69 100644 --- a/pyscf/dh/energy/dhbase.py +++ b/pyscf/dh/energy/dhbase.py @@ -149,10 +149,10 @@ class EngBase(lib.StreamObject, ABC): return frozen_core.mask @classmethod - def from_cls(cls, mf_dh: "EngBase", *args, **kwargs): + def from_cls(cls, mf_dh: "EngBase", *args, copy_all=False, **kwargs): """ Build a new (inherited) instance by a base EngBase instance. - Only some attributes are copied. Excluded attributes including + By default only some attributes are copied. Excluded attributes including - ``results`` - ``tensors`` - ``_tmpfile`` @@ -160,8 +160,12 @@ class EngBase(lib.StreamObject, ABC): """ mf_new = cls(*args, **kwargs) mf_new.__dict__.update(mf_dh.__dict__) - mf_new.results = dict() - mf_new.tensors = dict() - mf_new._tmpfile = lib.H5TmpFile() - mf_new.e_corr = NotImplemented + if not copy_all: + mf_new.results = dict() + mf_new.tensors = dict() + mf_new._tmpfile = lib.H5TmpFile() + mf_new.e_corr = NotImplemented return mf_new + + def to_resp(self): + raise NotImplementedError diff --git a/pyscf/dh/energy/hdft/rhdft.py b/pyscf/dh/energy/hdft/rhdft.py index 2cda438..a597819 100644 --- a/pyscf/dh/energy/hdft/rhdft.py +++ b/pyscf/dh/energy/hdft/rhdft.py @@ -456,6 +456,10 @@ class RHDFT(EngBase): self.scf.kernel(*args, **kwargs) return self.e_tot + def to_resp(self): + from pyscf.dh.response.hdft.rhdft import RHDFTResp + return RHDFTResp.from_cls(self, self.scf, copy_all=True) + get_energy_exactx = staticmethod(get_energy_restricted_exactx) get_energy_noxc = staticmethod(get_energy_restricted_noxc) get_energy_vv10 = staticmethod(get_energy_vv10) diff --git a/pyscf/dh/energy/mp2/rmp2.py b/pyscf/dh/energy/mp2/rmp2.py index 7404ebe..3fc90dd 100644 --- a/pyscf/dh/energy/mp2/rmp2.py +++ b/pyscf/dh/energy/mp2/rmp2.py @@ -401,6 +401,10 @@ class RMP2RI(MP2Base): self.results.update(results) return results + def to_resp(self): + from pyscf.dh.response.mp2.rmp2ri import RMP2RespRI + return RMP2RespRI.from_cls(self, self.scf, copy_all=True) + kernel_energy_mp2 = staticmethod(kernel_energy_rmp2_ri_incore) kernel = driver_eng_mp2 diff --git a/pyscf/dh/response/__init__.py b/pyscf/dh/response/__init__.py new file mode 100644 index 0000000..eea327d --- /dev/null +++ b/pyscf/dh/response/__init__.py @@ -0,0 +1,4 @@ +from .respbase import RespBase +from . import dh, hdft, mp2, respbase + + diff --git a/pyscf/dh/response/dh/rdh.py b/pyscf/dh/response/dh/rdh.py new file mode 100644 index 0000000..5e5ce89 --- /dev/null +++ b/pyscf/dh/response/dh/rdh.py @@ -0,0 +1,197 @@ +""" Doubly-Hybrid Response-Related Utilities. """ +import numpy as np + +from pyscf.dh.energy.dh import RDH +from pyscf.dh.response import RespBase +from pyscf.dh.response.hdft.rhdft import RHDFTResp +from pyscf.scf import cphf + + +class RDHResp(RDH, RespBase): + + def __init__(self, *args, **kwargs): + super().__init__(*args, **kwargs) + + # generate response instance for SCF (to obtain Ax0_Core) + HDFTResp = RHDFTResp if self.restricted else NotImplemented + self.scf_resp = HDFTResp(self.scf) + self._inherited_updated = False + + @property + def Ax0_Core(self): + """ Fock response of underlying SCF object in MO basis. """ + if self._Ax0_Core is NotImplemented: + self._Ax0_Core = self.scf_resp.Ax0_Core + return self._Ax0_Core + + @Ax0_Core.setter + def Ax0_Core(self, Ax0_Core): + self._Ax0_Core = Ax0_Core + + # in response instance, we first transfer all child instances into response + def to_scf(self, *args, **kwargs): + mf_resp = super().to_scf(*args, **kwargs).to_resp() + mf_resp.Ax0_Core = self.scf_resp.Ax0_Core + return mf_resp + + def to_mp2(self, *args, **kwargs): + assert len(args) == 0 + mf_resp = super().to_mp2(**kwargs).to_resp() + mf_resp.Ax0_Core = self.scf_resp.Ax0_Core + return mf_resp + + def to_iepa(self, *args, **kwargs): + raise NotImplementedError("Not implemented for response functions.") + + def to_ring_ccd(self, *args, **kwargs): + raise NotImplementedError("Not implemented for response functions.") + + def update_inherited(self): + """ Update inherited attribute, to transfer all child method instances to response. """ + + if self._inherited_updated: + return + + # if energy not evaluated, then evaluate energy first + if len(self.inherited) == 0: + self.kernel() + + # for energy evaluation, instance of low_rung may not be generated. + if len(self.inherited["low_rung"][1]) == 0: + HDFTResp = RHDFTResp if self.restricted else NotImplemented + self.inherited["low_rung"][1].append(HDFTResp(self.scf, xc=self.inherited["low_rung"][0])) + + # transform instances to response functions + # note that if Ax0_Core appears, then this object is already response, or have been inherited + for key in self.inherited: + for idx in range(len(self.inherited[key][1])): + instance = self.inherited[key][1][idx] + if not hasattr(instance, "Ax0_Core"): + instance = instance.to_resp() + self.inherited[key][1][idx] = instance + + self._inherited_updated = True + + def make_lag_vo(self): + if "lag_vo" in self.tensors: + return self.tensors["lag_vo"] + + self.update_inherited() + + nocc, nvir, nmo = self.nocc, self.nvir, self.nmo + lag_vo = np.zeros((nvir, nocc)) + for key in self.inherited: + for instance in self.inherited[key][1]: + lag_vo += instance.make_lag_vo() + + self.tensors["lag_vo"] = lag_vo + return lag_vo + + def make_rdm1_resp(self, ao=False): + if "rdm1_resp" in self.tensors: + return self.tensors["rdm1_resp"] + + self.update_inherited() + + rdm1_resp = np.diag(self.mo_occ) + for key in self.inherited: + for instance in self.inherited[key][1]: + if hasattr(instance, "make_rdm1_corr"): + rdm1_resp += instance.make_rdm1_corr() + + nocc, nvir, nmo = self.nocc, self.nvir, self.nmo + so, sv = slice(0, nocc), slice(nocc, nmo) + mo_energy = self.mo_energy + mo_occ = self.mo_occ + lag_vo = self.make_lag_vo() + max_cycle = self.max_cycle_cpks + tol = self.tol_cpks + + Ax0_Core = self.Ax0_Core + rdm1_resp_vo = cphf.solve( + Ax0_Core(sv, so, sv, so), mo_energy, mo_occ, lag_vo, + max_cycle=max_cycle, tol=tol)[0] + + rdm1_resp[sv, so] += rdm1_resp_vo + + self.tensors["rdm1_resp"] = rdm1_resp + + if ao: + rdm1_resp = self.mo_coeff @ rdm1_resp @ self.mo_coeff.T + return rdm1_resp + + def make_dipole(self): + # prepare input + mol = self.mol + rdm1_ao = self.make_rdm1_resp(ao=True) + int1e_r = mol.intor("int1e_r") + + dip_elec = - np.einsum("uv, tuv -> t", rdm1_ao, int1e_r) + dip_nuc = np.einsum("A, At -> t", mol.atom_charges(), mol.atom_coords()) + dip = dip_elec + dip_nuc + self.tensors["dipole"] = dip + return dip + + +if __name__ == '__main__': + def main_1(): + # test energy is the same for response and general DH + from pyscf import gto, scf + mol = gto.Mole(atom="O; H 1 0.94; H 1 0.94 2 104.5", basis="6-31G").build() + mf_resp = RDHResp(mol, xc="XYG3").run() + mf = RDH(mol, xc="XYG3").run() + print(np.allclose(mf_resp.e_tot, mf.e_tot)) + + def main_2(): + # test RSH functional + from pyscf import gto, scf + mol = gto.Mole(atom="O; H 1 0.94; H 1 0.94 2 104.5", basis="6-31G").build() + mf = RDH(mol, xc="RS-PBE-P86").run() + print(mf.results) + mf_resp = RDHResp(mol, xc="RS-PBE-P86").run() + print(mf_resp.results) + print(np.allclose(mf_resp.e_tot, mf.e_tot)) + + from pyscf import lib + lib.num_threads() + + def main_3(): + # test unbalanced OS contribution + from pyscf import gto, scf + mol = gto.Mole(atom="O; H 1 0.94; H 1 0.94 2 104.5", basis="6-31G").build() + mf = RDH(mol, xc="XYGJ-OS").run() + mf_resp = RDHResp(mol, xc="XYGJ-OS").run() + print(np.allclose(mf_resp.e_tot, mf.e_tot)) + + from pyscf import lib + + def main_4(): + # try response density and dipole + from pyscf import gto, scf, dft + + np.set_printoptions(8, suppress=True, linewidth=150) + + mol = gto.Mole(atom="O; H 1 0.94; H 1 0.94 2 104.5", basis="6-31G").build() + mf_resp = RDHResp(mol, xc="XYG3").run() + + def eng_with_dipole_field(t, h): + mf_scf = dft.RKS(mol, xc="B3LYPg").density_fit("cc-pVDZ-jkfit") + mf_scf.get_hcore = lambda *args, **kwargs: scf.hf.get_hcore(mol) - h * mol.intor("int1e_r")[t] + mf_scf.run(conv_tol=1e-12) + mf_mp = RDH(mf_scf, xc="XYG3").run() + return mf_mp.e_tot + + eng_array = np.zeros((2, 3)) + h = 1e-4 + for idx, h in [(0, h), [1, -h]]: + for t in (0, 1, 2): + eng_array[idx, t] = eng_with_dipole_field(t, h) + dip_elec_num = - (eng_array[0] - eng_array[1]) / (2 * h) + + dip_nuc = np.einsum("A, At -> t", mol.atom_charges(), mol.atom_coords()) + dip_anal = mf_resp.make_dipole() + dip_num = dip_elec_num + dip_nuc + print(dip_anal) + print(dip_num) + + main_4() diff --git a/pyscf/dh/response/hdft/rhdft.py b/pyscf/dh/response/hdft/rhdft.py index 799b664..d64806f 100644 --- a/pyscf/dh/response/hdft/rhdft.py +++ b/pyscf/dh/response/hdft/rhdft.py @@ -2,7 +2,8 @@ from pyscf.dh import RHDFT from pyscf.dh import util -from pyscf import gto, dft, lib, __config__ +from pyscf import gto, dft, lib, __config__, scf +from pyscf.dh.response import RespBase from pyscf.scf import _response_functions # this import is not unnecessary from pyscf.dh.energy.hdft.rhdft import get_rho import numpy as np @@ -236,7 +237,7 @@ def get_xc_integral(ni, mol, grids, xc, dm): return tensors -class RHDFTResp(RHDFT): +class RHDFTResp(RHDFT, RespBase): def __init__(self, *args, **kwargs): super().__init__(*args, **kwargs) @@ -284,7 +285,7 @@ class RHDFTResp(RHDFT): self.tensors.update(tensors) return cderi_uaa - def Ax0_Core(self, sp, sq, sr, ss): + def get_Ax0_Core(self, sp, sq, sr, ss): r""" Convenient function for evaluation of Fock response in MO basis :math:`\sum_{rs} A_{pq, rs} X_{rs}^\mathbb{A}`. @@ -305,7 +306,7 @@ class RHDFTResp(RHDFT): """ # if not RI, then use general Ax0_Core_resp if not hasattr(self.scf, "with_df") or not self.use_eri_cpks: - return self.Ax0_Core_resp(sp, sq, sr, ss) + return self.get_Ax0_Core_resp(sp, sq, sr, ss) # try if satisfies CPKS evaluation (vovo) lst_nmo = np.arange(self.nmo) @@ -315,15 +316,15 @@ class RHDFTResp(RHDFT): if ( np.all(lst_nmo[sp] == lst_vir) and np.all(lst_nmo[sq] == lst_occ) and np.all(lst_nmo[sr] == lst_vir) and np.all(lst_nmo[ss] == lst_occ)): - return self.Ax0_cpks() + return self.get_Ax0_cpks() else: # otherwise, use response by PySCF - return self.Ax0_Core_resp(sp, sq, sr, ss) + return self.get_Ax0_Core_resp(sp, sq, sr, ss) except ValueError: # dimension not match - return self.Ax0_Core_resp(sp, sq, sr, ss) + return self.get_Ax0_Core_resp(sp, sq, sr, ss) - def Ax0_Core_resp(self, sp, sq, sr, ss, vresp=None, mo_coeff=None): + def get_Ax0_Core_resp(self, sp, sq, sr, ss, vresp=None, mo_coeff=None): r""" Convenient function for evaluation of Fock response in MO basis :math:`\sum_{rs} A_{pq, rs} X_{rs}^\mathbb{A}` by PySCF's response function. @@ -394,7 +395,7 @@ class RHDFTResp(RHDFT): self.tensors["eri_cpks_vovo"] = eri_cpks_vovo return eri_cpks_vovo - def Ax0_cpks_HF(self): + def get_Ax0_cpks_HF(self): eri_cpks_vovo = self.tensors.get("eri_cpks_vovo", self.make_eri_cpks_vovo()) ax0_cpks_hf = Ax0_cpks_HF(eri_cpks_vovo, self.max_memory, self.verbose) return ax0_cpks_hf @@ -420,7 +421,7 @@ class RHDFTResp(RHDFT): self.tensors.update(tensors) return tensors - def Ax0_Core_KS(self, sp, sq, sr, ss): + def get_Ax0_Core_KS(self, sp, sq, sr, ss): ni = self.scf._numint mol = self.mol grids = self.grids_cpks @@ -443,17 +444,40 @@ class RHDFTResp(RHDFT): verbose=self.verbose) return ax0_core_ks - def Ax0_cpks(self): + def get_Ax0_cpks(self): nocc, nmo = self.nocc, self.nmo so, sv = slice(0, nocc), slice(nocc, nmo) - ax0_core_ks = self.Ax0_Core_KS(sv, so, sv, so) - ax0_cpks_hf = self.Ax0_cpks_HF() + ax0_core_ks = self.get_Ax0_Core_KS(sv, so, sv, so) + ax0_cpks_hf = self.get_Ax0_cpks_HF() def Ax0_cpks_inner(X): res = ax0_cpks_hf(X) + ax0_core_ks(X) return res return Ax0_cpks_inner + @property + def Ax0_Core(self): + """ Fock response of underlying SCF object in MO basis. """ + if self._Ax0_Core is NotImplemented: + self._Ax0_Core = self.get_Ax0_Core + return self._Ax0_Core + + @Ax0_Core.setter + def Ax0_Core(self, Ax0_Core): + self._Ax0_Core = Ax0_Core + + def make_lag_vo(self): + r""" Generate hybrid DFT contribution to Lagrangian vir-occ block :math:`L_{ai}`. """ + if "lag_vo" in self.tensors: + return self.tensors["lag_vo"] + + nocc, nmo = self.nocc, self.nmo + so, sv = slice(0, nocc), slice(nocc, nmo) + mo_coeff = self.mo_coeff + lag_vo = 4 * mo_coeff[:, sv].T @ self.hdft.get_fock(dm=self.scf.make_rdm1()) @ mo_coeff[:, so] + self.tensors["lag_vo"] = lag_vo + return lag_vo + if __name__ == '__main__': def main_1(): @@ -470,8 +494,8 @@ if __name__ == '__main__': nocc, nvir, nmo = mf_scf.nocc, mf_scf.nvir, mf_scf.nmo so, sv = slice(0, nocc), slice(nocc, nmo) X = np.random.randn(3, nvir, nocc) - ax_cpks = mf_scf.Ax0_cpks()(X) - ax_core = mf_scf.Ax0_Core_resp(sv, so, sv, so)(X) + ax_cpks = mf_scf.get_Ax0_cpks()(X) + ax_core = mf_scf.get_Ax0_Core_resp(sv, so, sv, so)(X) print(np.allclose(ax_cpks, ax_core)) # print(ax_core[0]) # print(ax_cpks[0]) @@ -490,8 +514,8 @@ if __name__ == '__main__': nocc, nvir, nmo = mf_scf.nocc, mf_scf.nvir, mf_scf.nmo so, sv = slice(0, nocc), slice(nocc, nmo) X = np.random.randn(3, nvir, nocc) - ax_cpks = mf_scf.Ax0_cpks()(X) - ax_core = mf_scf.Ax0_Core_resp(sv, so, sv, so)(X) + ax_cpks = mf_scf.get_Ax0_cpks()(X) + ax_core = mf_scf.get_Ax0_Core_resp(sv, so, sv, so)(X) print(np.allclose(ax_cpks, ax_core)) # print(ax_core[0]) # print(ax_cpks[0]) diff --git a/pyscf/dh/response/hdft/test/test_rhdft.py b/pyscf/dh/response/hdft/test/test_rhdft.py index ad48566..e1df99e 100644 --- a/pyscf/dh/response/hdft/test/test_rhdft.py +++ b/pyscf/dh/response/hdft/test/test_rhdft.py @@ -18,8 +18,8 @@ class TestRHDFTResp(unittest.TestCase): nocc, nvir, nmo = mf_scf.nocc, mf_scf.nvir, mf_scf.nmo so, sv = slice(0, nocc), slice(nocc, nmo) X = np.random.randn(3, nvir, nocc) - ax_cpks = mf_scf.Ax0_cpks()(X) - ax_core = mf_scf.Ax0_Core_resp(sv, so, sv, so)(X) + ax_cpks = mf_scf.get_Ax0_cpks()(X) + ax_core = mf_scf.get_Ax0_Core_resp(sv, so, sv, so)(X) self.assertTrue("eri_cpks_vovo" in mf_scf.tensors) self.assertTrue(np.allclose(ax_cpks, ax_core)) @@ -35,8 +35,8 @@ class TestRHDFTResp(unittest.TestCase): nocc, nvir, nmo = mf_scf.nocc, mf_scf.nvir, mf_scf.nmo so, sv = slice(0, nocc), slice(nocc, nmo) X = np.random.randn(3, nvir, nocc) - ax_cpks = mf_scf.Ax0_Core(sv, so, sv, so)(X) - ax_core = mf_scf.Ax0_Core_resp(sv, so, sv, so)(X) + ax_cpks = mf_scf.get_Ax0_Core(sv, so, sv, so)(X) + ax_core = mf_scf.get_Ax0_Core_resp(sv, so, sv, so)(X) self.assertTrue("eri_cpks_vovo" in mf_scf.tensors) self.assertTrue(np.allclose(ax_cpks, ax_core)) @@ -52,7 +52,7 @@ class TestRHDFTResp(unittest.TestCase): nocc, nvir, nmo = mf_scf.nocc, mf_scf.nvir, mf_scf.nmo so, sv = slice(0, nocc), slice(nocc, nmo) X = np.random.randn(3, nvir, nocc) - ax_cpks = mf_scf.Ax0_Core(sv, so, sv, so)(X) - ax_core = mf_scf.Ax0_Core_resp(sv, so, sv, so)(X) + ax_cpks = mf_scf.get_Ax0_Core(sv, so, sv, so)(X) + ax_core = mf_scf.get_Ax0_Core_resp(sv, so, sv, so)(X) self.assertTrue(np.allclose(ax_cpks, ax_core)) diff --git a/pyscf/dh/response/mp2/rmp2ri.py b/pyscf/dh/response/mp2/rmp2ri.py index 5e58876..b2d5b12 100644 --- a/pyscf/dh/response/mp2/rmp2ri.py +++ b/pyscf/dh/response/mp2/rmp2ri.py @@ -3,6 +3,7 @@ from pyscf.dh import RMP2RI from pyscf.dh import util from pyscf import scf, lib, __config__ +from pyscf.dh.response import RespBase from pyscf.scf import cphf import h5py import numpy as np @@ -303,7 +304,7 @@ def get_rdm1_corr_resp( return tensors -class RMP2RespRI(RMP2RI): +class RMP2RespRI(RMP2RI, RespBase): def __init__(self, *args, **kwargs): super().__init__(*args, **kwargs) @@ -411,6 +412,7 @@ class RMP2RespRI(RMP2RI): c_os, c_ss, incore_t_oovv, verbose=verbose, max_memory=max_memory, h5file=h5file) self.tensors.update(tensors) + results = {util.pad_omega(key, self.omega): val for (key, val) in results.items()} util.update_results(self.results, results) def make_t_oovv(self): @@ -558,22 +560,6 @@ class RMP2RespRI(RMP2RI): self.tensors["dipole"] = dip return dip - @property - def Ax0_Core(self): - """ Fock response of underlying SCF object in MO basis. """ - if self._Ax0_Core is NotImplemented: - restricted = isinstance(self.scf, scf.rhf.RHF) - from pyscf.dh.response.hdft import RHDFTResp - UHDFTResp = NotImplemented - HDFTResp = RHDFTResp if restricted else UHDFTResp - mf_scf = HDFTResp(self.scf) - self._Ax0_Core = mf_scf.Ax0_Core - return self._Ax0_Core - - @Ax0_Core.setter - def Ax0_Core(self, Ax0_Core): - self._Ax0_Core = Ax0_Core - kernel = driver_eng_mp2 get_mp2_integrals = staticmethod(get_mp2_integrals) diff --git a/pyscf/dh/response/respbase.py b/pyscf/dh/response/respbase.py new file mode 100644 index 0000000..153b24b --- /dev/null +++ b/pyscf/dh/response/respbase.py @@ -0,0 +1,33 @@ +""" Base class of response instances. """ + +from pyscf.dh.energy import EngBase +from pyscf import scf, __config__ + + +CONFIG_max_cycle_cpks = getattr(__config__, "max_cycle_cpks", 20) +CONFIG_tol_cpks = getattr(__config__, "tol_cpks", 1e-9) + + +class RespBase(EngBase): + def __init__(self, *args, **kwargs): + super().__init__(*args, **kwargs) + + self._Ax0_Core = NotImplemented + self.max_cycle_cpks = CONFIG_max_cycle_cpks + self.tol_cpks = CONFIG_tol_cpks + + @property + def Ax0_Core(self): + """ Fock response of underlying SCF object in MO basis. """ + if self._Ax0_Core is NotImplemented: + restricted = isinstance(self.scf, scf.rhf.RHF) + from pyscf.dh.response.hdft import RHDFTResp + UHDFTResp = NotImplemented + HDFTResp = RHDFTResp if restricted else UHDFTResp + mf_scf = HDFTResp(self.scf) + self._Ax0_Core = mf_scf.get_Ax0_Core + return self._Ax0_Core + + @Ax0_Core.setter + def Ax0_Core(self, Ax0_Core): + self._Ax0_Core = Ax0_Core -- Gitee From 27f6543802d99bcd1998e9813bd86e5bd5073d69 Mon Sep 17 00:00:00 2001 From: ajz34 Date: Wed, 10 May 2023 13:40:28 +0800 Subject: [PATCH 17/20] fix(dh resp): allow Ax0_cpks accept HF instance --- pyscf/dh/response/hdft/rhdft.py | 13 +++++++++++-- 1 file changed, 11 insertions(+), 2 deletions(-) diff --git a/pyscf/dh/response/hdft/rhdft.py b/pyscf/dh/response/hdft/rhdft.py index d64806f..03061f6 100644 --- a/pyscf/dh/response/hdft/rhdft.py +++ b/pyscf/dh/response/hdft/rhdft.py @@ -359,8 +359,13 @@ class RHDFTResp(RHDFT, RespBase): # This function will utilize mf.scf._numint to specify hybrid and omega coefficients - ni = self.scf._numint # type: dft.numint.NumInt - omega, alpha, hyb = ni.rsh_and_hybrid_coeff(self.scf.xc) + if not hasattr(self.scf, "xc"): + # not a DFT instance, configure as HF instance + omega, alpha, hyb = 0, 0, 1 + else: + ni = self.scf._numint # type: dft.numint.NumInt + omega, alpha, hyb = ni.rsh_and_hybrid_coeff(self.scf.xc) + nvir, nocc = self.nvir, self.nocc cderi_uaa = self.tensors.get("cderi_uaa", self.make_cderi_uaa()) @@ -422,6 +427,10 @@ class RHDFTResp(RHDFT, RespBase): return tensors def get_Ax0_Core_KS(self, sp, sq, sr, ss): + if not hasattr(self.scf, "xc"): + # not a DFT instance, KS contribution is zero + return lambda *args, **kwargs: 0 + ni = self.scf._numint mol = self.mol grids = self.grids_cpks -- Gitee From f15f55f09cf39405a3821795c21f93101d4358be Mon Sep 17 00:00:00 2001 From: ajz34 Date: Wed, 10 May 2023 14:28:28 +0800 Subject: [PATCH 18/20] feat(dh resp): rdhresp (wip) - abstract function get_rdm1_resp_vo --- pyscf/dh/response/dh/rdh.py | 26 ++++++++-------- pyscf/dh/response/mp2/rmp2ri.py | 11 +++---- pyscf/dh/response/respbase.py | 53 ++++++++++++++++++++++++++++++++- 3 files changed, 69 insertions(+), 21 deletions(-) diff --git a/pyscf/dh/response/dh/rdh.py b/pyscf/dh/response/dh/rdh.py index 5e5ce89..58c8cef 100644 --- a/pyscf/dh/response/dh/rdh.py +++ b/pyscf/dh/response/dh/rdh.py @@ -4,7 +4,7 @@ import numpy as np from pyscf.dh.energy.dh import RDH from pyscf.dh.response import RespBase from pyscf.dh.response.hdft.rhdft import RHDFTResp -from pyscf.scf import cphf +from pyscf.dh.response.respbase import get_rdm1_resp_vo_restricted class RDHResp(RDH, RespBase): @@ -106,11 +106,12 @@ class RDHResp(RDH, RespBase): lag_vo = self.make_lag_vo() max_cycle = self.max_cycle_cpks tol = self.tol_cpks - Ax0_Core = self.Ax0_Core - rdm1_resp_vo = cphf.solve( - Ax0_Core(sv, so, sv, so), mo_energy, mo_occ, lag_vo, - max_cycle=max_cycle, tol=tol)[0] + verbose = self.verbose + + rdm1_resp_vo = get_rdm1_resp_vo_restricted( + lag_vo, mo_energy, mo_occ, Ax0_Core, + max_cycle=max_cycle, tol=tol, verbose=verbose) rdm1_resp[sv, so] += rdm1_resp_vo @@ -140,7 +141,7 @@ if __name__ == '__main__': mol = gto.Mole(atom="O; H 1 0.94; H 1 0.94 2 104.5", basis="6-31G").build() mf_resp = RDHResp(mol, xc="XYG3").run() mf = RDH(mol, xc="XYG3").run() - print(np.allclose(mf_resp.e_tot, mf.e_tot)) + assert np.allclose(mf_resp.e_tot, mf.e_tot) def main_2(): # test RSH functional @@ -150,10 +151,7 @@ if __name__ == '__main__': print(mf.results) mf_resp = RDHResp(mol, xc="RS-PBE-P86").run() print(mf_resp.results) - print(np.allclose(mf_resp.e_tot, mf.e_tot)) - - from pyscf import lib - lib.num_threads() + assert np.allclose(mf_resp.e_tot, mf.e_tot) def main_3(): # test unbalanced OS contribution @@ -161,9 +159,7 @@ if __name__ == '__main__': mol = gto.Mole(atom="O; H 1 0.94; H 1 0.94 2 104.5", basis="6-31G").build() mf = RDH(mol, xc="XYGJ-OS").run() mf_resp = RDHResp(mol, xc="XYGJ-OS").run() - print(np.allclose(mf_resp.e_tot, mf.e_tot)) - - from pyscf import lib + assert np.allclose(mf_resp.e_tot, mf.e_tot) def main_4(): # try response density and dipole @@ -191,7 +187,11 @@ if __name__ == '__main__': dip_nuc = np.einsum("A, At -> t", mol.atom_charges(), mol.atom_coords()) dip_anal = mf_resp.make_dipole() dip_num = dip_elec_num + dip_nuc + np.allclose(dip_anal, dip_num, rtol=1e-5, atol=1e-6) print(dip_anal) print(dip_num) + main_1() + main_2() + main_3() main_4() diff --git a/pyscf/dh/response/mp2/rmp2ri.py b/pyscf/dh/response/mp2/rmp2ri.py index b2d5b12..f826374 100644 --- a/pyscf/dh/response/mp2/rmp2ri.py +++ b/pyscf/dh/response/mp2/rmp2ri.py @@ -4,7 +4,7 @@ from pyscf.dh import RMP2RI from pyscf.dh import util from pyscf import scf, lib, __config__ from pyscf.dh.response import RespBase -from pyscf.scf import cphf +from pyscf.dh.response.respbase import get_rdm1_resp_vo_restricted import h5py import numpy as np @@ -284,19 +284,16 @@ def get_rdm1_corr_resp( nocc = (mo_occ > 0).sum() nvir = (mo_occ == 0).sum() nmo = nocc + nvir - assert lag_vo.shape == (nvir, nocc) assert rdm1_corr.shape == (nmo, nmo) - assert len(mo_energy) == nmo - assert len(mo_occ) == nmo # prepare essential matrices and slices so, sv = slice(0, nocc), slice(nocc, nmo) # cp-ks evaluation rdm1_corr_resp = rdm1_corr.copy() - rdm1_corr_resp_vo = cphf.solve( - Ax0_Core(sv, so, sv, so), mo_energy, mo_occ, lag_vo, - max_cycle=max_cycle, tol=tol)[0] + rdm1_corr_resp_vo = get_rdm1_resp_vo_restricted( + lag_vo, mo_energy, mo_occ, Ax0_Core, + max_cycle=max_cycle, tol=tol, verbose=verbose) rdm1_corr_resp[sv, so] = rdm1_corr_resp_vo log.timer("get_rdm1_corr_resp", *time0) diff --git a/pyscf/dh/response/respbase.py b/pyscf/dh/response/respbase.py index 153b24b..569286f 100644 --- a/pyscf/dh/response/respbase.py +++ b/pyscf/dh/response/respbase.py @@ -1,7 +1,8 @@ """ Base class of response instances. """ from pyscf.dh.energy import EngBase -from pyscf import scf, __config__ +from pyscf import scf, __config__, lib +from pyscf.scf import cphf CONFIG_max_cycle_cpks = getattr(__config__, "max_cycle_cpks", 20) @@ -31,3 +32,53 @@ class RespBase(EngBase): @Ax0_Core.setter def Ax0_Core(self, Ax0_Core): self._Ax0_Core = Ax0_Core + + +def get_rdm1_resp_vo_restricted( + lag_vo, mo_energy, mo_occ, Ax0_Core, + max_cycle=20, tol=1e-9, verbose=lib.logger.NOTE): + r""" Solution of response 1-RDM vir-occ part by CP-KS equation (restricted). + + .. math:: + A'_{ai, bj} D_{bj} &= L_{ai} + + Parameters + ---------- + lag_vo : np.ndarray + mo_energy : np.ndarray + mo_occ : np.ndarray + Ax0_Core : callable + max_cycle : int + tol : float + verbose : int + + Returns + ------- + dict[str, np.ndarray] + + See Also + -------- + get_rdm1_corr + """ + + log = lib.logger.new_logger(verbose=verbose) + time0 = lib.logger.process_clock(), lib.logger.perf_counter() + + # dimension definition and check sanity + nocc = (mo_occ > 0).sum() + nvir = (mo_occ == 0).sum() + nmo = nocc + nvir + assert lag_vo.shape == (nvir, nocc) + assert len(mo_energy) == nmo + assert len(mo_occ) == nmo + + # prepare essential matrices and slices + so, sv = slice(0, nocc), slice(nocc, nmo) + + # cp-ks evaluation + rdm1_vo = cphf.solve( + Ax0_Core(sv, so, sv, so), mo_energy, mo_occ, lag_vo, + max_cycle=max_cycle, tol=tol)[0] + + log.timer("get_rdm1_resp_vo_restricted", *time0) + return rdm1_vo -- Gitee From 4c8c114aa976b41d627654632e4e557cbfb612c2 Mon Sep 17 00:00:00 2001 From: ajz34 Date: Wed, 10 May 2023 15:19:09 +0800 Subject: [PATCH 19/20] fix(dh eng): allow larger dimension in array_add_with_diff_rows --- pyscf/dh/energy/hdft/rhdft.py | 22 ++++++++++++++++------ 1 file changed, 16 insertions(+), 6 deletions(-) diff --git a/pyscf/dh/energy/hdft/rhdft.py b/pyscf/dh/energy/hdft/rhdft.py index 2cda438..0d4050f 100644 --- a/pyscf/dh/energy/hdft/rhdft.py +++ b/pyscf/dh/energy/hdft/rhdft.py @@ -282,13 +282,23 @@ def numint_customized(xc, _mol=None): def array_add_with_diff_rows(mat1, mat2): assert len(mat1.shape) == len(mat2.shape) - assert mat1.shape[1:] == mat2.shape[1:] - if len(mat1) >= len(mat2): - mat1[:len(mat2)] += mat2 - return mat1 + assert mat1.shape[-1] == mat2.shape[-1] + if len(mat1) < len(mat2): + mat1, mat2 = mat2, mat1 + + if mat1.ndim == 1: + mat1 += mat2 + elif mat1.ndim == 2: + mat1[:mat2.shape[0]] += mat2 + elif mat1.ndim == 3: + mat1[:mat2.shape[0], :mat2.shape[1]] += mat2 + elif mat1.ndim == 4: + mat1[:mat2.shape[0], :mat2.shape[1], :mat2.shape[2]] += mat2 + elif mat1.ndim == 5: + mat1[:mat2.shape[0], :mat2.shape[1], :mat2.shape[2], :mat2.shape[3]] += mat2 else: - mat2[:len(mat1)] += mat1 - return mat2 + raise NotImplementedError + return mat1 def eval_xc_eff(*args, **kwargs): exc, vxc, fxc, kxc = gen_lists[0](*args, **kwargs) -- Gitee From 6c90e43ae3ac259fef4807384c32fe08dc6001e0 Mon Sep 17 00:00:00 2001 From: ajz34 Date: Wed, 10 May 2023 15:21:16 +0800 Subject: [PATCH 20/20] feat(dh resp): rdhresp (wip) - add range separate function dipole (customize numint, various omega of MP2) --- pyscf/dh/response/dh/rdh.py | 40 +++++++++++++++++++++++++++++---- pyscf/dh/response/hdft/rhdft.py | 17 +++++++++++++- pyscf/dh/response/mp2/rmp2ri.py | 11 ++++++--- 3 files changed, 60 insertions(+), 8 deletions(-) diff --git a/pyscf/dh/response/dh/rdh.py b/pyscf/dh/response/dh/rdh.py index 58c8cef..d8bc107 100644 --- a/pyscf/dh/response/dh/rdh.py +++ b/pyscf/dh/response/dh/rdh.py @@ -191,7 +191,39 @@ if __name__ == '__main__': print(dip_anal) print(dip_num) - main_1() - main_2() - main_3() - main_4() + def main_5(): + # try response density and dipole + from pyscf import gto, scf, dft + + np.set_printoptions(8, suppress=True, linewidth=150) + + mol = gto.Mole(atom="O; H 1 0.94; H 1 0.94 2 104.5", basis="6-31G").build() + xc_code = "0.5*HF + 0.5*LR_HF(0.7) + 0.5*SSR(GGA_X_PBE, 0.7), 0.75*SSR(GGA_C_P86, 0.7) + MP2(0.5, 0.5) + RS_MP2(-0.7, -0.25, -0.25) + RS_MP2(0.7, 0.5, 0.5)" + mf_resp = RDHResp(mol, xc=xc_code).run() + + def eng_with_dipole_field(t, h): + mf_mp = RDH(mol, xc=xc_code) + mf_mp.scf.get_hcore = lambda *args, **kwargs: scf.hf.get_hcore(mol) - h * mol.intor("int1e_r")[t] + mf_mp.scf.conv_tol = 1e-12 + mf_mp.run() + return mf_mp.e_tot + + eng_array = np.zeros((2, 3)) + h = 1e-4 + for idx, h in [(0, h), [1, -h]]: + for t in (0, 1, 2): + eng_array[idx, t] = eng_with_dipole_field(t, h) + dip_elec_num = - (eng_array[0] - eng_array[1]) / (2 * h) + + dip_nuc = np.einsum("A, At -> t", mol.atom_charges(), mol.atom_coords()) + dip_anal = mf_resp.make_dipole() + dip_num = dip_elec_num + dip_nuc + np.allclose(dip_anal, dip_num, rtol=1e-5, atol=1e-6) + print(dip_anal) + print(dip_num) + + # main_1() + # main_2() + # main_3() + # main_4() + main_5() diff --git a/pyscf/dh/response/hdft/rhdft.py b/pyscf/dh/response/hdft/rhdft.py index 03061f6..9ce8aca 100644 --- a/pyscf/dh/response/hdft/rhdft.py +++ b/pyscf/dh/response/hdft/rhdft.py @@ -253,7 +253,22 @@ class RHDFTResp(RHDFT, RespBase): @cached_property def vresp(self): """ Fock response function (derivative w.r.t. molecular coefficient in AO basis). """ - return self.scf.gen_response() + try: + return self.scf.gen_response() + except ValueError: + # case that have customized xc + # should pass check of ni.libxc.test_deriv_order and ni.libxc.is_hybrid_xc + ni = self.scf._numint + omega, alpha, hyb = ni.rsh_and_hybrid_coeff(self.scf.xc, self.mol.spin) + if abs(hyb) > 1e-10 or abs(omega) > 1e-10: + fake_xc = "B3LYPg" + else: + fake_xc = "PBE" + actual_xc = self.scf.xc + self.scf.xc = fake_xc + resp = self.scf.gen_response() + self.scf.xc = actual_xc + return resp def make_cderi_uaa(self, omega=0): """ Generate cholesky decomposed ERI (all block, full orbitals, s1 symm, in memory/disk). """ diff --git a/pyscf/dh/response/mp2/rmp2ri.py b/pyscf/dh/response/mp2/rmp2ri.py index f826374..bf33436 100644 --- a/pyscf/dh/response/mp2/rmp2ri.py +++ b/pyscf/dh/response/mp2/rmp2ri.py @@ -385,9 +385,14 @@ class RMP2RespRI(RMP2RI, RespBase): if "cderi_uoo" in self.tensors: return self.tensors["cderi_uoo"] - mo_coeff_occ = self.mo_coeff[:, self.mask_occ] - with_df = self.with_df - cderi_uoo = util.get_cderi_mo(with_df, mo_coeff_occ, max_memory=self.max_memory) + # dimension and mask + mask_occ_act = self.mask_occ_act + + # In response evaluation, ERI of all orbitals are generally required. + # Thus, cderi_uaa is a must-evaluated tensor. + # note that cderi_uaa may be stored by h5py, so indexing by list should be done in following way + cderi_uaa = self.tensors.get("cderi_uaa", self.make_cderi_uaa()) + cderi_uoo = cderi_uaa[:, mask_occ_act, :][:, :, mask_occ_act] tensors = {"cderi_uoo": cderi_uoo} self.tensors.update(tensors) -- Gitee