Merge fde961a into 041cdcb

adcc/OperatorIntegrals.py

@@ -137,6 +137,26 @@ def nabla(self):
         """
         return self.import_dipole_like_operator("nabla", is_symmetric=False)
 
+    @property
+    def density_dependent_operators(self):
+        if not hasattr(self.provider_ao, "density_dependent_operators"):
+            return {}
+        ret = {}
+        for op in self.provider_ao.density_dependent_operators:
+            ao_callback = self.provider_ao.density_dependent_operators[op]
+
+            def process_operator(dm, callback=ao_callback):
+                dm_ao = sum(dm.to_ao_basis())
+                v_ao = callback(dm_ao)
+                v_bb = replicate_ao_block(self.mospaces, v_ao, is_symmetric=True)
+                v_ff = OneParticleOperator(self.mospaces, is_symmetric=True)
+                transform_operator_ao2mo(
+                    v_bb, v_ff, self.__coefficients, self.__conv_tol
+                )
+                return v_ff
+            ret[op] = process_operator
+        return ret
+
     @property
     def timer(self):
         ret = Timer()

adcc/backends/psi4.py

@@ -53,6 +53,17 @@ def magnetic_dipole(self):
     def nabla(self):
         return [-1.0 * np.asarray(comp) for comp in self.mints.ao_nabla()]
 
+    @property
+    def density_dependent_operators(self):
+        ret = {}
+        if hasattr(self.wfn, "pe_state"):
+            ret["pe_induction_elec"] = lambda dm: \
+                self.wfn.pe_state.get_pe_contribution(
+                    psi4.core.Matrix.from_array(dm.to_ndarray()),
+                    elec_only=True
+            )[1]
+        return ret
+
 
 class Psi4EriBuilder(EriBuilder):
     def __init__(self, wfn, n_orbs, n_orbs_alpha, n_alpha, n_beta, restricted):

adcc/backends/pyscf.py

@@ -21,6 +21,7 @@
 ##
 ## ---------------------------------------------------------------------
 import numpy as np
+import warnings
 
 from libadcc import HartreeFockProvider
 from adcc.misc import cached_property
@@ -38,7 +39,12 @@ def __init__(self, scfres):
 
     @cached_property
     def electric_dipole(self):
-        return list(self.scfres.mol.intor_symmetric('int1e_r', comp=3))
+        if hasattr(self.scfres, "with_solvent"):
+            if self.scfres.with_solvent.eef:
+                warnings.warn("Using modified dipole operator due to EEF.")
+            return self.scfres.with_solvent.effective_dipole_operator()
+        else:
+            return list(self.scfres.mol.intor_symmetric('int1e_r', comp=3))
 
     @cached_property
     def magnetic_dipole(self):
@@ -55,6 +61,16 @@ def nabla(self):
                 -1.0 * self.scfres.mol.intor('int1e_ipovlp', comp=3, hermi=2)
             )
 
+    @property
+    def density_dependent_operators(self):
+        ret = {}
+        if hasattr(self.scfres, "with_solvent"):
+            if isinstance(self.scfres.with_solvent, solvent.pol_embed.PolEmbed):
+                ret["pe_induction_elec"] = lambda dm: \
+                    self.scfres.with_solvent._exec_cppe(dm.to_ndarray(),
+                                                        elec_only=True)[1]
+        return ret
+
 
 # TODO: refactor ERI builder to be more general
 # IntegralBuilder would be good

adcc/backends/test_backends_polembed.py

@@ -22,9 +22,13 @@
 ## ---------------------------------------------------------------------
 import unittest
 import itertools
+import numpy as np
 import adcc
 import adcc.backends
 
+from adcc.adc_pp.matrix import AdcBlock
+from adcc import OneParticleOperator, AmplitudeVector
+
 from numpy.testing import assert_allclose
 
 import pytest
@@ -52,7 +56,7 @@
 @pytest.mark.skipif(len(backends) == 0, reason="No backend found.")
 @expand_test_templates(list(itertools.product(basissets, methods, backends)))
 class TestPolarizableEmbedding(unittest.TestCase):
-    def template_pe_formaldehyde(self, basis, method, backend):
+    def template_pe_formaldehyde_perturbative(self, basis, method, backend):
         basename = f"formaldehyde_{basis}_pe_{method}"
         qc_result = qchem_data[basename]
         pe_options = {"potfile": pe_potentials["fa_6w"]}
@@ -82,3 +86,39 @@ def template_pe_formaldehyde(self, basis, method, backend):
             state.pe_ptlr_correction,
             atol=1e-5
         )
+
+    def test_pe_formaldehyde_coupling_sto3g_pe_adc2(self):
+        basis = "sto-3g"
+        method = "adc2"
+        backend = "pyscf"
+        pe_options = {"potfile": pe_potentials["fa_6w"]}
+        scfres = cached_backend_hf(backend, "formaldehyde", basis,
+                                   pe_options=pe_options)
+        assert_allclose(scfres.energy_scf, -112.36904349555, atol=1e-8)
+
+        # construct a normal ADC matrix
+        matrix = adcc.AdcMatrix(method, scfres)
+
+        # additional matrix apply for PE
+        def block_ph_ph_0_pe(hf, mp, intermediates):
+            op = hf.operators
+
+            def apply(ampl):
+                tdm = OneParticleOperator(mp, is_symmetric=False)
+                tdm.vo = ampl.ph.transpose()
+                vpe = op.density_dependent_operators["pe_induction_elec"](tdm)
+                return AmplitudeVector(ph=vpe.ov)
+            return AdcBlock(apply, 0)
+
+        # register the additional function with the matrix
+        # this should be implemented in the AdcMatrix ctor
+        # (to also get all the diagonal setup right)
+        matrix.blocks_ph['ph_ph_pe'] = block_ph_ph_0_pe(
+            matrix.reference_state, matrix.ground_state, matrix.intermediates
+        )
+        assert_allclose(matrix.ground_state.energy(2), -112.4804685653, atol=1e-8)
+        exci_tm = np.array([
+            0.1733632144, 0.3815240343, 0.5097618218, 0.5189443358, 0.5670575778
+        ])
+        state = adcc.run_adc(matrix, n_singlets=5, conv_tol=1e-7)
+        assert_allclose(state.excitation_energy_uncorrected, exci_tm, atol=1e-6)