LR-PCM ADC for PSI4 (#132)

* lr_pcm. should work

* fix backend pcm detection for psi4

* some comments and cleanup

* rename to potential

* add test for PCM/PSI4 with PSI4 CIS reference

* introduced env variable for cashed_backend_hf

* some cleanup

* fix problem with PSI4 pe + pcm tests

* reviews

* fix yaml generator

* cleanup + psi4 specific tests

* some more cleanup

* first try implementing ptlr pcm

* add some tests - 1) versus LR - 2) against itself

* increased printout in yaml and some cleanup

* add LR and ptLR for pyscf + tests

* review + creating PCM and moving PE examples

* fix test (cc-pvdz only in full mode)

* documentation

* review

* adding to authors

Co-authored-by: Jonas Leitner <j.leitner@stud.uni-heidelberg.de>

adcc/OperatorIntegrals.py

@@ -102,7 +102,8 @@ def available(self):
             "electric_dipole",
             "magnetic_dipole",
             "nabla",
-            "pe_induction_elec"
+            "pe_induction_elec",
+            "pcm_potential_elec"
         )
         return [integral for integral in integrals
                 if hasattr(self.provider_ao, integral)]
@@ -187,6 +188,19 @@ def pe_induction_elec(self):
         callback = self.provider_ao.pe_induction_elec
         return self.__import_density_dependent_operator(callback)
 
+    @property
+    def pcm_potential_elec(self):
+        """
+        Returns a function to obtain the (density-dependent)
+        electronic PCM potential operator in the molecular orbital basis
+        """
+        if "pcm_potential_elec" not in self.available:
+            raise NotImplementedError("Electronic PCM potential operator "
+                                      "not implemented "
+                                      f"in {self.provider_ao.backend} backend.")
+        callback = self.provider_ao.pcm_potential_elec
+        return self.__import_density_dependent_operator(callback)
+
     @property
     def timer(self):
         ret = Timer()

adcc/__init__.py

@@ -69,7 +69,7 @@
 __version__ = "0.15.11"
 __license__ = "GPL v3"
 __url__ = "https://adc-connect.org"
-__authors__ = ["Michael F. Herbst", "Maximilian Scheurer"]
+__authors__ = ["Michael F. Herbst", "Maximilian Scheurer", "Jonas Leitner"]
 __email__ = "developers@adc-connect.org"
 __contributors__ = []
 

adcc/adc_pp/environment.py

@@ -41,3 +41,22 @@ def apply(ampl):
         vpe = op.pe_induction_elec(tdm)
         return AmplitudeVector(ph=vpe.ov)
     return AdcBlock(apply, 0)
+
+
+def block_ph_ph_0_pcm(hf, mp, intermediates):
+    """
+    Constructs an :py:class:`AdcBlock` that describes the
+    linear response coupling to the polarizable environment
+    from PCM via a CIS-like transition density as described
+    in 10.1021/ct300763v, eq 63. Since the contribution
+    depends on the input amplitude itself,
+    a diagonal term cannot be formulated.
+    """
+    op = hf.operators
+
+    def apply(ampl):
+        tdm = OneParticleOperator(mp, is_symmetric=False)
+        tdm.vo = ampl.ph.transpose()
+        vpcm = op.pcm_potential_elec(tdm)
+        return AmplitudeVector(ph=vpcm.ov)
+    return AdcBlock(apply, 0)

adcc/backends/psi4.py

@@ -64,6 +64,16 @@ def pe_induction_elec_ao(dm):
                 )[1]
             return pe_induction_elec_ao
 
+    @property
+    def pcm_potential_elec(self):
+        if self.wfn.PCM_enabled():
+            def pcm_potential_elec_ao(dm):
+                return psi4.core.PCM.compute_V(
+                    self.wfn.get_PCM(),
+                    psi4.core.Matrix.from_array(dm.to_ndarray())
+                )
+            return pcm_potential_elec_ao
+
 
 class Psi4EriBuilder(EriBuilder):
     def __init__(self, wfn, n_orbs, n_orbs_alpha, n_alpha, n_beta, restricted):
@@ -106,6 +116,14 @@ def pe_energy(self, dm, elec_only=True):
                                                         elec_only=elec_only)
         return e_pe
 
+    def pcm_energy(self, dm):
+        psi_dm = psi4.core.Matrix.from_array(dm.to_ndarray())
+        # computes the Fock matrix contribution.
+        # By contraction with the tdm, the electronic energy contribution is
+        # calculated.
+        V_pcm = psi4.core.PCM.compute_V(self.wfn.get_PCM(), psi_dm).to_array()
+        return np.einsum("uv,uv->", dm.to_ndarray(), V_pcm)
+
     @property
     def excitation_energy_corrections(self):
         ret = []
