diff --git a/adcc/OperatorIntegrals.py b/adcc/OperatorIntegrals.py index 270529eb..0bfb0b9a 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_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() diff --git a/adcc/adc_pp/environment.py b/adcc/adc_pp/environment.py index 1ee02d42..aad254e0 100644 --- a/adcc/adc_pp/environment.py +++ b/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) diff --git a/adcc/backends/psi4.py b/adcc/backends/psi4.py index 94b3aa24..1e74b828 100644 --- a/adcc/backends/psi4.py +++ b/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", 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..c1d00c41 --- /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, + 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) 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 d41261df..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): + env=(None, None)): """ Run the SCF for a backend and a particular test case (if not done) and return the result. @@ -189,13 +189,21 @@ 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 @@ -203,7 +211,11 @@ def payload(): 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: 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