From 69ddf100dec424a0cb78afdd457154ac6603be8c Mon Sep 17 00:00:00 2001 From: Maximilian Scheurer Date: Wed, 9 Sep 2020 21:51:17 +0200 Subject: [PATCH 01/21] PE ADC postSCF coupling --- adcc/OperatorIntegrals.py | 20 ++++++++++++++++++++ adcc/PeAdcMatrix.py | 22 ++++++++++++++++++++++ adcc/backends/psi4.py | 11 +++++++++++ adcc/backends/pyscf.py | 18 +++++++++++++++++- 4 files changed, 70 insertions(+), 1 deletion(-) create mode 100644 adcc/PeAdcMatrix.py diff --git a/adcc/OperatorIntegrals.py b/adcc/OperatorIntegrals.py index a3ae3501..1f0be057 100644 --- a/adcc/OperatorIntegrals.py +++ b/adcc/OperatorIntegrals.py @@ -137,6 +137,26 @@ def nabla(self): """ return self.import_dipole_like_operator("nabla", is_symmetric=False) + @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/PeAdcMatrix.py b/adcc/PeAdcMatrix.py new file mode 100644 index 00000000..5e8cc676 --- /dev/null +++ b/adcc/PeAdcMatrix.py @@ -0,0 +1,22 @@ +from adcc import AdcMatrix +from .OneParticleOperator import OneParticleOperator + + +class PeAdcMatrix(AdcMatrix): + def matvec(self, in_ampl): + out_ampl = super().matvec(in_ampl) + operators = self.reference_state.operators + tdm = OneParticleOperator(self.ground_state, is_symmetric=False) + tdm.vo = in_ampl.ph.transpose() + vpe = operators.density_dependent_operators["pe_induction_elec"](tdm) + out_ampl.ph += vpe.ov + return out_ampl + + +# TODO: replace with existing shifted matrix +class PeShiftedMat(PeAdcMatrix): + omega = 0 + + def matvec(self, in_ampl): + out_ampl = super().matvec(in_ampl) + return out_ampl - in_ampl * self.omega diff --git a/adcc/backends/psi4.py b/adcc/backends/psi4.py index c9714b93..66c861ca 100644 --- a/adcc/backends/psi4.py +++ b/adcc/backends/psi4.py @@ -53,6 +53,17 @@ def magnetic_dipole(self): def nabla(self): return [-1.0 * np.asarray(comp) for comp in self.mints.ao_nabla()] + @property + def density_dependent_operators(self): + ret = {} + if hasattr(self.wfn, "pe_state"): + ret["pe_induction_elec"] = lambda dm: \ + self.wfn.pe_state.get_pe_contribution( + psi4.core.Matrix.from_array(dm.to_ndarray()), + elec_only=True + )[1] + return ret + class Psi4EriBuilder(EriBuilder): def __init__(self, wfn, n_orbs, n_orbs_alpha, n_alpha, n_beta, restricted): diff --git a/adcc/backends/pyscf.py b/adcc/backends/pyscf.py index dcafc719..5afc462d 100644 --- a/adcc/backends/pyscf.py +++ b/adcc/backends/pyscf.py @@ -21,6 +21,7 @@ ## ## --------------------------------------------------------------------- import numpy as np +import warnings from libadcc import HartreeFockProvider from adcc.misc import cached_property @@ -38,7 +39,12 @@ def __init__(self, scfres): @cached_property def electric_dipole(self): - return list(self.scfres.mol.intor_symmetric('int1e_r', comp=3)) + if hasattr(self.scfres, "with_solvent"): + if self.scfres.with_solvent.eef: + warnings.warn("Using modified dipole operator due to EEF.") + return self.scfres.with_solvent.effective_dipole_operator() + else: + return list(self.scfres.mol.intor_symmetric('int1e_r', comp=3)) @cached_property def magnetic_dipole(self): @@ -55,6 +61,16 @@ def nabla(self): -1.0 * self.scfres.mol.intor('int1e_ipovlp', comp=3, hermi=2) ) + @property + def density_dependent_operators(self): + ret = {} + if hasattr(self.scfres, "with_solvent"): + if isinstance(self.scfres.with_solvent, solvent.pol_embed.PolEmbed): + ret["pe_induction_elec"] = lambda dm: \ + self.scfres.with_solvent._exec_cppe(dm.to_ndarray(), + elec_only=True)[1] + return ret + # TODO: refactor ERI builder to be more general # IntegralBuilder would be good From 0567ad69ac889c8124c2b20a1d80ec138816df0a Mon Sep 17 00:00:00 2001 From: Maximilian Scheurer Date: Thu, 7 Jan 2021 17:09:02 +0100 Subject: [PATCH 02/21] adapt to new amplitude vector --- adcc/PeAdcMatrix.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/adcc/PeAdcMatrix.py b/adcc/PeAdcMatrix.py index 5e8cc676..0e8ab27f 100644 --- a/adcc/PeAdcMatrix.py +++ b/adcc/PeAdcMatrix.py @@ -9,7 +9,8 @@ def matvec(self, in_ampl): tdm = OneParticleOperator(self.ground_state, is_symmetric=False) tdm.vo = in_ampl.ph.transpose() vpe = operators.density_dependent_operators["pe_induction_elec"](tdm) - out_ampl.ph += vpe.ov + # TODO: out_ampl.ph += does not seem to work?! + out_ampl['ph'] += vpe.ov return out_ampl From 9a8b6b2d08074e78b2854f85c680e36bc8587916 Mon Sep 17 00:00:00 2001 From: Maximilian Scheurer Date: Sun, 17 Jan 2021 14:18:03 +0100 Subject: [PATCH 03/21] refactor matrix apply --- adcc/PeAdcMatrix.py | 29 +++++++++++------------------ 1 file changed, 11 insertions(+), 18 deletions(-) diff --git a/adcc/PeAdcMatrix.py b/adcc/PeAdcMatrix.py index 0e8ab27f..3bd3d872 100644 --- a/adcc/PeAdcMatrix.py +++ b/adcc/PeAdcMatrix.py @@ -3,21 +3,14 @@ class PeAdcMatrix(AdcMatrix): - def matvec(self, in_ampl): - out_ampl = super().matvec(in_ampl) - operators = self.reference_state.operators - tdm = OneParticleOperator(self.ground_state, is_symmetric=False) - tdm.vo = in_ampl.ph.transpose() - vpe = operators.density_dependent_operators["pe_induction_elec"](tdm) - # TODO: out_ampl.ph += does not seem to work?! - out_ampl['ph'] += vpe.ov - return out_ampl - - -# TODO: replace with existing shifted matrix -class PeShiftedMat(PeAdcMatrix): - omega = 0 - - def matvec(self, in_ampl): - out_ampl = super().matvec(in_ampl) - return out_ampl - in_ampl * self.omega + def block_apply(self, block, in_vec): + ret = super().block_apply(block, in_vec) + if block == "ph_ph": + # CIS-like contribution only adds to ph_ph block + with self.timer.record("apply/pe_coupling"): + op = self.reference_state.operators + tdm = OneParticleOperator(self.ground_state, is_symmetric=False) + tdm.vo = in_vec.ph.transpose() + vpe = op.density_dependent_operators["pe_induction_elec"](tdm) + ret.ph += vpe.ov + return ret From fde961a1624abcb2d78b46fb77e3574b4f596a0f Mon Sep 17 00:00:00 2001 From: Maximilian Scheurer Date: Sun, 21 Feb 2021 12:40:32 +0100 Subject: [PATCH 04/21] nuke extra matrix, play with AdcBlock --- adcc/PeAdcMatrix.py | 16 ---------- adcc/backends/test_backends_polembed.py | 42 ++++++++++++++++++++++++- 2 files changed, 41 insertions(+), 17 deletions(-) delete mode 100644 adcc/PeAdcMatrix.py diff --git a/adcc/PeAdcMatrix.py b/adcc/PeAdcMatrix.py deleted file mode 100644 index 3bd3d872..00000000 --- a/adcc/PeAdcMatrix.py +++ /dev/null @@ -1,16 +0,0 @@ -from adcc import AdcMatrix -from .OneParticleOperator import OneParticleOperator - - -class PeAdcMatrix(AdcMatrix): - def block_apply(self, block, in_vec): - ret = super().block_apply(block, in_vec) - if block == "ph_ph": - # CIS-like contribution only adds to ph_ph block - with self.timer.record("apply/pe_coupling"): - op = self.reference_state.operators - tdm = OneParticleOperator(self.ground_state, is_symmetric=False) - tdm.vo = in_vec.ph.transpose() - vpe = op.density_dependent_operators["pe_induction_elec"](tdm) - ret.ph += vpe.ov - return ret diff --git a/adcc/backends/test_backends_polembed.py b/adcc/backends/test_backends_polembed.py index 3123f8bd..2adc2dd0 100644 --- a/adcc/backends/test_backends_polembed.py +++ b/adcc/backends/test_backends_polembed.py @@ -22,9 +22,13 @@ ## --------------------------------------------------------------------- import unittest import itertools +import numpy as np import adcc import adcc.backends +from adcc.adc_pp.matrix import AdcBlock +from adcc import OneParticleOperator, AmplitudeVector + from numpy.testing import assert_allclose import pytest @@ -52,7 +56,7 @@ @pytest.mark.skipif(len(backends) == 0, reason="No backend found.") @expand_test_templates(list(itertools.product(basissets, methods, backends))) class TestPolarizableEmbedding(unittest.TestCase): - def template_pe_formaldehyde(self, basis, method, backend): + def template_pe_formaldehyde_perturbative(self, basis, method, backend): basename = f"formaldehyde_{basis}_pe_{method}" qc_result = qchem_data[basename] pe_options = {"potfile": pe_potentials["fa_6w"]} @@ -82,3 +86,39 @@ def template_pe_formaldehyde(self, basis, method, backend): state.pe_ptlr_correction, atol=1e-5 ) + + def test_pe_formaldehyde_coupling_sto3g_pe_adc2(self): + basis = "sto-3g" + method = "adc2" + backend = "pyscf" + pe_options = {"potfile": pe_potentials["fa_6w"]} + scfres = cached_backend_hf(backend, "formaldehyde", basis, + pe_options=pe_options) + assert_allclose(scfres.energy_scf, -112.36904349555, atol=1e-8) + + # construct a normal ADC matrix + matrix = adcc.AdcMatrix(method, scfres) + + # additional matrix apply for PE + def block_ph_ph_0_pe(hf, mp, intermediates): + op = hf.operators + + def apply(ampl): + tdm = OneParticleOperator(mp, is_symmetric=False) + tdm.vo = ampl.ph.transpose() + vpe = op.density_dependent_operators["pe_induction_elec"](tdm) + return AmplitudeVector(ph=vpe.ov) + return AdcBlock(apply, 0) + + # register the additional function with the matrix + # this should be implemented in the AdcMatrix ctor + # (to also get all the diagonal setup right) + matrix.blocks_ph['ph_ph_pe'] = block_ph_ph_0_pe( + matrix.reference_state, matrix.ground_state, matrix.intermediates + ) + assert_allclose(matrix.ground_state.energy(2), -112.4804685653, atol=1e-8) + exci_tm = np.array([ + 0.1733632144, 0.3815240343, 0.5097618218, 0.5189443358, 0.5670575778 + ]) + state = adcc.run_adc(matrix, n_singlets=5, conv_tol=1e-7) + assert_allclose(state.excitation_energy_uncorrected, exci_tm, atol=1e-6) From de5bbb93d9a286aca75d8a2a40e409859c43d26e Mon Sep 17 00:00:00 2001 From: Maximilian Scheurer Date: Tue, 23 Feb 2021 10:19:31 +0100 Subject: [PATCH 05/21] add extra test case, yml file --- adcc/backends/pyscf.py | 7 +----- adcc/backends/test_backends_polembed.py | 31 ++++++++++++++----------- adcc/testdata/cache.py | 7 +++--- adcc/testdata/tmole_dump.yml | 14 +++++++++++ 4 files changed, 37 insertions(+), 22 deletions(-) create mode 100644 adcc/testdata/tmole_dump.yml diff --git a/adcc/backends/pyscf.py b/adcc/backends/pyscf.py index 5afc462d..5c9b3b61 100644 --- a/adcc/backends/pyscf.py +++ b/adcc/backends/pyscf.py @@ -39,12 +39,7 @@ def __init__(self, scfres): @cached_property def electric_dipole(self): - if hasattr(self.scfres, "with_solvent"): - if self.scfres.with_solvent.eef: - warnings.warn("Using modified dipole operator due to EEF.") - return self.scfres.with_solvent.effective_dipole_operator() - else: - return list(self.scfres.mol.intor_symmetric('int1e_r', comp=3)) + return list(self.scfres.mol.intor_symmetric('int1e_r', comp=3)) @cached_property def magnetic_dipole(self): diff --git a/adcc/backends/test_backends_polembed.py b/adcc/backends/test_backends_polembed.py index 2adc2dd0..1323b730 100644 --- a/adcc/backends/test_backends_polembed.py +++ b/adcc/backends/test_backends_polembed.py @@ -22,7 +22,6 @@ ## --------------------------------------------------------------------- import unittest import itertools -import numpy as np import adcc import adcc.backends @@ -35,7 +34,7 @@ from ..misc import expand_test_templates from .testing import cached_backend_hf -from ..testdata.cache import qchem_data +from ..testdata.cache import qchem_data, tmole_data from ..testdata.static_data import pe_potentials try: @@ -56,7 +55,7 @@ @pytest.mark.skipif(len(backends) == 0, reason="No backend found.") @expand_test_templates(list(itertools.product(basissets, methods, backends))) class TestPolarizableEmbedding(unittest.TestCase): - def template_pe_formaldehyde_perturbative(self, basis, method, backend): + def template_pe_perturbative_formaldehyde(self, basis, method, backend): basename = f"formaldehyde_{basis}_pe_{method}" qc_result = qchem_data[basename] pe_options = {"potfile": pe_potentials["fa_6w"]} @@ -87,14 +86,15 @@ def template_pe_formaldehyde_perturbative(self, basis, method, backend): atol=1e-5 ) - def test_pe_formaldehyde_coupling_sto3g_pe_adc2(self): - basis = "sto-3g" - method = "adc2" - backend = "pyscf" + def template_pe_coupling_formaldehyde(self, basis, method, backend): + if method != "adc2": + pytest.skip("") + basename = f"formaldehyde_{basis}_pe_{method}" + tm_result = tmole_data[basename] pe_options = {"potfile": pe_potentials["fa_6w"]} scfres = cached_backend_hf(backend, "formaldehyde", basis, pe_options=pe_options) - assert_allclose(scfres.energy_scf, -112.36904349555, atol=1e-8) + assert_allclose(scfres.energy_scf, tm_result["energy_scf"], atol=1e-8) # construct a normal ADC matrix matrix = adcc.AdcMatrix(method, scfres) @@ -116,9 +116,14 @@ def apply(ampl): matrix.blocks_ph['ph_ph_pe'] = block_ph_ph_0_pe( matrix.reference_state, matrix.ground_state, matrix.intermediates ) - assert_allclose(matrix.ground_state.energy(2), -112.4804685653, atol=1e-8) - exci_tm = np.array([ - 0.1733632144, 0.3815240343, 0.5097618218, 0.5189443358, 0.5670575778 - ]) + assert_allclose( + matrix.ground_state.energy(2), + tm_result["energy_mp2"], + atol=1e-8 + ) state = adcc.run_adc(matrix, n_singlets=5, conv_tol=1e-7) - assert_allclose(state.excitation_energy_uncorrected, exci_tm, atol=1e-6) + assert_allclose( + state.excitation_energy_uncorrected, + tm_result["excitation_energy"], + atol=1e-6 + ) diff --git a/adcc/testdata/cache.py b/adcc/testdata/cache.py index 447449fd..4dbfab89 100644 --- a/adcc/testdata/cache.py +++ b/adcc/testdata/cache.py @@ -248,12 +248,13 @@ def lists_to_ndarray(dictionary): return data -def read_qchem_data(): +def read_yaml_data(fname): thisdir = os.path.dirname(__file__) - yaml_file = os.path.join(thisdir, "qchem_dump.yml") + yaml_file = os.path.join(thisdir, fname) with open(yaml_file, "r") as f: data_raw = yaml.safe_load(f) return lists_to_ndarray(data_raw) -qchem_data = read_qchem_data() +qchem_data = read_yaml_data("qchem_dump.yml") +tmole_data = read_yaml_data("tmole_dump.yml") diff --git a/adcc/testdata/tmole_dump.yml b/adcc/testdata/tmole_dump.yml new file mode 100644 index 00000000..88e8660e --- /dev/null +++ b/adcc/testdata/tmole_dump.yml @@ -0,0 +1,14 @@ +formaldehyde_sto3g_pe_adc2: + basis: sto-3g + method: adc2 + molecule: formaldehyde + energy_scf: -112.36904349555 + energy_mp2: -112.4804685653 + excitation_energy: [0.1733632144, 0.3815240343, 0.5097618218, 0.5189443358, 0.5670575778] +formaldehyde_ccpvdz_pe_adc2: + basis: ccpvdz + method: adc2 + molecule: formaldehyde + energy_scf: -113.89976872224 + energy_mp2: -114.2156673890 + excitation_energy: [0.1596037963, 0.3123344811, 0.3621306427, 0.3795572252, 0.4089162982] \ No newline at end of file From 6ae427e6705ed598d2722e6a461b581886b3512f Mon Sep 17 00:00:00 2001 From: Maximilian Scheurer Date: Wed, 10 Mar 2021 10:03:51 +0100 Subject: [PATCH 06/21] solvent block in matrix --- adcc/AdcMatrix.py | 14 +++++++++++++- adcc/backends/pyscf.py | 1 - adcc/backends/test_backends_polembed.py | 10 +--------- 3 files changed, 14 insertions(+), 11 deletions(-) diff --git a/adcc/AdcMatrix.py b/adcc/AdcMatrix.py index d40d2058..7ba3e5e0 100644 --- a/adcc/AdcMatrix.py +++ b/adcc/AdcMatrix.py @@ -50,7 +50,8 @@ class AdcMatrix(AdcMatrixlike): "adc3": dict(ph_ph=3, ph_pphh=2, pphh_ph=2, pphh_pphh=1), # noqa: E501 } - def __init__(self, method, hf_or_mp, block_orders=None, intermediates=None): + def __init__(self, method, hf_or_mp, block_orders=None, intermediates=None, + solvent_block=None): """ Initialise an ADC matrix. @@ -65,6 +66,9 @@ def __init__(self, method, hf_or_mp, block_orders=None, intermediates=None): If not set, defaults according to the selected ADC method are chosen. intermediates : adcc.Intermediates or NoneType Allows to pass intermediates to re-use to this class. + solvent_block: callable + Function to construct an AdcBlock responsible for + coupling to a solvent model """ if isinstance(hf_or_mp, (libadcc.ReferenceState, libadcc.HartreeFockSolution_i)): @@ -125,6 +129,14 @@ def __init__(self, method, hf_or_mp, block_orders=None, intermediates=None): self.__diagonal.evaluate() self.__init_space_data(self.__diagonal) + if solvent_block is not None: + if not callable(solvent_block): + raise TypeError("solvent_block needs to be" + " a callable.") + self.blocks_ph["solvent"] = solvent_block( + self.reference_state, self.ground_state, self.intermediates + ) + def __init_space_data(self, diagonal): """Update the cached data regarding the spaces of the ADC matrix""" self.axis_spaces = {} diff --git a/adcc/backends/pyscf.py b/adcc/backends/pyscf.py index 5c9b3b61..ab62b486 100644 --- a/adcc/backends/pyscf.py +++ b/adcc/backends/pyscf.py @@ -21,7 +21,6 @@ ## ## --------------------------------------------------------------------- import numpy as np -import warnings from libadcc import HartreeFockProvider from adcc.misc import cached_property diff --git a/adcc/backends/test_backends_polembed.py b/adcc/backends/test_backends_polembed.py index 1323b730..4e9008a2 100644 --- a/adcc/backends/test_backends_polembed.py +++ b/adcc/backends/test_backends_polembed.py @@ -96,9 +96,6 @@ def template_pe_coupling_formaldehyde(self, basis, method, backend): pe_options=pe_options) assert_allclose(scfres.energy_scf, tm_result["energy_scf"], atol=1e-8) - # construct a normal ADC matrix - matrix = adcc.AdcMatrix(method, scfres) - # additional matrix apply for PE def block_ph_ph_0_pe(hf, mp, intermediates): op = hf.operators @@ -110,12 +107,7 @@ def apply(ampl): return AmplitudeVector(ph=vpe.ov) return AdcBlock(apply, 0) - # register the additional function with the matrix - # this should be implemented in the AdcMatrix ctor - # (to also get all the diagonal setup right) - matrix.blocks_ph['ph_ph_pe'] = block_ph_ph_0_pe( - matrix.reference_state, matrix.ground_state, matrix.intermediates - ) + matrix = adcc.AdcMatrix(method, scfres, solvent_block=block_ph_ph_0_pe) assert_allclose( matrix.ground_state.energy(2), tm_result["energy_mp2"], From 6361637f292a034ee6ee24a09d4bda0b639c1e46 Mon Sep 17 00:00:00 2001 From: Maximilian Scheurer Date: Tue, 23 Mar 2021 12:25:50 +0100 Subject: [PATCH 07/21] fix psi4 pe bug --- 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 66c861ca..d4a6e031 100644 --- a/adcc/backends/psi4.py +++ b/adcc/backends/psi4.py @@ -250,7 +250,7 @@ def run_hf(xyz, basis, charge=0, multiplicity=1, conv_tol=1e-11, 'reference': "RHF", }) if pe_options: - psi4.set_options({"pe", "true"}) + psi4.set_options({"pe": "true"}) psi4.set_module_options("pe", {"potfile": pe_options["potfile"]}) if multiplicity != 1: From 49dafeade3e39e4add9f4762a50d3502554bc361 Mon Sep 17 00:00:00 2001 From: Maximilian Scheurer Date: Wed, 9 Sep 2020 21:51:17 +0200 Subject: [PATCH 08/21] PE ADC postSCF coupling adapt to new amplitude vector refactor matrix apply nuke extra matrix, play with AdcBlock add extra test case, yml file solvent block in matrix fix psi4 pe bug use composable AdcBlock verbose exceptions, extra file with solvent matrix equations add style file --- .clang-format | 89 +++++++++++++++++++++++++ .github/workflows/ci.yaml | 2 +- adcc/AdcMatrix.py | 19 +++++- adcc/OperatorIntegrals.py | 20 ++++++ adcc/adc_pp/matrix.py | 61 ++++++++++++++--- adcc/adc_pp/solvent.py | 42 ++++++++++++ adcc/backends/psi4.py | 13 +++- adcc/backends/pyscf.py | 10 +++ adcc/backends/test_backends_polembed.py | 33 ++++++++- adcc/testdata/cache.py | 7 +- adcc/testdata/tmole_dump.yml | 14 ++++ 11 files changed, 291 insertions(+), 19 deletions(-) create mode 100644 .clang-format create mode 100644 adcc/adc_pp/solvent.py create mode 100644 adcc/testdata/tmole_dump.yml diff --git a/.clang-format b/.clang-format new file mode 100644 index 00000000..d6d616bc --- /dev/null +++ b/.clang-format @@ -0,0 +1,89 @@ +## --------------------------------------------------------------------- +## +## Dummy licence header +## +## --------------------------------------------------------------------- + +## DO NOT EDIT +## This file is automatically generated from a file in the repository "krims". +## Edit the original and call the script "update_from_sister_repos.sh" instead. + +--- +Language: Cpp +AccessModifierOffset: -1 +AlignAfterOpenBracket: Align +AlignConsecutiveAssignments: true +AlignConsecutiveDeclarations: false +AlignEscapedNewlinesLeft: true +AlignOperands: true +AlignTrailingComments: true +AllowAllParametersOfDeclarationOnNextLine: true +AllowShortBlocksOnASingleLine: false +AllowShortCaseLabelsOnASingleLine: false +AllowShortFunctionsOnASingleLine: All +AllowShortIfStatementsOnASingleLine: true +AllowShortLoopsOnASingleLine: true +AlwaysBreakAfterReturnType: None +AlwaysBreakBeforeMultilineStrings: true +AlwaysBreakTemplateDeclarations: true +BinPackArguments: true +BinPackParameters: true +BraceWrapping: + AfterClass: false + AfterControlStatement: false + AfterEnum: false + AfterFunction: false + AfterNamespace: false + AfterObjCDeclaration: false + AfterStruct: false + AfterUnion: false + BeforeCatch: false + BeforeElse: false + IndentBraces: false +BreakBeforeBinaryOperators: None +BreakBeforeBraces: Attach +BreakBeforeTernaryOperators: true +BreakConstructorInitializersBeforeComma: false +ColumnLimit: 90 +CommentPragmas: '^ IWYU pragma:' +ConstructorInitializerAllOnOneLineOrOnePerLine: true +ConstructorInitializerIndentWidth: 6 +ContinuationIndentWidth: 6 +Cpp11BracedListStyle: true +DerivePointerAlignment: false +DisableFormat: false +ExperimentalAutoDetectBinPacking: false +ForEachMacros: [ foreach, Q_FOREACH, BOOST_FOREACH ] +IndentCaseLabels: true +IndentWidth: 2 +IndentWrappedFunctionNames: false +KeepEmptyLinesAtTheStartOfBlocks: true +MaxEmptyLinesToKeep: 1 +NamespaceIndentation: None +ObjCBlockIndentWidth: 2 +ObjCSpaceAfterProperty: false +ObjCSpaceBeforeProtocolList: false +PenaltyBreakBeforeFirstCallParameter: 1 +PenaltyBreakComment: 300 +PenaltyBreakFirstLessLess: 120 +PenaltyBreakString: 1000 +PenaltyExcessCharacter: 1000000 +PenaltyReturnTypeOnItsOwnLine: 200 +PointerAlignment: Left +ReflowComments: true +SortIncludes: true +SpaceAfterCStyleCast: false +SpaceBeforeAssignmentOperators: true +SpaceBeforeParens: ControlStatements +SpaceInEmptyParentheses: false +SpacesBeforeTrailingComments: 2 +SpacesInAngles: false +SpacesInContainerLiterals: true +SpacesInCStyleCastParentheses: false +SpacesInParentheses: false +SpacesInSquareBrackets: false +Standard: Cpp11 +TabWidth: 8 +UseTab: Never +... + diff --git a/.github/workflows/ci.yaml b/.github/workflows/ci.yaml index 029a3852..520876c6 100644 --- a/.github/workflows/ci.yaml +++ b/.github/workflows/ci.yaml @@ -119,5 +119,5 @@ jobs: - name: Checking C++ files with clang-format working-directory: libadcc run: | - find . -type f \( -name '*.cc' -o -name '*.hh' \) -print0 | xargs -r0 clang-format -i + find . -type f \( -name '*.cc' -o -name '*.hh' \) -print0 | xargs -r0 clang-format --style=file -i git diff --exit-code diff --git a/adcc/AdcMatrix.py b/adcc/AdcMatrix.py index d40d2058..593dc589 100644 --- a/adcc/AdcMatrix.py +++ b/adcc/AdcMatrix.py @@ -50,7 +50,8 @@ class AdcMatrix(AdcMatrixlike): "adc3": dict(ph_ph=3, ph_pphh=2, pphh_ph=2, pphh_pphh=1), # noqa: E501 } - def __init__(self, method, hf_or_mp, block_orders=None, intermediates=None): + def __init__(self, method, hf_or_mp, block_orders=None, intermediates=None, + additional_blocks=None): """ Initialise an ADC matrix. @@ -120,6 +121,22 @@ def __init__(self, method, hf_or_mp, block_orders=None, intermediates=None): variant=variant) for block, order in block_orders.items() if order is not None } + if additional_blocks is not None: + if not isinstance(additional_blocks, dict): + raise TypeError("additional_blocks needs to " + " be a dict.") + for space in additional_blocks: + if space not in self.blocks_ph: + raise ValueError("Can only add blocks" + " to existing matrix blocks.") + block_fun = additional_blocks[space] + if not callable(block_fun): + raise TypeError("Items in additional_blocks " + "must be callable.") + block = block_fun( + self.reference_state, self.ground_state, self.intermediates + ) + self.blocks_ph[space].add_block(block) self.__diagonal = sum(bl.diagonal for bl in self.blocks_ph.values() if bl.diagonal) self.__diagonal.evaluate() diff --git a/adcc/OperatorIntegrals.py b/adcc/OperatorIntegrals.py index a3ae3501..1f0be057 100644 --- a/adcc/OperatorIntegrals.py +++ b/adcc/OperatorIntegrals.py @@ -137,6 +137,26 @@ def nabla(self): """ return self.import_dipole_like_operator("nabla", is_symmetric=False) + @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/matrix.py b/adcc/adc_pp/matrix.py index f4892a37..cbafba54 100644 --- a/adcc/adc_pp/matrix.py +++ b/adcc/adc_pp/matrix.py @@ -21,7 +21,6 @@ ## ## --------------------------------------------------------------------- from math import sqrt -from collections import namedtuple from adcc import block as b from adcc.functions import direct_sum, einsum, zeros_like @@ -39,15 +38,53 @@ # really so much our focus. -# -# Dispatch routine -# -""" -`apply` is a function mapping an AmplitudeVector to the contribution of this -block to the result of applying the ADC matrix. `diagonal` is an `AmplitudeVector` -containing the expression to the diagonal of the ADC matrix from this block. -""" -AdcBlock = namedtuple("AdcBlock", ["apply", "diagonal"]) +class AdcBlock: + def __init__(self, apply, diagonal): + """AdcBlock contains the matrix apply and diagonal + routines for a specific ADC matrix block and is used to compose + the dispatch routines for :py:class:`AdcMatrix` + + Parameters + ---------- + apply : callable + function mapping an AmplitudeVector to the contribution of + this block to the result of applying the ADC matrix + diagonal : AmplitudeVector, int, float + expression to the diagonal of the ADC matrix from this block + """ + if not callable(apply): + raise TypeError("apply needs to be callable.") + if not isinstance(diagonal, (AmplitudeVector, int, float)): + raise TypeError("diagonal needs to be an AmplitudeVector or a number.") + self._applies = [apply] + self._diagonals = [diagonal] + + def add_block(self, other): + """Adds another :py:class:`AdcBlock` to this AdcBlock. + + Parameters + ---------- + other : AdcBlock + block to be added + """ + if not isinstance(other, AdcBlock): + raise TypeError("other must be of type AdcBlock.") + self._applies.append(other.apply) + self._diagonals.append(other.diagonal) + + def apply(self, ampl): + """Applies the block to an input AmplitudeVector + + Parameters + ---------- + ampl : AmplitudeVector + Input AmplitudeVector + """ + return sum(app(ampl) for app in self._applies) + + @property + def diagonal(self): + return sum(self._diagonals) def block(ground_state, spaces, order, variant=None, intermediates=None): @@ -58,7 +95,9 @@ def block(ground_state, spaces, order, variant=None, intermediates=None): It is assumed largely, that CVS is equivalent to mp.has_core_occupied_space, while one would probably want in the long run that one can have an "o2" space, - but not do CVS + but not do CVS. + + Returns an instance of :py:class:`AdcBlock`. """ if isinstance(variant, str): variant = [variant] diff --git a/adcc/adc_pp/solvent.py b/adcc/adc_pp/solvent.py new file mode 100644 index 00000000..c7624c6d --- /dev/null +++ b/adcc/adc_pp/solvent.py @@ -0,0 +1,42 @@ +#!/usr/bin/env python3 +## vi: tabstop=4 shiftwidth=4 softtabstop=4 expandtab +## --------------------------------------------------------------------- +## +## Copyright (C) 2021 by the adcc authors +## +## This file is part of adcc. +## +## adcc is free software: you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published +## by the Free Software Foundation, either version 3 of the License, or +## (at your option) any later version. +## +## adcc is distributed in the hope that it will be useful, +## but WITHOUT ANY WARRANTY; without even the implied warranty of +## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +## GNU General Public License for more details. +## +## You should have received a copy of the GNU General Public License +## along with adcc. If not, see . +## +## --------------------------------------------------------------------- +from adcc import OneParticleOperator, AmplitudeVector +from .matrix import AdcBlock + + +def block_ph_ph_0_pe(hf, mp, intermediates): + """ + Constructs an :py:class:`AdcBlock` that describes the + coupling to the polarizable environment from PE via a CIS-like + coupling 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() + vpe = op.density_dependent_operators["pe_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 c9714b93..d4a6e031 100644 --- a/adcc/backends/psi4.py +++ b/adcc/backends/psi4.py @@ -53,6 +53,17 @@ def magnetic_dipole(self): def nabla(self): return [-1.0 * np.asarray(comp) for comp in self.mints.ao_nabla()] + @property + def density_dependent_operators(self): + ret = {} + if hasattr(self.wfn, "pe_state"): + ret["pe_induction_elec"] = lambda dm: \ + self.wfn.pe_state.get_pe_contribution( + psi4.core.Matrix.from_array(dm.to_ndarray()), + elec_only=True + )[1] + return ret + class Psi4EriBuilder(EriBuilder): def __init__(self, wfn, n_orbs, n_orbs_alpha, n_alpha, n_beta, restricted): @@ -239,7 +250,7 @@ def run_hf(xyz, basis, charge=0, multiplicity=1, conv_tol=1e-11, 'reference': "RHF", }) if pe_options: - psi4.set_options({"pe", "true"}) + psi4.set_options({"pe": "true"}) psi4.set_module_options("pe", {"potfile": pe_options["potfile"]}) if multiplicity != 1: diff --git a/adcc/backends/pyscf.py b/adcc/backends/pyscf.py index dcafc719..ab62b486 100644 --- a/adcc/backends/pyscf.py +++ b/adcc/backends/pyscf.py @@ -55,6 +55,16 @@ def nabla(self): -1.0 * self.scfres.mol.intor('int1e_ipovlp', comp=3, hermi=2) ) + @property + def density_dependent_operators(self): + ret = {} + if hasattr(self.scfres, "with_solvent"): + if isinstance(self.scfres.with_solvent, solvent.pol_embed.PolEmbed): + ret["pe_induction_elec"] = lambda dm: \ + self.scfres.with_solvent._exec_cppe(dm.to_ndarray(), + elec_only=True)[1] + return ret + # TODO: refactor ERI builder to be more general # IntegralBuilder would be good diff --git a/adcc/backends/test_backends_polembed.py b/adcc/backends/test_backends_polembed.py index 3123f8bd..ec02e27a 100644 --- a/adcc/backends/test_backends_polembed.py +++ b/adcc/backends/test_backends_polembed.py @@ -31,9 +31,11 @@ from ..misc import expand_test_templates from .testing import cached_backend_hf -from ..testdata.cache import qchem_data +from ..testdata.cache import qchem_data, tmole_data from ..testdata.static_data import pe_potentials +from ..adc_pp.solvent import block_ph_ph_0_pe + try: import cppe # noqa: F401 @@ -52,7 +54,7 @@ @pytest.mark.skipif(len(backends) == 0, reason="No backend found.") @expand_test_templates(list(itertools.product(basissets, methods, backends))) class TestPolarizableEmbedding(unittest.TestCase): - def template_pe_formaldehyde(self, basis, method, backend): + def template_pe_perturbative_formaldehyde(self, basis, method, backend): basename = f"formaldehyde_{basis}_pe_{method}" qc_result = qchem_data[basename] pe_options = {"potfile": pe_potentials["fa_6w"]} @@ -82,3 +84,30 @@ def template_pe_formaldehyde(self, basis, method, backend): state.pe_ptlr_correction, atol=1e-5 ) + + def template_pe_coupling_formaldehyde(self, basis, method, backend): + if method != "adc2": + pytest.skip("") + basename = f"formaldehyde_{basis}_pe_{method}" + tm_result = tmole_data[basename] + pe_options = {"potfile": pe_potentials["fa_6w"]} + scfres = cached_backend_hf(backend, "formaldehyde", basis, + pe_options=pe_options) + assert_allclose(scfres.energy_scf, tm_result["energy_scf"], atol=1e-8) + + matrix = adcc.AdcMatrix( + method, scfres, + additional_blocks={'ph_ph': block_ph_ph_0_pe} + ) + + assert_allclose( + matrix.ground_state.energy(2), + tm_result["energy_mp2"], + atol=1e-8 + ) + state = adcc.run_adc(matrix, n_singlets=5, conv_tol=1e-7) + assert_allclose( + state.excitation_energy_uncorrected, + tm_result["excitation_energy"], + atol=1e-6 + ) diff --git a/adcc/testdata/cache.py b/adcc/testdata/cache.py index 447449fd..4dbfab89 100644 --- a/adcc/testdata/cache.py +++ b/adcc/testdata/cache.py @@ -248,12 +248,13 @@ def lists_to_ndarray(dictionary): return data -def read_qchem_data(): +def read_yaml_data(fname): thisdir = os.path.dirname(__file__) - yaml_file = os.path.join(thisdir, "qchem_dump.yml") + yaml_file = os.path.join(thisdir, fname) with open(yaml_file, "r") as f: data_raw = yaml.safe_load(f) return lists_to_ndarray(data_raw) -qchem_data = read_qchem_data() +qchem_data = read_yaml_data("qchem_dump.yml") +tmole_data = read_yaml_data("tmole_dump.yml") diff --git a/adcc/testdata/tmole_dump.yml b/adcc/testdata/tmole_dump.yml new file mode 100644 index 00000000..88e8660e --- /dev/null +++ b/adcc/testdata/tmole_dump.yml @@ -0,0 +1,14 @@ +formaldehyde_sto3g_pe_adc2: + basis: sto-3g + method: adc2 + molecule: formaldehyde + energy_scf: -112.36904349555 + energy_mp2: -112.4804685653 + excitation_energy: [0.1733632144, 0.3815240343, 0.5097618218, 0.5189443358, 0.5670575778] +formaldehyde_ccpvdz_pe_adc2: + basis: ccpvdz + method: adc2 + molecule: formaldehyde + energy_scf: -113.89976872224 + energy_mp2: -114.2156673890 + excitation_energy: [0.1596037963, 0.3123344811, 0.3621306427, 0.3795572252, 0.4089162982] \ No newline at end of file From c24c9dfaa01102f043804085016ddfdaea7b03f1 Mon Sep 17 00:00:00 2001 From: Maximilian Scheurer Date: Mon, 19 Apr 2021 16:48:33 +0200 Subject: [PATCH 09/21] composable AdcExtraTerm stuff refactor AdcExtraTerm and EnergyCorrection cleanup and improve coverage --- .clang-format | 89 ------------------- .github/workflows/ci.yaml | 2 +- adcc/AdcMatrix.py | 111 +++++++++++++++++++----- adcc/ExcitedStates.py | 91 +++++++++++++------ adcc/backends/psi4.py | 19 ++-- adcc/backends/pyscf.py | 19 ++-- adcc/backends/test_backends_polembed.py | 16 +++- adcc/backends/veloxchem.py | 19 ++-- adcc/test_AdcMatrix.py | 64 +++++++++++++- adcc/test_ExcitedStates.py | 30 ++++--- adcc/workflow.py | 8 ++ setup.py | 2 +- 12 files changed, 298 insertions(+), 172 deletions(-) delete mode 100644 .clang-format diff --git a/.clang-format b/.clang-format deleted file mode 100644 index d6d616bc..00000000 --- a/.clang-format +++ /dev/null @@ -1,89 +0,0 @@ -## --------------------------------------------------------------------- -## -## Dummy licence header -## -## --------------------------------------------------------------------- - -## DO NOT EDIT -## This file is automatically generated from a file in the repository "krims". -## Edit the original and call the script "update_from_sister_repos.sh" instead. - ---- -Language: Cpp -AccessModifierOffset: -1 -AlignAfterOpenBracket: Align -AlignConsecutiveAssignments: true -AlignConsecutiveDeclarations: false -AlignEscapedNewlinesLeft: true -AlignOperands: true -AlignTrailingComments: true -AllowAllParametersOfDeclarationOnNextLine: true -AllowShortBlocksOnASingleLine: false -AllowShortCaseLabelsOnASingleLine: false -AllowShortFunctionsOnASingleLine: All -AllowShortIfStatementsOnASingleLine: true -AllowShortLoopsOnASingleLine: true -AlwaysBreakAfterReturnType: None -AlwaysBreakBeforeMultilineStrings: true -AlwaysBreakTemplateDeclarations: true -BinPackArguments: true -BinPackParameters: true -BraceWrapping: - AfterClass: false - AfterControlStatement: false - AfterEnum: false - AfterFunction: false - AfterNamespace: false - AfterObjCDeclaration: false - AfterStruct: false - AfterUnion: false - BeforeCatch: false - BeforeElse: false - IndentBraces: false -BreakBeforeBinaryOperators: None -BreakBeforeBraces: Attach -BreakBeforeTernaryOperators: true -BreakConstructorInitializersBeforeComma: false -ColumnLimit: 90 -CommentPragmas: '^ IWYU pragma:' -ConstructorInitializerAllOnOneLineOrOnePerLine: true -ConstructorInitializerIndentWidth: 6 -ContinuationIndentWidth: 6 -Cpp11BracedListStyle: true -DerivePointerAlignment: false -DisableFormat: false -ExperimentalAutoDetectBinPacking: false -ForEachMacros: [ foreach, Q_FOREACH, BOOST_FOREACH ] -IndentCaseLabels: true -IndentWidth: 2 -IndentWrappedFunctionNames: false -KeepEmptyLinesAtTheStartOfBlocks: true -MaxEmptyLinesToKeep: 1 -NamespaceIndentation: None -ObjCBlockIndentWidth: 2 -ObjCSpaceAfterProperty: false -ObjCSpaceBeforeProtocolList: false -PenaltyBreakBeforeFirstCallParameter: 1 -PenaltyBreakComment: 300 -PenaltyBreakFirstLessLess: 120 -PenaltyBreakString: 1000 -PenaltyExcessCharacter: 1000000 -PenaltyReturnTypeOnItsOwnLine: 200 -PointerAlignment: Left -ReflowComments: true -SortIncludes: true -SpaceAfterCStyleCast: false -SpaceBeforeAssignmentOperators: true -SpaceBeforeParens: ControlStatements -SpaceInEmptyParentheses: false -SpacesBeforeTrailingComments: 2 -SpacesInAngles: false -SpacesInContainerLiterals: true -SpacesInCStyleCastParentheses: false -SpacesInParentheses: false -SpacesInSquareBrackets: false -Standard: Cpp11 -TabWidth: 8 -UseTab: Never -... - diff --git a/.github/workflows/ci.yaml b/.github/workflows/ci.yaml index 520876c6..029a3852 100644 --- a/.github/workflows/ci.yaml +++ b/.github/workflows/ci.yaml @@ -119,5 +119,5 @@ jobs: - name: Checking C++ files with clang-format working-directory: libadcc run: | - find . -type f \( -name '*.cc' -o -name '*.hh' \) -print0 | xargs -r0 clang-format --style=file -i + find . -type f \( -name '*.cc' -o -name '*.hh' \) -print0 | xargs -r0 clang-format -i git diff --exit-code diff --git a/adcc/AdcMatrix.py b/adcc/AdcMatrix.py index 593dc589..68b886ad 100644 --- a/adcc/AdcMatrix.py +++ b/adcc/AdcMatrix.py @@ -32,6 +32,48 @@ from .AmplitudeVector import AmplitudeVector +class AdcExtraTerm: + # NOTE: currently requires the matrix + # to allow for block construction with the same + # interface as the usual adc_pp, i.e., + # reference_state, ground_state, intermediates + def __init__(self, matrix, blocks): + """Initialise an AdcExtraTerm. + This class can be used to add customs terms + to an existing :py:class:`AdcMatrix` + + Parameters + ---------- + matrix : AdcMatrix + The matrix for which the extra term + should be created. + blocks : dict + A dictionary where the key labels the matrix block + and the item denotes a callable to construct + an :py:class:`AdcBlock` + """ + self.ground_state = matrix.ground_state + self.reference_state = matrix.reference_state + self.intermediates = matrix.intermediates + self.blocks = {} + if not isinstance(blocks, dict): + raise TypeError("blocks needs to be a dict.") + for space in blocks: + block_fun = blocks[space] + if not callable(block_fun): + raise TypeError("Items in additional_blocks must be callable.") + block = block_fun( + self.reference_state, self.ground_state, self.intermediates + ) + self.blocks[space] = block + + def __iadd__(self, other): + # disable in-place addition getting accidentally + # short-circuited to AdcMatrix.__radd__ + raise NotImplementedError("In-place addition not implemented " + " for AdcExtraTerm") + + class AdcMatrixlike: """ Base class marker for all objects like ADC matrices. @@ -50,8 +92,7 @@ class AdcMatrix(AdcMatrixlike): "adc3": dict(ph_ph=3, ph_pphh=2, pphh_ph=2, pphh_pphh=1), # noqa: E501 } - def __init__(self, method, hf_or_mp, block_orders=None, intermediates=None, - additional_blocks=None): + def __init__(self, method, hf_or_mp, block_orders=None, intermediates=None): """ Initialise an ADC matrix. @@ -85,6 +126,7 @@ def __init__(self, method, hf_or_mp, block_orders=None, intermediates=None, self.mospaces = hf_or_mp.reference_state.mospaces self.is_core_valence_separated = method.is_core_valence_separated self.ndim = 2 + self.extra_terms = [] self.intermediates = intermediates if self.intermediates is None: @@ -113,35 +155,64 @@ def __init__(self, method, hf_or_mp, block_orders=None, intermediates=None, # Build the blocks and diagonals with self.timer.record("build"): variant = None - if method.is_core_valence_separated: + if self.is_core_valence_separated: variant = "cvs" self.blocks_ph = { # TODO Rename to self.block in 0.16.0 block: ppmatrix.block(self.ground_state, block.split("_"), order=order, intermediates=self.intermediates, variant=variant) - for block, order in block_orders.items() if order is not None + for block, order in self.block_orders.items() if order is not None } - if additional_blocks is not None: - if not isinstance(additional_blocks, dict): - raise TypeError("additional_blocks needs to " - " be a dict.") - for space in additional_blocks: - if space not in self.blocks_ph: - raise ValueError("Can only add blocks" - " to existing matrix blocks.") - block_fun = additional_blocks[space] - if not callable(block_fun): - raise TypeError("Items in additional_blocks " - "must be callable.") - block = block_fun( - self.reference_state, self.ground_state, self.intermediates - ) - self.blocks_ph[space].add_block(block) self.__diagonal = sum(bl.diagonal for bl in self.blocks_ph.values() if bl.diagonal) self.__diagonal.evaluate() self.__init_space_data(self.__diagonal) + def __iadd__(self, other): + """In-place addition of an :py:class:`AdcExtraTerm` + + Parameters + ---------- + other : AdcExtraTerm + the extra term to be added + """ + assert isinstance(other, AdcExtraTerm) + assert all(k in self.blocks_ph for k in other.blocks) + for sp in other.blocks: + ob = other.blocks[sp] + self.blocks_ph[sp].add_block(ob) + diag = sum(bl.diagonal for bl in other.blocks.values() + if bl.diagonal) + # iadd does not work with numbers + self.__diagonal = self.__diagonal + diag + self.__diagonal.evaluate() + self.extra_terms.append(other) + return self + + def __add__(self, other): + """Addition of an :py:class:`AdcExtraTerm`, creating + a copy of self and adding the term to the new matrix + + Parameters + ---------- + other : AdcExtraTerm + the extra term to be added + + Returns + ------- + AdcMatrix + a copy of the AdcMatrix with the extra term added + """ + # NOTE: re-computes the (expensive) diagonal... + ret = AdcMatrix(self.method, self.ground_state, + block_orders=self.block_orders, + intermediates=self.intermediates) + ret += other + return ret + + def __radd__(self, other): + return self.__add__(other) + def __init_space_data(self, diagonal): """Update the cached data regarding the spaces of the ADC matrix""" self.axis_spaces = {} diff --git a/adcc/ExcitedStates.py b/adcc/ExcitedStates.py index 10164c04..377c2e9a 100644 --- a/adcc/ExcitedStates.py +++ b/adcc/ExcitedStates.py @@ -38,6 +38,19 @@ from .FormatDominantElements import FormatDominantElements +class EnergyCorrection: + def __init__(self, name, function, description=None): + assert isinstance(name, str) + assert callable(function) + self.name = name + self.function = function + self.description = description + + def __call__(self, excitation): + assert isinstance(excitation, Excitation) + return self.function(excitation) + + class FormatExcitationVector: def __init__(self, matrix, tolerance=0.01, index_format=None): """ @@ -127,8 +140,7 @@ def format(self, vector): class ExcitedStates(ElectronicTransition): - def __init__(self, data, method=None, property_method=None, - excitation_energy_corrections={}): + def __init__(self, data, method=None, property_method=None): """Construct an ExcitedStates class from some data obtained from an interative solver or another :class:`ExcitedStates` object. @@ -153,33 +165,46 @@ def __init__(self, data, method=None, property_method=None, property_method : str, optional Provide an explicit method for property calculations to override the automatic selection. - excitation_energy_corrections : dict, optional - Provide a dictionary of functions to compute corrections for - excitation energies, called for each excitation individually """ super().__init__(data, method, property_method) - - if hasattr(data, "excitation_energy_corrections"): - merged = data.excitation_energy_corrections.copy() - merged.update(excitation_energy_corrections) - excitation_energy_corrections = merged - - if hasattr(self.reference_state, "excitation_energy_corrections"): - merged = self.reference_state.excitation_energy_corrections.copy() - merged.update(excitation_energy_corrections) - excitation_energy_corrections = merged - - self.excitation_energy_corrections = excitation_energy_corrections - for key in self.excitation_energy_corrections: - corr_function = self.excitation_energy_corrections[key] - if not callable(corr_function): - raise TypeError("Elements in excitation_energy_corrections " - "must be callable.") - # call the function for exc. energy correction - correction = np.array([corr_function(exci) - for exci in self.excitations]) - setattr(self, key, correction) - self._excitation_energy += correction + self._excitation_energy_corrections = [] + # copy energy corrections if possible + # and avoids re-computation of the corrections + if hasattr(data, "_excitation_energy_corrections"): + for eec in data._excitation_energy_corrections: + correction_energy = getattr(data, eec.name) + setattr(self, eec.name, correction_energy) + self._excitation_energy += correction_energy + self._excitation_energy_corrections.append(eec) + + def __add_energy_correction(self, correction): + assert isinstance(correction, EnergyCorrection) + if hasattr(self, correction.name): + raise ValueError("ExcitedStates already has a correction" + f" with the name '{correction.name}'") + correction_energy = np.array([correction(exci) + for exci in self.excitations]) + setattr(self, correction.name, correction_energy) + self._excitation_energy += correction_energy + self._excitation_energy_corrections.append(correction) + + def __iadd__(self, other): + if isinstance(other, EnergyCorrection): + self.__add_energy_correction(other) + elif isinstance(other, list): + for k in other: + assert isinstance(k, EnergyCorrection) + self.__add_energy_correction(k) + else: + raise TypeError("Can only add EnergyCorrection (or list" + " of EnergyCorrection) to" + f" ExcitedState, not '{type(other)}'") + return self + + def __add__(self, other): + ret = ExcitedStates(self, self.method, self.property_method) + ret += other + return ret @cached_property @mark_excitation_property(transform_to_ao=True) @@ -338,6 +363,16 @@ def describe(self, oscillator_strengths=True, rotatory_strengths=False, text += body.format(i=i, ene=self.excitation_energy[i], ev=self.excitation_energy[i] * eV, **fields) text += separator + "\n" + if len(self._excitation_energy_corrections): + head_corr = "| excitation energy corrections included:" + text += head_corr + nspace = len(separator) - len(head_corr) - 1 + text += nspace * " " + "|\n" + maxlen = len(separator) - 8 + for eec in self._excitation_energy_corrections: + label = f"| - {eec.name:{maxlen}.{maxlen}}|\n" + text += label + text += separator + "\n" return text def _repr_pretty_(self, pp, cycle): @@ -398,7 +433,7 @@ def to_dataframe(self): Atomic units are used for all values. """ propkeys = self.excitation_property_keys - propkeys.extend(self.excitation_energy_corrections.keys()) + propkeys.extend([k.name for k in self._excitation_energy_corrections]) data = { "excitation": np.arange(0, self.size, dtype=int), "kind": np.tile(self.kind, self.size) diff --git a/adcc/backends/psi4.py b/adcc/backends/psi4.py index d4a6e031..6729b64f 100644 --- a/adcc/backends/psi4.py +++ b/adcc/backends/psi4.py @@ -29,6 +29,7 @@ from .EriBuilder import EriBuilder from ..exceptions import InvalidReference +from ..ExcitedStates import EnergyCorrection class Psi4OperatorIntegralProvider: @@ -108,12 +109,20 @@ def pe_energy(self, dm, elec_only=True): @property def excitation_energy_corrections(self): - ret = {} + ret = [] if hasattr(self.wfn, "pe_state"): - ret["pe_ptlr_correction"] = lambda view: \ - 2.0 * self.pe_energy(view.transition_dm_ao, elec_only=True) - ret["pe_ptss_correction"] = lambda view: \ - self.pe_energy(view.state_diffdm_ao, elec_only=True) + ptlr = EnergyCorrection( + "pe_ptlr_correction", + lambda view: 2.0 * self.pe_energy(view.transition_dm_ao, + elec_only=True) + ) + ptss = EnergyCorrection( + "pe_ptss_correction", + lambda view: self.pe_energy(view.state_diffdm_ao, + elec_only=True) + ) + ret.append(ptlr) + ret.append(ptss) return ret def get_backend(self): diff --git a/adcc/backends/pyscf.py b/adcc/backends/pyscf.py index ab62b486..ede80d08 100644 --- a/adcc/backends/pyscf.py +++ b/adcc/backends/pyscf.py @@ -27,6 +27,7 @@ from .EriBuilder import EriBuilder from ..exceptions import InvalidReference +from ..ExcitedStates import EnergyCorrection from pyscf import ao2mo, gto, scf, solvent @@ -126,13 +127,21 @@ def pe_energy(self, dm, elec_only=True): @property def excitation_energy_corrections(self): - ret = {} + ret = [] if hasattr(self.scfres, "with_solvent"): if isinstance(self.scfres.with_solvent, solvent.pol_embed.PolEmbed): - ret["pe_ptlr_correction"] = lambda view: \ - 2.0 * self.pe_energy(view.transition_dm_ao, elec_only=True) - ret["pe_ptss_correction"] = lambda view: \ - self.pe_energy(view.state_diffdm_ao, elec_only=True) + ptlr = EnergyCorrection( + "pe_ptlr_correction", + lambda view: 2.0 * self.pe_energy(view.transition_dm_ao, + elec_only=True) + ) + ptss = EnergyCorrection( + "pe_ptss_correction", + lambda view: self.pe_energy(view.state_diffdm_ao, + elec_only=True) + ) + ret.append(ptlr) + ret.append(ptss) return ret def get_backend(self): diff --git a/adcc/backends/test_backends_polembed.py b/adcc/backends/test_backends_polembed.py index ec02e27a..36054d46 100644 --- a/adcc/backends/test_backends_polembed.py +++ b/adcc/backends/test_backends_polembed.py @@ -34,6 +34,7 @@ from ..testdata.cache import qchem_data, tmole_data from ..testdata.static_data import pe_potentials +from ..AdcMatrix import AdcExtraTerm from ..adc_pp.solvent import block_ph_ph_0_pe try: @@ -62,6 +63,9 @@ def template_pe_perturbative_formaldehyde(self, basis, method, backend): pe_options=pe_options) state = adcc.run_adc(scfres, method=method, n_singlets=5, conv_tol=1e-10) + corrs = state.reference_state.excitation_energy_corrections + state += corrs + assert_allclose( qc_result["excitation_energy"], state.excitation_energy_uncorrected, @@ -95,10 +99,14 @@ def template_pe_coupling_formaldehyde(self, basis, method, backend): pe_options=pe_options) assert_allclose(scfres.energy_scf, tm_result["energy_scf"], atol=1e-8) - matrix = adcc.AdcMatrix( - method, scfres, - additional_blocks={'ph_ph': block_ph_ph_0_pe} - ) + matrix = adcc.AdcMatrix(method, scfres) + solvent = AdcExtraTerm(matrix, {'ph_ph': block_ph_ph_0_pe}) + + with pytest.raises(NotImplementedError): + solvent += matrix + + matrix += solvent + assert len(matrix.extra_terms) assert_allclose( matrix.ground_state.energy(2), diff --git a/adcc/backends/veloxchem.py b/adcc/backends/veloxchem.py index 3725acc3..92ca8370 100644 --- a/adcc/backends/veloxchem.py +++ b/adcc/backends/veloxchem.py @@ -32,6 +32,7 @@ from .EriBuilder import EriBuilder from ..exceptions import InvalidReference +from ..ExcitedStates import EnergyCorrection from veloxchem.mpitask import MpiTask from veloxchem.veloxchemlib import (AngularMomentumIntegralsDriver, @@ -130,12 +131,20 @@ def pe_energy(self, dm, elec_only=True): @property def excitation_energy_corrections(self): - ret = {} + ret = [] if hasattr(self.scfdrv, "pe_drv"): - ret["pe_ptlr_correction"] = lambda view: \ - 2.0 * self.pe_energy(view.transition_dm_ao, elec_only=True) - ret["pe_ptss_correction"] = lambda view: \ - self.pe_energy(view.state_diffdm_ao, elec_only=True) + ptlr = EnergyCorrection( + "pe_ptlr_correction", + lambda view: 2.0 * self.pe_energy(view.transition_dm_ao, + elec_only=True) + ) + ptss = EnergyCorrection( + "pe_ptss_correction", + lambda view: self.pe_energy(view.state_diffdm_ao, + elec_only=True) + ) + ret.append(ptlr) + ret.append(ptss) return ret def get_backend(self): diff --git a/adcc/test_AdcMatrix.py b/adcc/test_AdcMatrix.py index f90f87f4..7b57ca01 100644 --- a/adcc/test_AdcMatrix.py +++ b/adcc/test_AdcMatrix.py @@ -22,12 +22,14 @@ ## --------------------------------------------------------------------- import adcc import unittest +import pytest import itertools import numpy as np from numpy.testing import assert_allclose -from adcc.AdcMatrix import AdcMatrixShifted +from adcc.AdcMatrix import AdcMatrixShifted, AdcExtraTerm +from adcc.adc_pp.matrix import AdcBlock from adcc.testdata.cache import cache from .misc import expand_test_templates @@ -214,6 +216,66 @@ def test_matvec_adc2(self): diffv = resv.ph - refv.ph assert diffv.dot(diffv) < 1e-12 + def test_extra_term(self): + ground_state = adcc.LazyMp(cache.refstate["h2o_sto3g"]) + matrix = adcc.AdcMatrix("adc2", ground_state) + shift = -0.3 + shifted = AdcMatrixShifted(matrix, shift) + # TODO: need to do this to differentiate between + # diagonals for ph and pphh + # if we just pass numbers, i.e., shift + # we get 2*shift on the diagonal :( + ones = matrix.diagonal().ones_like() + + with pytest.raises(TypeError): + AdcBlock(0, 0) + with pytest.raises(TypeError): + AdcBlock(lambda x: x, "fail") + with pytest.raises(TypeError): + abc = AdcBlock(lambda x: x, ones) + abc.add_block("fail") + + with pytest.raises(TypeError): + AdcExtraTerm(matrix, "fail") + with pytest.raises(TypeError): + AdcExtraTerm(matrix, {"fail": "not_callable"}) + + def __shift_ph(hf, mp, intermediates): + def apply(invec): + return adcc.AmplitudeVector(ph=shift * invec.ph) + diag = adcc.AmplitudeVector(ph=shift * ones.ph) + return AdcBlock(apply, diag) + + def __shift_pphh(hf, mp, intermediates): + def apply(invec): + return adcc.AmplitudeVector(pphh=shift * invec.pphh) + diag = adcc.AmplitudeVector(pphh=shift * ones.pphh) + return AdcBlock(apply, diag) + extra = AdcExtraTerm( + matrix, {'ph_ph': __shift_ph, 'pphh_pphh': __shift_pphh} + ) + shifted_2 = matrix + extra + shifted_3 = extra + matrix + for manual in [shifted_2, shifted_3]: + assert_allclose( + shifted.diagonal().ph.to_ndarray(), + manual.diagonal().ph.to_ndarray(), + atol=1e-12 + ) + assert_allclose( + shifted.diagonal().pphh.to_ndarray(), + manual.diagonal().pphh.to_ndarray(), + atol=1e-12 + ) + vec = adcc.guess_zero(matrix) + vec.set_random() + ref = shifted @ vec + ret = manual @ vec + diff_s = ref.ph - ret.ph + diff_d = ref.pphh - ret.pphh + assert np.max(np.abs(diff_s.to_ndarray())) < 1e-12 + assert np.max(np.abs(diff_d.to_ndarray())) < 1e-12 + @expand_test_templates(testcases) class TestAdcMatrixShifted(unittest.TestCase): diff --git a/adcc/test_ExcitedStates.py b/adcc/test_ExcitedStates.py index 8827985e..b4d1e0b3 100644 --- a/adcc/test_ExcitedStates.py +++ b/adcc/test_ExcitedStates.py @@ -21,11 +21,12 @@ ## ## --------------------------------------------------------------------- import unittest +import pytest from numpy.testing import assert_allclose from adcc.testdata.cache import cache from adcc.OneParticleOperator import OneParticleOperator -from adcc.ExcitedStates import ExcitedStates +from adcc.ExcitedStates import EnergyCorrection from .test_state_densities import Runners @@ -64,12 +65,13 @@ def base_test(self, system, method, kind): method = method.replace("_", "-") state = cache.adc_states[system][method][kind] - corrections = { - "custom_correction1": lambda exci: exci.excitation_energy ** 2, - "custom_correction2": lambda exci: 2.0, - } - state_corrected = ExcitedStates(state, - excitation_energy_corrections=corrections) + cc1 = EnergyCorrection("custom_correction1", + lambda exci: exci.excitation_energy ** 2) + cc2 = EnergyCorrection("custom_correction2", + lambda exci: 2.0) + cc3 = EnergyCorrection("custom_correction3", + lambda exci: -42.0) + state_corrected = state + [cc1, cc2] for i in range(state.size): assert hasattr(state_corrected, "custom_correction1") assert hasattr(state_corrected, "custom_correction2") @@ -79,20 +81,22 @@ def base_test(self, system, method, kind): assert_allclose(state.excitation_energy[i] + corr, state_corrected.excitation_energy[i]) - corrections2 = {"custom_correction2": lambda exci: 3.0, - "custom_correction3": lambda exci: -42.0} - state_corrected2 = ExcitedStates( - state_corrected, excitation_energy_corrections=corrections2 - ) + with pytest.raises(ValueError): + state_corrected += cc2 + with pytest.raises(TypeError): + state_corrected += 1 + + state_corrected2 = state_corrected + cc3 for i in range(state.size): assert hasattr(state_corrected2, "custom_correction1") assert hasattr(state_corrected2, "custom_correction2") assert hasattr(state_corrected2, "custom_correction3") assert_allclose(state.excitation_energy[i], state_corrected2.excitation_energy_uncorrected[i]) - corr = state.excitation_energy[i] ** 2 + 3.0 - 42.0 + corr = state.excitation_energy[i] ** 2 + 2.0 - 42.0 assert_allclose(state.excitation_energy[i] + corr, state_corrected2.excitation_energy[i]) + state_corrected2.describe() class TestDataFrameExport(unittest.TestCase, Runners): diff --git a/adcc/workflow.py b/adcc/workflow.py index 37d78774..af4d3ce6 100644 --- a/adcc/workflow.py +++ b/adcc/workflow.py @@ -191,6 +191,10 @@ def run_adc(data_or_matrix, n_states=None, kind="any", conv_tol=None, if eigensolver is None: eigensolver = "davidson" + # add solvent coupling terms to matrix? + # matrix += pe_coupling_term + # solvent_scheme: scf, pt, ptlr, ptss, ... + diagres = diagonalise_adcmatrix( matrix, n_states, kind, guesses=guesses, n_guesses=n_guesses, n_guesses_doubles=n_guesses_doubles, conv_tol=conv_tol, output=output, @@ -198,6 +202,10 @@ def run_adc(data_or_matrix, n_states=None, kind="any", conv_tol=None, exstates = ExcitedStates(diagres) exstates.kind = kind exstates.spin_change = spin_change + + # if solvent_scheme == "ptss": + # ptss_correction = PerturbativeSolventCorrection("state-specific") + # exstates += ptss_correction return exstates diff --git a/setup.py b/setup.py index d4064629..a00779ea 100755 --- a/setup.py +++ b/setup.py @@ -460,7 +460,7 @@ def libadcc_extension(): def is_conda_build(): return ( os.environ.get("CONDA_BUILD", None) == "1" - or os.environ.get("CONDA_EXE", None) + or os.environ.get("CONDA_EXE", None) is not None ) From c3765bc74c35ed6fd296ef2d3eafc19aac40f77d Mon Sep 17 00:00:00 2001 From: Maximilian Scheurer Date: Thu, 22 Apr 2021 19:26:40 +0200 Subject: [PATCH 10/21] user interface for solvent schemes make adcc usable again, some cleanup update pe examples minor cleanup --- adcc/AdcMatrix.py | 12 ++- adcc/ExcitedStates.py | 3 +- adcc/ReferenceState.py | 3 +- adcc/backends/psi4.py | 15 +++- adcc/backends/pyscf.py | 36 +++++--- adcc/backends/test_backends_polembed.py | 27 +++++- adcc/backends/veloxchem.py | 6 +- adcc/test_AdcMatrix.py | 6 ++ adcc/workflow.py | 98 +++++++++++++++++++-- examples/pna/psi4_adc2_pna_6w_pol_embed.py | 8 +- examples/pna/pyscf_adc2_pna_6w_pol_embed.py | 25 +++--- 11 files changed, 187 insertions(+), 52 deletions(-) diff --git a/adcc/AdcMatrix.py b/adcc/AdcMatrix.py index 68b886ad..84f67112 100644 --- a/adcc/AdcMatrix.py +++ b/adcc/AdcMatrix.py @@ -176,8 +176,16 @@ def __iadd__(self, other): other : AdcExtraTerm the extra term to be added """ - assert isinstance(other, AdcExtraTerm) - assert all(k in self.blocks_ph for k in other.blocks) + if isinstance(other, list): + for k in other: + self += k + return self + if not isinstance(other, AdcExtraTerm): + raise TypeError("Can only add AdcExtra term or" + " lists thereof to AdcMatrix.") + if not all(k in self.blocks_ph for k in other.blocks): + raise ValueError("Can only add to blocks of" + " AdcMatrix that already exist.") for sp in other.blocks: ob = other.blocks[sp] self.blocks_ph[sp].add_block(ob) diff --git a/adcc/ExcitedStates.py b/adcc/ExcitedStates.py index 377c2e9a..ef6d62b5 100644 --- a/adcc/ExcitedStates.py +++ b/adcc/ExcitedStates.py @@ -193,8 +193,7 @@ def __iadd__(self, other): self.__add_energy_correction(other) elif isinstance(other, list): for k in other: - assert isinstance(k, EnergyCorrection) - self.__add_energy_correction(k) + self += k else: raise TypeError("Can only add EnergyCorrection (or list" " of EnergyCorrection) to" diff --git a/adcc/ReferenceState.py b/adcc/ReferenceState.py index a74e79c8..14c2b343 100644 --- a/adcc/ReferenceState.py +++ b/adcc/ReferenceState.py @@ -155,7 +155,8 @@ def __init__(self, hfdata, core_orbitals=None, frozen_core=None, self.orbital_coefficients, self.conv_tol ) - for name in ["excitation_energy_corrections"]: + self.solvent = None # no solvent attached by default + for name in ["excitation_energy_corrections", "solvent"]: if hasattr(hfdata, name): setattr(self, name, getattr(hfdata, name)) diff --git a/adcc/backends/psi4.py b/adcc/backends/psi4.py index 6729b64f..2c7a1287 100644 --- a/adcc/backends/psi4.py +++ b/adcc/backends/psi4.py @@ -109,8 +109,8 @@ def pe_energy(self, dm, elec_only=True): @property def excitation_energy_corrections(self): - ret = [] - if hasattr(self.wfn, "pe_state"): + ret = {} + if self.solvent == "pe": ptlr = EnergyCorrection( "pe_ptlr_correction", lambda view: 2.0 * self.pe_energy(view.transition_dm_ao, @@ -121,8 +121,15 @@ def excitation_energy_corrections(self): lambda view: self.pe_energy(view.state_diffdm_ao, elec_only=True) ) - ret.append(ptlr) - ret.append(ptss) + ret["pe_ptlr_correction"] = ptlr + ret["pe_ptss_correction"] = ptss + return ret + + @property + def solvent(self): + ret = None + if hasattr(self.wfn, "pe_state"): + ret = "pe" return ret def get_backend(self): diff --git a/adcc/backends/pyscf.py b/adcc/backends/pyscf.py index ede80d08..f019b09f 100644 --- a/adcc/backends/pyscf.py +++ b/adcc/backends/pyscf.py @@ -127,21 +127,31 @@ def pe_energy(self, dm, elec_only=True): @property def excitation_energy_corrections(self): - ret = [] + ret = {} + if self.solvent == "pe": + ptlr = EnergyCorrection( + "pe_ptlr_correction", + lambda view: 2.0 * self.pe_energy(view.transition_dm_ao, + elec_only=True) + ) + ptss = EnergyCorrection( + "pe_ptss_correction", + lambda view: self.pe_energy(view.state_diffdm_ao, + elec_only=True) + ) + # NOTE: I don't like the duplicate 'name', + # but this way it's easier to exctract the corrections + # directly + ret["pe_ptlr_correction"] = ptlr + ret["pe_ptss_correction"] = ptss + return ret + + @property + def solvent(self): + ret = None if hasattr(self.scfres, "with_solvent"): if isinstance(self.scfres.with_solvent, solvent.pol_embed.PolEmbed): - ptlr = EnergyCorrection( - "pe_ptlr_correction", - lambda view: 2.0 * self.pe_energy(view.transition_dm_ao, - elec_only=True) - ) - ptss = EnergyCorrection( - "pe_ptss_correction", - lambda view: self.pe_energy(view.state_diffdm_ao, - elec_only=True) - ) - ret.append(ptlr) - ret.append(ptss) + ret = "pe" return ret def get_backend(self): diff --git a/adcc/backends/test_backends_polembed.py b/adcc/backends/test_backends_polembed.py index 36054d46..081b79fc 100644 --- a/adcc/backends/test_backends_polembed.py +++ b/adcc/backends/test_backends_polembed.py @@ -31,6 +31,7 @@ from ..misc import expand_test_templates from .testing import cached_backend_hf +from ..exceptions import InputError from ..testdata.cache import qchem_data, tmole_data from ..testdata.static_data import pe_potentials @@ -62,9 +63,8 @@ def template_pe_perturbative_formaldehyde(self, basis, method, backend): scfres = cached_backend_hf(backend, "formaldehyde", basis, pe_options=pe_options) state = adcc.run_adc(scfres, method=method, - n_singlets=5, conv_tol=1e-10) - corrs = state.reference_state.excitation_energy_corrections - state += corrs + n_singlets=5, conv_tol=1e-10, + solvent_scheme=["ptlr", "ptss"]) assert_allclose( qc_result["excitation_energy"], @@ -105,6 +105,7 @@ def template_pe_coupling_formaldehyde(self, basis, method, backend): with pytest.raises(NotImplementedError): solvent += matrix + # manually add the coupling term matrix += solvent assert len(matrix.extra_terms) @@ -113,7 +114,25 @@ def template_pe_coupling_formaldehyde(self, basis, method, backend): tm_result["energy_mp2"], atol=1e-8 ) - state = adcc.run_adc(matrix, n_singlets=5, conv_tol=1e-7) + state = adcc.run_adc(matrix, n_singlets=5, conv_tol=1e-7, + solvent_scheme="hf") + assert_allclose( + state.excitation_energy_uncorrected, + tm_result["excitation_energy"], + atol=1e-6 + ) + + # invalid combination + with pytest.raises(InputError): + adcc.run_adc(scfres, method=method, n_singlets=5, + solvent_scheme=["postscf", "ptlr"]) + # no scheme 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, solvent_scheme="postscf") assert_allclose( state.excitation_energy_uncorrected, tm_result["excitation_energy"], diff --git a/adcc/backends/veloxchem.py b/adcc/backends/veloxchem.py index 92ca8370..ef08fec6 100644 --- a/adcc/backends/veloxchem.py +++ b/adcc/backends/veloxchem.py @@ -131,7 +131,7 @@ def pe_energy(self, dm, elec_only=True): @property def excitation_energy_corrections(self): - ret = [] + ret = {} if hasattr(self.scfdrv, "pe_drv"): ptlr = EnergyCorrection( "pe_ptlr_correction", @@ -143,8 +143,8 @@ def excitation_energy_corrections(self): lambda view: self.pe_energy(view.state_diffdm_ao, elec_only=True) ) - ret.append(ptlr) - ret.append(ptss) + ret["pe_ptlr_correction"] = ptlr + ret["pe_ptss_correction"] = ptss return ret def get_backend(self): diff --git a/adcc/test_AdcMatrix.py b/adcc/test_AdcMatrix.py index 7b57ca01..c4fb8eb8 100644 --- a/adcc/test_AdcMatrix.py +++ b/adcc/test_AdcMatrix.py @@ -218,6 +218,9 @@ def test_matvec_adc2(self): def test_extra_term(self): ground_state = adcc.LazyMp(cache.refstate["h2o_sto3g"]) + matrix_adc1 = adcc.AdcMatrix("adc1", ground_state) + with pytest.raises(TypeError): + matrix_adc1 += 42 matrix = adcc.AdcMatrix("adc2", ground_state) shift = -0.3 shifted = AdcMatrixShifted(matrix, shift) @@ -254,6 +257,9 @@ def apply(invec): extra = AdcExtraTerm( matrix, {'ph_ph': __shift_ph, 'pphh_pphh': __shift_pphh} ) + # cannot add to 'pphh_pphh' in ADC(1) matrix + with pytest.raises(ValueError): + matrix_adc1 += extra shifted_2 = matrix + extra shifted_3 = extra + matrix for manual in [shifted_2, shifted_3]: diff --git a/adcc/workflow.py b/adcc/workflow.py index af4d3ce6..9f7db3ef 100644 --- a/adcc/workflow.py +++ b/adcc/workflow.py @@ -29,7 +29,7 @@ from .guess import (guesses_any, guesses_singlet, guesses_spin_flip, guesses_triplet) from .LazyMp import LazyMp -from .AdcMatrix import AdcMatrix, AdcMatrixlike +from .AdcMatrix import AdcMatrix, AdcMatrixlike, AdcExtraTerm from .AdcMethod import AdcMethod from .exceptions import InputError from .ExcitedStates import ExcitedStates @@ -47,6 +47,7 @@ def run_adc(data_or_matrix, n_states=None, kind="any", conv_tol=None, n_guesses_doubles=None, output=sys.stdout, core_orbitals=None, frozen_core=None, frozen_virtual=None, method=None, n_singlets=None, n_triplets=None, n_spin_flip=None, + solvent_scheme=None, **solverargs): """Run an ADC calculation. @@ -191,9 +192,11 @@ def run_adc(data_or_matrix, n_states=None, kind="any", conv_tol=None, if eigensolver is None: eigensolver = "davidson" - # add solvent coupling terms to matrix? - # matrix += pe_coupling_term - # solvent_scheme: scf, pt, ptlr, ptss, ... + # Setup solvent coupling terms and energy corrections + solvent_matrix_terms, solvent_energy_corrections = \ + setup_solvent(matrix, solvent_scheme) + # add terms to matrix + matrix += solvent_matrix_terms diagres = diagonalise_adcmatrix( matrix, n_states, kind, guesses=guesses, n_guesses=n_guesses, @@ -203,9 +206,8 @@ def run_adc(data_or_matrix, n_states=None, kind="any", conv_tol=None, exstates.kind = kind exstates.spin_change = spin_change - # if solvent_scheme == "ptss": - # ptss_correction = PerturbativeSolventCorrection("state-specific") - # exstates += ptss_correction + # add corrections to excited states + exstates += solvent_energy_corrections return exstates @@ -503,3 +505,85 @@ def setup_solver_printing(solmethod_name, matrix, kind, default_print, def inner_callback(state, identifier): default_print(state, identifier, output) return inner_callback + + +def setup_solvent(matrix, solvent_scheme): + """ + Setup solvent matrix terms and/or energy corrections. + Internal function called from run_adc. + """ + # TODO: move someplace meaningful & reasonable... + valid_solvent_schemes = { + "hf": """ + only couple via the 'solvated' orbitals of the HF reference + state, no additional matrix terms or perturbative corrections + are used automatically + """, + "ptss": """ + perturbative state-specific (ptSS) correction, computed based on + the difference density between the ground and excited state + """, + "ptlr": """ + perturbative linear-response (ptLR) correction, computed based on + the transition density between the ground and excited state + """, + # NOTE: could also be called 'lr'... + "postscf": """ + iterative coupling to the solvent via a CIS-like coupling + density matrix, the term is added to the ADC matrix + """, + } + valid_scheme_names = list(valid_solvent_schemes.keys()) + + hf = matrix.reference_state + if hf.solvent and not solvent_scheme: + raise InputError( + "Solvent found in reference state, but no solvent_scheme" + " specified. Please select one or more of the following" + f" schemes: {valid_solvent_schemes.keys()}." + ) + elif solvent_scheme and not hf.solvent: + raise InputError( + "solvent_scheme specified, but no solvent method" + " was found in reference state." + ) + elif solvent_scheme is None and hf.solvent is None: + # nothing to do + return [], [] + if not isinstance(solvent_scheme, list): + solvent_scheme = [solvent_scheme] + if any(scheme not in valid_solvent_schemes for scheme in solvent_scheme): + raise InputError("Invalid solvent_scheme given." + f" Valid schemes are: {valid_scheme_names}.") + if "hf" in solvent_scheme: + if len(solvent_scheme) > 1: + raise InputError("hf solvent_scheme cannot be combined" + " with other schemes.") + return [], [] + if "postscf" in solvent_scheme and "ptlr" in solvent_scheme: + raise InputError("Cannot combine ptlr correction with iterative" + " postscf linear response procedure.") + + extra_matrix_terms = [] + energy_corrections = [] + + pt_schemes = [pt for pt in ["ptss", "ptlr"] if pt in solvent_scheme] + for pt in pt_schemes: + hf_corr = hf.excitation_energy_corrections + eec_key = f"{hf.solvent}_{pt}_correction" + if eec_key not in hf_corr: + raise ValueError(f"{pt} correction requested, but could not find" + f" the needed function {eec_key} in" + f" reference state from backend {hf.backend}.") + energy_corrections.append(hf_corr[eec_key]) + if "postscf" in solvent_scheme: + from adcc.adc_pp import solvent as solvent_mat_terms + block_key = f"block_ph_ph_0_{hf.solvent}" + if not hasattr(solvent_mat_terms, block_key): + raise NotImplementedError("Matrix term for postscf coupling" + f" with solvent {hf.solvent}" + " not implemented.") + block_fun = getattr(solvent_mat_terms, block_key) + extra_matrix_terms.append(AdcExtraTerm(matrix, {'ph_ph': block_fun})) + + return extra_matrix_terms, energy_corrections diff --git a/examples/pna/psi4_adc2_pna_6w_pol_embed.py b/examples/pna/psi4_adc2_pna_6w_pol_embed.py index 080fbd26..6102b8f8 100644 --- a/examples/pna/psi4_adc2_pna_6w_pol_embed.py +++ b/examples/pna/psi4_adc2_pna_6w_pol_embed.py @@ -29,10 +29,7 @@ no_com """) -# set the number of cores equal to the auto-determined value from -# the adcc ThreadPool psi4.core.set_num_threads(4) -# psi4.core.be_quiet() psi4.set_options({'basis': "sto-3g", 'scf_type': 'pk', 'pe': 'true', @@ -42,5 +39,6 @@ scf_e, wfn = psi4.energy('SCF', return_wfn=True) # Run an adc2 calculation: -state = adcc.adc2(wfn, n_singlets=5, conv_tol=1e-8) -print(state.pe_ptss_energy_correction) +state = adcc.adc2(wfn, n_singlets=5, conv_tol=1e-8, + solvent_scheme=['ptss', 'ptlr']) +print(state.describe()) diff --git a/examples/pna/pyscf_adc2_pna_6w_pol_embed.py b/examples/pna/pyscf_adc2_pna_6w_pol_embed.py index 5dae2738..738f4d55 100644 --- a/examples/pna/pyscf_adc2_pna_6w_pol_embed.py +++ b/examples/pna/pyscf_adc2_pna_6w_pol_embed.py @@ -2,7 +2,6 @@ ## vi: tabstop=4 shiftwidth=4 softtabstop=4 expandtab import adcc -import cppe from pyscf import gto, scf from pyscf.solvent import PE @@ -30,18 +29,22 @@ """, basis='sto-3g', ) -pe_options = cppe.PeOptions() -pe_options.potfile = "pna_6w.pot" -# pe_options.iso_pol = True -scfres = PE(scf.RHF(mol), pe_options) -scfres.conv_tol = 1e-10 -scfres.conv_tol_grad = 1e-10 +scfres = PE(scf.RHF(mol), {"potfile": "pna_6w.pot"}) +scfres.conv_tol = 1e-8 +scfres.conv_tol_grad = 1e-6 scfres.max_cycle = 250 scfres.kernel() -state = adcc.adc2(scfres, n_singlets=5, conv_tol=1e-8) +# model the solvent through perturbative corrections +state_pt = adcc.adc2(scfres, n_singlets=5, conv_tol=1e-5, + solvent_scheme=['ptss', 'ptlr']) -ptlr = state.pe_ptlr_correction -ptss = state.pe_ptss_correction -print(ptlr, ptss) +# now model the solvent through linear-response/postscf coupling +# in the ADC matrix, re-using the matrix from previous run. +# This will modify state_pt.matrix +state_lr = adcc.run_adc(state_pt.matrix, n_singlets=5, conv_tol=1e-5, + solvent_scheme='postscf') + +print(state_pt.describe()) +print(state_lr.describe()) From c87aebd944175c11a6e190be78ea3e1656e0dcd3 Mon Sep 17 00:00:00 2001 From: Maximilian Scheurer Date: Tue, 27 Apr 2021 11:34:23 +0200 Subject: [PATCH 11/21] some revs refactor pe induction operator more revs revert adcblock, patch function --- adcc/AdcMatrix.py | 39 ++++++------ adcc/ExcitedStates.py | 21 +++++-- adcc/OperatorIntegrals.py | 84 +++++++++++++++++-------- adcc/adc_pp/matrix.py | 57 +++-------------- adcc/adc_pp/solvent.py | 2 +- adcc/backends/psi4.py | 11 ++-- adcc/backends/pyscf.py | 12 ++-- adcc/backends/test_backends_polembed.py | 5 +- adcc/test_AdcMatrix.py | 8 --- adcc/workflow.py | 13 ++-- 10 files changed, 122 insertions(+), 130 deletions(-) diff --git a/adcc/AdcMatrix.py b/adcc/AdcMatrix.py index 84f67112..1a8a10ee 100644 --- a/adcc/AdcMatrix.py +++ b/adcc/AdcMatrix.py @@ -67,12 +67,6 @@ def __init__(self, matrix, blocks): ) self.blocks[space] = block - def __iadd__(self, other): - # disable in-place addition getting accidentally - # short-circuited to AdcMatrix.__radd__ - raise NotImplementedError("In-place addition not implemented " - " for AdcExtraTerm") - class AdcMatrixlike: """ @@ -157,13 +151,15 @@ def __init__(self, method, hf_or_mp, block_orders=None, intermediates=None): variant = None if self.is_core_valence_separated: variant = "cvs" - self.blocks_ph = { # TODO Rename to self.block in 0.16.0 + blocks = { block: ppmatrix.block(self.ground_state, block.split("_"), order=order, intermediates=self.intermediates, variant=variant) for block, order in self.block_orders.items() if order is not None } - self.__diagonal = sum(bl.diagonal for bl in self.blocks_ph.values() + # TODO Rename to self.block in 0.16.0 + self.blocks_ph = {bl: blocks[bl].apply for bl in blocks} + self.__diagonal = sum(bl.diagonal for bl in blocks.values() if bl.diagonal) self.__diagonal.evaluate() self.__init_space_data(self.__diagonal) @@ -176,23 +172,22 @@ def __iadd__(self, other): other : AdcExtraTerm the extra term to be added """ - if isinstance(other, list): - for k in other: - self += k - return self if not isinstance(other, AdcExtraTerm): - raise TypeError("Can only add AdcExtra term or" - " lists thereof to AdcMatrix.") + return NotImplemented if not all(k in self.blocks_ph for k in other.blocks): raise ValueError("Can only add to blocks of" " AdcMatrix that already exist.") for sp in other.blocks: - ob = other.blocks[sp] - self.blocks_ph[sp].add_block(ob) - diag = sum(bl.diagonal for bl in other.blocks.values() - if bl.diagonal) + orig_app = self.blocks_ph[sp] + other_app = other.blocks[sp].apply + + def patched_apply(ampl, original=orig_app, other=other_app): + return sum(app(ampl) for app in (original, other)) + self.blocks_ph[sp] = patched_apply + other_diagonal = sum(bl.diagonal for bl in other.blocks.values() + if bl.diagonal) # iadd does not work with numbers - self.__diagonal = self.__diagonal + diag + self.__diagonal = self.__diagonal + other_diagonal self.__diagonal.evaluate() self.extra_terms.append(other) return self @@ -211,6 +206,8 @@ def __add__(self, other): AdcMatrix a copy of the AdcMatrix with the extra term added """ + if not isinstance(other, AdcExtraTerm): + return NotImplemented # NOTE: re-computes the (expensive) diagonal... ret = AdcMatrix(self.method, self.ground_state, block_orders=self.block_orders, @@ -304,7 +301,7 @@ def block_apply(self, block, tensor): with self.timer.record(f"apply/{block}"): outblock, inblock = block.split("_") ampl = AmplitudeVector(**{inblock: tensor}) - ret = self.blocks_ph[block].apply(ampl) + ret = self.blocks_ph[block](ampl) return getattr(ret, outblock) @timed_member_call() @@ -313,7 +310,7 @@ def matvec(self, v): Compute the matrix-vector product of the ADC matrix with an excitation amplitude and return the result. """ - return sum(block.apply(v) for block in self.blocks_ph.values()) + return sum(block(v) for block in self.blocks_ph.values()) def rmatvec(self, v): # ADC matrix is symmetric diff --git a/adcc/ExcitedStates.py b/adcc/ExcitedStates.py index ef6d62b5..3ce9e641 100644 --- a/adcc/ExcitedStates.py +++ b/adcc/ExcitedStates.py @@ -39,12 +39,25 @@ class EnergyCorrection: - def __init__(self, name, function, description=None): - assert isinstance(name, str) - assert callable(function) + def __init__(self, name, function): + """A helper class to represent excitation energy + corrections. + + Parameters + ---------- + name : str + descriptive name of the energy correction + function : callable + function that takes a :py:class:`Excitation` + as single argument and returns the energy + correction as a float + """ + if not isinstance(name, str): + raise TypeError("name needs to be a string.") + if not callable(function): + raise TypeError("function needs to be callable.") self.name = name self.function = function - self.description = description def __call__(self, excitation): assert isinstance(excitation, Excitation) diff --git a/adcc/OperatorIntegrals.py b/adcc/OperatorIntegrals.py index 1f0be057..270529eb 100644 --- a/adcc/OperatorIntegrals.py +++ b/adcc/OperatorIntegrals.py @@ -32,16 +32,21 @@ def transform_operator_ao2mo(tensor_bb, tensor_ff, coefficients, conv_tol=1e-14): - """ - Take a block-diagonal tensor in the atomic orbital basis + """Take a block-diagonal tensor in the atomic orbital basis and transform it into the molecular orbital basis in the convention used by adcc. - @param tensor_bb Block-diagonal tensor in the atomic orbital basis - @param tensor_ff Output tensor with the symmetry set-up to contain - the operator in the molecular orbital representation - @param coefficients Function providing coefficient blocks - @param conv_tol SCF convergence tolerance + Parameters + ---------- + tensor_bb : Tensor + Block-diagonal tensor in the atomic orbital basis + tensor_ff : Tensor + Output tensor with the symmetry set-up to contain + the operator in the molecular orbital representation + coefficients : callable + Function providing coefficient blocks + conv_tol : float, optional + SCF convergence tolerance, by default 1e-14 """ for blk in tensor_ff.blocks: assert len(blk) == 4 @@ -93,8 +98,13 @@ def provider_ao(self): @cached_property def available(self): """Which integrals are available in the underlying backend""" - return [integral - for integral in ("electric_dipole", "magnetic_dipole", "nabla") + integrals = ( + "electric_dipole", + "magnetic_dipole", + "nabla", + "pe_induction_elec" + ) + return [integral for integral in integrals if hasattr(self.provider_ao, integral)] def import_dipole_like_operator(self, integral, is_symmetric=True): @@ -137,25 +147,45 @@ def nabla(self): """ return self.import_dipole_like_operator("nabla", is_symmetric=False) + def __import_density_dependent_operator(self, ao_callback, is_symmetric=True): + """Returns a function that imports a density-dependent operator. + The returned function imports the operator and transforms it to the + molecular orbital basis. + + Parameters + ---------- + ao_callback : callable + Function that computes the operator in atomic orbitals using + a :py:class:`OneParticleOperator` (the density matrix + in atomic orbitals) as single argument + is_symmetric : bool, optional + if the imported operator is symmetric, by default True + """ + def process_operator(dm, callback=ao_callback, is_symmetric=is_symmetric): + dm_ao = sum(dm.to_ao_basis()) + v_ao = callback(dm_ao) + v_bb = replicate_ao_block( + self.mospaces, v_ao, is_symmetric=is_symmetric + ) + v_ff = OneParticleOperator(self.mospaces, is_symmetric=is_symmetric) + transform_operator_ao2mo( + v_bb, v_ff, self.__coefficients, self.__conv_tol + ) + return v_ff + return process_operator + @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 + def pe_induction_elec(self): + """ + Returns a function to obtain the (density-dependent) + PE electronic induction operator in the molecular orbital basis + """ + if "pe_induction_elec" not in self.available: + raise NotImplementedError("PE electronic induction operator " + "not implemented " + f"in {self.provider_ao.backend} backend.") + callback = self.provider_ao.pe_induction_elec + return self.__import_density_dependent_operator(callback) @property def timer(self): diff --git a/adcc/adc_pp/matrix.py b/adcc/adc_pp/matrix.py index cbafba54..5ef0b4b1 100644 --- a/adcc/adc_pp/matrix.py +++ b/adcc/adc_pp/matrix.py @@ -21,6 +21,7 @@ ## ## --------------------------------------------------------------------- from math import sqrt +from collections import namedtuple from adcc import block as b from adcc.functions import direct_sum, einsum, zeros_like @@ -38,53 +39,15 @@ # really so much our focus. -class AdcBlock: - def __init__(self, apply, diagonal): - """AdcBlock contains the matrix apply and diagonal - routines for a specific ADC matrix block and is used to compose - the dispatch routines for :py:class:`AdcMatrix` - - Parameters - ---------- - apply : callable - function mapping an AmplitudeVector to the contribution of - this block to the result of applying the ADC matrix - diagonal : AmplitudeVector, int, float - expression to the diagonal of the ADC matrix from this block - """ - if not callable(apply): - raise TypeError("apply needs to be callable.") - if not isinstance(diagonal, (AmplitudeVector, int, float)): - raise TypeError("diagonal needs to be an AmplitudeVector or a number.") - self._applies = [apply] - self._diagonals = [diagonal] - - def add_block(self, other): - """Adds another :py:class:`AdcBlock` to this AdcBlock. - - Parameters - ---------- - other : AdcBlock - block to be added - """ - if not isinstance(other, AdcBlock): - raise TypeError("other must be of type AdcBlock.") - self._applies.append(other.apply) - self._diagonals.append(other.diagonal) - - def apply(self, ampl): - """Applies the block to an input AmplitudeVector - - Parameters - ---------- - ampl : AmplitudeVector - Input AmplitudeVector - """ - return sum(app(ampl) for app in self._applies) - - @property - def diagonal(self): - return sum(self._diagonals) +# +# Dispatch routine +# +""" +`apply` is a function mapping an AmplitudeVector to the contribution of this +block to the result of applying the ADC matrix. `diagonal` is an `AmplitudeVector` +containing the expression to the diagonal of the ADC matrix from this block. +""" +AdcBlock = namedtuple("AdcBlock", ["apply", "diagonal"]) def block(ground_state, spaces, order, variant=None, intermediates=None): diff --git a/adcc/adc_pp/solvent.py b/adcc/adc_pp/solvent.py index c7624c6d..e8d084e1 100644 --- a/adcc/adc_pp/solvent.py +++ b/adcc/adc_pp/solvent.py @@ -37,6 +37,6 @@ def block_ph_ph_0_pe(hf, mp, intermediates): def apply(ampl): tdm = OneParticleOperator(mp, is_symmetric=False) tdm.vo = ampl.ph.transpose() - vpe = op.density_dependent_operators["pe_induction_elec"](tdm) + vpe = op.pe_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 2c7a1287..a832e7bf 100644 --- a/adcc/backends/psi4.py +++ b/adcc/backends/psi4.py @@ -55,15 +55,14 @@ def nabla(self): return [-1.0 * np.asarray(comp) for comp in self.mints.ao_nabla()] @property - def density_dependent_operators(self): - ret = {} + def pe_induction_elec(self): if hasattr(self.wfn, "pe_state"): - ret["pe_induction_elec"] = lambda dm: \ - self.wfn.pe_state.get_pe_contribution( + def pe_induction_elec_ao(dm): + return self.wfn.pe_state.get_pe_contribution( psi4.core.Matrix.from_array(dm.to_ndarray()), elec_only=True - )[1] - return ret + )[1] + return pe_induction_elec_ao class Psi4EriBuilder(EriBuilder): diff --git a/adcc/backends/pyscf.py b/adcc/backends/pyscf.py index f019b09f..44ba3fc5 100644 --- a/adcc/backends/pyscf.py +++ b/adcc/backends/pyscf.py @@ -57,14 +57,14 @@ def nabla(self): ) @property - def density_dependent_operators(self): - ret = {} + def pe_induction_elec(self): if hasattr(self.scfres, "with_solvent"): if isinstance(self.scfres.with_solvent, solvent.pol_embed.PolEmbed): - ret["pe_induction_elec"] = lambda dm: \ - self.scfres.with_solvent._exec_cppe(dm.to_ndarray(), - elec_only=True)[1] - return ret + 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 # TODO: refactor ERI builder to be more general diff --git a/adcc/backends/test_backends_polembed.py b/adcc/backends/test_backends_polembed.py index 081b79fc..deae214a 100644 --- a/adcc/backends/test_backends_polembed.py +++ b/adcc/backends/test_backends_polembed.py @@ -91,7 +91,7 @@ def template_pe_perturbative_formaldehyde(self, basis, method, backend): def template_pe_coupling_formaldehyde(self, basis, method, backend): if method != "adc2": - pytest.skip("") + pytest.skip("Reference only exists for adc2.") basename = f"formaldehyde_{basis}_pe_{method}" tm_result = tmole_data[basename] pe_options = {"potfile": pe_potentials["fa_6w"]} @@ -102,9 +102,6 @@ def template_pe_coupling_formaldehyde(self, basis, method, backend): matrix = adcc.AdcMatrix(method, scfres) solvent = AdcExtraTerm(matrix, {'ph_ph': block_ph_ph_0_pe}) - with pytest.raises(NotImplementedError): - solvent += matrix - # manually add the coupling term matrix += solvent assert len(matrix.extra_terms) diff --git a/adcc/test_AdcMatrix.py b/adcc/test_AdcMatrix.py index c4fb8eb8..523e333c 100644 --- a/adcc/test_AdcMatrix.py +++ b/adcc/test_AdcMatrix.py @@ -230,14 +230,6 @@ def test_extra_term(self): # we get 2*shift on the diagonal :( ones = matrix.diagonal().ones_like() - with pytest.raises(TypeError): - AdcBlock(0, 0) - with pytest.raises(TypeError): - AdcBlock(lambda x: x, "fail") - with pytest.raises(TypeError): - abc = AdcBlock(lambda x: x, ones) - abc.add_block("fail") - with pytest.raises(TypeError): AdcExtraTerm(matrix, "fail") with pytest.raises(TypeError): diff --git a/adcc/workflow.py b/adcc/workflow.py index 9f7db3ef..8fdaafb1 100644 --- a/adcc/workflow.py +++ b/adcc/workflow.py @@ -193,10 +193,11 @@ def run_adc(data_or_matrix, n_states=None, kind="any", conv_tol=None, eigensolver = "davidson" # Setup solvent coupling terms and energy corrections - solvent_matrix_terms, solvent_energy_corrections = \ - setup_solvent(matrix, solvent_scheme) + ret = setup_solvent(matrix, solvent_scheme) + solvent_matrix_term, solvent_energy_corrections = ret # add terms to matrix - matrix += solvent_matrix_terms + if solvent_matrix_term: + matrix += solvent_matrix_term diagres = diagonalise_adcmatrix( matrix, n_states, kind, guesses=guesses, n_guesses=n_guesses, @@ -564,7 +565,7 @@ def setup_solvent(matrix, solvent_scheme): raise InputError("Cannot combine ptlr correction with iterative" " postscf linear response procedure.") - extra_matrix_terms = [] + extra_matrix_term = None energy_corrections = [] pt_schemes = [pt for pt in ["ptss", "ptlr"] if pt in solvent_scheme] @@ -584,6 +585,6 @@ def setup_solvent(matrix, solvent_scheme): f" with solvent {hf.solvent}" " not implemented.") block_fun = getattr(solvent_mat_terms, block_key) - extra_matrix_terms.append(AdcExtraTerm(matrix, {'ph_ph': block_fun})) + extra_matrix_term = AdcExtraTerm(matrix, {'ph_ph': block_fun}) - return extra_matrix_terms, energy_corrections + return extra_matrix_term, energy_corrections From a2da43d88260d489d16605efae967fb5494abe6f Mon Sep 17 00:00:00 2001 From: Maximilian Scheurer Date: Fri, 30 Apr 2021 18:01:15 +0200 Subject: [PATCH 12/21] refactoring user interface more refactoring, make adcc usable again --- adcc/AdcMatrix.py | 24 +++- adcc/AmplitudeVector.py | 4 + adcc/ReferenceState.py | 4 +- adcc/adc_pp/{solvent.py => environment.py} | 0 adcc/adc_pp/matrix.py | 2 - adcc/backends/psi4.py | 4 +- adcc/backends/pyscf.py | 4 +- adcc/backends/test_backends_polembed.py | 12 +- adcc/test_AdcMatrix.py | 8 ++ adcc/workflow.py | 147 ++++++++++---------- examples/pna/psi4_adc2_pna_6w_pol_embed.py | 2 +- examples/pna/pyscf_adc2_pna_6w_pol_embed.py | 6 +- 12 files changed, 123 insertions(+), 94 deletions(-) rename adcc/adc_pp/{solvent.py => environment.py} (100%) diff --git a/adcc/AdcMatrix.py b/adcc/AdcMatrix.py index 1a8a10ee..ab1238fb 100644 --- a/adcc/AdcMatrix.py +++ b/adcc/AdcMatrix.py @@ -86,7 +86,8 @@ class AdcMatrix(AdcMatrixlike): "adc3": dict(ph_ph=3, ph_pphh=2, pphh_ph=2, pphh_pphh=1), # noqa: E501 } - def __init__(self, method, hf_or_mp, block_orders=None, intermediates=None): + def __init__(self, method, hf_or_mp, block_orders=None, intermediates=None, + diagonal_precomputed=None): """ Initialise an ADC matrix. @@ -101,6 +102,8 @@ def __init__(self, method, hf_or_mp, block_orders=None, intermediates=None): If not set, defaults according to the selected ADC method are chosen. intermediates : adcc.Intermediates or NoneType Allows to pass intermediates to re-use to this class. + diagonal_precomputed: adcc.AmplitudeVector + Allows to pass a pre-computed diagonal, for internal use only. """ if isinstance(hf_or_mp, (libadcc.ReferenceState, libadcc.HartreeFockSolution_i)): @@ -159,9 +162,18 @@ def __init__(self, method, hf_or_mp, block_orders=None, intermediates=None): } # TODO Rename to self.block in 0.16.0 self.blocks_ph = {bl: blocks[bl].apply for bl in blocks} - self.__diagonal = sum(bl.diagonal for bl in blocks.values() - if bl.diagonal) - self.__diagonal.evaluate() + if diagonal_precomputed: + if not isinstance(diagonal_precomputed, AmplitudeVector): + raise TypeError("diagonal_precomputed needs to be" + " an AmplitudeVector.") + if diagonal_precomputed.needs_evaluation: + raise ValueError("diagonal_precomputed must already" + " be evaluated.") + self.__diagonal = diagonal_precomputed + else: + self.__diagonal = sum(bl.diagonal for bl in blocks.values() + if bl.diagonal) + self.__diagonal.evaluate() self.__init_space_data(self.__diagonal) def __iadd__(self, other): @@ -208,10 +220,10 @@ def __add__(self, other): """ if not isinstance(other, AdcExtraTerm): return NotImplemented - # NOTE: re-computes the (expensive) diagonal... ret = AdcMatrix(self.method, self.ground_state, block_orders=self.block_orders, - intermediates=self.intermediates) + intermediates=self.intermediates, + diagonal_precomputed=self.diagonal()) ret += other return ret diff --git a/adcc/AmplitudeVector.py b/adcc/AmplitudeVector.py index 0e858ef3..bc661d39 100644 --- a/adcc/AmplitudeVector.py +++ b/adcc/AmplitudeVector.py @@ -110,6 +110,10 @@ def evaluate(self): t.evaluate() return self + @property + def needs_evaluation(self): + return any(t.needs_evaluation for k, t in self.items()) + def ones_like(self): """Return an empty AmplitudeVector of the same shape and symmetry""" return AmplitudeVector(**{k: t.ones_like() for k, t in self.items()}) diff --git a/adcc/ReferenceState.py b/adcc/ReferenceState.py index 14c2b343..a3274c61 100644 --- a/adcc/ReferenceState.py +++ b/adcc/ReferenceState.py @@ -155,8 +155,8 @@ def __init__(self, hfdata, core_orbitals=None, frozen_core=None, self.orbital_coefficients, self.conv_tol ) - self.solvent = None # no solvent attached by default - for name in ["excitation_energy_corrections", "solvent"]: + self.environment = None # no environment attached by default + for name in ["excitation_energy_corrections", "environment"]: if hasattr(hfdata, name): setattr(self, name, getattr(hfdata, name)) diff --git a/adcc/adc_pp/solvent.py b/adcc/adc_pp/environment.py similarity index 100% rename from adcc/adc_pp/solvent.py rename to adcc/adc_pp/environment.py diff --git a/adcc/adc_pp/matrix.py b/adcc/adc_pp/matrix.py index 5ef0b4b1..b7cf69b3 100644 --- a/adcc/adc_pp/matrix.py +++ b/adcc/adc_pp/matrix.py @@ -59,8 +59,6 @@ def block(ground_state, spaces, order, variant=None, intermediates=None): It is assumed largely, that CVS is equivalent to mp.has_core_occupied_space, while one would probably want in the long run that one can have an "o2" space, but not do CVS. - - Returns an instance of :py:class:`AdcBlock`. """ if isinstance(variant, str): variant = [variant] diff --git a/adcc/backends/psi4.py b/adcc/backends/psi4.py index a832e7bf..70b3ec0d 100644 --- a/adcc/backends/psi4.py +++ b/adcc/backends/psi4.py @@ -109,7 +109,7 @@ def pe_energy(self, dm, elec_only=True): @property def excitation_energy_corrections(self): ret = {} - if self.solvent == "pe": + if self.environment == "pe": ptlr = EnergyCorrection( "pe_ptlr_correction", lambda view: 2.0 * self.pe_energy(view.transition_dm_ao, @@ -125,7 +125,7 @@ def excitation_energy_corrections(self): return ret @property - def solvent(self): + def environment(self): ret = None if hasattr(self.wfn, "pe_state"): ret = "pe" diff --git a/adcc/backends/pyscf.py b/adcc/backends/pyscf.py index 44ba3fc5..d6fd744a 100644 --- a/adcc/backends/pyscf.py +++ b/adcc/backends/pyscf.py @@ -128,7 +128,7 @@ def pe_energy(self, dm, elec_only=True): @property def excitation_energy_corrections(self): ret = {} - if self.solvent == "pe": + if self.environment == "pe": ptlr = EnergyCorrection( "pe_ptlr_correction", lambda view: 2.0 * self.pe_energy(view.transition_dm_ao, @@ -147,7 +147,7 @@ def excitation_energy_corrections(self): return ret @property - def solvent(self): + def environment(self): ret = None if hasattr(self.scfres, "with_solvent"): if isinstance(self.scfres.with_solvent, solvent.pol_embed.PolEmbed): diff --git a/adcc/backends/test_backends_polembed.py b/adcc/backends/test_backends_polembed.py index deae214a..2ddf7783 100644 --- a/adcc/backends/test_backends_polembed.py +++ b/adcc/backends/test_backends_polembed.py @@ -36,7 +36,7 @@ from ..testdata.static_data import pe_potentials from ..AdcMatrix import AdcExtraTerm -from ..adc_pp.solvent import block_ph_ph_0_pe +from ..adc_pp.environment import block_ph_ph_0_pe try: import cppe # noqa: F401 @@ -64,7 +64,7 @@ def template_pe_perturbative_formaldehyde(self, basis, method, backend): pe_options=pe_options) state = adcc.run_adc(scfres, method=method, n_singlets=5, conv_tol=1e-10, - solvent_scheme=["ptlr", "ptss"]) + environment=["ptlr", "ptss"]) assert_allclose( qc_result["excitation_energy"], @@ -89,7 +89,7 @@ def template_pe_perturbative_formaldehyde(self, basis, method, backend): atol=1e-5 ) - def template_pe_coupling_formaldehyde(self, basis, method, backend): + def template_pe_linear_response_formaldehyde(self, basis, method, backend): if method != "adc2": pytest.skip("Reference only exists for adc2.") basename = f"formaldehyde_{basis}_pe_{method}" @@ -112,7 +112,7 @@ def template_pe_coupling_formaldehyde(self, basis, method, backend): atol=1e-8 ) state = adcc.run_adc(matrix, n_singlets=5, conv_tol=1e-7, - solvent_scheme="hf") + environment=False) assert_allclose( state.excitation_energy_uncorrected, tm_result["excitation_energy"], @@ -122,14 +122,14 @@ def template_pe_coupling_formaldehyde(self, basis, method, backend): # invalid combination with pytest.raises(InputError): adcc.run_adc(scfres, method=method, n_singlets=5, - solvent_scheme=["postscf", "ptlr"]) + environment={"linear_response": True, "ptlr": True}) # no scheme 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, solvent_scheme="postscf") + conv_tol=1e-7, environment="linear_response") assert_allclose( state.excitation_energy_uncorrected, tm_result["excitation_energy"], diff --git a/adcc/test_AdcMatrix.py b/adcc/test_AdcMatrix.py index 523e333c..62b02441 100644 --- a/adcc/test_AdcMatrix.py +++ b/adcc/test_AdcMatrix.py @@ -222,6 +222,14 @@ def test_extra_term(self): with pytest.raises(TypeError): matrix_adc1 += 42 matrix = adcc.AdcMatrix("adc2", ground_state) + + with pytest.raises(TypeError): + adcc.AdcMatrix("adc2", ground_state, + diagonal_precomputed=42) + with pytest.raises(ValueError): + adcc.AdcMatrix("adc2", ground_state, + diagonal_precomputed=matrix.diagonal() + 42) + shift = -0.3 shifted = AdcMatrixShifted(matrix, shift) # TODO: need to do this to differentiate between diff --git a/adcc/workflow.py b/adcc/workflow.py index 8fdaafb1..91152b01 100644 --- a/adcc/workflow.py +++ b/adcc/workflow.py @@ -47,8 +47,7 @@ def run_adc(data_or_matrix, n_states=None, kind="any", conv_tol=None, n_guesses_doubles=None, output=sys.stdout, core_orbitals=None, frozen_core=None, frozen_virtual=None, method=None, n_singlets=None, n_triplets=None, n_spin_flip=None, - solvent_scheme=None, - **solverargs): + environment=None, **solverargs): """Run an ADC calculation. Main entry point to run an ADC calculation. The reference to build the ADC @@ -192,12 +191,12 @@ def run_adc(data_or_matrix, n_states=None, kind="any", conv_tol=None, if eigensolver is None: eigensolver = "davidson" - # Setup solvent coupling terms and energy corrections - ret = setup_solvent(matrix, solvent_scheme) - solvent_matrix_term, solvent_energy_corrections = ret + # Setup environment coupling terms and energy corrections + ret = setup_environment(matrix, environment) + env_matrix_term, env_energy_corrections = ret # add terms to matrix - if solvent_matrix_term: - matrix += solvent_matrix_term + if env_matrix_term: + matrix += env_matrix_term diagres = diagonalise_adcmatrix( matrix, n_states, kind, guesses=guesses, n_guesses=n_guesses, @@ -207,8 +206,8 @@ def run_adc(data_or_matrix, n_states=None, kind="any", conv_tol=None, exstates.kind = kind exstates.spin_change = spin_change - # add corrections to excited states - exstates += solvent_energy_corrections + # add environment corrections to excited states + exstates += env_energy_corrections return exstates @@ -508,83 +507,91 @@ def inner_callback(state, identifier): return inner_callback -def setup_solvent(matrix, solvent_scheme): +def setup_environment(matrix, environment): """ - Setup solvent matrix terms and/or energy corrections. + Setup environment matrix terms and/or energy corrections. Internal function called from run_adc. """ - # TODO: move someplace meaningful & reasonable... - valid_solvent_schemes = { - "hf": """ - only couple via the 'solvated' orbitals of the HF reference - state, no additional matrix terms or perturbative corrections - are used automatically - """, - "ptss": """ - perturbative state-specific (ptSS) correction, computed based on - the difference density between the ground and excited state - """, - "ptlr": """ - perturbative linear-response (ptLR) correction, computed based on - the transition density between the ground and excited state - """, - # NOTE: could also be called 'lr'... - "postscf": """ - iterative coupling to the solvent via a CIS-like coupling - density matrix, the term is added to the ADC matrix - """, - } - valid_scheme_names = list(valid_solvent_schemes.keys()) - + valid_envs = ["ptss", "ptlr", "linear_response"] + # valid_solvent_schemes = { + # "hf": """ + # only couple via the 'solvated' orbitals of the HF reference + # state, no additional matrix terms or perturbative corrections + # are used automatically + # """, + # "ptss": """ + # perturbative state-specific (ptSS) correction, computed based on + # the difference density between the ground and excited state + # """, + # "ptlr": """ + # perturbative linear-response (ptLR) correction, computed based on + # the transition density between the ground and excited state + # """, + # # NOTE: could also be called 'lr'... + # "postscf": """ + # iterative coupling to the solvent via a CIS-like coupling + # density matrix, the term is added to the ADC matrix + # """, + # } hf = matrix.reference_state - if hf.solvent and not solvent_scheme: + if hf.environment and environment is None: raise InputError( - "Solvent found in reference state, but no solvent_scheme" - " specified. Please select one or more of the following" - f" schemes: {valid_solvent_schemes.keys()}." + "Environment found in reference state, but no environment" + " configuration specified. Please select from the following" + f" schemes: {valid_envs}." ) - elif solvent_scheme and not hf.solvent: + elif environment and not hf.environment: raise InputError( - "solvent_scheme specified, but no solvent method" + "Environment specified, but no environment" " was found in reference state." ) - elif solvent_scheme is None and hf.solvent is None: - # nothing to do - return [], [] - if not isinstance(solvent_scheme, list): - solvent_scheme = [solvent_scheme] - if any(scheme not in valid_solvent_schemes for scheme in solvent_scheme): - raise InputError("Invalid solvent_scheme given." - f" Valid schemes are: {valid_scheme_names}.") - if "hf" in solvent_scheme: - if len(solvent_scheme) > 1: - raise InputError("hf solvent_scheme cannot be combined" - " with other schemes.") - return [], [] - if "postscf" in solvent_scheme and "ptlr" in solvent_scheme: - raise InputError("Cannot combine ptlr correction with iterative" - " postscf linear response procedure.") - - extra_matrix_term = None + elif not hf.environment: + environment = False + + if isinstance(environment, bool): + environment = {"ptss": True, "ptlr": True} if environment else {} + elif isinstance(environment, list): + environment = {k: True for k in environment} + elif isinstance(environment, str): + environment = {environment: True} + elif not isinstance(environment, dict): + raise TypeError("Invalid type for environment parameter" + f"' {type(environment)}'.") + + if any(env not in valid_envs for env in environment): + raise InputError("Invalid key specified for environment." + f" Valid keys are '{valid_envs}'.") + + env_matrix_term = None energy_corrections = [] - pt_schemes = [pt for pt in ["ptss", "ptlr"] if pt in solvent_scheme] - for pt in pt_schemes: + forbidden_combinations = [ + ["ptlr", "linear_response"], + ] + for fbc in forbidden_combinations: + if all(environment.get(k, False) for k in fbc): + raise InputError("Combination of environment schemes" + f" '{fbc}' not allowed. Check the" + " adcc documentation for more details.") + + for pt in ["ptss", "ptlr"]: + if not environment.get(pt, False): + continue hf_corr = hf.excitation_energy_corrections - eec_key = f"{hf.solvent}_{pt}_correction" + eec_key = f"{hf.environment}_{pt}_correction" if eec_key not in hf_corr: raise ValueError(f"{pt} correction requested, but could not find" f" the needed function {eec_key} in" f" reference state from backend {hf.backend}.") energy_corrections.append(hf_corr[eec_key]) - if "postscf" in solvent_scheme: - from adcc.adc_pp import solvent as solvent_mat_terms - block_key = f"block_ph_ph_0_{hf.solvent}" - if not hasattr(solvent_mat_terms, block_key): - raise NotImplementedError("Matrix term for postscf coupling" - f" with solvent {hf.solvent}" + if environment.get("linear_response", False): + from adcc.adc_pp import environment as adcpp_env + block_key = f"block_ph_ph_0_{hf.environment}" + if not hasattr(adcpp_env, block_key): + raise NotImplementedError("Matrix term for linear response coupling" + f" with solvent {hf.environment}" " not implemented.") - block_fun = getattr(solvent_mat_terms, block_key) - extra_matrix_term = AdcExtraTerm(matrix, {'ph_ph': block_fun}) + block_fun = getattr(adcpp_env, block_key) + env_matrix_term = AdcExtraTerm(matrix, {'ph_ph': block_fun}) - return extra_matrix_term, energy_corrections + return env_matrix_term, energy_corrections diff --git a/examples/pna/psi4_adc2_pna_6w_pol_embed.py b/examples/pna/psi4_adc2_pna_6w_pol_embed.py index 6102b8f8..e8a3b2ff 100644 --- a/examples/pna/psi4_adc2_pna_6w_pol_embed.py +++ b/examples/pna/psi4_adc2_pna_6w_pol_embed.py @@ -40,5 +40,5 @@ # Run an adc2 calculation: state = adcc.adc2(wfn, n_singlets=5, conv_tol=1e-8, - solvent_scheme=['ptss', 'ptlr']) + environment=True) print(state.describe()) diff --git a/examples/pna/pyscf_adc2_pna_6w_pol_embed.py b/examples/pna/pyscf_adc2_pna_6w_pol_embed.py index 738f4d55..cd8b7be8 100644 --- a/examples/pna/pyscf_adc2_pna_6w_pol_embed.py +++ b/examples/pna/pyscf_adc2_pna_6w_pol_embed.py @@ -38,13 +38,13 @@ # model the solvent through perturbative corrections state_pt = adcc.adc2(scfres, n_singlets=5, conv_tol=1e-5, - solvent_scheme=['ptss', 'ptlr']) + environment=['ptss', 'ptlr']) -# now model the solvent through linear-response/postscf coupling +# now model the solvent through linear-response coupling # in the ADC matrix, re-using the matrix from previous run. # This will modify state_pt.matrix state_lr = adcc.run_adc(state_pt.matrix, n_singlets=5, conv_tol=1e-5, - solvent_scheme='postscf') + environment='linear_response') print(state_pt.describe()) print(state_lr.describe()) From 983b2b0bb07bb8a2ca3b45cfee1ec4175f9ece1e Mon Sep 17 00:00:00 2001 From: Maximilian Scheurer Date: Tue, 4 May 2021 10:00:39 +0200 Subject: [PATCH 13/21] some docs and cleanup more docs and cleanup cleanup --- adcc/AdcMatrix.py | 6 +-- adcc/adc_pp/environment.py | 7 +-- adcc/backends/pyscf.py | 2 +- adcc/workflow.py | 24 ++-------- docs/calculations.rst | 91 ++++++++++++++++++++++++++++++++++++++ docs/ref.bib | 55 +++++++++++++++++++++++ 6 files changed, 156 insertions(+), 29 deletions(-) diff --git a/adcc/AdcMatrix.py b/adcc/AdcMatrix.py index ab1238fb..e381b316 100644 --- a/adcc/AdcMatrix.py +++ b/adcc/AdcMatrix.py @@ -33,10 +33,6 @@ class AdcExtraTerm: - # NOTE: currently requires the matrix - # to allow for block construction with the same - # interface as the usual adc_pp, i.e., - # reference_state, ground_state, intermediates def __init__(self, matrix, blocks): """Initialise an AdcExtraTerm. This class can be used to add customs terms @@ -109,7 +105,7 @@ def __init__(self, method, hf_or_mp, block_orders=None, intermediates=None, libadcc.HartreeFockSolution_i)): hf_or_mp = LazyMp(hf_or_mp) if not isinstance(hf_or_mp, LazyMp): - raise TypeError("mp_results is not a valid object. It needs to be " + raise TypeError("hf_or_mp is not a valid object. It needs to be " "either a LazyMp, a ReferenceState or a " "HartreeFockSolution_i.") diff --git a/adcc/adc_pp/environment.py b/adcc/adc_pp/environment.py index e8d084e1..1ee02d42 100644 --- a/adcc/adc_pp/environment.py +++ b/adcc/adc_pp/environment.py @@ -27,9 +27,10 @@ def block_ph_ph_0_pe(hf, mp, intermediates): """ Constructs an :py:class:`AdcBlock` that describes the - coupling to the polarizable environment from PE via a CIS-like - coupling density as described in 10.1021/ct300763v, eq 63. - Since the contribution depends on the input amplitude itself, + linear response coupling to the polarizable environment + from PE 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 diff --git a/adcc/backends/pyscf.py b/adcc/backends/pyscf.py index d6fd744a..dcbfdc78 100644 --- a/adcc/backends/pyscf.py +++ b/adcc/backends/pyscf.py @@ -281,7 +281,7 @@ def run_hf(xyz, basis, charge=0, multiplicity=1, conv_tol=1e-11, ) if pe_options: from pyscf.solvent import PE - mf = PE(scf.HF(mol), pe_options["potfile"]) + mf = PE(scf.HF(mol), pe_options) else: mf = scf.HF(mol) mf.conv_tol = conv_tol diff --git a/adcc/workflow.py b/adcc/workflow.py index 91152b01..8165ceab 100644 --- a/adcc/workflow.py +++ b/adcc/workflow.py @@ -129,6 +129,10 @@ def run_adc(data_or_matrix, n_states=None, kind="any", conv_tol=None, virtuals for both the MP and ADC methods performed). For ways to define these see the description in :py:class:`adcc.ReferenceState`. + environment : bool or list or dict, optional + The keywords to specify how coupling to an environment model, + e.g. PE, is treated. For details see :ref:`environment`. + Other parameters ---------------- max_subspace : int, optional @@ -513,26 +517,6 @@ def setup_environment(matrix, environment): Internal function called from run_adc. """ valid_envs = ["ptss", "ptlr", "linear_response"] - # valid_solvent_schemes = { - # "hf": """ - # only couple via the 'solvated' orbitals of the HF reference - # state, no additional matrix terms or perturbative corrections - # are used automatically - # """, - # "ptss": """ - # perturbative state-specific (ptSS) correction, computed based on - # the difference density between the ground and excited state - # """, - # "ptlr": """ - # perturbative linear-response (ptLR) correction, computed based on - # the transition density between the ground and excited state - # """, - # # NOTE: could also be called 'lr'... - # "postscf": """ - # iterative coupling to the solvent via a CIS-like coupling - # density matrix, the term is added to the ADC matrix - # """, - # } hf = matrix.reference_state if hf.environment and environment is None: raise InputError( diff --git a/docs/calculations.rst b/docs/calculations.rst index e5f7d38d..71e60995 100644 --- a/docs/calculations.rst +++ b/docs/calculations.rst @@ -27,6 +27,8 @@ to be applied. Arbitrary combinations of these variants, e.g. applying **both** CVS and FC approximations are supported as well. See :ref:`frozen-spaces` for details. +Calculations with :ref:`environment` are also supported. + General ADC(n) calculations --------------------------- General ADC(n) calculations, @@ -685,6 +687,95 @@ In fact all three may be combined jointly with any available ADC method, if desired. +.. _`environment`: + +Polarisable Embedding +--------------------- + +ADC calculations with the Polarisable Embedding (PE) model are supported +for the PySCF and Psi4 backends via the `CPPE library `_ :cite:`Scheurer2019`. +In the PE model, interactions with the environment are represented by a +multi-center multipole expansion for electrostatics, and polarisation is modeled +via dipole polarizabilities located at the expansion sites. +For a general introduction of PE and a tutorial on how to set up calculations, please see the tutorial review :cite:`Steinmann2019`. +The embedding potential needed for PE can be generated using `PyFraME `_, which is installable +via ``pip install pyframe``. + +There are different options to include environment effects in ADC excited state calculations, summarised in +the following table: + ++------------------------------------------------+-----------------------+--------------------------------------------------------------------------+-----------------------------------------------+ +| Name | ``environment`` | Comment | Reference | ++================================================+=======================+==========================================================================+===============================================+ +| coupling through reference state only | ``False`` | only couple via the 'solvated' orbitals of the SCF reference state, | :cite:`Scheurer2018` | +| | | no additional matrix terms or corrections are used | | ++------------------------------------------------+-----------------------+--------------------------------------------------------------------------+-----------------------------------------------+ +| perturbative state-specific correction (ptSS) | ``"ptss"`` | computed from the difference density betweenthe ground and excited state | :cite:`Scheurer2018` | ++------------------------------------------------+-----------------------+--------------------------------------------------------------------------+-----------------------------------------------+ +| perturbative linear-response correction (ptLR) | ``"ptlr"`` | computed from the transition density between | :cite:`Scheurer2018` | +| | | the ground and excited state | | ++------------------------------------------------+-----------------------+--------------------------------------------------------------------------+-----------------------------------------------+ +| linear response iterative coupling | ``"linear_response"`` | iterative coupling to the solvent via a CIS-like coupling density, | :cite:`Lunkenheimer2013`, :cite:`Marefat2018` | +| | | the additional term is added to the ADC matrix | | ++------------------------------------------------+-----------------------+--------------------------------------------------------------------------+-----------------------------------------------+ + +The scheme can be selected with the ``environment`` parameter in :func:`adcc.run_adc` (and also in the short-hand method functions, e.g. :func:`adcc.adc2`). +If a PE-SCF ground state is found but no ``environment`` parameter is specified, an error will be thrown. +Specifying ``environment=True`` will enable both perturbative corrections, equivalent to ``environment=["ptss", "ptlr"]``. +Combining ``"ptlr"`` with ``"linear_response"`` is not allowed since both describe the same physical effect in a different manner. + +The following example computes PE-ADC(2) excited states of para-nitroaniline in the presence of six water molecules +a) with perturbative corrections and +b) with the linear response scheme. The results of both schemes are then printed out for comparison. + +.. code-block:: python + + import adcc + from pyscf import gto, scf + from pyscf.solvent import PE + + mol = gto.M( + atom=""" + C 8.64800 1.07500 -1.71100 + C 9.48200 0.43000 -0.80800 + C 9.39600 0.75000 0.53800 + C 8.48200 1.71200 0.99500 + C 7.65300 2.34500 0.05500 + C 7.73200 2.03100 -1.29200 + H 10.18300 -0.30900 -1.16400 + H 10.04400 0.25200 1.24700 + H 6.94200 3.08900 0.38900 + H 7.09700 2.51500 -2.01800 + N 8.40100 2.02500 2.32500 + N 8.73400 0.74100 -3.12900 + O 7.98000 1.33100 -3.90100 + O 9.55600 -0.11000 -3.46600 + H 7.74900 2.71100 2.65200 + H 8.99100 1.57500 2.99500 + """, + basis='sto-3g', + ) + + scfres = PE(scf.RHF(mol), {"potfile": "pna_6w.pot"}) + scfres.conv_tol = 1e-8 + scfres.conv_tol_grad = 1e-6 + scfres.max_cycle = 250 + scfres.kernel() + + # model the solvent through perturbative corrections + state_pt = adcc.adc2(scfres, n_singlets=5, conv_tol=1e-5, + environment=['ptss', 'ptlr']) + + # now model the solvent through linear-response coupling + # in the ADC matrix, re-using the matrix from previous run. + # This will modify state_pt.matrix + state_lr = adcc.run_adc(state_pt.matrix, n_singlets=5, conv_tol=1e-5, + environment='linear_response') + + print(state_pt.describe()) + print(state_lr.describe()) + + Further examples and details ---------------------------- Some further examples can be found in the ``examples`` folder diff --git a/docs/ref.bib b/docs/ref.bib index 4d3ab548..f3d520a1 100644 --- a/docs/ref.bib +++ b/docs/ref.bib @@ -297,3 +297,58 @@ @Article{Wormit2014 pages = {774-784}, doi = {10.1080/00268976.2013.859313}, } + +@article{Scheurer2018, +author = {Scheurer, Maximilian and Herbst, Michael F. and Reinholdt, Peter and Olsen, J{\'{o}}gvan Magnus Haugaard and Dreuw, Andreas and Kongsted, Jacob}, +doi = {10.1021/acs.jctc.8b00576}, +issn = {1549-9618}, +journal = {J. Chem. Theory Comput.}, +number = {9}, +pages = {4870−4883}, +pmid = {30086234}, +title = {{Polarizable Embedding Combined with the Algebraic Diagrammatic Construction: Tackling Excited States in Biomolecular Systems}}, +volume = {14}, +year = {2018} +} + +@article{Lunkenheimer2013, +author = {Lunkenheimer, Bernd and K{\"{o}}hn, Andreas}, +doi = {10.1021/ct300763v}, +issn = {15499626}, +journal = {J. Chem. Theory Comput.}, +title = {{Solvent effects on electronically excited states using the conductor-like screening model and the second-order correlated method ADC(2)}}, +year = {2013} +} + +@article{Marefat2018, +author = {{Marefat Khah}, Alireza and {Karbalaei Khani}, Sarah and H{\"{a}}ttig, Christof}, +doi = {10.1021/acs.jctc.8b00396}, +issn = {15499626}, +journal = {J. Chem. Theory Comput.}, +pmid = {30040882}, +title = {{Analytic Excited State Gradients for the QM/MM Polarizable Embedded Second-Order Algebraic Diagrammatic Construction for the Polarization Propagator PE-ADC(2)}}, +year = {2018} +} + +@article{Scheurer2019, + title={CPPE: An open-source C++ and Python library for polarizable embedding}, + author={Scheurer, Maximilian and Reinholdt, Peter and Kjellgren, Erik Rosendahl and Olsen, Jógvan Magnus Haugaard and Dreuw, Andreas and Kongsted, Jacob}, + journal={J. Chem. Theory Comput.}, + volume={15}, + number={11}, + pages={6154--6163}, + year={2019}, + doi={10.1021/acs.jctc.9b00758}, +} + +@ARTICLE{Steinmann2019, + title = {Response properties of embedded molecules through the polarizable embedding model}, + author = {Steinmann, Casper and Reinholdt, Peter and N{\o}rby, Morten + Steen and Kongsted, Jacob and Olsen, J{\'o}gvan Magnus Haugaard}, + journal = {Int. J. Quantum Chem.}, + volume = {119}, + number = {1}, + pages = {e25717}, + year = {2019}, + doi = {10.1002/qua.25717}, +} \ No newline at end of file From b913264488a203e061bc100c0b38deea7fcd355c Mon Sep 17 00:00:00 2001 From: Jonas Leitner Date: Sun, 9 May 2021 09:55:10 +0200 Subject: [PATCH 14/21] lr_pcm. should work --- adcc/OperatorIntegrals.py | 16 +++++++++++++++- adcc/adc_pp/environment.py | 10 ++++++++++ adcc/backends/psi4.py | 10 ++++++++++ 3 files changed, 35 insertions(+), 1 deletion(-) diff --git a/adcc/OperatorIntegrals.py b/adcc/OperatorIntegrals.py index 5972938a..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,19 @@ 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"): 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 70b3ec0d..ba975ce0 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): @@ -129,6 +137,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 05c6f1b2485ddfe63c76a0ed3a84f347c79bf9ae Mon Sep 17 00:00:00 2001 From: Jonas Leitner Date: Wed, 12 May 2021 18:14:15 +0200 Subject: [PATCH 15/21] 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 ba975ce0..56de0bff 100644 --- a/adcc/backends/psi4.py +++ b/adcc/backends/psi4.py @@ -137,7 +137,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 8b09746a77f5c418aea585028feed5e8ef6956fe Mon Sep 17 00:00:00 2001 From: Maximilian Scheurer Date: Wed, 19 May 2021 16:35:49 +0200 Subject: [PATCH 16/21] reviews --- adcc/AdcMatrix.py | 15 ++++++++------- adcc/ExcitedStates.py | 6 +++--- adcc/backends/psi4.py | 6 +++--- adcc/backends/pyscf.py | 9 +++------ adcc/backends/veloxchem.py | 6 +++--- adcc/test_AdcMatrix.py | 14 +++++++------- adcc/workflow.py | 4 ++-- docs/calculations.rst | 31 +++++++++++++++++++++++++++++++ 8 files changed, 60 insertions(+), 31 deletions(-) diff --git a/adcc/AdcMatrix.py b/adcc/AdcMatrix.py index e381b316..2028252e 100644 --- a/adcc/AdcMatrix.py +++ b/adcc/AdcMatrix.py @@ -112,6 +112,14 @@ def __init__(self, method, hf_or_mp, block_orders=None, intermediates=None, if not isinstance(method, AdcMethod): method = AdcMethod(method) + if diagonal_precomputed: + if not isinstance(diagonal_precomputed, AmplitudeVector): + raise TypeError("diagonal_precomputed needs to be" + " an AmplitudeVector.") + if diagonal_precomputed.needs_evaluation: + raise ValueError("diagonal_precomputed must already" + " be evaluated.") + self.timer = Timer() self.method = method self.ground_state = hf_or_mp @@ -159,12 +167,6 @@ def __init__(self, method, hf_or_mp, block_orders=None, intermediates=None, # TODO Rename to self.block in 0.16.0 self.blocks_ph = {bl: blocks[bl].apply for bl in blocks} if diagonal_precomputed: - if not isinstance(diagonal_precomputed, AmplitudeVector): - raise TypeError("diagonal_precomputed needs to be" - " an AmplitudeVector.") - if diagonal_precomputed.needs_evaluation: - raise ValueError("diagonal_precomputed must already" - " be evaluated.") self.__diagonal = diagonal_precomputed else: self.__diagonal = sum(bl.diagonal for bl in blocks.values() @@ -194,7 +196,6 @@ def patched_apply(ampl, original=orig_app, other=other_app): self.blocks_ph[sp] = patched_apply other_diagonal = sum(bl.diagonal for bl in other.blocks.values() if bl.diagonal) - # iadd does not work with numbers self.__diagonal = self.__diagonal + other_diagonal self.__diagonal.evaluate() self.extra_terms.append(other) diff --git a/adcc/ExcitedStates.py b/adcc/ExcitedStates.py index 3ce9e641..8307c9ad 100644 --- a/adcc/ExcitedStates.py +++ b/adcc/ExcitedStates.py @@ -208,12 +208,12 @@ def __iadd__(self, other): for k in other: self += k else: - raise TypeError("Can only add EnergyCorrection (or list" - " of EnergyCorrection) to" - f" ExcitedState, not '{type(other)}'") + return NotImplemented return self def __add__(self, other): + if not isinstance(other, (EnergyCorrection, list)): + return NotImplemented ret = ExcitedStates(self, self.method, self.property_method) ret += other return ret diff --git a/adcc/backends/psi4.py b/adcc/backends/psi4.py index 70b3ec0d..3b662e49 100644 --- a/adcc/backends/psi4.py +++ b/adcc/backends/psi4.py @@ -108,7 +108,7 @@ def pe_energy(self, dm, elec_only=True): @property def excitation_energy_corrections(self): - ret = {} + ret = [] if self.environment == "pe": ptlr = EnergyCorrection( "pe_ptlr_correction", @@ -120,8 +120,8 @@ def excitation_energy_corrections(self): lambda view: self.pe_energy(view.state_diffdm_ao, elec_only=True) ) - ret["pe_ptlr_correction"] = ptlr - ret["pe_ptss_correction"] = ptss + ret.extend([ptlr, ptss]) + ret = {ec.name: ec for ec in ret} return ret @property diff --git a/adcc/backends/pyscf.py b/adcc/backends/pyscf.py index dcbfdc78..ad71279c 100644 --- a/adcc/backends/pyscf.py +++ b/adcc/backends/pyscf.py @@ -127,7 +127,7 @@ def pe_energy(self, dm, elec_only=True): @property def excitation_energy_corrections(self): - ret = {} + ret = [] if self.environment == "pe": ptlr = EnergyCorrection( "pe_ptlr_correction", @@ -139,11 +139,8 @@ def excitation_energy_corrections(self): lambda view: self.pe_energy(view.state_diffdm_ao, elec_only=True) ) - # NOTE: I don't like the duplicate 'name', - # but this way it's easier to exctract the corrections - # directly - ret["pe_ptlr_correction"] = ptlr - ret["pe_ptss_correction"] = ptss + ret.extend([ptlr, ptss]) + ret = {ec.name: ec for ec in ret} return ret @property diff --git a/adcc/backends/veloxchem.py b/adcc/backends/veloxchem.py index ef08fec6..b6dc1128 100644 --- a/adcc/backends/veloxchem.py +++ b/adcc/backends/veloxchem.py @@ -131,7 +131,7 @@ def pe_energy(self, dm, elec_only=True): @property def excitation_energy_corrections(self): - ret = {} + ret = [] if hasattr(self.scfdrv, "pe_drv"): ptlr = EnergyCorrection( "pe_ptlr_correction", @@ -143,8 +143,8 @@ def excitation_energy_corrections(self): lambda view: self.pe_energy(view.state_diffdm_ao, elec_only=True) ) - ret["pe_ptlr_correction"] = ptlr - ret["pe_ptss_correction"] = ptss + ret.extend([ptlr, ptss]) + ret = {ec.name: ec for ec in ret} return ret def get_backend(self): diff --git a/adcc/test_AdcMatrix.py b/adcc/test_AdcMatrix.py index 62b02441..ea5d0b92 100644 --- a/adcc/test_AdcMatrix.py +++ b/adcc/test_AdcMatrix.py @@ -229,20 +229,19 @@ def test_extra_term(self): with pytest.raises(ValueError): adcc.AdcMatrix("adc2", ground_state, diagonal_precomputed=matrix.diagonal() + 42) + with pytest.raises(TypeError): + AdcExtraTerm(matrix, "fail") + with pytest.raises(TypeError): + AdcExtraTerm(matrix, {"fail": "not_callable"}) shift = -0.3 shifted = AdcMatrixShifted(matrix, shift) - # TODO: need to do this to differentiate between + # TODO: need to use AmplitudeVector to differentiate between # diagonals for ph and pphh # if we just pass numbers, i.e., shift - # we get 2*shift on the diagonal :( + # we get 2*shift on the diagonal ones = matrix.diagonal().ones_like() - with pytest.raises(TypeError): - AdcExtraTerm(matrix, "fail") - with pytest.raises(TypeError): - AdcExtraTerm(matrix, {"fail": "not_callable"}) - def __shift_ph(hf, mp, intermediates): def apply(invec): return adcc.AmplitudeVector(ph=shift * invec.ph) @@ -260,6 +259,7 @@ def apply(invec): # cannot add to 'pphh_pphh' in ADC(1) matrix with pytest.raises(ValueError): matrix_adc1 += extra + shifted_2 = matrix + extra shifted_3 = extra + matrix for manual in [shifted_2, shifted_3]: diff --git a/adcc/workflow.py b/adcc/workflow.py index 8165ceab..abb72f7f 100644 --- a/adcc/workflow.py +++ b/adcc/workflow.py @@ -522,7 +522,7 @@ def setup_environment(matrix, environment): raise InputError( "Environment found in reference state, but no environment" " configuration specified. Please select from the following" - f" schemes: {valid_envs}." + f" schemes: {valid_envs} or set to False." ) elif environment and not hf.environment: raise InputError( @@ -530,7 +530,7 @@ def setup_environment(matrix, environment): " was found in reference state." ) elif not hf.environment: - environment = False + environment = {} if isinstance(environment, bool): environment = {"ptss": True, "ptlr": True} if environment else {} diff --git a/docs/calculations.rst b/docs/calculations.rst index 71e60995..92275699 100644 --- a/docs/calculations.rst +++ b/docs/calculations.rst @@ -776,6 +776,37 @@ b) with the linear response scheme. The results of both schemes are then printed print(state_lr.describe()) +The output of the last two lines is:: + + +--------------------------------------------------------------+ + | adc2 singlet , converged | + +--------------------------------------------------------------+ + | # excitation energy osc str |v1|^2 |v2|^2 | + | (au) (eV) | + | 0 0.1434972 3.904756 0.0000 0.9187 0.08128 | + | 1 0.1554448 4.229869 0.0000 0.9179 0.08211 | + | 2 0.2102638 5.721569 0.0209 0.8977 0.1023 | + | 3 0.2375643 6.464453 0.6198 0.9033 0.09666 | + | 4 0.2699134 7.344718 0.0762 0.8975 0.1025 | + +--------------------------------------------------------------+ + | excitation energy corrections included: | + | - pe_ptss_correction | + | - pe_ptlr_correction | + +--------------------------------------------------------------+ + + +--------------------------------------------------------------+ + | adc2 singlet , converged | + +--------------------------------------------------------------+ + | # excitation energy osc str |v1|^2 |v2|^2 | + | (au) (eV) | + | 0 0.1435641 3.906577 0.0000 0.9187 0.08128 | + | 1 0.1555516 4.232775 0.0000 0.9179 0.08211 | + | 2 0.210272 5.721794 0.0212 0.8977 0.1023 | + | 3 0.2378427 6.47203 0.6266 0.9034 0.09663 | + | 4 0.2698889 7.34405 0.0805 0.898 0.102 | + +--------------------------------------------------------------+ + + Further examples and details ---------------------------- Some further examples can be found in the ``examples`` folder From 34d9821dfc70189f970af4dc0b4fed6db1c4d3d2 Mon Sep 17 00:00:00 2001 From: Jonas Leitner Date: Tue, 25 May 2021 16:46:44 +0200 Subject: [PATCH 17/21] 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 75f5cfba..4cbaf271 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 0ba777086c76bf2cd5b67a7343c3427b29a47982 Mon Sep 17 00:00:00 2001 From: Jonas Leitner Date: Fri, 25 Jun 2021 09:20:46 +0200 Subject: [PATCH 18/21] rename to potential --- adcc/OperatorIntegrals.py | 12 ++++++------ adcc/adc_pp/environment.py | 3 ++- adcc/backends/psi4.py | 6 +++--- 3 files changed, 11 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 c71cb2ac..c6c7d436 100644 --- a/adcc/adc_pp/environment.py +++ b/adcc/adc_pp/environment.py @@ -42,6 +42,7 @@ 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 @@ -56,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 1a203360737a4ce4b37295c62126cac4ebd87c60 Mon Sep 17 00:00:00 2001 From: Jonas Leitner Date: Sun, 11 Jul 2021 14:45:32 +0200 Subject: [PATCH 19/21] 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 ad4c50b0b563565dd26c481c2a0d2b64fc84fb77 Mon Sep 17 00:00:00 2001 From: Jonas Leitner Date: Mon, 12 Jul 2021 11:11:54 +0200 Subject: [PATCH 20/21] 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 9e8b82c4..69702122 100644 --- a/adcc/backends/test_backends_polembed.py +++ b/adcc/backends/test_backends_polembed.py @@ -64,7 +64,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"]) @@ -99,7 +99,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 2b343e0ce722b797ca1b06e9b86bc6487cd01524 Mon Sep 17 00:00:00 2001 From: Jonas Leitner Date: Mon, 12 Jul 2021 11:24:55 +0200 Subject: [PATCH 21/21] some cleanup --- adcc/OperatorIntegrals.py | 20 -------------------- adcc/backends/test_backends_polembed.py | 3 --- 2 files changed, 23 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() diff --git a/adcc/backends/test_backends_polembed.py b/adcc/backends/test_backends_polembed.py index 69702122..9339fc70 100644 --- a/adcc/backends/test_backends_polembed.py +++ b/adcc/backends/test_backends_polembed.py @@ -25,9 +25,6 @@ import adcc import adcc.backends -from adcc.adc_pp.matrix import AdcBlock -from adcc import OneParticleOperator, AmplitudeVector - from numpy.testing import assert_allclose import pytest