Merge 1a3627c 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,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/test_backends_pcm.py

@@ -0,0 +1,158 @@
+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
+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
+
+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)))
+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")
+
+        # compare ptLR result to LR data
+        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-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)
+
+    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]
+
+        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)
+
+        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})
+
+        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,
+            psi4_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
+        )
+
+        # 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,
+            psi4_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)
+
+
+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)

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")
+psi4_data = read_yaml_data("psi4_dump.yml")

adcc/testdata/generate_psi4_data.py

@@ -0,0 +1,110 @@
+import psi4
+import adcc
+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.get("weight", 0.3)}
+            }}
+            Medium {{
+                Solvertype = {pcm_options.get("pcm_method", "IEFPCM")}
+                Solvent = {pcm_options.get("solvent", "Water")}
+                Nonequilibrium = {pcm_options.get("neq", True)}
+            }}
+        """)
+    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 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:
+        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["lr_excitation_energy"] = [round(float(r["EXCITATION ENERGY"]), 9)
+                                   for r in res]
+    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), 9) for e in state.excitation_energy
+        ]
+    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}.")
+
+    with open("psi4_dump.yml", 'w') as yamlfile:
+        yaml.safe_dump(psi4_results, yamlfile,
+                       sort_keys=False, default_flow_style=None)
+
+
+if __name__ == "__main__":
+    main()