diff --git a/pyscf/dh/energy/__init__.py b/pyscf/dh/energy/__init__.py index 4edca96b8c5d8d2b94591e520a814d4bcfc742cc..c7d865482b4ce4e557246241a14e5de56d5cf2cb 100644 --- a/pyscf/dh/energy/__init__.py +++ b/pyscf/dh/energy/__init__.py @@ -1,13 +1,11 @@ from .dhbase import EngBase -from .rhdft import RHDFT -from .uhdft import UHDFT - -from . import mp2, iepa, ring_ccd +from . import mp2, iepa, ring_ccd, hdft from .mp2 import ( RMP2RIPySCF, RMP2ConvPySCF, UMP2ConvPySCF, RMP2Conv, RMP2RI, UMP2Conv, UMP2RI) from .iepa import (RIEPAConv, RIEPARI, UIEPAConv, UIEPARI) from .ring_ccd import (RRingCCDConv, URingCCDConv) +from .hdft import RHDFT, UHDFT from .dh import DH diff --git a/pyscf/dh/energy/dh.py b/pyscf/dh/energy/dh.py index 6025e65b669464793ea45268fdf0477912fcca46..dbe6761c1e5836d7cc19966ac077b56852ed92cc 100644 --- a/pyscf/dh/energy/dh.py +++ b/pyscf/dh/energy/dh.py @@ -2,9 +2,9 @@ from pyscf.dh.energy import EngBase from typing import Tuple, List from pyscf.dh import util from pyscf.dh.util import XCType, XCList, XCDH, update_results, pad_omega -from pyscf.dh.energy.rhdft import get_rho, numint_customized +from pyscf.dh.energy.hdft.rhdft import get_rho, numint_customized from pyscf import lib, scf, gto, dft, df, __config__ -from pyscf.dh.energy.rhdft import custom_mf +from pyscf.dh.energy.hdft.rhdft import custom_mf CONFIG_etb_first = getattr(__config__, "etb_first", False) @@ -614,15 +614,17 @@ class DH(EngBase): def e_tot(self): return self.results[f"eng_dh_{self.xc.xc_eng.token}"] - def to_scf(self, **kwargs): + def to_scf(self, xc=None, **kwargs): # import if self.restricted: from pyscf.dh import RHDFT as HDFT else: from pyscf.dh import UHDFT as HDFT + xc = xc if xc is not None else self.xc.xc_scf + # generate instance - mf = HDFT.from_rdh(self, self.scf, self.xc.xc_scf, **kwargs) + mf = HDFT.from_cls(self, self.scf, xc, **kwargs) return mf @@ -641,9 +643,9 @@ class DH(EngBase): # generate instance if route_mp2.lower().startswith("ri"): - mf = MP2RI.from_rdh(self, self.scf, **kwargs) + mf = MP2RI.from_cls(self, self.scf, **kwargs) elif route_mp2.lower().startswith("conv"): - mf = MP2Conv.from_rdh(self, self.scf, **kwargs) + mf = MP2Conv.from_cls(self, self.scf, **kwargs) else: assert False, "Not recognized route_mp2." @@ -670,9 +672,9 @@ class DH(EngBase): # generate instance if route_iepa.lower().startswith("ri"): - mf = IEPARI.from_rdh(self, self.scf, **kwargs) + mf = IEPARI.from_cls(self, self.scf, **kwargs) elif route_iepa.lower().startswith("conv"): - mf = IEPAConv.from_rdh(self, self.scf, **kwargs) + mf = IEPAConv.from_cls(self, self.scf, **kwargs) else: assert False, "Not recognized route_iepa." @@ -699,7 +701,7 @@ class DH(EngBase): max_cycle_ring_ccd = self.flags.get("max_cycle_ring_ccd", NotImplemented) # generated instance - mf = RingCCDConv.from_rdh(self, self.scf, **kwargs) + mf = RingCCDConv.from_cls(self, self.scf, **kwargs) # fill configurations if tol_eng_ring_ccd is not NotImplemented: diff --git a/pyscf/dh/energy/dhbase.py b/pyscf/dh/energy/dhbase.py index db818836f2348e2f94a0bf88033f9549763a8be9..4ef5b691d96ab573e7bb52b15eb8972f356c4522 100644 --- a/pyscf/dh/energy/dhbase.py +++ b/pyscf/dh/energy/dhbase.py @@ -149,19 +149,23 @@ class EngBase(lib.StreamObject, ABC): return frozen_core.mask @classmethod - def from_rdh(cls, mf_dh: "EngBase", *args, **kwargs): + def from_cls(cls, mf_dh: "EngBase", *args, copy_all=False, **kwargs): """ Build a new (inherited) instance by a base EngBase instance. - Only some attributes are copied. Excluded attributes including + By default only some attributes are copied. Excluded attributes including - ``results`` - ``tensors`` - ``_tmpfile`` - ``e_corr`` """ mf_new = cls(*args, **kwargs) - mf_new.with_df = mf_dh.with_df - mf_new.frozen = mf_dh.frozen - mf_new.verbose = mf_dh.verbose - mf_new.max_memory = mf_dh.max_memory - mf_new.base = mf_dh + mf_new.__dict__.update(mf_dh.__dict__) + if not copy_all: + mf_new.results = dict() + mf_new.tensors = dict() + mf_new._tmpfile = lib.H5TmpFile() + mf_new.e_corr = NotImplemented return mf_new + + def to_resp(self): + raise NotImplementedError diff --git a/pyscf/dh/energy/hdft/__init__.py b/pyscf/dh/energy/hdft/__init__.py new file mode 100644 index 0000000000000000000000000000000000000000..93bffda9978c3cadd9da44f1dc6c77ab407def03 --- /dev/null +++ b/pyscf/dh/energy/hdft/__init__.py @@ -0,0 +1,4 @@ +from . import rhdft, uhdft + +from .rhdft import RHDFT +from .uhdft import UHDFT diff --git a/pyscf/dh/energy/rhdft.py b/pyscf/dh/energy/hdft/rhdft.py similarity index 91% rename from pyscf/dh/energy/rhdft.py rename to pyscf/dh/energy/hdft/rhdft.py index 4a971e419d33bf0497d416bd09dd3a167fc2b5ff..6a6e4b241cc9c342afe6a429861609fd146ebfdb 100644 --- a/pyscf/dh/energy/rhdft.py +++ b/pyscf/dh/energy/hdft/rhdft.py @@ -1,8 +1,11 @@ +from functools import cached_property + 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__ import numpy as np +import copy from types import MethodType CONFIG_ssr_x_fr = getattr(__config__, "ssr_x_fr", "LDA_X") @@ -279,13 +282,23 @@ def numint_customized(xc, _mol=None): def array_add_with_diff_rows(mat1, mat2): assert len(mat1.shape) == len(mat2.shape) - assert mat1.shape[1:] == mat2.shape[1:] - if len(mat1) >= len(mat2): - mat1[:len(mat2)] += mat2 - return mat1 + assert mat1.shape[-1] == mat2.shape[-1] + if len(mat1) < len(mat2): + mat1, mat2 = mat2, mat1 + + if mat1.ndim == 1: + mat1 += mat2 + elif mat1.ndim == 2: + mat1[:mat2.shape[0]] += mat2 + elif mat1.ndim == 3: + mat1[:mat2.shape[0], :mat2.shape[1]] += mat2 + elif mat1.ndim == 4: + mat1[:mat2.shape[0], :mat2.shape[1], :mat2.shape[2]] += mat2 + elif mat1.ndim == 5: + mat1[:mat2.shape[0], :mat2.shape[1], :mat2.shape[2], :mat2.shape[3]] += mat2 else: - mat2[:len(mat1)] += mat1 - return mat2 + raise NotImplementedError + return mat1 def eval_xc_eff(*args, **kwargs): exc, vxc, fxc, kxc = gen_lists[0](*args, **kwargs) @@ -333,6 +346,7 @@ def custom_mf(mf, xc, auxbasis_or_with_df=None): Note that if no with_df object passed in, the density-fitting setting of an SCF object is left as is. So leaving option ``auxbasis_or_with_df=None`` will not convert density-fitting SCF to conventional SCF. """ + mf = copy.copy(mf) verbose = mf.verbose log = lib.logger.new_logger(verbose=verbose) restricted = isinstance(mf, scf.hf.RHF) @@ -394,9 +408,14 @@ class RHDFT(EngBase): It's better to initialize this object first, before actually running SCF iterations. """ - def __init__(self, mf, xc): + def __init__(self, mf, xc=None): super().__init__(mf) - self.xc = xc # type: XCList + if xc is not None: + self.xc = xc + elif not hasattr(mf, "xc"): + self.xc = "HF" + else: + self.xc = self.scf.xc if isinstance(self.xc, str): xc_scf = XCList(self.xc, code_scf=True) xc_eng = XCList(self.xc, code_scf=False) @@ -405,15 +424,16 @@ class RHDFT(EngBase): self.xc = xc_scf else: xc_scf = self.xc - self._scf = custom_mf(mf, xc_scf) - @property + self.hdft = custom_mf(mf, xc_scf) + + @cached_property def restricted(self): - return isinstance(self.scf, scf.rhf.RHF) + return isinstance(self.hdft, scf.rhf.RHF) - @property + @cached_property def e_tot(self) -> float: - return self.scf.e_tot + return self.hdft.energy_tot() def make_energy_purexc(self, xc_lists, numint=None, dm=None): """ Evaluate energy contributions of pure (DFT) exchange-correlation effects. @@ -431,9 +451,9 @@ class RHDFT(EngBase): -------- get_energy_purexc """ - grids = self.scf.grids + grids = self.hdft.grids if dm is None: - dm = self.scf.make_rdm1() + dm = self.hdft.make_rdm1() dm = np.asarray(dm) if (self.restricted and dm.ndim != 2) or (not self.restricted and (dm.ndim != 3 or dm.shape[0] != 2)): raise ValueError("Dimension of input density matrix is not correct.") @@ -442,7 +462,13 @@ class RHDFT(EngBase): xc_lists, rho, grids.weights, self.restricted, numint=numint) def kernel(self, *args, **kwargs): - return self.scf.kernel(*args, **kwargs) + if not self.scf.converged: + self.scf.kernel(*args, **kwargs) + return self.e_tot + + def to_resp(self): + from pyscf.dh.response.hdft.rhdft import RHDFTResp + return RHDFTResp.from_cls(self, self.scf, copy_all=True) get_energy_exactx = staticmethod(get_energy_restricted_exactx) get_energy_noxc = staticmethod(get_energy_restricted_noxc) diff --git a/pyscf/dh/energy/uhdft.py b/pyscf/dh/energy/hdft/uhdft.py similarity index 97% rename from pyscf/dh/energy/uhdft.py rename to pyscf/dh/energy/hdft/uhdft.py index ecab7194c54e6d717b520b094cd822d37aa7b0e5..039a5c2d24643d87b1dfdfc910cb02a1a85870b0 100644 --- a/pyscf/dh/energy/uhdft.py +++ b/pyscf/dh/energy/hdft/uhdft.py @@ -1,6 +1,5 @@ -from pyscf.dh.energy.rhdft import RHDFT +from pyscf.dh.energy.hdft.rhdft import RHDFT from pyscf.dh import util -from pyscf import dft import numpy as np diff --git a/pyscf/dh/energy/mp2/rmp2.py b/pyscf/dh/energy/mp2/rmp2.py index 7404ebe970071e008e79824b864ec42b526efbf0..3fc90dd961b934276b0153aa8f97cbb8de3738b1 100644 --- a/pyscf/dh/energy/mp2/rmp2.py +++ b/pyscf/dh/energy/mp2/rmp2.py @@ -401,6 +401,10 @@ class RMP2RI(MP2Base): self.results.update(results) return results + def to_resp(self): + from pyscf.dh.response.mp2.rmp2ri import RMP2RespRI + return RMP2RespRI.from_cls(self, self.scf, copy_all=True) + kernel_energy_mp2 = staticmethod(kernel_energy_rmp2_ri_incore) kernel = driver_eng_mp2 diff --git a/pyscf/dh/energy/rresp.py b/pyscf/dh/energy/rresp.py deleted file mode 100644 index d6ba04ec54044df8feb9f51509e342261ae694e7..0000000000000000000000000000000000000000 --- a/pyscf/dh/energy/rresp.py +++ /dev/null @@ -1,98 +0,0 @@ -import numpy as np -from pyscf import lib - -from pyscf.dh import util - - -def kernel_resp_eri(params, Y_mo, nocc, hyb, - Y_mo_lr=None, alpha=None, - max_memory=2000, verbose=lib.logger.NOTE): - """ Prepare eri evluated in response function - - Cholesky decomposed ERIs in this function should be jk form - (density-fitting object used in SCF object). - - Parameters - ---------- - params : util.Params - Y_mo : np.ndarray - nocc : int - hyb : float - Y_mo_lr : np.ndarray - alpha : float - max_memory : float - verbose : int - - Returns - ------- - np.ndarray - """ - log = lib.logger.new_logger(verbose=verbose) - naux, nmo, _ = Y_mo.shape - nvir = nmo - nocc - is_rsh = Y_mo_lr is not None - einsum = lib.einsum - so, sv = slice(0, nocc), slice(nocc, nmo) - # sanity check - assert Y_mo.dtype is not np.complex - assert Y_mo.shape == (naux, nmo, nmo) - if is_rsh: - assert Y_mo_lr.shape == (naux, nmo, nmo) - assert alpha is not None - # create space bulk - incore_resp_eri = params.flags["incore_resp_eri"] - incore_resp_eri = util.parse_incore_flag( - incore_resp_eri, nocc**2 * nvir**2, max_memory) - log.debug("Actual incore strategy of incore_resp_eri: {:}".format(incore_resp_eri)) - log.debug("Creating `resp_eri`, shape {:}.".format((nvir, nocc, nvir, nocc))) - resp_eri = params.tensors.create("resp_eri", shape=(nvir, nocc, nvir, nocc), chunk=(1, 1, nvir, nocc)) - # take some parts of cholesky eri into memory - Y_vo = np.asarray(Y_mo[:, sv, so]) - Y_oo = np.asarray(Y_mo[:, so, so]) - Y_vo_lr = Y_oo_lr = None - if is_rsh: - Y_vo_lr = np.asarray(Y_mo_lr[:, sv, so]) - Y_oo_lr = np.asarray(Y_mo_lr[:, so, so]) - - nbatch = util.calc_batch_size( - nvir * naux + 3 * nocc ** 2 * nvir, max_memory, 2 * Y_mo.size) - - def save_eri_cpks(sA, buf): - resp_eri[sA] = buf - - with lib.call_in_background(save_eri_cpks) as async_write: - for sA in util.gen_batch(nocc, nmo, nbatch): - sAvir = slice(sA.start - nocc, sA.stop - nocc) - print(sAvir) - resp_eri_buf = ( - + 4 * einsum("Pai, Pbj -> aibj", Y_vo[:, sAvir], Y_vo) - - hyb * einsum("Paj, Pbi -> aibj", Y_vo[:, sAvir], Y_vo) - - hyb * einsum("Pij, Pab -> aibj", Y_oo, Y_mo[:, sA, sv])) - if is_rsh: - resp_eri_buf += ( - - (alpha - hyb) * einsum("Paj, Pbi -> aibj", Y_vo_lr[:, sAvir], Y_vo_lr) - - (alpha - hyb) * einsum("Pij, Pab -> aibj", Y_oo_lr, Y_mo_lr[:, sA, sv])) - async_write(sAvir, resp_eri_buf) - - return dict() - - -def resp_Ax0_HF(si, sa, sj, sb, hyb, Y_mo, - alpha=None, Y_mo_lr=None, - max_memory=2000, verbose=lib.logger.NOTE): - naux, nmo, _ = Y_mo.shape - ni, na = si.stop - si.start, sa.stop - sa.start - # sanity check - is_rsh = Y_mo_lr is not None - - def resp_Ax0_HF_inner(X): - X_shape = X.shape - X = X.reshape((-1, X_shape[-2], X_shape[-1])) - res = np.zeros((X.shape[0], ni, na)) - nbatch = util.calc_batch_size(nmo**2, max_memory, X.size + res.size) - for saux in util.gen_batch(0, naux, nbatch): - pass - - - - diff --git a/pyscf/dh/response/__init__.py b/pyscf/dh/response/__init__.py new file mode 100644 index 0000000000000000000000000000000000000000..eea327d4431e8f1ddb129133c1d7b41e02c4f063 --- /dev/null +++ b/pyscf/dh/response/__init__.py @@ -0,0 +1,4 @@ +from .respbase import RespBase +from . import dh, hdft, mp2, respbase + + diff --git a/pyscf/dh/response/dh/rdh.py b/pyscf/dh/response/dh/rdh.py new file mode 100644 index 0000000000000000000000000000000000000000..d8bc1078e3bf9156294a24a4cf04dd41f84be41a --- /dev/null +++ b/pyscf/dh/response/dh/rdh.py @@ -0,0 +1,229 @@ +""" Doubly-Hybrid Response-Related Utilities. """ +import numpy as np + +from pyscf.dh.energy.dh import RDH +from pyscf.dh.response import RespBase +from pyscf.dh.response.hdft.rhdft import RHDFTResp +from pyscf.dh.response.respbase import get_rdm1_resp_vo_restricted + + +class RDHResp(RDH, RespBase): + + def __init__(self, *args, **kwargs): + super().__init__(*args, **kwargs) + + # generate response instance for SCF (to obtain Ax0_Core) + HDFTResp = RHDFTResp if self.restricted else NotImplemented + self.scf_resp = HDFTResp(self.scf) + self._inherited_updated = False + + @property + def Ax0_Core(self): + """ Fock response of underlying SCF object in MO basis. """ + if self._Ax0_Core is NotImplemented: + self._Ax0_Core = self.scf_resp.Ax0_Core + return self._Ax0_Core + + @Ax0_Core.setter + def Ax0_Core(self, Ax0_Core): + self._Ax0_Core = Ax0_Core + + # in response instance, we first transfer all child instances into response + def to_scf(self, *args, **kwargs): + mf_resp = super().to_scf(*args, **kwargs).to_resp() + mf_resp.Ax0_Core = self.scf_resp.Ax0_Core + return mf_resp + + def to_mp2(self, *args, **kwargs): + assert len(args) == 0 + mf_resp = super().to_mp2(**kwargs).to_resp() + mf_resp.Ax0_Core = self.scf_resp.Ax0_Core + return mf_resp + + def to_iepa(self, *args, **kwargs): + raise NotImplementedError("Not implemented for response functions.") + + def to_ring_ccd(self, *args, **kwargs): + raise NotImplementedError("Not implemented for response functions.") + + def update_inherited(self): + """ Update inherited attribute, to transfer all child method instances to response. """ + + if self._inherited_updated: + return + + # if energy not evaluated, then evaluate energy first + if len(self.inherited) == 0: + self.kernel() + + # for energy evaluation, instance of low_rung may not be generated. + if len(self.inherited["low_rung"][1]) == 0: + HDFTResp = RHDFTResp if self.restricted else NotImplemented + self.inherited["low_rung"][1].append(HDFTResp(self.scf, xc=self.inherited["low_rung"][0])) + + # transform instances to response functions + # note that if Ax0_Core appears, then this object is already response, or have been inherited + for key in self.inherited: + for idx in range(len(self.inherited[key][1])): + instance = self.inherited[key][1][idx] + if not hasattr(instance, "Ax0_Core"): + instance = instance.to_resp() + self.inherited[key][1][idx] = instance + + self._inherited_updated = True + + def make_lag_vo(self): + if "lag_vo" in self.tensors: + return self.tensors["lag_vo"] + + self.update_inherited() + + nocc, nvir, nmo = self.nocc, self.nvir, self.nmo + lag_vo = np.zeros((nvir, nocc)) + for key in self.inherited: + for instance in self.inherited[key][1]: + lag_vo += instance.make_lag_vo() + + self.tensors["lag_vo"] = lag_vo + return lag_vo + + def make_rdm1_resp(self, ao=False): + if "rdm1_resp" in self.tensors: + return self.tensors["rdm1_resp"] + + self.update_inherited() + + rdm1_resp = np.diag(self.mo_occ) + for key in self.inherited: + for instance in self.inherited[key][1]: + if hasattr(instance, "make_rdm1_corr"): + rdm1_resp += instance.make_rdm1_corr() + + nocc, nvir, nmo = self.nocc, self.nvir, self.nmo + so, sv = slice(0, nocc), slice(nocc, nmo) + mo_energy = self.mo_energy + mo_occ = self.mo_occ + lag_vo = self.make_lag_vo() + max_cycle = self.max_cycle_cpks + tol = self.tol_cpks + Ax0_Core = self.Ax0_Core + verbose = self.verbose + + rdm1_resp_vo = get_rdm1_resp_vo_restricted( + lag_vo, mo_energy, mo_occ, Ax0_Core, + max_cycle=max_cycle, tol=tol, verbose=verbose) + + rdm1_resp[sv, so] += rdm1_resp_vo + + self.tensors["rdm1_resp"] = rdm1_resp + + if ao: + rdm1_resp = self.mo_coeff @ rdm1_resp @ self.mo_coeff.T + return rdm1_resp + + def make_dipole(self): + # prepare input + mol = self.mol + rdm1_ao = self.make_rdm1_resp(ao=True) + int1e_r = mol.intor("int1e_r") + + dip_elec = - np.einsum("uv, tuv -> t", rdm1_ao, int1e_r) + dip_nuc = np.einsum("A, At -> t", mol.atom_charges(), mol.atom_coords()) + dip = dip_elec + dip_nuc + self.tensors["dipole"] = dip + return dip + + +if __name__ == '__main__': + def main_1(): + # test energy is the same for response and general DH + from pyscf import gto, scf + mol = gto.Mole(atom="O; H 1 0.94; H 1 0.94 2 104.5", basis="6-31G").build() + mf_resp = RDHResp(mol, xc="XYG3").run() + mf = RDH(mol, xc="XYG3").run() + assert np.allclose(mf_resp.e_tot, mf.e_tot) + + def main_2(): + # test RSH functional + from pyscf import gto, scf + mol = gto.Mole(atom="O; H 1 0.94; H 1 0.94 2 104.5", basis="6-31G").build() + mf = RDH(mol, xc="RS-PBE-P86").run() + print(mf.results) + mf_resp = RDHResp(mol, xc="RS-PBE-P86").run() + print(mf_resp.results) + assert np.allclose(mf_resp.e_tot, mf.e_tot) + + def main_3(): + # test unbalanced OS contribution + from pyscf import gto, scf + mol = gto.Mole(atom="O; H 1 0.94; H 1 0.94 2 104.5", basis="6-31G").build() + mf = RDH(mol, xc="XYGJ-OS").run() + mf_resp = RDHResp(mol, xc="XYGJ-OS").run() + assert np.allclose(mf_resp.e_tot, mf.e_tot) + + def main_4(): + # try response density and dipole + from pyscf import gto, scf, dft + + np.set_printoptions(8, suppress=True, linewidth=150) + + mol = gto.Mole(atom="O; H 1 0.94; H 1 0.94 2 104.5", basis="6-31G").build() + mf_resp = RDHResp(mol, xc="XYG3").run() + + def eng_with_dipole_field(t, h): + mf_scf = dft.RKS(mol, xc="B3LYPg").density_fit("cc-pVDZ-jkfit") + mf_scf.get_hcore = lambda *args, **kwargs: scf.hf.get_hcore(mol) - h * mol.intor("int1e_r")[t] + mf_scf.run(conv_tol=1e-12) + mf_mp = RDH(mf_scf, xc="XYG3").run() + return mf_mp.e_tot + + eng_array = np.zeros((2, 3)) + h = 1e-4 + for idx, h in [(0, h), [1, -h]]: + for t in (0, 1, 2): + eng_array[idx, t] = eng_with_dipole_field(t, h) + dip_elec_num = - (eng_array[0] - eng_array[1]) / (2 * h) + + dip_nuc = np.einsum("A, At -> t", mol.atom_charges(), mol.atom_coords()) + dip_anal = mf_resp.make_dipole() + dip_num = dip_elec_num + dip_nuc + np.allclose(dip_anal, dip_num, rtol=1e-5, atol=1e-6) + print(dip_anal) + print(dip_num) + + def main_5(): + # try response density and dipole + from pyscf import gto, scf, dft + + np.set_printoptions(8, suppress=True, linewidth=150) + + mol = gto.Mole(atom="O; H 1 0.94; H 1 0.94 2 104.5", basis="6-31G").build() + xc_code = "0.5*HF + 0.5*LR_HF(0.7) + 0.5*SSR(GGA_X_PBE, 0.7), 0.75*SSR(GGA_C_P86, 0.7) + MP2(0.5, 0.5) + RS_MP2(-0.7, -0.25, -0.25) + RS_MP2(0.7, 0.5, 0.5)" + mf_resp = RDHResp(mol, xc=xc_code).run() + + def eng_with_dipole_field(t, h): + mf_mp = RDH(mol, xc=xc_code) + mf_mp.scf.get_hcore = lambda *args, **kwargs: scf.hf.get_hcore(mol) - h * mol.intor("int1e_r")[t] + mf_mp.scf.conv_tol = 1e-12 + mf_mp.run() + return mf_mp.e_tot + + eng_array = np.zeros((2, 3)) + h = 1e-4 + for idx, h in [(0, h), [1, -h]]: + for t in (0, 1, 2): + eng_array[idx, t] = eng_with_dipole_field(t, h) + dip_elec_num = - (eng_array[0] - eng_array[1]) / (2 * h) + + dip_nuc = np.einsum("A, At -> t", mol.atom_charges(), mol.atom_coords()) + dip_anal = mf_resp.make_dipole() + dip_num = dip_elec_num + dip_nuc + np.allclose(dip_anal, dip_num, rtol=1e-5, atol=1e-6) + print(dip_anal) + print(dip_num) + + # main_1() + # main_2() + # main_3() + # main_4() + main_5() diff --git a/pyscf/dh/response/hdft/__init__.py b/pyscf/dh/response/hdft/__init__.py new file mode 100644 index 0000000000000000000000000000000000000000..253f044a9f1f538e0346898695b48d7771c26799 --- /dev/null +++ b/pyscf/dh/response/hdft/__init__.py @@ -0,0 +1 @@ +from pyscf.dh.response.hdft.rhdft import RHDFTResp diff --git a/pyscf/dh/response/hdft/rhdft.py b/pyscf/dh/response/hdft/rhdft.py new file mode 100644 index 0000000000000000000000000000000000000000..9ce8aca6e54ebc91b259f81ff8cb61bf757e388d --- /dev/null +++ b/pyscf/dh/response/hdft/rhdft.py @@ -0,0 +1,548 @@ +""" Hybrid-DFT Response-Related Utilities. """ + +from pyscf.dh import RHDFT +from pyscf.dh import util +from pyscf import gto, dft, lib, __config__, scf +from pyscf.dh.response import RespBase +from pyscf.scf import _response_functions # this import is not unnecessary +from pyscf.dh.energy.hdft.rhdft import get_rho +import numpy as np +from functools import cached_property + + +CONFIG_dh_verbose = getattr(__config__, "dh_verbose", lib.logger.NOTE) +CONFIG_incore_cderi_uaa_hdft = getattr(__config__, "incore_cderi_uaa_hdft", "auto") +CONFIG_incore_eri_cpks_vovo = getattr(__config__, "incore_eri_cpks_vovo", "auto") +CONFIG_dft_gen_grid_Grids_level = getattr(__config__, 'dft_gen_grid_Grids_level', 3) +CONFIG_use_eri_cpks = getattr(__config__, "use_eri_cpks", True) + + +def get_eri_cpks_vovo( + cderi_uaa, mo_occ, cj, ck, eri_cpks_vovo=None, + max_memory=2000, verbose=CONFIG_dh_verbose): + r""" Obtain ERI term :math:`A_{ai, bj}^\mathrm{HF}`. + + For given + + Parameters + ---------- + cderi_uaa : np.ndarray + mo_occ : np.ndarray + cj : float + ck : float + eri_cpks_vovo : np.ndarray or None + max_memory : float + verbose : int + + Returns + ------- + np.ndarray + """ + log = lib.logger.new_logger(verbose=verbose) + time0 = lib.logger.process_clock(), lib.logger.perf_counter() + + # dimension and sanity check + naux = cderi_uaa.shape[0] + nmo = len(mo_occ) + nocc = (mo_occ != 0).sum() + nvir = (mo_occ == 0).sum() + so, sv = slice(0, nocc), slice(nocc, nmo) + assert cderi_uaa.shape == (naux, nmo, nmo) + + # allocate eri_cpks_vovo if necessary + if eri_cpks_vovo is None: + eri_cpks_vovo = np.zeros((nvir, nocc, nvir, nocc)) + else: + assert eri_cpks_vovo.shape == (nvir, nocc, nvir, nocc) + + # copy some tensors to memory + cderi_uvo = np.asarray(cderi_uaa[:, sv, so]) + cderi_uoo = np.asarray(cderi_uaa[:, so, so]) + + # prepare for async and batch + mem_avail = max_memory - lib.current_memory()[0] + nbatch = util.calc_batch_size(nvir*naux + 2*nocc**2*nvir, mem_avail) + + def load_cderi_uaa(slc): + return np.asarray(cderi_uaa[:, slc, sv]) + + def write_eri_cpks_vovo(slc, blk): + eri_cpks_vovo[slc] += blk + + batches = util.gen_batch(nocc, nmo, nbatch) + with lib.call_in_background(write_eri_cpks_vovo) as async_write_eri_cpks_vovo: + for sA, cderi_uAa in zip(batches, lib.map_with_prefetch(load_cderi_uaa, batches)): + log.debug(f"[DEBUG] in get_eri_cpks_vovo, slice {sA} of virtual orbitals {nocc, nmo}") + sAvir = slice(sA.start - nocc, sA.stop - nocc) + blk = np.zeros((sA.stop - sA.start, nocc, nvir, nocc)) + if abs(cj) > 1e-10: + blk += 4 * cj * lib.einsum("Pai, Pbj -> aibj", cderi_uvo[:, sAvir], cderi_uvo) + if abs(ck) > 1e-10: + blk -= ck * ( + + lib.einsum("Paj, Pbi -> aibj", cderi_uvo[:, sAvir], cderi_uvo) + + lib.einsum("Pij, Pab -> aibj", cderi_uoo, cderi_uAa)) + async_write_eri_cpks_vovo(sAvir, blk) + + log.timer("get_eri_cpks_vovo", *time0) + return eri_cpks_vovo + + +def Ax0_cpks_HF(eri_cpks_vovo, max_memory=2000, verbose=CONFIG_dh_verbose): + r""" Convenient function for evaluation of HF contribution of Fock response in MO basis + :math:`\sum_{rs} A_{ai, bj} X_{bj}^\mathbb{A}` by explicitly contraction to MO ERI :math:`(ai|bj)`. + + Parameters + ---------- + eri_cpks_vovo : np.ndarray + max_memory : float + verbose : int + + Returns + ------- + callable + """ + nvir, nocc = eri_cpks_vovo.shape[:2] + + def Ax0_cpks_HF_inner(X): + log = lib.logger.new_logger(verbose=verbose) + time0 = lib.logger.process_clock(), lib.logger.perf_counter() + + X_shape = X.shape + X = X.reshape((-1, X_shape[-2], X_shape[-1])) + res = np.zeros_like(X) + + def load_eri_cpks_vovo(slc): + return eri_cpks_vovo[slc] + + mem_avail = max_memory - lib.current_memory()[0] + nbatch = util.calc_batch_size(nocc**2 * nvir, mem_avail) + batches = util.gen_batch(0, nvir, nbatch) + + for sA, eri_cpks_Vovo in zip(batches, lib.map_with_prefetch(load_eri_cpks_vovo, batches)): + res[:, sA] = lib.einsum("aibj, Abj -> Aai", eri_cpks_Vovo, X) + + res.shape = list(X_shape[:-2]) + [res.shape[-2], res.shape[-1]] + + log.timer("Ax0_cpks_HF_inner", *time0) + return res + + return Ax0_cpks_HF_inner + + +def Ax0_Core_KS( + sp, sq, sr, ss, + mo_coeff, xc_setting, xc_kernel, + max_memory=2000, verbose=CONFIG_dh_verbose): + r""" Convenient function for evaluation of pure DFT contribution of Fock response in MO basis + :math:`\sum_{rs} A_{pq, rs} X_{rs}^\mathbb{A}`. + + Parameters + ---------- + sp, sq, sr, ss : slice or list + mo_coeff : np.ndarray + xc_setting : tuple[dft.numint.NumInt, gto.Mole, dft.Grids, str, np.ndarray] + xc_kernel : tuple[np.ndarray, np.ndarray, np.ndarray] + max_memory : float + verbose : int + + Returns + ------- + callable + """ + C = mo_coeff + ni, mol, grids, xc, dm = xc_setting + rho, _, fxc = xc_kernel + + def Ax0_Core_KS_inner(X): + log = lib.logger.new_logger(verbose=verbose) + time0 = lib.logger.process_clock(), lib.logger.perf_counter() + mem_avail = max_memory - lib.current_memory()[0] + + X_shape = X.shape + X = X.reshape((-1, X_shape[-2], X_shape[-1])) + dmX = C[:, sr] @ X @ C[:, ss].T + dmX += dmX.swapaxes(-1, -2) + ax_ao = ni.nr_rks_fxc(mol, grids, xc, dm, dmX, hermi=1, rho0=rho, fxc=fxc, max_memory=mem_avail) + res = 2 * C[:, sp].T @ ax_ao @ C[:, sq] + res.shape = list(X_shape[:-2]) + [res.shape[-2], res.shape[-1]] + + log.timer("Ax0_Core_KS_inner", *time0) + return res + + return Ax0_Core_KS_inner + + +def Ax0_Core_resp( + sp, sq, sr, ss, vresp, mo_coeff, + verbose=CONFIG_dh_verbose): + r""" Convenient function for evaluation of Fock response in MO basis + :math:`\sum_{rs} A_{pq, rs} X_{rs}^\mathbb{A}` by PySCF's response function. + + 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. + verbose : int + Print verbosity. + + 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): + log = lib.logger.new_logger(verbose=verbose) + time0 = lib.logger.process_clock(), lib.logger.perf_counter() + + X_shape = X.shape + X = X.reshape((-1, X_shape[-2], X_shape[-1])) + dmX = C[:, sr] @ X @ C[:, ss].T + 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]] + + log.timer("Ax0_Core_resp_inner", *time0) + return res + + return Ax0_Core_resp_inner + + +def get_xc_integral(ni, mol, grids, xc, dm): + rho = get_rho(mol, grids, dm) + try: + exc, vxc, fxc, kxc = ni.eval_xc_eff(xc, rho, deriv=3) + tensors = { + "rho": rho, + f"exc_{xc}": exc, + f"vxc_{xc}": vxc, + f"fxc_{xc}": fxc, + f"kxc_{xc}": kxc, + } + except NotImplementedError: + exc, vxc, fxc, _ = ni.eval_xc_eff(xc, rho, deriv=2) + tensors = { + "rho": rho, + f"exc_{xc}": exc, + f"vxc_{xc}": vxc, + f"fxc_{xc}": fxc, + } + return tensors + + +class RHDFTResp(RHDFT, RespBase): + + def __init__(self, *args, **kwargs): + super().__init__(*args, **kwargs) + self.incore_cderi_uaa = CONFIG_incore_cderi_uaa_hdft + self.incore_eri_cpks_vovo = CONFIG_incore_eri_cpks_vovo + self.use_eri_cpks = CONFIG_use_eri_cpks + + grid_level_cpks = max(CONFIG_dft_gen_grid_Grids_level - 1, 1) + self.grids_cpks = dft.Grids(self.mol) + self.grids_cpks.level = grid_level_cpks + self.grids_cpks.build() + + @cached_property + def vresp(self): + """ Fock response function (derivative w.r.t. molecular coefficient in AO basis). """ + try: + return self.scf.gen_response() + except ValueError: + # case that have customized xc + # should pass check of ni.libxc.test_deriv_order and ni.libxc.is_hybrid_xc + ni = self.scf._numint + omega, alpha, hyb = ni.rsh_and_hybrid_coeff(self.scf.xc, self.mol.spin) + if abs(hyb) > 1e-10 or abs(omega) > 1e-10: + fake_xc = "B3LYPg" + else: + fake_xc = "PBE" + actual_xc = self.scf.xc + self.scf.xc = fake_xc + resp = self.scf.gen_response() + self.scf.xc = actual_xc + return resp + + def make_cderi_uaa(self, omega=0): + """ Generate cholesky decomposed ERI (all block, full orbitals, s1 symm, in memory/disk). """ + if util.pad_omega("cderi_uaa", omega) in self.tensors: + return self.tensors[util.pad_omega("cderi_uaa", omega)] + + log = lib.logger.new_logger(verbose=self.verbose) + + # dimension and mask + mo_coeff = self.mo_coeff + nmo = self.nmo + + # density fitting preparation + with_df = util.get_with_df_omega(self.scf.with_df, omega) + naux = with_df.get_naoaux() + + # array allocation + mem_avail = self.max_memory - lib.current_memory()[0] + incore_cderi_uaa = self.incore_cderi_uaa + cderi_uaa = util.allocate_array( + incore_cderi_uaa, (naux, nmo, nmo), mem_avail, + h5file=self._tmpfile, name=util.pad_omega("cderi_uaa", omega), zero_init=False) + log.info(f"[INFO] Store type of cderi_uaa: {type(cderi_uaa)}") + + # generate array + util.get_cderi_mo(with_df, mo_coeff, cderi_uaa, max_memory=self.max_memory) + + tensors = {util.pad_omega("cderi_uaa", omega): cderi_uaa} + self.tensors.update(tensors) + return cderi_uaa + + def get_Ax0_Core(self, sp, sq, sr, ss): + r""" Convenient function for evaluation of Fock response in MO basis + :math:`\sum_{rs} A_{pq, rs} X_{rs}^\mathbb{A}`. + + 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. + """ + # if not RI, then use general Ax0_Core_resp + if not hasattr(self.scf, "with_df") or not self.use_eri_cpks: + return self.get_Ax0_Core_resp(sp, sq, sr, ss) + + # try if satisfies CPKS evaluation (vovo) + lst_nmo = np.arange(self.nmo) + lst_occ = lst_nmo[:self.nocc] + lst_vir = lst_nmo[self.nocc:] + try: + if ( + np.all(lst_nmo[sp] == lst_vir) and np.all(lst_nmo[sq] == lst_occ) and + np.all(lst_nmo[sr] == lst_vir) and np.all(lst_nmo[ss] == lst_occ)): + return self.get_Ax0_cpks() + else: + # otherwise, use response by PySCF + return self.get_Ax0_Core_resp(sp, sq, sr, ss) + except ValueError: + # dimension not match + return self.get_Ax0_Core_resp(sp, sq, sr, ss) + + def get_Ax0_Core_resp(self, sp, sq, sr, ss, vresp=None, mo_coeff=None): + r""" Convenient function for evaluation of Fock response in MO basis + :math:`\sum_{rs} A_{pq, rs} X_{rs}^\mathbb{A}` by PySCF's response function. + + 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) + + def make_eri_cpks_vovo(self): + r""" Generate ERI for CP-KS evaluation :math:`(ai, bj)` for current exch-corr setting. """ + if "eri_cpks_vovo" in self.tensors: + return self.tensors["eri_cpks_vovo"] + + # This function will utilize mf.scf._numint to specify hybrid and omega coefficients + + if not hasattr(self.scf, "xc"): + # not a DFT instance, configure as HF instance + omega, alpha, hyb = 0, 0, 1 + else: + ni = self.scf._numint # type: dft.numint.NumInt + omega, alpha, hyb = ni.rsh_and_hybrid_coeff(self.scf.xc) + + nvir, nocc = self.nvir, self.nocc + cderi_uaa = self.tensors.get("cderi_uaa", self.make_cderi_uaa()) + + eri_cpks_vovo = util.allocate_array( + incore=self.incore_eri_cpks_vovo, + shape=(nvir, nocc, nvir, nocc), + max_memory=self.max_memory, + h5file=self._tmpfile, + name="eri_cpks_vovo", + chunk=(1, nocc, nvir, nocc)) + + get_eri_cpks_vovo( + cderi_uaa=cderi_uaa, + mo_occ=self.mo_occ, + cj=1, ck=hyb, + eri_cpks_vovo=eri_cpks_vovo, + max_memory=self.max_memory, + verbose=self.verbose + ) + + if abs(omega) > 1e-10: + cderi_uaa = self.make_cderi_uaa(omega=omega) + get_eri_cpks_vovo( + cderi_uaa=cderi_uaa, + mo_occ=self.mo_occ, + cj=0, ck=alpha-hyb, + eri_cpks_vovo=eri_cpks_vovo, + max_memory=self.max_memory, + verbose=self.verbose + ) + + self.tensors["eri_cpks_vovo"] = eri_cpks_vovo + return eri_cpks_vovo + + def get_Ax0_cpks_HF(self): + eri_cpks_vovo = self.tensors.get("eri_cpks_vovo", self.make_eri_cpks_vovo()) + ax0_cpks_hf = Ax0_cpks_HF(eri_cpks_vovo, self.max_memory, self.verbose) + return ax0_cpks_hf + + def make_xc_integral(self): + ni = self.scf._numint # type: dft.numint.NumInt + mol = self.mol + grids = self.scf.grids + xc_token = self.xc.token + dm = self.scf.make_rdm1() + tensors = get_xc_integral(ni, mol, grids, xc_token, dm) + self.tensors.update(tensors) + return tensors + + def make_xc_integral_cpks(self): + ni = self.scf._numint # type: dft.numint.NumInt + mol = self.mol + grids = self.grids_cpks + xc_token = self.xc.token + dm = self.scf.make_rdm1() + tensors = get_xc_integral(ni, mol, grids, xc_token, dm) + tensors = {key + "/cpks": val for key, val in tensors.items()} + self.tensors.update(tensors) + return tensors + + def get_Ax0_Core_KS(self, sp, sq, sr, ss): + if not hasattr(self.scf, "xc"): + # not a DFT instance, KS contribution is zero + return lambda *args, **kwargs: 0 + + ni = self.scf._numint + mol = self.mol + grids = self.grids_cpks + xc_token = self.xc.token + dm = self.scf.make_rdm1() + if f"fxc_{xc_token}" not in self.tensors: + self.make_xc_integral_cpks() + rho = self.tensors[f"rho/cpks"] + vxc = self.tensors[f"vxc_{xc_token}/cpks"] + fxc = self.tensors[f"fxc_{xc_token}/cpks"] + + xc_setting = ni, mol, grids, xc_token, dm + xc_kernel = rho, vxc, fxc + mo_coeff = self.mo_coeff + + ax0_core_ks = Ax0_Core_KS( + sp, sq, sr, ss, + mo_coeff, xc_setting, xc_kernel, + max_memory=self.max_memory, + verbose=self.verbose) + return ax0_core_ks + + def get_Ax0_cpks(self): + nocc, nmo = self.nocc, self.nmo + so, sv = slice(0, nocc), slice(nocc, nmo) + ax0_core_ks = self.get_Ax0_Core_KS(sv, so, sv, so) + ax0_cpks_hf = self.get_Ax0_cpks_HF() + + def Ax0_cpks_inner(X): + res = ax0_cpks_hf(X) + ax0_core_ks(X) + return res + return Ax0_cpks_inner + + @property + def Ax0_Core(self): + """ Fock response of underlying SCF object in MO basis. """ + if self._Ax0_Core is NotImplemented: + self._Ax0_Core = self.get_Ax0_Core + return self._Ax0_Core + + @Ax0_Core.setter + def Ax0_Core(self, Ax0_Core): + self._Ax0_Core = Ax0_Core + + def make_lag_vo(self): + r""" Generate hybrid DFT contribution to Lagrangian vir-occ block :math:`L_{ai}`. """ + if "lag_vo" in self.tensors: + return self.tensors["lag_vo"] + + nocc, nmo = self.nocc, self.nmo + so, sv = slice(0, nocc), slice(nocc, nmo) + mo_coeff = self.mo_coeff + lag_vo = 4 * mo_coeff[:, sv].T @ self.hdft.get_fock(dm=self.scf.make_rdm1()) @ mo_coeff[:, so] + self.tensors["lag_vo"] = lag_vo + return lag_vo + + +if __name__ == '__main__': + def main_1(): + # test that RSH functional is correct for Ax0_cpks + from pyscf import gto, scf + np.set_printoptions(5, suppress=True) + np.random.seed(0) + + mol = gto.Mole(atom="O; H 1 0.94; H 1 0.94 2 104.5", basis="6-31G").build() + mf = dft.RKS(mol, xc="CAM-B3LYP").density_fit().run() + mf_scf = RHDFTResp(mf) + mf_scf.grids_cpks = mf_scf.scf.grids + + nocc, nvir, nmo = mf_scf.nocc, mf_scf.nvir, mf_scf.nmo + so, sv = slice(0, nocc), slice(nocc, nmo) + X = np.random.randn(3, nvir, nocc) + ax_cpks = mf_scf.get_Ax0_cpks()(X) + ax_core = mf_scf.get_Ax0_Core_resp(sv, so, sv, so)(X) + print(np.allclose(ax_cpks, ax_core)) + # print(ax_core[0]) + # print(ax_cpks[0]) + + def main_2(): + # test that meta-GGA functional is correct for Ax0_cpks + from pyscf import gto, scf + np.set_printoptions(5, suppress=True) + np.random.seed(0) + + mol = gto.Mole(atom="O; H 1 0.94; H 1 0.94 2 104.5", basis="6-31G").build() + mf = dft.RKS(mol, xc="TPSS0").density_fit().run() + mf_scf = RHDFTResp(mf) + mf_scf.grids_cpks = mf_scf.scf.grids + + nocc, nvir, nmo = mf_scf.nocc, mf_scf.nvir, mf_scf.nmo + so, sv = slice(0, nocc), slice(nocc, nmo) + X = np.random.randn(3, nvir, nocc) + ax_cpks = mf_scf.get_Ax0_cpks()(X) + ax_core = mf_scf.get_Ax0_Core_resp(sv, so, sv, so)(X) + print(np.allclose(ax_cpks, ax_core)) + # print(ax_core[0]) + # print(ax_cpks[0]) + + main_1() + main_2() diff --git a/pyscf/dh/response/hdft/test/test_rhdft.py b/pyscf/dh/response/hdft/test/test_rhdft.py new file mode 100644 index 0000000000000000000000000000000000000000..e1df99e0b8eda1108d50645d96addec32a618ed7 --- /dev/null +++ b/pyscf/dh/response/hdft/test/test_rhdft.py @@ -0,0 +1,58 @@ +import unittest +from pyscf.dh.response.hdft.rhdft import RHDFTResp +from pyscf import gto, dft +import numpy as np + + +class TestRHDFTResp(unittest.TestCase): + + def test_rsh_ax0_cpks(self): + # test that RSH functional is correct for Ax0_cpks + np.random.seed(0) + + mol = gto.Mole(atom="O; H 1 0.94; H 1 0.94 2 104.5", basis="6-31G").build() + mf = dft.RKS(mol, xc="CAM-B3LYP").density_fit().run() + mf_scf = RHDFTResp(mf) + mf_scf.grids_cpks = mf_scf.scf.grids + + nocc, nvir, nmo = mf_scf.nocc, mf_scf.nvir, mf_scf.nmo + so, sv = slice(0, nocc), slice(nocc, nmo) + X = np.random.randn(3, nvir, nocc) + ax_cpks = mf_scf.get_Ax0_cpks()(X) + ax_core = mf_scf.get_Ax0_Core_resp(sv, so, sv, so)(X) + self.assertTrue("eri_cpks_vovo" in mf_scf.tensors) + self.assertTrue(np.allclose(ax_cpks, ax_core)) + + def test_meta_ax0_cpks(self): + # test that meta-GGA functional is correct for Ax0_cpks + np.random.seed(0) + + mol = gto.Mole(atom="O; H 1 0.94; H 1 0.94 2 104.5", basis="6-31G").build() + mf = dft.RKS(mol, xc="TPSS0").density_fit().run() + mf_scf = RHDFTResp(mf) + mf_scf.grids_cpks = mf_scf.scf.grids + + nocc, nvir, nmo = mf_scf.nocc, mf_scf.nvir, mf_scf.nmo + so, sv = slice(0, nocc), slice(nocc, nmo) + X = np.random.randn(3, nvir, nocc) + ax_cpks = mf_scf.get_Ax0_Core(sv, so, sv, so)(X) + ax_core = mf_scf.get_Ax0_Core_resp(sv, so, sv, so)(X) + self.assertTrue("eri_cpks_vovo" in mf_scf.tensors) + self.assertTrue(np.allclose(ax_cpks, ax_core)) + + def test_no_density_fit(self): + # test that RSH functional is correct for Ax0_cpks if no density fitting + np.random.seed(0) + + mol = gto.Mole(atom="O; H 1 0.94; H 1 0.94 2 104.5", basis="6-31G").build() + mf = dft.RKS(mol, xc="CAM-B3LYP").run() + mf_scf = RHDFTResp(mf) + mf_scf.grids_cpks = mf_scf.scf.grids + + nocc, nvir, nmo = mf_scf.nocc, mf_scf.nvir, mf_scf.nmo + so, sv = slice(0, nocc), slice(nocc, nmo) + X = np.random.randn(3, nvir, nocc) + ax_cpks = mf_scf.get_Ax0_Core(sv, so, sv, so)(X) + ax_core = mf_scf.get_Ax0_Core_resp(sv, so, sv, so)(X) + self.assertTrue(np.allclose(ax_cpks, ax_core)) + diff --git a/pyscf/dh/response/mp2/rmp2ri.py b/pyscf/dh/response/mp2/rmp2ri.py new file mode 100644 index 0000000000000000000000000000000000000000..bf33436fd272c6681ba7baff38af3122bcf8ea64 --- /dev/null +++ b/pyscf/dh/response/mp2/rmp2ri.py @@ -0,0 +1,642 @@ +""" RI-MP2 Response-Related Utilities. """ + +from pyscf.dh import RMP2RI +from pyscf.dh import util +from pyscf import scf, lib, __config__ +from pyscf.dh.response import RespBase +from pyscf.dh.response.respbase import get_rdm1_resp_vo_restricted +import h5py +import numpy as np + +CONFIG_incore_t_oovv_mp2 = getattr(__config__, "incore_t_oovv_mp2", "auto") +CONFIG_incore_cderi_uaa_mp2 = getattr(__config__, "incore_cderi_uaa_mp2", "auto") +CONFIG_max_cycle_cpks = getattr(__config__, "max_cycle_cpks", 20) +CONFIG_tol_cpks = getattr(__config__, "tol_cpks", 1e-9) + + +def get_mp2_integrals( + cderi_uov, mo_occ, mo_energy, mask_act, + c_os, c_ss, incore_t_oovv, + verbose=lib.logger.NOTE, max_memory=2000, h5file=NotImplemented): + r""" Get 3-index transformed MP2 amplitude. + + .. math:: + t_{ij}^{ab} &= Y_{ia, P} Y_{jb, P} / D_{ij}^{ab} \\ + T_{ij}^{ab} &= 2 c_\mathrm{os} t_{ij}^{ab} - c_\mathrm{ss} t_{ij}^{ba} \\ + \Gamma_{ia, P} &= T_{ij}^{ab} Y_{jb, P} \\ + + Parameters + ---------- + cderi_uov : np.ndarray + Cholesky decomposed ERI in MO basis (occ-vir block) :math:`Y_{ia, P}`. + dim: (naux, nocc_act, nvir_act). + mo_occ : np.ndarray + Occupation number. + dim: (nmo, ). + mo_energy : np.ndarray + Molecular orbital energy. + dim: (nmo, ). + mask_act : np.ndarray + Active molecular orbital mask. + dim: (nmo, ); dtype: bool. + + c_os : float + Oppo-spin coefficient. + c_ss : float + Same-spin coefficient. + incore_t_oovv : bool or str or float or int + Whether save :math:`t_{ij}^{ab}` into memory (True) or disk (False). + + verbose : int + Print verbosity. + max_memory : int or float + Maximum memory available for molecule object. + h5file : h5py.File + HDF5 file, used when ``t_oovv`` is stored in disk. + + Returns + ------- + dict[str, np.ndarray], dict[str, float] + + Tensors and results are returned. Tensors included: + + - `t_oovv`: MP2 amplitude :math:`t_{ij}^{ab}`. dim: (nocc_act, nocc_act, nvir_act, nvir_act). + - `G_uov`: MP2 amplitude contracted by ERI :math:`\Gamma_{ia, P}`. dim: (naux, nocc_act, nvir_act). + - `rdm1_corr`: Response density matrix of MP2 correlation contribution in MO (non-response) + :math:`D_{pq}^{(2)}`. dim: (nmo, nmo). + + Results are the MP2 energy contributions (including oppo-spin and same-spin). + """ + log = lib.logger.new_logger(verbose=verbose) + time0 = lib.logger.process_clock(), lib.logger.perf_counter() + + # dimension definition and check sanity + mask_occ_act = (mo_occ > 0) & mask_act + mask_vir_act = (mo_occ == 0) & mask_act + nocc_act = mask_occ_act.sum() + nvir_act = mask_vir_act.sum() + nmo = len(mask_act) + naux = cderi_uov.shape[0] + assert len(mo_occ) == nmo + assert (naux, nocc_act, nvir_act) == cderi_uov.shape + assert incore_t_oovv is not None, "t_oovv must be stored for evaluation of MP2 response" + + # preparation + eo = mo_energy[mask_occ_act] + ev = mo_energy[mask_vir_act] + D_ovv = lib.direct_sum("j - a - b -> jab", eo, ev, ev) + + block_oo_act = np.ix_(mask_occ_act, mask_occ_act) + block_vv_act = np.ix_(mask_vir_act, mask_vir_act) + + # allocate results + rdm1_corr = np.zeros((nmo, nmo)) + G_uov = np.zeros((naux, nocc_act, nvir_act)) + t_oovv = util.allocate_array( + incore_t_oovv, (nocc_act, nocc_act, nvir_act, nvir_act), max_memory, + h5file=h5file, name="t_oovv", zero_init=False, chunk=(1, 1, nvir_act, nvir_act)) + eng_bi1 = eng_bi2 = 0 + + # for async write t_oovv + def write_t_oovv(slc, buf): + t_oovv[slc] = buf + + # prepare batch + mem_avail = max_memory - lib.current_memory()[0] + nbatch = util.calc_batch_size(4 * nocc_act * nvir_act**2 + naux * nvir_act, mem_avail) + + with lib.call_in_background(write_t_oovv) as async_write_t_oovv: + for sI in util.gen_batch(0, nocc_act, nbatch): + log.debug(f"[DEBUG] Loop in get_mp2_integrals, slice {sI} of {nocc_act} orbitals.") + D_Oovv = lib.direct_sum("I + jab -> Ijab", eo[sI], D_ovv) + g_Oovv = lib.einsum("PIa, Pjb -> Ijab", cderi_uov[:, sI], cderi_uov) + t_Oovv = g_Oovv / D_Oovv + async_write_t_oovv(sI, t_Oovv) + T_Oovv = util.restricted_biorthogonalize(t_Oovv, 1, c_os, c_ss) + eng_bi1 += lib.einsum("Ijab, Ijab ->", t_Oovv, g_Oovv) + eng_bi2 += lib.einsum("Ijab, Ijba ->", t_Oovv, g_Oovv) + rdm1_corr[block_vv_act] += 2 * lib.einsum("ijac, ijbc -> ab", T_Oovv, t_Oovv) + rdm1_corr[block_oo_act] -= 2 * lib.einsum("ijab, ikab -> jk", T_Oovv, t_Oovv) + G_uov[:, sI] += lib.einsum("ijab, Pjb -> Pia", T_Oovv, cderi_uov) + + # results + eng_os = eng_bi1 + eng_ss = eng_bi1 - eng_bi2 + eng_mp2 = c_os * eng_os + c_ss * eng_ss + results = dict() + results["eng_corr_MP2_bi1"] = eng_bi1 + results["eng_corr_MP2_bi2"] = eng_bi2 + results["eng_corr_MP2_OS"] = eng_os + results["eng_corr_MP2_SS"] = eng_ss + results["eng_corr_MP2"] = eng_mp2 + log.note(f"[RESULT] Energy corr MP2 of same-spin: {eng_ss :18.10f}") + log.note(f"[RESULT] Energy corr MP2 of oppo-spin: {eng_os :18.10f}") + log.note(f"[RESULT] Energy corr MP2 of total: {eng_mp2:18.10f}") + + # tensors + tensors = { + "t_oovv": t_oovv, + "G_uov": G_uov, + "rdm1_corr": rdm1_corr} + + log.timer("get_mp2_integrals", *time0) + return tensors, results + + +def get_W_I(cderi_uov, cderi_uoo, G_uov, verbose=lib.logger.NOTE): + r""" Part I of MO-energy-weighted density matrix. + + .. math:: + W_{ij} [\mathrm{I}] &= -2 \Gamma_{ia, P} Y_{ja, P} \\ + W_{ab} [\mathrm{I}] &= -2 \Gamma_{ia, P} Y_{ib, P} \\ + W_{ai} [\mathrm{I}] &= -4 \Gamma_{ja, P} Y_{ij, P} + + This energy-weighted density matrix is not symmetric by current implementation. + + Parameters + ---------- + cderi_uov : np.ndarray + cderi_uoo : np.ndarray + G_uov : np.ndarray + verbose : int + + Returns + ------- + dict[str, np.ndarray] + """ + + log = lib.logger.new_logger(verbose=verbose) + time0 = lib.logger.process_clock(), lib.logger.perf_counter() + + # dimension definition and check sanity + naux, nocc, nvir = cderi_uov.shape + nmo = nocc + nvir + assert cderi_uoo.shape == (naux, nocc, nocc) + assert G_uov.shape == (naux, nocc, nvir) + + # prepare essential matrices and slices + so, sv = slice(0, nocc), slice(nocc, nmo) + + W_I = np.zeros((nmo, nmo)) + W_I[so, so] = - 2 * lib.einsum("Pia, Pja -> ij", G_uov, cderi_uov) + W_I[sv, sv] = - 2 * lib.einsum("Pia, Pib -> ab", G_uov, cderi_uov) + W_I[sv, so] = - 4 * lib.einsum("Pja, Pij -> ai", G_uov, cderi_uoo) + + log.timer("get_W_I", *time0) + tensors = {"W_I": W_I} + return tensors + + +def get_lag_vo( + G_uov, cderi_uaa, W_I, rdm1_corr, Ax0_Core, + max_memory=2000, verbose=lib.logger.NOTE): + r""" MP2 contribution to Lagrangian vir-occ block. + + .. math:: + L_{ai} = - 4 \Gamma_{ja, P} Y_{ij, P} + 4 \Gamma_{ib, P} Y_{ab, P} + A_{ai, pq} D_{pq}^\mathrm{(2)} + + where :math:`- 4 \Gamma_{ja, P} Y_{ij, P}` is evaluated by :math:`W_{ai} [\mathrm{I}]`. + + Some pratical meaning of lagrangian, is that for vir-occ block, it can be defined by MP2 energy derivative wrt + MO coefficients: + + .. math:: + \mathscr{L}_{pq} &= C_{\mu p} \frac{\partial E^\mathrm{(2)}}{\partial C_{\mu q}} \\ + L_{ai} &= \mathscr{L}_{ai} - \mathscr{L}_{ia} + + Parameters + ---------- + G_uov : np.ndarray + cderi_uaa : np.ndarray + W_I : np.ndarray + rdm1_corr : np.ndarray + Ax0_Core : callable + max_memory : float + verbose : int + + Returns + ------- + dict[str, np.ndarray] + """ + + log = lib.logger.new_logger(verbose=verbose) + time0 = lib.logger.process_clock(), lib.logger.perf_counter() + + # dimension definition and check sanity + naux, nocc, nvir = G_uov.shape + nmo = nocc + nvir + assert cderi_uaa.shape == (naux, nmo, nmo) + if W_I is not None: + assert W_I.shape == (nmo, nmo) + + # prepare essential matrices and slices + so, sv, sa = slice(0, nocc), slice(nocc, nmo), slice(0, nmo) + + # generate lagrangian occ-vir block + lag_vo = np.zeros((nvir, nocc)) + lag_vo += W_I[sv, so] + lag_vo += Ax0_Core(sv, so, sa, sa)(rdm1_corr) + mem_avail = max_memory - lib.current_memory()[0] + nbatch = util.calc_batch_size(nvir ** 2 + nocc * nvir, mem_avail) + for saux in util.gen_batch(0, naux, nbatch): + lag_vo += 4 * lib.einsum("Pib, Pab -> ai", G_uov[saux], cderi_uaa[saux, sv, sv]) + + log.timer("get_lag_vo", *time0) + tensors = {"lag_vo": lag_vo} + return tensors + + +def get_rdm1_corr_resp( + rdm1_corr, lag_vo, + mo_energy, mo_occ, Ax0_Core, + max_cycle=20, tol=1e-9, verbose=lib.logger.NOTE): + r""" 1-RDM of response MP2 correlation by MO basis. + + For other parts, response density matrix is the same to the usual 1-RDM. + + .. math:: + A'_{ai, bj} D_{bj}^\mathrm{(2)} &= L_{ai} + + Parameters + ---------- + rdm1_corr : np.ndarray + lag_vo : np.ndarray + mo_energy : np.ndarray + mo_occ : np.ndarray + Ax0_Core : callable + max_cycle : int + tol : float + verbose : int + + Returns + ------- + dict[str, np.ndarray] + + See Also + -------- + get_rdm1_corr + """ + + log = lib.logger.new_logger(verbose=verbose) + time0 = lib.logger.process_clock(), lib.logger.perf_counter() + + # dimension definition and check sanity + nocc = (mo_occ > 0).sum() + nvir = (mo_occ == 0).sum() + nmo = nocc + nvir + assert rdm1_corr.shape == (nmo, 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 = get_rdm1_resp_vo_restricted( + lag_vo, mo_energy, mo_occ, Ax0_Core, + max_cycle=max_cycle, tol=tol, verbose=verbose) + rdm1_corr_resp[sv, so] = rdm1_corr_resp_vo + + log.timer("get_rdm1_corr_resp", *time0) + tensors = {"rdm1_corr_resp": rdm1_corr_resp} + return tensors + + +class RMP2RespRI(RMP2RI, RespBase): + + def __init__(self, *args, **kwargs): + super().__init__(*args, **kwargs) + self.c_os = kwargs.get("c_os", 1) + self.c_ss = kwargs.get("c_ss", 1) + self.incore_cderi_uaa = CONFIG_incore_cderi_uaa_mp2 + self.incore_t_oovv = kwargs.get("incore_t_oovv", CONFIG_incore_t_oovv_mp2) + self.max_cycle_cpks = CONFIG_max_cycle_cpks + self.tol_cpks = CONFIG_tol_cpks + self._Ax0_Core = NotImplemented + + @property + def mask_occ(self): + return self.mo_occ != 0 + + @property + def mask_occ_act(self): + return self.get_frozen_mask() & (self.mo_occ != 0) + + @property + def mask_vir(self): + return self.mo_occ == 0 + + @property + def mask_vir_act(self): + return self.get_frozen_mask() & (self.mo_occ == 0) + + def make_cderi_uaa(self): + """ Generate cholesky decomposed ERI (all block, full orbitals, s1 symm, in memory/disk). """ + if "cderi_uaa" in self.tensors: + return self.tensors["cderi_uaa"] + + log = lib.logger.new_logger(verbose=self.verbose) + + # dimension and mask + mo_coeff = self.mo_coeff + nmo = self.nmo + + # 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, nmo, nmo), mem_avail, + h5file=self._tmpfile, name="cderi_uaa", zero_init=False) + log.info(f"[INFO] Store type of cderi_uaa: {type(cderi_uaa)}") + + # generate array + util.get_cderi_mo(with_df, mo_coeff, cderi_uaa, max_memory=self.max_memory) + + tensors = {"cderi_uaa": cderi_uaa} + self.tensors.update(tensors) + return cderi_uaa + + def make_cderi_uov(self): + """ Generate cholesky decomposed ERI (occ-vir block, active only, in memory). """ + if "cderi_uov" in self.tensors: + return self.tensors["cderi_uov"] + + # dimension and mask + mask_occ_act = self.mask_occ_act + mask_vir_act = self.mask_vir_act + + # In response evaluation, ERI of all orbitals are generally required. + # Thus, cderi_uaa is a must-evaluated tensor. + # note that cderi_uaa may be stored by h5py, so indexing by list should be done in following way + cderi_uaa = self.tensors.get("cderi_uaa", self.make_cderi_uaa()) + cderi_uov = cderi_uaa[:, mask_occ_act, :][:, :, mask_vir_act] + + tensors = {"cderi_uov": cderi_uov} + self.tensors.update(tensors) + return cderi_uov + + def make_cderi_uoo(self): + """ Generate cholesky decomposed ERI (occ-occ block, in memory). """ + # todo: frozen core not applied + if "cderi_uoo" in self.tensors: + return self.tensors["cderi_uoo"] + + # dimension and mask + mask_occ_act = self.mask_occ_act + + # In response evaluation, ERI of all orbitals are generally required. + # Thus, cderi_uaa is a must-evaluated tensor. + # note that cderi_uaa may be stored by h5py, so indexing by list should be done in following way + cderi_uaa = self.tensors.get("cderi_uaa", self.make_cderi_uaa()) + cderi_uoo = cderi_uaa[:, mask_occ_act, :][:, :, mask_occ_act] + + tensors = {"cderi_uoo": cderi_uoo} + self.tensors.update(tensors) + return cderi_uoo + + def _make_mp2_integrals(self): + cderi_uov = self.tensors.get("cderi_uov", self.make_cderi_uov()) + mo_occ = self.mo_occ + mo_energy = self.mo_energy + mask_act = self.get_frozen_mask() + c_os = self.c_os + c_ss = self.c_ss + incore_t_oovv = self.incore_t_oovv + verbose = self.verbose + max_memory = self.max_memory + h5file = self._tmpfile + tensors, results = get_mp2_integrals( + cderi_uov, mo_occ, mo_energy, mask_act, + c_os, c_ss, incore_t_oovv, + verbose=verbose, max_memory=max_memory, h5file=h5file) + self.tensors.update(tensors) + results = {util.pad_omega(key, self.omega): val for (key, val) in results.items()} + util.update_results(self.results, results) + + def make_t_oovv(self): + r""" MP2 amplitude :math:`t_{ij}^{ab}`. + + dim: (nocc_act, nocc_act, nvir_act, nvir_act) + """ + if "t_oovv" in self.tensors: + return self.tensors["t_oovv"] + + self._make_mp2_integrals() + return self.tensors["t_oovv"] + + def make_G_uov(self): + r""" MP2 amplitude contracted by ERI :math:`\Gamma_{ia, P}`. + + dim: (naux, nocc_act, nvir_act) + """ + if "G_uov" in self.tensors: + return self.tensors["G_uov"] + + self._make_mp2_integrals() + return self.tensors["G_uov"] + + def make_rdm1_corr(self): + r""" Response density matrix of MP2 correlation contribution in MO (non-response) :math:`D_{pq}^{(2)}`. + + dim: (nmo, nmo) + """ + if "rdm1_corr" in self.tensors: + return self.tensors["rdm1_corr"] + + self._make_mp2_integrals() + return self.tensors["rdm1_corr"] + + def driver_eng_mp2(self, **kwargs): + """ Driver of MP2 energy. + + For response computation of MP2, energy evaluation is accompanied with evaluation of MP2 integrals. + """ + if "eng_corr_MP2" not in self.results: + self._make_mp2_integrals() + + return self.results + + def make_W_I(self): + r""" Part I of MO-energy-weighted density matrix W_{pq} [\mathrm{I}]. + + dim: (nmo, nmo) + """ + if "W_I" in self.tensors: + return self.tensors["W_I"] + + # prepare tensors + cderi_uov = self.tensors.get("cderi_uov", self.make_cderi_uov()) + cderi_uoo = self.tensors.get("cderi_uoo", self.make_cderi_uoo()) + G_uov = self.tensors.get("G_uov", self.make_G_uov()) + + # main function + tensors = get_W_I(cderi_uov, cderi_uoo, G_uov, verbose=self.verbose) + self.tensors.update(tensors) + return tensors["W_I"] + + def make_lag_vo(self): + r""" Generate MP2 contribution to Lagrangian vir-occ block :math:`L_{ai}`. """ + if "lag_vo" in self.tensors: + return self.tensors["lag_vo"] + + # prepare tensors + W_I = self.tensors.get("W_I", self.make_W_I()) + cderi_uaa = self.tensors.get("cderi_uaa", self.make_cderi_uaa()) + G_uov = self.tensors["G_uov"] + rdm1_corr = self.tensors["rdm1_corr"] + Ax0_Core = self.Ax0_Core + + # main function + tensors = get_lag_vo( + G_uov, cderi_uaa, W_I, rdm1_corr, Ax0_Core, + max_memory=self.max_memory, verbose=self.verbose) + + self.tensors.update(tensors) + return tensors["lag_vo"] + + def make_rdm1_corr_resp(self): + """ Generate 1-RDM (response) of MP2 correlation :math:`D_{pq}^{(2)}`. """ + if "rdm1_corr_resp" in self.tensors: + return self.tensors["rdm1_corr_resp"] + + # prepare input + lag_vo = self.tensors.get("lag_vo", self.make_lag_vo()) + rdm1_corr = self.tensors["rdm1_corr"] + mask = self.get_frozen_mask() + mo_energy = self.mo_energy[mask] + mo_occ = self.mo_occ[mask] + Ax0_Core = self.Ax0_Core + max_cycle = self.max_cycle_cpks + tol = self.tol_cpks + verbose = self.verbose + + # main function + tensors = get_rdm1_corr_resp( + rdm1_corr, lag_vo, mo_energy, mo_occ, Ax0_Core, + max_cycle=max_cycle, tol=tol, verbose=verbose) + + self.tensors.update(tensors) + return tensors["rdm1_corr_resp"] + + def make_rdm1(self, ao=False): + r""" Generate 1-RDM (non-response) of MP2 :math:`D_{pq}^{MP2}` in MO or :math:`D_{\mu \nu}^{MP2}` in AO. """ + # prepare input + rdm1_corr = self.tensors.get("rdm1_corr", self.make_rdm1_corr()) + + rdm1 = np.diag(self.mo_occ) + mask = self.get_frozen_mask() + ix_act = np.ix_(mask, mask) + rdm1[ix_act] += rdm1_corr + self.tensors["rdm1"] = rdm1 + if ao: + rdm1 = self.mo_coeff @ rdm1 @ self.mo_coeff.T + return rdm1 + + def make_rdm1_resp(self, ao=False): + r""" Generate 1-RDM (response) of MP2 :math:`D_{pq}^{MP2}` in MO or :math:`D_{\mu \nu}^{MP2}` in AO. """ + # prepare input + rdm1_corr_resp = self.tensors.get("rdm1_corr_resp", self.make_rdm1_corr_resp()) + + rdm1 = np.diag(self.mo_occ) + mask = self.get_frozen_mask() + ix_act = np.ix_(mask, mask) + rdm1[ix_act] += rdm1_corr_resp + self.tensors["rdm1_resp"] = rdm1 + if ao: + rdm1 = self.mo_coeff @ rdm1 @ self.mo_coeff.T + return rdm1 + + def make_dipole(self): + # prepare input + mol = self.mol + rdm1_ao = self.make_rdm1_resp(ao=True) + int1e_r = mol.intor("int1e_r") + + dip_elec = - np.einsum("uv, tuv -> t", rdm1_ao, int1e_r) + dip_nuc = np.einsum("A, At -> t", mol.atom_charges(), mol.atom_coords()) + dip = dip_elec + dip_nuc + self.tensors["dipole"] = dip + return dip + + kernel = driver_eng_mp2 + get_mp2_integrals = staticmethod(get_mp2_integrals) + + +if __name__ == '__main__': + def main_1(): + from pyscf import gto, scf + mol = gto.Mole(atom="O; H 1 0.94; H 1 0.94 2 104.5", basis="6-31G").build() + mf = scf.RHF(mol).run() + mf_mp2 = RMP2RespRI(mf) + print(mf_mp2.make_dipole()) + + def main_2(): + # test numerical dipole + from pyscf import gto, scf, mp, df + np.set_printoptions(8, suppress=True, linewidth=150) + mol = gto.Mole(atom="O; H 1 0.94; H 1 0.94 2 104.5", basis="6-31G").build() + + def eng_with_dipole_field(t, h): + mf_scf = scf.RHF(mol).density_fit("cc-pVDZ-jkfit") + mf_scf.get_hcore = lambda *args, **kwargs: scf.hf.get_hcore(mol) - h * mol.intor("int1e_r")[t] + mf_scf.run(conv_tol=1e-12) + mf_mp = RMP2RI(mf_scf).run() + # mf_mp = mp.dfmp2.DFMP2(mf_scf) + # mf_mp.with_df = df.DF(mol, df.make_auxbasis(mol, mp2fit=True)) + # mf_mp.run() + return mf_mp.e_tot + + mf_scf = scf.RHF(mol).density_fit("cc-pVDZ-jkfit").run() + mf_mp = RMP2RespRI(mf_scf).run() + dip_anal = mf_mp.make_dipole() + + eng_array = np.zeros((2, 3)) + h = 1e-4 + for idx, h in [(0, h), [1, -h]]: + for t in (0, 1, 2): + eng_array[idx, t] = eng_with_dipole_field(t, h) + dip_elec_num = - (eng_array[0] - eng_array[1]) / (2 * h) + dip_nuc = np.einsum("A, At -> t", mol.atom_charges(), mol.atom_coords()) + dip_num = dip_elec_num + dip_nuc + + print(mf_scf.dip_moment(unit="AU")) + print(dip_num) + print(dip_anal) + + def main_3(): + # test numerical dipole for DFT PT2 + from pyscf import gto, scf, dft + np.set_printoptions(8, suppress=True, linewidth=150) + mol = gto.Mole(atom="O; H 1 0.94; H 1 0.94 2 104.5", basis="6-31G").build() + + c_os, c_ss = 1.3, 0.6 + + def eng_with_dipole_field(t, h): + mf_scf = dft.RKS(mol, xc="B3LYPg").density_fit("cc-pVDZ-jkfit") + mf_scf.get_hcore = lambda *args, **kwargs: scf.hf.get_hcore(mol) - h * mol.intor("int1e_r")[t] + mf_scf.run(conv_tol=1e-12) + mf_mp = RMP2RI(mf_scf).run() + return mf_mp.scf.e_tot + c_os * mf_mp.results["eng_corr_MP2_OS"] + c_ss * mf_mp.results["eng_corr_MP2_SS"] + + mf_scf = dft.RKS(mol, xc="B3LYPg").density_fit("cc-pVDZ-jkfit").run() + mf_mp = RMP2RespRI(mf_scf).run(c_os=c_os, c_ss=c_ss) + dip_anal = mf_mp.make_dipole() + + eng_array = np.zeros((2, 3)) + h = 1e-4 + for idx, h in [(0, h), [1, -h]]: + for t in (0, 1, 2): + eng_array[idx, t] = eng_with_dipole_field(t, h) + dip_elec_num = - (eng_array[0] - eng_array[1]) / (2 * h) + dip_nuc = np.einsum("A, At -> t", mol.atom_charges(), mol.atom_coords()) + dip_num = dip_elec_num + dip_nuc + + print(mf_scf.dip_moment(unit="AU")) + print(dip_num) + print(dip_anal) + + main_3() diff --git a/pyscf/dh/response/mp2/test/test_rmp2ri.py b/pyscf/dh/response/mp2/test/test_rmp2ri.py new file mode 100644 index 0000000000000000000000000000000000000000..f9fbce0b39fdd97ce11d9fe133f2f9b277858e05 --- /dev/null +++ b/pyscf/dh/response/mp2/test/test_rmp2ri.py @@ -0,0 +1,62 @@ +import unittest +from pyscf.dh.response.mp2.rmp2ri import RMP2RespRI +from pyscf.dh.energy.mp2.rmp2 import RMP2RI +import numpy as np + + +class TestRMP2RI(unittest.TestCase): + def test_num_dipole_ref_hf(self): + from pyscf import gto, scf + mol = gto.Mole(atom="O; H 1 0.94; H 1 0.94 2 104.5", basis="6-31G").build() + + def eng_with_dipole_field(t, h): + mf_scf = scf.RHF(mol).density_fit("cc-pVDZ-jkfit") + mf_scf.get_hcore = lambda *args, **kwargs: scf.hf.get_hcore(mol) - h * mol.intor("int1e_r")[t] + mf_scf.run(conv_tol=1e-12) + mf_mp = RMP2RI(mf_scf).run() + return mf_mp.e_tot + + mf_scf = scf.RHF(mol).density_fit("cc-pVDZ-jkfit").run() + mf_mp = RMP2RespRI(mf_scf).run() + dip_anal = mf_mp.make_dipole() + + eng_array = np.zeros((2, 3)) + h = 1e-4 + for idx, h in [(0, h), [1, -h]]: + for t in (0, 1, 2): + eng_array[idx, t] = eng_with_dipole_field(t, h) + dip_elec_num = - (eng_array[0] - eng_array[1]) / (2 * h) + dip_nuc = np.einsum("A, At -> t", mol.atom_charges(), mol.atom_coords()) + dip_num = dip_elec_num + dip_nuc + + self.assertTrue(np.allclose(dip_num, dip_anal, atol=1e-5, rtol=1e-7)) + + def test_num_dipole_ref_b3lyp(self): + from pyscf import gto, scf, dft + mol = gto.Mole(atom="O; H 1 0.94; H 1 0.94 2 104.5", basis="6-31G").build() + + c_os, c_ss = 1.3, 0.6 + + def eng_with_dipole_field(t, h): + mf_scf = dft.RKS(mol, xc="B3LYPg").density_fit("cc-pVDZ-jkfit") + mf_scf.get_hcore = lambda *args, **kwargs: scf.hf.get_hcore(mol) - h * mol.intor("int1e_r")[t] + mf_scf.run(conv_tol=1e-12) + mf_mp = RMP2RI(mf_scf).run() + return mf_mp.scf.e_tot + c_os * mf_mp.results["eng_corr_MP2_OS"] + c_ss * mf_mp.results["eng_corr_MP2_SS"] + + mf_scf = dft.RKS(mol, xc="B3LYPg").density_fit("cc-pVDZ-jkfit").run() + mf_mp = RMP2RespRI(mf_scf).run(c_os=c_os, c_ss=c_ss) + dip_anal = mf_mp.make_dipole() + + eng_array = np.zeros((2, 3)) + h = 1e-4 + for idx, h in [(0, h), [1, -h]]: + for t in (0, 1, 2): + eng_array[idx, t] = eng_with_dipole_field(t, h) + dip_elec_num = - (eng_array[0] - eng_array[1]) / (2 * h) + dip_nuc = np.einsum("A, At -> t", mol.atom_charges(), mol.atom_coords()) + dip_num = dip_elec_num + dip_nuc + + self.assertTrue(np.allclose(dip_num, dip_anal, atol=1e-5, rtol=1e-7)) + + diff --git a/pyscf/dh/response/respbase.py b/pyscf/dh/response/respbase.py new file mode 100644 index 0000000000000000000000000000000000000000..569286fa91c36e24f4adc916c31ad5c1d8ffa14e --- /dev/null +++ b/pyscf/dh/response/respbase.py @@ -0,0 +1,84 @@ +""" Base class of response instances. """ + +from pyscf.dh.energy import EngBase +from pyscf import scf, __config__, lib +from pyscf.scf import cphf + + +CONFIG_max_cycle_cpks = getattr(__config__, "max_cycle_cpks", 20) +CONFIG_tol_cpks = getattr(__config__, "tol_cpks", 1e-9) + + +class RespBase(EngBase): + def __init__(self, *args, **kwargs): + super().__init__(*args, **kwargs) + + self._Ax0_Core = NotImplemented + self.max_cycle_cpks = CONFIG_max_cycle_cpks + self.tol_cpks = CONFIG_tol_cpks + + @property + def Ax0_Core(self): + """ Fock response of underlying SCF object in MO basis. """ + if self._Ax0_Core is NotImplemented: + restricted = isinstance(self.scf, scf.rhf.RHF) + from pyscf.dh.response.hdft import RHDFTResp + UHDFTResp = NotImplemented + HDFTResp = RHDFTResp if restricted else UHDFTResp + mf_scf = HDFTResp(self.scf) + self._Ax0_Core = mf_scf.get_Ax0_Core + return self._Ax0_Core + + @Ax0_Core.setter + def Ax0_Core(self, Ax0_Core): + self._Ax0_Core = Ax0_Core + + +def get_rdm1_resp_vo_restricted( + lag_vo, mo_energy, mo_occ, Ax0_Core, + max_cycle=20, tol=1e-9, verbose=lib.logger.NOTE): + r""" Solution of response 1-RDM vir-occ part by CP-KS equation (restricted). + + .. math:: + A'_{ai, bj} D_{bj} &= L_{ai} + + Parameters + ---------- + lag_vo : np.ndarray + mo_energy : np.ndarray + mo_occ : np.ndarray + Ax0_Core : callable + max_cycle : int + tol : float + verbose : int + + Returns + ------- + dict[str, np.ndarray] + + See Also + -------- + get_rdm1_corr + """ + + log = lib.logger.new_logger(verbose=verbose) + time0 = lib.logger.process_clock(), lib.logger.perf_counter() + + # dimension definition and check sanity + nocc = (mo_occ > 0).sum() + nvir = (mo_occ == 0).sum() + nmo = nocc + nvir + assert lag_vo.shape == (nvir, nocc) + assert len(mo_energy) == nmo + assert len(mo_occ) == nmo + + # prepare essential matrices and slices + so, sv = slice(0, nocc), slice(nocc, nmo) + + # cp-ks evaluation + rdm1_vo = cphf.solve( + Ax0_Core(sv, so, sv, so), mo_energy, mo_occ, lag_vo, + max_cycle=max_cycle, tol=tol)[0] + + log.timer("get_rdm1_resp_vo_restricted", *time0) + return rdm1_vo diff --git a/pyscf/dh/util/df_addon.py b/pyscf/dh/util/df_addon.py index b05f610b515590966720944b77f7ef674fa98201..f1e7879a9268a6bc3ff7cc438fc5b6355eb13a4d 100644 --- a/pyscf/dh/util/df_addon.py +++ b/pyscf/dh/util/df_addon.py @@ -5,7 +5,9 @@ from pyscf.ao2mo import _ao2mo from pyscf.dh.util import calc_batch_size -def get_cderi_mo(with_df, mo_coeff, Y_mo=None, pqslice=None, max_memory=2000): +def get_cderi_mo( + with_df, mo_coeff, Y_mo=None, pqslice=None, + max_memory=2000, verbose=lib.logger.NOTE): """ Get Cholesky decomposed ERI in MO basis. Parameters @@ -20,12 +22,17 @@ def get_cderi_mo(with_df, mo_coeff, Y_mo=None, pqslice=None, max_memory=2000): Slice of orbitals to be transformed. max_memory : float Maximum memory available. + verbose : int + Print verbosity. Returns ------- np.ndarray Cholesky decomposed ERI in MO basis. """ + log = lib.logger.new_logger(verbose=verbose) + time0 = lib.logger.process_clock(), lib.logger.perf_counter() + naux = with_df.get_naoaux() nmo = mo_coeff.shape[-1] if pqslice is None: @@ -55,6 +62,8 @@ def get_cderi_mo(with_df, mo_coeff, Y_mo=None, pqslice=None, max_memory=2000): Y_ao.reshape(naux, nao, nao), mo_coeff[:, pqslice[0]:pqslice[1]].conj(), mo_coeff[:, pqslice[2]:pqslice[3]]) p0 = p1 + + log.timer("get_mp2_integrals", *time0) return Y_mo