Calculate pdft_feff (#38)

* initial stuff * adding pdft_feff stuff, nothing new yet, just some interface stuff * small edit * add lazy feff kernel and extract pdft_eff funcs * fix imports * update, make real kernel * add IDE to gitignore * fix call sig * add test for ao2mo feff * feff with state-average density contraction * trying to understand how to test feff * not working but ohwell * got it to agree for ncore = 0 case... * add doc strings * add more doc strings for pdft_feff * fix? idk it still doesn't want to work :( * works ! * working again * clean up test * fix linter * update pdft_feff interface for general c_dm1s and c_cascm2, and cache veff1, veff2 in lpdft object * update reference for L-PDFT in example * update reference for L-PDFT * fix lint
pyscf · Sep 16, 2023 · f80d249 · f80d249
1 parent 5e4236c
commit f80d249
Show file tree

Hide file tree

Showing 16 changed files with 1,971 additions and 1,158 deletions.
diff --git a/.gitignore b/.gitignore
@@ -127,3 +127,7 @@ dmypy.json
 
 # Pyre type checker
 .pyre/
+
+# IDEs
+.idea
+.vscode
diff --git a/doc/mcpdft/README.md b/doc/mcpdft/README.md
@@ -37,7 +37,7 @@ Features
       and/or point groups)
     - Extended multi-state MC-PDFT (XMS-PDFT): [*Faraday Discuss* **2020**, 224, 348-372]
     - Compressed multi-state MC-PDFT (CMS-PDFT): [*JCTC* **2020**, *16*, 7444]
-    - Linearized PDFT (L-PDFT): [*JCTC* **2023**]
+    - Linearized PDFT (L-PDFT): [*JCTC* **2023**, *19*, 3172]
 * On-the-fly generation of on-top density functionals from underlying KS-DFT
   'LDA' or 'GGA' exchange-correlation functionals as defined in Libxc.
     - Translated functionals: [*JCTC* **2014**, *10*, 3669]
@@ -80,7 +80,7 @@ Features
 [*JCTC* **2021**, *17*, 2775]: http://dx.doi.org/10.1021/acs.jctc.0c01346
 [*Mol Phys* **2022**, 120]: http://dx.doi.org/10.1080/00268976.2022.2110534
 [*Faraday Discuss* **2020**, 224, 348-372]: http://dx.doi.org/10.1039/D0FD00037J
-[*JCTC* **2023**]: https://dx.doi.org/10.1021/acs.jctc.3c00207
+[*JCTC* **2023**, *19*, 3172]: https://dx.doi.org/10.1021/acs.jctc.3c00207
 
 [comment]: <> (Code hyperlinks)
 [examples/mcpdft/02-hybrid_functionals.py]: examples/mcpdft/02-hybrid_functionals.py

diff --git a/examples/mcpdft/42-linearized-pdft.py b/examples/mcpdft/42-linearized-pdft.py
@@ -3,9 +3,9 @@
 '''
 linearized pair-density functional theory
 
-A multi-states extension of MC-PDFT for states close in energy.
-Generates an effective L-PDFT Hamiltonian through a Taylor expansion of
-the MC-PDFT energy expression [ChemRxiv: 10.26434/chemrxiv-2023-7d0gv].
+A multi-states extension of MC-PDFT for states close in energy. Generates an
+effective L-PDFT Hamiltonian through a Taylor expansion of the MC-PDFT energy
+expression [J Chem Theory Comput: 10.1021/acs.jctc.3c00207].
 '''
 
 from pyscf import gto, scf, mcpdft

