Merge a36c658 into e4ca2ed

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/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,15 @@ 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):
@@ -128,6 +137,8 @@ def environment(self):
         ret = None
         if hasattr(self.wfn, "pe_state"):
             ret = "pe"
+        if self.wfn.PCM_enabled():
+            ret = "pcm"
         return ret
 
     def get_backend(self):
@@ -238,13 +249,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):
+           conv_tol_grad=1e-8, max_iter=150, pe_options=None, pcm=None):
     basissets = {
         "sto3g": "sto-3g",
         "def2tzvp": "def2-tzvp",
         "ccpvdz": "cc-pvdz",
     }
 
+    psi4.core.clean_options()
     mol = psi4.geometry(f"""
         {charge} {multiplicity}
         {xyz}
@@ -267,6 +279,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",

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,

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,
+                                   env=("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)

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)

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):
+                      env=(None, None)):
     """
     Run the SCF for a backend and a particular test case (if not done)
     and return the result.
@@ -189,21 +189,33 @@ 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)
+                                     pe_options=pe_options, pcm=pcm)
         return adcc.backends.import_scf_results(hfres)
 
     # For reasons not clear to me (mfh), caching does not work
     # with pyscf
     if backend == "pyscf":
         return payload()
 
-    key = (backend, molecule, basis, str(multiplicity))
+    # 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))
     try:
         return __cache_cached_backend_hf[key]
     except NameError:

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")

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]