From 31e1bec6184ebc8021658e8c73866f70c9170f42 Mon Sep 17 00:00:00 2001 From: ajz34 Date: Fri, 21 Apr 2023 11:28:42 +0800 Subject: [PATCH 01/12] docs(dh eng): add comment on custom_mf --- pyscf/dh/energy/rhdft.py | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/pyscf/dh/energy/rhdft.py b/pyscf/dh/energy/rhdft.py index af61d01..4a971e4 100644 --- a/pyscf/dh/energy/rhdft.py +++ b/pyscf/dh/energy/rhdft.py @@ -327,6 +327,11 @@ def custom_mf(mf, xc, auxbasis_or_with_df=None): Returns ------- dft.rks.RKS or dft.uks.UKS + + Notes + ----- + 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. """ verbose = mf.verbose log = lib.logger.new_logger(verbose=verbose) -- Gitee From 66ee971a429884ebebda80fe338078fea64c42b4 Mon Sep 17 00:00:00 2001 From: ajz34 Date: Fri, 21 Apr 2023 15:51:54 +0800 Subject: [PATCH 02/12] feat(dh eng): add Ax0_Core_resp --- pyscf/dh/energy/rhdft.py | 93 ++++++++++++++++++++++++++++++++++++++-- 1 file changed, 90 insertions(+), 3 deletions(-) diff --git a/pyscf/dh/energy/rhdft.py b/pyscf/dh/energy/rhdft.py index 4a971e4..d79cf01 100644 --- a/pyscf/dh/energy/rhdft.py +++ b/pyscf/dh/energy/rhdft.py @@ -2,8 +2,10 @@ from pyscf.dh import util from pyscf.dh.util import XCList, XCType, XCInfo from pyscf.dh.energy import EngBase from pyscf import dft, lib, scf, df, __config__ +from pyscf.scf import _response_functions # this import is not unnecessary import numpy as np from types import MethodType +from functools import cache, cached_property CONFIG_ssr_x_fr = getattr(__config__, "ssr_x_fr", "LDA_X") CONFIG_ssr_x_sr = getattr(__config__, "ssr_x_sr", "LDA_X_ERF") @@ -108,7 +110,7 @@ def get_energy_purexc(xc_lists, rho, weights, restricted, numint=None): """ Evaluate energy contributions of pure (DFT) exchange-correlation effects. Note that this kernel does not count HF, LR_HF and advanced correlation into account. - To evaluate exact exchange (HF or LR_HF), use ``kernel_energy_restricted_exactx``. + To evaluate exact exchange (HF or LR_HF), use ``get_energy_restricted_exactx``. Parameters ---------- @@ -131,7 +133,7 @@ def get_energy_purexc(xc_lists, rho, weights, restricted, numint=None): See Also -------- - kernel_energy_restricted_exactx + get_energy_restricted_exactx """ if not isinstance(xc_lists, list): xc_lists = [xc_lists] @@ -309,6 +311,40 @@ def numint_customized(xc, _mol=None): return ni_custom +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}`. + + 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 + + def custom_mf(mf, xc, auxbasis_or_with_df=None): """ Customize options of PySCF's mf object. @@ -444,6 +480,39 @@ class RHDFT(EngBase): def kernel(self, *args, **kwargs): return self.scf.kernel(*args, **kwargs) + @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_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}`. + + 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) + get_energy_exactx = staticmethod(get_energy_restricted_exactx) get_energy_noxc = staticmethod(get_energy_restricted_noxc) get_energy_vv10 = staticmethod(get_energy_vv10) @@ -460,4 +529,22 @@ if __name__ == '__main__': res = mf.make_energy_purexc([", LYP", "B88, ", "HF", "LR_HF(0.5)", "SSR(GGA_X_B88, 0.5), "]) print(res) - main_1() + def main_2(): + from pyscf import gto, scf + mol = gto.Mole(atom="O; H 1 0.94; H 1 0.94 2 104.5").build() + mf_s = scf.RHF(mol) + mf = RHDFT(mf_s, xc="B3LYPg").run() + np.random.seed(0) + nocc, nvir = mf.nocc, mf.nvir + dmX = np.random.randn(3, nocc, nvir) + mo_occ = mf.mo_occ + idx_occ = mo_occ > 0 + idx_vir = ~idx_occ + log = lib.logger.new_logger(verbose=5) + # first epoch should be a lot slower, since mf.scf.gen_response will cache xc kernel first + for epoch in range(3): + t0, t1 = lib.logger.process_clock(), lib.logger.perf_counter() + ax = mf.Ax0_Core_resp(idx_occ, idx_vir, idx_occ, idx_vir)(dmX) + lib.logger.timer(log, f"Epoch {epoch} of Ax0_Core_resp", t0, t1) + + main_2() -- Gitee From a555ccaf1cad6a44ee465c06e76a027b993c8b0e Mon Sep 17 00:00:00 2001 From: ajz34 Date: Thu, 27 Apr 2023 14:33:44 +0800 Subject: [PATCH 03/12] feat(dh eng): add dipole (wip) --- pyscf/dh/__init__.py | 2 +- pyscf/dh/energy/__init__.py | 2 +- pyscf/dh/energy/dh.py | 31 +- pyscf/dh/energy/mp2/__init__.py | 3 +- pyscf/dh/energy/mp2/rmp2.py | 186 +--------- pyscf/dh/energy/mp2/rmp2ri.py | 633 ++++++++++++++++++++++++++++++++ pyscf/dh/energy/mp2/ump2.py | 24 +- pyscf/dh/energy/rhdft.py | 38 +- pyscf/dh/util/df_addon.py | 20 +- pyscf/dh/util/frozen_core.py | 7 +- pyscf/dh/util/helper.py | 10 +- 11 files changed, 747 insertions(+), 209 deletions(-) create mode 100644 pyscf/dh/energy/mp2/rmp2ri.py diff --git a/pyscf/dh/__init__.py b/pyscf/dh/__init__.py index 22b0475..0fd3963 100644 --- a/pyscf/dh/__init__.py +++ b/pyscf/dh/__init__.py @@ -3,7 +3,7 @@ from . import energy from .energy import ( RHDFT, UHDFT, - RMP2RIPySCF, RMP2ConvPySCF, UMP2ConvPySCF, RMP2Conv, RMP2RI, UMP2Conv, UMP2RI, + RMP2RIPySCF, RMP2ConvPySCF, UMP2ConvPySCF, RMP2Conv, UMP2Conv, UMP2RI, RIEPARI, RIEPAConv, UIEPAConv, UIEPARI, RRingCCDConv ) diff --git a/pyscf/dh/energy/__init__.py b/pyscf/dh/energy/__init__.py index 843baee..04b1d25 100644 --- a/pyscf/dh/energy/__init__.py +++ b/pyscf/dh/energy/__init__.py @@ -6,7 +6,7 @@ from .uhdft import UHDFT from . import mp2, iepa from .mp2 import ( RMP2RIPySCF, RMP2ConvPySCF, UMP2ConvPySCF, - RMP2Conv, RMP2RI, UMP2Conv, UMP2RI) + RMP2Conv, UMP2Conv, RMP2RI, UMP2RI) from .iepa import (RIEPAConv, RIEPARI, UIEPAConv, UIEPARI) from .ring_ccd import RRingCCDConv diff --git a/pyscf/dh/energy/dh.py b/pyscf/dh/energy/dh.py index ef45357..e14e2f1 100644 --- a/pyscf/dh/energy/dh.py +++ b/pyscf/dh/energy/dh.py @@ -86,16 +86,6 @@ def _process_energy_mp2(mf_dh, xc_list, force_evaluate=False): log.info(f"[INFO] XCList extracted by process_energy_mp2: {xc_mp2.token}") log.info(f"[INFO] XCList remains by process_energy_mp2: {xc_extracted.token}") log.info("[INFO] MP2 detected") - # generate omega list - # parameter of RSMP2: omega, c_os, c_ss - omega_list = [] - for info in xc_mp2: - if XCType.MP2 in info.type: - omega_list.append(0) - else: - assert XCType.RSMP2 in info.type - omega_list.append(info.parameters[0]) - assert len(set(omega_list)) == len(omega_list) def comput_mp2(): eng_tot = 0 @@ -106,18 +96,21 @@ def _process_energy_mp2(mf_dh, xc_list, force_evaluate=False): omega = 0 else: omega, c_os, c_ss = info.parameters - eng = info.fac * ( - + c_os * results[pad_omega("eng_corr_MP2_OS", omega)] - + c_ss * results[pad_omega("eng_corr_MP2_SS", omega)]) + eng = info.fac * results[pad_omega("eng_corr_MP2", omega)] log.note(f"[RESULT] Energy of correlation {info.token}: {eng:20.12f}") eng_tot += eng return eng_tot def force_comput_mp2(): - for omega in omega_list: + for info in xc_mp2: + if XCType.MP2 in info.type: + c_os, c_ss = info.parameters + omega = 0 + else: + omega, c_os, c_ss = info.parameters if pad_omega("eng_corr_MP2_SS", omega) in mf_dh.results and not force_evaluate: continue # results have been evaluated - mf_mp2 = mf_dh.to_mp2(omega=omega).run() + mf_mp2 = mf_dh.to_mp2(omega=omega, c_os=c_os, c_ss=c_ss).run() update_results(mf_dh.results, mf_mp2.results) mf_dh.inherited.append((mf_mp2, xc_mp2)) eng_tot = comput_mp2() @@ -596,12 +589,12 @@ class DH(EngBase): def to_scf(self, **kwargs): # import if self.restricted: - from pyscf.dh import RHDFT as DHSCF + from pyscf.dh import RHDFT as HDFT else: - from pyscf.dh import UHDFT as DHSCF + from pyscf.dh import UHDFT as HDFT # generate instance - mf = DHSCF.from_rdh(self, self.scf, self.xc.xc_scf, **kwargs) + mf = HDFT.from_rdh(self, self.scf, self.xc.xc_scf, **kwargs) return mf @@ -609,7 +602,7 @@ class DH(EngBase): # import if self.restricted: from pyscf.dh import RMP2Conv as MP2Conv - from pyscf.dh import RMP2RI as MP2RI + from pyscf.dh.energy.mp2 import RMP2RI as MP2RI else: from pyscf.dh import UMP2Conv as MP2Conv from pyscf.dh import UMP2RI as MP2RI diff --git a/pyscf/dh/energy/mp2/__init__.py b/pyscf/dh/energy/mp2/__init__.py index 5ed9b16..cbe73a7 100644 --- a/pyscf/dh/energy/mp2/__init__.py +++ b/pyscf/dh/energy/mp2/__init__.py @@ -1,3 +1,4 @@ from . import rmp2, ump2 -from .rmp2 import RMP2ConvPySCF, RMP2Conv, RMP2RI, RMP2RIPySCF +from .rmp2 import RMP2ConvPySCF, RMP2Conv, RMP2RIPySCF +from .rmp2ri import RMP2RI from .ump2 import UMP2ConvPySCF, UMP2Conv, UMP2RI diff --git a/pyscf/dh/energy/mp2/rmp2.py b/pyscf/dh/energy/mp2/rmp2.py index 7404ebe..40cd88b 100644 --- a/pyscf/dh/energy/mp2/rmp2.py +++ b/pyscf/dh/energy/mp2/rmp2.py @@ -56,9 +56,11 @@ CONFIG_etb_first = getattr(__config__, "etb_first", False) class MP2Base(EngBase): """ Restricted MP2 base class. """ - def __init__(self, mf, frozen=None, omega=0, with_df=None, **kwargs): + def __init__(self, mf, frozen=None, omega=0, with_df=None, c_os=1.0, c_ss=1.0, **kwargs): super().__init__(mf) self.omega = omega + self.c_os = c_os + self.c_ss = c_ss self.incore_t_oovv_mp2 = CONFIG_incore_t_oovv_mp2 self.frozen = frozen if frozen is not None else 0 self.frac_num = None @@ -76,8 +78,10 @@ class MP2Base(EngBase): class RMP2ConvPySCF(mp.mp2.RMP2, EngBase): """ Restricted MP2 class of doubly hybrid with conventional integral evaluated by PySCF. """ - def __init__(self, mf, *args, **kwargs): + def __init__(self, mf, c_os=1.0, c_ss=1.0, *args, **kwargs): EngBase.__init__(self, mf) + self.c_os = c_os + self.c_ss = c_ss super().__init__(mf, *args, **kwargs) def get_frozen_mask(self): # type: () -> np.ndarray @@ -87,7 +91,7 @@ class RMP2ConvPySCF(mp.mp2.RMP2, EngBase): kernel_output = super().kernel(*args, **kwargs) self.results["eng_corr_MP2_OS"] = self.e_corr_os self.results["eng_corr_MP2_SS"] = self.e_corr_ss - self.results["eng_corr_MP2"] = self.e_corr + self.results["eng_corr_MP2"] = self.c_os * self.e_corr_os + self.c_ss * self.e_corr_ss return kernel_output kernel = driver_eng_mp2 @@ -106,8 +110,10 @@ class RMP2RIPySCF(mp.dfmp2.DFMP2, EngBase): def update_amps(self, t2, eris): raise NotImplementedError - def __init__(self, mf, *args, **kwargs): + def __init__(self, mf, c_os=1.0, c_ss=1.0, *args, **kwargs): EngBase.__init__(self, mf) + self.c_os = c_os + self.c_ss = c_ss super().__init__(mf, *args, **kwargs) def get_frozen_mask(self): # type: () -> np.ndarray @@ -117,7 +123,7 @@ class RMP2RIPySCF(mp.dfmp2.DFMP2, EngBase): kernel_output = super().kernel(*args, **kwargs) self.results["eng_corr_MP2_OS"] = self.e_corr_os self.results["eng_corr_MP2_SS"] = self.e_corr_ss - self.results["eng_corr_MP2"] = self.e_corr + self.results["eng_corr_MP2"] = self.c_os * self.e_corr_os + self.c_ss * self.e_corr_ss return kernel_output kernel = driver_eng_mp2 @@ -128,7 +134,7 @@ class RMP2RIPySCF(mp.dfmp2.DFMP2, EngBase): # region RMP2Conv def kernel_energy_rmp2_conv_full_incore( - mo_energy, mo_coeff, eri_or_mol, mo_occ, + mo_energy, mo_coeff, eri_or_mol, mo_occ, c_os, c_ss, t_oovv=None, frac_num=None, verbose=lib.logger.NOTE, **_kwargs): """ Kernel of restricted MP2 energy by conventional method. @@ -142,6 +148,10 @@ def kernel_energy_rmp2_conv_full_incore( ERI that is recognized by ``pyscf.ao2mo.general``. mo_occ : np.ndarray Molecular orbital occupation numbers. + c_os : float + Coefficient of oppo-spin contribution. + c_ss : float + Coefficient of same-spin contribution. t_oovv : np.ndarray Store space for ``t_oovv`` @@ -203,7 +213,7 @@ def kernel_energy_rmp2_conv_full_incore( # report eng_os = eng_bi1 eng_ss = eng_bi1 - eng_bi2 - eng_mp2 = eng_os + eng_ss + eng_mp2 = c_os * eng_os + c_ss * eng_ss # results results = dict() results["eng_corr_MP2_bi1"] = eng_bi1 @@ -244,6 +254,7 @@ class RMP2Conv(MP2Base): eri_or_mol = eri_or_mol if eri_or_mol is not None else self.mol results = self.kernel_energy_mp2( mo_energy_act, mo_coeff_act, eri_or_mol, mo_occ_act, + self.c_os, self.c_ss, t_oovv=t_oovv, frac_num=frac_num_act, verbose=self.verbose, **kwargs) self.e_corr = results["eng_corr_MP2"] # pad omega @@ -257,159 +268,10 @@ class RMP2Conv(MP2Base): # endregion -# region RMP2RI - -def kernel_energy_rmp2_ri_incore( - mo_energy, cderi_uov, - t_oovv=None, frac_num=None, verbose=lib.logger.NOTE, max_memory=2000, cderi_uov_2=None): - """ Kernel of MP2 energy by RI integral. - - Parameters - ---------- - mo_energy : np.ndarray - Molecular orbital energy levels. - cderi_uov : np.ndarray - Cholesky decomposed 3c2e ERI in MO basis (occ-vir part). - - t_oovv : np.ndarray - Store space for ``t_oovv`` - frac_num : np.ndarray - Fractional occupation number list. - verbose : int - Verbose level for PySCF. - max_memory : float - Allocatable memory in MB. - cderi_uov_2 : np.ndarray - Another part of 3c2e ERI in MO basis (occ-vir part). This is mostly used in magnetic computations. - - Notes - ----- - - For RI approximation, ERI integral is set to be - - .. math:: - g_{ij}^{ab} = (ia|jb) = Y_{ia, P} Y_{jb, P} - - For energy of fractional occupation system, computation method is chosen - by eq (7) from Su2016 (10.1021/acs.jctc.6b00197). - """ - log = lib.logger.new_logger(verbose=verbose) - log.info("[INFO] Start unrestricted RI-MP2") - - naux, nocc, nvir = cderi_uov.shape - if frac_num is not None: - frac_occ, frac_vir = frac_num[:nocc], frac_num[nocc:] - else: - frac_occ = frac_vir = None - eo = mo_energy[:nocc] - ev = mo_energy[nocc:] - - # loops - log.info("[INFO] Start RI-MP2 loop") - nbatch = util.calc_batch_size(4 * nocc * nvir ** 2, max_memory, dtype=cderi_uov.dtype) - eng_bi1 = eng_bi2 = 0 - for sI in util.gen_batch(0, nocc, nbatch): - log.info(f"[INFO] MP2 loop i: {sI}") - if cderi_uov_2 is None: - g_Iajb = lib.einsum("PIa, Pjb -> Iajb", cderi_uov[:, sI], cderi_uov) - else: - g_Iajb = 0.5 * lib.einsum("PIa, Pjb -> Iajb", cderi_uov[:, sI], cderi_uov_2) - g_Iajb += 0.5 * lib.einsum("PIa, Pjb -> Iajb", cderi_uov_2[:, sI], cderi_uov) - D_Ijab = lib.direct_sum("i + j - a - b -> ijab", eo[sI], eo, ev, ev) - t_Ijab = lib.einsum("Iajb, Ijab -> Ijab", g_Iajb, 1 / D_Ijab) - if t_oovv is not None: - t_oovv[sI] = t_Ijab - if frac_num is not None: - n_Ijab = lib.einsum("i, j, a, b -> ijab", frac_occ[sI], frac_occ, 1 - frac_vir, 1 - frac_vir) - eng_bi1 += lib.einsum("Ijab, Ijab, Iajb ->", n_Ijab, t_Ijab.conj(), g_Iajb) - eng_bi2 += lib.einsum("Ijab, Ijab, Ibja ->", n_Ijab, t_Ijab.conj(), g_Iajb) - else: - eng_bi1 += lib.einsum("Ijab, Iajb ->", t_Ijab.conj(), g_Iajb) - eng_bi2 += lib.einsum("Ijab, Ibja ->", t_Ijab.conj(), g_Iajb) - eng_bi1 = util.check_real(eng_bi1) - eng_bi2 = util.check_real(eng_bi2) - log.info("[INFO] MP2 energy computation finished.") - # report - eng_os = eng_bi1 - eng_ss = eng_bi1 - eng_bi2 - eng_mp2 = eng_os + eng_ss - # results - 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}") - return results - - -class RMP2RI(MP2Base): - """ Restricted MP2 class of doubly hybrid with RI integral. """ - - def __init__(self, mf, frozen=None, omega=0, with_df=None, **kwargs): - super().__init__(mf, frozen=frozen, omega=omega, with_df=with_df, **kwargs) - self.with_df_2 = None - - def driver_eng_mp2(self, **kwargs): - 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_energy_act = self.mo_energy[mask] - frac_num_act = self.frac_num - if frac_num_act is not None: - frac_num_act = frac_num_act[mask] - # prepare t_oovv - max_memory = self.max_memory - lib.current_memory()[0] - t_oovv = util.allocate_array( - self.incore_t_oovv_mp2, (nocc_act, nocc_act, nvir_act, nvir_act), max_memory, - h5file=self._tmpfile, name="t_oovv_mp2", zero_init=False, chunk=(1, 1, nvir_act, nvir_act), - dtype=mo_coeff_act.dtype) - # generate cderi_uov - omega = self.omega - with_df = util.get_with_df_omega(self.with_df, omega) - max_memory = self.max_memory - lib.current_memory()[0] - cderi_uov = self.tensors.get("cderi_uov", None) - if cderi_uov is None: - cderi_uov = util.get_cderi_mo( - with_df, mo_coeff_act, None, (0, nocc_act, nocc_act, nact), max_memory) - self.tensors["cderi_uov"] = cderi_uov - # cderi_uov_2 is rarely called, so do not try to build omega for this special case - cderi_uov_2 = None - max_memory = self.max_memory - lib.current_memory()[0] - if self.with_df_2 is not None: - cderi_uov_2 = util.get_cderi_mo( - self.with_df_2, mo_coeff_act, None, (0, nocc_act, nocc_act, nact), max_memory) - # kernel - max_memory = self.max_memory - lib.current_memory()[0] - results = self.kernel_energy_mp2( - mo_energy_act, cderi_uov, - t_oovv=t_oovv, - frac_num=frac_num_act, - verbose=self.verbose, - max_memory=max_memory, - cderi_uov_2=cderi_uov_2, - **kwargs) - - self.e_corr = results["eng_corr_MP2"] - # pad omega - results = {pad_omega(key, self.omega): val for (key, val) in results.items()} - self.results.update(results) - return results - - kernel_energy_mp2 = staticmethod(kernel_energy_rmp2_ri_incore) - kernel = driver_eng_mp2 - -# endregion - - if __name__ == '__main__': def main_1(): + # test RMP2ConvPySCF 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 = scf.RHF(mol).run() @@ -418,6 +280,7 @@ if __name__ == '__main__': print(mf_mp.results) def main_2(): + # test RMP2Conv 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 = scf.RHF(mol).run() @@ -425,14 +288,5 @@ if __name__ == '__main__': print(mf_mp.e_tot) print(mf_mp.results) - def main_3(): - 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 = scf.RHF(mol).run() - mf_mp = RMP2RI(mf_scf, frozen=[1, 2]).run() - print(mf_mp.e_tot) - print(mf_mp.results) - main_1() main_2() - main_3() diff --git a/pyscf/dh/energy/mp2/rmp2ri.py b/pyscf/dh/energy/mp2/rmp2ri.py new file mode 100644 index 0000000..ed7d168 --- /dev/null +++ b/pyscf/dh/energy/mp2/rmp2ri.py @@ -0,0 +1,633 @@ +import numpy as np +from pyscf import lib, scf, __config__ +from pyscf.dh import util +from pyscf.dh.energy.mp2.rmp2 import MP2Base +from pyscf.dh.util import pad_omega +from pyscf.scf import cphf + +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 kernel_energy_rmp2_ri_incore( + mo_energy, cderi_uov, + c_os, c_ss, + t_oovv=None, frac_num=None, verbose=lib.logger.NOTE, max_memory=2000, cderi_uov_2=None): + """ Kernel of MP2 energy by RI integral. + + Parameters + ---------- + mo_energy : np.ndarray + Molecular orbital energy levels. + cderi_uov : np.ndarray + Cholesky decomposed 3c2e ERI in MO basis (occ-vir part). + c_os : float + Coefficient of oppo-spin contribution. + c_ss : float + Coefficient of same-spin contribution. + + t_oovv : np.ndarray + Store space for ``t_oovv`` + frac_num : np.ndarray + Fractional occupation number list. + verbose : int + Verbose level for PySCF. + max_memory : float + Allocatable memory in MB. + cderi_uov_2 : np.ndarray + Another part of 3c2e ERI in MO basis (occ-vir part). This is mostly used in magnetic computations. + + Notes + ----- + + For RI approximation, ERI integral is set to be + + .. math:: + g_{ij}^{ab} = (ia|jb) = Y_{ia, P} Y_{jb, P} + + For energy of fractional occupation system, computation method is chosen + by eq (7) from Su2016 (10.1021/acs.jctc.6b00197). + """ + log = lib.logger.new_logger(verbose=verbose) + time0 = lib.logger.process_clock(), lib.logger.perf_counter() + log.info("[INFO] Start unrestricted RI-MP2") + + naux, nocc, nvir = cderi_uov.shape + if frac_num is not None: + frac_occ, frac_vir = frac_num[:nocc], frac_num[nocc:] + else: + frac_occ = frac_vir = None + eo = mo_energy[:nocc] + ev = mo_energy[nocc:] + + # loops + log.info("[INFO] Start RI-MP2 loop") + mem_avail = max_memory - lib.current_memory()[0] + nbatch = util.calc_batch_size(4 * nocc * nvir ** 2, mem_avail, dtype=cderi_uov.dtype) + eng_bi1 = eng_bi2 = 0 + for sI in util.gen_batch(0, nocc, nbatch): + log.info(f"[INFO] MP2 loop i: {sI}") + if cderi_uov_2 is None: + g_Iajb = lib.einsum("PIa, Pjb -> Iajb", cderi_uov[:, sI], cderi_uov) + else: + g_Iajb = 0.5 * lib.einsum("PIa, Pjb -> Iajb", cderi_uov[:, sI], cderi_uov_2) + g_Iajb += 0.5 * lib.einsum("PIa, Pjb -> Iajb", cderi_uov_2[:, sI], cderi_uov) + D_Ijab = lib.direct_sum("i + j - a - b -> ijab", eo[sI], eo, ev, ev) + t_Ijab = lib.einsum("Iajb, Ijab -> Ijab", g_Iajb, 1 / D_Ijab) + if t_oovv is not None: + t_oovv[sI] = t_Ijab + if frac_num is not None: + n_Ijab = lib.einsum("i, j, a, b -> ijab", frac_occ[sI], frac_occ, 1 - frac_vir, 1 - frac_vir) + eng_bi1 += lib.einsum("Ijab, Ijab, Iajb ->", n_Ijab, t_Ijab.conj(), g_Iajb) + eng_bi2 += lib.einsum("Ijab, Ijab, Ibja ->", n_Ijab, t_Ijab.conj(), g_Iajb) + else: + eng_bi1 += lib.einsum("Ijab, Iajb ->", t_Ijab.conj(), g_Iajb) + eng_bi2 += lib.einsum("Ijab, Ibja ->", t_Ijab.conj(), g_Iajb) + eng_bi1 = util.check_real(eng_bi1) + eng_bi2 = util.check_real(eng_bi2) + log.info("[INFO] MP2 energy computation finished.") + # report + eng_os = eng_bi1 + eng_ss = eng_bi1 - eng_bi2 + eng_mp2 = c_os * eng_os + c_ss * eng_ss + # results + 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}") + log.timer("kernel_energy_rmp2_ri_incore", *time0) + return results + + +def get_rdm1_corr( + mo_energy, cderi_uov, t_oovv, + c_os, c_ss, + verbose=lib.logger.NOTE, max_memory=2000): + r""" MP2 correlation part of density matrix. + + .. math:: + D_{jk} &= T_{ij}^{ab} t_{ik}^{ab} \\ + D_{ab} &= T_{ij}^{ac} t_{ij}^{bc} + + Parameters + ---------- + mo_energy : np.ndarray + cderi_uov : np.ndarray + t_oovv : np.ndarray + c_os : float + c_ss : float + verbose : int + max_memory : int or float + + 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 len(mo_energy) == nmo + assert t_oovv.shape == (nocc, nocc, nvir, nvir) + + # prepare essential matrices + so, sv = slice(0, nocc), slice(nocc, nmo) + + # allocate results + rdm1_corr = np.zeros((nmo, nmo)) + + def load_slice(slc): + return slc, t_oovv[slc] + + mem_avail = max_memory - lib.current_memory()[0] + nbatch = util.calc_batch_size(2 * nocc * nvir**2, mem_avail) + for sI, t_Oovv in lib.map_with_prefetch(load_slice, util.gen_batch(0, nocc, nbatch)): + log.debug(f"[INFO] Generation of MP2 rdm1 oo vv and G_uov, {sI} of total {nocc} occ orbitals") + T_Oovv = util.restricted_biorthogonalize(t_Oovv, 1, c_os, c_ss) + rdm1_corr[sv, sv] += 2 * lib.einsum("ijac, ijbc -> ab", T_Oovv, t_Oovv) + rdm1_corr[so, so] -= 2 * lib.einsum("ijab, ikab -> jk", T_Oovv, t_Oovv) + + log.timer("get_rdm1_corr", *time0) + tensors = {"rdm1_corr": rdm1_corr} + return tensors + + +def get_G_uov( + mo_energy, cderi_uov, t_oovv, + c_os, c_ss, + verbose=lib.logger.NOTE, max_memory=2000): + r""" Get 3-index transformed MP2 amplitude. + + .. math:: + \Gamma_{ia, P} = T_{ij}^{ab} Y_{jb, P} + + Parameters + ---------- + mo_energy : np.ndarray + cderi_uov : np.ndarray + t_oovv : np.ndarray + c_os : float + c_ss : float + verbose : int + max_memory : int or float + + 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 len(mo_energy) == nmo + assert t_oovv.shape == (nocc, nocc, nvir, nvir) + + # prepare essential matrices + so, sv = slice(0, nocc), slice(nocc, nmo) + + # allocate results + rdm1_corr = np.zeros((nmo, nmo)) + G_uov = np.zeros((naux, nocc, nvir)) + + def load_slice(slc): + return slc, t_oovv[slc] + + mem_avail = max_memory - lib.current_memory()[0] + nbatch = util.calc_batch_size(2 * nocc * nvir**2 + naux * nvir, mem_avail) + for sI, t_Oovv in lib.map_with_prefetch(load_slice, util.gen_batch(0, nocc, nbatch)): + log.debug(f"[INFO] Generation of MP2 rdm1 oo vv and G_uov, {sI} of total {nocc} occ orbitals") + T_Oovv = util.restricted_biorthogonalize(t_Oovv, 1, c_os, c_ss) + rdm1_corr[sv, sv] += 2 * lib.einsum("ijac, ijbc -> ab", T_Oovv, t_Oovv) + rdm1_corr[so, so] -= 2 * lib.einsum("ijab, ikab -> jk", T_Oovv, t_Oovv) + G_uov[:, sI] += lib.einsum("ijab, Pjb -> Pia", T_Oovv, cderi_uov) + + log.timer("get_G_uov", *time0) + tensors = { + "rdm1_corr": rdm1_corr, + "G_uov": G_uov, + } + return tensors + + +def get_W_I(cderi_uov, cderi_uoo, G_uov, verbose=lib.logger.NOTE): + + 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): + + 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): + + 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 + + # symmetrize + rdm1_corr_resp = 0.5 * (rdm1_corr_resp + rdm1_corr_resp.T) + + log.timer("get_rdm1_corr_resp", *time0) + tensors = {"rdm1_corr_resp": rdm1_corr_resp} + return tensors + + +class RMP2RI(MP2Base): + """ Restricted MP2 class of doubly hybrid with RI integral. """ + + def __init__(self, mf, **kwargs): + super().__init__(mf, **kwargs) + self.with_df_2 = None + self.incore_cderi_uaa = CONFIG_incore_cderi_uaa_mp2 + self.max_cycle_cpks = CONFIG_max_cycle_cpks + self.tol_cpks = CONFIG_tol_cpks + self._Ax0_Core = NotImplemented + + def make_cderi_uov(self): + """ Generate cholesky decomposed ERI (in memory, occ-vir block). """ + 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] + omega = self.omega + with_df = util.get_with_df_omega(self.with_df, omega) + cderi_uov = util.get_cderi_mo( + with_df, mo_coeff_act, None, (0, nocc_act, nocc_act, nact), + max_memory=self.max_memory, verbose=self.verbose) + + tensors = {"cderi_uov": cderi_uov} + self.tensors.update(tensors) + return cderi_uov + + def make_cderi_uoo(self): + """ Generate cholesky decomposed ERI (in memory, occ-occ block). """ + mask = self.get_frozen_mask() + idx_occ = (self.mo_occ > 0) & mask + mo_occ_act = self.mo_coeff[:, idx_occ] + nocc_act = idx_occ.sum() + omega = self.omega + with_df = util.get_with_df_omega(self.with_df, omega) + cderi_uoo = util.get_cderi_mo( + with_df, mo_occ_act, None, (0, nocc_act, 0, nocc_act), + max_memory=self.max_memory, verbose=self.verbose) + + tensors = {"cderi_uoo": cderi_uoo} + self.tensors.update(tensors) + return cderi_uoo + + def make_cderi_uaa(self): + """ Generate cholesky decomposed ERI (all block, s1 symm, in memory/disk). """ + 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 driver_eng_mp2(self, **kwargs): + """ Driver of MP2 energy. """ + 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_energy_act = self.mo_energy[mask] + frac_num_act = self.frac_num + if frac_num_act is not None: + frac_num_act = frac_num_act[mask] + # prepare t_oovv + max_memory = self.max_memory - lib.current_memory()[0] + t_oovv = util.allocate_array( + self.incore_t_oovv_mp2, (nocc_act, nocc_act, nvir_act, nvir_act), max_memory, + h5file=self._tmpfile, name="t_oovv_mp2", zero_init=False, chunk=(1, 1, nvir_act, nvir_act), + dtype=mo_coeff_act.dtype) + if t_oovv is not None: + self.tensors["t_oovv"] = t_oovv + # generate cderi_uov + cderi_uov = self.tensors.get("cderi_uov", self.make_cderi_uov()) + # cderi_uov_2 is rarely called, so do not try to build omega for this special case + cderi_uov_2 = None + if self.with_df_2 is not None: + cderi_uov_2 = util.get_cderi_mo( + self.with_df_2, mo_coeff_act, None, (0, nocc_act, nocc_act, nact), + max_memory=self.max_memory, verbose=self.verbose) + # kernel + results = self.kernel_energy_mp2( + mo_energy_act, cderi_uov, + self.c_os, self.c_ss, + t_oovv=t_oovv, + frac_num=frac_num_act, + verbose=self.verbose, + max_memory=self.max_memory, + cderi_uov_2=cderi_uov_2, + **kwargs) + + self.e_corr = results["eng_corr_MP2"] + # pad omega + results = {pad_omega(key, self.omega): val for (key, val) in results.items()} + self.results.update(results) + return results + + @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 import RHDFT, UHDFT + HDFT = RHDFT if restricted else UHDFT + mf_scf = HDFT(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 + + def make_rdm1_corr(self): + log = lib.logger.new_logger(verbose=self.verbose) + mask = self.get_frozen_mask() + mo_energy = self.mo_energy[mask] + cderi_uov = self.tensors["cderi_uov"] + t_oovv = self.tensors.get("t_oovv", None) + if t_oovv is None: + log.warn("t_oovv for MP2 has not been stored. " + "To perform G_uov, t_oovv and MP2 energy is to be re-evaluated.") + self.incore_t_oovv_mp2 = True + self.driver_eng_mp2() + t_oovv = self.tensors["t_oovv"] + tensors = get_rdm1_corr( + mo_energy, cderi_uov, t_oovv, self.c_ss, self.c_os, + verbose=self.verbose, max_memory=self.max_memory) + self.tensors.update(tensors) + return tensors["rdm1_corr"] + + def make_G_uov(self): + log = lib.logger.new_logger(verbose=self.verbose) + mask = self.get_frozen_mask() + mo_energy = self.mo_energy[mask] + cderi_uov = self.tensors["cderi_uov"] + t_oovv = self.tensors.get("t_oovv", None) + if t_oovv is None: + log.warn("t_oovv for MP2 has not been stored. " + "To perform G_uov, t_oovv and MP2 energy is to be re-evaluated.") + self.incore_t_oovv_mp2 = True + self.driver_eng_mp2() + t_oovv = self.tensors["t_oovv"] + tensors = get_G_uov( + mo_energy, cderi_uov, t_oovv, self.c_ss, self.c_os, + verbose=self.verbose, max_memory=self.max_memory) + self.tensors.update(tensors) + return tensors["G_uov"] + + def make_W_I(self): + log = lib.logger.new_logger(verbose=self.verbose) + + # prepare tensors + if any([key not in self.tensors for key in ["rdm1_corr", "G_uov", "t_oovv"]]): + log.warn("Some tensors (rdm1_corr, G_uov, t_oovv) has not been generated. " + "Perform make_mp2_integrals first.") + self.make_G_uov() + cderi_uov = self.tensors["cderi_uov"] + cderi_uoo = self.tensors.get("cderi_uoo", self.make_cderi_uoo()) + G_uov = self.tensors["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): + # 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): + # 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): + # 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): + # 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 + + kernel_energy_mp2 = staticmethod(kernel_energy_rmp2_ri_incore) + kernel = driver_eng_mp2 + + +if __name__ == '__main__': + + def main_1(): + # test RMP2RI + 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 = scf.RHF(mol).run() + mf_mp = RMP2RI(mf_scf, frozen=[1, 2]).run() + print(mf_mp.e_tot) + print(mf_mp.results) + + def main_2(): + # test response rdm1 generation + from pyscf import gto, scf + np.set_printoptions(4, 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_scf = scf.RHF(mol).run() + mf_mp = RMP2RI(mf_scf, frozen=[2]).run(incore_t_oovv_mp2=True) + mf_mp.make_rdm1() + for key, val in mf_mp.tensors.items(): + print(key, val.shape) + print(mf_mp.tensors["rdm1"]) + + def main_3(): + # 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 = RMP2RI(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_3() diff --git a/pyscf/dh/energy/mp2/ump2.py b/pyscf/dh/energy/mp2/ump2.py index dd0d739..af1820a 100644 --- a/pyscf/dh/energy/mp2/ump2.py +++ b/pyscf/dh/energy/mp2/ump2.py @@ -2,7 +2,8 @@ r""" Unrestricted MP2. """ from pyscf.dh import util -from pyscf.dh.energy.mp2.rmp2 import RMP2ConvPySCF, RMP2Conv, RMP2RI +from pyscf.dh.energy.mp2.rmp2 import RMP2ConvPySCF, RMP2Conv +from pyscf.dh.energy.mp2.rmp2ri import RMP2RI from pyscf import ao2mo, lib, mp, __config__, df import numpy as np @@ -39,6 +40,7 @@ class UMP2ConvPySCF(mp.ump2.UMP2, RMP2ConvPySCF): def kernel_energy_ump2_conv_full_incore( mo_energy, mo_coeff, eri_or_mol, mo_occ, + c_os, c_ss, t_oovv=None, frac_num=None, verbose=lib.logger.NOTE): """ Kernel of unrestricted MP2 energy by conventional method. @@ -50,6 +52,10 @@ def kernel_energy_ump2_conv_full_incore( Molecular coefficients. eri_or_mol : np.ndarray or gto.Mole ERI that is recognized by ``pyscf.ao2mo.general``. + c_os : float + Coefficient of oppo-spin contribution. + c_ss : float + Coefficient of same-spin contribution. t_oovv : list[np.ndarray] Store space for ``t_oovv`` @@ -120,7 +126,7 @@ def kernel_energy_ump2_conv_full_incore( eng_spin[0] *= 0.25 eng_spin[2] *= 0.25 eng_spin = util.check_real(eng_spin) - eng_mp2 = eng_spin[1] + (eng_spin[0] + eng_spin[2]) + eng_mp2 = c_os * eng_spin[1] + c_ss * (eng_spin[0] + eng_spin[2]) # finalize results results = dict() results["eng_corr_MP2_aa"] = eng_spin[0] @@ -161,7 +167,7 @@ class UMP2Conv(RMP2Conv): for s0, s1, ss, ssn in ((0, 0, 0, "aa"), (0, 1, 1, "ab"), (1, 1, 2, "bb")): t_oovv[ss] = util.allocate_array( incore_t_oovv, shape=(nocc_act[s0], nocc_act[s1], nvir_act[s0], nvir_act[s1]), - max_memory=max_memory, + mem_avail=max_memory, h5file=self._tmpfile, name=f"t_oovv_{ssn}", dtype=mo_coeff_act[0].dtype) @@ -172,6 +178,7 @@ class UMP2Conv(RMP2Conv): eri_or_mol = eri_or_mol if eri_or_mol is not None else self.mol results = self.kernel_energy_mp2( mo_energy_act, mo_coeff_act, eri_or_mol, mo_occ_act, + self.c_os, self.c_ss, t_oovv=t_oovv, frac_num=frac_num_act, verbose=self.verbose, **kwargs) self.e_corr = results["eng_corr_MP2"] # pad omega @@ -188,7 +195,7 @@ class UMP2Conv(RMP2Conv): # region UMP2RI def kernel_energy_ump2_ri_incore( - mo_energy, cderi_uov, + mo_energy, cderi_uov, c_os, c_ss, t_oovv=None, frac_num=None, verbose=lib.logger.NOTE, max_memory=2000, cderi_uov_2=None): """ Kernel of unrestricted MP2 energy by RI integral. @@ -203,6 +210,10 @@ def kernel_energy_ump2_ri_incore( Molecular orbital energy levels. cderi_uov : list[np.ndarray] Cholesky decomposed 3c2e ERI in MO basis (occ-vir part). Spin in (aa, bb). + c_os : float + Coefficient of oppo-spin contribution. + c_ss : float + Coefficient of same-spin contribution. t_oovv : list[np.ndarray] Store space for ``t_oovv`` @@ -269,7 +280,7 @@ def kernel_energy_ump2_ri_incore( eng_spin[0] *= 0.25 eng_spin[2] *= 0.25 eng_spin = util.check_real(eng_spin) - eng_mp2 = eng_spin[1] + (eng_spin[0] + eng_spin[2]) + eng_mp2 = c_os * eng_spin[1] + c_ss * (eng_spin[0] + eng_spin[2]) # finalize results results = dict() results["eng_corr_MP2_aa"] = eng_spin[0] @@ -310,7 +321,7 @@ class UMP2RI(RMP2RI): for s0, s1, ss, ssn in ((0, 0, 0, "aa"), (0, 1, 1, "ab"), (1, 1, 2, "bb")): t_oovv[ss] = util.allocate_array( incore_t_oovv, shape=(nocc_act[s0], nocc_act[s1], nvir_act[s0], nvir_act[s1]), - max_memory=max_memory, + mem_avail=max_memory, h5file=self._tmpfile, name=f"t_oovv_{ssn}", dtype=mo_coeff_act[0].dtype) @@ -338,6 +349,7 @@ class UMP2RI(RMP2RI): max_memory = self.max_memory - lib.current_memory()[0] results = self.kernel_energy_mp2( mo_energy_act, cderi_uov, + self.c_os, self.c_ss, t_oovv=t_oovv, frac_num=frac_num_act, verbose=self.verbose, diff --git a/pyscf/dh/energy/rhdft.py b/pyscf/dh/energy/rhdft.py index d79cf01..3cdce2c 100644 --- a/pyscf/dh/energy/rhdft.py +++ b/pyscf/dh/energy/rhdft.py @@ -313,7 +313,7 @@ def numint_customized(xc, _mol=None): 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}`. + :math:`\sum_{rs} A_{pq, rs} X_{rs}^\mathbb{A}` by PySCF's response function. Parameters ---------- @@ -430,9 +430,15 @@ 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 None: + if hasattr(mf, "xc"): + self.xc = xc + else: + self.xc = "HF" # not KS object + else: + self.xc = xc # type: XCList if isinstance(self.xc, str): xc_scf = XCList(self.xc, code_scf=True) xc_eng = XCList(self.xc, code_scf=False) @@ -485,10 +491,34 @@ class RHDFT(EngBase): """ Fock response function (derivative w.r.t. molecular coefficient in AO basis). """ return self.scf.gen_response() - def Ax0_Core_resp(self, sp, sq, sr, ss, vresp=None, mo_coeff=None): + 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 diff --git a/pyscf/dh/util/df_addon.py b/pyscf/dh/util/df_addon.py index b05f610..c7e6ad2 100644 --- a/pyscf/dh/util/df_addon.py +++ b/pyscf/dh/util/df_addon.py @@ -1,11 +1,13 @@ import numpy as np import copy -from pyscf import df, lib +from pyscf import df, lib, __config__ 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 @@ -19,13 +21,18 @@ def get_cderi_mo(with_df, mo_coeff, Y_mo=None, pqslice=None, max_memory=2000): pqslice : tuple[int] Slice of orbitals to be transformed. max_memory : float - Maximum memory available. + Maximum memory of molecular object. + verbose : int + Print level 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: @@ -38,13 +45,16 @@ def get_cderi_mo(with_df, mo_coeff, Y_mo=None, pqslice=None, max_memory=2000): p0, p1 = 0, 0 preflop = 0 if not isinstance(Y_mo, np.ndarray) else Y_mo.size - nbatch = calc_batch_size(2*nump*numq, max_memory, preflop) + mem_avail = max_memory - lib.current_memory()[0] + nbatch = calc_batch_size(2*nump*numq, mem_avail, preflop) + log.debug(f"[DEBUG] Number of available batch size in get_cderi_mo: {nbatch}") if hasattr(with_df._cderi, "shape"): # array alike aosym = "s2" if len(with_df._cderi.shape) == 2 else "s1" else: # assert pyscf when memory is low aosym = "s2" for Y_ao in with_df.loop(nbatch): p1 = p0 + Y_ao.shape[0] + log.debug(f"[DEBUG] cderi transformation of auxiliary basis: {p0, p1}") if Y_ao.dtype == np.double and mo_coeff.dtype == np.double: Y_mo[p0:p1] = _ao2mo.nr_e2(Y_ao, mo_coeff, pqslice, aosym=aosym, mosym="s1").reshape(p1 - p0, nump, numq) else: @@ -55,6 +65,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_cderi_mo", *time0) return Y_mo diff --git a/pyscf/dh/util/frozen_core.py b/pyscf/dh/util/frozen_core.py index 5621431..420cedb 100644 --- a/pyscf/dh/util/frozen_core.py +++ b/pyscf/dh/util/frozen_core.py @@ -875,10 +875,13 @@ class FrozenCore: mask = np.ones((nset, nmo), dtype=bool) if isinstance(self.rule, list or np.ndarray): + arr = self.rule + # 0. frozen set is empty + if len(arr) == 0: + mask[:] = True # 1. case of frozen/active list defined, highest priority # first perform mask; if self.act, then reverse mask in order to make active orbitals to be True - arr = self.rule - if not hasattr(arr[0], "__iter__"): + elif not hasattr(arr[0], "__iter__"): # one list of frozen orbitals: mask[:, arr] = False else: diff --git a/pyscf/dh/util/helper.py b/pyscf/dh/util/helper.py index b47938d..c1afd2a 100644 --- a/pyscf/dh/util/helper.py +++ b/pyscf/dh/util/helper.py @@ -86,7 +86,7 @@ def parse_incore_flag(flag, unit_flop, mem_avail, pre_flop=0, dtype=float): Memory available in MB. pre_flop : int Number of data preserved in memory. Unit in number. - dtype : type + dtype : np.dtype Type of data. Such as np.float64, complex, etc. Returns @@ -217,7 +217,7 @@ def pad_omega(s, omega): return f"{s}_omega({omega:.6f})" -def allocate_array(incore, shape, max_memory, +def allocate_array(incore, shape, mem_avail, h5file=None, name=None, dtype=float, zero_init=True, chunk=None, **kwargs): """ Allocate an array with given memory estimation and incore stragety. @@ -227,8 +227,8 @@ def allocate_array(incore, shape, max_memory, Incore flag. Also see :class:`parse_incore_flag`. shape : tuple Shape of allocated array. - max_memory : int or float - Maximum memory in MB. + mem_avail : int or float + Memory available for allocation. h5file : h5py.File HDF5 file instance. Array may save into this file if disk is required. name : str @@ -244,7 +244,7 @@ def allocate_array(incore, shape, max_memory, ------- np.array or h5py.Dataset """ - incore = parse_incore_flag(incore, int(np.prod(shape)), max_memory, dtype=dtype) + incore = parse_incore_flag(incore, int(np.prod(shape)), mem_avail, dtype=dtype) if incore is None: return None elif incore is True: -- Gitee From a85bfc2879040abbd0bea6aa893af40759fdad0a Mon Sep 17 00:00:00 2001 From: ajz34 Date: Thu, 27 Apr 2023 18:43:41 +0800 Subject: [PATCH 04/12] docs(dh eng): add docstring --- pyscf/dh/energy/mp2/rmp2ri.py | 104 ++++++++++++++++++++++++++++++---- 1 file changed, 93 insertions(+), 11 deletions(-) diff --git a/pyscf/dh/energy/mp2/rmp2ri.py b/pyscf/dh/energy/mp2/rmp2ri.py index ed7d168..dd9ddcb 100644 --- a/pyscf/dh/energy/mp2/rmp2ri.py +++ b/pyscf/dh/energy/mp2/rmp2ri.py @@ -109,11 +109,13 @@ def get_rdm1_corr( mo_energy, cderi_uov, t_oovv, c_os, c_ss, verbose=lib.logger.NOTE, max_memory=2000): - r""" MP2 correlation part of density matrix. + r""" 1-RDM of MP2 correlation by MO basis. .. math:: - D_{jk} &= T_{ij}^{ab} t_{ik}^{ab} \\ - D_{ab} &= T_{ij}^{ac} t_{ij}^{bc} + D_{jk}^\mathrm{(2)} &= 2 T_{ij}^{ab} t_{ik}^{ab} \\ + D_{ab}^\mathrm{(2)} &= 2 T_{ij}^{ac} t_{ij}^{bc} + + By definition of this program, 1-RDM is not response density. Parameters ---------- @@ -166,6 +168,8 @@ def get_G_uov( verbose=lib.logger.NOTE, max_memory=2000): r""" Get 3-index transformed MP2 amplitude. + Evaluation of this 3-index amplitude is performed together with non-response MP2 correlation 1-RDM. + .. math:: \Gamma_{ia, P} = T_{ij}^{ab} Y_{jb, P} @@ -220,6 +224,26 @@ def get_G_uov( 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() @@ -246,6 +270,34 @@ def get_W_I(cderi_uov, cderi_uoo, G_uov, verbose=lib.logger.NOTE): 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 + cderi_uaa + W_I + rdm1_corr + Ax0_Core + max_memory + verbose + + Returns + ------- + dict[str, np.ndarray] + """ log = lib.logger.new_logger(verbose=verbose) time0 = lib.logger.process_clock(), lib.logger.perf_counter() @@ -278,6 +330,32 @@ 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() @@ -301,9 +379,6 @@ def get_rdm1_corr_resp( max_cycle=max_cycle, tol=tol)[0] rdm1_corr_resp[sv, so] = rdm1_corr_resp_vo - # symmetrize - rdm1_corr_resp = 0.5 * (rdm1_corr_resp + rdm1_corr_resp.T) - log.timer("get_rdm1_corr_resp", *time0) tensors = {"rdm1_corr_resp": rdm1_corr_resp} return tensors @@ -321,7 +396,7 @@ class RMP2RI(MP2Base): self._Ax0_Core = NotImplemented def make_cderi_uov(self): - """ Generate cholesky decomposed ERI (in memory, occ-vir block). """ + r""" Generate cholesky decomposed ERI :math:`Y_{ia, P}` (in memory, occ-vir block). """ mask = self.get_frozen_mask() nocc_act = (mask & (self.mo_occ != 0)).sum() nvir_act = (mask & (self.mo_occ == 0)).sum() @@ -338,7 +413,7 @@ class RMP2RI(MP2Base): return cderi_uov def make_cderi_uoo(self): - """ Generate cholesky decomposed ERI (in memory, occ-occ block). """ + r""" Generate cholesky decomposed ERI math:`Y_{ij, P}` (in memory, occ-occ block). """ mask = self.get_frozen_mask() idx_occ = (self.mo_occ > 0) & mask mo_occ_act = self.mo_coeff[:, idx_occ] @@ -354,7 +429,7 @@ class RMP2RI(MP2Base): return cderi_uoo def make_cderi_uaa(self): - """ Generate cholesky decomposed ERI (all block, s1 symm, in memory/disk). """ + r""" Generate cholesky decomposed ERI math:`Y_{pq, P}` (all block, s1 symm, in memory/disk). """ log = lib.logger.new_logger(verbose=self.verbose) # dimension and mask @@ -385,7 +460,7 @@ class RMP2RI(MP2Base): return cderi_uaa def driver_eng_mp2(self, **kwargs): - """ Driver of MP2 energy. """ + r""" Driver of MP2 energy. """ mask = self.get_frozen_mask() nocc_act = (mask & (self.mo_occ != 0)).sum() nvir_act = (mask & (self.mo_occ == 0)).sum() @@ -430,7 +505,7 @@ class RMP2RI(MP2Base): @property def Ax0_Core(self): - """ Fock response of underlying SCF object in MO basis. """ + r""" 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 import RHDFT, UHDFT @@ -444,6 +519,7 @@ class RMP2RI(MP2Base): self._Ax0_Core = Ax0_Core def make_rdm1_corr(self): + r""" Generate 1-RDM (non-response) of MP2 correlation :math:`D_{pq}^{(2)}`. """ log = lib.logger.new_logger(verbose=self.verbose) mask = self.get_frozen_mask() mo_energy = self.mo_energy[mask] @@ -462,6 +538,7 @@ class RMP2RI(MP2Base): return tensors["rdm1_corr"] def make_G_uov(self): + r""" Generate 3-index transformed MP2 amplitude :math:`\Gamma_{ia, P}`. """ log = lib.logger.new_logger(verbose=self.verbose) mask = self.get_frozen_mask() mo_energy = self.mo_energy[mask] @@ -480,6 +557,7 @@ class RMP2RI(MP2Base): return tensors["G_uov"] def make_W_I(self): + r""" Generate part I of MO-energy-weighted density matrix. :math:`W_{pq} [\mathrm{I}]`. """ log = lib.logger.new_logger(verbose=self.verbose) # prepare tensors @@ -497,6 +575,7 @@ class RMP2RI(MP2Base): return tensors["W_I"] def make_lag_vo(self): + r""" Generate MP2 contribution to Lagrangian vir-occ block :math:`L_{ai}`. """ # prepare tensors W_I = self.tensors.get("W_I", self.make_W_I()) cderi_uaa = self.tensors.get("cderi_uaa", self.make_cderi_uaa()) @@ -513,6 +592,7 @@ class RMP2RI(MP2Base): return tensors["lag_vo"] def make_rdm1_corr_resp(self): + """ Generate 1-RDM (response) of MP2 correlation :math:`D_{pq}^{(2)}`. """ # prepare input lag_vo = self.tensors.get("lag_vo", self.make_lag_vo()) rdm1_corr = self.tensors["rdm1_corr"] @@ -533,6 +613,7 @@ class RMP2RI(MP2Base): 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()) @@ -546,6 +627,7 @@ class RMP2RI(MP2Base): 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()) -- Gitee From c7d89727a07769decaefe57331e696dee71bc5a1 Mon Sep 17 00:00:00 2001 From: ajz34 Date: Wed, 3 May 2023 09:53:50 +0800 Subject: [PATCH 05/12] docs(dh eng): minor edit --- pyscf/dh/energy/mp2/rmp2ri.py | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/pyscf/dh/energy/mp2/rmp2ri.py b/pyscf/dh/energy/mp2/rmp2ri.py index dd9ddcb..8b21db3 100644 --- a/pyscf/dh/energy/mp2/rmp2ri.py +++ b/pyscf/dh/energy/mp2/rmp2ri.py @@ -286,13 +286,13 @@ def get_lag_vo( Parameters ---------- - G_uov - cderi_uaa - W_I - rdm1_corr - Ax0_Core - max_memory - verbose + 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 ------- -- Gitee From 14dc4c885d0c5c4b743835325ffe332df6edfe66 Mon Sep 17 00:00:00 2001 From: ajz34 Date: Wed, 3 May 2023 16:33:34 +0800 Subject: [PATCH 06/12] refactor(dh eng): change logic of energy evaluation process and xc splitting --- pyscf/dh/energy/dh.py | 179 +++++++++++++++++++++++++----------------- 1 file changed, 106 insertions(+), 73 deletions(-) diff --git a/pyscf/dh/energy/dh.py b/pyscf/dh/energy/dh.py index e14e2f1..a009a7a 100644 --- a/pyscf/dh/energy/dh.py +++ b/pyscf/dh/energy/dh.py @@ -15,6 +15,48 @@ CONFIG_route_iepa = getattr(__config__, "route_iepa", CONFIG_route) CONFIG_frozen = getattr(__config__, "frozen", 0) +def _process_xc_split(xc_list): + """ Split exch-corr list by certain type of xc. + + Exch-corr list of input should be energy functional. + + Parameters + ---------- + xc_list : XCList + + Returns + ------- + dict[str, XCList] + """ + + result = dict() + # 1. low-rung + xc_extracted = xc_list.copy() + xc_low_rung = xc_extracted.extract_by_xctype(XCType.RUNG_LOW) + xc_extracted = xc_extracted.remove(xc_low_rung) + result["low_rung"] = xc_low_rung + # 2. advanced correlation + # 2.1 IEPA (may also resolve evaluation of MP2) + xc_iepa = xc_extracted.extract_by_xctype(XCType.IEPA) + if len(xc_iepa) != 0: + xc_iepa = xc_extracted.extract_by_xctype(XCType.IEPA | XCType.MP2) + xc_extracted = xc_extracted.remove(xc_iepa) + result["iepa"] = xc_iepa + # 2.2 MP2 + xc_mp2 = xc_extracted.extract_by_xctype(XCType.MP2 | XCType.RSMP2) + xc_extracted = xc_extracted.remove(xc_mp2) + result["mp2"] = xc_mp2 + # 2.3 Ring CCD + xc_ring_ccd = xc_extracted.extract_by_xctype(XCType.RS_RING_CCD) + xc_extracted = xc_extracted.remove(xc_ring_ccd) + result["ring_ccd"] = xc_ring_ccd + + # finalize + if len(xc_extracted) > 0: + raise RuntimeError(f"Some xc terms not evaluated! Possibly bg of program.\nXC not evaluated: {xc_extracted}") + return result + + def _process_energy_exx(mf_dh, xc_list, force_evaluate=False): """ Evaluate exact exchange energy. @@ -61,31 +103,27 @@ def _process_energy_exx(mf_dh, xc_list, force_evaluate=False): eng_tot += eng update_results(mf_dh.results, results) - mf_dh.inherited.append((mf_scf, xc_exx)) return xc_extracted, eng_tot -def _process_energy_mp2(mf_dh, xc_list, force_evaluate=False): +def _process_energy_mp2(mf_dh, xc_mp2, force_evaluate=False): """ Evaluate MP2 correlation energy. Parameters ---------- mf_dh : DH - xc_list : XCList + xc_mp2 : XCList force_evaluate : bool Returns ------- tuple[XCList, float] """ - log = lib.logger.new_logger(verbose=mf_dh.verbose) - xc_mp2 = xc_list.extract_by_xctype(XCType.MP2 | XCType.RSMP2) if len(xc_mp2) == 0: - return xc_list, 0 - xc_extracted = xc_list.remove(xc_mp2, inplace=False) - log.info(f"[INFO] XCList extracted by process_energy_mp2: {xc_mp2.token}") - log.info(f"[INFO] XCList remains by process_energy_mp2: {xc_extracted.token}") - log.info("[INFO] MP2 detected") + return 0 + + log = lib.logger.new_logger(verbose=mf_dh.verbose) + log.info(f"[INFO] XCList to be evaluated by process_energy_mp2: {xc_mp2.token}") def comput_mp2(): eng_tot = 0 @@ -96,7 +134,9 @@ def _process_energy_mp2(mf_dh, xc_list, force_evaluate=False): omega = 0 else: omega, c_os, c_ss = info.parameters - eng = info.fac * results[pad_omega("eng_corr_MP2", omega)] + eng = info.fac * ( + + c_os * results[pad_omega("eng_corr_MP2_OS", omega)] + + c_ss * results[pad_omega("eng_corr_MP2_SS", omega)]) log.note(f"[RESULT] Energy of correlation {info.token}: {eng:20.12f}") eng_tot += eng return eng_tot @@ -112,7 +152,7 @@ def _process_energy_mp2(mf_dh, xc_list, force_evaluate=False): continue # results have been evaluated mf_mp2 = mf_dh.to_mp2(omega=omega, c_os=c_os, c_ss=c_ss).run() update_results(mf_dh.results, mf_mp2.results) - mf_dh.inherited.append((mf_mp2, xc_mp2)) + mf_dh.inherited["mp2"][1].append(mf_mp2) eng_tot = comput_mp2() return eng_tot @@ -124,32 +164,28 @@ def _process_energy_mp2(mf_dh, xc_list, force_evaluate=False): except KeyError: eng_tot = force_comput_mp2() - return xc_extracted, eng_tot + return eng_tot -def _process_energy_iepa(mf_dh, xc_list, force_evaluate=False): +def _process_energy_iepa(mf_dh, xc_iepa, force_evaluate=False): """ Evaluate IEPA-like correlation energy. Parameters ---------- mf_dh : DH - xc_list : XCList + xc_iepa : XCList force_evaluate : bool Returns ------- - tuple[XCList, float] + float """ log = lib.logger.new_logger(verbose=mf_dh.verbose) - xc_iepa = xc_list.extract_by_xctype(XCType.IEPA) if len(xc_iepa) == 0: - return xc_list, 0 + return 0 # run MP2 while in IEPA if found - xc_iepa = xc_list.extract_by_xctype(XCType.IEPA | XCType.MP2) - xc_extracted = xc_list.remove(xc_iepa, inplace=False) - log.info(f"[INFO] XCList extracted by process_energy_iepa: {xc_iepa.token}") - log.info(f"[INFO] XCList remains by process_energy_iepa: {xc_extracted.token}") + log.info(f"[INFO] XCList to be evaluated by process_energy_iepa: {xc_iepa.token}") # prepare IEPA iepa_schemes = [info.name for info in xc_iepa] @@ -170,7 +206,7 @@ def _process_energy_iepa(mf_dh, xc_list, force_evaluate=False): mf_iepa = mf_dh.to_iepa().run(iepa_schemes=iepa_schemes) update_results(mf_dh.results, mf_iepa.results) eng_tot = comput_iepa() - mf_dh.inherited.append((mf_iepa, xc_iepa)) + mf_dh.inherited["iepa"][1].append(mf_iepa) return eng_tot if force_evaluate: @@ -181,33 +217,30 @@ def _process_energy_iepa(mf_dh, xc_list, force_evaluate=False): except KeyError: eng_tot = force_comput_iepa() - return xc_extracted, eng_tot + return eng_tot -def _process_energy_drpa(mf_dh, xc_list, force_evaluate=False): - """ Evaluate RPA-like correlation energy. +def _process_energy_ring_ccd(mf_dh, xc_ring_ccd, force_evaluate=False): + """ Evaluate Ring-CCD-like correlation energy. Parameters ---------- mf_dh : DH - xc_list : XCList + xc_ring_ccd : XCList force_evaluate : bool Returns ------- tuple[XCList, float] """ - log = lib.logger.new_logger(verbose=mf_dh.verbose) - xc_ring_ccd = xc_list.extract_by_xctype(XCType.RS_RING_CCD) - xc_extracted = xc_list.remove(xc_ring_ccd, inplace=False) if len(xc_ring_ccd) == 0: - return xc_list, 0 - log.info(f"[INFO] XCList extracted by process_energy_drpa (RS_RING_CCD): {xc_ring_ccd.token}") - log.info(f"[INFO] XCList remains by process_energy_drpa (RS_RING_CCD): {xc_extracted.token}") - log.info(f"[INFO] Ring-CCD detected") + return 0 + + log = lib.logger.new_logger(verbose=mf_dh.verbose) + log.info(f"[INFO] XCList to be evaluated by process_energy_drpa: {xc_ring_ccd.token}") # generate omega list - # parameter of RSMP2: omega, c_os, c_ss + # parameter of RS-Ring-CCD: omega, c_os, c_ss omega_list = [] for xc_info in xc_ring_ccd: omega_list.append(xc_info.parameters[0]) @@ -230,7 +263,7 @@ def _process_energy_drpa(mf_dh, xc_list, force_evaluate=False): continue # results have been evaluated mf_ring_ccd = mf_dh.to_ring_ccd(omega=omega).run() update_results(mf_dh.results, mf_ring_ccd.results) - mf_dh.inherited.append((mf_ring_ccd, xc_ring_ccd)) + mf_dh.inherited["ring_ccd"][1].append(mf_ring_ccd) eng_tot = comput_ring_ccd() return eng_tot @@ -242,21 +275,20 @@ def _process_energy_drpa(mf_dh, xc_list, force_evaluate=False): except KeyError: eng_tot = force_comput_ring_ccd() - return xc_extracted, eng_tot + return eng_tot -def _process_energy_low_rung(mf_dh, xc_list, xc_to_parse=None): +def _process_energy_low_rung(mf_dh, xc_list): """ Evaluate low-rung pure DFT exchange/correlation energy. Parameters ---------- mf_dh : DH xc_list : XCList - xc_to_parse : XCList Returns ------- - tuple[XCList, float] + float """ # logic of this function # 1. first try the whole low_rung token @@ -268,13 +300,10 @@ def _process_energy_low_rung(mf_dh, xc_list, xc_to_parse=None): log = lib.logger.new_logger(verbose=mf_dh.verbose) mf_scf = mf_dh.to_scf() ni = dft.numint.NumInt() - if xc_to_parse is None: - xc_to_parse = xc_list.extract_by_xctype(XCType.RUNG_LOW) - xc_extracted = xc_list.remove(xc_to_parse, inplace=False) - log.info(f"[INFO] XCList to be parsed in process_energy_low_rung: {xc_to_parse.token}") + log.info(f"[INFO] XCList to be parsed in process_energy_low_rung: {xc_list.token}") - if len(xc_to_parse) == 0: - return xc_list, 0 + if len(xc_list) == 0: + return 0 def handle_multiple_omega(): log.warn( @@ -284,8 +313,8 @@ def _process_energy_low_rung(mf_dh, xc_list, xc_to_parse=None): "and our program currently could not handle it.\n" "Anyway, this is purely experimental and use with caution.") eng_tot = 0 - xc_exx = xc_to_parse.extract_by_xctype(XCType.EXX) - xc_non_exx = xc_to_parse.remove(xc_exx, inplace=False) + xc_exx = xc_list.extract_by_xctype(XCType.EXX) + xc_non_exx = xc_list.remove(xc_exx, inplace=False) log.info(f"[INFO] XCList extracted by process_energy_low_rung (handle_multiple_omega, " f"exx): {xc_exx.token}") log.info(f"[INFO] XCList extracted by process_energy_low_rung (handle_multiple_omega, " @@ -300,7 +329,7 @@ def _process_energy_low_rung(mf_dh, xc_list, xc_to_parse=None): f"This is probably bug.") eng_tot += eng_exx + eng_non_exx log.note(f"[RESULT] Energy of process_energy_low_rung (handle_multiple_omega): {eng_tot}") - return xc_extracted, eng_tot + return eng_tot else: log.warn( "Seems pure-DFT part has different values of omega.\n" @@ -312,12 +341,12 @@ def _process_energy_low_rung(mf_dh, xc_list, xc_to_parse=None): assert len(xc_item_extracted) == 0 eng_tot += eng_item log.note(f"[RESULT] Energy of process_energy_low_rung (handle_multiple_omega): {eng_tot}") - return xc_extracted, eng_tot + return eng_tot numint = None # perform the logic of function try: - ni._xc_type(xc_to_parse.token) + ni._xc_type(xc_list.token) except (ValueError, KeyError) as err: if isinstance(err, ValueError) and "Different values of omega" in err.args[0]: return handle_multiple_omega() @@ -325,7 +354,7 @@ def _process_energy_low_rung(mf_dh, xc_list, xc_to_parse=None): or isinstance(err, ValueError) and "too many values to unpack" in err.args[0]: log.info("[INFO] Unknown functional to PySCF. Try build custom numint.") try: - numint = numint_customized(xc_to_parse) + numint = numint_customized(xc_list) except ValueError as err_numint: if "Different values of omega" in err_numint.args[0]: log.warn("We first resolve different omegas, and we later move on numint.") @@ -339,8 +368,8 @@ def _process_energy_low_rung(mf_dh, xc_list, xc_to_parse=None): if numint is None: numint = dft.numint.NumInt() # exx part - xc_exx = xc_to_parse.extract_by_xctype(XCType.EXX) - xc_non_exx = xc_to_parse.remove(xc_exx, inplace=False) + xc_exx = xc_list.extract_by_xctype(XCType.EXX) + xc_non_exx = xc_list.remove(xc_exx, inplace=False) log.info(f"[INFO] XCList extracted by process_energy_low_rung (exx): {xc_exx.token}") log.info(f"[INFO] XCList extracted by process_energy_low_rung (non_exx): {xc_non_exx.token}") xc_exx_extracted, eng_exx = _process_energy_exx(mf_dh, xc_exx) @@ -369,7 +398,7 @@ def _process_energy_low_rung(mf_dh, xc_list, xc_to_parse=None): eng_tot += results["eng_purexc_" + xc_non_exx.token] log.note(f"[RESULT] Energy of process_energy_low_rung (non_exx): {results['eng_purexc_' + xc_non_exx.token]}") log.note(f"[RESULT] Energy of process_energy_low_rung (total): {eng_tot}") - return xc_extracted, eng_tot + return eng_tot def driver_energy_dh(mf_dh, xc=None, force_evaluate=False): @@ -383,6 +412,8 @@ def driver_energy_dh(mf_dh, xc=None, force_evaluate=False): Token of exchange and correlation. """ log = mf_dh.log + + # prepare xc tokens mf_dh.build() if xc is None: xc = mf_dh.xc.xc_eng @@ -390,43 +421,46 @@ def driver_energy_dh(mf_dh, xc=None, force_evaluate=False): xc = XCList(xc, code_scf=False) assert isinstance(xc, XCList) xc = xc.copy() + xc_splitted = _process_xc_split(xc) + for key, xc_list in xc_splitted.items(): + mf_dh.inherited[key] = (xc_list, []) + result = dict() eng_tot = 0. # 1. general xc - xc_extracted = xc.copy() - xc_low_rung = xc_extracted.extract_by_xctype(XCType.RUNG_LOW) - xc_extracted = xc_extracted.remove(xc_low_rung) + xc_low_rung = xc_splitted["low_rung"] if xc_low_rung == mf_dh.xc.xc_scf and not force_evaluate: # If low-rung of xc_eng and xc_scf shares the same xc formula, # then we just use the SCF energy to evaluate low-rung part of total doubly hybrid energy. log.info("[INFO] xc of SCF is the same to xc of energy in rung-low part. Add SCF energy to total energy.") eng_tot += mf_dh.scf.e_tot else: - xc_low_rung_extracted, eng_low_rung = _process_energy_low_rung(mf_dh, xc_low_rung) - assert len(xc_low_rung_extracted) == 0 + eng_low_rung = _process_energy_low_rung(mf_dh, xc_low_rung) eng_tot += eng_low_rung result.update(mf_dh.to_scf().get_energy_noxc(mf_dh.scf, mf_dh.scf.make_rdm1())) eng_tot += result["eng_noxc"] # 2. other correlation # 2.1 IEPA (may also resolve evaluation of MP2) - xc_extracted, eng_iepa = _process_energy_iepa(mf_dh, xc_extracted) + xc_iepa = xc_splitted["iepa"] + eng_iepa = _process_energy_iepa(mf_dh, xc_iepa) eng_tot += eng_iepa # 2.2 MP2 - xc_extracted, eng_mp2 = _process_energy_mp2(mf_dh, xc_extracted) + xc_mp2 = xc_splitted["mp2"] + eng_mp2 = _process_energy_mp2(mf_dh, xc_mp2) eng_tot += eng_mp2 - # 2.3 dRPA - xc_extracted, eng_drpa = _process_energy_drpa(mf_dh, xc_extracted) - eng_tot += eng_drpa + # 2.3 Ring CCD + xc_ring_ccd = xc_splitted["ring_ccd"] + eng_ring_ccd = _process_energy_ring_ccd(mf_dh, xc_ring_ccd) + eng_tot += eng_ring_ccd # # 2.4 VV10 # xc_extracted, eng_vdw = _process_energy_vdw(mf_dh, xc_extracted) # eng_tot += eng_vdw - # finalize - if len(xc_extracted) > 0: - raise RuntimeError("Some xc terms not evaluated! Possibly bug of program.") + result[f"eng_dh_{xc.token}"] = eng_tot + update_results(mf_dh.results, result) log.note(f"[RESULT] Energy of {xc.token}: {eng_tot:20.12f}") - return result + return eng_tot class DH(EngBase): @@ -452,7 +486,7 @@ class DH(EngBase): # cheat to generate someting by __init__ from base class mol = mf_or_mol if isinstance(mf_or_mol, gto.Mole) else mf_or_mol.mol super().__init__(scf.HF(mol)) - self.inherited = [] # type: List[Tuple[EngBase, XCList]] + self.inherited = dict() # type: dict[str, tuple[XCList, list[EngBase]]] self.xc = NotImplemented # type: XCDH self.log = NotImplemented # type: lib.logger.Logger self.flags = flags if flags is not None else dict() @@ -578,9 +612,8 @@ class DH(EngBase): self.flags.update(kwargs) force_evaluate = self.flags.get("force_evaluate", False) - results = driver_energy_dh(self, xc, force_evaluate) - update_results(self.results, results) - return self.results[f"eng_dh_{xc.token}"] + eng_tot = driver_energy_dh(self, xc, force_evaluate) + return eng_tot @property def e_tot(self): -- Gitee From e037bc54a624aa04da4e45f27b67b3bbebc2809b Mon Sep 17 00:00:00 2001 From: ajz34 Date: Fri, 5 May 2023 14:07:56 +0800 Subject: [PATCH 07/12] Revert "docs(dh eng): minor edit" This reverts commit c7d89727a07769decaefe57331e696dee71bc5a1. --- pyscf/dh/energy/mp2/rmp2ri.py | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/pyscf/dh/energy/mp2/rmp2ri.py b/pyscf/dh/energy/mp2/rmp2ri.py index 8b21db3..dd9ddcb 100644 --- a/pyscf/dh/energy/mp2/rmp2ri.py +++ b/pyscf/dh/energy/mp2/rmp2ri.py @@ -286,13 +286,13 @@ def get_lag_vo( 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 + G_uov + cderi_uaa + W_I + rdm1_corr + Ax0_Core + max_memory + verbose Returns ------- -- Gitee From 5fd85dead19bdee86a2c47e3f416b362fdc8a818 Mon Sep 17 00:00:00 2001 From: ajz34 Date: Fri, 5 May 2023 14:07:56 +0800 Subject: [PATCH 08/12] Revert "docs(dh eng): add docstring" This reverts commit a85bfc2879040abbd0bea6aa893af40759fdad0a. --- pyscf/dh/energy/mp2/rmp2ri.py | 104 ++++------------------------------ 1 file changed, 11 insertions(+), 93 deletions(-) diff --git a/pyscf/dh/energy/mp2/rmp2ri.py b/pyscf/dh/energy/mp2/rmp2ri.py index dd9ddcb..ed7d168 100644 --- a/pyscf/dh/energy/mp2/rmp2ri.py +++ b/pyscf/dh/energy/mp2/rmp2ri.py @@ -109,13 +109,11 @@ def get_rdm1_corr( mo_energy, cderi_uov, t_oovv, c_os, c_ss, verbose=lib.logger.NOTE, max_memory=2000): - r""" 1-RDM of MP2 correlation by MO basis. + r""" MP2 correlation part of density matrix. .. math:: - D_{jk}^\mathrm{(2)} &= 2 T_{ij}^{ab} t_{ik}^{ab} \\ - D_{ab}^\mathrm{(2)} &= 2 T_{ij}^{ac} t_{ij}^{bc} - - By definition of this program, 1-RDM is not response density. + D_{jk} &= T_{ij}^{ab} t_{ik}^{ab} \\ + D_{ab} &= T_{ij}^{ac} t_{ij}^{bc} Parameters ---------- @@ -168,8 +166,6 @@ def get_G_uov( verbose=lib.logger.NOTE, max_memory=2000): r""" Get 3-index transformed MP2 amplitude. - Evaluation of this 3-index amplitude is performed together with non-response MP2 correlation 1-RDM. - .. math:: \Gamma_{ia, P} = T_{ij}^{ab} Y_{jb, P} @@ -224,26 +220,6 @@ def get_G_uov( 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() @@ -270,34 +246,6 @@ def get_W_I(cderi_uov, cderi_uoo, G_uov, verbose=lib.logger.NOTE): 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 - cderi_uaa - W_I - rdm1_corr - Ax0_Core - max_memory - verbose - - Returns - ------- - dict[str, np.ndarray] - """ log = lib.logger.new_logger(verbose=verbose) time0 = lib.logger.process_clock(), lib.logger.perf_counter() @@ -330,32 +278,6 @@ 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() @@ -379,6 +301,9 @@ def get_rdm1_corr_resp( max_cycle=max_cycle, tol=tol)[0] rdm1_corr_resp[sv, so] = rdm1_corr_resp_vo + # symmetrize + rdm1_corr_resp = 0.5 * (rdm1_corr_resp + rdm1_corr_resp.T) + log.timer("get_rdm1_corr_resp", *time0) tensors = {"rdm1_corr_resp": rdm1_corr_resp} return tensors @@ -396,7 +321,7 @@ class RMP2RI(MP2Base): self._Ax0_Core = NotImplemented def make_cderi_uov(self): - r""" Generate cholesky decomposed ERI :math:`Y_{ia, P}` (in memory, occ-vir block). """ + """ Generate cholesky decomposed ERI (in memory, occ-vir block). """ mask = self.get_frozen_mask() nocc_act = (mask & (self.mo_occ != 0)).sum() nvir_act = (mask & (self.mo_occ == 0)).sum() @@ -413,7 +338,7 @@ class RMP2RI(MP2Base): return cderi_uov def make_cderi_uoo(self): - r""" Generate cholesky decomposed ERI math:`Y_{ij, P}` (in memory, occ-occ block). """ + """ Generate cholesky decomposed ERI (in memory, occ-occ block). """ mask = self.get_frozen_mask() idx_occ = (self.mo_occ > 0) & mask mo_occ_act = self.mo_coeff[:, idx_occ] @@ -429,7 +354,7 @@ class RMP2RI(MP2Base): return cderi_uoo def make_cderi_uaa(self): - r""" Generate cholesky decomposed ERI math:`Y_{pq, P}` (all block, s1 symm, in memory/disk). """ + """ Generate cholesky decomposed ERI (all block, s1 symm, in memory/disk). """ log = lib.logger.new_logger(verbose=self.verbose) # dimension and mask @@ -460,7 +385,7 @@ class RMP2RI(MP2Base): return cderi_uaa def driver_eng_mp2(self, **kwargs): - r""" Driver of MP2 energy. """ + """ Driver of MP2 energy. """ mask = self.get_frozen_mask() nocc_act = (mask & (self.mo_occ != 0)).sum() nvir_act = (mask & (self.mo_occ == 0)).sum() @@ -505,7 +430,7 @@ class RMP2RI(MP2Base): @property def Ax0_Core(self): - r""" Fock response of underlying SCF object in MO basis. """ + """ 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 import RHDFT, UHDFT @@ -519,7 +444,6 @@ class RMP2RI(MP2Base): self._Ax0_Core = Ax0_Core def make_rdm1_corr(self): - r""" Generate 1-RDM (non-response) of MP2 correlation :math:`D_{pq}^{(2)}`. """ log = lib.logger.new_logger(verbose=self.verbose) mask = self.get_frozen_mask() mo_energy = self.mo_energy[mask] @@ -538,7 +462,6 @@ class RMP2RI(MP2Base): return tensors["rdm1_corr"] def make_G_uov(self): - r""" Generate 3-index transformed MP2 amplitude :math:`\Gamma_{ia, P}`. """ log = lib.logger.new_logger(verbose=self.verbose) mask = self.get_frozen_mask() mo_energy = self.mo_energy[mask] @@ -557,7 +480,6 @@ class RMP2RI(MP2Base): return tensors["G_uov"] def make_W_I(self): - r""" Generate part I of MO-energy-weighted density matrix. :math:`W_{pq} [\mathrm{I}]`. """ log = lib.logger.new_logger(verbose=self.verbose) # prepare tensors @@ -575,7 +497,6 @@ class RMP2RI(MP2Base): return tensors["W_I"] def make_lag_vo(self): - r""" Generate MP2 contribution to Lagrangian vir-occ block :math:`L_{ai}`. """ # prepare tensors W_I = self.tensors.get("W_I", self.make_W_I()) cderi_uaa = self.tensors.get("cderi_uaa", self.make_cderi_uaa()) @@ -592,7 +513,6 @@ class RMP2RI(MP2Base): return tensors["lag_vo"] def make_rdm1_corr_resp(self): - """ Generate 1-RDM (response) of MP2 correlation :math:`D_{pq}^{(2)}`. """ # prepare input lag_vo = self.tensors.get("lag_vo", self.make_lag_vo()) rdm1_corr = self.tensors["rdm1_corr"] @@ -613,7 +533,6 @@ class RMP2RI(MP2Base): 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()) @@ -627,7 +546,6 @@ class RMP2RI(MP2Base): 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()) -- Gitee From 9e6a8ac3c9b1612032b4972c46b6c2849dba4638 Mon Sep 17 00:00:00 2001 From: ajz34 Date: Fri, 5 May 2023 14:11:59 +0800 Subject: [PATCH 09/12] Revert "feat(dh eng): add dipole (wip)" This reverts commit a555ccaf --- pyscf/dh/__init__.py | 2 +- pyscf/dh/energy/__init__.py | 2 +- pyscf/dh/energy/mp2/__init__.py | 3 +- pyscf/dh/energy/mp2/rmp2.py | 186 +++++++++- pyscf/dh/energy/mp2/rmp2ri.py | 633 -------------------------------- pyscf/dh/energy/mp2/ump2.py | 24 +- pyscf/dh/energy/rhdft.py | 38 +- pyscf/dh/util/df_addon.py | 20 +- pyscf/dh/util/frozen_core.py | 7 +- pyscf/dh/util/helper.py | 10 +- 10 files changed, 190 insertions(+), 735 deletions(-) delete mode 100644 pyscf/dh/energy/mp2/rmp2ri.py diff --git a/pyscf/dh/__init__.py b/pyscf/dh/__init__.py index 0fd3963..22b0475 100644 --- a/pyscf/dh/__init__.py +++ b/pyscf/dh/__init__.py @@ -3,7 +3,7 @@ from . import energy from .energy import ( RHDFT, UHDFT, - RMP2RIPySCF, RMP2ConvPySCF, UMP2ConvPySCF, RMP2Conv, UMP2Conv, UMP2RI, + RMP2RIPySCF, RMP2ConvPySCF, UMP2ConvPySCF, RMP2Conv, RMP2RI, UMP2Conv, UMP2RI, RIEPARI, RIEPAConv, UIEPAConv, UIEPARI, RRingCCDConv ) diff --git a/pyscf/dh/energy/__init__.py b/pyscf/dh/energy/__init__.py index 04b1d25..843baee 100644 --- a/pyscf/dh/energy/__init__.py +++ b/pyscf/dh/energy/__init__.py @@ -6,7 +6,7 @@ from .uhdft import UHDFT from . import mp2, iepa from .mp2 import ( RMP2RIPySCF, RMP2ConvPySCF, UMP2ConvPySCF, - RMP2Conv, UMP2Conv, RMP2RI, UMP2RI) + RMP2Conv, RMP2RI, UMP2Conv, UMP2RI) from .iepa import (RIEPAConv, RIEPARI, UIEPAConv, UIEPARI) from .ring_ccd import RRingCCDConv diff --git a/pyscf/dh/energy/mp2/__init__.py b/pyscf/dh/energy/mp2/__init__.py index cbe73a7..5ed9b16 100644 --- a/pyscf/dh/energy/mp2/__init__.py +++ b/pyscf/dh/energy/mp2/__init__.py @@ -1,4 +1,3 @@ from . import rmp2, ump2 -from .rmp2 import RMP2ConvPySCF, RMP2Conv, RMP2RIPySCF -from .rmp2ri import RMP2RI +from .rmp2 import RMP2ConvPySCF, RMP2Conv, RMP2RI, RMP2RIPySCF from .ump2 import UMP2ConvPySCF, UMP2Conv, UMP2RI diff --git a/pyscf/dh/energy/mp2/rmp2.py b/pyscf/dh/energy/mp2/rmp2.py index 40cd88b..7404ebe 100644 --- a/pyscf/dh/energy/mp2/rmp2.py +++ b/pyscf/dh/energy/mp2/rmp2.py @@ -56,11 +56,9 @@ CONFIG_etb_first = getattr(__config__, "etb_first", False) class MP2Base(EngBase): """ Restricted MP2 base class. """ - def __init__(self, mf, frozen=None, omega=0, with_df=None, c_os=1.0, c_ss=1.0, **kwargs): + def __init__(self, mf, frozen=None, omega=0, with_df=None, **kwargs): super().__init__(mf) self.omega = omega - self.c_os = c_os - self.c_ss = c_ss self.incore_t_oovv_mp2 = CONFIG_incore_t_oovv_mp2 self.frozen = frozen if frozen is not None else 0 self.frac_num = None @@ -78,10 +76,8 @@ class MP2Base(EngBase): class RMP2ConvPySCF(mp.mp2.RMP2, EngBase): """ Restricted MP2 class of doubly hybrid with conventional integral evaluated by PySCF. """ - def __init__(self, mf, c_os=1.0, c_ss=1.0, *args, **kwargs): + def __init__(self, mf, *args, **kwargs): EngBase.__init__(self, mf) - self.c_os = c_os - self.c_ss = c_ss super().__init__(mf, *args, **kwargs) def get_frozen_mask(self): # type: () -> np.ndarray @@ -91,7 +87,7 @@ class RMP2ConvPySCF(mp.mp2.RMP2, EngBase): kernel_output = super().kernel(*args, **kwargs) self.results["eng_corr_MP2_OS"] = self.e_corr_os self.results["eng_corr_MP2_SS"] = self.e_corr_ss - self.results["eng_corr_MP2"] = self.c_os * self.e_corr_os + self.c_ss * self.e_corr_ss + self.results["eng_corr_MP2"] = self.e_corr return kernel_output kernel = driver_eng_mp2 @@ -110,10 +106,8 @@ class RMP2RIPySCF(mp.dfmp2.DFMP2, EngBase): def update_amps(self, t2, eris): raise NotImplementedError - def __init__(self, mf, c_os=1.0, c_ss=1.0, *args, **kwargs): + def __init__(self, mf, *args, **kwargs): EngBase.__init__(self, mf) - self.c_os = c_os - self.c_ss = c_ss super().__init__(mf, *args, **kwargs) def get_frozen_mask(self): # type: () -> np.ndarray @@ -123,7 +117,7 @@ class RMP2RIPySCF(mp.dfmp2.DFMP2, EngBase): kernel_output = super().kernel(*args, **kwargs) self.results["eng_corr_MP2_OS"] = self.e_corr_os self.results["eng_corr_MP2_SS"] = self.e_corr_ss - self.results["eng_corr_MP2"] = self.c_os * self.e_corr_os + self.c_ss * self.e_corr_ss + self.results["eng_corr_MP2"] = self.e_corr return kernel_output kernel = driver_eng_mp2 @@ -134,7 +128,7 @@ class RMP2RIPySCF(mp.dfmp2.DFMP2, EngBase): # region RMP2Conv def kernel_energy_rmp2_conv_full_incore( - mo_energy, mo_coeff, eri_or_mol, mo_occ, c_os, c_ss, + mo_energy, mo_coeff, eri_or_mol, mo_occ, t_oovv=None, frac_num=None, verbose=lib.logger.NOTE, **_kwargs): """ Kernel of restricted MP2 energy by conventional method. @@ -148,10 +142,6 @@ def kernel_energy_rmp2_conv_full_incore( ERI that is recognized by ``pyscf.ao2mo.general``. mo_occ : np.ndarray Molecular orbital occupation numbers. - c_os : float - Coefficient of oppo-spin contribution. - c_ss : float - Coefficient of same-spin contribution. t_oovv : np.ndarray Store space for ``t_oovv`` @@ -213,7 +203,7 @@ def kernel_energy_rmp2_conv_full_incore( # report eng_os = eng_bi1 eng_ss = eng_bi1 - eng_bi2 - eng_mp2 = c_os * eng_os + c_ss * eng_ss + eng_mp2 = eng_os + eng_ss # results results = dict() results["eng_corr_MP2_bi1"] = eng_bi1 @@ -254,7 +244,6 @@ class RMP2Conv(MP2Base): eri_or_mol = eri_or_mol if eri_or_mol is not None else self.mol results = self.kernel_energy_mp2( mo_energy_act, mo_coeff_act, eri_or_mol, mo_occ_act, - self.c_os, self.c_ss, t_oovv=t_oovv, frac_num=frac_num_act, verbose=self.verbose, **kwargs) self.e_corr = results["eng_corr_MP2"] # pad omega @@ -268,10 +257,159 @@ class RMP2Conv(MP2Base): # endregion +# region RMP2RI + +def kernel_energy_rmp2_ri_incore( + mo_energy, cderi_uov, + t_oovv=None, frac_num=None, verbose=lib.logger.NOTE, max_memory=2000, cderi_uov_2=None): + """ Kernel of MP2 energy by RI integral. + + Parameters + ---------- + mo_energy : np.ndarray + Molecular orbital energy levels. + cderi_uov : np.ndarray + Cholesky decomposed 3c2e ERI in MO basis (occ-vir part). + + t_oovv : np.ndarray + Store space for ``t_oovv`` + frac_num : np.ndarray + Fractional occupation number list. + verbose : int + Verbose level for PySCF. + max_memory : float + Allocatable memory in MB. + cderi_uov_2 : np.ndarray + Another part of 3c2e ERI in MO basis (occ-vir part). This is mostly used in magnetic computations. + + Notes + ----- + + For RI approximation, ERI integral is set to be + + .. math:: + g_{ij}^{ab} = (ia|jb) = Y_{ia, P} Y_{jb, P} + + For energy of fractional occupation system, computation method is chosen + by eq (7) from Su2016 (10.1021/acs.jctc.6b00197). + """ + log = lib.logger.new_logger(verbose=verbose) + log.info("[INFO] Start unrestricted RI-MP2") + + naux, nocc, nvir = cderi_uov.shape + if frac_num is not None: + frac_occ, frac_vir = frac_num[:nocc], frac_num[nocc:] + else: + frac_occ = frac_vir = None + eo = mo_energy[:nocc] + ev = mo_energy[nocc:] + + # loops + log.info("[INFO] Start RI-MP2 loop") + nbatch = util.calc_batch_size(4 * nocc * nvir ** 2, max_memory, dtype=cderi_uov.dtype) + eng_bi1 = eng_bi2 = 0 + for sI in util.gen_batch(0, nocc, nbatch): + log.info(f"[INFO] MP2 loop i: {sI}") + if cderi_uov_2 is None: + g_Iajb = lib.einsum("PIa, Pjb -> Iajb", cderi_uov[:, sI], cderi_uov) + else: + g_Iajb = 0.5 * lib.einsum("PIa, Pjb -> Iajb", cderi_uov[:, sI], cderi_uov_2) + g_Iajb += 0.5 * lib.einsum("PIa, Pjb -> Iajb", cderi_uov_2[:, sI], cderi_uov) + D_Ijab = lib.direct_sum("i + j - a - b -> ijab", eo[sI], eo, ev, ev) + t_Ijab = lib.einsum("Iajb, Ijab -> Ijab", g_Iajb, 1 / D_Ijab) + if t_oovv is not None: + t_oovv[sI] = t_Ijab + if frac_num is not None: + n_Ijab = lib.einsum("i, j, a, b -> ijab", frac_occ[sI], frac_occ, 1 - frac_vir, 1 - frac_vir) + eng_bi1 += lib.einsum("Ijab, Ijab, Iajb ->", n_Ijab, t_Ijab.conj(), g_Iajb) + eng_bi2 += lib.einsum("Ijab, Ijab, Ibja ->", n_Ijab, t_Ijab.conj(), g_Iajb) + else: + eng_bi1 += lib.einsum("Ijab, Iajb ->", t_Ijab.conj(), g_Iajb) + eng_bi2 += lib.einsum("Ijab, Ibja ->", t_Ijab.conj(), g_Iajb) + eng_bi1 = util.check_real(eng_bi1) + eng_bi2 = util.check_real(eng_bi2) + log.info("[INFO] MP2 energy computation finished.") + # report + eng_os = eng_bi1 + eng_ss = eng_bi1 - eng_bi2 + eng_mp2 = eng_os + eng_ss + # results + 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}") + return results + + +class RMP2RI(MP2Base): + """ Restricted MP2 class of doubly hybrid with RI integral. """ + + def __init__(self, mf, frozen=None, omega=0, with_df=None, **kwargs): + super().__init__(mf, frozen=frozen, omega=omega, with_df=with_df, **kwargs) + self.with_df_2 = None + + def driver_eng_mp2(self, **kwargs): + 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_energy_act = self.mo_energy[mask] + frac_num_act = self.frac_num + if frac_num_act is not None: + frac_num_act = frac_num_act[mask] + # prepare t_oovv + max_memory = self.max_memory - lib.current_memory()[0] + t_oovv = util.allocate_array( + self.incore_t_oovv_mp2, (nocc_act, nocc_act, nvir_act, nvir_act), max_memory, + h5file=self._tmpfile, name="t_oovv_mp2", zero_init=False, chunk=(1, 1, nvir_act, nvir_act), + dtype=mo_coeff_act.dtype) + # generate cderi_uov + omega = self.omega + with_df = util.get_with_df_omega(self.with_df, omega) + max_memory = self.max_memory - lib.current_memory()[0] + cderi_uov = self.tensors.get("cderi_uov", None) + if cderi_uov is None: + cderi_uov = util.get_cderi_mo( + with_df, mo_coeff_act, None, (0, nocc_act, nocc_act, nact), max_memory) + self.tensors["cderi_uov"] = cderi_uov + # cderi_uov_2 is rarely called, so do not try to build omega for this special case + cderi_uov_2 = None + max_memory = self.max_memory - lib.current_memory()[0] + if self.with_df_2 is not None: + cderi_uov_2 = util.get_cderi_mo( + self.with_df_2, mo_coeff_act, None, (0, nocc_act, nocc_act, nact), max_memory) + # kernel + max_memory = self.max_memory - lib.current_memory()[0] + results = self.kernel_energy_mp2( + mo_energy_act, cderi_uov, + t_oovv=t_oovv, + frac_num=frac_num_act, + verbose=self.verbose, + max_memory=max_memory, + cderi_uov_2=cderi_uov_2, + **kwargs) + + self.e_corr = results["eng_corr_MP2"] + # pad omega + results = {pad_omega(key, self.omega): val for (key, val) in results.items()} + self.results.update(results) + return results + + kernel_energy_mp2 = staticmethod(kernel_energy_rmp2_ri_incore) + kernel = driver_eng_mp2 + +# endregion + + if __name__ == '__main__': def main_1(): - # test RMP2ConvPySCF 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 = scf.RHF(mol).run() @@ -280,7 +418,6 @@ if __name__ == '__main__': print(mf_mp.results) def main_2(): - # test RMP2Conv 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 = scf.RHF(mol).run() @@ -288,5 +425,14 @@ if __name__ == '__main__': print(mf_mp.e_tot) print(mf_mp.results) + def main_3(): + 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 = scf.RHF(mol).run() + mf_mp = RMP2RI(mf_scf, frozen=[1, 2]).run() + print(mf_mp.e_tot) + print(mf_mp.results) + main_1() main_2() + main_3() diff --git a/pyscf/dh/energy/mp2/rmp2ri.py b/pyscf/dh/energy/mp2/rmp2ri.py deleted file mode 100644 index ed7d168..0000000 --- a/pyscf/dh/energy/mp2/rmp2ri.py +++ /dev/null @@ -1,633 +0,0 @@ -import numpy as np -from pyscf import lib, scf, __config__ -from pyscf.dh import util -from pyscf.dh.energy.mp2.rmp2 import MP2Base -from pyscf.dh.util import pad_omega -from pyscf.scf import cphf - -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 kernel_energy_rmp2_ri_incore( - mo_energy, cderi_uov, - c_os, c_ss, - t_oovv=None, frac_num=None, verbose=lib.logger.NOTE, max_memory=2000, cderi_uov_2=None): - """ Kernel of MP2 energy by RI integral. - - Parameters - ---------- - mo_energy : np.ndarray - Molecular orbital energy levels. - cderi_uov : np.ndarray - Cholesky decomposed 3c2e ERI in MO basis (occ-vir part). - c_os : float - Coefficient of oppo-spin contribution. - c_ss : float - Coefficient of same-spin contribution. - - t_oovv : np.ndarray - Store space for ``t_oovv`` - frac_num : np.ndarray - Fractional occupation number list. - verbose : int - Verbose level for PySCF. - max_memory : float - Allocatable memory in MB. - cderi_uov_2 : np.ndarray - Another part of 3c2e ERI in MO basis (occ-vir part). This is mostly used in magnetic computations. - - Notes - ----- - - For RI approximation, ERI integral is set to be - - .. math:: - g_{ij}^{ab} = (ia|jb) = Y_{ia, P} Y_{jb, P} - - For energy of fractional occupation system, computation method is chosen - by eq (7) from Su2016 (10.1021/acs.jctc.6b00197). - """ - log = lib.logger.new_logger(verbose=verbose) - time0 = lib.logger.process_clock(), lib.logger.perf_counter() - log.info("[INFO] Start unrestricted RI-MP2") - - naux, nocc, nvir = cderi_uov.shape - if frac_num is not None: - frac_occ, frac_vir = frac_num[:nocc], frac_num[nocc:] - else: - frac_occ = frac_vir = None - eo = mo_energy[:nocc] - ev = mo_energy[nocc:] - - # loops - log.info("[INFO] Start RI-MP2 loop") - mem_avail = max_memory - lib.current_memory()[0] - nbatch = util.calc_batch_size(4 * nocc * nvir ** 2, mem_avail, dtype=cderi_uov.dtype) - eng_bi1 = eng_bi2 = 0 - for sI in util.gen_batch(0, nocc, nbatch): - log.info(f"[INFO] MP2 loop i: {sI}") - if cderi_uov_2 is None: - g_Iajb = lib.einsum("PIa, Pjb -> Iajb", cderi_uov[:, sI], cderi_uov) - else: - g_Iajb = 0.5 * lib.einsum("PIa, Pjb -> Iajb", cderi_uov[:, sI], cderi_uov_2) - g_Iajb += 0.5 * lib.einsum("PIa, Pjb -> Iajb", cderi_uov_2[:, sI], cderi_uov) - D_Ijab = lib.direct_sum("i + j - a - b -> ijab", eo[sI], eo, ev, ev) - t_Ijab = lib.einsum("Iajb, Ijab -> Ijab", g_Iajb, 1 / D_Ijab) - if t_oovv is not None: - t_oovv[sI] = t_Ijab - if frac_num is not None: - n_Ijab = lib.einsum("i, j, a, b -> ijab", frac_occ[sI], frac_occ, 1 - frac_vir, 1 - frac_vir) - eng_bi1 += lib.einsum("Ijab, Ijab, Iajb ->", n_Ijab, t_Ijab.conj(), g_Iajb) - eng_bi2 += lib.einsum("Ijab, Ijab, Ibja ->", n_Ijab, t_Ijab.conj(), g_Iajb) - else: - eng_bi1 += lib.einsum("Ijab, Iajb ->", t_Ijab.conj(), g_Iajb) - eng_bi2 += lib.einsum("Ijab, Ibja ->", t_Ijab.conj(), g_Iajb) - eng_bi1 = util.check_real(eng_bi1) - eng_bi2 = util.check_real(eng_bi2) - log.info("[INFO] MP2 energy computation finished.") - # report - eng_os = eng_bi1 - eng_ss = eng_bi1 - eng_bi2 - eng_mp2 = c_os * eng_os + c_ss * eng_ss - # results - 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}") - log.timer("kernel_energy_rmp2_ri_incore", *time0) - return results - - -def get_rdm1_corr( - mo_energy, cderi_uov, t_oovv, - c_os, c_ss, - verbose=lib.logger.NOTE, max_memory=2000): - r""" MP2 correlation part of density matrix. - - .. math:: - D_{jk} &= T_{ij}^{ab} t_{ik}^{ab} \\ - D_{ab} &= T_{ij}^{ac} t_{ij}^{bc} - - Parameters - ---------- - mo_energy : np.ndarray - cderi_uov : np.ndarray - t_oovv : np.ndarray - c_os : float - c_ss : float - verbose : int - max_memory : int or float - - 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 len(mo_energy) == nmo - assert t_oovv.shape == (nocc, nocc, nvir, nvir) - - # prepare essential matrices - so, sv = slice(0, nocc), slice(nocc, nmo) - - # allocate results - rdm1_corr = np.zeros((nmo, nmo)) - - def load_slice(slc): - return slc, t_oovv[slc] - - mem_avail = max_memory - lib.current_memory()[0] - nbatch = util.calc_batch_size(2 * nocc * nvir**2, mem_avail) - for sI, t_Oovv in lib.map_with_prefetch(load_slice, util.gen_batch(0, nocc, nbatch)): - log.debug(f"[INFO] Generation of MP2 rdm1 oo vv and G_uov, {sI} of total {nocc} occ orbitals") - T_Oovv = util.restricted_biorthogonalize(t_Oovv, 1, c_os, c_ss) - rdm1_corr[sv, sv] += 2 * lib.einsum("ijac, ijbc -> ab", T_Oovv, t_Oovv) - rdm1_corr[so, so] -= 2 * lib.einsum("ijab, ikab -> jk", T_Oovv, t_Oovv) - - log.timer("get_rdm1_corr", *time0) - tensors = {"rdm1_corr": rdm1_corr} - return tensors - - -def get_G_uov( - mo_energy, cderi_uov, t_oovv, - c_os, c_ss, - verbose=lib.logger.NOTE, max_memory=2000): - r""" Get 3-index transformed MP2 amplitude. - - .. math:: - \Gamma_{ia, P} = T_{ij}^{ab} Y_{jb, P} - - Parameters - ---------- - mo_energy : np.ndarray - cderi_uov : np.ndarray - t_oovv : np.ndarray - c_os : float - c_ss : float - verbose : int - max_memory : int or float - - 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 len(mo_energy) == nmo - assert t_oovv.shape == (nocc, nocc, nvir, nvir) - - # prepare essential matrices - so, sv = slice(0, nocc), slice(nocc, nmo) - - # allocate results - rdm1_corr = np.zeros((nmo, nmo)) - G_uov = np.zeros((naux, nocc, nvir)) - - def load_slice(slc): - return slc, t_oovv[slc] - - mem_avail = max_memory - lib.current_memory()[0] - nbatch = util.calc_batch_size(2 * nocc * nvir**2 + naux * nvir, mem_avail) - for sI, t_Oovv in lib.map_with_prefetch(load_slice, util.gen_batch(0, nocc, nbatch)): - log.debug(f"[INFO] Generation of MP2 rdm1 oo vv and G_uov, {sI} of total {nocc} occ orbitals") - T_Oovv = util.restricted_biorthogonalize(t_Oovv, 1, c_os, c_ss) - rdm1_corr[sv, sv] += 2 * lib.einsum("ijac, ijbc -> ab", T_Oovv, t_Oovv) - rdm1_corr[so, so] -= 2 * lib.einsum("ijab, ikab -> jk", T_Oovv, t_Oovv) - G_uov[:, sI] += lib.einsum("ijab, Pjb -> Pia", T_Oovv, cderi_uov) - - log.timer("get_G_uov", *time0) - tensors = { - "rdm1_corr": rdm1_corr, - "G_uov": G_uov, - } - return tensors - - -def get_W_I(cderi_uov, cderi_uoo, G_uov, verbose=lib.logger.NOTE): - - 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): - - 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): - - 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 - - # symmetrize - rdm1_corr_resp = 0.5 * (rdm1_corr_resp + rdm1_corr_resp.T) - - log.timer("get_rdm1_corr_resp", *time0) - tensors = {"rdm1_corr_resp": rdm1_corr_resp} - return tensors - - -class RMP2RI(MP2Base): - """ Restricted MP2 class of doubly hybrid with RI integral. """ - - def __init__(self, mf, **kwargs): - super().__init__(mf, **kwargs) - self.with_df_2 = None - self.incore_cderi_uaa = CONFIG_incore_cderi_uaa_mp2 - self.max_cycle_cpks = CONFIG_max_cycle_cpks - self.tol_cpks = CONFIG_tol_cpks - self._Ax0_Core = NotImplemented - - def make_cderi_uov(self): - """ Generate cholesky decomposed ERI (in memory, occ-vir block). """ - 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] - omega = self.omega - with_df = util.get_with_df_omega(self.with_df, omega) - cderi_uov = util.get_cderi_mo( - with_df, mo_coeff_act, None, (0, nocc_act, nocc_act, nact), - max_memory=self.max_memory, verbose=self.verbose) - - tensors = {"cderi_uov": cderi_uov} - self.tensors.update(tensors) - return cderi_uov - - def make_cderi_uoo(self): - """ Generate cholesky decomposed ERI (in memory, occ-occ block). """ - mask = self.get_frozen_mask() - idx_occ = (self.mo_occ > 0) & mask - mo_occ_act = self.mo_coeff[:, idx_occ] - nocc_act = idx_occ.sum() - omega = self.omega - with_df = util.get_with_df_omega(self.with_df, omega) - cderi_uoo = util.get_cderi_mo( - with_df, mo_occ_act, None, (0, nocc_act, 0, nocc_act), - max_memory=self.max_memory, verbose=self.verbose) - - tensors = {"cderi_uoo": cderi_uoo} - self.tensors.update(tensors) - return cderi_uoo - - def make_cderi_uaa(self): - """ Generate cholesky decomposed ERI (all block, s1 symm, in memory/disk). """ - 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 driver_eng_mp2(self, **kwargs): - """ Driver of MP2 energy. """ - 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_energy_act = self.mo_energy[mask] - frac_num_act = self.frac_num - if frac_num_act is not None: - frac_num_act = frac_num_act[mask] - # prepare t_oovv - max_memory = self.max_memory - lib.current_memory()[0] - t_oovv = util.allocate_array( - self.incore_t_oovv_mp2, (nocc_act, nocc_act, nvir_act, nvir_act), max_memory, - h5file=self._tmpfile, name="t_oovv_mp2", zero_init=False, chunk=(1, 1, nvir_act, nvir_act), - dtype=mo_coeff_act.dtype) - if t_oovv is not None: - self.tensors["t_oovv"] = t_oovv - # generate cderi_uov - cderi_uov = self.tensors.get("cderi_uov", self.make_cderi_uov()) - # cderi_uov_2 is rarely called, so do not try to build omega for this special case - cderi_uov_2 = None - if self.with_df_2 is not None: - cderi_uov_2 = util.get_cderi_mo( - self.with_df_2, mo_coeff_act, None, (0, nocc_act, nocc_act, nact), - max_memory=self.max_memory, verbose=self.verbose) - # kernel - results = self.kernel_energy_mp2( - mo_energy_act, cderi_uov, - self.c_os, self.c_ss, - t_oovv=t_oovv, - frac_num=frac_num_act, - verbose=self.verbose, - max_memory=self.max_memory, - cderi_uov_2=cderi_uov_2, - **kwargs) - - self.e_corr = results["eng_corr_MP2"] - # pad omega - results = {pad_omega(key, self.omega): val for (key, val) in results.items()} - self.results.update(results) - return results - - @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 import RHDFT, UHDFT - HDFT = RHDFT if restricted else UHDFT - mf_scf = HDFT(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 - - def make_rdm1_corr(self): - log = lib.logger.new_logger(verbose=self.verbose) - mask = self.get_frozen_mask() - mo_energy = self.mo_energy[mask] - cderi_uov = self.tensors["cderi_uov"] - t_oovv = self.tensors.get("t_oovv", None) - if t_oovv is None: - log.warn("t_oovv for MP2 has not been stored. " - "To perform G_uov, t_oovv and MP2 energy is to be re-evaluated.") - self.incore_t_oovv_mp2 = True - self.driver_eng_mp2() - t_oovv = self.tensors["t_oovv"] - tensors = get_rdm1_corr( - mo_energy, cderi_uov, t_oovv, self.c_ss, self.c_os, - verbose=self.verbose, max_memory=self.max_memory) - self.tensors.update(tensors) - return tensors["rdm1_corr"] - - def make_G_uov(self): - log = lib.logger.new_logger(verbose=self.verbose) - mask = self.get_frozen_mask() - mo_energy = self.mo_energy[mask] - cderi_uov = self.tensors["cderi_uov"] - t_oovv = self.tensors.get("t_oovv", None) - if t_oovv is None: - log.warn("t_oovv for MP2 has not been stored. " - "To perform G_uov, t_oovv and MP2 energy is to be re-evaluated.") - self.incore_t_oovv_mp2 = True - self.driver_eng_mp2() - t_oovv = self.tensors["t_oovv"] - tensors = get_G_uov( - mo_energy, cderi_uov, t_oovv, self.c_ss, self.c_os, - verbose=self.verbose, max_memory=self.max_memory) - self.tensors.update(tensors) - return tensors["G_uov"] - - def make_W_I(self): - log = lib.logger.new_logger(verbose=self.verbose) - - # prepare tensors - if any([key not in self.tensors for key in ["rdm1_corr", "G_uov", "t_oovv"]]): - log.warn("Some tensors (rdm1_corr, G_uov, t_oovv) has not been generated. " - "Perform make_mp2_integrals first.") - self.make_G_uov() - cderi_uov = self.tensors["cderi_uov"] - cderi_uoo = self.tensors.get("cderi_uoo", self.make_cderi_uoo()) - G_uov = self.tensors["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): - # 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): - # 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): - # 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): - # 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 - - kernel_energy_mp2 = staticmethod(kernel_energy_rmp2_ri_incore) - kernel = driver_eng_mp2 - - -if __name__ == '__main__': - - def main_1(): - # test RMP2RI - 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 = scf.RHF(mol).run() - mf_mp = RMP2RI(mf_scf, frozen=[1, 2]).run() - print(mf_mp.e_tot) - print(mf_mp.results) - - def main_2(): - # test response rdm1 generation - from pyscf import gto, scf - np.set_printoptions(4, 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_scf = scf.RHF(mol).run() - mf_mp = RMP2RI(mf_scf, frozen=[2]).run(incore_t_oovv_mp2=True) - mf_mp.make_rdm1() - for key, val in mf_mp.tensors.items(): - print(key, val.shape) - print(mf_mp.tensors["rdm1"]) - - def main_3(): - # 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 = RMP2RI(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_3() diff --git a/pyscf/dh/energy/mp2/ump2.py b/pyscf/dh/energy/mp2/ump2.py index af1820a..dd0d739 100644 --- a/pyscf/dh/energy/mp2/ump2.py +++ b/pyscf/dh/energy/mp2/ump2.py @@ -2,8 +2,7 @@ r""" Unrestricted MP2. """ from pyscf.dh import util -from pyscf.dh.energy.mp2.rmp2 import RMP2ConvPySCF, RMP2Conv -from pyscf.dh.energy.mp2.rmp2ri import RMP2RI +from pyscf.dh.energy.mp2.rmp2 import RMP2ConvPySCF, RMP2Conv, RMP2RI from pyscf import ao2mo, lib, mp, __config__, df import numpy as np @@ -40,7 +39,6 @@ class UMP2ConvPySCF(mp.ump2.UMP2, RMP2ConvPySCF): def kernel_energy_ump2_conv_full_incore( mo_energy, mo_coeff, eri_or_mol, mo_occ, - c_os, c_ss, t_oovv=None, frac_num=None, verbose=lib.logger.NOTE): """ Kernel of unrestricted MP2 energy by conventional method. @@ -52,10 +50,6 @@ def kernel_energy_ump2_conv_full_incore( Molecular coefficients. eri_or_mol : np.ndarray or gto.Mole ERI that is recognized by ``pyscf.ao2mo.general``. - c_os : float - Coefficient of oppo-spin contribution. - c_ss : float - Coefficient of same-spin contribution. t_oovv : list[np.ndarray] Store space for ``t_oovv`` @@ -126,7 +120,7 @@ def kernel_energy_ump2_conv_full_incore( eng_spin[0] *= 0.25 eng_spin[2] *= 0.25 eng_spin = util.check_real(eng_spin) - eng_mp2 = c_os * eng_spin[1] + c_ss * (eng_spin[0] + eng_spin[2]) + eng_mp2 = eng_spin[1] + (eng_spin[0] + eng_spin[2]) # finalize results results = dict() results["eng_corr_MP2_aa"] = eng_spin[0] @@ -167,7 +161,7 @@ class UMP2Conv(RMP2Conv): for s0, s1, ss, ssn in ((0, 0, 0, "aa"), (0, 1, 1, "ab"), (1, 1, 2, "bb")): t_oovv[ss] = util.allocate_array( incore_t_oovv, shape=(nocc_act[s0], nocc_act[s1], nvir_act[s0], nvir_act[s1]), - mem_avail=max_memory, + max_memory=max_memory, h5file=self._tmpfile, name=f"t_oovv_{ssn}", dtype=mo_coeff_act[0].dtype) @@ -178,7 +172,6 @@ class UMP2Conv(RMP2Conv): eri_or_mol = eri_or_mol if eri_or_mol is not None else self.mol results = self.kernel_energy_mp2( mo_energy_act, mo_coeff_act, eri_or_mol, mo_occ_act, - self.c_os, self.c_ss, t_oovv=t_oovv, frac_num=frac_num_act, verbose=self.verbose, **kwargs) self.e_corr = results["eng_corr_MP2"] # pad omega @@ -195,7 +188,7 @@ class UMP2Conv(RMP2Conv): # region UMP2RI def kernel_energy_ump2_ri_incore( - mo_energy, cderi_uov, c_os, c_ss, + mo_energy, cderi_uov, t_oovv=None, frac_num=None, verbose=lib.logger.NOTE, max_memory=2000, cderi_uov_2=None): """ Kernel of unrestricted MP2 energy by RI integral. @@ -210,10 +203,6 @@ def kernel_energy_ump2_ri_incore( Molecular orbital energy levels. cderi_uov : list[np.ndarray] Cholesky decomposed 3c2e ERI in MO basis (occ-vir part). Spin in (aa, bb). - c_os : float - Coefficient of oppo-spin contribution. - c_ss : float - Coefficient of same-spin contribution. t_oovv : list[np.ndarray] Store space for ``t_oovv`` @@ -280,7 +269,7 @@ def kernel_energy_ump2_ri_incore( eng_spin[0] *= 0.25 eng_spin[2] *= 0.25 eng_spin = util.check_real(eng_spin) - eng_mp2 = c_os * eng_spin[1] + c_ss * (eng_spin[0] + eng_spin[2]) + eng_mp2 = eng_spin[1] + (eng_spin[0] + eng_spin[2]) # finalize results results = dict() results["eng_corr_MP2_aa"] = eng_spin[0] @@ -321,7 +310,7 @@ class UMP2RI(RMP2RI): for s0, s1, ss, ssn in ((0, 0, 0, "aa"), (0, 1, 1, "ab"), (1, 1, 2, "bb")): t_oovv[ss] = util.allocate_array( incore_t_oovv, shape=(nocc_act[s0], nocc_act[s1], nvir_act[s0], nvir_act[s1]), - mem_avail=max_memory, + max_memory=max_memory, h5file=self._tmpfile, name=f"t_oovv_{ssn}", dtype=mo_coeff_act[0].dtype) @@ -349,7 +338,6 @@ class UMP2RI(RMP2RI): max_memory = self.max_memory - lib.current_memory()[0] results = self.kernel_energy_mp2( mo_energy_act, cderi_uov, - self.c_os, self.c_ss, t_oovv=t_oovv, frac_num=frac_num_act, verbose=self.verbose, diff --git a/pyscf/dh/energy/rhdft.py b/pyscf/dh/energy/rhdft.py index 3cdce2c..d79cf01 100644 --- a/pyscf/dh/energy/rhdft.py +++ b/pyscf/dh/energy/rhdft.py @@ -313,7 +313,7 @@ def numint_customized(xc, _mol=None): 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. + :math:`\sum_{rs} A_{pq, rs} X_{rs}^\mathbb{A}`. Parameters ---------- @@ -430,15 +430,9 @@ class RHDFT(EngBase): It's better to initialize this object first, before actually running SCF iterations. """ - def __init__(self, mf, xc=None): + def __init__(self, mf, xc): super().__init__(mf) - if xc is None: - if hasattr(mf, "xc"): - self.xc = xc - else: - self.xc = "HF" # not KS object - else: - self.xc = xc # type: XCList + self.xc = xc # type: XCList if isinstance(self.xc, str): xc_scf = XCList(self.xc, code_scf=True) xc_eng = XCList(self.xc, code_scf=False) @@ -491,33 +485,9 @@ class RHDFT(EngBase): """ 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. + :math:`\sum_{rs} A_{pq, rs} X_{rs}^\mathbb{A}`. Parameters ---------- diff --git a/pyscf/dh/util/df_addon.py b/pyscf/dh/util/df_addon.py index c7e6ad2..b05f610 100644 --- a/pyscf/dh/util/df_addon.py +++ b/pyscf/dh/util/df_addon.py @@ -1,13 +1,11 @@ import numpy as np import copy -from pyscf import df, lib, __config__ +from pyscf import df, lib 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, verbose=lib.logger.NOTE): +def get_cderi_mo(with_df, mo_coeff, Y_mo=None, pqslice=None, max_memory=2000): """ Get Cholesky decomposed ERI in MO basis. Parameters @@ -21,18 +19,13 @@ def get_cderi_mo( pqslice : tuple[int] Slice of orbitals to be transformed. max_memory : float - Maximum memory of molecular object. - verbose : int - Print level verbosity. + Maximum memory available. 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: @@ -45,16 +38,13 @@ def get_cderi_mo( p0, p1 = 0, 0 preflop = 0 if not isinstance(Y_mo, np.ndarray) else Y_mo.size - mem_avail = max_memory - lib.current_memory()[0] - nbatch = calc_batch_size(2*nump*numq, mem_avail, preflop) - log.debug(f"[DEBUG] Number of available batch size in get_cderi_mo: {nbatch}") + nbatch = calc_batch_size(2*nump*numq, max_memory, preflop) if hasattr(with_df._cderi, "shape"): # array alike aosym = "s2" if len(with_df._cderi.shape) == 2 else "s1" else: # assert pyscf when memory is low aosym = "s2" for Y_ao in with_df.loop(nbatch): p1 = p0 + Y_ao.shape[0] - log.debug(f"[DEBUG] cderi transformation of auxiliary basis: {p0, p1}") if Y_ao.dtype == np.double and mo_coeff.dtype == np.double: Y_mo[p0:p1] = _ao2mo.nr_e2(Y_ao, mo_coeff, pqslice, aosym=aosym, mosym="s1").reshape(p1 - p0, nump, numq) else: @@ -65,8 +55,6 @@ def get_cderi_mo( Y_ao.reshape(naux, nao, nao), mo_coeff[:, pqslice[0]:pqslice[1]].conj(), mo_coeff[:, pqslice[2]:pqslice[3]]) p0 = p1 - - log.timer("get_cderi_mo", *time0) return Y_mo diff --git a/pyscf/dh/util/frozen_core.py b/pyscf/dh/util/frozen_core.py index 420cedb..5621431 100644 --- a/pyscf/dh/util/frozen_core.py +++ b/pyscf/dh/util/frozen_core.py @@ -875,13 +875,10 @@ class FrozenCore: mask = np.ones((nset, nmo), dtype=bool) if isinstance(self.rule, list or np.ndarray): - arr = self.rule - # 0. frozen set is empty - if len(arr) == 0: - mask[:] = True # 1. case of frozen/active list defined, highest priority # first perform mask; if self.act, then reverse mask in order to make active orbitals to be True - elif not hasattr(arr[0], "__iter__"): + arr = self.rule + if not hasattr(arr[0], "__iter__"): # one list of frozen orbitals: mask[:, arr] = False else: diff --git a/pyscf/dh/util/helper.py b/pyscf/dh/util/helper.py index c1afd2a..b47938d 100644 --- a/pyscf/dh/util/helper.py +++ b/pyscf/dh/util/helper.py @@ -86,7 +86,7 @@ def parse_incore_flag(flag, unit_flop, mem_avail, pre_flop=0, dtype=float): Memory available in MB. pre_flop : int Number of data preserved in memory. Unit in number. - dtype : np.dtype + dtype : type Type of data. Such as np.float64, complex, etc. Returns @@ -217,7 +217,7 @@ def pad_omega(s, omega): return f"{s}_omega({omega:.6f})" -def allocate_array(incore, shape, mem_avail, +def allocate_array(incore, shape, max_memory, h5file=None, name=None, dtype=float, zero_init=True, chunk=None, **kwargs): """ Allocate an array with given memory estimation and incore stragety. @@ -227,8 +227,8 @@ def allocate_array(incore, shape, mem_avail, Incore flag. Also see :class:`parse_incore_flag`. shape : tuple Shape of allocated array. - mem_avail : int or float - Memory available for allocation. + max_memory : int or float + Maximum memory in MB. h5file : h5py.File HDF5 file instance. Array may save into this file if disk is required. name : str @@ -244,7 +244,7 @@ def allocate_array(incore, shape, mem_avail, ------- np.array or h5py.Dataset """ - incore = parse_incore_flag(incore, int(np.prod(shape)), mem_avail, dtype=dtype) + incore = parse_incore_flag(incore, int(np.prod(shape)), max_memory, dtype=dtype) if incore is None: return None elif incore is True: -- Gitee From 578e762e7d2d70bdc7672ff33a192c23534503fe Mon Sep 17 00:00:00 2001 From: ajz34 Date: Fri, 5 May 2023 14:12:04 +0800 Subject: [PATCH 10/12] Revert "feat(dh eng): add Ax0_Core_resp" This reverts commit 66ee971a429884ebebda80fe338078fea64c42b4. --- pyscf/dh/energy/rhdft.py | 93 ++-------------------------------------- 1 file changed, 3 insertions(+), 90 deletions(-) diff --git a/pyscf/dh/energy/rhdft.py b/pyscf/dh/energy/rhdft.py index d79cf01..4a971e4 100644 --- a/pyscf/dh/energy/rhdft.py +++ b/pyscf/dh/energy/rhdft.py @@ -2,10 +2,8 @@ from pyscf.dh import util from pyscf.dh.util import XCList, XCType, XCInfo from pyscf.dh.energy import EngBase from pyscf import dft, lib, scf, df, __config__ -from pyscf.scf import _response_functions # this import is not unnecessary import numpy as np from types import MethodType -from functools import cache, cached_property CONFIG_ssr_x_fr = getattr(__config__, "ssr_x_fr", "LDA_X") CONFIG_ssr_x_sr = getattr(__config__, "ssr_x_sr", "LDA_X_ERF") @@ -110,7 +108,7 @@ def get_energy_purexc(xc_lists, rho, weights, restricted, numint=None): """ Evaluate energy contributions of pure (DFT) exchange-correlation effects. Note that this kernel does not count HF, LR_HF and advanced correlation into account. - To evaluate exact exchange (HF or LR_HF), use ``get_energy_restricted_exactx``. + To evaluate exact exchange (HF or LR_HF), use ``kernel_energy_restricted_exactx``. Parameters ---------- @@ -133,7 +131,7 @@ def get_energy_purexc(xc_lists, rho, weights, restricted, numint=None): See Also -------- - get_energy_restricted_exactx + kernel_energy_restricted_exactx """ if not isinstance(xc_lists, list): xc_lists = [xc_lists] @@ -311,40 +309,6 @@ def numint_customized(xc, _mol=None): return ni_custom -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}`. - - 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 - - def custom_mf(mf, xc, auxbasis_or_with_df=None): """ Customize options of PySCF's mf object. @@ -480,39 +444,6 @@ class RHDFT(EngBase): def kernel(self, *args, **kwargs): return self.scf.kernel(*args, **kwargs) - @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_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}`. - - 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) - get_energy_exactx = staticmethod(get_energy_restricted_exactx) get_energy_noxc = staticmethod(get_energy_restricted_noxc) get_energy_vv10 = staticmethod(get_energy_vv10) @@ -529,22 +460,4 @@ if __name__ == '__main__': res = mf.make_energy_purexc([", LYP", "B88, ", "HF", "LR_HF(0.5)", "SSR(GGA_X_B88, 0.5), "]) print(res) - def main_2(): - from pyscf import gto, scf - mol = gto.Mole(atom="O; H 1 0.94; H 1 0.94 2 104.5").build() - mf_s = scf.RHF(mol) - mf = RHDFT(mf_s, xc="B3LYPg").run() - np.random.seed(0) - nocc, nvir = mf.nocc, mf.nvir - dmX = np.random.randn(3, nocc, nvir) - mo_occ = mf.mo_occ - idx_occ = mo_occ > 0 - idx_vir = ~idx_occ - log = lib.logger.new_logger(verbose=5) - # first epoch should be a lot slower, since mf.scf.gen_response will cache xc kernel first - for epoch in range(3): - t0, t1 = lib.logger.process_clock(), lib.logger.perf_counter() - ax = mf.Ax0_Core_resp(idx_occ, idx_vir, idx_occ, idx_vir)(dmX) - lib.logger.timer(log, f"Epoch {epoch} of Ax0_Core_resp", t0, t1) - - main_2() + main_1() -- Gitee From d22b838384a1892e5b7e775f96d15f667d2bffdf Mon Sep 17 00:00:00 2001 From: ajz34 Date: Fri, 5 May 2023 14:34:46 +0800 Subject: [PATCH 11/12] fix(dh eng): fix energy evaluation for multiple omega --- pyscf/dh/energy/dh.py | 9 ++------- 1 file changed, 2 insertions(+), 7 deletions(-) diff --git a/pyscf/dh/energy/dh.py b/pyscf/dh/energy/dh.py index a009a7a..6025e65 100644 --- a/pyscf/dh/energy/dh.py +++ b/pyscf/dh/energy/dh.py @@ -321,12 +321,8 @@ def _process_energy_low_rung(mf_dh, xc_list): f"non_exx): {xc_non_exx.token}") if len(xc_exx) != 0: xc_exx_extracted, eng_exx = _process_energy_exx(mf_dh, xc_exx) - xc_non_exx_extracted, eng_non_exx = _process_energy_low_rung(mf_dh, xc_non_exx) + eng_non_exx = _process_energy_low_rung(mf_dh, xc_non_exx) assert len(xc_exx_extracted) == 0 - if len(xc_non_exx_extracted) != 0: - raise RuntimeError( - f"Finally left some xc not parsed: {xc_non_exx_extracted.token}. " - f"This is probably bug.") eng_tot += eng_exx + eng_non_exx log.note(f"[RESULT] Energy of process_energy_low_rung (handle_multiple_omega): {eng_tot}") return eng_tot @@ -337,8 +333,7 @@ def _process_energy_low_rung(mf_dh, xc_list): "We evaluate each term one-by-one.") xc_non_exx_list = [XCList().build_from_list([xc_info]) for xc_info in xc_non_exx] for xc_non_exx_item in xc_non_exx_list: - xc_item_extracted, eng_item = _process_energy_low_rung(mf_dh, xc_non_exx_item) - assert len(xc_item_extracted) == 0 + eng_item = _process_energy_low_rung(mf_dh, xc_non_exx_item) eng_tot += eng_item log.note(f"[RESULT] Energy of process_energy_low_rung (handle_multiple_omega): {eng_tot}") return eng_tot -- Gitee From 1f8059feaae26d72e8b22453ef935729ff9e7761 Mon Sep 17 00:00:00 2001 From: ajz34 Date: Fri, 5 May 2023 22:44:28 +0800 Subject: [PATCH 12/12] feat(dh eng): unrestricted dRPA by ring-CCD --- pyscf/dh/__init__.py | 2 +- pyscf/dh/energy/__init__.py | 4 +- pyscf/dh/energy/ring_ccd/__init__.py | 3 +- pyscf/dh/energy/ring_ccd/rring_ccd.py | 35 ++-- .../dh/energy/ring_ccd/test/test_uring_ccd.py | 61 ++++++ pyscf/dh/energy/ring_ccd/uring_ccd.py | 186 ++++++++++++++++++ 6 files changed, 276 insertions(+), 15 deletions(-) create mode 100644 pyscf/dh/energy/ring_ccd/test/test_uring_ccd.py create mode 100644 pyscf/dh/energy/ring_ccd/uring_ccd.py diff --git a/pyscf/dh/__init__.py b/pyscf/dh/__init__.py index 22b0475..878bd37 100644 --- a/pyscf/dh/__init__.py +++ b/pyscf/dh/__init__.py @@ -5,7 +5,7 @@ from .energy import ( RHDFT, UHDFT, RMP2RIPySCF, RMP2ConvPySCF, UMP2ConvPySCF, RMP2Conv, RMP2RI, UMP2Conv, UMP2RI, RIEPARI, RIEPAConv, UIEPAConv, UIEPARI, - RRingCCDConv + RRingCCDConv, URingCCDConv ) from .energy import DH diff --git a/pyscf/dh/energy/__init__.py b/pyscf/dh/energy/__init__.py index 843baee..4edca96 100644 --- a/pyscf/dh/energy/__init__.py +++ b/pyscf/dh/energy/__init__.py @@ -3,11 +3,11 @@ from .dhbase import EngBase from .rhdft import RHDFT from .uhdft import UHDFT -from . import mp2, iepa +from . import mp2, iepa, ring_ccd from .mp2 import ( RMP2RIPySCF, RMP2ConvPySCF, UMP2ConvPySCF, RMP2Conv, RMP2RI, UMP2Conv, UMP2RI) from .iepa import (RIEPAConv, RIEPARI, UIEPAConv, UIEPARI) -from .ring_ccd import RRingCCDConv +from .ring_ccd import (RRingCCDConv, URingCCDConv) from .dh import DH diff --git a/pyscf/dh/energy/ring_ccd/__init__.py b/pyscf/dh/energy/ring_ccd/__init__.py index d1939c3..bec96ce 100644 --- a/pyscf/dh/energy/ring_ccd/__init__.py +++ b/pyscf/dh/energy/ring_ccd/__init__.py @@ -1,3 +1,4 @@ -from . import rring_ccd +from . import rring_ccd, uring_ccd from .rring_ccd import RRingCCDConv +from .uring_ccd import URingCCDConv diff --git a/pyscf/dh/energy/ring_ccd/rring_ccd.py b/pyscf/dh/energy/ring_ccd/rring_ccd.py index 1c677e6..88ea5f6 100644 --- a/pyscf/dh/energy/ring_ccd/rring_ccd.py +++ b/pyscf/dh/energy/ring_ccd/rring_ccd.py @@ -25,6 +25,8 @@ from pyscf.dh.util import pad_omega CONFIG_tol_eng_ring_ccd = getattr(__config__, "tol_eng_ring_ccd", 1e-8) CONFIG_tol_amp_ring_ccd = getattr(__config__, "tol_amp_ring_ccd", 1e-6) CONFIG_max_cycle_ring_ccd = getattr(__config__, "max_cycle_ring_ccd", 64) +CONFIG_diis_start_ring_ccd = getattr(__config__, "diis_start_ring_ccd", 3) +CONFIG_diis_space_ring_ccd = getattr(__config__, "diis_space_ring_ccd", 6) CONFIG_etb_first = getattr(__config__, "etb_first", False) @@ -36,6 +38,8 @@ class RingCCDBase(EngBase): self.conv_tol = CONFIG_tol_eng_ring_ccd self.conv_tol_amp = CONFIG_tol_amp_ring_ccd self.max_cycle = CONFIG_max_cycle_ring_ccd + self.diis_start = CONFIG_diis_start_ring_ccd + self.diis_space = CONFIG_diis_space_ring_ccd if with_df is None: with_df = getattr(self.scf, "with_df", None) if with_df is None: @@ -48,6 +52,7 @@ class RingCCDBase(EngBase): def kernel_energy_rring_ccd_conv( mo_energy, mo_coeff, eri_or_mol, mo_occ, tol_e=1e-8, tol_amp=1e-6, max_cycle=64, + diis_start=3, diis_space=6, verbose=lib.logger.NOTE): """ dRPA evaluation by ring-CCD with conventional integral. @@ -68,13 +73,17 @@ def kernel_energy_rring_ccd_conv( Threshold of L2 norm of ring-CCD amplitude while in DIIS update. max_cycle : int Maximum iteration of ring-CCD iteration. + diis_start : int + Start iteration number of DIIS. + diis_space : int + Space of DIIS. verbose : int Verbose level for PySCF. """ log = lib.logger.new_logger(verbose=verbose) - log.warn("Conventional integral of MP2 is not recommended!\n" - "Use density fitting approximation is recommended.") - log.info(f"[INFO] dRPA (ring-CCD) iteration, tol_e {tol_e:9.4e}, tol_amp {tol_amp:9.4e}, max_cycle {max_cycle:3d}") + # log.warn("Conventional integral of MP2 is not recommended!\n" + # "Use density fitting approximation is recommended.") + log.info(f"[INFO] Ring-CCD iteration, tol_e {tol_e:9.4e}, tol_amp {tol_amp:9.4e}, max_cycle {max_cycle:3d}") mask_occ = mo_occ != 0 mask_vir = mo_occ == 0 @@ -95,33 +104,35 @@ def kernel_energy_rring_ccd_conv( D = D_iajb.reshape(dim_ov, dim_ov) def update_T(T, B, D): - return - 1 / D * (B - 2 * B @ T - 2 * T @ B + 4 * T @ B @ T) + # return - 1 / D * (B - 2 * B @ T - 2 * T @ B + 4 * T @ B @ T) + BT = B @ T + return - 1 / D * (B - 2 * BT - 2 * BT.T + 4 * T @ BT) - # begin diis, start from third iteration, space of diis is 6 + # begin diis # T_old = np.zeros_like(B) T_new = np.zeros_like(B) eng_os = eng_ss = eng_drpa = 0 diis = lib.diis.DIIS() - diis.space = 6 + diis.space = diis_space converged = False for epoch in range(max_cycle): T_old, T_new = T_new, update_T(T_new, B, D) eng_old = eng_drpa - if epoch > 3: + if epoch > diis_start: T_new = diis.update(T_new) eng_os = eng_ss = - np.einsum("AB, AB ->", T_new, B) eng_drpa = eng_os + eng_ss err_amp = np.linalg.norm(T_new - T_old) err_e = abs(eng_drpa - eng_old) log.info( - f"[INFO] dRPA (ring-CCD) energy in iter {epoch:3d}: {eng_drpa:20.12f}, " + f"[INFO] Ring-CCD energy in iter {epoch:3d}: {eng_drpa:20.12f}, " f"amplitude L2 error {err_amp:9.4e}, eng_err {err_e:9.4e}") if err_amp < tol_amp and err_e < tol_e: converged = True - log.info("[INFO] dRPA (ring-CCD) amplitude converges.") + log.info("[INFO] Ring-CCD amplitude converges.") break if not converged: - log.warn(f"dRPA (ring-CCD) not converged in {max_cycle} iterations!") + log.warn(f"Ring-CCD not converged in {max_cycle} iterations!") results = dict() results["eng_corr_RING_CCD_OS"] = eng_os @@ -156,7 +167,7 @@ class RRingCCDConv(RingCCDBase): eri_or_mol = self.scf._eri if self.omega == 0 else mol eri_or_mol = eri_or_mol if eri_or_mol is not None else mol with mol.with_range_coulomb(self.omega): - results = kernel_energy_rring_ccd_conv( + results = self.kernel_energy_ring_ccd( mo_energy=mo_energy_act, mo_coeff=mo_coeff_act, eri_or_mol=eri_or_mol, @@ -164,6 +175,8 @@ class RRingCCDConv(RingCCDBase): tol_e=self.conv_tol, tol_amp=self.conv_tol_amp, max_cycle=self.max_cycle, + diis_start=self.diis_start, + diis_space=self.diis_space, verbose=self.verbose ) # pad omega diff --git a/pyscf/dh/energy/ring_ccd/test/test_uring_ccd.py b/pyscf/dh/energy/ring_ccd/test/test_uring_ccd.py new file mode 100644 index 0000000..d1f236b --- /dev/null +++ b/pyscf/dh/energy/ring_ccd/test/test_uring_ccd.py @@ -0,0 +1,61 @@ +import unittest +from pyscf import gto, scf, df, dft +from pyscf.dh import URingCCDConv + + +def get_mf_oh_hf_mrcc(): + coord = """ + H + O 1 R1 + """.replace("R1", "2.0") + + mol = gto.Mole(atom=coord, basis="cc-pVTZ", unit="AU", spin=1).build() + return scf.UHF(mol).density_fit(auxbasis="cc-pVTZ-jkfit").run() + + +def get_mf_oh_drpa75_mrcc(): + coord = """ + H + O 1 R1 + """.replace("R1", "2.0") + + mol = gto.Mole(atom=coord, basis="aug-cc-pVTZ", unit="AU", spin=1).build() + return dft.UKS(mol, xc="0.75*HF + 0.25*PBE, PBE").density_fit(auxbasis="aug-cc-pVTZ-jkfit").run() + + +class TestEngURingCCD(unittest.TestCase): + @classmethod + def setUpClass(cls) -> None: + cls.mf_oh_hf_mrcc = get_mf_oh_hf_mrcc() + cls.mf_oh_drpa75_mrcc = get_mf_oh_drpa75_mrcc() + + def test_eng_ring_ccd_known(self): + # mrcc: MINP_OH_cc-pVTZ_dRPA + mf_s = self.mf_oh_hf_mrcc + mol = mf_s.mol + # cheat to get ri-eri and evaluate by conv + mf_s._eri = df.DF(mol, auxbasis="cc-pVTZ-ri").get_ao_eri() + mf_dh = URingCCDConv(mf_s).run(frozen="FreezeNobleGasCore") + REF_MRCC = -0.254148360718 + self.assertAlmostEqual(mf_dh.results["eng_corr_RING_CCD"], REF_MRCC, 5) + + def test_eng_drpa75_known(self): + # mrcc: MINP_OH_aug-cc-pVTZ_dRPA75 + # a test case of how to evaluate doubly hybrid without handling xc code in package + mf_s = self.mf_oh_drpa75_mrcc + mol = mf_s.mol + # cheat to get ri-eri and evaluate by conv + mf_s._eri = df.DF(mol, auxbasis="aug-cc-pVTZ-ri").get_ao_eri() + mf_dh = URingCCDConv(mf_s).run(frozen="FreezeNobleGasCore") + mf_n = dft.RKS(mol).set(mo_coeff=mf_s.mo_coeff, mo_occ=mf_s.mo_occ, xc="0.75*HF + 0.25*PBE, ") + eng_low_rung = mf_n.energy_tot(dm=mf_s.make_rdm1()) + eng_rring_ccd = mf_dh.results["eng_corr_RING_CCD"] + eng_rring_ccd_os = mf_dh.results["eng_corr_RING_CCD_OS"] + eng_rring_ccd_ss = mf_dh.results["eng_corr_RING_CCD_SS"] + eng_tot = eng_low_rung + eng_rring_ccd + eng_scs = eng_low_rung + 1.5 * eng_rring_ccd_os + 0.5 * eng_rring_ccd_ss + REF_MRCC_dRPA75 = -75.683803588262 + REF_MRCC_SCS_dRPA75 = -75.678981497985 + self.assertAlmostEqual(eng_tot, REF_MRCC_dRPA75, 5) + self.assertAlmostEqual(eng_scs, REF_MRCC_SCS_dRPA75, 5) + diff --git a/pyscf/dh/energy/ring_ccd/uring_ccd.py b/pyscf/dh/energy/ring_ccd/uring_ccd.py new file mode 100644 index 0000000..c106952 --- /dev/null +++ b/pyscf/dh/energy/ring_ccd/uring_ccd.py @@ -0,0 +1,186 @@ +""" Unrestricted Ring-CCD. """ + + +from pyscf.dh.energy.ring_ccd.rring_ccd import RRingCCDConv +from pyscf.dh.util import pad_omega +from pyscf import ao2mo, lib +import numpy as np + +α, β = 0, 1 +αα, αβ, ββ = 0, 1, 2 + + +def kernel_energy_uring_ccd_conv( + mo_energy, mo_coeff, eri_or_mol, mo_occ, + tol_e=1e-8, tol_amp=1e-6, max_cycle=64, + diis_start=3, diis_space=6, + verbose=lib.logger.NOTE): + """ dRPA evaluation by ring-CCD with conventional integral (unrestricted). + + Parameters + ---------- + mo_energy : list[np.ndarray] + Molecular orbital energy levels. + mo_coeff : list[np.ndarray] + Molecular coefficients. + eri_or_mol : np.ndarray or gto.Mole + ERI that is recognized by ``pyscf.ao2mo.general``. + mo_occ : list[np.ndarray] + Molecular orbital occupation numbers. + + tol_e : float + Threshold of ring-CCD energy difference while in DIIS update. + tol_amp : float + Threshold of L2 norm of ring-CCD amplitude while in DIIS update. + max_cycle : int + Maximum iteration of ring-CCD iteration. + diis_start : int + Start iteration number of DIIS. + diis_space : int + Space of DIIS. + verbose : int + Verbose level for PySCF. + """ + log = lib.logger.new_logger(verbose=verbose) + # log.warn("Conventional integral of MP2 is not recommended!\n" + # "Use density fitting approximation is recommended.") + log.info(f"[INFO] Ring-CCD iteration, tol_e {tol_e:9.4e}, tol_amp {tol_amp:9.4e}, max_cycle {max_cycle:3d}") + + # preparation + mask_occ = [mo_occ[σ] != 0 for σ in (α, β)] + mask_vir = [mo_occ[σ] == 0 for σ in (α, β)] + eo = [mo_energy[σ][mask_occ[σ]] for σ in (α, β)] + ev = [mo_energy[σ][mask_vir[σ]] for σ in (α, β)] + Co = [mo_coeff[σ][:, mask_occ[σ]] for σ in (α, β)] + Cv = [mo_coeff[σ][:, mask_vir[σ]] for σ in (α, β)] + nocc = tuple([mask_occ[σ].sum() for σ in (α, β)]) + nvir = tuple([mask_vir[σ].sum() for σ in (α, β)]) + + # ao2mo preparation + log.info("[INFO] Start ao2mo") + g_ovov = [ + ao2mo.general(eri_or_mol, (Co[σ], Cv[σ], Co[ς], Cv[ς])).reshape(nocc[σ], nvir[σ], nocc[ς], nvir[ς]) + for (σ, ς) in ((α, α), (α, β), (β, β))] + D_ovov = [ + lib.direct_sum("i - a + j - b -> iajb", eo[σ], ev[σ], eo[ς], ev[ς]) + for (σ, ς) in ((α, α), (α, β), (β, β))] + log.info("[INFO] Finish ao2mo") + + # convert to dim_ov + dim_ov = [nocc[σ] * nvir[σ] for σ in (α, β)] + B = [g_ovov[σς].reshape(dim_ov[σ], dim_ov[ς]) for (σ, ς, σς) in ((α, α, αα), (α, β, αβ), (β, β, ββ))] + D = [D_ovov[σς].reshape(dim_ov[σ], dim_ov[ς]) for (σ, ς, σς) in ((α, α, αα), (α, β, αβ), (β, β, ββ))] + + def update_T(T, B, D): + B_αα, B_αβ, B_ββ = B + D_αα, D_αβ, D_ββ = D + T_αα, T_αβ, T_ββ = T + B_βα, D_βα, T_βα = B_αβ.T, D_αβ.T, T_αβ.T + BT_αα = B_αα @ T_αα + B_αβ @ T_βα + BT_αβ = B_αα @ T_αβ + B_αβ @ T_ββ + BT_βα = B_βα @ T_αα + B_ββ @ T_βα + BT_ββ = B_βα @ T_αβ + B_ββ @ T_ββ + # T_new_αα = - 1 / D_αα * ( + # + B_αα - B_αα @ T_αα - B_αβ @ T_βα - T_αα @ B_αα - T_αβ @ B_βα + # + T_αα @ B_αα @ T_αα + T_αα @ B_αβ @ T_βα + T_αβ @ B_βα @ T_αα + T_αβ @ B_ββ @ T_βα) + # T_new_αβ = - 1 / D_αβ * ( + # + B_αβ - B_αα @ T_αβ - B_αβ @ T_ββ - T_αα @ B_αβ - T_αβ @ B_ββ + # + T_αα @ B_αα @ T_αβ + T_αα @ B_αβ @ T_ββ + T_αβ @ B_βα @ T_αβ + T_αβ @ B_ββ @ T_ββ) + # T_new_ββ = - 1 / D_ββ * ( + # + B_ββ - B_ββ @ T_ββ - B_βα @ T_αβ - T_ββ @ B_ββ - T_βα @ B_αβ + # + T_ββ @ B_ββ @ T_ββ + T_ββ @ B_βα @ T_αβ + T_βα @ B_αβ @ T_ββ + T_βα @ B_αα @ T_αβ) + T_new_αα = - 1 / D_αα * (B_αα - BT_αα - BT_αα.T + T_αα @ BT_αα + T_αβ @ BT_βα) + T_new_αβ = - 1 / D_αβ * (B_αβ - BT_αβ - BT_βα.T + T_αα @ BT_αβ + T_αβ @ BT_ββ) + T_new_ββ = - 1 / D_ββ * (B_ββ - BT_ββ - BT_ββ.T + T_βα @ BT_αβ + T_ββ @ BT_ββ) + return [T_new_αα, T_new_αβ, T_new_ββ] + + # begin diis + T_new = [np.zeros_like(B[σς]) for σς in (αα, αβ, ββ)] + eng_os = eng_ss = eng_drpa = 0 + diis = [lib.diis.DIIS() for _ in (αα, αβ, ββ)] + for σς in (αα, αβ, ββ): + diis[σς].space = diis_space + converged = False + for epoch in range(max_cycle): + T_old, T_new = T_new, update_T(T_new, B, D) + eng_old = eng_drpa + if epoch > diis_start: + T_new = [diis[σς].update(T_new[σς]) for σς in (αα, αβ, ββ)] + eng_os = - np.einsum("AB, AB -> ", T_new[αβ], B[αβ]) + eng_ss = ( + - 0.5 * np.einsum("AB, AB -> ", T_new[αα], B[αα]) + - 0.5 * np.einsum("AB, AB -> ", T_new[ββ], B[ββ])) + eng_drpa = eng_os + eng_ss + err_amp = [np.linalg.norm(T_new[σς] - T_old[σς]) for σς in (αα, αβ, ββ)] + err_e = abs(eng_drpa - eng_old) + log.info( + f"[INFO] Ring-CCD energy in iter {epoch:3d}: {eng_drpa:20.12f}, " + f"amplitude L2 error {err_amp}, eng_err {err_e:9.4e}") + if max(err_amp) < tol_amp and err_e < tol_e: + converged = True + log.info("[INFO] Ring-CCD amplitude converges.") + break + if not converged: + log.warn(f"Ring-CCD not converged in {max_cycle} iterations!") + + results = dict() + results["eng_corr_RING_CCD_OS"] = eng_os + results["eng_corr_RING_CCD_SS"] = eng_ss + results["eng_corr_RING_CCD"] = eng_drpa + results["converged_RING_CCD"] = converged + log.info(f"[RESULT] Energy corr ring-CCD of same-spin: {eng_os :18.10f}") + log.info(f"[RESULT] Energy corr ring-CCD of oppo-spin: {eng_ss :18.10f}") + log.info(f"[RESULT] Energy corr ring-CCD of total : {eng_drpa:18.10f}") + return results + + +class URingCCDConv(RRingCCDConv): + + def driver_eng_ring_ccd(self, **kwargs): + mask = self.get_frozen_mask() + mol = self.mol + mo_coeff_act = [self.mo_coeff[σ][:, mask[σ]] for σ in (α, β)] + mo_energy_act = [self.mo_energy[σ][mask[σ]] for σ in (α, β)] + mo_occ_act = [self.mo_occ[σ][mask[σ]] for σ in (α, β)] + # eri generator + eri_or_mol = self.scf._eri if self.omega == 0 else mol + eri_or_mol = eri_or_mol if eri_or_mol is not None else mol + results = self.kernel_energy_ring_ccd( + mo_energy=mo_energy_act, + mo_coeff=mo_coeff_act, + eri_or_mol=eri_or_mol, + mo_occ=mo_occ_act, + tol_e=self.conv_tol, + tol_amp=self.conv_tol_amp, + max_cycle=self.max_cycle, + diis_start=self.diis_start, + diis_space=self.diis_space, + verbose=self.verbose + ) + # pad omega + results = {pad_omega(key, self.omega): val for (key, val) in results.items()} + self.results.update(results) + return results + + kernel_energy_ring_ccd = staticmethod(kernel_energy_uring_ccd_conv) + kernel = driver_eng_ring_ccd + + +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_r = RRingCCDConv(mf).run() + print(mf_r.results) + + mf = mf.to_uhf() + mf_u = URingCCDConv(mf).run() + print(mf_u.results) + + mol = gto.Mole(atom="O; H 1 0.94; H 1 0.94 2 104.5", basis="6-31G", charge=1, spin=1).build() + mf = scf.UHF(mol).run() + mf_u = URingCCDConv(mf).run(frozen=[[0], [1, 2]]) + print(mf_u.results) + + main_1() -- Gitee