From 82ba87a299919c1df85138c1c4ed5b068a4f5551 Mon Sep 17 00:00:00 2001 From: hebrewsnabla Date: Thu, 2 Feb 2023 22:38:43 +0800 Subject: [PATCH 01/22] include cint.h directory for cmake --- pyscf/lib/CMakeLists.txt | 1 + 1 file changed, 1 insertion(+) diff --git a/pyscf/lib/CMakeLists.txt b/pyscf/lib/CMakeLists.txt index 89158ec9..ce6f20ea 100644 --- a/pyscf/lib/CMakeLists.txt +++ b/pyscf/lib/CMakeLists.txt @@ -87,6 +87,7 @@ if(NOT PYSCF_SOURCE_DIR) endif() message(STATUS "Include pyscf source dir: ${PYSCF_SOURCE_DIR}") include_directories(${PYSCF_SOURCE_DIR}/lib ${PYSCF_SOURCE_DIR}/lib/deps/include) +include_directories(${CINT_DIR}/include) configure_file( "${PYSCF_SOURCE_DIR}/lib/config.h.in" "${PYSCF_SOURCE_DIR}/lib/config.h") From 5e5a4c698e3e28e5bb9273ddfe0b00c5c6f97794 Mon Sep 17 00:00:00 2001 From: JangidBhavnesh Date: Sun, 8 Dec 2024 17:20:40 -0600 Subject: [PATCH 02/22] Translated MGGA activated --- pyscf/mcpdft/otfnal.py | 13 ++++++++++++- pyscf/mcpdft/tfnal_derivs.py | 20 ++++++++++++++++---- 2 files changed, 28 insertions(+), 5 deletions(-) diff --git a/pyscf/mcpdft/otfnal.py b/pyscf/mcpdft/otfnal.py index a9581be6..4a4b3f88 100644 --- a/pyscf/mcpdft/otfnal.py +++ b/pyscf/mcpdft/otfnal.py @@ -253,11 +253,21 @@ def get_ratio (self, Pi, rho_avg): def get_rho_translated (self, Pi, rho, _fn_deriv=0): r''' Compute the "translated" alpha and beta densities: + For the unrestricted case, + rho = [rho^a, rho^b] + Here: + rho^a will have dim of 1,4 or 6 depends on the functional. For MGGA, + rho^a = [rho_u,grad_xu, grad_yu, grad_zu, laplacian_u, tau_u] + Similar for rho_b. + + The translation is done as follows: rho_t^a = (rho/2) * (1 + zeta) rho_t^b = (rho/2) * (1 - zeta) rho'_t^a = (rho'/2) * (1 + zeta) rho'_t^b = (rho'/2) * (1 - zeta) + laplacian_t^a = (laplacian/2) * (1 + zeta) + tau_t^a = (tau/2) * (1 + zeta) See "get_zeta" for the meaning of "zeta" @@ -278,7 +288,8 @@ def get_rho_translated (self, Pi, rho, _fn_deriv=0): Returns: rho_t : ndarray of shape (2,*,ngrids) - Translated spin density (and derivatives) + Translated spin density (and derivatives) in case of LDA or GGAs + Translated spin density, derivatives, and kinetic energy density in case of MGGA ''' # For nonzero charge & pair density, set alpha dens = beta dens diff --git a/pyscf/mcpdft/tfnal_derivs.py b/pyscf/mcpdft/tfnal_derivs.py index 199794e7..4b729692 100644 --- a/pyscf/mcpdft/tfnal_derivs.py +++ b/pyscf/mcpdft/tfnal_derivs.py @@ -106,8 +106,7 @@ def eval_ot(otfnal, rho, Pi, dderiv=1, weights=None, _unpack_vot=True): assert (rho.shape[0] == 2) nderiv = rho.shape[1] nderiv_Pi = Pi.shape[0] - if nderiv > 4: - raise NotImplementedError("Translation of meta-GGA functionals") + rho_t = otfnal.get_rho_translated(Pi, rho) # LDA in libxc has a special numerical problem with zero-valued densities # in one spin @@ -117,8 +116,21 @@ def eval_ot(otfnal, rho, Pi, dderiv=1, weights=None, _unpack_vot=True): idx = (rho_t[0, 0] < 1e-15) & (rho_t[1, 0] > 1e-15) rho_t[0, 0, idx] = 1e-15 rho_tot = rho.sum(0) - rho_deriv = rho_tot[1:4, :] if nderiv > 1 else None - Pi_deriv = Pi[1:4, :] if nderiv_Pi > 1 else None + + if 1 < nderiv <= 4: + rho_deriv = rho_tot[1:4, :] + elif nderiv > 4: + rho_deriv = rho_tot[1:6, :] + else: + rho_deriv = None + + if 1 < nderiv_Pi <= 4: + Pi_deriv = Pi[1:4, :] + elif nderiv > 4: + Pi_deriv = Pi[1:6, :] + else: + Pi_deriv = None + xc_grid = otfnal._numint.eval_xc(otfnal.otxc, (rho_t[0, :, :], rho_t[1, :, :]), spin=1, relativity=0, deriv=dderiv, verbose=otfnal.verbose)[:dderiv + 1] From 4c4914e8b5ccafa3cc654810bac0bb9f8e37f816 Mon Sep 17 00:00:00 2001 From: JangidBhavnesh Date: Sun, 8 Dec 2024 17:23:51 -0600 Subject: [PATCH 03/22] CI Flake Error Fixed --- pyscf/mcpdft/tfnal_derivs.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyscf/mcpdft/tfnal_derivs.py b/pyscf/mcpdft/tfnal_derivs.py index 4b729692..da907213 100644 --- a/pyscf/mcpdft/tfnal_derivs.py +++ b/pyscf/mcpdft/tfnal_derivs.py @@ -120,7 +120,7 @@ def eval_ot(otfnal, rho, Pi, dderiv=1, weights=None, _unpack_vot=True): if 1 < nderiv <= 4: rho_deriv = rho_tot[1:4, :] elif nderiv > 4: - rho_deriv = rho_tot[1:6, :] + rho_deriv = rho_tot[1:6, :] else: rho_deriv = None From 4f360bcb636fafdd2630bd5277e0d33ae08ec43a Mon Sep 17 00:00:00 2001 From: JangidBhavnesh Date: Mon, 9 Dec 2024 11:21:04 -0600 Subject: [PATCH 04/22] Example added for MetaGGAs --- examples/mcpdft/03-metaGGA_functionals.py | 18 ++++++++++++++++++ 1 file changed, 18 insertions(+) create mode 100644 examples/mcpdft/03-metaGGA_functionals.py diff --git a/examples/mcpdft/03-metaGGA_functionals.py b/examples/mcpdft/03-metaGGA_functionals.py new file mode 100644 index 00000000..71f60d43 --- /dev/null +++ b/examples/mcpdft/03-metaGGA_functionals.py @@ -0,0 +1,18 @@ +#!/usr/bin/env/python +from pyscf import gto, scf, mcpdft + +mol = gto.M ( + atom = 'O 0 0 0; O 0 0 1.2', + basis = 'ccpvdz', + spin = 2) + +mf = scf.RHF (mol).run () + +# The translation of Meta-GGAs and hybrid-Meta-GGAs [ChemRxiv. 2024; doi:10.26434/chemrxiv-2024-5hc8g-v2] + +# Translated-Meta-GGA +mc = mcpdft.CASCI(mf, 'tM06L', 6, 8).run () + +# Hybrid-Translated-Meta-GGA +tM06L0 = 't' + mcpdft.hyb('M06L',0.25, hyb_type='average') +mc = mcpdft.CASCI(mf, tM06L0, 6, 8).run () From 11efbc0815aeb5ad9140a1429de939761e94931d Mon Sep 17 00:00:00 2001 From: JangidBhavnesh Date: Thu, 19 Dec 2024 21:51:37 -0600 Subject: [PATCH 05/22] MC23 Functional added --- examples/mcpdft/03-metaGGA_functionals.py | 23 +++++++++ pyscf/mcpdft/otfnal.py | 63 ++++++++++++++++++++++- 2 files changed, 85 insertions(+), 1 deletion(-) diff --git a/examples/mcpdft/03-metaGGA_functionals.py b/examples/mcpdft/03-metaGGA_functionals.py index 71f60d43..5d242925 100644 --- a/examples/mcpdft/03-metaGGA_functionals.py +++ b/examples/mcpdft/03-metaGGA_functionals.py @@ -16,3 +16,26 @@ # Hybrid-Translated-Meta-GGA tM06L0 = 't' + mcpdft.hyb('M06L',0.25, hyb_type='average') mc = mcpdft.CASCI(mf, tM06L0, 6, 8).run () + +# MC23 on-top functional [ChemRxiv. 2024; doi:10.26434/chemrxiv-2024-5hc8g-v2] + +# MC23 = { '0.2952*HF + (1-0.2952)*rep-M06L, 0.2952*HF + (1-0.2952)*rep-M06L'}} +# To calculate the MC23 energies, +# First: Calculate the MCPDFT Energies using repM06L +# Second: Add the CAS Contributation + +# State-Specific +mc = mcpdft.CASCI(mf, 'trepM06L', 6, 8) +e_pdft = mc.kernel()[0] +emc23 = 0.2952*mc.e_mcscf + (1-0.2952)*e_pdft +print("MC-PDFT state %d Energy at %s OT Functional: %d" % (0, "MC23", emc23)) + +# State-average +nroots=2 +mc = mcpdft.CASCI(mf, 'trepM06L', 6, 8) +mc.fcisolver.nroots=nroots +e_pdft = mc.kernel()[0] + +for i in range(nroots): + emc23 = 0.2952*mc.e_mcscf[i] + (1-0.2952)*e_pdft[i] + print("MC-PDFT state %d Energy at %s OT Functional: %d" % (i, "MC23", emc23)) diff --git a/pyscf/mcpdft/otfnal.py b/pyscf/mcpdft/otfnal.py index 4a4b3f88..09067411 100644 --- a/pyscf/mcpdft/otfnal.py +++ b/pyscf/mcpdft/otfnal.py @@ -32,7 +32,58 @@ FT_C = getattr(__config__, 'mcpdft_otfnal_ftransfnal_C', -85.38149682) OT_HYB_ALIAS = {'PBE0' : '0.25*HF + 0.75*PBE, 0.25*HF + 0.75*PBE'} +REG_FUNCTIONALS={} +def register_ot_libxc(xc_base): + ''' + This is a wrapper function to registers the new ot-functionals if they haven't been + registered previously. + Args: + xc_base: str + The name of the ot-functional to be registered. + ''' + xc_base = xc_base.replace("-", "") # In case, the functional name has a hyphen + + def _register_repM06L(): + ''' + This function register the reparametrized + M06-L functional (used for MC23 Functional) in the libxc library. + ''' + from pyscf import dft2 + dft.libxc = dft2.libxc + dft.numint.libxc = dft2.libxc + dft.numint.LibXCMixin.libxc = dft2.libxc + + # Reparametrized-M06L: rep-M06L + # MC23 = { '0.2952*HF + (1-0.2952)*rep-M06L, 0.2952*HF + (1-0.2952)*rep-M06L'}} + # Here, I am registering the rep-M06L functional only. + + # repM06L_C has 27 parameters + MC23_C = np.array([0.06, 0.0031, 0.00515088, 0.00304966, 2.427648, 3.707473, + -7.943377, -2.521466, 2.658691, 2.932276, -8.832841e-01, + -1.895247, -2.899644, -5.068570e-01, -2.712838, 9.416102e-02, + -3.485860e-03, -5.811240e-04, 6.668814e-04, 0.0, 2.669169e-01, + -7.563289e-02, 7.036292e-02, 3.493904e-04, 6.360837e-04, 0.0, 1e-10]) + + # repM06L_X has 18 parameters + MC23_X = np.array([3.352197, 6.332929e-01, -9.469553e-01, 2.030835e-01, + 2.503819, 8.085354e-01, -3.619144, -5.572321e-01, + -4.506606, 9.614774e-01, 6.977048, -1.309337, -2.426371, + -7.896540e-03, 1.364510e-02, -1.714252e-06, -4.698672e-05, 0.0]) + + # See: pyscf/pyscf/dft/libxc_funcs.txt + XC_ID_MGGA_C_M06_L = 233 + XC_ID_MGGA_X_M06_L = 203 + + libxc_register_code = 'repM06L'.lower () + + dft2.libxc.register_custom_functional_(libxc_register_code, 'M06L', + ext_params={XC_ID_MGGA_X_M06_L: MC23_X, + XC_ID_MGGA_C_M06_L: MC23_C}) + + if xc_base.upper() == 'REPM06L' and REG_FUNCTIONALS.get(xc_base.upper()) is None: + _register_repM06L() + def energy_ot (ot, casdm1s, casdm2, mo_coeff, ncore, max_memory=2000, hermi=1): '''Compute the on-top energy - the last term in @@ -747,8 +798,18 @@ def get_transfnal (mol, otxc): 'On-top pair-density functional names other than "translated" (t) or ' '"fully-translated (ft).' ) + xc_base = OT_HYB_ALIAS.get (xc_base.upper (), xc_base) - if ',' not in xc_base and _libxc.is_hybrid_or_rsh (xc_base): + + if xc_base.replace("-","").upper() == 'REPM06L': + # Register the repM06L functional with libxc, if not already done + register_ot_libxc (xc_base) + + # Interface only takes the lower case name for the + # reparametrized functionals. + xc_base = xc_base.lower () + + elif ',' not in xc_base and _libxc.is_hybrid_or_rsh (xc_base): raise NotImplementedError ( 'Aliased or built-in translated hybrid or range-separated ' 'functionals\nother than those listed in otfnal.OT_HYB_ALIAS. ' From e30dcc418ead44efe84d44149977bb4c067ba9af Mon Sep 17 00:00:00 2001 From: JangidBhavnesh Date: Fri, 20 Dec 2024 10:42:41 -0600 Subject: [PATCH 06/22] Flake8 errors fixed --- pyscf/mcpdft/otfnal.py | 21 +++++++++------------ 1 file changed, 9 insertions(+), 12 deletions(-) diff --git a/pyscf/mcpdft/otfnal.py b/pyscf/mcpdft/otfnal.py index 09067411..6683a910 100644 --- a/pyscf/mcpdft/otfnal.py +++ b/pyscf/mcpdft/otfnal.py @@ -46,18 +46,18 @@ def register_ot_libxc(xc_base): def _register_repM06L(): ''' - This function register the reparametrized + This function register the reparametrized M06-L functional (used for MC23 Functional) in the libxc library. ''' from pyscf import dft2 dft.libxc = dft2.libxc dft.numint.libxc = dft2.libxc dft.numint.LibXCMixin.libxc = dft2.libxc - + # Reparametrized-M06L: rep-M06L # MC23 = { '0.2952*HF + (1-0.2952)*rep-M06L, 0.2952*HF + (1-0.2952)*rep-M06L'}} # Here, I am registering the rep-M06L functional only. - + # repM06L_C has 27 parameters MC23_C = np.array([0.06, 0.0031, 0.00515088, 0.00304966, 2.427648, 3.707473, -7.943377, -2.521466, 2.658691, 2.932276, -8.832841e-01, @@ -70,20 +70,20 @@ def _register_repM06L(): 2.503819, 8.085354e-01, -3.619144, -5.572321e-01, -4.506606, 9.614774e-01, 6.977048, -1.309337, -2.426371, -7.896540e-03, 1.364510e-02, -1.714252e-06, -4.698672e-05, 0.0]) - + # See: pyscf/pyscf/dft/libxc_funcs.txt XC_ID_MGGA_C_M06_L = 233 XC_ID_MGGA_X_M06_L = 203 libxc_register_code = 'repM06L'.lower () - dft2.libxc.register_custom_functional_(libxc_register_code, 'M06L', + dft2.libxc.register_custom_functional_(libxc_register_code, 'M06L', ext_params={XC_ID_MGGA_X_M06_L: MC23_X, XC_ID_MGGA_C_M06_L: MC23_C}) - + if xc_base.upper() == 'REPM06L' and REG_FUNCTIONALS.get(xc_base.upper()) is None: _register_repM06L() - + def energy_ot (ot, casdm1s, casdm2, mo_coeff, ncore, max_memory=2000, hermi=1): '''Compute the on-top energy - the last term in @@ -798,16 +798,15 @@ def get_transfnal (mol, otxc): 'On-top pair-density functional names other than "translated" (t) or ' '"fully-translated (ft).' ) - xc_base = OT_HYB_ALIAS.get (xc_base.upper (), xc_base) if xc_base.replace("-","").upper() == 'REPM06L': # Register the repM06L functional with libxc, if not already done register_ot_libxc (xc_base) - + # Interface only takes the lower case name for the # reparametrized functionals. - xc_base = xc_base.lower () + xc_base = xc_base.lower () elif ',' not in xc_base and _libxc.is_hybrid_or_rsh (xc_base): raise NotImplementedError ( @@ -821,8 +820,6 @@ def get_transfnal (mol, otxc): ks.xc = xc_base return fnal_class (ks) - - class colle_salvetti_corr (otfnal): From d4fe6f32345a972628b9e0dc4fa16a5138c5af94 Mon Sep 17 00:00:00 2001 From: JangidBhavnesh Date: Sun, 22 Dec 2024 21:26:12 -0600 Subject: [PATCH 07/22] Few modification and Unit test added --- examples/mcpdft/03-metaGGA_functionals.py | 22 +-- pyscf/mcpdft/_libxc.py | 4 +- pyscf/mcpdft/otfnal.py | 166 +++++++++++-------- pyscf/mcpdft/test/test_mgga.py | 188 ++++++++++++++++++++++ pyscf/mcpdft/tfnal_derivs.py | 12 +- 5 files changed, 298 insertions(+), 94 deletions(-) create mode 100644 pyscf/mcpdft/test/test_mgga.py diff --git a/examples/mcpdft/03-metaGGA_functionals.py b/examples/mcpdft/03-metaGGA_functionals.py index 5d242925..37e5327e 100644 --- a/examples/mcpdft/03-metaGGA_functionals.py +++ b/examples/mcpdft/03-metaGGA_functionals.py @@ -17,25 +17,17 @@ tM06L0 = 't' + mcpdft.hyb('M06L',0.25, hyb_type='average') mc = mcpdft.CASCI(mf, tM06L0, 6, 8).run () -# MC23 on-top functional [ChemRxiv. 2024; doi:10.26434/chemrxiv-2024-5hc8g-v2] +# MC23: meta-hybrid on-top functional [ChemRxiv. 2024; doi:10.26434/chemrxiv-2024-5hc8g-v2] +# To calculate the MC23 energies, you can define functional as 'tMC23'. -# MC23 = { '0.2952*HF + (1-0.2952)*rep-M06L, 0.2952*HF + (1-0.2952)*rep-M06L'}} -# To calculate the MC23 energies, -# First: Calculate the MCPDFT Energies using repM06L -# Second: Add the CAS Contributation +# Note: It's 'fully-translated' version 'ftMC23' is not defined as of now. # State-Specific -mc = mcpdft.CASCI(mf, 'trepM06L', 6, 8) -e_pdft = mc.kernel()[0] -emc23 = 0.2952*mc.e_mcscf + (1-0.2952)*e_pdft -print("MC-PDFT state %d Energy at %s OT Functional: %d" % (0, "MC23", emc23)) +mc = mcpdft.CASCI(mf, 'tMC23', 6, 8) +mc.kernel() # State-average nroots=2 -mc = mcpdft.CASCI(mf, 'trepM06L', 6, 8) +mc = mcpdft.CASCI(mf, 'tMC23', 6, 8) mc.fcisolver.nroots=nroots -e_pdft = mc.kernel()[0] - -for i in range(nroots): - emc23 = 0.2952*mc.e_mcscf[i] + (1-0.2952)*e_pdft[i] - print("MC-PDFT state %d Energy at %s OT Functional: %d" % (i, "MC23", emc23)) +mc.kernel()[0] diff --git a/pyscf/mcpdft/_libxc.py b/pyscf/mcpdft/_libxc.py index 9c35d312..7d951009 100644 --- a/pyscf/mcpdft/_libxc.py +++ b/pyscf/mcpdft/_libxc.py @@ -13,8 +13,8 @@ # See the License for the specific language governing permissions and # limitations under the License. # -from pyscf.dft.libxc import XC_ALIAS, XC_CODES, XC_KEYS -from pyscf.dft.libxc import hybrid_coeff, rsh_coeff +from pyscf.dft2.libxc import XC_ALIAS, XC_CODES, XC_KEYS +from pyscf.dft2.libxc import hybrid_coeff, rsh_coeff from pyscf import lib XC_ALIAS_KEYS = set (XC_ALIAS.keys ()) diff --git a/pyscf/mcpdft/otfnal.py b/pyscf/mcpdft/otfnal.py index 6683a910..a7cb2d6e 100644 --- a/pyscf/mcpdft/otfnal.py +++ b/pyscf/mcpdft/otfnal.py @@ -19,70 +19,92 @@ from scipy import linalg from pyscf import lib, dft from pyscf.lib import logger +from pyscf.dft2 import libxc from pyscf.dft.gen_grid import Grids from pyscf.dft.numint import _NumInt, NumInt from pyscf.mcpdft import pdft_veff, tfnal_derivs, _libxc, _dms, pdft_feff, pdft_eff from pyscf.mcpdft.otpd import get_ontop_pair_density from pyscf import __config__ +# Points to the libxc of the dft2 module +dft.numint.libxc = libxc +dft.numint.LibXCMixin.libxc = libxc + FT_R0 = getattr(__config__, 'mcpdft_otfnal_ftransfnal_R0', 0.9) FT_R1 = getattr(__config__, 'mcpdft_otfnal_ftransfnal_R1', 1.15) FT_A = getattr(__config__, 'mcpdft_otfnal_ftransfnal_A', -475.60656009) FT_B = getattr(__config__, 'mcpdft_otfnal_ftransfnal_B', -379.47331922) FT_C = getattr(__config__, 'mcpdft_otfnal_ftransfnal_C', -85.38149682) -OT_HYB_ALIAS = {'PBE0' : '0.25*HF + 0.75*PBE, 0.25*HF + 0.75*PBE'} -REG_FUNCTIONALS={} - -def register_ot_libxc(xc_base): +OT_HYB_ALIAS = {'PBE0' : '0.25*HF + 0.75*PBE, 0.25*HF + 0.75*PBE', + 'MC23' : 'mc23'} # Note: mc23 is hybrid Meta-GGA OT Functional. + +REG_OT_FUNCTIONALS={} + +# ALIAS for the preset on-top functional +OT_PRESET={ + # Reparametrized-M06L: rep-M06L + # MC23 = { '0.2952*HF + (1-0.2952)*rep-M06L, 0.2952*HF + (1-0.2952)*rep-M06L'}} + # XC_ID_MGGA_C_M06_L = 233 + # XC_ID_MGGA_X_M06_L = 203 + 'MC23':{ + 'xc_code':'M06L', + 'ext_params':{203: np.array([3.352197, 6.332929e-01, -9.469553e-01, 2.030835e-01, + 2.503819, 8.085354e-01, -3.619144, -5.572321e-01, + -4.506606, 9.614774e-01, 6.977048, -1.309337, -2.426371, + -7.896540e-03, 1.364510e-02, -1.714252e-06, -4.698672e-05, 0.0]), + 233: np.array([0.06, 0.0031, 0.00515088, 0.00304966, 2.427648, 3.707473, + -7.943377, -2.521466, 2.658691, 2.932276, -8.832841e-01, + -1.895247, -2.899644, -5.068570e-01, -2.712838, 9.416102e-02, + -3.485860e-03, -5.811240e-04, 6.668814e-04, 0.0, 2.669169e-01, + -7.563289e-02, 7.036292e-02, 3.493904e-04, 6.360837e-04, 0.0, 1e-10])}, + 'hyb':(0.2952,0,0), + 'facs':(0.7048,0.7048)} + } + +def register_otfnal(xc_code, preset): ''' - This is a wrapper function to registers the new ot-functionals if they haven't been + This function registers the new on-top functional if it hasn't been registered previously. Args: - xc_base: str - The name of the ot-functional to be registered. + xc_code: str + The name of the on-top functional to be registered. + preset: dict + The dictionary containing the information about the on-top functional + to be registered. + xc_base: str + The name of the underylying KS-functional in the libxc library. + ext_params: dict, with LibXC exchange and correlation functional integer ID as key, and + an array-like object containing the functional parameters as value. + hyb: tuple + The hybrid functional parameters. + facs: tuple + The mixing factors. + kwargs: dict + The additional keyword arguments. ''' - xc_base = xc_base.replace("-", "") # In case, the functional name has a hyphen - - def _register_repM06L(): - ''' - This function register the reparametrized - M06-L functional (used for MC23 Functional) in the libxc library. - ''' - from pyscf import dft2 - dft.libxc = dft2.libxc - dft.numint.libxc = dft2.libxc - dft.numint.LibXCMixin.libxc = dft2.libxc - - # Reparametrized-M06L: rep-M06L - # MC23 = { '0.2952*HF + (1-0.2952)*rep-M06L, 0.2952*HF + (1-0.2952)*rep-M06L'}} - # Here, I am registering the rep-M06L functional only. - - # repM06L_C has 27 parameters - MC23_C = np.array([0.06, 0.0031, 0.00515088, 0.00304966, 2.427648, 3.707473, - -7.943377, -2.521466, 2.658691, 2.932276, -8.832841e-01, - -1.895247, -2.899644, -5.068570e-01, -2.712838, 9.416102e-02, - -3.485860e-03, -5.811240e-04, 6.668814e-04, 0.0, 2.669169e-01, - -7.563289e-02, 7.036292e-02, 3.493904e-04, 6.360837e-04, 0.0, 1e-10]) - - # repM06L_X has 18 parameters - MC23_X = np.array([3.352197, 6.332929e-01, -9.469553e-01, 2.030835e-01, - 2.503819, 8.085354e-01, -3.619144, -5.572321e-01, - -4.506606, 9.614774e-01, 6.977048, -1.309337, -2.426371, - -7.896540e-03, 1.364510e-02, -1.714252e-06, -4.698672e-05, 0.0]) - - # See: pyscf/pyscf/dft/libxc_funcs.txt - XC_ID_MGGA_C_M06_L = 233 - XC_ID_MGGA_X_M06_L = 203 - - libxc_register_code = 'repM06L'.lower () - - dft2.libxc.register_custom_functional_(libxc_register_code, 'M06L', - ext_params={XC_ID_MGGA_X_M06_L: MC23_X, - XC_ID_MGGA_C_M06_L: MC23_C}) - - if xc_base.upper() == 'REPM06L' and REG_FUNCTIONALS.get(xc_base.upper()) is None: - _register_repM06L() + libxc_register_code = xc_code.lower () + libxc_base_code = preset['xc_base'] + ext_params = preset['ext_params'] + hyb = preset['hyb'] + facs = preset['facs'] + libxc.register_custom_functional_(libxc_register_code, libxc_base_code, + ext_params=ext_params, hyb=hyb, facs=facs) + +def _get_regsitered_ot_functional(xc_code, mol): + ''' + This function returns the on-top functional if it has been registered + previously. + Args: + xc_code: str + The name of the on-top functional to be registered. + ''' + if (xc_code.upper() not in REG_OT_FUNCTIONALS) and (xc_code.upper() in OT_PRESET): + preset = OT_PRESET[xc_code.upper()] + register_otfnal(xc_code, preset) + REG_OT_FUNCTIONALS[xc_code.upper()] = {'hyb_x':preset.get('hyb',0)[0], 'hyb_c':preset.get('hyb',0)[0]} + logger.info(mol, 'Registered the on-top functional: %s', xc_code) + return xc_code.upper() def energy_ot (ot, casdm1s, casdm2, mo_coeff, ncore, max_memory=2000, hermi=1): '''Compute the on-top energy - the last term in @@ -318,7 +340,9 @@ def get_rho_translated (self, Pi, rho, _fn_deriv=0): rho'_t^a = (rho'/2) * (1 + zeta) rho'_t^b = (rho'/2) * (1 - zeta) laplacian_t^a = (laplacian/2) * (1 + zeta) + laplacian_t^b = (laplacian/2) * (1 - zeta) tau_t^a = (tau/2) * (1 + zeta) + tau_t^b = (tau/2) * (1 - zeta) See "get_zeta" for the meaning of "zeta" @@ -798,17 +822,14 @@ def get_transfnal (mol, otxc): 'On-top pair-density functional names other than "translated" (t) or ' '"fully-translated (ft).' ) - xc_base = OT_HYB_ALIAS.get (xc_base.upper (), xc_base) + # Try to register the functional with libxc, if not already done + xc_base = _get_regsitered_ot_functional (xc_base, mol) - if xc_base.replace("-","").upper() == 'REPM06L': - # Register the repM06L functional with libxc, if not already done - register_ot_libxc (xc_base) - - # Interface only takes the lower case name for the - # reparametrized functionals. - xc_base = xc_base.lower () + xc_base = OT_HYB_ALIAS.get (xc_base.upper (), xc_base) - elif ',' not in xc_base and _libxc.is_hybrid_or_rsh (xc_base): + if ',' not in xc_base and \ + (xc_base.upper() not in REG_OT_FUNCTIONALS) and \ + (_libxc.is_hybrid_or_rsh (xc_base)): raise NotImplementedError ( 'Aliased or built-in translated hybrid or range-separated ' 'functionals\nother than those listed in otfnal.OT_HYB_ALIAS. ' @@ -870,19 +891,24 @@ def _hybrid_2c_coeff (ni, xc_code, spin=0): exchange and correlation components of the hybrid coefficent separately ''' - hyb_tot = _NumInt.hybrid_coeff (ni, xc_code, spin=spin) - if hyb_tot == 0: return [0, 0] - - # For exchange-only functionals, hyb_c = hyb_x - x_code, c_code = _libxc.split_x_c_comma (xc_code) - x_code = x_code + ',' - c_code = ',' + c_code - - # All factors of 'HF' are summed by default. Therefore just run the same - # code for the exchange and correlation parts of the string separately - hyb_x = _NumInt.hybrid_coeff(ni, x_code, spin=spin) if len (x_code) else 0 - hyb_c = _NumInt.hybrid_coeff(ni, c_code, spin=spin) if len (c_code) else 0 - return [hyb_x, hyb_c] + if xc_code.upper() in REG_OT_FUNCTIONALS: + hyb_x = REG_OT_FUNCTIONALS[xc_code.upper ()].get('hyb_x', 0) + hyb_c = REG_OT_FUNCTIONALS[xc_code.upper ()].get('hyb_c', 0) + return [hyb_x, hyb_c] + else: + hyb_tot = _NumInt.hybrid_coeff (ni, xc_code, spin=spin) + if hyb_tot == 0: return [0, 0] + + # For exchange-only functionals, hyb_c = hyb_x + x_code, c_code = _libxc.split_x_c_comma (xc_code) + x_code = x_code + ',' + c_code = ',' + c_code + + # All factors of 'HF' are summed by default. Therefore just run the same + # code for the exchange and correlation parts of the string separately + hyb_x = _NumInt.hybrid_coeff(ni, x_code, spin=spin) if len (x_code) else 0 + hyb_c = _NumInt.hybrid_coeff(ni, c_code, spin=spin) if len (c_code) else 0 + return [hyb_x, hyb_c] def make_scaled_fnal (xc_code, hyb_x = 0, hyb_c = 0, fnal_x = None, fnal_c = None): diff --git a/pyscf/mcpdft/test/test_mgga.py b/pyscf/mcpdft/test/test_mgga.py new file mode 100644 index 00000000..9394df95 --- /dev/null +++ b/pyscf/mcpdft/test/test_mgga.py @@ -0,0 +1,188 @@ +#!/usr/bin/env python +# Copyright 2014-2024 The PySCF Developers. All Rights Reserved. +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +# +# Author: Bhavnesh Jangid + +import tempfile +import numpy as np +from pyscf import gto, scf, dft, fci +from pyscf import mcpdft +import unittest + +''' +In this unit-test, test the MCPDFT energies calculated for the LiH +molecule at the state-specific and state-average (2-states) using +1. Meta-GGA functional (M06L) +2. Hybrid-meta-GGA functional M06L0 +3. MC23 Functional +4. Customized Functional + +Test the MCPDFT energies calculated for the triplet water molecule at the +5. Meta-GGA functional (M06L) +6. MC23 Functional + +Note: The reference values from OpenMolcas v22.02, tag 177-gc48a1862b +''' + +def get_lih (r, stateaverage=False, functional='tM06L', basis='sto3g'): + mol = gto.M (atom='Li 0 0 0\nH {} 0 0'.format (r), basis=basis, + output='/dev/null', verbose=0) + mf = scf.RHF (mol).run () + if stateaverage: + mc = mcpdft.CASSCF (mf, functional, 2, 2, grids_level=6) + mc = mc.state_average_([0.5, 0.5]) + else: + mc = mcpdft.CASSCF(mf, functional, 5, 2, grids_level=6) + + mc.fix_spin_(ss=0) + mc = mc.run() + return mc + +def get_water_triplet(functional='tM06L', basis='6-31G'): + mol = gto.M(atom=''' + O 0. 0.000 0.1174 + H 0. 0.757 -0.4696 + H 0. -0.757 -0.4696 + ''',basis=basis, output='/dev/null', verbose=0) + + mf = scf.RHF(mol).run() + + mc = mcpdft.CASSCF(mf, functional, 4, 4, grids_level=6) + solver1 = fci.addons.fix_spin(solver1, ss=2) + solver1.spin = 2 + mc.fcisolver = solver1 + mc.run() + return mc + +def setUpModule(): + global original_grids + global get_lih, lih_tm06l, lih_tmc23, lih_tm06l_sa2, lih_tmc23_sa2 + global get_lih_tm06l0, lih_tm06l0 + global get_lih_custom, lih_custom + global get_water_triplet, water_tm06l, water_tmc23 + + original_grids = dft.radi.ATOM_SPECIFIC_TREUTLER_GRIDS + dft.radi.ATOM_SPECIFIC_TREUTLER_GRIDS = False + lih_tm06l = get_lih(1.5, functional='tM06L') + lih_tmc23 = get_lih(1.5, functional='tMC23') + lih_tm06l_sa2 = get_lih(1.5, stateaverage=True, functional='tM06L') + lih_tmc23_sa2 = get_lih(1.5, stateaverage=True, functional='tMC23') + lih_tm06l0 = get_lih_tm06l0(1.5) + lih_custom = get_lih_custom(1.5) + water_tm06l = get_water_triplet() + water_tmc23 = get_water_triplet(functional='tMC23') + +def tearDownModule(): + global original_grids, lih_tm06l, lih_tmc23, lih_tm06l_sa2, lih_tmc23_sa2 + global lih_tm06l0, get_lih_custom, lih_custom, water_tm06l, water_tmc23 + dft.radi.ATOM_SPECIFIC_TREUTLER_GRIDS = original_grids + lih_tm06l.mol.stdout.close()#Done + lih_tmc23.mol.stdout.close()#Done + lih_tm06l_sa2.mol.stdout.close()#Done + lih_tmc23_sa2.mol.stdout.close()#Done + lih_tm06l0.mol.stdout.close()#Done + lih_custom.mol.stdout.close()#Done + water_tm06l.mol.stdout.close()#Done + water_tmc23.mol.stdout.close()#Done + del original_grids, lih_tm06l, lih_tmc23, lih_tm06l_sa2, lih_tmc23_sa2 + del lih_tm06l0, get_lih_custom, lih_custom, water_tm06l, water_tmc23 + +class KnownValues(unittest.TestCase): + + def assertListAlmostEqual(self, first_list, second_list, expected): + self.assertTrue(len(first_list) == len(second_list)) + for first, second in zip(first_list, second_list): + self.assertAlmostEqual(first, second, expected) + + def test_tmgga(self): + e_mcscf = lih_tm06l.e_mcscf + epdft = lih_tm06l.e_pdft + + sa_e_mcscf = lih_tm06l_sa2.e_mcscf + sa_epdft = lih_tm06l_sa2.e_states + + E_CASSCF_EXPECTED = -7.88112386 + E_MCPDFT_EXPECTED = -7.72800554 + SA_E_CASSCF_EXPECTED = [-7.88112386] + SA_E_MCPDFT_EXPECTED = [-7.88112386] + + self.assertAlmostEqual(e_mcscf, E_CASSCF_EXPECTED, 7) + self.assertAlmostEqual(epdft, E_MCPDFT_EXPECTED, 7) + self.assertListAlmostEqual(sa_e_mcscf, SA_E_CASSCF_EXPECTED, 7) + self.assertListAlmostEqual(sa_epdft, SA_E_MCPDFT_EXPECTED, 7) + + def test_t_hyb_mgga(self): + e_mcscf = lih_tm06l0.e_mcscf + epdft = lih_tm06l0.e_pdft + + E_CASSCF_EXPECTED = -7.88112386 + E_MCPDFT_EXPECTED = -7.72800554 + + self.assertAlmostEqual(e_mcscf, E_CASSCF_EXPECTED, 7) + self.assertAlmostEqual(epdft, E_MCPDFT_EXPECTED, 7) + + def test_tmc23(self): + e_mcscf = lih_tmc23.e_mcscf + epdft = lih_tmc23.e_pdft + + sa_e_mcscf = lih_tmc23_sa2.e_mcscf + sa_epdft = lih_tmc23_sa2.e_states + + E_CASSCF_EXPECTED = -7.88112386 + E_MCPDFT_EXPECTED = -7.72800554 + SA_E_CASSCF_EXPECTED = [-7.88112386] + SA_E_MCPDFT_EXPECTED = [-7.88112386] + + self.assertAlmostEqual(e_mcscf, E_CASSCF_EXPECTED, 7) + self.assertAlmostEqual(epdft, E_MCPDFT_EXPECTED, 7) + self.assertListAlmostEqual(sa_e_mcscf, SA_E_CASSCF_EXPECTED, 7) + self.assertListAlmostEqual(sa_epdft, SA_E_MCPDFT_EXPECTED, 7) + + + def test_water_triplet_tm06l(self): + e_mcscf = water_tm06l.e_mcscf + epdft = water_tm06l.e_pdft + + E_CASSCF_EXPECTED = -7.88112386 + E_MCPDFT_EXPECTED = -7.72800554 + + self.assertAlmostEqual(e_mcscf, E_CASSCF_EXPECTED, 7) + self.assertAlmostEqual(epdft, E_MCPDFT_EXPECTED, 7) + + def test_water_triplet_tmc23(self): + e_mcscf = water_tmc23.e_mcscf + epdft = water_tmc23.e_pdft + + E_CASSCF_EXPECTED = -7.88112386 + E_MCPDFT_EXPECTED = -7.72800554 + + self.assertAlmostEqual(e_mcscf, E_CASSCF_EXPECTED, 7) + self.assertAlmostEqual(epdft, E_MCPDFT_EXPECTED, 7) + + def test_custom_ot_functional(self): + e_mcscf = lih_custom.e_mcscf + epdft = lih_custom.e_pdft + + E_CASSCF_EXPECTED = -7.88112386 + E_MCPDFT_EXPECTED = -7.72800554 + + self.assertAlmostEqual(e_mcscf, E_CASSCF_EXPECTED, 7) + self.assertAlmostEqual(epdft, E_MCPDFT_EXPECTED, 7) + + + +if __name__ == "__main__": + print("Full Tests for MGGAs and MC23") + unittest.main() diff --git a/pyscf/mcpdft/tfnal_derivs.py b/pyscf/mcpdft/tfnal_derivs.py index da907213..4b891ffd 100644 --- a/pyscf/mcpdft/tfnal_derivs.py +++ b/pyscf/mcpdft/tfnal_derivs.py @@ -117,19 +117,17 @@ def eval_ot(otfnal, rho, Pi, dderiv=1, weights=None, _unpack_vot=True): rho_t[0, 0, idx] = 1e-15 rho_tot = rho.sum(0) + if nderiv > 4 and dderiv > 0: + raise NotImplementedError("Meta-GGA functional derivatives") + if 1 < nderiv <= 4: rho_deriv = rho_tot[1:4, :] - elif nderiv > 4: + elif 4 < nderiv <= 6: rho_deriv = rho_tot[1:6, :] else: rho_deriv = None - if 1 < nderiv_Pi <= 4: - Pi_deriv = Pi[1:4, :] - elif nderiv > 4: - Pi_deriv = Pi[1:6, :] - else: - Pi_deriv = None + Pi_deriv = Pi[1:4, :] if nderiv_Pi > 1 else None xc_grid = otfnal._numint.eval_xc(otfnal.otxc, (rho_t[0, :, :], rho_t[1, :, :]), spin=1, relativity=0, deriv=dderiv, From 55b9e96df40485834c0a81c94160da7fe6dcf465 Mon Sep 17 00:00:00 2001 From: JangidBhavnesh Date: Mon, 23 Dec 2024 13:32:54 -0600 Subject: [PATCH 08/22] Unit-Tests modified --- pyscf/mcpdft/otfnal.py | 26 ++++- pyscf/mcpdft/test/test_mgga.py | 190 ++++++++++++++++++--------------- 2 files changed, 127 insertions(+), 89 deletions(-) diff --git a/pyscf/mcpdft/otfnal.py b/pyscf/mcpdft/otfnal.py index a7cb2d6e..7da50ce5 100644 --- a/pyscf/mcpdft/otfnal.py +++ b/pyscf/mcpdft/otfnal.py @@ -48,7 +48,7 @@ # XC_ID_MGGA_C_M06_L = 233 # XC_ID_MGGA_X_M06_L = 203 'MC23':{ - 'xc_code':'M06L', + 'xc_base':'M06L', 'ext_params':{203: np.array([3.352197, 6.332929e-01, -9.469553e-01, 2.030835e-01, 2.503819, 8.085354e-01, -3.619144, -5.572321e-01, -4.506606, 9.614774e-01, 6.977048, -1.309337, -2.426371, @@ -86,10 +86,29 @@ def register_otfnal(xc_code, preset): libxc_register_code = xc_code.lower () libxc_base_code = preset['xc_base'] ext_params = preset['ext_params'] - hyb = preset['hyb'] - facs = preset['facs'] + hyb = preset.get('hyb', None) + facs = preset.get('facs', None) libxc.register_custom_functional_(libxc_register_code, libxc_base_code, ext_params=ext_params, hyb=hyb, facs=facs) + REG_OT_FUNCTIONALS[xc_code.upper()] = {'hyb_x':preset.get('hyb',[0])[0], + 'hyb_c':preset.get('hyb',[0])[0]} + +def unregister_otfnal(xc_code): + ''' + This function unregisters the on-top functional if it has been registered + previously. + Args: + xc_code: str + The name of the on-top functional to be unregistered. + ''' + try: + if xc_code.upper() in REG_OT_FUNCTIONALS: + libxc_unregister_code = xc_code.lower() + libxc.unregister_custom_functional_(libxc_unregister_code) + del REG_OT_FUNCTIONALS[xc_code.upper()] + + except Exception as e: + raise RuntimeError(f"Failed to unregister functional '{xc_code}': {e}") from e def _get_regsitered_ot_functional(xc_code, mol): ''' @@ -102,7 +121,6 @@ def _get_regsitered_ot_functional(xc_code, mol): if (xc_code.upper() not in REG_OT_FUNCTIONALS) and (xc_code.upper() in OT_PRESET): preset = OT_PRESET[xc_code.upper()] register_otfnal(xc_code, preset) - REG_OT_FUNCTIONALS[xc_code.upper()] = {'hyb_x':preset.get('hyb',0)[0], 'hyb_c':preset.get('hyb',0)[0]} logger.info(mol, 'Registered the on-top functional: %s', xc_code) return xc_code.upper() diff --git a/pyscf/mcpdft/test/test_mgga.py b/pyscf/mcpdft/test/test_mgga.py index 9394df95..94803e28 100644 --- a/pyscf/mcpdft/test/test_mgga.py +++ b/pyscf/mcpdft/test/test_mgga.py @@ -15,7 +15,6 @@ # # Author: Bhavnesh Jangid -import tempfile import numpy as np from pyscf import gto, scf, dft, fci from pyscf import mcpdft @@ -24,27 +23,70 @@ ''' In this unit-test, test the MCPDFT energies calculated for the LiH molecule at the state-specific and state-average (2-states) using -1. Meta-GGA functional (M06L) -2. Hybrid-meta-GGA functional M06L0 -3. MC23 Functional -4. Customized Functional +1. Meta-GGA functional (tM06L) +2. Hybrid-meta-GGA functional tM06L0 +3. tMC23 Functional Test the MCPDFT energies calculated for the triplet water molecule at the -5. Meta-GGA functional (M06L) -6. MC23 Functional - -Note: The reference values from OpenMolcas v22.02, tag 177-gc48a1862b +4. Meta-GGA functional (M06L) +5. MC23 Functional + +Note: The reference values from OpenMolcas v24.10, tag 682-gf74be507d +The OpenMolcas results were obtained with this grid settings +&SEWARD +Grid Input +RQuad=TA +NR=100 +LMAX=41 +NOPrun +NOSCreening ''' +# To be consistent with OpenMolcas Grid Settings. Grids_att is defined as below +# Source: pyscf/mcpdft/test/test_diatomic_energies.py + +om_ta_alpha = [0.8, 0.9, # H, He + 1.8, 1.4, # Li, Be + 1.3, 1.1, 0.9, 0.9, 0.9, 0.9, # B - Ne + 1.4, 1.3, # Na, Mg + 1.3, 1.2, 1.1, 1.0, 1.0, 1.0, # Al - Ar + 1.5, 1.4, # K, Ca + 1.3, 1.2, 1.2, 1.2, 1.2, 1.2, 1.2, 1.1, 1.1, 1.1, # Sc - Zn + 1.1, 1.0, 0.9, 0.9, 0.9, 0.9] # Ga - Kr + +def om_treutler_ahlrichs(n, chg, *args, **kwargs): + ''' + "Treutler-Ahlrichs" as implemented in OpenMolcas + ''' + r = np.empty(n) + dr = np.empty(n) + alpha = om_ta_alpha[chg-1] + step = 2.0 / (n+1) # = numpy.pi / (n+1) + ln2 = alpha / np.log(2) + for i in range(n): + x = (i+1)*step - 1 # = numpy.cos((i+1)*step) + r [i] = -ln2*(1+x)**.6 * np.log((1-x)/2) + dr[i] = (step #* numpy.sin((i+1)*step) + * ln2*(1+x)**.6 *(-.6/(1+x)*np.log((1-x)/2)+1/(1-x))) + return r[::-1], dr[::-1] + +my_grids = {'atom_grid': (99,590), + 'radi_method': om_treutler_ahlrichs, + 'prune': False, + 'radii_adjust': None} + def get_lih (r, stateaverage=False, functional='tM06L', basis='sto3g'): mol = gto.M (atom='Li 0 0 0\nH {} 0 0'.format (r), basis=basis, output='/dev/null', verbose=0) mf = scf.RHF (mol).run () + if functional == 'tM06L0': + tM06L0 = 't' + mcpdft.hyb('M06L',0.25, hyb_type='average') + mc = mcpdft.CASSCF(mf, tM06L0, 5, 2, grids_attr=my_grids) + else: + mc = mcpdft.CASSCF(mf, functional, 5, 2, grids_attr=my_grids) + if stateaverage: - mc = mcpdft.CASSCF (mf, functional, 2, 2, grids_level=6) mc = mc.state_average_([0.5, 0.5]) - else: - mc = mcpdft.CASSCF(mf, functional, 5, 2, grids_level=6) mc.fix_spin_(ss=0) mc = mc.run() @@ -55,49 +97,40 @@ def get_water_triplet(functional='tM06L', basis='6-31G'): O 0. 0.000 0.1174 H 0. 0.757 -0.4696 H 0. -0.757 -0.4696 - ''',basis=basis, output='/dev/null', verbose=0) + ''',basis=basis, spin=2,output='/dev/null', verbose=0) mf = scf.RHF(mol).run() - mc = mcpdft.CASSCF(mf, functional, 4, 4, grids_level=6) + mc = mcpdft.CASSCF(mf, functional, 2, 2, grids_attr=my_grids) + solver1 = fci.direct_spin1.FCI(mol) solver1 = fci.addons.fix_spin(solver1, ss=2) - solver1.spin = 2 mc.fcisolver = solver1 - mc.run() + mc = mc.run() return mc def setUpModule(): - global original_grids - global get_lih, lih_tm06l, lih_tmc23, lih_tm06l_sa2, lih_tmc23_sa2 - global get_lih_tm06l0, lih_tm06l0 - global get_lih_custom, lih_custom + global get_lih, lih_tm06l, lih_tmc23, lih_tm06l_sa2, lih_tmc23_sa2,lih_tm06l0 global get_water_triplet, water_tm06l, water_tmc23 - - original_grids = dft.radi.ATOM_SPECIFIC_TREUTLER_GRIDS - dft.radi.ATOM_SPECIFIC_TREUTLER_GRIDS = False lih_tm06l = get_lih(1.5, functional='tM06L') lih_tmc23 = get_lih(1.5, functional='tMC23') lih_tm06l_sa2 = get_lih(1.5, stateaverage=True, functional='tM06L') lih_tmc23_sa2 = get_lih(1.5, stateaverage=True, functional='tMC23') - lih_tm06l0 = get_lih_tm06l0(1.5) - lih_custom = get_lih_custom(1.5) + lih_tm06l0 = get_lih(1.5, functional='tM06L0') water_tm06l = get_water_triplet() water_tmc23 = get_water_triplet(functional='tMC23') def tearDownModule(): - global original_grids, lih_tm06l, lih_tmc23, lih_tm06l_sa2, lih_tmc23_sa2 - global lih_tm06l0, get_lih_custom, lih_custom, water_tm06l, water_tmc23 - dft.radi.ATOM_SPECIFIC_TREUTLER_GRIDS = original_grids - lih_tm06l.mol.stdout.close()#Done - lih_tmc23.mol.stdout.close()#Done - lih_tm06l_sa2.mol.stdout.close()#Done - lih_tmc23_sa2.mol.stdout.close()#Done - lih_tm06l0.mol.stdout.close()#Done - lih_custom.mol.stdout.close()#Done - water_tm06l.mol.stdout.close()#Done - water_tmc23.mol.stdout.close()#Done - del original_grids, lih_tm06l, lih_tmc23, lih_tm06l_sa2, lih_tmc23_sa2 - del lih_tm06l0, get_lih_custom, lih_custom, water_tm06l, water_tmc23 + global lih_tm06l, lih_tmc23, lih_tm06l_sa2, lih_tmc23_sa2 + global lih_tm06l0, water_tm06l, water_tmc23 + lih_tm06l.mol.stdout.close() + lih_tmc23.mol.stdout.close() + lih_tm06l_sa2.mol.stdout.close() + lih_tmc23_sa2.mol.stdout.close() + lih_tm06l0.mol.stdout.close() + water_tm06l.mol.stdout.close() + water_tmc23.mol.stdout.close() + del lih_tm06l, lih_tmc23, lih_tm06l_sa2, lih_tmc23_sa2 + del lih_tm06l0, water_tm06l, water_tmc23 class KnownValues(unittest.TestCase): @@ -108,81 +141,68 @@ def assertListAlmostEqual(self, first_list, second_list, expected): def test_tmgga(self): e_mcscf = lih_tm06l.e_mcscf - epdft = lih_tm06l.e_pdft + epdft = lih_tm06l.e_tot sa_e_mcscf = lih_tm06l_sa2.e_mcscf sa_epdft = lih_tm06l_sa2.e_states - E_CASSCF_EXPECTED = -7.88112386 - E_MCPDFT_EXPECTED = -7.72800554 - SA_E_CASSCF_EXPECTED = [-7.88112386] - SA_E_MCPDFT_EXPECTED = [-7.88112386] + E_CASSCF_EXPECTED = -7.88214917 + E_MCPDFT_EXPECTED = -7.95814186 + SA_E_CASSCF_EXPECTED = [-7.88205449, -7.74391704] + SA_E_MCPDFT_EXPECTED = [-7.95807682, -7.79920022] - self.assertAlmostEqual(e_mcscf, E_CASSCF_EXPECTED, 7) - self.assertAlmostEqual(epdft, E_MCPDFT_EXPECTED, 7) - self.assertListAlmostEqual(sa_e_mcscf, SA_E_CASSCF_EXPECTED, 7) - self.assertListAlmostEqual(sa_epdft, SA_E_MCPDFT_EXPECTED, 7) + self.assertAlmostEqual(e_mcscf, E_CASSCF_EXPECTED, 6) + self.assertAlmostEqual(epdft, E_MCPDFT_EXPECTED, 6) + self.assertListAlmostEqual(sa_e_mcscf, SA_E_CASSCF_EXPECTED, 6) + self.assertListAlmostEqual(sa_epdft, SA_E_MCPDFT_EXPECTED, 6) def test_t_hyb_mgga(self): e_mcscf = lih_tm06l0.e_mcscf - epdft = lih_tm06l0.e_pdft + epdft = lih_tm06l0.e_tot - E_CASSCF_EXPECTED = -7.88112386 - E_MCPDFT_EXPECTED = -7.72800554 + E_CASSCF_EXPECTED = -7.88214917 + E_MCPDFT_EXPECTED = -7.93914369 - self.assertAlmostEqual(e_mcscf, E_CASSCF_EXPECTED, 7) - self.assertAlmostEqual(epdft, E_MCPDFT_EXPECTED, 7) + self.assertAlmostEqual(e_mcscf, E_CASSCF_EXPECTED, 6) + self.assertAlmostEqual(epdft, E_MCPDFT_EXPECTED, 6) def test_tmc23(self): e_mcscf = lih_tmc23.e_mcscf - epdft = lih_tmc23.e_pdft + epdft = lih_tmc23.e_tot sa_e_mcscf = lih_tmc23_sa2.e_mcscf sa_epdft = lih_tmc23_sa2.e_states - E_CASSCF_EXPECTED = -7.88112386 - E_MCPDFT_EXPECTED = -7.72800554 - SA_E_CASSCF_EXPECTED = [-7.88112386] - SA_E_MCPDFT_EXPECTED = [-7.88112386] - - self.assertAlmostEqual(e_mcscf, E_CASSCF_EXPECTED, 7) - self.assertAlmostEqual(epdft, E_MCPDFT_EXPECTED, 7) - self.assertListAlmostEqual(sa_e_mcscf, SA_E_CASSCF_EXPECTED, 7) - self.assertListAlmostEqual(sa_epdft, SA_E_MCPDFT_EXPECTED, 7) + E_CASSCF_EXPECTED = -7.88214917 + E_MCPDFT_EXPECTED = -7.95098727 + SA_E_CASSCF_EXPECTED = [-7.88205449, -7.74391704] + SA_E_MCPDFT_EXPECTED = [-7.95093826, -7.80604012] + self.assertAlmostEqual(e_mcscf, E_CASSCF_EXPECTED, 6) + self.assertAlmostEqual(epdft, E_MCPDFT_EXPECTED, 6) + self.assertListAlmostEqual(sa_e_mcscf, SA_E_CASSCF_EXPECTED, 6) + self.assertListAlmostEqual(sa_epdft, SA_E_MCPDFT_EXPECTED, 6) def test_water_triplet_tm06l(self): e_mcscf = water_tm06l.e_mcscf - epdft = water_tm06l.e_pdft + epdft = water_tm06l.e_tot - E_CASSCF_EXPECTED = -7.88112386 - E_MCPDFT_EXPECTED = -7.72800554 + E_CASSCF_EXPECTED = -75.72365496 + E_MCPDFT_EXPECTED = -76.07686505 - self.assertAlmostEqual(e_mcscf, E_CASSCF_EXPECTED, 7) - self.assertAlmostEqual(epdft, E_MCPDFT_EXPECTED, 7) + self.assertAlmostEqual(e_mcscf, E_CASSCF_EXPECTED, 6) + self.assertAlmostEqual(epdft, E_MCPDFT_EXPECTED, 6) def test_water_triplet_tmc23(self): e_mcscf = water_tmc23.e_mcscf - epdft = water_tmc23.e_pdft - - E_CASSCF_EXPECTED = -7.88112386 - E_MCPDFT_EXPECTED = -7.72800554 - - self.assertAlmostEqual(e_mcscf, E_CASSCF_EXPECTED, 7) - self.assertAlmostEqual(epdft, E_MCPDFT_EXPECTED, 7) + epdft = water_tmc23.e_tot - def test_custom_ot_functional(self): - e_mcscf = lih_custom.e_mcscf - epdft = lih_custom.e_pdft - - E_CASSCF_EXPECTED = -7.88112386 - E_MCPDFT_EXPECTED = -7.72800554 + E_CASSCF_EXPECTED = -75.72365496 + E_MCPDFT_EXPECTED = -76.02630019 - self.assertAlmostEqual(e_mcscf, E_CASSCF_EXPECTED, 7) - self.assertAlmostEqual(epdft, E_MCPDFT_EXPECTED, 7) - - + self.assertAlmostEqual(e_mcscf, E_CASSCF_EXPECTED, 6) + self.assertAlmostEqual(epdft, E_MCPDFT_EXPECTED, 6) if __name__ == "__main__": - print("Full Tests for MGGAs and MC23") + print("Full Tests for MGGAs, Hybrid-MGGAs, and MC23") unittest.main() From cc68eee667b293fa6364e24890988a0314ecea24 Mon Sep 17 00:00:00 2001 From: JangidBhavnesh Date: Mon, 23 Dec 2024 13:35:16 -0600 Subject: [PATCH 09/22] Flake8 errors fixed --- pyscf/mcpdft/test/test_mgga.py | 24 ++++++++++++------------ 1 file changed, 12 insertions(+), 12 deletions(-) diff --git a/pyscf/mcpdft/test/test_mgga.py b/pyscf/mcpdft/test/test_mgga.py index 94803e28..3789bafb 100644 --- a/pyscf/mcpdft/test/test_mgga.py +++ b/pyscf/mcpdft/test/test_mgga.py @@ -131,7 +131,7 @@ def tearDownModule(): water_tmc23.mol.stdout.close() del lih_tm06l, lih_tmc23, lih_tm06l_sa2, lih_tmc23_sa2 del lih_tm06l0, water_tm06l, water_tmc23 - + class KnownValues(unittest.TestCase): def assertListAlmostEqual(self, first_list, second_list, expected): @@ -142,7 +142,7 @@ def assertListAlmostEqual(self, first_list, second_list, expected): def test_tmgga(self): e_mcscf = lih_tm06l.e_mcscf epdft = lih_tm06l.e_tot - + sa_e_mcscf = lih_tm06l_sa2.e_mcscf sa_epdft = lih_tm06l_sa2.e_states @@ -155,21 +155,21 @@ def test_tmgga(self): self.assertAlmostEqual(epdft, E_MCPDFT_EXPECTED, 6) self.assertListAlmostEqual(sa_e_mcscf, SA_E_CASSCF_EXPECTED, 6) self.assertListAlmostEqual(sa_epdft, SA_E_MCPDFT_EXPECTED, 6) - + def test_t_hyb_mgga(self): e_mcscf = lih_tm06l0.e_mcscf epdft = lih_tm06l0.e_tot - + E_CASSCF_EXPECTED = -7.88214917 E_MCPDFT_EXPECTED = -7.93914369 - + self.assertAlmostEqual(e_mcscf, E_CASSCF_EXPECTED, 6) self.assertAlmostEqual(epdft, E_MCPDFT_EXPECTED, 6) - + def test_tmc23(self): e_mcscf = lih_tmc23.e_mcscf epdft = lih_tmc23.e_tot - + sa_e_mcscf = lih_tmc23_sa2.e_mcscf sa_epdft = lih_tmc23_sa2.e_states @@ -186,20 +186,20 @@ def test_tmc23(self): def test_water_triplet_tm06l(self): e_mcscf = water_tm06l.e_mcscf epdft = water_tm06l.e_tot - + E_CASSCF_EXPECTED = -75.72365496 E_MCPDFT_EXPECTED = -76.07686505 - + self.assertAlmostEqual(e_mcscf, E_CASSCF_EXPECTED, 6) self.assertAlmostEqual(epdft, E_MCPDFT_EXPECTED, 6) - + def test_water_triplet_tmc23(self): e_mcscf = water_tmc23.e_mcscf epdft = water_tmc23.e_tot - + E_CASSCF_EXPECTED = -75.72365496 E_MCPDFT_EXPECTED = -76.02630019 - + self.assertAlmostEqual(e_mcscf, E_CASSCF_EXPECTED, 6) self.assertAlmostEqual(epdft, E_MCPDFT_EXPECTED, 6) From f662246700ea2ca9c7d4cd7e29568cf553bf281c Mon Sep 17 00:00:00 2001 From: Bhavnesh Jangid Date: Mon, 23 Dec 2024 17:59:42 -0600 Subject: [PATCH 10/22] Sanity check added, and few other fix --- pyscf/mcpdft/otfnal.py | 20 ++++++++++++++------ pyscf/mcpdft/test/test_mgga.py | 4 +++- 2 files changed, 17 insertions(+), 7 deletions(-) diff --git a/pyscf/mcpdft/otfnal.py b/pyscf/mcpdft/otfnal.py index 7da50ce5..98c19d87 100644 --- a/pyscf/mcpdft/otfnal.py +++ b/pyscf/mcpdft/otfnal.py @@ -26,10 +26,6 @@ from pyscf.mcpdft.otpd import get_ontop_pair_density from pyscf import __config__ -# Points to the libxc of the dft2 module -dft.numint.libxc = libxc -dft.numint.LibXCMixin.libxc = libxc - FT_R0 = getattr(__config__, 'mcpdft_otfnal_ftransfnal_R0', 0.9) FT_R1 = getattr(__config__, 'mcpdft_otfnal_ftransfnal_R1', 1.15) FT_A = getattr(__config__, 'mcpdft_otfnal_ftransfnal_A', -475.60656009) @@ -110,7 +106,7 @@ def unregister_otfnal(xc_code): except Exception as e: raise RuntimeError(f"Failed to unregister functional '{xc_code}': {e}") from e -def _get_regsitered_ot_functional(xc_code, mol): +def _get_registered_ot_functional(xc_code, mol): ''' This function returns the on-top functional if it has been registered previously. @@ -296,6 +292,7 @@ def __init__ (self, ks, **kwargs): otfnal.__init__(self, ks.mol, **kwargs) self.otxc = 't' + ks.xc self._numint = copy.copy (ks._numint) + self._numint.libxc = libxc self.grids = copy.copy (ks.grids) self._numint.hybrid_coeff = t_hybrid_coeff.__get__(self._numint) self._numint.nlc_coeff = t_nlc_coeff.__get__(self._numint) @@ -672,6 +669,7 @@ def __init__ (self, ks, **kwargs): self.C=FT_C self.otxc = 'ft' + ks.xc self._numint = copy.copy (ks._numint) + self._numint.libxc = libxc self.grids = copy.copy (ks.grids) self._numint.hybrid_coeff = ft_hybrid_coeff.__get__(self._numint) self._numint.nlc_coeff = ft_nlc_coeff.__get__(self._numint) @@ -827,6 +825,15 @@ def d_jT_op (self, x, rho, Pi, **kwargs): _CS_c_DEFAULT = 0.2533 _CS_d_DEFAULT = 0.349 +def _sanity_check_ftot(xc_code): + ''' + This function will check the functional type and will + raise the warning for fully-translated MGGAs or custom functionals. + ''' + xc_type = libxc.xc_type(xc_code) + if xc_type not in ['LDA', 'GGA']: + msg = f"fully-translated {xc_type} on-top functionals are not defined" + raise NotImplementedError(msg) def get_transfnal (mol, otxc): if otxc.upper ().startswith ('T'): @@ -834,6 +841,7 @@ def get_transfnal (mol, otxc): fnal_class = transfnal elif otxc.upper ().startswith ('FT'): xc_base = otxc[2:] + _sanity_check_ftot(xc_base) fnal_class = ftransfnal else: raise NotImplementedError ( @@ -841,7 +849,7 @@ def get_transfnal (mol, otxc): '"fully-translated (ft).' ) # Try to register the functional with libxc, if not already done - xc_base = _get_regsitered_ot_functional (xc_base, mol) + xc_base = _get_registered_ot_functional (xc_base, mol) xc_base = OT_HYB_ALIAS.get (xc_base.upper (), xc_base) diff --git a/pyscf/mcpdft/test/test_mgga.py b/pyscf/mcpdft/test/test_mgga.py index 3789bafb..1c7b5510 100644 --- a/pyscf/mcpdft/test/test_mgga.py +++ b/pyscf/mcpdft/test/test_mgga.py @@ -31,7 +31,9 @@ 4. Meta-GGA functional (M06L) 5. MC23 Functional -Note: The reference values from OpenMolcas v24.10, tag 682-gf74be507d +Note: The reference values are calculated using same geometries and basis set +from OpenMolcas v24.10, tag 682-gf74be507d. + The OpenMolcas results were obtained with this grid settings &SEWARD Grid Input From dbc0b91b33572aa81c9b42d71b86e82d4a202aec Mon Sep 17 00:00:00 2001 From: JangidBhavnesh Date: Tue, 24 Dec 2024 09:17:41 -0600 Subject: [PATCH 11/22] Few sanity checks added --- pyscf/mcpdft/test/test_mgga.py | 14 +++++++++++++- pyscf/mcpdft/tfnal_derivs.py | 2 ++ 2 files changed, 15 insertions(+), 1 deletion(-) diff --git a/pyscf/mcpdft/test/test_mgga.py b/pyscf/mcpdft/test/test_mgga.py index 3789bafb..f5d59cbf 100644 --- a/pyscf/mcpdft/test/test_mgga.py +++ b/pyscf/mcpdft/test/test_mgga.py @@ -31,7 +31,9 @@ 4. Meta-GGA functional (M06L) 5. MC23 Functional -Note: The reference values from OpenMolcas v24.10, tag 682-gf74be507d +Note: The reference values are generated from +OpenMolcas v24.10, tag 682-gf74be507d + The OpenMolcas results were obtained with this grid settings &SEWARD Grid Input @@ -146,6 +148,8 @@ def test_tmgga(self): sa_e_mcscf = lih_tm06l_sa2.e_mcscf sa_epdft = lih_tm06l_sa2.e_states + # The CAS and MCPDFT reference values are generated using + # OpenMolcas v24.10, tag 682-gf74be507d E_CASSCF_EXPECTED = -7.88214917 E_MCPDFT_EXPECTED = -7.95814186 SA_E_CASSCF_EXPECTED = [-7.88205449, -7.74391704] @@ -160,6 +164,8 @@ def test_t_hyb_mgga(self): e_mcscf = lih_tm06l0.e_mcscf epdft = lih_tm06l0.e_tot + # The CAS and MCPDFT reference values are generated using + # OpenMolcas v24.10, tag 682-gf74be507d E_CASSCF_EXPECTED = -7.88214917 E_MCPDFT_EXPECTED = -7.93914369 @@ -173,6 +179,8 @@ def test_tmc23(self): sa_e_mcscf = lih_tmc23_sa2.e_mcscf sa_epdft = lih_tmc23_sa2.e_states + # The CAS and MCPDFT reference values are generated using + # OpenMolcas v24.10, tag 682-gf74be507d E_CASSCF_EXPECTED = -7.88214917 E_MCPDFT_EXPECTED = -7.95098727 SA_E_CASSCF_EXPECTED = [-7.88205449, -7.74391704] @@ -187,6 +195,8 @@ def test_water_triplet_tm06l(self): e_mcscf = water_tm06l.e_mcscf epdft = water_tm06l.e_tot + # The CAS and MCPDFT reference values are generated using + # OpenMolcas v24.10, tag 682-gf74be507d E_CASSCF_EXPECTED = -75.72365496 E_MCPDFT_EXPECTED = -76.07686505 @@ -197,6 +207,8 @@ def test_water_triplet_tmc23(self): e_mcscf = water_tmc23.e_mcscf epdft = water_tmc23.e_tot + # The CAS and MCPDFT reference values are generated using + # OpenMolcas v24.10, tag 682-gf74be507d E_CASSCF_EXPECTED = -75.72365496 E_MCPDFT_EXPECTED = -76.02630019 diff --git a/pyscf/mcpdft/tfnal_derivs.py b/pyscf/mcpdft/tfnal_derivs.py index 4b891ffd..9383ceab 100644 --- a/pyscf/mcpdft/tfnal_derivs.py +++ b/pyscf/mcpdft/tfnal_derivs.py @@ -104,6 +104,8 @@ def eval_ot(otfnal, rho, Pi, dderiv=1, weights=None, _unpack_vot=True): if rho.ndim == 2: rho = rho[:, None, :] if Pi.ndim == 1: Pi = Pi[None, :] assert (rho.shape[0] == 2) + assert (rho.shape[1] <= 6), "Undefined behavior for this function" + nderiv = rho.shape[1] nderiv_Pi = Pi.shape[0] From bb743d16ecf997d935e32a3740f41991ad57882a Mon Sep 17 00:00:00 2001 From: Qiming Sun Date: Tue, 24 Dec 2024 09:36:12 -0800 Subject: [PATCH 12/22] Add MSDFT method (#77) * Add MSDFT method * Using generalized Slater-Condon rules for HF integrals in NOCI * Add comments * typo * Lint * Fix bug in nuclear energy treatments * Error in sm_t * Update tests * Update tests * Adjust threshold in tests * Adjust tests * Numerical noises in mcpdft tests --- NOTICE | 2 + examples/msdft/01-simple-noci.py | 34 ++ pyscf/mcpdft/test/test_diatomic_energies.py | 2 +- pyscf/msdft/__init__.py | 1 + pyscf/msdft/noci.py | 417 ++++++++++++++++++++ pyscf/msdft/tests/test_noci.py | 126 ++++++ 6 files changed, 581 insertions(+), 1 deletion(-) create mode 100644 examples/msdft/01-simple-noci.py create mode 100644 pyscf/msdft/__init__.py create mode 100644 pyscf/msdft/noci.py create mode 100644 pyscf/msdft/tests/test_noci.py diff --git a/NOTICE b/NOTICE index dda9b3f1..18ac35d4 100644 --- a/NOTICE +++ b/NOTICE @@ -1,6 +1,7 @@ The MC-PDFT module within this package was developed by: (in chronological order of first commit) +Qiming Sun Matthew R Hermes (University of Chicago) Dayou Zhang (University of Minnesota) Aleksandr Lykhin (University of Chicago) @@ -12,3 +13,4 @@ Bhavnesh Jangid Shirong Wang Jiachen Li Jincheng Yu +Peng Bao \ No newline at end of file diff --git a/examples/msdft/01-simple-noci.py b/examples/msdft/01-simple-noci.py new file mode 100644 index 00000000..35e7c899 --- /dev/null +++ b/examples/msdft/01-simple-noci.py @@ -0,0 +1,34 @@ +#!/usr/bin/env/python + +# Author: Peng Bao +# Edited by: Qiming Sun + +from pyscf import gto, msdft + +mol = gto.M(atom=''' +H 1.080977 -2.558832 0.000000 +H -1.080977 2.558832 0.000000 +H 2.103773 -1.017723 0.000000 +H -2.103773 1.017723 0.000000 +H -0.973565 -1.219040 0.000000 +H 0.973565 1.219040 0.000000 +C 0.000000 0.728881 0.000000 +C 0.000000 -0.728881 0.000000 +C 1.117962 -1.474815 0.000000 +C -1.117962 1.474815 0.000000 +''', basis='sto-3g') + +mf = msdft.NOCI(mol) +mf.xc = 'pbe0' + +h = homo = mol.nelec[0] - 1 +l = h + 1 +# Single excitation orbital pair +mf.s = [[h,l],[h-1,l],[h,l+1],[h-1,l+1]] +# Double excitation orbital pair +mf.d = [[h,l]] + +mf.run() +# reference: +#[-153.93158107 -153.8742658 -153.82198958 -153.69666086 -153.59511111 +# -153.53734913 -153.5155775 -153.47367943 -153.40221993 -153.37353437] diff --git a/pyscf/mcpdft/test/test_diatomic_energies.py b/pyscf/mcpdft/test/test_diatomic_energies.py index 11c45942..33dafb05 100644 --- a/pyscf/mcpdft/test/test_diatomic_energies.py +++ b/pyscf/mcpdft/test/test_diatomic_energies.py @@ -126,7 +126,7 @@ def test_h2_cms3ftlda22_631g (self): # commit: bd596f6cabd6da0301f3623af2de6a14082b34b5 for i in range (3): with self.subTest (state=i): - self.assertAlmostEqual (e[i], e_ref[i], 5) + self.assertAlmostEqual (e[i], e_ref[i], 4) def test_h2_cms2ftlda22_631g (self): e = diatomic ('H', 'H', 1.3, 'ftLDA,VWN3', '6-31G', 2, 2, 2) diff --git a/pyscf/msdft/__init__.py b/pyscf/msdft/__init__.py new file mode 100644 index 00000000..0da86f75 --- /dev/null +++ b/pyscf/msdft/__init__.py @@ -0,0 +1 @@ +from .noci import NOCI diff --git a/pyscf/msdft/noci.py b/pyscf/msdft/noci.py new file mode 100644 index 00000000..4df737d1 --- /dev/null +++ b/pyscf/msdft/noci.py @@ -0,0 +1,417 @@ +#!/usr/bin/env python +# +# Copyright 2024 The PySCF Developers. All Rights Reserved. +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +# +# Author: Peng Bao +# Edited by: Qiming Sun + +''' +Multistate Density Functional Theory (MSDFT) + +References: +[1] Block-Localized Excitation for Excimer Complex and Diabatic Coupling + Peng Bao, Christian P. Hettich, Qiang Shi, and Jiali Gao + J. Chem. Theory Comput. 2021, 17, 240-254 +[2] Block-Localized Density Functional Theory (BLDFT), Diabatic Coupling, and + Their Use in Valence Bond Theory for Representing Reactive Potential Energy + Surfaces + Alessandro Cembran, Lingchun Song, Yirong Mo and Jiali Gao + J. Chem. Theory Comput. 2009, 5, 2702-2716 +[3] Beyond Kohn-Sham Approximation: Hybrid Multistate Wave Function and + Density Functional Theory + Jiali Gao, Adam Grofe, Haisheng Ren, Peng Bao + J. Phys. Chem. Lett. 2016, 7, 5143-5149 +[4] Spin-Multiplet Components and Energy Splittings by Multistate Density + Functional Theory + Adam Grofe, Xin Chen, Wenjian Liu, Jiali Gao + J. Phys. Chem. Lett. 2017, 8, 4838-4845 +''' + +import numpy as np +import scipy.linalg +from pyscf import lib +from pyscf import scf +from pyscf import dft +from pyscf.lib import logger +from pyscf.data.nist import HARTREE2EV as au2ev + +__all__ = ['NOCI'] + +def hf_det_ovlp(msks, mfs): + '''Compute the standard interaction between two non-orthogonal + determinants I and J + ''' + log = logger.new_logger(msks) + mol = msks.mol + _mf = mfs[0].copy() + ovlp = _mf.get_ovlp() + neleca, nelecb = _mf.nelec + occ_mos = [] + for mf in mfs: + mo_coeff_a = mf.mo_coeff[0] + mo_coeff_b = mf.mo_coeff[1] + mo_occ_a = mf.mo_occ[0] + mo_occ_b = mf.mo_occ[1] + occ_mo_a = mo_coeff_a[:,mo_occ_a>0] + occ_mo_b = mo_coeff_b[:,mo_occ_b>0] + occ_mos.append([occ_mo_a, occ_mo_b]) + if occ_mo_a.shape[1] != neleca or occ_mo_b.shape[1] != nelecb: + raise RuntimeError('Electron numbers must be equal') + assert mo_coeff_a.dtype == np.float64 + + # can be evaluated using + # * the generalized Slater-Condon rule (J. Chem. Phys. 131, 124113, 2009) or + # * the density matrix limitation (dml) method (J. Am. Chem. SOC. 1990, 112, 4214). + # E = Tr(dm_12, H[dm_12]) where dm_12 = T_I S_{IJ}^{-1} T_J. + #def dml(mo1_a, mo1_b, mo2_a, mo2_b, f): + # '''See also the det_ovlp function in test_noci.py''' + # o_a = mo1_a.conj().T.dot(s).dot(mo2_a) + # o_b = mo1_b.conj().T.dot(s).dot(mo2_b) + # u_a, s_a, vt_a = scipy.linalg.svd(o_a) + # u_b, s_b, vt_b = scipy.linalg.svd(o_b) + # s_a = np.where(abs(s_a) > 1e-11, s_a, 1e-11) + # s_b = np.where(abs(s_b) > 1e-11, s_b, 1e-11) + # x_a = (u_a/s_a).dot(vt_a) + # x_b = (u_b/s_b).dot(vt_b) + # phase = (np.linalg.det(u_a) * np.linalg.det(u_b) * + # np.linalg.det(vt_a) * np.linalg.det(vt_b)) + # det_ovlp = phase * np.prod(s_a)*np.prod(s_b) + # # One-particle asymmetric density matrix. See also pyscf.scf.uhf.make_asym_dm + # dm_a = mo1_a.dot(x_a).dot(mo2_a.conj().T) + # dm_b = mo1_b.dot(x_b).dot(mo2_b.conj().T) + # dm_01 = (dm_a, dm_b) + # return scf.uhf.UHF.energy_elec(mol, dm_01) * det_ovlp + # + # when I and J differ by symmetry and is strictly zero, the density + # matrix limitation method may encounter numerical issues since |S_{IJ}| is + # strictly zero. + # Here, citing the discussions in https://github.com/pyscf/pyscf-forge/pull/77: + # * If there is only one zero singular value due to symmetry, then the symmetry of + # the two-electron integrals will ultimately cancel the crazy (AO-basis) JK + # matrix. (ii'|jj') and (ij'|ji') are only nonzero if j\to j' corresponds to the + # same symmetry element as i\to i', but this would imply that the jth singular + # value must also be zero, and the i==j case is canceled by exchange. So the + # two-electron part of the energy can't contribute, and neither can the + # one-electron part because it will be zero by symmetry. + # * If there are two singular values corresponding to the same symmetry change, + # then the two states do not in fact have different symmetries. The artificial + # 1e-11 floored singular values cancel between the two factors of the density + # matrix and the final multiplication by det_ovlp, and all is well. + # * If there are more than two zero singular values, then neither of the two terms + # of the Hamiltonian can cancel all of the 1e-11 factors in det_ovlp, and the + # whole thing is at most ~1e-11, which is close enough to zero in most cases. + # + # Below, is evaluated using dml when zero singular values in the + # orbital overlap. Otherwise, generalized Slater-Condon rule is applied. + + # Five elements are cached in det_ovlp_cache: + # * Determinants overlap + # * Asymmetric density matrix for evaluating JK + # * The second density matrix for computing HF energy: Tr(dm, hcore+VHF/2) + # * A factor for generalized Slater-Condon integral + # * Number of different spin-orbitals + det_ovlp_cache = {} + svd_threshold = msks.svd_threshold + for i, mf_bra in enumerate(mfs): + mo1_a, mo1_b = occ_mos[i] + for j, mf_ket in enumerate(mfs[:i]): + mo2_a, mo2_b = occ_mos[j] + o_a = mo1_a.conj().T.dot(ovlp).dot(mo2_a) + o_b = mo1_b.conj().T.dot(ovlp).dot(mo2_b) + + u_a, s_a, vt_a = scipy.linalg.svd(o_a) + u_b, s_b, vt_b = scipy.linalg.svd(o_b) + s_a_overlapped = s_a[s_a > svd_threshold] + s_b_overlapped = s_b[s_b > svd_threshold] + differs_a = s_a.size - s_a_overlapped.size + differs_b = s_b.size - s_b_overlapped.size + differs = differs_a + differs_b + log.debug1('states = (%d %d), GSC differs = %d', i, j, differs) + + if differs == 0: + # Evaluate using the density matrix limitation method + det_ovlp = np.linalg.det(o_a) * np.linalg.det(o_b) + # One-particle asymmetric density matrix. See also pyscf.scf.uhf.make_asym_dm + dm_a = mo1_a.dot(np.linalg.solve(o_a.T, mo2_a.conj().T)) + dm_b = mo1_b.dot(np.linalg.solve(o_b.T, mo2_b.conj().T)) + dm_01 = (dm_a, dm_b) + det_ovlp_cache[i, j] = (det_ovlp, dm_01, dm_01, det_ovlp, differs) + + elif differs == 1: # Generalized Slater-Condon rule + det_ovlp = np.linalg.det(o_a) * np.linalg.det(o_b) + c1_a = mo1_a.dot(u_a[:,s_asvd_threshold].dot(vt_a[s_a>svd_threshold]) + x_b = u_b[:,s_b>svd_threshold].dot(vt_b[s_b>svd_threshold]) + dm_a = mo1_a.dot(x_a).dot(mo2_a.conj().T) + dm_b = mo1_b.dot(x_b).dot(mo2_b.conj().T) + dm_r = (dm_a, dm_b) + + phase = (np.linalg.det(u_a) * np.linalg.det(u_b) * + np.linalg.det(vt_a) * np.linalg.det(vt_b)) + fac = phase * np.prod(s_a_overlapped)*np.prod(s_b_overlapped) + + det_ovlp_cache[i, j] = (det_ovlp, dm_r, dm_01, fac, differs) + + elif differs == 2: # Generalized Slater-Condon rule + det_ovlp = np.linalg.det(o_a) * np.linalg.det(o_b) + c1_a = mo1_a.dot(u_a[:,s_a', hcore, dm_01[0]).real + e1 += np.einsum('ij,ji->', hcore, dm_01[1]).real + e2 = np.einsum('ij,ji->', vhf_01[0], dm_01[0]).real + e2 += np.einsum('ij,ji->', vhf_01[1], dm_01[1]).real + e2 *= .5 + if differs == 2: + # When differed by 2 spin-orbitals, only the two-electron part + # contributes to Slater-Condon integrals + h[idx] = e2 * fac + else: + h[idx] = (e1 + e2) * fac + + h = lib.hermi_triu(h, inplace=True) + s = lib.hermi_triu(s, inplace=True) + return h, s + +def multi_states_scf(msks, ground_ks=None): + '''Construct multiple Kohn-Sham instances for states specified in MSDFT''' + log = logger.new_logger(msks) + if ground_ks is None: + ground_ks = dft.UKS(msks.mol, xc=msks.xc).run() + else: + assert isinstance(ground_ks, dft.uks.UKS) + + neleca, nelecb = ground_ks.nelec + assert neleca == nelecb + + mfs_s = [] + mfs_t = [] + for n, (i, a) in enumerate(msks.s): + log.debug('KS for single excitation %s->%s', i, a) + occ = ground_ks.mo_occ.copy() + occ[0][i] = 0 + occ[0][a] = 1 + mf = ground_ks.copy() + mf = scf.addons.mom_occ(mf, ground_ks.mo_coeff, occ) + dm_init = mf.make_rdm1(ground_ks.mo_coeff, occ) + mf.kernel(dm0=dm_init) + mfs_s.append(mf) + # single excitation for beta electrons + mf = mf.copy() + mf.mo_coeff = mf.mo_coeff[::-1] + mf.mo_occ = mf.mo_occ[::-1] + mf.mo_energy = mf.mo_energy[::-1] + mfs_s.append(mf) + + # spin-flip excitation + log.debug('KS for spin-flip single excitation %s->%s', i, a) + occ = ground_ks.mo_occ.copy() + occ[1][i] = 0 + occ[0][a] = 1 + mf = ground_ks.copy() + mf.nelec = neleca+1, nelecb-1 + mf = scf.addons.mom_occ(mf, ground_ks.mo_coeff, occ) + dm_init = mf.make_rdm1(ground_ks.mo_coeff, occ) + mf.kernel(dm0=dm_init) + mfs_t.append(mf) + + mfs_d = [] + for n, (i, a) in enumerate(msks.d): + log.debug('KS for double excitation (%s,%s)->(%s,%s)', i, i, a, a) + occ = ground_ks.mo_occ.copy() + occ[0][i] = 0 + occ[0][a] = 1 + occ[1][i] = 0 + occ[1][a] = 1 + mf = ground_ks.copy() + mf = scf.addons.mom_occ(mf, ground_ks.mo_coeff, occ) + dm_init = mf.make_rdm1(ground_ks.mo_coeff, occ) + mf.kernel(dm0=dm_init) + mfs_d.append(mf) + + e_g = ground_ks.e_tot + log.info('Ground state KS energy = %g', e_g) + log.info('Doubly excited energy:') + for i, mf in enumerate(mfs_d): + e_d = mf.e_tot + log.info('%-2d %18.15g AU %12.6g eV', i+1, e_d, (e_d-e_g)*au2ev) + + log.info('Single and triple excitation:') + log.info(' E(S) E(T) dEt dEs') + for i, (mf_s, mf_t) in enumerate(zip(mfs_s[::2], mfs_t)): + dEt = (mf_t.e_tot - e_g) * au2ev + e_split = (mf_s.e_tot - mf_t.e_tot) * au2ev + log.info('%-2d %18.15g AU %18.15g AU %15.9g eV %15.9g eV', + i+1, mf_s.e_tot, mf_t.e_tot, dEt, e_split) + return [ground_ks], mfs_s, mfs_d, mfs_t + +class NOCI(lib.StreamObject): + ''' + Nonorthogonal Configuration Interaction (NOCI) of Multistate Density Functional Theory (MSDFT) + + Attributes: + xc : str + Name of exchange-correlation functional + s : + A list of singly excited orbital pairs. Each pair [i, a] means an + excitation from occupied orbital i to a. + d : + A list of doubly excited orbital pairs. Each pair [i, a] means both + alpha and beta electrons at orbital i are excited to orbital a. + coup : int + How to compute the electronic coupling between diabatic states (Bao, JCTC, 17, 240). + * 0: geometric average over diagonal terms. + * 1: determinant-weighted average of correlation. + * 2 (default): overlap-scaled average of correlation. + ci_g : bool + Whether to compute the adiabatic ground-state energy. True by default. + sm_t: bool + Use the energy difference between mix state and Ms=1 triplet state + as the coupling between two symmetry-adapted mix state. This can be + more accurate than the approximate HF coupling. True by default. + + Saved results: + e_tot : float + Total HF energy (electronic energy plus nuclear repulsion) + csfvec : array + CI coefficients + mfs : + KS instances of the underlying diabatic states) + ''' + _keys = { + 'mol', 'verbose', 'stdout', 'xc', 'coup', 'ci_g', 's', 'd', 'sm_t', + 'e_tot', 'csfvec', 'mfs', + } + + coup = 2 + ci_g = True + sm_t = True + svd_threshold = 1e-7 + + def __init__(self, mol, xc=None): + self.mol = mol + self.verbose = mol.verbose + self.stdout = mol.stdout + self.xc = xc + self.s = [] + self.d = [] +################################################## +# don't modify the following attributes, they are not input options + self.e_tot = 0 + self.csfvec = None + self.mfs = None + + def dump_flags(self, verbose=None): + log = logger.new_logger(self, verbose) + if log.verbose < logger.INFO: + return self + + log.info('\n') + log.info('******** %s ********', self.__class__) + log.info('xc = %s', self.xc) + log.info('coup = %s', self.coup) + log.info('ci_g = %s', self.ci_g) + log.info('sm_t = %s', self.sm_t) + log.info('single excitation = %s', self.s) + log.info('double excitation = %s', self.d) + log.info('Overlap svd threshold = %s', self.svd_threshold) + return self + + def kernel(self, ground_ks=None): + log = logger.new_logger(self) + self.dump_flags(log) + self.check_sanity() + + mf_gs, mfs_s, mfs_d, mfs_t = multi_states_scf(self, ground_ks) + mfs = mfs_s + mfs_d + if self.ci_g: + mfs = mfs + mf_gs + self.mfs = mfs + + e_hf, s_csf = hf_det_ovlp(self, mfs) + Enuc = self.mol.energy_nuc() + e_ks = np.array([mf.e_tot for mf in mfs]) + log.debug1('KS energies %s', e_ks) + e_ks -= Enuc + + # Compute transition density functional energy + if self.coup == 0: + # geometric average over diagonal terms. + d = e_ks / e_hf.diagonal() + h_tdf = e_hf * (d[:,None] * d)**.5 + s_csf * Enuc + elif self.coup == 1: + d = e_hf.diagonal() + # determinant-weighted average of correlation. + h_tdf = e_hf * (e_ks[:,None]+e_ks) / (d[:,None]+d) + s_csf * Enuc + elif self.coup == 2: + # overlap-scaled average of correlation. + d = e_ks - e_hf.diagonal() + h_tdf = e_hf + s_csf * ((d[:,None] + d) / 2 + Enuc) + + if self.sm_t: + n_triplets = len(mfs_t) + assert n_triplets * 2 == len(mfs_s) + e_t = np.array([mf.e_tot for mf in mfs_t]) + e_s = e_ks[:n_triplets*2] + Enuc + log.debug1('KS singlet energies %s', e_s) + log.debug1('KS triplet energies %s', e_t) + + for i in range(n_triplets): + j = 2*i + s_t_coupling = e_s[j] + (s_csf[j,j+1] - 1.) * e_t[i] + h_tdf[j,j+1] = h_tdf[j+1,j] = s_t_coupling + self.e_tot, self.csfvec = scipy.linalg.eigh(h_tdf, s_csf) + log.note('MSDFT eigs %s', self.e_tot) + return self.e_tot + + @property + def converged(self): + return all(mf.converged for mf in self.mfs) + + to_gpu = lib.to_gpu diff --git a/pyscf/msdft/tests/test_noci.py b/pyscf/msdft/tests/test_noci.py new file mode 100644 index 00000000..9b25ea96 --- /dev/null +++ b/pyscf/msdft/tests/test_noci.py @@ -0,0 +1,126 @@ +from functools import reduce +import numpy +import numpy as np +import scipy.linalg +from pyscf import gto, scf, lib, fci, ao2mo +from pyscf.msdft import noci + +def test_hf_det_ovlp(): + mol = gto.M(atom=''' +O 0. 0. 0. +H 0. -.757 .587 +H 0. .757 .587 +H 0.5 0.1 -0.2 +''', basis='6-31g', spin=1) + ms_ks = noci.NOCI(mol) + # Reduce iterations to prevent numerical instablity + mf0 = mol.UKS(xc='b3lyp').run(max_cycle=1) + mf1 = mf0.copy() + occ = mf0.mo_occ.copy() + occ[0][mf0.nelec[0]-1] = 0 + occ[0][mf0.nelec[0]+1] = 1 + mf1 = scf.addons.mom_occ(mf1, mf0.mo_coeff, occ).run(max_cycle=1) + h, s = noci.hf_det_ovlp(ms_ks, [mf0, mf1]) + ref = np.array([[-9.35176786e+01, -6.82503177e-02], + [-6.82503177e-02, -9.33368874e+01]]) + assert abs(h/ref - 1.).max() < 1e-7 + +def test_noci_e_tot(): + mol = gto.M(atom=''' +N 0. 0. 0. +H 0. -1.51 1.17 +H 0. 1.51 1.17 +H 1.5 0.1 -0.2 +''', basis='6-31g') + mf = noci.NOCI(mol) + mf.xc = 'pbe0' + mf.s = [[4,5], [4,6]] + mf.d = [[4,5]] + mf.sm_t = False + mf.run() + assert abs(mf.e_tot[0] - -56.161179917474) < 1e-8 + mf.sm_t = True + mf.run() + assert abs(mf.e_tot[0] - -56.161253460503) < 1e-8 + + mol = gto.M(atom=''' +O 0. 0. 0. +H 0. -1.51 1.17 +H 0. 1.51 1.17 +''', basis='6-31g') + mf = noci.NOCI(mol) + mf.xc = 'pbe0' + mf.s = [[4,5]] + mf.d = [[4,5]] + mf.run() + assert abs(mf.e_tot[0] - -76.0190855601) < 1e-7 + +def det_ovlp(mo1, mo2, occ1, occ2, ovlp): + if numpy.sum(occ1) !=numpy.sum(occ2): + raise RuntimeError('Electron numbers are not equal. Electronic coupling does not exist.') + c1_a = mo1[0][:, occ1[0]>0] + c1_b = mo1[1][:, occ1[1]>0] + c2_a = mo2[0][:, occ2[0]>0] + c2_b = mo2[1][:, occ2[1]>0] + o_a = numpy.asarray(reduce(numpy.dot, (c1_a.conj().T, ovlp, c2_a))) + o_b = numpy.asarray(reduce(numpy.dot, (c1_b.conj().T, ovlp, c2_b))) + u_a, s_a, vt_a = scipy.linalg.svd(o_a) + u_b, s_b, vt_b = scipy.linalg.svd(o_b) + s_a = numpy.where(abs(s_a) > 1.0e-11, s_a, 1.0e-11) + s_b = numpy.where(abs(s_b) > 1.0e-11, s_b, 1.0e-11) + OV = numpy.linalg.det(u_a)*numpy.linalg.det(u_b) \ + * numpy.prod(s_a)*numpy.prod(s_b) \ + * numpy.linalg.det(vt_a)*numpy.linalg.det(vt_b) + x_a = reduce(numpy.dot, (u_a*numpy.reciprocal(s_a), vt_a)) + x_b = reduce(numpy.dot, (u_b*numpy.reciprocal(s_b), vt_b)) + return OV, numpy.array((x_a, x_b)) + +def scoup_dml(mol, mo0, mo1, occ0, occ1): + mf = scf.UHF(mol) + # Calculate overlap between two determiant + s, x = det_ovlp(mo0, mo1, occ0, occ1, mf.get_ovlp()) + # Construct density matrix + dm_01 = mf.make_asym_dm(mo0, mo1, occ0, occ1, x) + # One-electron part contrbution + h1e = mf.get_hcore(mol) + # Two-electron part contrbution. D_{IF} is asymmetric + #vhf_01 = get_veff(mf, dm_01, hermi=0) + vj, vk = mf.get_jk(mol, dm_01, hermi=0) + vhf_01 = vj[0] + vj[1] - vk + # New total energy + e_01 = mf.energy_elec(dm_01, h1e, vhf_01) + return s, s * e_01[0], dm_01 + +def test_scoup_vs_fci(): + numpy.random.seed(4) + coords = numpy.random.rand(6, 3) + mol = gto.M(atom=[('H', r) for r in coords], verbose=1) + mf = mol.UHF().run() + nmo = mf.mo_coeff[0].shape[1] + nelec = (3,3) + u = np.linalg.svd(np.random.rand(nmo,nmo))[0] + mo = mf.mo_coeff[0], mf.mo_coeff[1].dot(u) + + eri_aa = ao2mo.kernel(mf._eri, mo[0]) + eri_bb = ao2mo.kernel(mf._eri, mo[1]) + eri_ab = ao2mo.kernel(mf._eri, (mo[0], mo[0], mo[1], mo[1])) + eri = eri_aa, eri_ab, eri_bb + h1e_a = mo[0].T.dot(mf.get_hcore()).dot(mo[0]) + h1e_b = mo[1].T.dot(mf.get_hcore()).dot(mo[1]) + h1e = h1e_a, h1e_b + h2e = fci.direct_uhf.absorb_h1e(h1e, eri, nmo, nelec, .5) + s1e_a = mf.mo_coeff[0].T.dot(mf.get_ovlp()).dot(mo[0]) + s1e_b = mf.mo_coeff[1].T.dot(mf.get_ovlp()).dot(mo[1]) + s = (s1e_a, s1e_b) + + linki = fci.direct_spin1._unpack(nmo, nelec, None) + na = linki[0].shape[0] + nb = linki[1].shape[0] + ket = np.zeros((na, nb)) + ket[0,0] = 1. + bra = ket.copy() + + ci1 = fci.direct_uhf.contract_2e(h2e, ket, nmo, nelec, linki) + ref = fci.addons.overlap(bra, ci1, nmo, nelec, s) + scoup_out = scoup_dml(mol, mf.mo_coeff, mo, mf.mo_occ, mf.mo_occ)[1] + assert abs(ref - scoup_out) < 1e-12 From bbb1f67c0e2ee56f90f8726b7e4b3df16553c4bd Mon Sep 17 00:00:00 2001 From: baopengbp <46704539+baopengbp@users.noreply.github.com> Date: Fri, 3 Jan 2025 10:21:46 +0800 Subject: [PATCH 13/22] =?UTF-8?q?Update=20noci.py=EF=BC=8Cfix=20output=20o?= =?UTF-8?q?f=20dEs=20and=20dSplit?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- pyscf/msdft/noci.py | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/pyscf/msdft/noci.py b/pyscf/msdft/noci.py index 4df737d1..46ca1b62 100644 --- a/pyscf/msdft/noci.py +++ b/pyscf/msdft/noci.py @@ -283,12 +283,13 @@ def multi_states_scf(msks, ground_ks=None): log.info('%-2d %18.15g AU %12.6g eV', i+1, e_d, (e_d-e_g)*au2ev) log.info('Single and triple excitation:') - log.info(' E(S) E(T) dEt dEs') + log.info(' E(S) E(T) dEt dEs dSplit') for i, (mf_s, mf_t) in enumerate(zip(mfs_s[::2], mfs_t)): dEt = (mf_t.e_tot - e_g) * au2ev e_split = (mf_s.e_tot - mf_t.e_tot) * au2ev - log.info('%-2d %18.15g AU %18.15g AU %15.9g eV %15.9g eV', - i+1, mf_s.e_tot, mf_t.e_tot, dEt, e_split) + dEs = dEt + 2 * e_split + log.info('%-2d %18.15g AU %18.15g AU %15.9g eV %15.9g eV %15.9g eV', + i+1, mf_s.e_tot, mf_t.e_tot, dEt, dEs, e_split) return [ground_ks], mfs_s, mfs_d, mfs_t class NOCI(lib.StreamObject): From 18900a5d1c2cda1e10ae3f477d2ee17d656ae8b0 Mon Sep 17 00:00:00 2001 From: hebrewsnabla Date: Fri, 3 Jan 2025 16:15:49 +0800 Subject: [PATCH 14/22] fix unnecessary pi deriv --- pyscf/mcpdft/otfnal.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/pyscf/mcpdft/otfnal.py b/pyscf/mcpdft/otfnal.py index a9581be6..5414b302 100644 --- a/pyscf/mcpdft/otfnal.py +++ b/pyscf/mcpdft/otfnal.py @@ -68,6 +68,7 @@ def energy_ot (ot, casdm1s, casdm2, mo_coeff, ncore, max_memory=2000, hermi=1): ni, xctype = ot._numint, ot.xctype if xctype=='HF': return E_ot dens_deriv = ot.dens_deriv + pi_deriv = 0 nao = mo_coeff.shape[0] ncas = casdm2.shape[0] @@ -84,7 +85,7 @@ def energy_ot (ot, casdm1s, casdm2, mo_coeff, ncore, max_memory=2000, hermi=1): rho = np.asarray ([m[0] (0, ao, mask, xctype) for m in make_rho]) t0 = logger.timer (ot, 'untransformed density', *t0) Pi = get_ontop_pair_density (ot, rho, ao, cascm2, mo_cas, - dens_deriv, mask) + pi_deriv, mask) t0 = logger.timer (ot, 'on-top pair density calculation', *t0) if rho.ndim == 2: rho = np.expand_dims (rho, 1) From 0887d485df80217a022bd69ea7e09a10dc9d7093 Mon Sep 17 00:00:00 2001 From: hebrewsnabla Date: Fri, 3 Jan 2025 16:35:09 +0800 Subject: [PATCH 15/22] add XC_DIR --- pyscf/lib/CMakeLists.txt | 1 + 1 file changed, 1 insertion(+) diff --git a/pyscf/lib/CMakeLists.txt b/pyscf/lib/CMakeLists.txt index 15ffa28d..bb588cc7 100644 --- a/pyscf/lib/CMakeLists.txt +++ b/pyscf/lib/CMakeLists.txt @@ -96,6 +96,7 @@ configure_file( # to find config.h link_directories(${PYSCF_SOURCE_DIR}/lib/deps/lib ${PYSCF_SOURCE_DIR}/lib/deps/lib64) link_directories(${PYSCF_SOURCE_DIR}/lib) +link_directories(${XC_DIR}/lib) message(STATUS "${PYSCF_SOURCE_DIR}/lib may need to be put in the environment LD_LIBRARY_PATH") # See also https://gitlab.kitware.com/cmake/community/wikis/doc/cmake/RPATH-handling From c730061bfc62efdd0527f84db7e8eeb2452bbcd8 Mon Sep 17 00:00:00 2001 From: hebrewsnabla Date: Fri, 3 Jan 2025 17:20:33 +0800 Subject: [PATCH 16/22] fix ftrans --- pyscf/mcpdft/otfnal.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pyscf/mcpdft/otfnal.py b/pyscf/mcpdft/otfnal.py index 5414b302..095bb29f 100644 --- a/pyscf/mcpdft/otfnal.py +++ b/pyscf/mcpdft/otfnal.py @@ -68,7 +68,7 @@ def energy_ot (ot, casdm1s, casdm2, mo_coeff, ncore, max_memory=2000, hermi=1): ni, xctype = ot._numint, ot.xctype if xctype=='HF': return E_ot dens_deriv = ot.dens_deriv - pi_deriv = 0 + Pi_deriv = ot.Pi_deriv nao = mo_coeff.shape[0] ncas = casdm2.shape[0] @@ -85,7 +85,7 @@ def energy_ot (ot, casdm1s, casdm2, mo_coeff, ncore, max_memory=2000, hermi=1): rho = np.asarray ([m[0] (0, ao, mask, xctype) for m in make_rho]) t0 = logger.timer (ot, 'untransformed density', *t0) Pi = get_ontop_pair_density (ot, rho, ao, cascm2, mo_cas, - pi_deriv, mask) + Pi_deriv, mask) t0 = logger.timer (ot, 'on-top pair density calculation', *t0) if rho.ndim == 2: rho = np.expand_dims (rho, 1) From 520e057aa08130705fb922a1874d74e04cb8ca09 Mon Sep 17 00:00:00 2001 From: Matthew R Hermes Date: Fri, 3 Jan 2025 11:09:52 -0600 Subject: [PATCH 17/22] 'MC23': 'tMC23' alias in otfnal Since "MC23" is the name of the functional, not "tMC23". --- examples/mcpdft/03-metaGGA_functionals.py | 7 ++----- pyscf/mcpdft/otfnal.py | 3 +++ pyscf/mcpdft/test/test_mgga.py | 8 ++++---- 3 files changed, 9 insertions(+), 9 deletions(-) diff --git a/examples/mcpdft/03-metaGGA_functionals.py b/examples/mcpdft/03-metaGGA_functionals.py index 37e5327e..394436b4 100644 --- a/examples/mcpdft/03-metaGGA_functionals.py +++ b/examples/mcpdft/03-metaGGA_functionals.py @@ -18,16 +18,13 @@ mc = mcpdft.CASCI(mf, tM06L0, 6, 8).run () # MC23: meta-hybrid on-top functional [ChemRxiv. 2024; doi:10.26434/chemrxiv-2024-5hc8g-v2] -# To calculate the MC23 energies, you can define functional as 'tMC23'. - -# Note: It's 'fully-translated' version 'ftMC23' is not defined as of now. # State-Specific -mc = mcpdft.CASCI(mf, 'tMC23', 6, 8) +mc = mcpdft.CASCI(mf, 'MC23', 6, 8) mc.kernel() # State-average nroots=2 -mc = mcpdft.CASCI(mf, 'tMC23', 6, 8) +mc = mcpdft.CASCI(mf, 'MC23', 6, 8) mc.fcisolver.nroots=nroots mc.kernel()[0] diff --git a/pyscf/mcpdft/otfnal.py b/pyscf/mcpdft/otfnal.py index 98c19d87..9700827e 100644 --- a/pyscf/mcpdft/otfnal.py +++ b/pyscf/mcpdft/otfnal.py @@ -32,6 +32,7 @@ FT_B = getattr(__config__, 'mcpdft_otfnal_ftransfnal_B', -379.47331922) FT_C = getattr(__config__, 'mcpdft_otfnal_ftransfnal_C', -85.38149682) +OT_ALIAS = {'MC23': 'tMC23'} OT_HYB_ALIAS = {'PBE0' : '0.25*HF + 0.75*PBE, 0.25*HF + 0.75*PBE', 'MC23' : 'mc23'} # Note: mc23 is hybrid Meta-GGA OT Functional. @@ -836,6 +837,8 @@ def _sanity_check_ftot(xc_code): raise NotImplementedError(msg) def get_transfnal (mol, otxc): + if otxc.upper () in OT_ALIAS: + otxc = OT_ALIAS[otxc.upper ()] if otxc.upper ().startswith ('T'): xc_base = otxc[1:] fnal_class = transfnal diff --git a/pyscf/mcpdft/test/test_mgga.py b/pyscf/mcpdft/test/test_mgga.py index f5d59cbf..b798d2a5 100644 --- a/pyscf/mcpdft/test/test_mgga.py +++ b/pyscf/mcpdft/test/test_mgga.py @@ -25,7 +25,7 @@ molecule at the state-specific and state-average (2-states) using 1. Meta-GGA functional (tM06L) 2. Hybrid-meta-GGA functional tM06L0 -3. tMC23 Functional +3. MC23 Functional Test the MCPDFT energies calculated for the triplet water molecule at the 4. Meta-GGA functional (M06L) @@ -114,12 +114,12 @@ def setUpModule(): global get_lih, lih_tm06l, lih_tmc23, lih_tm06l_sa2, lih_tmc23_sa2,lih_tm06l0 global get_water_triplet, water_tm06l, water_tmc23 lih_tm06l = get_lih(1.5, functional='tM06L') - lih_tmc23 = get_lih(1.5, functional='tMC23') + lih_tmc23 = get_lih(1.5, functional='MC23') lih_tm06l_sa2 = get_lih(1.5, stateaverage=True, functional='tM06L') - lih_tmc23_sa2 = get_lih(1.5, stateaverage=True, functional='tMC23') + lih_tmc23_sa2 = get_lih(1.5, stateaverage=True, functional='MC23') lih_tm06l0 = get_lih(1.5, functional='tM06L0') water_tm06l = get_water_triplet() - water_tmc23 = get_water_triplet(functional='tMC23') + water_tmc23 = get_water_triplet(functional='MC23') def tearDownModule(): global lih_tm06l, lih_tmc23, lih_tm06l_sa2, lih_tmc23_sa2 From 1625e2eba7608655881a76b8125605126a09d560 Mon Sep 17 00:00:00 2001 From: Bhavnesh Jangid <85997686+JangidBhavnesh@users.noreply.github.com> Date: Fri, 3 Jan 2025 11:30:22 -0600 Subject: [PATCH 18/22] Updated the citation. --- examples/mcpdft/03-metaGGA_functionals.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/examples/mcpdft/03-metaGGA_functionals.py b/examples/mcpdft/03-metaGGA_functionals.py index 394436b4..b61c0197 100644 --- a/examples/mcpdft/03-metaGGA_functionals.py +++ b/examples/mcpdft/03-metaGGA_functionals.py @@ -8,7 +8,7 @@ mf = scf.RHF (mol).run () -# The translation of Meta-GGAs and hybrid-Meta-GGAs [ChemRxiv. 2024; doi:10.26434/chemrxiv-2024-5hc8g-v2] +# The translation of Meta-GGAs and hybrid-Meta-GGAs [PNAS, 122, 1, 2025, e2419413121; https://doi.org/10.1073/pnas.2419413121] # Translated-Meta-GGA mc = mcpdft.CASCI(mf, 'tM06L', 6, 8).run () @@ -17,7 +17,7 @@ tM06L0 = 't' + mcpdft.hyb('M06L',0.25, hyb_type='average') mc = mcpdft.CASCI(mf, tM06L0, 6, 8).run () -# MC23: meta-hybrid on-top functional [ChemRxiv. 2024; doi:10.26434/chemrxiv-2024-5hc8g-v2] +# MC23: meta-hybrid on-top functional [PNAS, 122, 1, 2025, e2419413121; https://doi.org/10.1073/pnas.2419413121] # State-Specific mc = mcpdft.CASCI(mf, 'MC23', 6, 8) From 1a2db3c1029e5038fc286f43fe4202422a8c13a4 Mon Sep 17 00:00:00 2001 From: hebrewsnabla Date: Sat, 4 Jan 2025 14:37:30 +0800 Subject: [PATCH 19/22] copy pyscf solution --- pyscf/lib/CMakeLists.txt | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pyscf/lib/CMakeLists.txt b/pyscf/lib/CMakeLists.txt index bb588cc7..18fa0b98 100644 --- a/pyscf/lib/CMakeLists.txt +++ b/pyscf/lib/CMakeLists.txt @@ -89,14 +89,14 @@ if(NOT PYSCF_SOURCE_DIR) endif() message(STATUS "Include pyscf source dir: ${PYSCF_SOURCE_DIR}") include_directories(${PYSCF_SOURCE_DIR}/lib ${PYSCF_SOURCE_DIR}/lib/deps/include) -include_directories(${CINT_DIR}/include) +include_directories(${CMAKE_INSTALL_PREFIX}/include) configure_file( "${PYSCF_SOURCE_DIR}/lib/config.h.in" "${PYSCF_SOURCE_DIR}/lib/config.h") # to find config.h link_directories(${PYSCF_SOURCE_DIR}/lib/deps/lib ${PYSCF_SOURCE_DIR}/lib/deps/lib64) link_directories(${PYSCF_SOURCE_DIR}/lib) -link_directories(${XC_DIR}/lib) +link_directories(${CMAKE_INSTALL_PREFIX}/lib ${CMAKE_INSTALL_PREFIX}/lib64) message(STATUS "${PYSCF_SOURCE_DIR}/lib may need to be put in the environment LD_LIBRARY_PATH") # See also https://gitlab.kitware.com/cmake/community/wikis/doc/cmake/RPATH-handling From ba533d7a9471ac0bd79226bca50dbc8aa319540b Mon Sep 17 00:00:00 2001 From: JangidBhavnesh Date: Wed, 8 Jan 2025 16:09:08 -0600 Subject: [PATCH 20/22] Sanity check added for cell object --- pyscf/mcpdft/__init__.py | 18 ++++++++++++++++++ 1 file changed, 18 insertions(+) diff --git a/pyscf/mcpdft/__init__.py b/pyscf/mcpdft/__init__.py index 9462373a..1821bcd0 100644 --- a/pyscf/mcpdft/__init__.py +++ b/pyscf/mcpdft/__init__.py @@ -21,11 +21,29 @@ from pyscf.lib import logger from pyscf.mcscf import mc1step, casci +def _sanity_check__of_mol(mc_or_mf_or_mol): + ''' + Sanity check to ensure input is a mol object, not a cell object. + Raises an error for cell objects. + ''' + from pyscf.pbc import gto as pbcgto + if isinstance(mc_or_mf_or_mol, (mc1step.CASSCF, casci.CASCI)): + mol = mc_or_mf_or_mol._scf.mol + elif isinstance(mc_or_mf_or_mol, gto.Mole): + mol = mc_or_mf_or_mol + else: + mol = mc_or_mf_or_mol.mol + + if isinstance(mol, pbcgto.cell.Cell): + raise NotImplementedError("MCPDFT not implemented for PBC") + # NOTE: As of 02/06/2022, initializing PySCF mcscf classes with a symmetry-enabled molecule # doesn't work. def _MCPDFT (mc_class, mc_or_mf_or_mol, ot, ncas, nelecas, ncore=None, frozen=None, **kwargs): + # Raise an error if mol is a cell object. + _sanity_check__of_mol(mc_or_mf_or_mol) if isinstance (mc_or_mf_or_mol, (mc1step.CASSCF, casci.CASCI)): mc0 = mc_or_mf_or_mol mf_or_mol = mc_or_mf_or_mol._scf From 18ea1b121b2cec11eb01da9401b9f353ce0dc841 Mon Sep 17 00:00:00 2001 From: JangidBhavnesh Date: Thu, 9 Jan 2025 17:13:30 -0600 Subject: [PATCH 21/22] Typo fixed --- pyscf/mcpdft/__init__.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pyscf/mcpdft/__init__.py b/pyscf/mcpdft/__init__.py index 1821bcd0..edd4350f 100644 --- a/pyscf/mcpdft/__init__.py +++ b/pyscf/mcpdft/__init__.py @@ -21,7 +21,7 @@ from pyscf.lib import logger from pyscf.mcscf import mc1step, casci -def _sanity_check__of_mol(mc_or_mf_or_mol): +def _sanity_check_of_mol(mc_or_mf_or_mol): ''' Sanity check to ensure input is a mol object, not a cell object. Raises an error for cell objects. @@ -43,7 +43,7 @@ def _sanity_check__of_mol(mc_or_mf_or_mol): def _MCPDFT (mc_class, mc_or_mf_or_mol, ot, ncas, nelecas, ncore=None, frozen=None, **kwargs): # Raise an error if mol is a cell object. - _sanity_check__of_mol(mc_or_mf_or_mol) + _sanity_check_of_mol(mc_or_mf_or_mol) if isinstance (mc_or_mf_or_mol, (mc1step.CASSCF, casci.CASCI)): mc0 = mc_or_mf_or_mol mf_or_mol = mc_or_mf_or_mol._scf From 0c2586886d890bae40eb4833151f9527c42f90bd Mon Sep 17 00:00:00 2001 From: Matthew R Hermes Date: Fri, 10 Jan 2025 10:41:03 -0600 Subject: [PATCH 22/22] PySCF v2.8.0 compatibility msdft/test_noci: SCF behavior change required msdft test change. lrdf/test_lrdf: Lower numerical precision for unknown reason. --- pyscf/lrdf/test/test_lrdf.py | 4 ++-- pyscf/msdft/tests/test_noci.py | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/pyscf/lrdf/test/test_lrdf.py b/pyscf/lrdf/test/test_lrdf.py index 5f70520d..caa4d978 100644 --- a/pyscf/lrdf/test/test_lrdf.py +++ b/pyscf/lrdf/test/test_lrdf.py @@ -54,8 +54,8 @@ def test_get_jk_sr(self): vj, vk = with_df._get_jk_sr(dm[None]) with mol.with_range_coulomb(-0.15): jref, kref = hf.get_jk(mol, dm) - self.assertAlmostEqual(abs(jref-vj[0]).max(), 0, 11) - self.assertAlmostEqual(abs(kref-vk[0]).max(), 0, 11) + self.assertAlmostEqual(abs(jref-vj[0]).max(), 0, 8) + self.assertAlmostEqual(abs(kref-vk[0]).max(), 0, 8) if __name__ == "__main__": print('Full Tests for LRDF') diff --git a/pyscf/msdft/tests/test_noci.py b/pyscf/msdft/tests/test_noci.py index 9b25ea96..5c44ff58 100644 --- a/pyscf/msdft/tests/test_noci.py +++ b/pyscf/msdft/tests/test_noci.py @@ -19,7 +19,7 @@ def test_hf_det_ovlp(): occ = mf0.mo_occ.copy() occ[0][mf0.nelec[0]-1] = 0 occ[0][mf0.nelec[0]+1] = 1 - mf1 = scf.addons.mom_occ(mf1, mf0.mo_coeff, occ).run(max_cycle=1) + mf1 = scf.addons.mom_occ(mf1, mf0.mo_coeff, occ).run(max_cycle=1, mo_coeff=None) h, s = noci.hf_det_ovlp(ms_ks, [mf0, mf1]) ref = np.array([[-9.35176786e+01, -6.82503177e-02], [-6.82503177e-02, -9.33368874e+01]])