diff --git a/examples/mcpdft/43-multi_state_mix.py b/examples/mcpdft/43-multi_state_mix.py
@@ -2,8 +2,8 @@
 
 
 '''
-Perform multi-state PDFT averaging over states of different spins
-and/or spatial symmetry
+Perform multi-state PDFT averaging over states of different spins and/or
+spatial symmetry
 
 The mcpdft.multi_state function maybe not generate the right spin or spatial
 symmetry as one needs.  This example shows how to put states with different
@@ -68,4 +68,4 @@
 
 mc = mcpdft.CASSCF(mf, "tPBE", 4, 4, grids_level=1)
 mc = mc.multi_state_mix([solver1, solver2], weights, "lin")
-mc.kernel()
+mc.kernel()
diff --git a/pyscf/grad/mcpdft.py b/pyscf/grad/mcpdft.py
@@ -17,17 +17,17 @@
 from pyscf.grad import rks as rks_grad
 from pyscf.dft import gen_grid
 from pyscf.lib import logger, pack_tril, current_memory, einsum, tag_array
-#from pyscf.grad import sacasscf
 from pyscf.grad import sacasscf
 from pyscf.mcscf.casci import cas_natorb
+
+from pyscf.mcpdft.pdft_eff import _contract_eff_rho
 from pyscf.mcpdft.otpd import get_ontop_pair_density, _grid_ao2mo
-from pyscf.mcpdft.pdft_veff import _contract_vot_rho, _contract_ao_vao
 from pyscf.mcpdft import _dms
 from functools import reduce
 from itertools import product
 from scipy import linalg
 import numpy as np
-import time, gc
+import gc
 
 BLKSIZE = gen_grid.BLKSIZE
 
@@ -220,7 +220,7 @@ def mcpdft_HellmanFeynman_grad (mc, ot, veff1, veff2, mo_coeff=None, ci=None,
 
             # Vpq + Vpqrs * Drs ; I'm not sure why the list comprehension down
             # there doesn't break ao's stride order but I'm not complaining
-            vrho = _contract_vot_rho (vPi, rho.sum (0), add_vrho=vrho)
+            vrho = _contract_eff_rho (vPi, rho.sum (0), add_eff_rho=vrho)
             tmp_dv = np.stack ([ot.get_veff_1body (rho, Pi, [ao_i, moval_occ],
                 w0[ip0:ip1], kern=vrho) for ao_i in aoval], axis=0)
             tmp_dv = (tmp_dv * mo_occ[None,:,:]
@@ -261,7 +261,7 @@ def mcpdft_HellmanFeynman_grad (mc, ot, veff1, veff2, mo_coeff=None, ci=None,
             gc.collect ()
 
     for k, ia in enumerate(atmlst):
-        shl0, shl1, p0, p1 = aoslices[ia]
+        p0, p1 = aoslices[ia][2:]
         h1ao = hcore_deriv(ia) # MRH: this should be the TRUE hcore
         de_hcore[k] += np.tensordot(h1ao, dm1)
         de_renorm[k] -= np.tensordot(s1[:,p0:p1], dme0[p0:p1]) * 2
@@ -386,8 +386,6 @@ def get_ham_response (self, state=None, atmlst=None, verbose=None, mo=None,
         if ci is None: ci = self.base.ci
         if (veff1 is None) or (veff2 is None):
             assert (False), kwargs
-            veff1, veff2 = self.base.get_pdft_veff (mo, ci[state],
-                incl_coul=True, paaa_only=True)
         fcasscf = self.make_fcasscf (state)
         fcasscf.mo_coeff = mo
         fcasscf.ci = ci[state]

diff --git a/pyscf/mcpdft/_libxc.py b/pyscf/mcpdft/_libxc.py
@@ -13,7 +13,6 @@
 # See the License for the specific language governing permissions and
 # limitations under the License.
 #
-import re
 from pyscf.dft.libxc import XC_ALIAS, XC_CODES, XC_KEYS
 from pyscf.dft.libxc import hybrid_coeff, rsh_coeff
 
@@ -52,7 +51,7 @@ def split_x_c_comma (xc):
     xc = xc.upper ()
     myerr = XCSplitError (xc)
     max_recurse = 5
-    for i in range (max_recurse):
+    for _ in range (max_recurse):
         if ',' in xc:
             break
         elif xc in XC_ALIAS_KEYS:
@@ -84,7 +83,7 @@ def split_x_c_comma (xc):
 
 def is_hybrid_or_rsh (xc_code):
     hyb = hybrid_coeff (xc_code)
-    omega, alpha, beta = rsh_coeff (xc_code)
+    omega = rsh_coeff (xc_code)[0]
     non0 = [abs (x)>1e-10 for x in (hyb, omega)]
     return any (non0)
 

diff --git a/pyscf/mcpdft/lpdft.py b/pyscf/mcpdft/lpdft.py
@@ -21,7 +21,7 @@
 
 from pyscf.lib import logger
 from pyscf.fci import direct_spin1
-from pyscf import mcpdft, lib
+from pyscf import mcpdft
 from pyscf.mcpdft import _dms
 from pyscf.mcscf.addons import StateAverageMCSCFSolver, \
     StateAverageMixFCISolver
@@ -49,7 +49,7 @@ def weighted_average_densities(mc, ci=None, weights=None):
     return _dms.make_weighted_casdm1s(mc, ci=ci, weights=weights), _dms.make_weighted_casdm2(mc, ci=ci, weights=weights)
 
 
-def get_lpdfthconst(mc, E_ot, veff1_0, veff2_0, casdm1s_0, casdm2_0, hyb=1.0, ncas=None, ncore=None):
+def get_lpdfthconst(mc, E_ot, casdm1s_0, casdm2_0, hyb=1.0, ncas=None, ncore=None):
     ''' Compute h_const for the L-PDFT Hamiltonian
 
     Args:
