diff --git a/pyscf/dh/__init__.py b/pyscf/dh/__init__.py index 22b0475601892a82430c79c27abbd7e93dcbc400..878bd37e6f947b3c946473b60dd978d4de2fcb51 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 843baee4b344edd466ef92168089525e85e5fd9c..4edca96b8c5d8d2b94591e520a814d4bcfc742cc 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/dh.py b/pyscf/dh/energy/dh.py index ef4535713edbb0008e5b708846c4251356c6ba0f..6025e65b669464793ea45268fdf0477912fcca46 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,41 +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") - # 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) + 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 @@ -114,12 +142,17 @@ def _process_energy_mp2(mf_dh, xc_list, force_evaluate=False): 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)) + mf_dh.inherited["mp2"][1].append(mf_mp2) eng_tot = comput_mp2() return eng_tot @@ -131,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] @@ -177,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: @@ -188,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]) @@ -237,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 @@ -249,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 @@ -275,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( @@ -291,23 +313,19 @@ 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, " 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 xc_extracted, eng_tot + return eng_tot else: log.warn( "Seems pure-DFT part has different values of omega.\n" @@ -315,16 +333,15 @@ def _process_energy_low_rung(mf_dh, xc_list, xc_to_parse=None): "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 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() @@ -332,7 +349,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.") @@ -346,8 +363,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) @@ -376,7 +393,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): @@ -390,6 +407,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 @@ -397,43 +416,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): @@ -459,7 +481,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() @@ -585,9 +607,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): @@ -596,12 +617,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 +630,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/rhdft.py b/pyscf/dh/energy/rhdft.py index af61d01da18991ebdb50dc6c18d42bb5cc82bb52..4a971e419d33bf0497d416bd09dd3a167fc2b5ff 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) diff --git a/pyscf/dh/energy/ring_ccd/__init__.py b/pyscf/dh/energy/ring_ccd/__init__.py index d1939c35ad5d1bd5478217f81e9c2abd61cb780a..bec96ce87ee3ae9152a0f058ab2dab7187862267 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 1c677e6be56db0875057d6f028ea565a1fd39fef..88ea5f68afa83e0565d0b310518ce8c1081e118e 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 0000000000000000000000000000000000000000..d1f236b142f882eb6f4e2e87039d8c1e7c96fb80 --- /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 0000000000000000000000000000000000000000..c106952e4c9f9719d79c2cd85c8b71e7d060fa19 --- /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()