@@ -121,13 +139,21 @@ def excitation_energy_corrections(self):
                                             elec_only=True)
             )
             ret.extend([ptlr, ptss])
+        if self.environment == "pcm":
+            ptlr = EnergyCorrection(
+                "pcm_ptlr_correction",
+                lambda view: self.pcm_energy(view.transition_dm_ao)
+            )
+            ret.extend([ptlr])
         return {ec.name: ec for ec in ret}
 
     @property
     def environment(self):
         ret = None
         if hasattr(self.wfn, "pe_state"):
             ret = "pe"
+        elif self.wfn.PCM_enabled():
+            ret = "pcm"
         return ret
 
     def get_backend(self):
@@ -245,6 +271,8 @@ def run_hf(xyz, basis, charge=0, multiplicity=1, conv_tol=1e-11,
         "ccpvdz": "cc-pvdz",
     }
 
+    # needed for PE and PCM tests
+    psi4.core.clean_options()
     mol = psi4.geometry(f"""
         {charge} {multiplicity}
         {xyz}

adcc/backends/pyscf.py

@@ -29,7 +29,8 @@
 from ..exceptions import InvalidReference
 from ..ExcitedStates import EnergyCorrection
 
-from pyscf import ao2mo, gto, scf, solvent
+from pyscf import ao2mo, gto, scf
+from pyscf.solvent import pol_embed, ddcosmo
 
 
 class PyScfOperatorIntegralProvider:
@@ -59,13 +60,28 @@ def nabla(self):
     @property
     def pe_induction_elec(self):
         if hasattr(self.scfres, "with_solvent"):
-            if isinstance(self.scfres.with_solvent, solvent.pol_embed.PolEmbed):
+            if isinstance(self.scfres.with_solvent, pol_embed.PolEmbed):
                 def pe_induction_elec_ao(dm):
                     return self.scfres.with_solvent._exec_cppe(
                         dm.to_ndarray(), elec_only=True
                     )[1]
                 return pe_induction_elec_ao
 
+    @property
+    def pcm_potential_elec(self):
+        if hasattr(self.scfres, "with_solvent"):
+            if isinstance(self.scfres.with_solvent, ddcosmo.DDCOSMO):
+                def pcm_potential_elec_ao(dm):
+                    # Since eps (dielectric constant) is the only solvent parameter
+                    # in pyscf and there is no solvent data available in the
+                    # program, the user needs to adjust scfres.with_solvent.eps
+                    # manually to the optical dielectric constant (if non
+                    # equilibrium solvation is desired).
+                    return self.scfres.with_solvent._B_dot_x(
+                        dm.to_ndarray()
+                    )
+                return pcm_potential_elec_ao
+
 
 # TODO: refactor ERI builder to be more general
 # IntegralBuilder would be good
@@ -125,6 +141,15 @@ def pe_energy(self, dm, elec_only=True):
         e_pe, _ = pe_state.kernel(dm.to_ndarray(), elec_only=elec_only)
         return e_pe
 
+    def pcm_energy(self, dm):
+        # Since eps (dielectric constant) is the only solvent parameter
+        # in pyscf and there is no solvent data available in the
+        # program, the user needs to adjust scfres.with_solvent.eps
+        # manually to the optical dielectric constant (if non
+        # equilibrium solvation is desired).
+        V_pcm = self.scfres.with_solvent._B_dot_x(dm.to_ndarray())
+        return np.einsum("uv,uv->", dm.to_ndarray(), V_pcm)
+
     @property
     def excitation_energy_corrections(self):
         ret = []
@@ -140,14 +165,22 @@ def excitation_energy_corrections(self):
                                             elec_only=True)
             )
             ret.extend([ptlr, ptss])
+        elif self.environment == "pcm":
+            ptlr = EnergyCorrection(
+                "pcm_ptlr_correction",
+                lambda view: self.pcm_energy(view.transition_dm_ao)
+            )
+            ret.extend([ptlr])
         return {ec.name: ec for ec in ret}
 
     @property
     def environment(self):
         ret = None
         if hasattr(self.scfres, "with_solvent"):
-            if isinstance(self.scfres.with_solvent, solvent.pol_embed.PolEmbed):
+            if isinstance(self.scfres.with_solvent, pol_embed.PolEmbed):
                 ret = "pe"
+            elif isinstance(self.scfres.with_solvent, ddcosmo.DDCOSMO):
+                ret = "pcm"
         return ret
 
     def get_backend(self):