@@ -58,22 +58,13 @@ def get_lpdfthconst(mc, E_ot, veff1_0, veff2_0, casdm1s_0, casdm2_0, hyb=1.0, nc
         E_ot : float
             On-top energy
 
-        veff1_0 : ndarray with shape (nao, nao)
-            1-body effective potential in the AO basis.
-            Should not include classical Coulomb potential term.
-            Generated from expansion density
-
-        veff2_0 : pyscf.mcscf.mc_ao2mo._ERIS instance
-            Relevant 2-body effective potential in the MO basis.
-            Generated from expansion density.
-
         casdm1s_0 : ndarray of shape (2, ncas, ncas)
-            Spin-separated 1-RDM in the active space generated
-            from expansion density.
+            Spin-separated 1-RDM in the active space generated from expansion
+            density.
 
         casdm2_0 : ndarray of shape (ncas, ncas, ncas, ncas)
-            Spin-summed 2-RDM in the active space generated
-            from expansion density.
+            Spin-summed 2-RDM in the active space generated from expansion
+            density.
 
         hyb : float
             Hybridization constant (lambda term)
@@ -101,20 +92,19 @@ def get_lpdfthconst(mc, E_ot, veff1_0, veff2_0, casdm1s_0, casdm2_0, hyb=1.0, nc
     vj = mc._scf.get_j(dm=dm1)
     e_j = np.tensordot(vj, dm1) / 2
 
-    # One-electron on-top potential energy
-    e_veff1 = np.tensordot(veff1_0, dm1)
+    e_veff1 = np.tensordot(mc.veff1, dm1)
 
     # Deal with 2-electron on-top potential energy
-    e_veff2 = veff2_0.energy_core
-    e_veff2 += np.tensordot(veff2_0.vhf_c[ncore:nocc, ncore:nocc], casdm1_0)
-    e_veff2 += 0.5 * np.tensordot(mc.get_h2lpdft(veff2_0), casdm2_0, axes=4)
+    e_veff2 = mc.veff2.energy_core
+    e_veff2 += np.tensordot(mc.veff2.vhf_c[ncore:nocc, ncore:nocc], casdm1_0)
+    e_veff2 += 0.5 * np.tensordot(mc.get_h2lpdft(), casdm2_0, axes=4)
 
     # h_nuc + E_ot - 1/2 g_pqrs D_pq D_rs - V_pq D_pq - 1/2 v_pqrs d_pqrs
     energy_core = hyb * mc.energy_nuc() + E_ot - hyb * e_j - e_veff1 - e_veff2
     return energy_core
 
 
-def transformed_h1e_for_cas(mc, E_ot, veff1_0, veff2_0, casdm1s_0, casdm2_0, hyb=1.0,
+def transformed_h1e_for_cas(mc, E_ot, casdm1s_0, casdm2_0, hyb=1.0,
                             mo_coeff=None, ncas=None, ncore=None):
     '''Compute the CAS one-particle L-PDFT Hamiltonian
 
@@ -124,29 +114,20 @@ def transformed_h1e_for_cas(mc, E_ot, veff1_0, veff2_0, casdm1s_0, casdm2_0, hyb
         E_ot : float
             On-top energy
 
-        veff1_0 : ndarray with shape (nao, nao)
-            1-body effective potential in the AO basis.
-            Should not include classical Coulomb potential term.
-            Generated from expansion density.
-
-        veff2_0 : pyscf.mcscf.mc_ao2mo._ERIS instance
-            Relevant 2-body effecive potential in the MO basis.
-            Generated from expansion density.
-
         casdm1s_0 : ndarray of shape (2,ncas,ncas)
-            Spin-separated 1-RDM in the active space generated
-            from expansion density
+            Spin-separated 1-RDM in the active space generated from expansion
+            density
 
         casdm2_0 : ndarray of shape (ncas,ncas,ncas,ncas)
-            Spin-summed 2-RDM in the active space generated
-            from expansion density
+            Spin-summed 2-RDM in the active space generated from expansion
+            density
 
         hyb : float
             Hybridization constant (lambda term)
 
         mo_coeff : ndarray of shape (nao,nmo)
-            A full set of molecular orbital coefficients. Taken from
-            self if not provided.
+            A full set of molecular orbital coefficients. Taken from self if
+            not provided.
 
         ncas : int
             Number of active space molecular orbitals
@@ -155,8 +136,9 @@ def transformed_h1e_for_cas(mc, E_ot, veff1_0, veff2_0, casdm1s_0, casdm2_0, hyb
             Number of core molecular orbitals
 
     Returns:
-        A tuple, the first is the effective one-electron linear PDFT Hamiltonian
-        defined in CAS space, the second is the modified core energy.
+        A tuple, the first is the effective one-electron linear PDFT
+        Hamiltonian defined in CAS space, the second is the modified core
+        energy.
     '''
     if mo_coeff is None: mo_coeff = mc.mo_coeff
     if ncas is None: ncas = mc.ncas
@@ -171,31 +153,26 @@ def transformed_h1e_for_cas(mc, E_ot, veff1_0, veff2_0, casdm1s_0, casdm2_0, hyb
     v_j = mc._scf.get_j(dm=dm1)
 
     # h_pq + V_pq + J_pq all in AO integrals
-    hcore_eff = hyb * mc.get_hcore() + veff1_0 + hyb * v_j
-    energy_core = mc.get_lpdfthconst(E_ot, veff1_0, veff2_0, casdm1s_0,
-                                     casdm2_0, hyb)
+    hcore_eff = hyb * mc.get_hcore() + mc.veff1 + hyb * v_j
+    energy_core = mc.get_lpdfthconst(E_ot, casdm1s_0, casdm2_0, hyb)
 
     if mo_core.size != 0:
         core_dm = np.dot(mo_core, mo_core.conj().T) * 2
         # This is precomputed in MRH's ERIS object
-        energy_core += veff2_0.energy_core
+        energy_core += mc.veff2.energy_core
         energy_core += np.tensordot(core_dm, hcore_eff).real
 
     h1eff = reduce(np.dot, (mo_cas.conj().T, hcore_eff, mo_cas))
     # Add in the 2-electron portion that acts as a 1-electron operator
-    h1eff += veff2_0.vhf_c[ncore:nocc, ncore:nocc]
+    h1eff += mc.veff2.vhf_c[ncore:nocc, ncore:nocc]
 
     return h1eff, energy_core
 
 
-def get_transformed_h2eff_for_cas(mc, veff2_0, ncore=None, ncas=None):
+def get_transformed_h2eff_for_cas(mc, ncore=None, ncas=None):
     '''Compute the CAS two-particle linear PDFT Hamiltonian
 
     Args:
-        veff2_0 : pyscf.mcscf.mc_ao2mo._ERIS instance
-            Relevant 2-body effecive potential in the MO basis.
-            Generated from expansion density.
-
         ncore : int
             Number of core MOs
 
@@ -208,16 +185,16 @@ def get_transformed_h2eff_for_cas(mc, veff2_0, ncore=None, ncas=None):
     if ncore is None: ncore = mc.ncore
     if ncas is None: ncas = mc.ncas
     nocc = ncore + ncas
-    return veff2_0.papa[ncore:nocc, :, ncore:nocc, :]
+    return mc.veff2.papa[ncore:nocc, :, ncore:nocc, :]
 
 
 def make_lpdft_ham_(mc, mo_coeff=None, ci=None, ot=None):
     '''Compute the L-PDFT Hamiltonian
 
     Args:
         mo_coeff : ndarray of shape (nao, nmo)
-            A full set of molecular orbital coefficients. Taken from
-            self if not provided.
+            A full set of molecular orbital coefficients. Taken from self if
+            not provided.
 
         ci : list of ndarrays of length nroots
             CI vectors should be from a converged CASSCF/CASCI calculation
@@ -251,12 +228,12 @@ def make_lpdft_ham_(mc, mo_coeff=None, ci=None, ot=None):
     ncas = mc.ncas
     casdm1s_0, casdm2_0 = mc.get_casdm12_0()
 
-    veff1_0, veff2_0, E_ot = mc.get_pdft_veff(mo=mo_coeff, casdm1s=casdm1s_0,
+    mc.veff1, mc.veff2, E_ot = mc.get_pdft_veff(mo=mo_coeff, casdm1s=casdm1s_0,
                                         casdm2=casdm2_0, drop_mcwfn=True, incl_energy=True)
 
     # This is all standard procedure for generating the hamiltonian in PySCF
-    h1, h0 = mc.get_h1lpdft(E_ot, veff1_0, veff2_0, casdm1s_0, casdm2_0, hyb=1.0-cas_hyb)
-    h2 = mc.get_h2lpdft(veff2_0)
+    h1, h0 = mc.get_h1lpdft(E_ot, casdm1s_0, casdm2_0, hyb=1.0-cas_hyb)
+    h2 = mc.get_h2lpdft()
     h2eff = direct_spin1.absorb_h1e(h1, h2, ncas, mc.nelecas, 0.5)
 
     def construct_ham_slice(solver, slice, nelecas):
@@ -274,23 +251,21 @@ def construct_ham_slice(solver, slice, nelecas):
         return construct_ham_slice(direct_spin1, slice(0, len(ci)), mc.nelecas)
 
     # We have a StateAverageMix Solver
-    # Todo, should maybe initialize this in a more obvious way?
     mc._irrep_slices = []
     start = 0
     for solver in mc.fcisolver.fcisolvers:
         end = start + solver.nroots
         mc._irrep_slices.append(slice(start, end))
         start = end
 
-    return [construct_ham_slice(s, irrep, mc.fcisolver._get_nelec (s, mc.nelecas))
+    return [construct_ham_slice(s, irrep, mc.fcisolver._get_nelec(s, mc.nelecas))
             for s, irrep in zip(mc.fcisolver.fcisolvers, mc._irrep_slices)]
 
 
-def kernel(mc, mo_coeff=None, ci0=None, ot=None, verbose=logger.NOTE):
+def kernel(mc, mo_coeff=None, ci0=None, ot=None, **kwargs):
     if ot is None: ot = mc.otfnal
     if mo_coeff is None: mo_coeff = mc.mo_coeff
 
-    #log = logger.new_logger(mc, verbose)
     mc.optimize_mcscf_(mo_coeff=mo_coeff, ci0=ci0)
     mc.lpdft_ham = mc.make_lpdft_ham_(ot=ot)
 
@@ -319,22 +294,31 @@ class _LPDFT(mcpdft.MultiStateMCPDFTSolver):
         e_states : ndarray of shape (nroots)
             L-PDFT final energies of the adiabatic states
         ci : list of length (nroots) of ndarrays
-            CI vectors in the optimized adiabatic basis of MC-SCF. Related to the
-            L-PDFT adiabat CI vectors by the expansion coefficients ``si_pdft''.
+            CI vectors in the optimized adiabatic basis of MC-SCF. Related to
+            the L-PDFT adiabat CI vectors by the expansion coefficients
+            ``si_pdft''.
         si_pdft : ndarray of shape (nroots, nroots)
-            Expansion coefficients of the L-PDFT adiabats in terms of the optimized
+            Expansion coefficients of the L-PDFT adiabats in terms of the
+            optimized
             MC-SCF adiabats
         e_mcscf : ndarray of shape (nroots)
             Energies of the MC-SCF adiabatic states
         lpdft_ham : ndarray of shape (nroots, nroots)
             L-PDFT Hamiltonian in the MC-SCF adiabatic basis
+        veff1 : ndarray of shape (nao, nao)
+            1-body effective potential in the AO basis computed using the
+            zeroth-order densities.
+        veff2 : pyscf.mcscf.mc_ao2mo._ERIS instance
+            Relevant 2-body effective potential in the MO basis.
     '''
 
     def __init__(self, mc):
         self.__dict__.update(mc.__dict__)
-        keys = set(('lpdft_ham', 'si_pdft'))
+        keys = set(('lpdft_ham', 'si_pdft', 'veff1', 'veff2'))
         self.lpdft_ham = None
         self.si_pdft = None
+        self.veff1 = None
+        self.veff2 = None
         self._keys = set((self.__dict__.keys())).union(keys)
 
     @property