From aa3b6c93905f7d54e67f255b0ae1da8adb5158c3 Mon Sep 17 00:00:00 2001
From: Jonas Leitner <j.leitner@stud.uni-heidelberg.de>
Date: Sun, 9 May 2021 09:55:10 +0200
Subject: [PATCH 01/16] lr_pcm. should work

---
 adcc/OperatorIntegrals.py  | 36 +++++++++++++++++++++++++++++++++++-
 adcc/adc_pp/environment.py | 10 ++++++++++
 adcc/backends/psi4.py      | 10 ++++++++++
 3 files changed, 55 insertions(+), 1 deletion(-)

diff --git a/adcc/OperatorIntegrals.py b/adcc/OperatorIntegrals.py
index 270529eb..f5877531 100644
--- a/adcc/OperatorIntegrals.py
+++ b/adcc/OperatorIntegrals.py
@@ -102,7 +102,8 @@ def available(self):
             "electric_dipole",
             "magnetic_dipole",
             "nabla",
-            "pe_induction_elec"
+            "pe_induction_elec",
+            "pcm_induction_elec"
         )
         return [integral for integral in integrals
                 if hasattr(self.provider_ao, integral)]
@@ -187,6 +188,39 @@ def pe_induction_elec(self):
         callback = self.provider_ao.pe_induction_elec
         return self.__import_density_dependent_operator(callback)
 
+    @property
+    def pcm_induction_elec(self):
+        """
+        Returns a function to obtain the (density-dependent)
+        PCM electronic induction operator in the molecular orbital basis
+        """
+        if "pcm_induction_elec" not in self.available:
+            raise NotImplementedError("PCM electronic induction operator "
+                                      "not implemented "
+                                      f"in {self.provider_ao.backend} backend.")
+        callback = self.provider_ao.pcm_induction_elec
+        return self.__import_density_dependent_operator(callback)
+
+    @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()
diff --git a/adcc/adc_pp/environment.py b/adcc/adc_pp/environment.py
index 1ee02d42..6467b6a8 100644
--- a/adcc/adc_pp/environment.py
+++ b/adcc/adc_pp/environment.py
@@ -41,3 +41,13 @@ 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):
+    op = hf.operators
+
+    def apply(ampl):
+        tdm = OneParticleOperator(mp, is_symmetric=False)
+        tdm.vo = ampl.ph.transpose()
+        vpe = op.pcm_induction_elec(tdm)
+        return AmplitudeVector(ph=vpe.ov)
+    return AdcBlock(apply, 0)
diff --git a/adcc/backends/psi4.py b/adcc/backends/psi4.py
index 94b3aa24..c1a9588b 100644
--- a/adcc/backends/psi4.py
+++ b/adcc/backends/psi4.py
@@ -64,6 +64,14 @@ def pe_induction_elec_ao(dm):
                 )[1]
             return pe_induction_elec_ao
 
+    @property
+    def pcm_induction_elec(self):
+        if self.wfn.PCM_enabled():
+            def pcm_induction_elec_ao(dm):
+                return psi4.core.PCM.compute_V(self.wfn.get_PCM(),
+                    psi4.core.Matrix.from_array(dm.to_ndarray()))
+            return pcm_induction_elec_ao
+
 
 class Psi4EriBuilder(EriBuilder):
     def __init__(self, wfn, n_orbs, n_orbs_alpha, n_alpha, n_beta, restricted):
@@ -128,6 +136,8 @@ def environment(self):
         ret = None
         if hasattr(self.wfn, "pe_state"):
             ret = "pe"
+        if hasattr(self.wfn, "PCM_enabled"):
+            ret = "pcm"
         return ret
 
     def get_backend(self):

From 1685605254e9aa6777cfcda6253c5461973dba5a Mon Sep 17 00:00:00 2001
From: Jonas Leitner <j.leitner@stud.uni-heidelberg.de>
Date: Wed, 12 May 2021 18:14:15 +0200
Subject: [PATCH 02/16] fix backend pcm detection for psi4

---
 adcc/backends/psi4.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/adcc/backends/psi4.py b/adcc/backends/psi4.py
index c1a9588b..5a5f55ab 100644
--- a/adcc/backends/psi4.py
+++ b/adcc/backends/psi4.py
@@ -136,7 +136,7 @@ def environment(self):
         ret = None
         if hasattr(self.wfn, "pe_state"):
             ret = "pe"
-        if hasattr(self.wfn, "PCM_enabled"):
+        if self.wfn.PCM_enabled():
             ret = "pcm"
         return ret
 

From 7859e4c6d97a7d3115aa3854398be4f9f379a156 Mon Sep 17 00:00:00 2001
From: Jonas Leitner <j.leitner@stud.uni-heidelberg.de>
Date: Tue, 25 May 2021 16:46:44 +0200
Subject: [PATCH 03/16] some comments and cleanup

---
 adcc/adc_pp/environment.py | 9 +++++++++
 adcc/backends/psi4.py      | 3 ++-
 2 files changed, 11 insertions(+), 1 deletion(-)

diff --git a/adcc/adc_pp/environment.py b/adcc/adc_pp/environment.py
index 6467b6a8..4160399f 100644
--- a/adcc/adc_pp/environment.py
+++ b/adcc/adc_pp/environment.py
@@ -42,7 +42,16 @@ def apply(ampl):
         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):
diff --git a/adcc/backends/psi4.py b/adcc/backends/psi4.py
index 5a5f55ab..dee3e922 100644
--- a/adcc/backends/psi4.py
+++ b/adcc/backends/psi4.py
@@ -68,7 +68,8 @@ def pe_induction_elec_ao(dm):
     def pcm_induction_elec(self):
         if self.wfn.PCM_enabled():
             def pcm_induction_elec_ao(dm):
-                return psi4.core.PCM.compute_V(self.wfn.get_PCM(),
+                return psi4.core.PCM.compute_V(
+                    self.wfn.get_PCM(),
                     psi4.core.Matrix.from_array(dm.to_ndarray()))
             return pcm_induction_elec_ao
 

From 9d94ce3bc45ba330a6a924d373a562fb8fd7b653 Mon Sep 17 00:00:00 2001
From: Jonas Leitner <j.leitner@stud.uni-heidelberg.de>
Date: Fri, 25 Jun 2021 09:20:46 +0200
Subject: [PATCH 04/16] rename to potential

---
 adcc/OperatorIntegrals.py  | 12 ++++++------
 adcc/adc_pp/environment.py |  2 +-
 adcc/backends/psi4.py      |  6 +++---
 3 files changed, 10 insertions(+), 10 deletions(-)

diff --git a/adcc/OperatorIntegrals.py b/adcc/OperatorIntegrals.py
index f5877531..65da20ed 100644
--- a/adcc/OperatorIntegrals.py
+++ b/adcc/OperatorIntegrals.py
@@ -103,7 +103,7 @@ def available(self):
             "magnetic_dipole",
             "nabla",
             "pe_induction_elec",
-            "pcm_induction_elec"
+            "pcm_potential_elec"
         )
         return [integral for integral in integrals
                 if hasattr(self.provider_ao, integral)]
@@ -189,16 +189,16 @@ def pe_induction_elec(self):
         return self.__import_density_dependent_operator(callback)
 
     @property
-    def pcm_induction_elec(self):
+    def pcm_potential_elec(self):
         """
         Returns a function to obtain the (density-dependent)
-        PCM electronic induction operator in the molecular orbital basis
+        electronic PCM potential operator in the molecular orbital basis
         """
