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/mcpdft/03-metaGGA_functionals.py b/examples/mcpdft/03-metaGGA_functionals.py new file mode 100644 index 00000000..b61c0197 --- /dev/null +++ b/examples/mcpdft/03-metaGGA_functionals.py @@ -0,0 +1,30 @@ +#!/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 [PNAS, 122, 1, 2025, e2419413121; https://doi.org/10.1073/pnas.2419413121] + +# 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 () + +# 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) +mc.kernel() + +# State-average +nroots=2 +mc = mcpdft.CASCI(mf, 'MC23', 6, 8) +mc.fcisolver.nroots=nroots +mc.kernel()[0] 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/lib/CMakeLists.txt b/pyscf/lib/CMakeLists.txt index 5a4b9f69..1fb802c7 100644 --- a/pyscf/lib/CMakeLists.txt +++ b/pyscf/lib/CMakeLists.txt @@ -89,13 +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(${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 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/mcpdft/__init__.py b/pyscf/mcpdft/__init__.py index 9462373a..edd4350f 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 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 a9581be6..e12e38b9 100644 --- a/pyscf/mcpdft/otfnal.py +++ b/pyscf/mcpdft/otfnal.py @@ -19,6 +19,7 @@ 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 @@ -31,7 +32,94 @@ 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'} +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. + +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_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, + -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 function registers the new on-top functional if it hasn't been + registered previously. + Args: + 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. + ''' + libxc_register_code = xc_code.lower () + libxc_base_code = preset['xc_base'] + ext_params = preset['ext_params'] + 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_registered_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) + 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 @@ -68,6 +156,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 = ot.Pi_deriv nao = mo_coeff.shape[0] ncas = casdm2.shape[0] @@ -84,7 +173,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) @@ -205,6 +294,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) @@ -253,11 +343,23 @@ 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) + 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" @@ -278,7 +380,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 @@ -568,6 +671,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) @@ -723,21 +827,39 @@ 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 () in OT_ALIAS: + otxc = OT_ALIAS[otxc.upper ()] if otxc.upper ().startswith ('T'): xc_base = otxc[1:] fnal_class = transfnal elif otxc.upper ().startswith ('FT'): xc_base = otxc[2:] + _sanity_check_ftot(xc_base) fnal_class = ftransfnal else: raise NotImplementedError ( 'On-top pair-density functional names other than "translated" (t) or ' '"fully-translated (ft).' ) + # Try to register the functional with libxc, if not already done + xc_base = _get_registered_ot_functional (xc_base, mol) + 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 ',' 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. ' @@ -749,8 +871,6 @@ def get_transfnal (mol, otxc): ks.xc = xc_base return fnal_class (ks) - - class colle_salvetti_corr (otfnal): @@ -801,19 +921,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_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/mcpdft/test/test_mgga.py b/pyscf/mcpdft/test/test_mgga.py new file mode 100644 index 00000000..b798d2a5 --- /dev/null +++ b/pyscf/mcpdft/test/test_mgga.py @@ -0,0 +1,220 @@ +#!/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 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 (tM06L) +2. Hybrid-meta-GGA functional tM06L0 +3. MC23 Functional + +Test the MCPDFT energies calculated for the triplet water molecule at the +4. Meta-GGA functional (M06L) +5. MC23 Functional + +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 +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 = mc.state_average_([0.5, 0.5]) + + 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, spin=2,output='/dev/null', verbose=0) + + mf = scf.RHF(mol).run() + + 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) + mc.fcisolver = solver1 + mc = mc.run() + return mc + +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='MC23') + lih_tm06l_sa2 = get_lih(1.5, stateaverage=True, functional='tM06L') + 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='MC23') + +def tearDownModule(): + 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): + + 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_tot + + 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] + SA_E_MCPDFT_EXPECTED = [-7.95807682, -7.79920022] + + 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_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 + + 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 + + # 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] + 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_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 + + 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 + + # 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 + + self.assertAlmostEqual(e_mcscf, E_CASSCF_EXPECTED, 6) + self.assertAlmostEqual(epdft, E_MCPDFT_EXPECTED, 6) + +if __name__ == "__main__": + print("Full Tests for MGGAs, Hybrid-MGGAs, and MC23") + unittest.main() diff --git a/pyscf/mcpdft/tfnal_derivs.py b/pyscf/mcpdft/tfnal_derivs.py index 199794e7..9383ceab 100644 --- a/pyscf/mcpdft/tfnal_derivs.py +++ b/pyscf/mcpdft/tfnal_derivs.py @@ -104,10 +104,11 @@ 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] - 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 +118,19 @@ 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 + + if nderiv > 4 and dderiv > 0: + raise NotImplementedError("Meta-GGA functional derivatives") + + if 1 < nderiv <= 4: + rho_deriv = rho_tot[1:4, :] + elif 4 < nderiv <= 6: + rho_deriv = rho_tot[1:6, :] + else: + rho_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, verbose=otfnal.verbose)[:dderiv + 1] 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..46ca1b62 --- /dev/null +++ b/pyscf/msdft/noci.py @@ -0,0 +1,418 @@ +#!/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 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 + 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): + ''' + 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..5c44ff58 --- /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, 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]]) + 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