Merge 0416b30 into fe3d42c

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

adcc/backends/pyscf.py

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

adcc/backends/test_backends_pcm.py

@@ -0,0 +1,204 @@
+import pytest
+import unittest
+import itertools
+import adcc
+import adcc.backends
+import os
+
+from adcc.misc import expand_test_templates
+from adcc.testdata import static_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"]
+
+
+@pytest.mark.skipif(len(backends) == 0, reason="No backend for PCM available.")
+@expand_test_templates(list(itertools.product(basissets, methods, backends)))
+class TestPCM(unittest.TestCase):
+    def template_pcm_ptlr_formaldehyde(self, basis, method, backend):
+        if method != "adc1":
+            pytest.skip("Data only available for adc1.")
+
+        c = config[backend]
+        basename = f"formaldehyde_{basis}_pcm_{method}"
+        result = c["data"][basename]
+
+        run_hf = c["run_hf"]
+        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=c["pcm_options"])
+
+        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,
+                        result["lr_excitation_energy"], atol=5e-3)
+
+        # Consistency check with values obtained with ADCc
+        assert_allclose(state.excitation_energy,
+                        result["ptlr_adcc_excitation_energy"], atol=1e-6)
+
+        if backend == "psi4":
+            # remove cavity files from PSI4 PCM calculations
+            remove_cavity_psi4()
+
+    def template_pcm_linear_response_formaldehyde(self, basis, method, backend):
+        if method != "adc1":
+            pytest.skip("Reference only exists for adc1.")
+
+        c = config[backend]
+        basename = f"formaldehyde_{basis}_pcm_{method}"
+        result = c["data"][basename]
+
+        run_hf = c["run_hf"]
+        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=c["pcm_options"])
+
+        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})
+
+        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,
+            result["lr_excitation_energy"],
+            atol=1e-5
+        )
+
+        state_cis = adcc.ExcitedStates(state, property_method="adc0")
+        assert_allclose(
+            state_cis.oscillator_strength,
+            result["lr_osc_strength"], atol=1e-3
+        )
+
+        # 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,
+            result["lr_excitation_energy"],
+            atol=1e-5
+        )
+
+        if backend == "psi4":
+            # remove cavity files from PSI4 PCM calculations
+            remove_cavity_psi4()
+
+
+def remove_cavity_psi4():
+    # removes 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,
+                    conv_tol_grad=1e-11, max_iter=150, pcm_options=None):
+    import psi4
+
+    # 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': 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.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({
+            '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)
+
+
+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)
+
+
+config = {
+    "psi4": {"data": psi4_data, "run_hf": psi4_run_pcm_hf,
+             "pcm_options": {"weight": 0.3, "pcm_method": "IEFPCM",
+                             "neq": True, "solvent": "Water"}},
+    "pyscf": {"data": pyscf_data, "run_hf": pyscf_run_pcm_hf,
+              "pcm_options": {"eps": 78.3553, "eps_opt": 1.78}}
+}