-        if "pcm_induction_elec" not in self.available:
-            raise NotImplementedError("PCM electronic induction operator "
+        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_induction_elec
+        callback = self.provider_ao.pcm_potential_elec
         return self.__import_density_dependent_operator(callback)
 
     @property
diff --git a/adcc/adc_pp/environment.py b/adcc/adc_pp/environment.py
index 4160399f..c6c7d436 100644
--- a/adcc/adc_pp/environment.py
+++ b/adcc/adc_pp/environment.py
@@ -57,6 +57,6 @@ def block_ph_ph_0_pcm(hf, mp, intermediates):
     def apply(ampl):
         tdm = OneParticleOperator(mp, is_symmetric=False)
         tdm.vo = ampl.ph.transpose()
-        vpe = op.pcm_induction_elec(tdm)
+        vpe = op.pcm_potential_elec(tdm)
         return AmplitudeVector(ph=vpe.ov)
     return AdcBlock(apply, 0)
diff --git a/adcc/backends/psi4.py b/adcc/backends/psi4.py
index dee3e922..8ee97d72 100644
--- a/adcc/backends/psi4.py
+++ b/adcc/backends/psi4.py
@@ -65,13 +65,13 @@ def pe_induction_elec_ao(dm):
             return pe_induction_elec_ao
 
     @property
-    def pcm_induction_elec(self):
+    def pcm_potential_elec(self):
         if self.wfn.PCM_enabled():
-            def pcm_induction_elec_ao(dm):
+            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_induction_elec_ao
+            return pcm_potential_elec_ao
 
 
 class Psi4EriBuilder(EriBuilder):

From 6d78bdd68a4fea01d4e76f9e600cefeffc8e568e Mon Sep 17 00:00:00 2001
From: Jonas Leitner <j.leitner@stud.uni-heidelberg.de>
Date: Sun, 11 Jul 2021 14:45:32 +0200
Subject: [PATCH 05/16] add test for PCM/PSI4 with PSI4 CIS reference

---
 adcc/adc_pp/environment.py         |  4 +-
 adcc/backends/psi4.py              | 12 ++++-
 adcc/backends/pyscf.py             |  2 +-
 adcc/backends/test_backends_pcm.py | 71 ++++++++++++++++++++++++++++++
 adcc/backends/testing.py           | 10 +++--
 adcc/testdata/cache.py             |  1 +
 adcc/testdata/psi_dump.yml         | 12 +++++
 7 files changed, 105 insertions(+), 7 deletions(-)
 create mode 100644 adcc/backends/test_backends_pcm.py
 create mode 100644 adcc/testdata/psi_dump.yml

diff --git a/adcc/adc_pp/environment.py b/adcc/adc_pp/environment.py
index c6c7d436..aad254e0 100644
--- a/adcc/adc_pp/environment.py
+++ b/adcc/adc_pp/environment.py
@@ -57,6 +57,6 @@ def block_ph_ph_0_pcm(hf, mp, intermediates):
     def apply(ampl):
         tdm = OneParticleOperator(mp, is_symmetric=False)
         tdm.vo = ampl.ph.transpose()
-        vpe = op.pcm_potential_elec(tdm)
-        return AmplitudeVector(ph=vpe.ov)
+        vpcm = op.pcm_potential_elec(tdm)
+        return AmplitudeVector(ph=vpcm.ov)
     return AdcBlock(apply, 0)
diff --git a/adcc/backends/psi4.py b/adcc/backends/psi4.py
index 8ee97d72..e3afa58a 100644
--- a/adcc/backends/psi4.py
+++ b/adcc/backends/psi4.py
@@ -249,7 +249,7 @@ def import_scf(wfn):
 
 
 def run_hf(xyz, basis, charge=0, multiplicity=1, conv_tol=1e-11,
-           conv_tol_grad=1e-8, max_iter=150, pe_options=None):
+           conv_tol_grad=1e-8, max_iter=150, pe_options=None, pcm=None):
     basissets = {
         "sto3g": "sto-3g",
         "def2tzvp": "def2-tzvp",
@@ -278,6 +278,16 @@ def run_hf(xyz, basis, charge=0, multiplicity=1, conv_tol=1e-11,
         psi4.set_options({"pe": "true"})
         psi4.set_module_options("pe", {"potfile": pe_options["potfile"]})
 
+    if pcm:
+        psi4.set_options({"pcm": True, "pcm_scf_type": "total"})
+        psi4.pcm_helper("""
+        Units = AU
+        Cavity {Type = Gepol
+                Area = 0.3}
+        Medium {Solvertype = IEFPCM
+                Nonequilibrium = True
+                Solvent = Water}""")
+
     if multiplicity != 1:
         psi4.set_options({
             'reference': "UHF",
diff --git a/adcc/backends/pyscf.py b/adcc/backends/pyscf.py
index 90128750..3926a7b5 100644
--- a/adcc/backends/pyscf.py
+++ b/adcc/backends/pyscf.py
@@ -262,7 +262,7 @@ def import_scf(scfres):
 
 
 def run_hf(xyz, basis, charge=0, multiplicity=1, conv_tol=1e-11,
-           conv_tol_grad=1e-9, max_iter=150, pe_options=None):
+           conv_tol_grad=1e-9, max_iter=150, pe_options=None, pcm=None):
     mol = gto.M(
         atom=xyz,
         basis=basis,
diff --git a/adcc/backends/test_backends_pcm.py b/adcc/backends/test_backends_pcm.py
new file mode 100644
index 00000000..f8855123
--- /dev/null
+++ b/adcc/backends/test_backends_pcm.py
@@ -0,0 +1,71 @@
+import pytest
+import unittest
+import itertools
+import adcc
+import adcc.backends
+import os
+
+from adcc.misc import expand_test_templates
+from adcc.testdata.cache import psi_data
+from .testing import cached_backend_hf
+from numpy.testing import assert_allclose
+from adcc.exceptions import InputError
+
+from adcc.adc_pp.environment import block_ph_ph_0_pcm
+from adcc.AdcMatrix import AdcExtraTerm
+
+# remove pyscf until implemented
+backends = [b for b in adcc.backends.available()
+            if b not in ["molsturm", "veloxchem", "pyscf"]]
+basissets = ["sto3g", "ccpvdz"]
+methods = ["adc1", "adc2", "adc3"]
+
+
+@pytest.mark.skipif(len(backends) == 0, reason="No backend found.")
+@expand_test_templates(list(itertools.product(basissets, methods, backends)))
+class TestPCM(unittest.TestCase):
+    def template_pcm_linear_response_formaldehyde(self, basis, method, backend):
+        if method != "adc1":
+            pytest.skip("Reference only exists for adc1.")
+        basename = f"formaldehyde_{basis}_pcm_{method}"
+        psi_result = psi_data[basename]
+        scfres = cached_backend_hf(backend, "formaldehyde", basis,
+                                   pcm=True)
+        assert_allclose(scfres.energy_scf, psi_result["energy_scf"], atol=1e-8)
+
+        matrix = adcc.AdcMatrix(method, scfres)
+        solvent = AdcExtraTerm(matrix, {'ph_ph': block_ph_ph_0_pcm})
+
+        matrix += solvent
+        assert len(matrix.extra_terms)
+
+        state = adcc.run_adc(matrix, n_singlets=5, conv_tol=1e-7,
+                             environment=False)
+        assert_allclose(
+            state.excitation_energy_uncorrected,
+            psi_result["excitation_energy"],
+            atol=1e-5
+        )
+
+        # invalid combination
+        with pytest.raises(InputError):
+            adcc.run_adc(scfres, method=method, n_singlets=5,
+                         environment={"linear_response": True, "ptlr": True})
+
+        # no environment specified
+        with pytest.raises(InputError):
+            adcc.run_adc(scfres, method=method, n_singlets=5)
+
+        # automatically add coupling term
+        state = adcc.run_adc(scfres, method=method, n_singlets=5,
+                             conv_tol=1e-7, environment="linear_response")
+        assert_allclose(
+            state.excitation_energy_uncorrected,
+            psi_result["excitation_energy"],
+            atol=1e-5
+        )
+
+        # remove cavity files from PSI4 PCM calculations
+        for cavityfile in os.listdir(os.getcwd()):
+            if cavityfile.startswith(("cavity.off_", "PEDRA.OUT_")):
+                os.remove(cavityfile)
diff --git a/adcc/backends/testing.py b/adcc/backends/testing.py
index d41261df..5ceac43f 100644
--- a/adcc/backends/testing.py
+++ b/adcc/backends/testing.py
@@ -178,7 +178,7 @@ def operator_import_from_ao_test(scfres, ao_dict, operator="electric_dipole"):
 
 
 def cached_backend_hf(backend, molecule, basis, multiplicity=1, conv_tol=1e-12,
-                      pe_options=None):
+                      pe_options=None, pcm=None):
     """
     Run the SCF for a backend and a particular test case (if not done)
     and return the result.
@@ -195,7 +195,7 @@ def payload():
                                      basis=basis, conv_tol=conv_tol,
                                      multiplicity=multiplicity,
                                      conv_tol_grad=conv_tol_grad,
-                                     pe_options=pe_options)
+                                     pe_options=pe_options, pcm=pcm)
         return adcc.backends.import_scf_results(hfres)
 
     # For reasons not clear to me (mfh), caching does not work
@@ -203,7 +203,11 @@ def payload():
     if backend == "pyscf":
         return payload()
 
-    key = (backend, molecule, basis, str(multiplicity))
+    # otherwise PE and PCM share the same key
+    if pcm:
+        key = (backend, molecule, basis, str(multiplicity), "pcm")
+    else:
+        key = (backend, molecule, basis, str(multiplicity))
     try:
         return __cache_cached_backend_hf[key]
     except NameError:
diff --git a/adcc/testdata/cache.py b/adcc/testdata/cache.py
index 4dbfab89..bb9c2cd2 100644
--- a/adcc/testdata/cache.py
+++ b/adcc/testdata/cache.py
@@ -258,3 +258,4 @@ def read_yaml_data(fname):
 
 qchem_data = read_yaml_data("qchem_dump.yml")
 tmole_data = read_yaml_data("tmole_dump.yml")
+psi_data = read_yaml_data("psi_dump.yml")
diff --git a/adcc/testdata/psi_dump.yml b/adcc/testdata/psi_dump.yml
new file mode 100644
index 00000000..e5a730f3
--- /dev/null
+++ b/adcc/testdata/psi_dump.yml
@@ -0,0 +1,12 @@
+formaldehyde_sto3g_pcm_adc1:
+  basis: sto-3g
+  method: adc1
+  molecule: formaldehyde
+  energy_scf: -112.3545830298461254
+  excitation_energy: [0.16467, 0.36090, 0.46549, 0.50874, 0.69386]
+formaldehyde_ccpvdz_pcm_adc1:
+  basis: ccpvdz
+  method: adc1
+  molecule: formaldehyde
+  energy_scf: -113.8842706706941641
+  excitation_energy: [0.17900, 0.38177, 0.38802, 0.39216, 0.43544]
\ No newline at end of file

From bbd4a3f10b6b1e600a6aeae34b767b06b6b225d0 Mon Sep 17 00:00:00 2001
From: Jonas Leitner <j.leitner@stud.uni-heidelberg.de>
Date: Mon, 12 Jul 2021 11:11:54 +0200
Subject: [PATCH 06/16] introduced env variable for cashed_backend_hf

---
 adcc/backends/test_backends_pcm.py      |  2 +-
 adcc/backends/test_backends_polembed.py |  4 ++--
 adcc/backends/testing.py                | 14 +++++++++++---
 3 files changed, 14 insertions(+), 6 deletions(-)

diff --git a/adcc/backends/test_backends_pcm.py b/adcc/backends/test_backends_pcm.py
index f8855123..c1d00c41 100644
--- a/adcc/backends/test_backends_pcm.py
+++ b/adcc/backends/test_backends_pcm.py
@@ -30,7 +30,7 @@ def template_pcm_linear_response_formaldehyde(self, basis, method, backend):
         basename = f"formaldehyde_{basis}_pcm_{method}"
         psi_result = psi_data[basename]
         scfres = cached_backend_hf(backend, "formaldehyde", basis,
-                                   pcm=True)
+                                   env=("pcm", True))
         assert_allclose(scfres.energy_scf, psi_result["energy_scf"], atol=1e-8)
 
         matrix = adcc.AdcMatrix(method, scfres)
diff --git a/adcc/backends/test_backends_polembed.py b/adcc/backends/test_backends_polembed.py
index 2ddf7783..9339fc70 100644
--- a/adcc/backends/test_backends_polembed.py
+++ b/adcc/backends/test_backends_polembed.py
@@ -61,7 +61,7 @@ def template_pe_perturbative_formaldehyde(self, basis, method, backend):
         qc_result = qchem_data[basename]
         pe_options = {"potfile": pe_potentials["fa_6w"]}
         scfres = cached_backend_hf(backend, "formaldehyde", basis,
-                                   pe_options=pe_options)
+                                   env=("pe", pe_options))
         state = adcc.run_adc(scfres, method=method,
                              n_singlets=5, conv_tol=1e-10,
                              environment=["ptlr", "ptss"])
@@ -96,7 +96,7 @@ def template_pe_linear_response_formaldehyde(self, basis, method, backend):
         tm_result = tmole_data[basename]
         pe_options = {"potfile": pe_potentials["fa_6w"]}
         scfres = cached_backend_hf(backend, "formaldehyde", basis,
-                                   pe_options=pe_options)
+                                   env=("pe", pe_options))
         assert_allclose(scfres.energy_scf, tm_result["energy_scf"], atol=1e-8)
 
         matrix = adcc.AdcMatrix(method, scfres)
diff --git a/adcc/backends/testing.py b/adcc/backends/testing.py
index 5ceac43f..83adecf9 100644
--- a/adcc/backends/testing.py
+++ b/adcc/backends/testing.py
@@ -178,7 +178,7 @@ def operator_import_from_ao_test(scfres, ao_dict, operator="electric_dipole"):
 
 
 def cached_backend_hf(backend, molecule, basis, multiplicity=1, conv_tol=1e-12,
-                      pe_options=None, pcm=None):
+                      env=(None, None)):
     """
     Run the SCF for a backend and a particular test case (if not done)
     and return the result.
@@ -189,6 +189,14 @@ def cached_backend_hf(backend, molecule, basis, multiplicity=1, conv_tol=1e-12,
 
     global __cache_cached_backend_hf
 
+    # so you can only assign one of them.
+    pe_options = None
+    pcm = None
+    if env[0] == "pcm":
+        pcm = env[1]
+    elif env[0] == "pe":
+        pe_options = env[1]
+
     def payload():
         conv_tol_grad = 10 * conv_tol
         hfres = adcc.backends.run_hf(backend, xyz=static_data.xyz[molecule],
@@ -204,8 +212,8 @@ def payload():
         return payload()
 
     # otherwise PE and PCM share the same key
-    if pcm:
-        key = (backend, molecule, basis, str(multiplicity), "pcm")
+    if env[0]:
+        key = (backend, molecule, basis, str(multiplicity), env[0])
     else:
         key = (backend, molecule, basis, str(multiplicity))
     try:

From 4c495a88c2c9b6545d84cd6c704c406f54292b2a Mon Sep 17 00:00:00 2001
From: Jonas Leitner <j.leitner@stud.uni-heidelberg.de>
Date: Mon, 12 Jul 2021 11:24:55 +0200
Subject: [PATCH 07/16] some cleanup

---
 adcc/OperatorIntegrals.py | 20 --------------------
 1 file changed, 20 deletions(-)

diff --git a/adcc/OperatorIntegrals.py b/adcc/OperatorIntegrals.py
index 65da20ed..0bfb0b9a 100644
--- a/adcc/OperatorIntegrals.py
+++ b/adcc/OperatorIntegrals.py
@@ -201,26 +201,6 @@ def pcm_potential_elec(self):
         callback = self.provider_ao.pcm_potential_elec
         return self.__import_density_dependent_operator(callback)
 
-    @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()

From a36c658589c479b9a50d060a1708825f59ac6d44 Mon Sep 17 00:00:00 2001
From: Jonas Leitner <j.leitner@stud.uni-heidelberg.de>
Date: Mon, 12 Jul 2021 17:15:46 +0200
Subject: [PATCH 08/16] fix problem with PSI4 pe + pcm tests

---
 adcc/backends/psi4.py | 1 +
 1 file changed, 1 insertion(+)

diff --git a/adcc/backends/psi4.py b/adcc/backends/psi4.py
index e3afa58a..1e74b828 100644
--- a/adcc/backends/psi4.py
+++ b/adcc/backends/psi4.py
@@ -256,6 +256,7 @@ def run_hf(xyz, basis, charge=0, multiplicity=1, conv_tol=1e-11,
         "ccpvdz": "cc-pvdz",
     }
 
+    psi4.core.clean_options()
     mol = psi4.geometry(f"""
         {charge} {multiplicity}
         {xyz}

From 078ed34d8ebdf88ec388b94d93040bfa388f4dc5 Mon Sep 17 00:00:00 2001
From: Jonas Leitner <j.leitner@stud.uni-heidelberg.de>
Date: Thu, 15 Jul 2021 10:10:12 +0200
Subject: [PATCH 09/16] reviews

---
 adcc/backends/psi4.py                         | 18 +---
 adcc/backends/pyscf.py                        |  2 +-
 adcc/backends/test_backends_pcm.py            | 87 +++++++++++++++--
 adcc/backends/test_backends_polembed.py       |  4 +-
 adcc/backends/testing.py                      | 18 +---
 adcc/testdata/cache.py                        |  2 +-
 adcc/testdata/generate_psi4_data.py           | 97 +++++++++++++++++++
 adcc/testdata/generate_psi4_data_test.yml     | 34 +++++++
 adcc/testdata/{psi_dump.yml => psi4_dump.yml} |  4 +-
 9 files changed, 225 insertions(+), 41 deletions(-)
 create mode 100644 adcc/testdata/generate_psi4_data.py
 create mode 100644 adcc/testdata/generate_psi4_data_test.yml
 rename adcc/testdata/{psi_dump.yml => psi4_dump.yml} (78%)

diff --git a/adcc/backends/psi4.py b/adcc/backends/psi4.py
index 1e74b828..b6ec0ed5 100644
--- a/adcc/backends/psi4.py
+++ b/adcc/backends/psi4.py
@@ -70,7 +70,8 @@ def pcm_potential_elec(self):
             def pcm_potential_elec_ao(dm):
                 return psi4.core.PCM.compute_V(
                     self.wfn.get_PCM(),
-                    psi4.core.Matrix.from_array(dm.to_ndarray()))
+                    psi4.core.Matrix.from_array(dm.to_ndarray())
+                )
             return pcm_potential_elec_ao
 
 
@@ -137,7 +138,7 @@ def environment(self):
         ret = None
         if hasattr(self.wfn, "pe_state"):
             ret = "pe"
-        if self.wfn.PCM_enabled():
+        elif self.wfn.PCM_enabled():
             ret = "pcm"
         return ret
 
@@ -249,13 +250,14 @@ def import_scf(wfn):
 
 
 def run_hf(xyz, basis, charge=0, multiplicity=1, conv_tol=1e-11,
-           conv_tol_grad=1e-8, max_iter=150, pe_options=None, pcm=None):
+           conv_tol_grad=1e-8, max_iter=150, pe_options=None):
     basissets = {
         "sto3g": "sto-3g",
         "def2tzvp": "def2-tzvp",
         "ccpvdz": "cc-pvdz",
     }
 
+    # needed for PE and PCM tests
     psi4.core.clean_options()
     mol = psi4.geometry(f"""
         {charge} {multiplicity}
@@ -279,16 +281,6 @@ def run_hf(xyz, basis, charge=0, multiplicity=1, conv_tol=1e-11,
         psi4.set_options({"pe": "true"})
         psi4.set_module_options("pe", {"potfile": pe_options["potfile"]})
 
-    if pcm:
-        psi4.set_options({"pcm": True, "pcm_scf_type": "total"})
-        psi4.pcm_helper("""
-        Units = AU
-        Cavity {Type = Gepol
-                Area = 0.3}
-        Medium {Solvertype = IEFPCM
-                Nonequilibrium = True
-                Solvent = Water}""")
-
     if multiplicity != 1:
         psi4.set_options({
             'reference': "UHF",
diff --git a/adcc/backends/pyscf.py b/adcc/backends/pyscf.py
index 3926a7b5..90128750 100644
--- a/adcc/backends/pyscf.py
+++ b/adcc/backends/pyscf.py
@@ -262,7 +262,7 @@ def import_scf(scfres):
 
 
 def run_hf(xyz, basis, charge=0, multiplicity=1, conv_tol=1e-11,
-           conv_tol_grad=1e-9, max_iter=150, pe_options=None, pcm=None):
+           conv_tol_grad=1e-9, max_iter=150, pe_options=None):
     mol = gto.M(
         atom=xyz,
         basis=basis,
diff --git a/adcc/backends/test_backends_pcm.py b/adcc/backends/test_backends_pcm.py
index c1d00c41..1a49eaeb 100644
--- a/adcc/backends/test_backends_pcm.py
+++ b/adcc/backends/test_backends_pcm.py
@@ -3,11 +3,12 @@
 import itertools
 import adcc
 import adcc.backends
+import psi4
 import os
 
 from adcc.misc import expand_test_templates
-from adcc.testdata.cache import psi_data
-from .testing import cached_backend_hf
+from adcc.testdata import static_data
+from adcc.testdata.cache import psi4_data
 from numpy.testing import assert_allclose
 from adcc.exceptions import InputError
 
@@ -28,10 +29,21 @@ def template_pcm_linear_response_formaldehyde(self, basis, method, backend):
         if method != "adc1":
             pytest.skip("Reference only exists for adc1.")
         basename = f"formaldehyde_{basis}_pcm_{method}"
-        psi_result = psi_data[basename]
-        scfres = cached_backend_hf(backend, "formaldehyde", basis,
-                                   env=("pcm", True))
-        assert_allclose(scfres.energy_scf, psi_result["energy_scf"], atol=1e-8)
+        psi4_result = psi4_data[basename]
+
+        # so you can adapt some options more easily if needed
+        psi4_pcm_options = {"weight": 0.3, "pcm_method": "IEFPCM", "neq": True,
+                            "solvent": "Water"}
+        pcm_options = {"psi4": psi4_pcm_options}
+
+        hf_name = f"{backend}_run_pcm_hf"
+        hf_function = globals()[hf_name]
+        scfres = hf_function(static_data.xyz["formaldehyde"], basis, charge=0,
+                             multiplicity=1, conv_tol=1e-12,
+                             conv_tol_grad=1e-11, max_iter=150,
+                             pcm_options=pcm_options[backend])
+
+        assert_allclose(scfres.energy_scf, psi4_result["energy_scf"], atol=1e-8)
 
         matrix = adcc.AdcMatrix(method, scfres)
         solvent = AdcExtraTerm(matrix, {'ph_ph': block_ph_ph_0_pcm})
@@ -43,10 +55,16 @@ def template_pcm_linear_response_formaldehyde(self, basis, method, backend):
                              environment=False)
         assert_allclose(
             state.excitation_energy_uncorrected,
-            psi_result["excitation_energy"],
+            psi4_result["excitation_energy"],
             atol=1e-5
         )
 
+        state_cis = adcc.ExcitedStates(state, property_method="adc0")
+        assert_allclose(
+            state_cis.oscillator_strength,
+            psi4_result["osc_strength"], atol=1e-3
+        )
+
         # invalid combination
         with pytest.raises(InputError):
             adcc.run_adc(scfres, method=method, n_singlets=5,
@@ -61,7 +79,7 @@ def template_pcm_linear_response_formaldehyde(self, basis, method, backend):
                              conv_tol=1e-7, environment="linear_response")
         assert_allclose(
             state.excitation_energy_uncorrected,
-            psi_result["excitation_energy"],
+            psi4_result["excitation_energy"],
             atol=1e-5
         )
 
@@ -69,3 +87,56 @@ def template_pcm_linear_response_formaldehyde(self, basis, method, backend):
         for cavityfile in os.listdir(os.getcwd()):
             if cavityfile.startswith(("cavity.off_", "PEDRA.OUT_")):
                 os.remove(cavityfile)
+
+
+def psi4_run_pcm_hf(xyz, basis, charge=0, multiplicity=1, conv_tol=1e-12,
+                    conv_tol_grad=1e-11, max_iter=150, pcm_options=None):
+
+    # add something to check if pcm_options is defined and a dict?
+
+    basissets = {
+        "sto3g": "sto-3g",
+        "def2tzvp": "def2-tzvp",
+        "ccpvdz": "cc-pvdz",
+    }
+
+    # needed for PE and PCM tests
+    psi4.core.clean_options()
+    mol = psi4.geometry(f"""
+        {charge} {multiplicity}
+        {xyz}
+        symmetry c1
+        units au
+        no_reorient
+        no_com
+    """)
+
+    psi4.core.be_quiet()
+    psi4.set_options({
+        'basis': basissets.get(basis, basis),
+        'scf_type': 'pk',
+        'e_convergence': conv_tol,
+        'd_convergence': conv_tol_grad,
+        'maxiter': max_iter,
+        'reference': "RHF",
+        "pcm": True,
+        "pcm_scf_type": "total",
+    })
+    psi4.pcm_helper(f"""
+        Units = AU
+        Cavity {{Type = Gepol
+                Area = {pcm_options["weight"]}}}
+        Medium {{Solvertype = {pcm_options["pcm_method"]}
+                Nonequilibrium = {pcm_options["neq"]}
+                Solvent = {pcm_options["solvent"]}}}""")
+
+    if multiplicity != 1:
+        psi4.set_options({
+            'reference': "UHF",
+            'maxiter': max_iter + 500,
+            'soscf': 'true'
+        })
+
+    _, wfn = psi4.energy('SCF', return_wfn=True, molecule=mol)
+    psi4.core.clean()
+    return adcc.backends.import_scf_results(wfn)
diff --git a/adcc/backends/test_backends_polembed.py b/adcc/backends/test_backends_polembed.py
index 9339fc70..2ddf7783 100644
--- a/adcc/backends/test_backends_polembed.py
+++ b/adcc/backends/test_backends_polembed.py
@@ -61,7 +61,7 @@ def template_pe_perturbative_formaldehyde(self, basis, method, backend):
         qc_result = qchem_data[basename]
         pe_options = {"potfile": pe_potentials["fa_6w"]}
         scfres = cached_backend_hf(backend, "formaldehyde", basis,
-                                   env=("pe", pe_options))
+                                   pe_options=pe_options)
         state = adcc.run_adc(scfres, method=method,
                              n_singlets=5, conv_tol=1e-10,
                              environment=["ptlr", "ptss"])
@@ -96,7 +96,7 @@ def template_pe_linear_response_formaldehyde(self, basis, method, backend):
         tm_result = tmole_data[basename]
         pe_options = {"potfile": pe_potentials["fa_6w"]}
         scfres = cached_backend_hf(backend, "formaldehyde", basis,
-                                   env=("pe", pe_options))
+                                   pe_options=pe_options)
         assert_allclose(scfres.energy_scf, tm_result["energy_scf"], atol=1e-8)
 
         matrix = adcc.AdcMatrix(method, scfres)
diff --git a/adcc/backends/testing.py b/adcc/backends/testing.py
index 83adecf9..d41261df 100644
--- a/adcc/backends/testing.py
+++ b/adcc/backends/testing.py
@@ -178,7 +178,7 @@ def operator_import_from_ao_test(scfres, ao_dict, operator="electric_dipole"):
 
 
 def cached_backend_hf(backend, molecule, basis, multiplicity=1, conv_tol=1e-12,
-                      env=(None, None)):
+                      pe_options=None):
     """
     Run the SCF for a backend and a particular test case (if not done)
     and return the result.
@@ -189,21 +189,13 @@ def cached_backend_hf(backend, molecule, basis, multiplicity=1, conv_tol=1e-12,
 
     global __cache_cached_backend_hf
 
-    # so you can only assign one of them.
-    pe_options = None
-    pcm = None
-    if env[0] == "pcm":
-        pcm = env[1]
-    elif env[0] == "pe":
-        pe_options = env[1]
-
     def payload():
         conv_tol_grad = 10 * conv_tol
         hfres = adcc.backends.run_hf(backend, xyz=static_data.xyz[molecule],
                                      basis=basis, conv_tol=conv_tol,
                                      multiplicity=multiplicity,
                                      conv_tol_grad=conv_tol_grad,
-                                     pe_options=pe_options, pcm=pcm)
+                                     pe_options=pe_options)
         return adcc.backends.import_scf_results(hfres)
 
     # For reasons not clear to me (mfh), caching does not work
@@ -211,11 +203,7 @@ def payload():
     if backend == "pyscf":
         return payload()
 
-    # otherwise PE and PCM share the same key
-    if env[0]:
-        key = (backend, molecule, basis, str(multiplicity), env[0])
-    else:
-        key = (backend, molecule, basis, str(multiplicity))
+    key = (backend, molecule, basis, str(multiplicity))
     try:
         return __cache_cached_backend_hf[key]
     except NameError:
diff --git a/adcc/testdata/cache.py b/adcc/testdata/cache.py
index bb9c2cd2..87e5d136 100644
--- a/adcc/testdata/cache.py
+++ b/adcc/testdata/cache.py
@@ -258,4 +258,4 @@ def read_yaml_data(fname):
 
 qchem_data = read_yaml_data("qchem_dump.yml")
 tmole_data = read_yaml_data("tmole_dump.yml")
-psi_data = read_yaml_data("psi_dump.yml")
+psi4_data = read_yaml_data("psi4_dump.yml")
diff --git a/adcc/testdata/generate_psi4_data.py b/adcc/testdata/generate_psi4_data.py
new file mode 100644
index 00000000..99eb9821
--- /dev/null
+++ b/adcc/testdata/generate_psi4_data.py
@@ -0,0 +1,97 @@
+import psi4
+from psi4.driver.procrouting.response.scf_response import tdscf_excitations
+from adcc.testdata import static_data
+import yaml
+import os
+
+
+def run_psi4_tdscf(xyz, basis, charge=0, multiplicity=1,
+                   conv_tol=1e-12, conv_tol_grad=1e-11, max_iter=150,
+                   pcm_options=None):
+    mol = psi4.geometry(f"""
+        {charge} {multiplicity}
+        {xyz}
+        units au
+        symmetry c1
+    """)
+    psi4.set_options({
+        'basis': basis,
+        'scf_type': "pK",
+        'e_convergence': conv_tol,
+        'd_convergence': conv_tol_grad,
+        'maxiter': max_iter,
+        'reference': "RHF",
+        'save_jk': True,
+    })
+    if pcm_options:
+        # think its better to use a dict for the pcm options, because there
+        # are already enough function arguments. Of course one could just
+        # fix the options... would be sufficient for now too.
+        psi4.set_options({'pcm': True, 'pcm_scf_type': "total"})
+        psi4.pcm_helper(f"""
+            Units = AU
+            Cavity {{
+                Type = Gepol
+                Area = {pcm_options["weight"]}
+            }}
+            Medium {{
+                Solvertype = {pcm_options["pcm_method"]}
+                Solvent = {pcm_options["solvent"]}
+                Nonequilibrium = {pcm_options["neq"]}
+            }}
+        """)
+    psi4.core.be_quiet()
+    _, wfn = psi4.energy('scf', return_wfn=True, molecule=mol)
+
+    res = tdscf_excitations(wfn, states=5, triplets="NONE", tda=True,
+                            r_convergence=1e-7)
+
+    # remove cavity files from PSI4 PCM calculations
+    if pcm_options:
+        for cavityfile in os.listdir(os.getcwd()):
+            if cavityfile.startswith(("cavity.off_", "PEDRA.OUT_")):
+                os.remove(cavityfile)
+    return wfn, res
+
+
+def dump_results(molecule, basis, **kwargs):
+    pcm_options = kwargs.get("pcm_options", None)
+    if pcm_options:
+        name = f"{molecule}_{basis}_pcm_adc1"
+    else:
+        name = f"{molecule}_{basis}_adc1"
+
+    geom = static_data.xyz[molecule]
+
+    ret = {"basis": basis,
+           "method": "adc1",
+           "molecule": molecule}
+
+    wfn, res = run_psi4_tdscf(geom, basis, pcm_options=pcm_options)
+
+    ret["energy_scf"] = wfn.energy()
+    # yaml safe_dump doesn't like np.floats
+    ret["excitation_energy"] = [round(float(r["EXCITATION ENERGY"]), 5)
+                                for r in res]
+    ret["osc_strength"] = [round(float(r["OSCILLATOR STRENGTH (LEN)"]), 4)
+                           for r in res]
+    return name, ret
+
+
+def main():
+    basis_set = ["sto-3g", "cc-pvdz"]
+    pcm_options = {"weight": 0.3, "pcm_method": "IEFPCM", "neq": True,
+                   "solvent": "Water"}
+    psi4_results = {}
+    for basis in basis_set:
+        key, ret = dump_results("formaldehyde", basis, pcm_options=pcm_options)
+        psi4_results[key] = ret
+        print(f"Dumped {key}.")
+    # how to print the lists as lists, while keeping the rest of the output as is?
+    # print just to some test file atm.
+    with open("generate_psi4_data_test.yml", 'w') as yamlfile:
+        yaml.safe_dump(psi4_results, yamlfile)
+
+
+if __name__ == "__main__":
+    main()
diff --git a/adcc/testdata/generate_psi4_data_test.yml b/adcc/testdata/generate_psi4_data_test.yml
new file mode 100644
index 00000000..e20bf4fb
--- /dev/null
+++ b/adcc/testdata/generate_psi4_data_test.yml
@@ -0,0 +1,34 @@
+formaldehyde_cc-pvdz_pcm_adc1:
+  basis: cc-pvdz
+  energy_scf: -113.88427067069416
+  excitation_energy:
+  - 0.179
+  - 0.38177
+  - 0.38802
+  - 0.39216
+  - 0.43544
+  method: adc1
+  molecule: formaldehyde
+  osc_strength:
+  - 0.0
+  - 0.0001
+  - 0.2487
+  - 0.2349
+  - 0.0
+formaldehyde_sto-3g_pcm_adc1:
+  basis: sto-3g
+  energy_scf: -112.35458302984613
+  excitation_energy:
+  - 0.16467
+  - 0.3609
+  - 0.46549
+  - 0.50874
+  - 0.69386
+  method: adc1
+  molecule: formaldehyde
+  osc_strength:
+  - 0.0
+  - 0.0107
+  - 0.3434
+  - 0.0
+  - 0.3586
diff --git a/adcc/testdata/psi_dump.yml b/adcc/testdata/psi4_dump.yml
similarity index 78%
rename from adcc/testdata/psi_dump.yml
rename to adcc/testdata/psi4_dump.yml
index e5a730f3..0c19547c 100644
--- a/adcc/testdata/psi_dump.yml
+++ b/adcc/testdata/psi4_dump.yml
@@ -4,9 +4,11 @@ formaldehyde_sto3g_pcm_adc1:
   molecule: formaldehyde
   energy_scf: -112.3545830298461254
   excitation_energy: [0.16467, 0.36090, 0.46549, 0.50874, 0.69386]
+  osc_strength: [0.0, 0.0107, 0.3434, 0.0, 0.3586]
 formaldehyde_ccpvdz_pcm_adc1:
   basis: ccpvdz
   method: adc1
   molecule: formaldehyde
   energy_scf: -113.8842706706941641
-  excitation_energy: [0.17900, 0.38177, 0.38802, 0.39216, 0.43544]
\ No newline at end of file
+  excitation_energy: [0.17900, 0.38177, 0.38802, 0.39216, 0.43544]
+  osc_strength: [0.0, 0.0001, 0.2487, 0.2349, 0.0]
\ No newline at end of file

From 98afca93f54b91eedf2bc885c636ce04b4765248 Mon Sep 17 00:00:00 2001
From: Jonas Leitner <j.leitner@stud.uni-heidelberg.de>
Date: Wed, 21 Jul 2021 10:34:09 +0200
Subject: [PATCH 10/16] fix yaml generator

---
 adcc/backends/test_backends_pcm.py        |  2 +-
 adcc/testdata/generate_psi4_data.py       |  5 ++--
 adcc/testdata/generate_psi4_data_test.yml | 34 -----------------------
 adcc/testdata/psi4_dump.yml               | 16 +++++------
 4 files changed, 12 insertions(+), 45 deletions(-)
 delete mode 100644 adcc/testdata/generate_psi4_data_test.yml

diff --git a/adcc/backends/test_backends_pcm.py b/adcc/backends/test_backends_pcm.py
index 1a49eaeb..e3edc9d4 100644
--- a/adcc/backends/test_backends_pcm.py
+++ b/adcc/backends/test_backends_pcm.py
@@ -18,7 +18,7 @@
 # remove pyscf until implemented
 backends = [b for b in adcc.backends.available()
             if b not in ["molsturm", "veloxchem", "pyscf"]]
-basissets = ["sto3g", "ccpvdz"]
+basissets = ["sto-3g", "cc-pvdz"]
 methods = ["adc1", "adc2", "adc3"]
 
 
diff --git a/adcc/testdata/generate_psi4_data.py b/adcc/testdata/generate_psi4_data.py
index 99eb9821..a5f4f9cd 100644
--- a/adcc/testdata/generate_psi4_data.py
+++ b/adcc/testdata/generate_psi4_data.py
@@ -89,8 +89,9 @@ def main():
         print(f"Dumped {key}.")
     # how to print the lists as lists, while keeping the rest of the output as is?
     # print just to some test file atm.
-    with open("generate_psi4_data_test.yml", 'w') as yamlfile:
-        yaml.safe_dump(psi4_results, yamlfile)
+    with open("psi4_dump.yml", 'w') as yamlfile:
+        yaml.safe_dump(psi4_results, yamlfile,
+                       sort_keys=False, default_flow_style=None)
 
 
 if __name__ == "__main__":
diff --git a/adcc/testdata/generate_psi4_data_test.yml b/adcc/testdata/generate_psi4_data_test.yml
deleted file mode 100644
index e20bf4fb..00000000
--- a/adcc/testdata/generate_psi4_data_test.yml
+++ /dev/null
@@ -1,34 +0,0 @@
-formaldehyde_cc-pvdz_pcm_adc1:
-  basis: cc-pvdz
-  energy_scf: -113.88427067069416
-  excitation_energy:
-  - 0.179
-  - 0.38177
-  - 0.38802
-  - 0.39216
-  - 0.43544
-  method: adc1
-  molecule: formaldehyde
-  osc_strength:
-  - 0.0
-  - 0.0001
-  - 0.2487
-  - 0.2349
-  - 0.0
-formaldehyde_sto-3g_pcm_adc1:
-  basis: sto-3g
-  energy_scf: -112.35458302984613
-  excitation_energy:
-  - 0.16467
-  - 0.3609
-  - 0.46549
-  - 0.50874
-  - 0.69386
-  method: adc1
-  molecule: formaldehyde
-  osc_strength:
-  - 0.0
-  - 0.0107
-  - 0.3434
-  - 0.0
-  - 0.3586
diff --git a/adcc/testdata/psi4_dump.yml b/adcc/testdata/psi4_dump.yml
index 0c19547c..db675255 100644
--- a/adcc/testdata/psi4_dump.yml
+++ b/adcc/testdata/psi4_dump.yml
@@ -1,14 +1,14 @@
-formaldehyde_sto3g_pcm_adc1:
+formaldehyde_sto-3g_pcm_adc1:
   basis: sto-3g
   method: adc1
   molecule: formaldehyde
-  energy_scf: -112.3545830298461254
-  excitation_energy: [0.16467, 0.36090, 0.46549, 0.50874, 0.69386]
+  energy_scf: -112.35458302984613
+  excitation_energy: [0.16467, 0.3609, 0.46549, 0.50874, 0.69386]
   osc_strength: [0.0, 0.0107, 0.3434, 0.0, 0.3586]
-formaldehyde_ccpvdz_pcm_adc1:
-  basis: ccpvdz
+formaldehyde_cc-pvdz_pcm_adc1:
+  basis: cc-pvdz
   method: adc1
   molecule: formaldehyde
-  energy_scf: -113.8842706706941641
-  excitation_energy: [0.17900, 0.38177, 0.38802, 0.39216, 0.43544]
-  osc_strength: [0.0, 0.0001, 0.2487, 0.2349, 0.0]
\ No newline at end of file
+  energy_scf: -113.88427067069416
+  excitation_energy: [0.179, 0.38177, 0.38802, 0.39216, 0.43544]
+  osc_strength: [0.0, 0.0001, 0.2487, 0.2349, 0.0]

From 49ec8cc25472c6150c2f8b2bfd45c67e60aaad55 Mon Sep 17 00:00:00 2001
From: Jonas Leitner <j.leitner@stud.uni-heidelberg.de>
Date: Mon, 9 Aug 2021 12:05:25 +0200
Subject: [PATCH 11/16] cleanup + psi4 specific tests

---
 adcc/backends/test_backends_pcm.py  | 19 ++++++++-----------
 adcc/testdata/generate_psi4_data.py |  3 +--
 2 files changed, 9 insertions(+), 13 deletions(-)

diff --git a/adcc/backends/test_backends_pcm.py b/adcc/backends/test_backends_pcm.py
index e3edc9d4..3e876535 100644
--- a/adcc/backends/test_backends_pcm.py
+++ b/adcc/backends/test_backends_pcm.py
@@ -3,7 +3,6 @@
 import itertools
 import adcc
 import adcc.backends
-import psi4
 import os
 
 from adcc.misc import expand_test_templates
@@ -23,9 +22,9 @@
 
 
 @pytest.mark.skipif(len(backends) == 0, reason="No backend found.")
-@expand_test_templates(list(itertools.product(basissets, methods, backends)))
-class TestPCM(unittest.TestCase):
-    def template_pcm_linear_response_formaldehyde(self, basis, method, backend):
+@expand_test_templates(list(itertools.product(basissets, methods)))
+class TestPCM_psi4(unittest.TestCase):
+    def template_pcm_linear_response_formaldehyde(self, basis, method):
         if method != "adc1":
             pytest.skip("Reference only exists for adc1.")
         basename = f"formaldehyde_{basis}_pcm_{method}"
@@ -34,14 +33,11 @@ def template_pcm_linear_response_formaldehyde(self, basis, method, backend):
         # so you can adapt some options more easily if needed
         psi4_pcm_options = {"weight": 0.3, "pcm_method": "IEFPCM", "neq": True,
                             "solvent": "Water"}
-        pcm_options = {"psi4": psi4_pcm_options}
 
-        hf_name = f"{backend}_run_pcm_hf"
-        hf_function = globals()[hf_name]
-        scfres = hf_function(static_data.xyz["formaldehyde"], basis, charge=0,
-                             multiplicity=1, conv_tol=1e-12,
-                             conv_tol_grad=1e-11, max_iter=150,
-                             pcm_options=pcm_options[backend])
+        scfres = psi4_run_pcm_hf(static_data.xyz["formaldehyde"], basis, charge=0,
+                                 multiplicity=1, conv_tol=1e-12,
+                                 conv_tol_grad=1e-11, max_iter=150,
+                                 pcm_options=psi4_pcm_options)
 
         assert_allclose(scfres.energy_scf, psi4_result["energy_scf"], atol=1e-8)
 
@@ -93,6 +89,7 @@ def psi4_run_pcm_hf(xyz, basis, charge=0, multiplicity=1, conv_tol=1e-12,
                     conv_tol_grad=1e-11, max_iter=150, pcm_options=None):
 
     # add something to check if pcm_options is defined and a dict?
+    import psi4
 
     basissets = {
         "sto3g": "sto-3g",
diff --git a/adcc/testdata/generate_psi4_data.py b/adcc/testdata/generate_psi4_data.py
index a5f4f9cd..1ef507a9 100644
--- a/adcc/testdata/generate_psi4_data.py
+++ b/adcc/testdata/generate_psi4_data.py
@@ -87,8 +87,7 @@ def main():
         key, ret = dump_results("formaldehyde", basis, pcm_options=pcm_options)
         psi4_results[key] = ret
         print(f"Dumped {key}.")
-    # how to print the lists as lists, while keeping the rest of the output as is?
-    # print just to some test file atm.
+
     with open("psi4_dump.yml", 'w') as yamlfile:
         yaml.safe_dump(psi4_results, yamlfile,
                        sort_keys=False, default_flow_style=None)

From 1326ea4fa924e88eaa7c5cd1a167277db422a85a Mon Sep 17 00:00:00 2001
From: Jonas Leitner <j.leitner@stud.uni-heidelberg.de>
Date: Mon, 9 Aug 2021 16:15:19 +0200
Subject: [PATCH 12/16] some more cleanup

---
 adcc/backends/test_backends_pcm.py | 10 ++++------
 1 file changed, 4 insertions(+), 6 deletions(-)

diff --git a/adcc/backends/test_backends_pcm.py b/adcc/backends/test_backends_pcm.py
index 3e876535..8eb0d81b 100644
--- a/adcc/backends/test_backends_pcm.py
+++ b/adcc/backends/test_backends_pcm.py
@@ -14,17 +14,15 @@
 from adcc.adc_pp.environment import block_ph_ph_0_pcm
 from adcc.AdcMatrix import AdcExtraTerm
 
-# remove pyscf until implemented
-backends = [b for b in adcc.backends.available()
-            if b not in ["molsturm", "veloxchem", "pyscf"]]
 basissets = ["sto-3g", "cc-pvdz"]
 methods = ["adc1", "adc2", "adc3"]
 
 
-@pytest.mark.skipif(len(backends) == 0, reason="No backend found.")
+@pytest.mark.skipif("psi4" not in adcc.backends.available(),
+                    reason="Psi4 backend not found.")
 @expand_test_templates(list(itertools.product(basissets, methods)))
-class TestPCM_psi4(unittest.TestCase):
-    def template_pcm_linear_response_formaldehyde(self, basis, method):
+class TestPCM(unittest.TestCase):
+    def template_pcm_psi4_linear_response_formaldehyde(self, basis, method):
         if method != "adc1":
             pytest.skip("Reference only exists for adc1.")
         basename = f"formaldehyde_{basis}_pcm_{method}"

From 9c1baffa15734c86f16bba92292847bdb8c7af2d Mon Sep 17 00:00:00 2001
From: Jonas Leitner <j.leitner@stud.uni-heidelberg.de>
Date: Sun, 7 Nov 2021 19:21:47 +0100
Subject: [PATCH 13/16] first try implementing ptlr pcm

---
 adcc/backends/lr_pcm_test.py | 53 ++++++++++++++++++++++++++++++++++++
 adcc/backends/psi4.py        | 16 +++++++++++
 2 files changed, 69 insertions(+)
 create mode 100644 adcc/backends/lr_pcm_test.py

diff --git a/adcc/backends/lr_pcm_test.py b/adcc/backends/lr_pcm_test.py
new file mode 100644
index 00000000..712d1229
--- /dev/null
+++ b/adcc/backends/lr_pcm_test.py
@@ -0,0 +1,53 @@
+import psi4
+import adcc
+
+
+mol = psi4.geometry("""
+0 1
+    C 2.0092420208996 3.8300915804899 0.8199294419789
+    O 2.1078857690998 2.0406638776593 2.1812021228452
+    H 2.0682421748693 5.7438044586615 1.5798996515014
+    H 1.8588483602149 3.6361694243085 -1.2192956060942
+units au
+symmetry c1
+""")
+psi4.set_options({
+    'basis': "6-31g",
+    'scf_type': "pK",
+    'e_convergence': 1e-12,
+    'd_convergence': 1e-11,
+    'reference': "RHF",
+    'pcm': True,
+    'pcm_scf_type': "total"
+    }
+)
+psi4.pcm_helper(""" 
+    Units = AU
+    Cavity {
+        Type = Gepol
+        Area = 0.3
+        Radiiset = Bondi
+        Scaling = True
+    }
+    Medium {
+        Solvertype = IEFPCM
+        Solvent = Water
+        Nonequilibrium = True
+    }
+""")
+
+# psi4.core.be_quiet()
+psi4.core.set_output_file('adcc.out', False)
+scf_e, scf_wfn = psi4.energy('scf', return_wfn=True)
+state = adcc.run_adc(scf_wfn, method="adc1", n_singlets=3,
+                     conv_tol=1e-7, environment="ptlr")
+
+print(state.describe())
+
+for i, energy in enumerate(state.excitation_energy):
+    print(f"\nstate: {i+1}")
+    print("uncorrected {:.7f}, corrected {:.7f}, diff(correction): {:.7f}"
+          .format(state.excitation_energy_uncorrected[i], energy,
+                  state.excitation_energy_uncorrected[i] - energy
+                  )
+          )
diff --git a/adcc/backends/psi4.py b/adcc/backends/psi4.py
index b6ec0ed5..4a55140f 100644
--- a/adcc/backends/psi4.py
+++ b/adcc/backends/psi4.py
@@ -116,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())
+        V_pcm = psi4.core.PCM.compute_V(self.wfn.get_PCM(), psi_dm).to_array()
+        # either tranform both to np array and use np.einsum
+        # or tranform V_pcm to libadcc.Tensor
+        # numpy seems to be easier for now
+        return np.einsum("uv,uv->", dm.to_ndarray(), V_pcm)
+
     @property
     def excitation_energy_corrections(self):
         ret = []
@@ -131,6 +139,14 @@ def excitation_energy_corrections(self):
                                             elec_only=True)
             )
             ret.extend([ptlr, ptss])
+        if self.environment == "pcm":
+            # for the correction the correct transition_dm has to be used?
+            # <I|V|0><J|Qf|0>
+            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

From 3d17c98145f6513aed361ca6f2cc2c256061c4bf Mon Sep 17 00:00:00 2001
From: Jonas Leitner <j.leitner@stud.uni-heidelberg.de>
Date: Thu, 25 Nov 2021 11:35:35 +0100
Subject: [PATCH 14/16] add some tests - 1) versus LR - 2) against itself

---
 adcc/backends/lr_pcm_test.py        | 53 -------------------------
 adcc/backends/psi4.py               |  7 +---
 adcc/backends/test_backends_pcm.py  | 60 +++++++++++++++++++++--------
 adcc/testdata/generate_psi4_data.py | 29 ++++++++++----
 adcc/testdata/psi4_dump.yml         | 10 +++--
 5 files changed, 72 insertions(+), 87 deletions(-)
 delete mode 100644 adcc/backends/lr_pcm_test.py

diff --git a/adcc/backends/lr_pcm_test.py b/adcc/backends/lr_pcm_test.py
deleted file mode 100644
index 712d1229..00000000
--- a/adcc/backends/lr_pcm_test.py
+++ /dev/null
@@ -1,53 +0,0 @@
-import psi4
-import adcc
-
-
-mol = psi4.geometry("""
-0 1
-    C 2.0092420208996 3.8300915804899 0.8199294419789
-    O 2.1078857690998 2.0406638776593 2.1812021228452
-    H 2.0682421748693 5.7438044586615 1.5798996515014
-    H 1.8588483602149 3.6361694243085 -1.2192956060942
-units au
-symmetry c1
-""")
-psi4.set_options({
-    'basis': "6-31g",
-    'scf_type': "pK",
-    'e_convergence': 1e-12,
-    'd_convergence': 1e-11,
-    'reference': "RHF",
-    'pcm': True,
-    'pcm_scf_type': "total"
-    }
-)
-psi4.pcm_helper(""" 
-    Units = AU
-    Cavity {
-        Type = Gepol
-        Area = 0.3
-        Radiiset = Bondi
-        Scaling = True
-    }
-    Medium {
-        Solvertype = IEFPCM
-        Solvent = Water
-        Nonequilibrium = True
-    }
-""")
-
-# psi4.core.be_quiet()
-psi4.core.set_output_file('adcc.out', False)
-scf_e, scf_wfn = psi4.energy('scf', return_wfn=True)
-state = adcc.run_adc(scf_wfn, method="adc1", n_singlets=3,
-                     conv_tol=1e-7, environment="ptlr")
-
-print(state.describe())
-
-for i, energy in enumerate(state.excitation_energy):
-    print(f"\nstate: {i+1}")
-    print("uncorrected {:.7f}, corrected {:.7f}, diff(correction): {:.7f}"
-          .format(state.excitation_energy_uncorrected[i], energy,
-                  state.excitation_energy_uncorrected[i] - energy
-                  )
-          )
diff --git a/adcc/backends/psi4.py b/adcc/backends/psi4.py
index 4a55140f..4d72844f 100644
--- a/adcc/backends/psi4.py
+++ b/adcc/backends/psi4.py
@@ -118,10 +118,9 @@ def pe_energy(self, dm, elec_only=True):
 
     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 contribution is obtained
         V_pcm = psi4.core.PCM.compute_V(self.wfn.get_PCM(), psi_dm).to_array()
-        # either tranform both to np array and use np.einsum
-        # or tranform V_pcm to libadcc.Tensor
-        # numpy seems to be easier for now
         return np.einsum("uv,uv->", dm.to_ndarray(), V_pcm)
 
     @property
@@ -140,8 +139,6 @@ def excitation_energy_corrections(self):
             )
             ret.extend([ptlr, ptss])
         if self.environment == "pcm":
-            # for the correction the correct transition_dm has to be used?
-            # <I|V|0><J|Qf|0>
             ptlr = EnergyCorrection(
                 "pcm_ptlr_correction",
                 lambda view: self.pcm_energy(view.transition_dm_ao)
diff --git a/adcc/backends/test_backends_pcm.py b/adcc/backends/test_backends_pcm.py
index 8eb0d81b..1e004ada 100644
--- a/adcc/backends/test_backends_pcm.py
+++ b/adcc/backends/test_backends_pcm.py
@@ -22,13 +22,47 @@
                     reason="Psi4 backend not found.")
 @expand_test_templates(list(itertools.product(basissets, methods)))
 class TestPCM(unittest.TestCase):
+    def template_pcm_psi4_ptlr_formaldehyde(self, basis, method):
+        if method != "adc1":
+            pytest.skip("Data only available for adc1.")
+        basename = f"formaldehyde_{basis}_pcm_{method}"
+        psi4_pcm_options = {"weight": 0.3, "pcm_method": "IEFPCM", "neq": True,
+                            "solvent": "Water"}
+        scfres = psi4_run_pcm_hf(static_data.xyz["formaldehyde"], basis, charge=0,
+                                 multiplicity=1, conv_tol=1e-12,
+                                 conv_tol_grad=1e-11, max_iter=150,
+                                 pcm_options=psi4_pcm_options)
+
+        psi4_result = psi4_data[basename]
+        assert_allclose(scfres.energy_scf, psi4_result["energy_scf"], atol=1e-8)
+
+        state = adcc.run_adc(scfres, method=method, n_singlets=5,
+                             conv_tol=1e-7, environment="ptlr")
+        # is it possible to perform a PTE calculation with PSI4? Don't find anything
+        # environment related in the tdscf code
+        # Add comparison of uncorrected zeroth order excitation energies.
+
+        # compare ptLR result to LR data
+        # usually agrees up to third significant figure for energies in Hartree
+        # However formaldehyde seems to show a slightly larger discrepancy
+        assert_allclose(state.excitation_energy,
+                        psi4_result["lr_excitation_energy"], atol=5 * 1e-3)
+
+        # Consistency check with values obtained with ADCc
+        assert_allclose(state.excitation_energy,
+                        psi4_result["ptlr_adcc_excitation_energy"], atol=1e-5)
+
+        # remove cavity files from PSI4 PCM calculations
+        for cavityfile in os.listdir(os.getcwd()):
+            if cavityfile.startswith(("cavity.off_", "PEDRA.OUT_")):
+                os.remove(cavityfile)
+
     def template_pcm_psi4_linear_response_formaldehyde(self, basis, method):
         if method != "adc1":
             pytest.skip("Reference only exists for adc1.")
         basename = f"formaldehyde_{basis}_pcm_{method}"
         psi4_result = psi4_data[basename]
 
-        # so you can adapt some options more easily if needed
         psi4_pcm_options = {"weight": 0.3, "pcm_method": "IEFPCM", "neq": True,
                             "solvent": "Water"}
 
@@ -49,14 +83,14 @@ def template_pcm_psi4_linear_response_formaldehyde(self, basis, method):
                              environment=False)
         assert_allclose(
             state.excitation_energy_uncorrected,
-            psi4_result["excitation_energy"],
+            psi4_result["lr_excitation_energy"],
             atol=1e-5
         )
 
         state_cis = adcc.ExcitedStates(state, property_method="adc0")
         assert_allclose(
             state_cis.oscillator_strength,
-            psi4_result["osc_strength"], atol=1e-3
+            psi4_result["lr_osc_strength"], atol=1e-3
         )
 
         # invalid combination
@@ -73,7 +107,7 @@ def template_pcm_psi4_linear_response_formaldehyde(self, basis, method):
                              conv_tol=1e-7, environment="linear_response")
         assert_allclose(
             state.excitation_energy_uncorrected,
-            psi4_result["excitation_energy"],
+            psi4_result["lr_excitation_energy"],
             atol=1e-5
         )
 
@@ -85,16 +119,8 @@ def template_pcm_psi4_linear_response_formaldehyde(self, basis, method):
 
 def psi4_run_pcm_hf(xyz, basis, charge=0, multiplicity=1, conv_tol=1e-12,
                     conv_tol_grad=1e-11, max_iter=150, pcm_options=None):
-
-    # add something to check if pcm_options is defined and a dict?
     import psi4
 
-    basissets = {
-        "sto3g": "sto-3g",
-        "def2tzvp": "def2-tzvp",
-        "ccpvdz": "cc-pvdz",
-    }
-
     # needed for PE and PCM tests
     psi4.core.clean_options()
     mol = psi4.geometry(f"""
@@ -108,7 +134,7 @@ def psi4_run_pcm_hf(xyz, basis, charge=0, multiplicity=1, conv_tol=1e-12,
 
     psi4.core.be_quiet()
     psi4.set_options({
-        'basis': basissets.get(basis, basis),
+        'basis': basis,
         'scf_type': 'pk',
         'e_convergence': conv_tol,
         'd_convergence': conv_tol_grad,
@@ -120,10 +146,10 @@ def psi4_run_pcm_hf(xyz, basis, charge=0, multiplicity=1, conv_tol=1e-12,
     psi4.pcm_helper(f"""
         Units = AU
         Cavity {{Type = Gepol
-                Area = {pcm_options["weight"]}}}
-        Medium {{Solvertype = {pcm_options["pcm_method"]}
-                Nonequilibrium = {pcm_options["neq"]}
-                Solvent = {pcm_options["solvent"]}}}""")
+                Area = {pcm_options.get("weight", 0.3)}}}
+        Medium {{Solvertype = {pcm_options.get("pcm_method", "IEFPCM")}
+                Nonequilibrium = {pcm_options.get("neq", True)}
+                Solvent = {pcm_options.get("solvent", "Water")}}}""")
 
     if multiplicity != 1:
         psi4.set_options({
diff --git a/adcc/testdata/generate_psi4_data.py b/adcc/testdata/generate_psi4_data.py
index 1ef507a9..835e561d 100644
--- a/adcc/testdata/generate_psi4_data.py
+++ b/adcc/testdata/generate_psi4_data.py
@@ -1,4 +1,5 @@
 import psi4
+import adcc
 from psi4.driver.procrouting.response.scf_response import tdscf_excitations
 from adcc.testdata import static_data
 import yaml
@@ -32,12 +33,12 @@ def run_psi4_tdscf(xyz, basis, charge=0, multiplicity=1,
             Units = AU
             Cavity {{
                 Type = Gepol
-                Area = {pcm_options["weight"]}
+                Area = {pcm_options.get("weight", 0.3)}
             }}
             Medium {{
-                Solvertype = {pcm_options["pcm_method"]}
-                Solvent = {pcm_options["solvent"]}
-                Nonequilibrium = {pcm_options["neq"]}
+                Solvertype = {pcm_options.get("pcm_method", "IEFPCM")}
+                Solvent = {pcm_options.get("solvent", "Water")}
+                Nonequilibrium = {pcm_options.get("neq", True)}
             }}
         """)
     psi4.core.be_quiet()
@@ -54,6 +55,12 @@ def run_psi4_tdscf(xyz, basis, charge=0, multiplicity=1,
     return wfn, res
 
 
+def run_adcc_ptlr(wfn):
+    state = adcc.run_adc(wfn, method="adc1", n_singlets=5,
+                         conv_tol=1e-7, environment="ptlr")
+    return state
+
+
 def dump_results(molecule, basis, **kwargs):
     pcm_options = kwargs.get("pcm_options", None)
     if pcm_options:
@@ -71,10 +78,16 @@ def dump_results(molecule, basis, **kwargs):
 
     ret["energy_scf"] = wfn.energy()
     # yaml safe_dump doesn't like np.floats
-    ret["excitation_energy"] = [round(float(r["EXCITATION ENERGY"]), 5)
-                                for r in res]
-    ret["osc_strength"] = [round(float(r["OSCILLATOR STRENGTH (LEN)"]), 4)
-                           for r in res]
+    ret["lr_excitation_energy"] = [round(float(r["EXCITATION ENERGY"]), 5)
+                                   for r in res]
+    ret["lr_osc_strength"] = [round(float(r["OSCILLATOR STRENGTH (LEN)"]), 4)
+                              for r in res]
+
+    if pcm_options:
+        state = run_adcc_ptlr(wfn)
+        ret["ptlr_adcc_excitation_energy"] = [
+            round(float(e), 5) for e in state.excitation_energy
+        ]
     return name, ret
 
 
diff --git a/adcc/testdata/psi4_dump.yml b/adcc/testdata/psi4_dump.yml
index db675255..108a050f 100644
--- a/adcc/testdata/psi4_dump.yml
+++ b/adcc/testdata/psi4_dump.yml
@@ -3,12 +3,14 @@ formaldehyde_sto-3g_pcm_adc1:
   method: adc1
   molecule: formaldehyde
   energy_scf: -112.35458302984613
-  excitation_energy: [0.16467, 0.3609, 0.46549, 0.50874, 0.69386]
-  osc_strength: [0.0, 0.0107, 0.3434, 0.0, 0.3586]
+  lr_excitation_energy: [0.16467, 0.3609, 0.46549, 0.50874, 0.69386]
+  lr_osc_strength: [0.0, 0.0107, 0.3434, 0.0, 0.3586]
+  ptlr_adcc_excitation_energy: [0.16468, 0.36093, 0.46758, 0.50874, 0.69447]
 formaldehyde_cc-pvdz_pcm_adc1:
   basis: cc-pvdz
   method: adc1
   molecule: formaldehyde
   energy_scf: -113.88427067069416
-  excitation_energy: [0.179, 0.38177, 0.38802, 0.39216, 0.43544]
-  osc_strength: [0.0, 0.0001, 0.2487, 0.2349, 0.0]
+  lr_excitation_energy: [0.179, 0.38177, 0.38802, 0.39216, 0.43544]
+  lr_osc_strength: [0.0, 0.0001, 0.2487, 0.2349, 0.0]
+  ptlr_adcc_excitation_energy: [0.17902, 0.38185, 0.38958, 0.39259, 0.43546]

From 1a3627c0758fc527ba1736e8ae3e6823d945d63d Mon Sep 17 00:00:00 2001
From: Jonas Leitner <j.leitner@stud.uni-heidelberg.de>
Date: Mon, 29 Nov 2021 19:51:11 +0100
Subject: [PATCH 15/16] increased printout in yaml and some cleanup

---
 adcc/backends/test_backends_pcm.py  |  7 +------
 adcc/testdata/generate_psi4_data.py |  6 +++---
 adcc/testdata/psi4_dump.yml         | 14 ++++++++------
 3 files changed, 12 insertions(+), 15 deletions(-)

diff --git a/adcc/backends/test_backends_pcm.py b/adcc/backends/test_backends_pcm.py
index 1e004ada..b51df88f 100644
--- a/adcc/backends/test_backends_pcm.py
+++ b/adcc/backends/test_backends_pcm.py
@@ -38,19 +38,14 @@ def template_pcm_psi4_ptlr_formaldehyde(self, basis, method):
 
         state = adcc.run_adc(scfres, method=method, n_singlets=5,
                              conv_tol=1e-7, environment="ptlr")
-        # is it possible to perform a PTE calculation with PSI4? Don't find anything
-        # environment related in the tdscf code
-        # Add comparison of uncorrected zeroth order excitation energies.
 
         # compare ptLR result to LR data
-        # usually agrees up to third significant figure for energies in Hartree
-        # However formaldehyde seems to show a slightly larger discrepancy
         assert_allclose(state.excitation_energy,
                         psi4_result["lr_excitation_energy"], atol=5 * 1e-3)
 
         # Consistency check with values obtained with ADCc
         assert_allclose(state.excitation_energy,
-                        psi4_result["ptlr_adcc_excitation_energy"], atol=1e-5)
+                        psi4_result["ptlr_adcc_excitation_energy"], atol=1e-6)
 
         # remove cavity files from PSI4 PCM calculations
         for cavityfile in os.listdir(os.getcwd()):
diff --git a/adcc/testdata/generate_psi4_data.py b/adcc/testdata/generate_psi4_data.py
index 835e561d..12d512f7 100644
--- a/adcc/testdata/generate_psi4_data.py
+++ b/adcc/testdata/generate_psi4_data.py
@@ -78,15 +78,15 @@ def dump_results(molecule, basis, **kwargs):
 
     ret["energy_scf"] = wfn.energy()
     # yaml safe_dump doesn't like np.floats
-    ret["lr_excitation_energy"] = [round(float(r["EXCITATION ENERGY"]), 5)
+    ret["lr_excitation_energy"] = [round(float(r["EXCITATION ENERGY"]), 9)
                                    for r in res]
-    ret["lr_osc_strength"] = [round(float(r["OSCILLATOR STRENGTH (LEN)"]), 4)
+    ret["lr_osc_strength"] = [round(float(r["OSCILLATOR STRENGTH (LEN)"]), 5)
                               for r in res]
 
     if pcm_options:
         state = run_adcc_ptlr(wfn)
         ret["ptlr_adcc_excitation_energy"] = [
-            round(float(e), 5) for e in state.excitation_energy
+            round(float(e), 9) for e in state.excitation_energy
         ]
     return name, ret
 
diff --git a/adcc/testdata/psi4_dump.yml b/adcc/testdata/psi4_dump.yml
index 108a050f..fb506bfc 100644
--- a/adcc/testdata/psi4_dump.yml
+++ b/adcc/testdata/psi4_dump.yml
@@ -3,14 +3,16 @@ formaldehyde_sto-3g_pcm_adc1:
   method: adc1
   molecule: formaldehyde
   energy_scf: -112.35458302984613
-  lr_excitation_energy: [0.16467, 0.3609, 0.46549, 0.50874, 0.69386]
-  lr_osc_strength: [0.0, 0.0107, 0.3434, 0.0, 0.3586]
-  ptlr_adcc_excitation_energy: [0.16468, 0.36093, 0.46758, 0.50874, 0.69447]
+  lr_excitation_energy: [0.164670987, 0.360903879, 0.465492571, 0.50874055, 0.693857844]
+  lr_osc_strength: [0.0, 0.01075, 0.34344, 0.0, 0.35858]
+  ptlr_adcc_excitation_energy: [0.164676366, 0.360930626, 0.467575499, 0.508742077,
+    0.694474733]
 formaldehyde_cc-pvdz_pcm_adc1:
   basis: cc-pvdz
   method: adc1
   molecule: formaldehyde
   energy_scf: -113.88427067069416
-  lr_excitation_energy: [0.179, 0.38177, 0.38802, 0.39216, 0.43544]
-  lr_osc_strength: [0.0, 0.0001, 0.2487, 0.2349, 0.0]
-  ptlr_adcc_excitation_energy: [0.17902, 0.38185, 0.38958, 0.39259, 0.43546]
+  lr_excitation_energy: [0.178995043, 0.38176984, 0.388019451, 0.392159876, 0.435439339]
+  lr_osc_strength: [0.0, 7.0e-05, 0.24871, 0.23489, 0.0]
+  ptlr_adcc_excitation_energy: [0.179023803, 0.381851509, 0.389582948, 0.392594025,
+    0.435463213]

From 7a0e0e7f3722769255e26bcd4992d5f2cfed06b2 Mon Sep 17 00:00:00 2001
From: Jonas Leitner <j.leitner@stud.uni-heidelberg.de>
Date: Thu, 20 Jan 2022 08:49:02 +0100
Subject: [PATCH 16/16] add LR and ptLR for pyscf + tests

---
 adcc/backends/pyscf.py               |  39 ++++++++-
 adcc/backends/test_backends_pcm.py   | 118 +++++++++++++++++++--------
 adcc/testdata/cache.py               |   1 +
 adcc/testdata/generate_psi4_data.py  |   5 +-
 adcc/testdata/generate_pyscf_data.py |  94 +++++++++++++++++++++
 adcc/testdata/pyscf_dump.yml         |  18 ++++
 6 files changed, 234 insertions(+), 41 deletions(-)
 create mode 100644 adcc/testdata/generate_pyscf_data.py
 create mode 100644 adcc/testdata/pyscf_dump.yml

diff --git a/adcc/backends/pyscf.py b/adcc/backends/pyscf.py
index 90128750..f3d947e4 100644
--- a/adcc/backends/pyscf.py
+++ b/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):
diff --git a/adcc/backends/test_backends_pcm.py b/adcc/backends/test_backends_pcm.py
index b51df88f..f98090b2 100644
--- a/adcc/backends/test_backends_pcm.py
+++ b/adcc/backends/test_backends_pcm.py
@@ -7,66 +7,82 @@
 
 from adcc.misc import expand_test_templates
 from adcc.testdata import static_data
-from adcc.testdata.cache import psi4_data
+from adcc.testdata.cache import psi4_data, pyscf_data
 from numpy.testing import assert_allclose
 from adcc.exceptions import InputError
 
 from adcc.adc_pp.environment import block_ph_ph_0_pcm
 from adcc.AdcMatrix import AdcExtraTerm
 
+backends = [b for b in adcc.backends.available()
+            if b in ["psi4", "pyscf"]]
 basissets = ["sto-3g", "cc-pvdz"]
 methods = ["adc1", "adc2", "adc3"]
 
 
-@pytest.mark.skipif("psi4" not in adcc.backends.available(),
-                    reason="Psi4 backend not found.")
-@expand_test_templates(list(itertools.product(basissets, methods)))
+@pytest.mark.skipif(len(backends) == 0, reason="Psi4 and Pyscf backends not found.")
+@expand_test_templates(list(itertools.product(basissets, methods, backends)))
 class TestPCM(unittest.TestCase):
-    def template_pcm_psi4_ptlr_formaldehyde(self, basis, method):
+    def template_pcm_ptlr_formaldehyde(self, basis, method, backend):
         if method != "adc1":
             pytest.skip("Data only available for adc1.")
+
+        if backend == "psi4":
+            data = psi4_data
+            run_hf = psi4_run_pcm_hf
+            pcm_options = {"weight": 0.3, "pcm_method": "IEFPCM", "neq": True,
+                           "solvent": "Water"}
+        elif backend == "pyscf":
+            data = pyscf_data
+            run_hf = pyscf_run_pcm_hf
+            pcm_options = {"eps": 78.3553, "eps_opt": 1.78}
         basename = f"formaldehyde_{basis}_pcm_{method}"
-        psi4_pcm_options = {"weight": 0.3, "pcm_method": "IEFPCM", "neq": True,
-                            "solvent": "Water"}
-        scfres = psi4_run_pcm_hf(static_data.xyz["formaldehyde"], basis, charge=0,
-                                 multiplicity=1, conv_tol=1e-12,
-                                 conv_tol_grad=1e-11, max_iter=150,
-                                 pcm_options=psi4_pcm_options)
+        result = data[basename]
+
+        scfres = run_hf(static_data.xyz["formaldehyde"], basis, charge=0,
+                        multiplicity=1, conv_tol=1e-12, conv_tol_grad=1e-11,
+                        max_iter=150, pcm_options=pcm_options)
 
-        psi4_result = psi4_data[basename]
-        assert_allclose(scfres.energy_scf, psi4_result["energy_scf"], atol=1e-8)
+        assert_allclose(scfres.energy_scf, result["energy_scf"], atol=1e-8)
 
         state = adcc.run_adc(scfres, method=method, n_singlets=5,
                              conv_tol=1e-7, environment="ptlr")
 
         # compare ptLR result to LR data
         assert_allclose(state.excitation_energy,
-                        psi4_result["lr_excitation_energy"], atol=5 * 1e-3)
+                        result["lr_excitation_energy"], atol=5 * 1e-3)
 
         # Consistency check with values obtained with ADCc
         assert_allclose(state.excitation_energy,
-                        psi4_result["ptlr_adcc_excitation_energy"], atol=1e-6)
+                        result["ptlr_adcc_excitation_energy"], atol=1e-6)
 
         # remove cavity files from PSI4 PCM calculations
-        for cavityfile in os.listdir(os.getcwd()):
-            if cavityfile.startswith(("cavity.off_", "PEDRA.OUT_")):
-                os.remove(cavityfile)
+        if backend == "psi4":
+            for cavityfile in os.listdir(os.getcwd()):
+                if cavityfile.startswith(("cavity.off_", "PEDRA.OUT_")):
+                    os.remove(cavityfile)
 
-    def template_pcm_psi4_linear_response_formaldehyde(self, basis, method):
+    def template_pcm_linear_response_formaldehyde(self, basis, method, backend):
         if method != "adc1":
             pytest.skip("Reference only exists for adc1.")
-        basename = f"formaldehyde_{basis}_pcm_{method}"
-        psi4_result = psi4_data[basename]
 
-        psi4_pcm_options = {"weight": 0.3, "pcm_method": "IEFPCM", "neq": True,
-                            "solvent": "Water"}
+        if backend == "psi4":
+            data = psi4_data
+            run_hf = psi4_run_pcm_hf
+            pcm_options = {"weight": 0.3, "pcm_method": "IEFPCM", "neq": True,
+                           "solvent": "Water"}
+        elif backend == "pyscf":
+            data = pyscf_data
+            run_hf = pyscf_run_pcm_hf
+            pcm_options = {"eps": 78.3553, "eps_opt": 1.78}
+        basename = f"formaldehyde_{basis}_pcm_{method}"
+        result = data[basename]
 
-        scfres = psi4_run_pcm_hf(static_data.xyz["formaldehyde"], basis, charge=0,
-                                 multiplicity=1, conv_tol=1e-12,
-                                 conv_tol_grad=1e-11, max_iter=150,
-                                 pcm_options=psi4_pcm_options)
+        scfres = run_hf(static_data.xyz["formaldehyde"], basis, charge=0,
+                        multiplicity=1, conv_tol=1e-12, conv_tol_grad=1e-11,
+                        max_iter=150, pcm_options=pcm_options)
 
-        assert_allclose(scfres.energy_scf, psi4_result["energy_scf"], atol=1e-8)
+        assert_allclose(scfres.energy_scf, result["energy_scf"], atol=1e-8)
 
         matrix = adcc.AdcMatrix(method, scfres)
         solvent = AdcExtraTerm(matrix, {'ph_ph': block_ph_ph_0_pcm})
@@ -78,14 +94,14 @@ def template_pcm_psi4_linear_response_formaldehyde(self, basis, method):
                              environment=False)
         assert_allclose(
             state.excitation_energy_uncorrected,
-            psi4_result["lr_excitation_energy"],
+            result["lr_excitation_energy"],
             atol=1e-5
         )
 
         state_cis = adcc.ExcitedStates(state, property_method="adc0")
         assert_allclose(
             state_cis.oscillator_strength,
-            psi4_result["lr_osc_strength"], atol=1e-3
+            result["lr_osc_strength"], atol=1e-3
         )
 
         # invalid combination
@@ -102,14 +118,15 @@ def template_pcm_psi4_linear_response_formaldehyde(self, basis, method):
                              conv_tol=1e-7, environment="linear_response")
         assert_allclose(
             state.excitation_energy_uncorrected,
-            psi4_result["lr_excitation_energy"],
+            result["lr_excitation_energy"],
             atol=1e-5
         )
 
-        # remove cavity files from PSI4 PCM calculations
-        for cavityfile in os.listdir(os.getcwd()):
-            if cavityfile.startswith(("cavity.off_", "PEDRA.OUT_")):
-                os.remove(cavityfile)
+        if backend == "psi4":
+            # remove cavity files from PSI4 PCM calculations
+            for cavityfile in os.listdir(os.getcwd()):
+                if cavityfile.startswith(("cavity.off_", "PEDRA.OUT_")):
+                    os.remove(cavityfile)
 
 
 def psi4_run_pcm_hf(xyz, basis, charge=0, multiplicity=1, conv_tol=1e-12,
@@ -156,3 +173,34 @@ def psi4_run_pcm_hf(xyz, basis, charge=0, multiplicity=1, conv_tol=1e-12,
     _, wfn = psi4.energy('SCF', return_wfn=True, molecule=mol)
     psi4.core.clean()
     return adcc.backends.import_scf_results(wfn)
+
+
+def pyscf_run_pcm_hf(xyz, basis, charge=0, multiplicity=1, conv_tol=1e-11,
+                     conv_tol_grad=1e-10, max_iter=150, pcm_options=None):
+    from pyscf import scf, gto
+    from pyscf.solvent import ddCOSMO
+
+    mol = gto.M(
+        atom=xyz,
+        unit="Bohr",
+        basis=basis,
+        # spin in the pyscf world is 2S
+        spin=multiplicity - 1,
+        charge=charge,
+        # Disable commandline argument parsing in pyscf
+        parse_arg=False,
+        dump_input=False,
+        verbose=0,
+    )
+
+    mf = ddCOSMO(scf.RHF(mol))
+    # default eps
+    mf.with_solvent.eps = pcm_options.get("eps")
+    mf.conv_tol = conv_tol
+    mf.conv_tol_grad = conv_tol_grad
+    mf.max_cycle = max_iter
+
+    mf.kernel()
+    # replace eps with eps_opt for the ADC calculation
+    mf.with_solvent.eps = pcm_options.get("eps_opt")
+    return adcc.backends.import_scf_results(mf)
diff --git a/adcc/testdata/cache.py b/adcc/testdata/cache.py
index 87e5d136..95887f32 100644
--- a/adcc/testdata/cache.py
+++ b/adcc/testdata/cache.py
@@ -259,3 +259,4 @@ def read_yaml_data(fname):
 qchem_data = read_yaml_data("qchem_dump.yml")
 tmole_data = read_yaml_data("tmole_dump.yml")
 psi4_data = read_yaml_data("psi4_dump.yml")
+pyscf_data = read_yaml_data("pyscf_dump.yml")
diff --git a/adcc/testdata/generate_psi4_data.py b/adcc/testdata/generate_psi4_data.py
index 12d512f7..42de4394 100644
--- a/adcc/testdata/generate_psi4_data.py
+++ b/adcc/testdata/generate_psi4_data.py
@@ -56,9 +56,8 @@ def run_psi4_tdscf(xyz, basis, charge=0, multiplicity=1,
 
 
 def run_adcc_ptlr(wfn):
-    state = adcc.run_adc(wfn, method="adc1", n_singlets=5,
-                         conv_tol=1e-7, environment="ptlr")
-    return state
+    return adcc.run_adc(wfn, method="adc1", n_singlets=5,
+                        conv_tol=1e-7, environment="ptlr")
 
 
 def dump_results(molecule, basis, **kwargs):
diff --git a/adcc/testdata/generate_pyscf_data.py b/adcc/testdata/generate_pyscf_data.py
new file mode 100644
index 00000000..83cafb6c
--- /dev/null
+++ b/adcc/testdata/generate_pyscf_data.py
@@ -0,0 +1,94 @@
+import adcc
+from pyscf import gto, scf, tdscf
+from pyscf.solvent import ddCOSMO
+from adcc.testdata import static_data
+import yaml
+
+
+def run_pyscf_tdscf(xyz, basis, charge=0, multiplicity=1, conv_tol=1e-12,
+                    conv_tol_grad=1e-11, max_iter=150, pcm_options=None):
+    mol = gto.M(
+        atom=xyz,
+        basis=basis,
+        unit="Bohr",
+        charge=charge,
+        # spin in the pyscf world is 2S
+        spin=multiplicity - 1,
+        verbose=0,
+        # Disable commandline argument parsing in pyscf
+        parse_arg=False,
+        dump_input=False,
+    )
+
+    if pcm_options:
+        mf = ddCOSMO(scf.RHF(mol))
+        mf.with_solvent.eps = pcm_options.get("eps")
+    else:
+        mf = scf.RHF(mol)
+    mf.conv_tol = conv_tol
+    mf.conv_tol_grad = conv_tol_grad
+    mf.max_cycle = max_iter
+
+    mf.kernel()
+
+    if pcm_options:
+        mf.with_solvent.eps = pcm_options.get("eps_opt")
+        # for n_eq solvation only PTE implemented
+        mf.with_solvent.equilibrium_solvation = True
+        cis = ddCOSMO(tdscf.TDA(mf))
+    else:
+        cis = tdscf.TDA(mf)
+    cis.nstates = 5
+    cis.conv_tol = 1e-7
+    cis.kernel()
+    return mf, cis
+
+
+def run_adcc_ptlr(wfn):
+    return adcc.run_adc(wfn, method="adc1", n_singlets=5,
+                        conv_tol=1e-7, environment="ptlr")
+
+
+def dump_results(molecule, basis, **kwargs):
+    pcm_options = kwargs.get("pcm_options", None)
+    if pcm_options:
+        name = f"{molecule}_{basis}_pcm_adc1"
+    else:
+        name = f"{molecule}_{basis}_adc1"
+
+    geom = static_data.xyz[molecule]
+
+    ret = {"basis": basis,
+           "method": "adc1",
+           "molecule": molecule}
+
+    hf, cis = run_pyscf_tdscf(geom, basis, pcm_options=pcm_options)
+
+    ret["energy_scf"] = float(hf.e_tot)
+    ret["lr_excitation_energy"] = [float(round(s, 9)) for s in cis.e]
+    ret["lr_osc_strength"] = [float(round(s, 5)) for s in cis.oscillator_strength()]
+
+    if pcm_options:
+        state = run_adcc_ptlr(hf)
+        ret["ptlr_adcc_excitation_energy"] = [
+            round(float(e), 9) for e in state.excitation_energy
+        ]
+    return name, ret
+
+
+def main():
+    basis_set = ["sto-3g", "cc-pvdz"]
+    pcm_options = {"eps": 78.3553, "eps_opt": 1.78}
+    pyscf_results = {}
+    for basis in basis_set:
+        key, ret = dump_results("formaldehyde", basis, pcm_options=pcm_options)
+        pyscf_results[key] = ret
+        print(f"Dumped {key}.")
+
+    with open("pyscf_dump.yml", 'w') as yamlfile:
+        yaml.safe_dump(pyscf_results, yamlfile,
+                       sort_keys=False, default_flow_style=None)
+
+
+if __name__ == "__main__":
+    main()
diff --git a/adcc/testdata/pyscf_dump.yml b/adcc/testdata/pyscf_dump.yml
new file mode 100644
index 00000000..5788c519
--- /dev/null
+++ b/adcc/testdata/pyscf_dump.yml
@@ -0,0 +1,18 @@
+formaldehyde_sto-3g_pcm_adc1:
+  basis: sto-3g
+  method: adc1
+  molecule: formaldehyde
+  energy_scf: -112.35408506197172
+  lr_excitation_energy: [0.163950251, 0.360005329, 0.465135225, 0.509146753, 0.693115282]
+  lr_osc_strength: [0.0, 0.01098, 0.34455, 0.0, 0.35987]
+  ptlr_adcc_excitation_energy: [0.163957725, 0.360034834, 0.467270151, 0.509148683,
+    0.693769922]
+formaldehyde_cc-pvdz_pcm_adc1:
+  basis: cc-pvdz
+  method: adc1
+  molecule: formaldehyde
+  energy_scf: -113.88287611072602
+  lr_excitation_energy: [0.177403519, 0.379894734, 0.386588779, 0.390625944, 0.435452412]
+  lr_osc_strength: [0.0, 0.00012, 0.24985, 0.2374, 0.0]
+  ptlr_adcc_excitation_energy: [0.177443175, 0.379990015, 0.388208087, 0.391106802,
+    0.435484664]