From a959f8da6ecaa16f1a12b55578157696085f047e Mon Sep 17 00:00:00 2001 From: Oliver Backhouse Date: Thu, 11 Dec 2025 09:48:48 +0000 Subject: [PATCH 01/17] Fix non-dyson categorisation and remove central moments --- dyson/expressions/adc.py | 2 +- dyson/expressions/expression.py | 7 +---- dyson/expressions/gw.py | 2 +- dyson/expressions/hf.py | 45 --------------------------------- tests/conftest.py | 10 -------- tests/test_corrvec.py | 2 -- tests/test_cpgf.py | 2 -- tests/test_davidson.py | 4 --- tests/test_density.py | 2 -- tests/test_downfolded.py | 2 -- tests/test_exact.py | 2 -- tests/test_mblgf.py | 4 --- tests/test_mblse.py | 4 --- 13 files changed, 3 insertions(+), 85 deletions(-) diff --git a/dyson/expressions/adc.py b/dyson/expressions/adc.py index 8eb7c72..6aa0883 100644 --- a/dyson/expressions/adc.py +++ b/dyson/expressions/adc.py @@ -131,7 +131,7 @@ def mol(self) -> Mole: @property def non_dyson(self) -> bool: """Whether the expression produces a non-Dyson Green's function.""" - return False + return True class BaseADC_1h(BaseADC): diff --git a/dyson/expressions/expression.py b/dyson/expressions/expression.py index c866748..f85b299 100644 --- a/dyson/expressions/expression.py +++ b/dyson/expressions/expression.py @@ -464,10 +464,6 @@ def __getattr__(cls, key: str) -> type[BaseExpression]: if cls._particle is None: raise ValueError("Particle expression is not set.") return cls._particle - elif key in {"central", "dyson"}: - if cls._dyson is None: - raise ValueError("Central (Dyson) expression is not set.") - return cls._dyson elif key in {"neutral", "ee", "ph"}: if cls._neutral is None: raise ValueError("Neutral expression is not set.") @@ -481,7 +477,7 @@ def __getattr__(cls, key: str) -> type[BaseExpression]: def _classes(cls) -> set[type[BaseExpression]]: """Get all classes in the collection.""" return { - cls for cls in [cls._hole, cls._particle, cls._dyson, cls._neutral] if cls is not None + cls for cls in [cls._hole, cls._particle, cls._neutral] if cls is not None } def __contains__(cls, key: str) -> bool: @@ -498,7 +494,6 @@ class ExpressionCollection(metaclass=_ExpressionCollectionMeta): _hole: type[BaseExpression] | None = None _particle: type[BaseExpression] | None = None - _dyson: type[BaseExpression] | None = None _neutral: type[BaseExpression] | None = None _name: str | None = None diff --git a/dyson/expressions/gw.py b/dyson/expressions/gw.py index e23abf9..e61c939 100644 --- a/dyson/expressions/gw.py +++ b/dyson/expressions/gw.py @@ -145,7 +145,7 @@ def eris(self) -> Array: @property def non_dyson(self) -> bool: """Whether the expression produces a non-Dyson Green's function.""" - return False + return True class TDAGW_Dyson(BaseGW_Dyson): diff --git a/dyson/expressions/hf.py b/dyson/expressions/hf.py index f212ba3..7f8e420 100644 --- a/dyson/expressions/hf.py +++ b/dyson/expressions/hf.py @@ -188,54 +188,9 @@ def nsingle(self) -> int: return self.nvir -class HF_Dyson(BaseHF): # pylint: disable=invalid-name - """HF expressions for the Dyson Green's function.""" - - def diagonal(self) -> Array: - """Get the diagonal of the Hamiltonian. - - Returns: - Diagonal of the Hamiltonian. - """ - return self.mo_energy - - def get_excitation_vector(self, orbital: int) -> Array: - r"""Obtain the vector corresponding to a fermionic operator acting on the ground state. - - This vector is a generalisation of - - .. math:: - f_i^{\pm} \left| \Psi_0 \right> - - where :math:`f_i^{\pm}` is the fermionic creation or annihilation operator, or a product - thereof, depending on the particular expression and what Green's function it corresponds to. - - The vector defines the excitaiton manifold probed by the Green's function corresponding to - the expression. - - Args: - orbital: Orbital index. - - Returns: - Excitation vector. - """ - return util.unit_vector(self.shape[0], orbital) - - @property - def nsingle(self) -> int: - """Number of configurations in the singles sector.""" - return self.nocc + self.nvir - - @property - def non_dyson(self) -> bool: - """Whether the expression produces a non-Dyson Green's function.""" - return False - - class HF(ExpressionCollection): """Collection of HF expressions for different parts of the Green's function.""" _hole = HF_1h _particle = HF_1p - _dyson = HF_Dyson _name = "HF" diff --git a/tests/conftest.py b/tests/conftest.py index 5fbd1ad..832a8ea 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -166,16 +166,6 @@ def _get_central_result( allow_hermitian: bool = True, ) -> Spectral: """Get the central result for the given mean-field method.""" - if "dyson" in expression_method: - expression = expression_method.dyson.from_mf(mf) - if expression.nconfig > 1024: - pytest.skip("Skipping test for large Hamiltonian") - if not expression.hermitian and not allow_hermitian: - pytest.skip("Skipping test for non-Hermitian Hamiltonian with negative weights") - exact = exact_cache(mf, expression_method.dyson) - assert exact.result is not None - return exact.result - # Combine hole and particle results expression_h = expression_method.h.from_mf(mf) expression_p = expression_method.p.from_mf(mf) diff --git a/tests/test_corrvec.py b/tests/test_corrvec.py index a39171a..e827315 100644 --- a/tests/test_corrvec.py +++ b/tests/test_corrvec.py @@ -28,8 +28,6 @@ def test_vs_exact_solver( expression = expression_cls.from_mf(mf) if expression.nconfig > 1024: # TODO: Make larger for CI runs? pytest.skip("Skipping test for large Hamiltonian") - if expression.nsingle == (expression.nocc + expression.nvir): - pytest.skip("Skipping test for central Hamiltonian") grid = RealFrequencyGrid.from_uniform(-2, 2, 16, 0.1) # Solve the Hamiltonian exactly diff --git a/tests/test_cpgf.py b/tests/test_cpgf.py index 0f0ba11..d42674c 100644 --- a/tests/test_cpgf.py +++ b/tests/test_cpgf.py @@ -29,8 +29,6 @@ def test_vs_exact_solver( expression = expression_cls.from_mf(mf) if expression.nconfig > 1024: # TODO: Make larger for CI runs? pytest.skip("Skipping test for large Hamiltonian") - if expression.nsingle == (expression.nocc + expression.nvir): - pytest.skip("Skipping test for central Hamiltonian") if isinstance(expression, BaseHF): pytest.skip("Skipping test for HF Hamiltonian") grid = RealFrequencyGrid.from_uniform(-2, 2, 16, 0.1) diff --git a/tests/test_davidson.py b/tests/test_davidson.py index 1dc347c..375d23f 100644 --- a/tests/test_davidson.py +++ b/tests/test_davidson.py @@ -30,8 +30,6 @@ def test_vs_exact_solver( expression = expression_cls.from_mf(mf) if expression.nconfig > 1024: # TODO: Make larger for CI runs? pytest.skip("Skipping test for large Hamiltonian") - if expression.nsingle == (expression.nocc + expression.nvir): - pytest.skip("Skipping test for central Hamiltonian") bra: Array = np.array(expression.get_excitation_bras()) ket: Array = np.array(expression.get_excitation_kets()) @@ -81,8 +79,6 @@ def test_vs_exact_solver_central( ) -> None: """Test the exact solver for central moments.""" # Get the quantities required from the expressions - if "o" not in expression_method or "p" not in expression_method: - pytest.skip("Skipping test for Dyson only expression") expression_h = expression_method.h.from_mf(mf) expression_p = expression_method.p.from_mf(mf) if expression_h.nconfig > 1024 or expression_p.nconfig > 1024: diff --git a/tests/test_density.py b/tests/test_density.py index b87682d..dafb6e7 100644 --- a/tests/test_density.py +++ b/tests/test_density.py @@ -26,8 +26,6 @@ def test_vs_exact_solver( exact_cache: ExactGetter, ) -> None: """Test DensityRelaxation compared to the exact solver.""" - if "h" not in expression_method or "p" not in expression_method: - pytest.skip("Skipping test for Dyson only expression") expression_h = expression_method.h.from_mf(mf) expression_p = expression_method.p.from_mf(mf) if expression_h.nconfig > 1024 or expression_p.nconfig > 1024: diff --git a/tests/test_downfolded.py b/tests/test_downfolded.py index 32be3d7..57c4de5 100644 --- a/tests/test_downfolded.py +++ b/tests/test_downfolded.py @@ -26,8 +26,6 @@ def test_vs_exact_solver( ) -> None: """Test Downfolded compared to the exact solver.""" # Get the quantities required from the expressions - if "h" not in expression_method or "p" not in expression_method: - pytest.skip("Skipping test for Dyson only expression") expression_h = expression_method.h.from_mf(mf) expression_p = expression_method.p.from_mf(mf) if expression_h.nconfig > 1024 or expression_p.nconfig > 1024: diff --git a/tests/test_exact.py b/tests/test_exact.py index e3d4993..dd514d6 100644 --- a/tests/test_exact.py +++ b/tests/test_exact.py @@ -69,8 +69,6 @@ def test_vs_exact_solver_central( ) -> None: """Test the exact solver for central moments.""" # Get the quantities required from the expressions - if "h" not in expression_method or "p" not in expression_method: - pytest.skip("Skipping test for Dyson only expression") expression_h = expression_method.h.from_mf(mf) expression_p = expression_method.p.from_mf(mf) if expression_h.nconfig > 1024 or expression_p.nconfig > 1024: diff --git a/tests/test_mblgf.py b/tests/test_mblgf.py index ad96789..611b906 100644 --- a/tests/test_mblgf.py +++ b/tests/test_mblgf.py @@ -27,8 +27,6 @@ def test_central_moments( ) -> None: """Test the recovery of the exact central moments from the MBLGF solver.""" # Get the quantities required from the expression - if "h" not in expression_method or "p" not in expression_method: - pytest.skip("Skipping test for Dyson only expression") expression_h = expression_method.h.from_mf(mf) expression_p = expression_method.p.from_mf(mf) nmom_gf = max_cycle * 2 + 2 @@ -68,8 +66,6 @@ def test_vs_exact_solver_central( ) -> None: """Test the MBLGF solver for central moments.""" # Get the quantities required from the expressions - if "h" not in expression_method or "p" not in expression_method: - pytest.skip("Skipping test for Dyson only expression") expression_h = expression_method.h.from_mf(mf) expression_p = expression_method.p.from_mf(mf) if expression_h.nconfig > 1024 or expression_p.nconfig > 1024: diff --git a/tests/test_mblse.py b/tests/test_mblse.py index 30f4ca8..f0f7468 100644 --- a/tests/test_mblse.py +++ b/tests/test_mblse.py @@ -27,8 +27,6 @@ def test_central_moments( ) -> None: """Test the recovery of the exact central moments from the MBLSE solver.""" # Get the quantities required from the expression - if "h" not in expression_method or "p" not in expression_method: - pytest.skip("Skipping test for Dyson only expression") expression_h = expression_method.h.from_mf(mf) expression_p = expression_method.p.from_mf(mf) nmom_gf = max_cycle * 2 + 4 @@ -61,8 +59,6 @@ def test_vs_exact_solver_central( max_cycle: int, ) -> None: # Get the quantities required from the expressions - if "h" not in expression_method or "p" not in expression_method: - pytest.skip("Skipping test for Dyson only expression") expression_h = expression_method.h.from_mf(mf) expression_p = expression_method.p.from_mf(mf) if expression_h.nconfig > 1024 or expression_p.nconfig > 1024: From cc5e05cecd073a943984f628c5d83d5d220d72cb Mon Sep 17 00:00:00 2001 From: Oliver Backhouse Date: Thu, 11 Dec 2025 09:51:07 +0000 Subject: [PATCH 02/17] Remove GW for now --- dyson/__init__.py | 5 +- dyson/expressions/__init__.py | 2 - dyson/expressions/gw.py | 252 ---------------------------------- tests/conftest.py | 6 +- tests/test_expressions.py | 26 +--- 5 files changed, 5 insertions(+), 286 deletions(-) delete mode 100644 dyson/expressions/gw.py diff --git a/dyson/__init__.py b/dyson/__init__.py index 6f3ab4c..e6b7aa5 100644 --- a/dyson/__init__.py +++ b/dyson/__init__.py @@ -98,9 +98,6 @@ * - :data:`~dyson.expressions.adc.ADC2x` - Algebraic diagrammatic construction extended second order excited states, based on a mean-field ground state. - * - :data:`~dyson.expressions.gw.TDAGW` - - GW theory with the Tamm--Dancoff approximation for the excited states, based on a - mean-field ground state. * - :data:`~dyson.expressions.hamiltonian.Hamiltonian` - General Hamiltonian expression, which accepts an array representing the supermatrix of the self-energy, and supports :mod:`scipy.sparse` matrices. @@ -139,4 +136,4 @@ CorrectionVector, CPGF, ) -from dyson.expressions import HF, CCSD, FCI, ADC2, ADC2x, TDAGW, Hamiltonian +from dyson.expressions import HF, CCSD, FCI, ADC2, ADC2x, Hamiltonian diff --git a/dyson/expressions/__init__.py b/dyson/expressions/__init__.py index 24d7b00..c6c5c43 100644 --- a/dyson/expressions/__init__.py +++ b/dyson/expressions/__init__.py @@ -81,7 +81,6 @@ ccsd fci adc - gw hamiltonian """ @@ -90,5 +89,4 @@ from dyson.expressions.ccsd import CCSD from dyson.expressions.fci import FCI from dyson.expressions.adc import ADC2, ADC2x -from dyson.expressions.gw import TDAGW from dyson.expressions.hamiltonian import Hamiltonian diff --git a/dyson/expressions/gw.py b/dyson/expressions/gw.py deleted file mode 100644 index e61c939..0000000 --- a/dyson/expressions/gw.py +++ /dev/null @@ -1,252 +0,0 @@ -"""GW approximation expressions [1]_ [2]_ [3]_. - -.. [1] Hedin, L. (1965). New Method for Calculating the One-Particle Green’s Function with - Application to the Electron-Gas Problem. Physical Review, 139(3A), A796–A823. - https://doi.org/10.1103/physrev.139.a796 - -.. [2] Aryasetiawan, F., & Gunnarsson, O. (1998). The GW method. Reports on Progress - in Physics, 61(3), 237–312. https://doi.org/10.1088/0034-4885/61/3/002 - -.. [3] Zhu, T., & Chan, G. K. (2021). All-Electron Gaussian-Based G0W0 for valence and core - excitation energies of periodic systems. Journal of Chemical Theory and Computation, 17(2), - 727–741. https://doi.org/10.1021/acs.jctc.0c00704 - -""" - -from __future__ import annotations - -from typing import TYPE_CHECKING - -from pyscf import gw, lib - -from dyson import numpy as np -from dyson import util -from dyson.expressions.expression import BaseExpression, ExpressionCollection -from dyson.representations.enums import Reduction - -if TYPE_CHECKING: - from pyscf.gto.mole import Mole - from pyscf.scf.hf import RHF - - from dyson.typing import Array - - -class BaseGW_Dyson(BaseExpression): - """Base class for GW expressions for the Dyson Green's function.""" - - hermitian_downfolded = True - hermitian_upfolded = False - - def __init__( - self, - mol: Mole, - gw_obj: gw.GW, - eris: Array | None = None, - ) -> None: - """Initialise the expression. - - Args: - mol: Molecule object. - gw_obj: GW object from PySCF. - eris: Density fitted electron repulsion integrals from PySCF. - """ - self._mol = mol - self._gw = gw_obj - self._eris = eris if eris is not None else gw_obj.ao2mo() - - if getattr(self._gw._scf, "xc", "hf") != "hf": - raise NotImplementedError( - "GW expressions currently only support Hartree--Fock mean-field objects." - ) - - @classmethod - def from_mf(cls, mf: RHF) -> BaseGW_Dyson: - """Create an expression from a mean-field object. - - Args: - mf: Mean-field object. - - Returns: - Expression object. - """ - return cls.from_gw(gw.GW(mf)) - - @classmethod - def from_gw(cls, gw: gw.GW) -> BaseGW_Dyson: - """Create an expression from a GW object. - - Args: - gw: GW object. - - Returns: - Expression object. - """ - return cls(gw._scf.mol, gw) - - def build_se_moments(self, nmom: int, reduction: Reduction = Reduction.NONE) -> Array: - """Build the self-energy moments. - - Args: - nmom: Number of moments to compute. - reduction: Reduction method to apply to the moments. - - Returns: - Moments of the self-energy. - """ - raise NotImplementedError("Self-energy moments not implemented for GW.") - - def get_excitation_vector(self, orbital: int) -> Array: - r"""Obtain the vector corresponding to a fermionic operator acting on the ground state. - - This vector is a generalisation of - - .. math:: - f_i^{\pm} \left| \Psi_0 \right> - - where :math:`f_i^{\pm}` is the fermionic creation or annihilation operator, or a product - thereof, depending on the particular expression and what Green's function it corresponds to. - - The vector defines the excitaiton manifold probed by the Green's function corresponding to - the expression. - - Args: - orbital: Orbital index. - - Returns: - Excitation vector. - """ - return util.unit_vector(self.shape[0], orbital) - - @property - def nsingle(self) -> int: - """Number of configurations in the singles sector.""" - return self.nocc + self.nvir - - @property - def nconfig(self) -> int: - """Number of configurations.""" - return self.nocc * self.nocc * self.nvir + self.nvir * self.nvir * self.nocc - - @property - def mol(self) -> Mole: - """Molecule object.""" - return self._mol - - @property - def gw(self) -> gw.GW: - """GW object.""" - return self._gw - - @property - def eris(self) -> Array: - """Density fitted electron repulsion integrals.""" - return self._eris - - @property - def non_dyson(self) -> bool: - """Whether the expression produces a non-Dyson Green's function.""" - return True - - -class TDAGW_Dyson(BaseGW_Dyson): - """GW expressions with Tamm--Dancoff (TDA) approximation for the Dyson Green's function.""" - - def apply_hamiltonian(self, vector: Array) -> Array: - """Apply the Hamiltonian to a vector. - - Args: - vector: Vector to apply Hamiltonian to. - - Returns: - Output vector. - """ - # Get the slices for each sector - o1 = slice(None, self.nocc) - v1 = slice(self.nocc, self.nocc + self.nvir) - o2 = slice(self.nocc + self.nvir, self.nocc + self.nvir + self.nocc * self.nocc * self.nvir) - v2 = slice(self.nocc + self.nvir + self.nocc * self.nocc * self.nvir, None) - - # Get the blocks of the ERIs - Lia = self.eris[:, o1, v1] - Lai = self.eris[:, v1, o1] - Lij = self.eris[:, o1, o1] - Lab = self.eris[:, v1, v1] - - # Get the blocks of the vector - vector_o1 = vector[o1] - vector_v1 = vector[v1] - vector_o2 = vector[o2].reshape(self.nocc, self.nocc, self.nvir) - vector_v2 = vector[v2].reshape(self.nocc, self.nvir, self.nvir) - - # Get the energy denominators - mo_energy = self.gw._scf.mo_energy if self.gw.mo_energy is None else self.gw.mo_energy - e_ija = lib.direct_sum("i+j-a->ija", mo_energy[o1], mo_energy[o1], mo_energy[v1]) - e_iab = lib.direct_sum("a+b-i->iab", mo_energy[v1], mo_energy[v1], mo_energy[o1]) - - # Perform the contractions - r_o1 = mo_energy[o1] * vector_o1 - r_o1 += util.einsum("Qik,Qcl,klc->i", Lij, Lai, vector_o2) * 2 - r_o1 += util.einsum("Qid,Qkc,kcd->i", Lia.conj(), Lia.conj(), vector_v2) * 2 - - r_v1 = mo_energy[v1] * vector_v1 - r_v1 += util.einsum("Qak,Qcl,klc->a", Lai, Lai, vector_o2) * 2 - r_v1 += util.einsum("Qad,Qkc,kcd->a", Lab.conj(), Lia.conj(), vector_v2) * 2 - - r_o2 = util.einsum("Qki,Qaj,k->ija", Lij.conj(), Lai.conj(), vector_o1) - r_o2 += util.einsum("Qbi,Qaj,b->ija", Lai.conj(), Lai.conj(), vector_v1) - r_o2 += util.einsum("ija,ija->ija", e_ija, vector_o2) - r_o2 -= util.einsum("Qja,Qlc,ilc->ija", Lia, Lia, vector_o2) * 2 - - r_v2 = util.einsum("Qjb,Qia,j->iab", Lia, Lia, vector_o1) - r_v2 += util.einsum("Qcb,Qia,c->iab", Lab, Lia, vector_v1) - r_v2 += util.einsum("iab,iab->iab", e_iab, vector_v2) - r_v2 += util.einsum("Qia,Qkc,kcb->iab", Lia, Lia, vector_v2) * 2 - - return np.concatenate([r_o1, r_v1, r_o2.ravel(), r_v2.ravel()]) - - def apply_hamiltonian_left(self, vector: Array) -> Array: - """Apply the Hamiltonian to a vector on the left. - - Args: - vector: Vector to apply Hamiltonian to. - - Returns: - Output vector. - """ - raise NotImplementedError("Left application of Hamiltonian is not implemented for TDA-GW.") - - def diagonal(self) -> Array: - """Get the diagonal of the Hamiltonian. - - Returns: - Diagonal of the Hamiltonian. - """ - # Get the slices for each sector - o1 = slice(None, self.nocc) - v1 = slice(self.nocc, None) - - # Get the blocks of the ERIs - Lia = self.eris[:, o1, v1] - Lai = self.eris[:, v1, o1] - - # Get the energy denominators - mo_energy = self.gw._scf.mo_energy if self.gw.mo_energy is None else self.gw.mo_energy - e_ija = lib.direct_sum("i+j-a->ija", mo_energy[o1], mo_energy[o1], mo_energy[v1]) - e_iab = lib.direct_sum("a+b-i->iab", mo_energy[v1], mo_energy[v1], mo_energy[o1]) - - # Build the diagonal - diag_o1 = mo_energy[o1].copy() - diag_v1 = mo_energy[v1].copy() - diag_o2 = e_ija.ravel() - diag_o2 -= util.einsum("Qja,Qaj,ii->ija", Lia, Lai, np.eye(self.nocc)).ravel() - diag_v2 = e_iab.ravel() - diag_v2 += util.einsum("Qai,Qia,bb->iab", Lai, Lia, np.eye(self.nvir)).ravel() - - return np.concatenate([diag_o1, diag_v1, diag_o2, diag_v2]) - - -class TDAGW(ExpressionCollection): - """Collection of TDAGW expressions for different parts of the Green's function.""" - - _dyson = TDAGW_Dyson - _name = "TDA-GW" diff --git a/tests/conftest.py b/tests/conftest.py index 832a8ea..147558c 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -8,7 +8,7 @@ from pyscf import gto, scf from dyson import numpy as np -from dyson.expressions import ADC2, CCSD, FCI, HF, TDAGW, ADC2x +from dyson.expressions import ADC2, CCSD, FCI, HF, ADC2x from dyson.representations.lehmann import Lehmann from dyson.representations.spectral import Spectral from dyson.solvers import Exact @@ -51,8 +51,8 @@ dm = mf.make_rdm1(mo, mf.mo_occ) mf = mf.run(dm) -METHODS = [HF, CCSD, FCI, ADC2, ADC2x, TDAGW] -METHOD_NAMES = ["HF", "CCSD", "FCI", "ADC2", "ADC2x", "TDAGW"] +METHODS = [HF, CCSD, FCI, ADC2, ADC2x] +METHOD_NAMES = ["HF", "CCSD", "FCI", "ADC2", "ADC2x"] def pytest_generate_tests(metafunc): # type: ignore diff --git a/tests/test_expressions.py b/tests/test_expressions.py index ec521fd..f5d24ef 100644 --- a/tests/test_expressions.py +++ b/tests/test_expressions.py @@ -10,7 +10,7 @@ import pytest from dyson import util -from dyson.expressions import ADC2, CCSD, FCI, HF, TDAGW, ADC2x +from dyson.expressions import ADC2, CCSD, FCI, HF, ADC2x from dyson.solvers import Davidson, Exact if TYPE_CHECKING: @@ -283,27 +283,3 @@ def test_adc2x(mf: scf.hf.RHF) -> None: gf_moment_0 = gf_moments_h[0] + gf_moments_p[0] assert np.allclose(gf_moment_0, np.eye(mf.mol.nao)) - - -def test_tdagw(mf: scf.hf.RHF, exact_cache: ExactGetter) -> None: - """Test the TDAGW expression.""" - tdagw = TDAGW["dyson"].from_mf(mf) - dft = mf.to_rks() - dft.xc = "hf" - - td = pyscf.tdscf.dTDA(dft) - td.nstates = np.sum(mf.mo_occ > 0) * np.sum(mf.mo_occ == 0) - td.kernel() - td.xy = np.array([(x, np.zeros_like(x)) for x, y in td.xy]) - gw_obj = pyscf.gw.GW(dft, tdmf=td, freq_int="exact") - gw_obj.kernel() - - # Get the IPs and EAs from the Exact solver - solver = exact_cache(mf, TDAGW["dyson"]) - assert solver.result is not None - gf = solver.result.get_greens_function() - mo_energy = gf.as_perturbed_mo_energy() - - # No diagonal approximation in TDAGW so large error - assert np.abs(mo_energy[tdagw.nocc - 1] - gw_obj.mo_energy[tdagw.nocc - 1]) < 1e-3 - assert np.abs(mo_energy[tdagw.nocc] - gw_obj.mo_energy[tdagw.nocc]) < 1e-3 From 5896053a6da5833048809e26329b888bb2ae4471 Mon Sep 17 00:00:00 2001 From: Oliver Backhouse Date: Wed, 17 Dec 2025 18:24:05 +0000 Subject: [PATCH 03/17] Refactor Spectral.combine --- dyson/representations/spectral.py | 55 +++--- dyson/solvers/static/davidson.py | 6 +- dyson/solvers/static/exact.py | 66 +------ dyson/util/__init__.py | 1 + dyson/util/linalg.py | 68 ++++++++ tests/test_expressions.py | 280 +++++++++++++++++++++++++----- 6 files changed, 339 insertions(+), 137 deletions(-) diff --git a/dyson/representations/spectral.py b/dyson/representations/spectral.py index 815cba7..881561c 100644 --- a/dyson/representations/spectral.py +++ b/dyson/representations/spectral.py @@ -249,23 +249,13 @@ def combine(self, *args: Spectral, chempot: float | None = None) -> Spectral: Returns: Combined spectral representation. """ - # TODO: just concatenate the eigenvectors...? args = (self, *args) if len(set(arg.nphys for arg in args)) != 1: raise ValueError( "All Spectral objects must have the same number of physical degrees of freedom." ) nphys = args[0].nphys - - # Sum the overlap and static self-energy matrices - static = sum([arg.get_static_self_energy() for arg in args], np.zeros((nphys, nphys))) - overlap = sum([arg.get_overlap() for arg in args], np.zeros((nphys, nphys))) - - # Orthogonalise under the metric of the overlap matrix hermitian = all(arg.hermitian for arg in args) - orth, _ = util.matrix_power(overlap, -0.5, hermitian=hermitian, return_error=False) - static = orth @ static @ orth - overlap = orth @ overlap @ orth # Check the chemical potentials if chempot is None: @@ -273,32 +263,43 @@ def combine(self, *args: Spectral, chempot: float | None = None) -> Spectral: chempots = [arg.chempot for arg in args if arg.chempot is not None] if not all(np.isclose(chempots[0], part) for part in chempots[1:]): raise ValueError( - "If not chempot is passed to combine, all chemical potentials must be " + "If chempot is not passed to combine, all chemical potentials must be " "equal in the inputs." ) chempot = chempots[0] - # Get the auxiliaries - energies = np.zeros((0)) + # Get the eigenvalues and dyson orbitals + eigvals = np.concatenate([arg.eigvals for arg in args]) left = np.zeros((nphys, 0)) right = np.zeros((nphys, 0)) for arg in args: - energies_i, couplings_i = arg.get_auxiliaries() - energies = np.concatenate([energies, energies_i]) - if arg.hermitian: - left = np.concatenate([left, couplings_i], axis=1) - else: - left_i, right_i = util.unpack_vectors(couplings_i) - left = np.concatenate([left, left_i], axis=1) - right = np.concatenate([right, right_i], axis=1) - couplings = np.array([left, right]) if not args[0].hermitian else left - - # Solve the eigenvalue problem - self_energy = Lehmann(energies, couplings) - result = Spectral( - *self_energy.diagonalise_matrix(static, overlap=overlap), nphys, chempot=chempot + orbitals = arg.get_dyson_orbitals()[1] + left_i, right_i = util.unpack_vectors(orbitals) + left = np.concatenate([left, left_i], axis=1) + right = np.concatenate([right, right_i], axis=1) + + # Orthogonalise under the metric of the overlap matrix + overlap = sum([arg.get_overlap() for arg in args], np.zeros((nphys, nphys))) + orth, _ = util.matrix_power(overlap, -0.5, hermitian=hermitian, return_error=False) + left = util.rotate_subspace(left, orth) + right = util.rotate_subspace(right, orth.T.conj()) + + # Construct the eigenvectors + eigvecs = util.project_eigenvectors( + None, + left, + right if not hermitian else None, ) + # Unorthogonalise the eigenvectors + unorth, _ = util.matrix_power(overlap, 0.5, hermitian=hermitian, return_error=False) + left = util.rotate_subspace(left, unorth.T.conj()) + right = util.rotate_subspace(right, unorth) + eigvecs = np.array([left, right]) + + # Construct the result + result = Spectral(eigvals, eigvecs, nphys, chempot=chempot, sort=True) + return result @cached_property diff --git a/dyson/solvers/static/davidson.py b/dyson/solvers/static/davidson.py index 929a746..5ae1343 100644 --- a/dyson/solvers/static/davidson.py +++ b/dyson/solvers/static/davidson.py @@ -21,7 +21,7 @@ from dyson.representations.lehmann import Lehmann from dyson.representations.spectral import Spectral from dyson.solvers.solver import StaticSolver -from dyson.solvers.static.exact import orthogonalise_self_energy, project_eigenvectors +from dyson.solvers.static.exact import orthogonalise_self_energy if TYPE_CHECKING: from typing import Any, Callable @@ -299,7 +299,9 @@ def _callback(env: dict[str, Any]) -> None: converged = converged[mask] # Get the full map onto physical + auxiliary and rotate the eigenvectors - eigvecs = project_eigenvectors(eigvecs, self.bra, self.ket if not self.hermitian else None) + eigvecs = util.project_eigenvectors( + eigvecs, self.bra, self.ket if not self.hermitian else None + ) # Store the results self.result = Spectral(eigvals, eigvecs, self.nphys) diff --git a/dyson/solvers/static/exact.py b/dyson/solvers/static/exact.py index 4111644..3a2c091 100644 --- a/dyson/solvers/static/exact.py +++ b/dyson/solvers/static/exact.py @@ -17,68 +17,6 @@ from dyson.typing import Array -def project_eigenvectors( - eigvecs: Array, - bra: Array, - ket: Array | None = None, -) -> Array: - """Project eigenvectors onto the physical plus auxiliary space. - - Args: - eigvecs: Eigenvectors to be projected. - bra: Bra state vector mapping the supermatrix to the physical space. - ket: Ket state vector mapping the supermatrix to the physical space. If ``None``, use the - same vectors as ``bra``. - - Returns: - Projected eigenvectors. - - Notes: - The physical space is defined by the ``bra`` and ``ket`` vectors, while the auxiliary part - is defined by the null space of the projector formed by the outer product of these vectors. - """ - hermitian = ket is None - nphys = bra.shape[0] - if not hermitian and eigvecs.ndim == 2: - raise ValueError( - "bra and ket both passed implying a non-hermitian system, but eigvecs is 2D." - ) - if ket is None: - ket = bra - - # Find a basis for the null space of the bra and ket vectors - projector = ket.T @ bra.conj() - vectors = util.null_space_basis(projector, hermitian=hermitian) - - # If the system is hermitian, the rotation is trivial - if hermitian: - rotation = np.concatenate([bra.T, vectors[0]], axis=1) - return rotation.T.conj() @ eigvecs - - # If the system is not hermitian, we need to ensure biorthonormality - overlap = ket.conj() @ bra.T - orth, orth_error = util.matrix_power(overlap, -0.5, hermitian=hermitian, return_error=True) - unorth, unorth_error = util.matrix_power(overlap, 0.5, hermitian=hermitian, return_error=True) - - # Work in an orthonormal physical basis - bra = bra.T @ orth - ket = ket.T @ orth.T.conj() - - # Biorthonormalise the physical plus auxiliary vectors - left = np.concatenate([ket, vectors[0]], axis=1) - right = np.concatenate([bra, vectors[1]], axis=1) - left, right = util.biorthonormalise(left, right) - - # Return the physical vectors to the original basis - left[:, :nphys] = left[:, :nphys] @ unorth.T.conj() - right[:, :nphys] = right[:, :nphys] @ unorth - - # Rotate the eigenvectors - eigvecs = np.array([left.T.conj() @ eigvecs[0], right.T.conj() @ eigvecs[1]]) - - return eigvecs - - def orthogonalise_self_energy( static: Array, self_energy: Lehmann, @@ -251,7 +189,9 @@ def kernel(self) -> Spectral: eigvecs = np.array([left, right]) # Get the full map onto physical + auxiliary and rotate the eigenvectors - eigvecs = project_eigenvectors(eigvecs, self.bra, self.ket if not self.hermitian else None) + eigvecs = util.project_eigenvectors( + eigvecs, self.bra, self.ket if not self.hermitian else None + ) # Store the result self.result = Spectral(eigvals, eigvecs, self.nphys) diff --git a/dyson/util/__init__.py b/dyson/util/__init__.py index d968eef..13e96ac 100644 --- a/dyson/util/__init__.py +++ b/dyson/util/__init__.py @@ -28,6 +28,7 @@ as_trace, unit_vector, null_space_basis, + project_eigenvectors, concatenate_paired_vectors, unpack_vectors, block_diag, diff --git a/dyson/util/linalg.py b/dyson/util/linalg.py index 7f8d9c6..f9a214f 100644 --- a/dyson/util/linalg.py +++ b/dyson/util/linalg.py @@ -299,6 +299,74 @@ def null_space_basis( return (left, right) if hermitian else (left, left) +@cache_by_id +def project_eigenvectors( + eigvecs: Array | None, + bra: Array, + ket: Array | None = None, +) -> Array: + """Project eigenvectors onto the physical plus auxiliary space. + + Args: + eigvecs: Eigenvectors to be projected. If ``None``, assume identity. + bra: Bra state vector mapping the supermatrix to the physical space. + ket: Ket state vector mapping the supermatrix to the physical space. If ``None``, use the + same vectors as ``bra``. + + Returns: + Projected eigenvectors. + + Notes: + The physical space is defined by the ``bra`` and ``ket`` vectors, while the auxiliary part + is defined by the null space of the projector formed by the outer product of these vectors. + """ + hermitian = ket is None + nphys = bra.shape[0] + if not hermitian and eigvecs is not None and eigvecs.ndim == 2: + raise ValueError( + "bra and ket both passed implying a non-hermitian system, but eigvecs is 2D." + ) + if ket is None: + ket = bra + + # Find a basis for the null space of the bra and ket vectors + projector = ket.T @ bra.conj() + vectors = null_space_basis(projector, hermitian=hermitian) + + # If the system is hermitian, the rotation is trivial + if hermitian: + rotation = np.concatenate([bra.T, vectors[0]], axis=1) + return rotation.T.conj() @ eigvecs if eigvecs is not None else rotation.T.conj() + + # If the system is not hermitian, we need to ensure biorthonormality + overlap = ket.conj() @ bra.T + orth, orth_error = matrix_power(overlap, -0.5, hermitian=hermitian, return_error=True) + unorth, unorth_error = matrix_power(overlap, 0.5, hermitian=hermitian, return_error=True) + + # Work in an orthonormal physical basis + bra = bra.T @ orth + ket = ket.T @ orth.T.conj() + + # Biorthonormalise the physical plus auxiliary vectors + left = np.concatenate([ket, vectors[0]], axis=1) + right = np.concatenate([bra, vectors[1]], axis=1) + left, right = biorthonormalise(left, right) + + # Return the physical vectors to the original basis + left[:, :nphys] = left[:, :nphys] @ unorth.T.conj() + right[:, :nphys] = right[:, :nphys] @ unorth + + # Rotate the eigenvectors + eigvecs = np.array( + [ + left.T.conj() @ eigvecs[0] if eigvecs is not None else left.T.conj(), + right.T.conj() @ eigvecs[1] if eigvecs is not None else right.T.conj(), + ] + ) + + return eigvecs + + @cache_by_id def matrix_power( matrix: Array, diff --git a/tests/test_expressions.py b/tests/test_expressions.py index f5d24ef..e93ff77 100644 --- a/tests/test_expressions.py +++ b/tests/test_expressions.py @@ -14,9 +14,12 @@ from dyson.solvers import Davidson, Exact if TYPE_CHECKING: + from typing import Any + from pyscf import scf from dyson.expressions.expression import BaseExpression + from dyson.solvers.solver import BaseSolver from .conftest import ExactGetter, Helper @@ -102,41 +105,89 @@ def test_static( assert np.allclose(static, gf_moments[1]) -def test_hf(mf: scf.hf.RHF) -> None: +@pytest.mark.parametrize( + "solver_cls, solver_kwargs", + [ + (Exact, dict()), + (Davidson, dict(nroots=3)), + ], +) +def test_hf( + helper: Helper, mf: scf.hf.RHF, solver_cls: type[BaseSolver], solver_kwargs: dict[str, Any] +) -> None: """Test the HF expression.""" hf_h = HF.h.from_mf(mf) hf_p = HF.p.from_mf(mf) - hf_dyson = HF["dyson"].from_mf(mf) - gf_h_moments = hf_h.build_gf_moments(2) - gf_p_moments = hf_p.build_gf_moments(2) - gf_dyson_moments = hf_dyson.build_gf_moments(2) + gf_moments_h = hf_h.build_gf_moments(2) + gf_moments_p = hf_p.build_gf_moments(2) # Get the energy from the hole moments h1e = np.einsum("pq,pi,qj->ij", mf.get_hcore(), mf.mo_coeff, mf.mo_coeff) - energy = util.gf_moments_galitskii_migdal(gf_h_moments, h1e, factor=1.0) + energy = util.gf_moments_galitskii_migdal(gf_moments_h, h1e, factor=1.0) assert np.abs(energy - mf.energy_elec()[0]) < 1e-8 # Get the Fock matrix Fock matrix from the moments fock_ref = np.einsum("pq,pi,qj->ij", mf.get_fock(), mf.mo_coeff, mf.mo_coeff) - fock = gf_h_moments[1] + gf_p_moments[1] + fock = gf_moments_h[1] + gf_moments_p[1] assert np.allclose(fock, fock_ref) - assert np.allclose(gf_dyson_moments[1], fock) + assert np.allclose(gf_moments_h[1] + gf_moments_p[1], fock) + + # Get the hole Green's function from the Exact solver + solver_h = solver_cls.from_expression(hf_h, **solver_kwargs) + solver_h.kernel() + nocc = mf.mol.nelectron // 2 + + assert solver_h.result is not None + assert np.allclose( + solver_h.result.get_greens_function().as_perturbed_mo_energy()[nocc-1], mf.mo_energy[nocc-1] + ) + if solver_cls is Exact: + assert helper.have_equal_moments(solver_h.result.get_greens_function(), gf_moments_h, 2) + + # Get the particle Green's function + solver_p = solver_cls.from_expression(hf_p, **solver_kwargs) + solver_p.kernel() + + assert solver_p.result is not None + assert np.allclose( + solver_p.result.get_greens_function().as_perturbed_mo_energy()[nocc], mf.mo_energy[nocc] + ) + if solver_cls is Exact: + assert helper.have_equal_moments(solver_p.result.get_greens_function(), gf_moments_p, 2) + + # Get the combined Green's function + result = solver_h.result.combine(solver_p.result) + + assert np.allclose( + result.get_greens_function().as_perturbed_mo_energy()[nocc-1], mf.mo_energy[nocc-1] + ) + assert np.allclose( + result.get_greens_function().as_perturbed_mo_energy()[nocc], mf.mo_energy[nocc] + ) + if solver_cls is Exact: + assert np.allclose(result.get_greens_function().as_perturbed_mo_energy(), mf.mo_energy) + assert helper.have_equal_moments( + result.get_greens_function(), gf_moments_h + gf_moments_p, 2 + ) - # Get the Green's function from the Exact solver - exact_h = Exact.from_expression(hf_h) - exact_h.kernel() - exact_p = Exact.from_expression(hf_p) - exact_p.kernel() - assert exact_h.result is not None - assert exact_p.result is not None - result = exact_h.result.combine(exact_p.result) + # Check the zeroth moments add to identity + gf_moment_0 = gf_moments_h[0] + gf_moments_p[0] - assert np.allclose(result.get_greens_function().as_perturbed_mo_energy(), mf.mo_energy) + assert np.allclose(gf_moment_0, np.eye(mf.mol.nao)) -def test_ccsd(mf: scf.hf.RHF) -> None: +@pytest.mark.parametrize( + "solver_cls, solver_kwargs", + [ + (Exact, dict()), + (Davidson, dict(nroots=3)), + ], +) +def test_ccsd( + helper: Helper, mf: scf.hf.RHF, solver_cls: type[BaseSolver], solver_kwargs: dict[str, Any] +) -> None: """Test the CCSD expression.""" ccsd_h = CCSD.h.from_mf(mf) ccsd_p = CCSD.p.from_mf(mf) @@ -158,13 +209,35 @@ def test_ccsd(mf: scf.hf.RHF) -> None: else: assert correct_energy - # Get the Green's function from the Davidson solver - davidson = Davidson.from_expression(ccsd_h, nroots=3) - davidson.kernel() + # Get the hole Green's function from + solver_h = solver_cls.from_expression(ccsd_h, **solver_kwargs) + solver_h.kernel() ip_ref, _ = pyscf_ccsd.ipccsd(nroots=3) - assert davidson.result is not None - assert np.allclose(davidson.result.eigvals[0], -ip_ref[-1]) + assert solver_h.result is not None + assert np.allclose(solver_h.result.eigvals[-1], -ip_ref[0]) + if solver_cls is Exact: + assert helper.have_equal_moments(solver_h.result.get_greens_function(), gf_moments_h, 2) + + # Get the particle Green's function + solver_p = solver_cls.from_expression(ccsd_p, **solver_kwargs) + solver_p.kernel() + ea_ref, _ = pyscf_ccsd.eaccsd(nroots=3) + + assert solver_p.result is not None + assert np.allclose(solver_p.result.eigvals[0], ea_ref[0]) + if solver_cls is Exact: + assert helper.have_equal_moments(solver_p.result.get_greens_function(), gf_moments_p, 2) + + # Get the combined Green's function + result = solver_h.result.combine(solver_p.result) + + assert np.allclose(result.get_greens_function().physical().occupied().energies[-1], -ip_ref[0]) + assert np.allclose(result.get_greens_function().physical().virtual().energies[0], ea_ref[0]) + if solver_cls is Exact: + assert helper.have_equal_moments( + result.get_greens_function(), gf_moments_h + gf_moments_p, 2 + ) # Check the RDM rdm1 = gf_moments_h[0].copy() @@ -179,7 +252,16 @@ def test_ccsd(mf: scf.hf.RHF) -> None: assert np.allclose(gf_moment_0, np.eye(mf.mol.nao)) -def test_fci(mf: scf.hf.RHF) -> None: +@pytest.mark.parametrize( + "solver_cls, solver_kwargs", + [ + (Exact, dict()), + (Davidson, dict(nroots=3)), + ], +) +def test_fci( + helper: Helper, mf: scf.hf.RHF, solver_cls: type[BaseSolver], solver_kwargs: dict[str, Any] +) -> None: """Test the FCI expression.""" fci_h = FCI.h.from_mf(mf) fci_p = FCI.p.from_mf(mf) @@ -190,10 +272,51 @@ def test_fci(mf: scf.hf.RHF) -> None: # Get the energy from the hole moments h1e = np.einsum("pq,pi,qj->ij", mf.get_hcore(), mf.mo_coeff, mf.mo_coeff) energy = util.gf_moments_galitskii_migdal(gf_moments_h, h1e, factor=1.0) - energy_ref = pyscf_fci.kernel()[0] - mf.mol.energy_nuc() + energy_fci = pyscf_fci.kernel()[0] + energy_ref = energy_fci - mf.mol.energy_nuc() assert np.abs(energy - energy_ref) < 1e-8 + # Get the hole Green's function + solver_h = solver_cls.from_expression(fci_h, **solver_kwargs) + solver_h.kernel() + mol_h = mf.mol.copy() + mol_h.nelec = (mf.mol.nelec[0] - 1, mf.mol.nelec[1]) + pyscf_fci_h = pyscf.fci.FCI(mf.__class__(mol_h).run()) + energy_fci_h = pyscf_fci_h.kernel()[0] + energy_ref_h = energy_fci_h - mf.mol.energy_nuc() + ip_ref = energy_ref - energy_ref_h + + assert solver_h.result is not None + assert np.allclose(solver_h.result.eigvals[-1], ip_ref) + if solver_cls is Exact: + assert helper.have_equal_moments(solver_h.result.get_greens_function(), gf_moments_h, 2) + + # Get the particle Green's function + solver_p = solver_cls.from_expression(fci_p, **solver_kwargs) + solver_p.kernel() + mol_p = mf.mol.copy() + mol_p.nelec = (mf.mol.nelec[0] + 1, mf.mol.nelec[1]) + pyscf_fci_p = pyscf.fci.FCI(mf.__class__(mol_p).run()) + energy_fci_p = pyscf_fci_p.kernel()[0] + energy_ref_p = energy_fci_p - mf.mol.energy_nuc() + ea_ref = energy_ref_p - energy_ref + + assert solver_p.result is not None + assert np.allclose(solver_p.result.eigvals[0], ea_ref) + if solver_cls is Exact: + assert helper.have_equal_moments(solver_p.result.get_greens_function(), gf_moments_p, 2) + + # Get the combined Green's function + result = solver_h.result.combine(solver_p.result) + + assert np.allclose(result.get_greens_function().physical().occupied().energies[-1], ip_ref) + assert np.allclose(result.get_greens_function().physical().virtual().energies[0], ea_ref) + if solver_cls is Exact: + assert helper.have_equal_moments( + result.get_greens_function(), gf_moments_h + gf_moments_p, 2 + ) + # Check the RDM rdm1 = fci_h.build_gf_moments(1)[0] * 2 rdm1_ref = pyscf_fci.make_rdm1(pyscf_fci.ci, mf.mol.nao, mf.mol.nelectron) @@ -214,28 +337,61 @@ def test_fci(mf: scf.hf.RHF) -> None: assert np.allclose(gf_moments_ccsd[1], gf_moments_h[1]) -def test_adc2(mf: scf.hf.RHF) -> None: +@pytest.mark.parametrize( + "solver_cls, solver_kwargs", + [ + (Exact, dict()), + (Davidson, dict(nroots=3)), + ], +) +def test_adc2( + helper: Helper, mf: scf.hf.RHF, solver_cls: type[BaseSolver], solver_kwargs: dict[str, Any] +) -> None: """Test the ADC(2) expression.""" adc_h = ADC2.h.from_mf(mf) adc_p = ADC2.p.from_mf(mf) - pyscf_adc = pyscf.adc.ADC(mf) + pyscf_adc_ip = pyscf.adc.ADC(mf) + pyscf_adc_ea = pyscf.adc.ADC(mf) + pyscf_adc_ea.method_type = "ea" gf_moments_h = adc_h.build_gf_moments(2) gf_moments_p = adc_p.build_gf_moments(2) # Get the energy from the hole moments h1e = np.einsum("pq,pi,qj->ij", mf.get_hcore(), mf.mo_coeff, mf.mo_coeff) energy = util.gf_moments_galitskii_migdal(gf_moments_h, h1e, factor=1.0) - energy_ref = mf.energy_elec()[0] + pyscf_adc.kernel_gs()[0] + energy_ref = mf.energy_elec()[0] + pyscf_adc_ip.kernel_gs()[0] assert np.abs(energy - energy_ref) < 1e-8 - # Get the Green's function from the Davidson solver - davidson = Davidson.from_expression(adc_h, nroots=3) - davidson.kernel() - ip_ref, _, _, _ = pyscf_adc.kernel(nroots=3) + # Get the hole Green's function from the Davidson solver + solver_h = solver_cls.from_expression(adc_h, **solver_kwargs) + solver_h.kernel() + ip_ref, _, _, _ = pyscf_adc_ip.kernel(nroots=3) + + assert solver_h.result is not None + assert np.allclose(solver_h.result.eigvals[-1], -ip_ref[0]) + if solver_cls is Exact: + assert helper.have_equal_moments(solver_h.result.get_greens_function(), gf_moments_h, 2) - assert davidson.result is not None - assert np.allclose(davidson.result.eigvals[0], -ip_ref[-1]) + # Get the particle Green's function from the Davidson solver + solver_p = solver_cls.from_expression(adc_p, **solver_kwargs) + solver_p.kernel() + ea_ref, _, _, _ = pyscf_adc_ea.kernel(nroots=3) + + assert solver_p.result is not None + assert np.allclose(solver_p.result.eigvals[0], ea_ref[0]) + if solver_cls is Exact: + assert helper.have_equal_moments(solver_p.result.get_greens_function(), gf_moments_p, 2) + + # Get the combined Green's function from the Davidson solver + result = solver_h.result.combine(solver_p.result) + + assert np.allclose(result.get_greens_function().physical().occupied().energies[-1], -ip_ref[0]) + assert np.allclose(result.get_greens_function().physical().virtual().energies[0], ea_ref[0]) + if solver_cls is Exact: + assert helper.have_equal_moments( + result.get_greens_function(), gf_moments_h + gf_moments_p, 2 + ) # Check the RDM rdm1 = adc_h.build_gf_moments(1)[0] * 2 @@ -249,29 +405,63 @@ def test_adc2(mf: scf.hf.RHF) -> None: assert np.allclose(gf_moment_0, np.eye(mf.mol.nao)) -def test_adc2x(mf: scf.hf.RHF) -> None: +@pytest.mark.parametrize( + "solver_cls, solver_kwargs", + [ + (Exact, dict()), + (Davidson, dict(nroots=3)), + ], +) +def test_adc2x( + helper: Helper, mf: scf.hf.RHF, solver_cls: type[BaseSolver], solver_kwargs: dict[str, Any] +) -> None: """Test the ADC(2)-x expression.""" adc_h = ADC2x.h.from_mf(mf) adc_p = ADC2x.p.from_mf(mf) - pyscf_adc = pyscf.adc.ADC(mf) - pyscf_adc.method = "adc(2)-x" + pyscf_adc_ip = pyscf.adc.ADC(mf) + pyscf_adc_ip.method = "adc(2)-x" + pyscf_adc_ea = pyscf.adc.ADC(mf) + pyscf_adc_ea.method = "adc(2)-x" + pyscf_adc_ea.method_type = "ea" gf_moments_h = adc_h.build_gf_moments(2) gf_moments_p = adc_p.build_gf_moments(2) # Get the energy from the hole moments h1e = np.einsum("pq,pi,qj->ij", mf.get_hcore(), mf.mo_coeff, mf.mo_coeff) energy = util.gf_moments_galitskii_migdal(gf_moments_h, h1e, factor=1.0) - energy_ref = mf.energy_elec()[0] + pyscf_adc.kernel_gs()[0] + energy_ref = mf.energy_elec()[0] + pyscf_adc_ip.kernel_gs()[0] assert np.abs(energy - energy_ref) < 1e-8 - # Get the Green's function from the Davidson solver - davidson = Davidson.from_expression(adc_h, nroots=3) - davidson.kernel() - ip_ref, _, _, _ = pyscf_adc.kernel(nroots=3) - - assert davidson.result is not None - assert np.allclose(davidson.result.eigvals[0], -ip_ref[-1]) + # Get the hole Green's function from the Davidson solver + solver_h = solver_cls.from_expression(adc_h, **solver_kwargs) + solver_h.kernel() + ip_ref, _, _, _ = pyscf_adc_ip.kernel(nroots=3) + + assert solver_h.result is not None + assert np.allclose(solver_h.result.eigvals[-1], -ip_ref[0]) + if solver_cls is Exact: + assert helper.have_equal_moments(solver_h.result.get_greens_function(), gf_moments_h, 2) + + # Get the particle Green's function from the Davidson solver + solver_p = solver_cls.from_expression(adc_p, **solver_kwargs) + solver_p.kernel() + ea_ref, _, _, _ = pyscf_adc_ea.kernel(nroots=3) + + assert solver_p.result is not None + assert np.allclose(solver_p.result.eigvals[0], ea_ref[0]) + if solver_cls is Exact: + assert helper.have_equal_moments(solver_p.result.get_greens_function(), gf_moments_p, 2) + + # Get the combined Green's function from the Davidson solver + result = solver_h.result.combine(solver_p.result) + + assert np.allclose(result.get_greens_function().physical().occupied().energies[-1], -ip_ref[0]) + assert np.allclose(result.get_greens_function().physical().virtual().energies[0], ea_ref[0]) + if solver_cls is Exact: + assert helper.have_equal_moments( + result.get_greens_function(), gf_moments_h + gf_moments_p, 2 + ) # Check the RDM rdm1 = adc_h.build_gf_moments(1)[0] * 2 From fe944eff73a6d5da22fcb6cf04c25f0d29efbe93 Mon Sep 17 00:00:00 2001 From: Oliver Backhouse Date: Wed, 17 Dec 2025 18:29:25 +0000 Subject: [PATCH 04/17] Formatting, introduce typing-extensions for dev --- dyson/expressions/expression.py | 7 +++---- dyson/solvers/solver.py | 9 +++++---- dyson/solvers/static/mblse.py | 4 ++-- pyproject.toml | 1 + tests/test_expressions.py | 17 +++++++++-------- 5 files changed, 20 insertions(+), 18 deletions(-) diff --git a/dyson/expressions/expression.py b/dyson/expressions/expression.py index f85b299..a2c0bf8 100644 --- a/dyson/expressions/expression.py +++ b/dyson/expressions/expression.py @@ -12,6 +12,7 @@ if TYPE_CHECKING: from typing import Callable + from typing_extensions import Self from pyscf.gto.mole import Mole from pyscf.scf.hf import RHF @@ -33,7 +34,7 @@ class BaseExpression(ABC): @classmethod @abstractmethod - def from_mf(cls, mf: RHF) -> BaseExpression: + def from_mf(cls, mf: RHF) -> Self: """Create an expression from a mean-field object. Args: @@ -476,9 +477,7 @@ def __getattr__(cls, key: str) -> type[BaseExpression]: @property def _classes(cls) -> set[type[BaseExpression]]: """Get all classes in the collection.""" - return { - cls for cls in [cls._hole, cls._particle, cls._neutral] if cls is not None - } + return {cls for cls in [cls._hole, cls._particle, cls._neutral] if cls is not None} def __contains__(cls, key: str) -> bool: """Check if an expression exists by its name.""" diff --git a/dyson/solvers/solver.py b/dyson/solvers/solver.py index d5df977..6b233fe 100644 --- a/dyson/solvers/solver.py +++ b/dyson/solvers/solver.py @@ -15,6 +15,7 @@ if TYPE_CHECKING: from typing import Any + from typing_extensions import Self from dyson.expressions.expression import BaseExpression from dyson.representations.dynamic import Dynamic @@ -32,7 +33,7 @@ def __init_subclass__(cls, *args: Any, **kwargs: Any) -> None: def wrap_init(init: Any) -> Any: """Wrapper to call __post_init__ after __init__.""" - def wrapped_init(self: BaseSolver, *args: Any, **kwargs: Any) -> None: + def wrapped_init(self: Self, *args: Any, **kwargs: Any) -> None: init(self, *args, **kwargs) if init.__name__ == "__init__": self.__log_init__() @@ -43,7 +44,7 @@ def wrapped_init(self: BaseSolver, *args: Any, **kwargs: Any) -> None: def wrap_kernel(kernel: Any) -> Any: """Wrapper to call __post_kernel__ after kernel.""" - def wrapped_kernel(self: BaseSolver, *args: Any, **kwargs: Any) -> Any: + def wrapped_kernel(self: Self, *args: Any, **kwargs: Any) -> Any: result = kernel(self, *args, **kwargs) if kernel.__name__ == "kernel": self.__post_kernel__() @@ -112,7 +113,7 @@ def from_self_energy( self_energy: Lehmann, overlap: Array | None = None, **kwargs: Any, - ) -> BaseSolver: + ) -> Self: """Create a solver from a self-energy. Args: @@ -132,7 +133,7 @@ def from_self_energy( @classmethod @abstractmethod - def from_expression(cls, expression: BaseExpression, **kwargs: Any) -> BaseSolver: + def from_expression(cls, expression: BaseExpression, **kwargs: Any) -> Self: """Create a solver from an expression. Args: diff --git a/dyson/solvers/static/mblse.py b/dyson/solvers/static/mblse.py index e2108c8..266b024 100644 --- a/dyson/solvers/static/mblse.py +++ b/dyson/solvers/static/mblse.py @@ -501,7 +501,7 @@ def from_self_energy( self_energy: Lehmann, overlap: Array | None = None, **kwargs: Any, - ) -> MBLSE: + ) -> MLSE: """Create a solver from a self-energy. Args: @@ -516,7 +516,7 @@ def from_self_energy( raise NotImplementedError("Cannot instantiate MLSE from a self-energy.") @classmethod - def from_expression(cls, expression: BaseExpression, **kwargs: Any) -> MBLSE: + def from_expression(cls, expression: BaseExpression, **kwargs: Any) -> MLSE: """Create a solver from an expression. Args: diff --git a/pyproject.toml b/pyproject.toml index d5b864e..6839af0 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -37,6 +37,7 @@ dependencies = [ dev = [ "ruff>=0.10.0", "mypy>=1.5.0", + "typing-extensions>=4.0.0", "coverage[toml]>=5.5.0", "pytest>=6.2.4", "pytest-cov>=4.0.0", diff --git a/tests/test_expressions.py b/tests/test_expressions.py index e93ff77..f4c7989 100644 --- a/tests/test_expressions.py +++ b/tests/test_expressions.py @@ -19,7 +19,7 @@ from pyscf import scf from dyson.expressions.expression import BaseExpression - from dyson.solvers.solver import BaseSolver + from dyson.solvers.solver import StaticSolver from .conftest import ExactGetter, Helper @@ -113,7 +113,7 @@ def test_static( ], ) def test_hf( - helper: Helper, mf: scf.hf.RHF, solver_cls: type[BaseSolver], solver_kwargs: dict[str, Any] + helper: Helper, mf: scf.hf.RHF, solver_cls: type[StaticSolver], solver_kwargs: dict[str, Any] ) -> None: """Test the HF expression.""" hf_h = HF.h.from_mf(mf) @@ -141,7 +141,8 @@ def test_hf( assert solver_h.result is not None assert np.allclose( - solver_h.result.get_greens_function().as_perturbed_mo_energy()[nocc-1], mf.mo_energy[nocc-1] + solver_h.result.get_greens_function().as_perturbed_mo_energy()[nocc - 1], + mf.mo_energy[nocc - 1], ) if solver_cls is Exact: assert helper.have_equal_moments(solver_h.result.get_greens_function(), gf_moments_h, 2) @@ -161,7 +162,7 @@ def test_hf( result = solver_h.result.combine(solver_p.result) assert np.allclose( - result.get_greens_function().as_perturbed_mo_energy()[nocc-1], mf.mo_energy[nocc-1] + result.get_greens_function().as_perturbed_mo_energy()[nocc - 1], mf.mo_energy[nocc - 1] ) assert np.allclose( result.get_greens_function().as_perturbed_mo_energy()[nocc], mf.mo_energy[nocc] @@ -186,7 +187,7 @@ def test_hf( ], ) def test_ccsd( - helper: Helper, mf: scf.hf.RHF, solver_cls: type[BaseSolver], solver_kwargs: dict[str, Any] + helper: Helper, mf: scf.hf.RHF, solver_cls: type[StaticSolver], solver_kwargs: dict[str, Any] ) -> None: """Test the CCSD expression.""" ccsd_h = CCSD.h.from_mf(mf) @@ -260,7 +261,7 @@ def test_ccsd( ], ) def test_fci( - helper: Helper, mf: scf.hf.RHF, solver_cls: type[BaseSolver], solver_kwargs: dict[str, Any] + helper: Helper, mf: scf.hf.RHF, solver_cls: type[StaticSolver], solver_kwargs: dict[str, Any] ) -> None: """Test the FCI expression.""" fci_h = FCI.h.from_mf(mf) @@ -345,7 +346,7 @@ def test_fci( ], ) def test_adc2( - helper: Helper, mf: scf.hf.RHF, solver_cls: type[BaseSolver], solver_kwargs: dict[str, Any] + helper: Helper, mf: scf.hf.RHF, solver_cls: type[StaticSolver], solver_kwargs: dict[str, Any] ) -> None: """Test the ADC(2) expression.""" adc_h = ADC2.h.from_mf(mf) @@ -413,7 +414,7 @@ def test_adc2( ], ) def test_adc2x( - helper: Helper, mf: scf.hf.RHF, solver_cls: type[BaseSolver], solver_kwargs: dict[str, Any] + helper: Helper, mf: scf.hf.RHF, solver_cls: type[StaticSolver], solver_kwargs: dict[str, Any] ) -> None: """Test the ADC(2)-x expression.""" adc_h = ADC2x.h.from_mf(mf) From 13c3f94a2b177f2978b2191d37501ce041f9bb31 Mon Sep 17 00:00:00 2001 From: Oliver Backhouse Date: Wed, 17 Dec 2025 18:38:31 +0000 Subject: [PATCH 05/17] ruff --- dyson/expressions/expression.py | 2 +- dyson/solvers/solver.py | 1 + 2 files changed, 2 insertions(+), 1 deletion(-) diff --git a/dyson/expressions/expression.py b/dyson/expressions/expression.py index a2c0bf8..0e0972e 100644 --- a/dyson/expressions/expression.py +++ b/dyson/expressions/expression.py @@ -12,10 +12,10 @@ if TYPE_CHECKING: from typing import Callable - from typing_extensions import Self from pyscf.gto.mole import Mole from pyscf.scf.hf import RHF + from typing_extensions import Self from dyson.typing import Array diff --git a/dyson/solvers/solver.py b/dyson/solvers/solver.py index 6b233fe..21945de 100644 --- a/dyson/solvers/solver.py +++ b/dyson/solvers/solver.py @@ -15,6 +15,7 @@ if TYPE_CHECKING: from typing import Any + from typing_extensions import Self from dyson.expressions.expression import BaseExpression From 82c074cc2a1fab4997e1ce593a0a6fb60ff4028a Mon Sep 17 00:00:00 2001 From: Oliver Backhouse Date: Fri, 23 Jan 2026 17:52:08 +0000 Subject: [PATCH 06/17] Fix spectral combination with target observables --- dyson/expressions/expression.py | 19 +++ dyson/representations/spectral.py | 235 +++++++++++++++++++++++------- dyson/solvers/static/mblse.py | 1 + tests/conftest.py | 2 +- tests/test_davidson.py | 11 +- tests/test_density.py | 4 +- tests/test_downfolded.py | 4 +- tests/test_exact.py | 2 +- tests/test_expressions.py | 32 ++-- tests/test_mblgf.py | 6 +- tests/test_mblse.py | 16 +- 11 files changed, 248 insertions(+), 84 deletions(-) diff --git a/dyson/expressions/expression.py b/dyson/expressions/expression.py index 0e0972e..d7cf5d4 100644 --- a/dyson/expressions/expression.py +++ b/dyson/expressions/expression.py @@ -399,6 +399,25 @@ def build_se_moments(self, nmom: int, reduction: Reduction = Reduction.NONE) -> """ pass + def build_overlap(self, reduction: Reduction = Reduction.NONE) -> Array: + """Build the matrix representing the overlap of the excitation vectors. + + Args: + reduction: Reduction to apply to the overlap. + + Returns: + Overlap matrix. + + Notes: + The overlap matrix is equal to the zeroth moment of the Green's function. + """ + return self.build_gf_moments( + 1, + store_vectors=True, + left=False, + reduction=reduction, + )[0] + @property @abstractmethod def mol(self) -> Mole: diff --git a/dyson/representations/spectral.py b/dyson/representations/spectral.py index 881561c..5d24e2c 100644 --- a/dyson/representations/spectral.py +++ b/dyson/representations/spectral.py @@ -96,6 +96,7 @@ def from_self_energy( static: Array, self_energy: Lehmann, overlap: Array | None = None, + sort: bool = False, ) -> Spectral: """Create a spectrum from a self-energy. @@ -103,6 +104,7 @@ def from_self_energy( static: Static part of the self-energy. self_energy: Self-energy. overlap: Overlap matrix for the physical space. + sort: Sort the eigenfunctions by eigenvalue. Returns: Spectrum object. @@ -111,6 +113,7 @@ def from_self_energy( *self_energy.diagonalise_matrix(static, overlap=overlap), self_energy.nphys, chempot=self_energy.chempot, + sort=sort, ) def sort_(self) -> None: @@ -238,69 +241,49 @@ def get_greens_function(self, chempot: float | None = None) -> Lehmann: chempot = 0.0 return Lehmann(*self.get_dyson_orbitals(), chempot=chempot) - def combine(self, *args: Spectral, chempot: float | None = None) -> Spectral: - """Combine multiple spectral representations. + @classmethod + def combine_for_greens_function( + cls, + *spectrals: Spectral, + overlap: Array | None = None, + chempot: float | None = None, + ) -> Spectral: + """Combine multiple spectral representations to conserve the Green's functions. Args: - args: Spectral representations to combine. - chempot: Chemical potential to be used in the Lehmann representations of the self-energy - and Green's function. + *spectrals: Spectral representations to combine. + overlap: Overlap matrix for the physical space. If not passed, the overlap matrices from the + individual spectral representations will be summed. + chempot: Chemical potential to use for the combined result. If not provided, the + chemical potential from the individual results will be used, raising an exception + if they are inconsistent. Returns: Combined spectral representation. """ - args = (self, *args) - if len(set(arg.nphys for arg in args)) != 1: - raise ValueError( - "All Spectral objects must have the same number of physical degrees of freedom." - ) - nphys = args[0].nphys - hermitian = all(arg.hermitian for arg in args) - - # Check the chemical potentials - if chempot is None: - if any(arg.chempot is not None for arg in args): - chempots = [arg.chempot for arg in args if arg.chempot is not None] - if not all(np.isclose(chempots[0], part) for part in chempots[1:]): - raise ValueError( - "If chempot is not passed to combine, all chemical potentials must be " - "equal in the inputs." - ) - chempot = chempots[0] - - # Get the eigenvalues and dyson orbitals - eigvals = np.concatenate([arg.eigvals for arg in args]) - left = np.zeros((nphys, 0)) - right = np.zeros((nphys, 0)) - for arg in args: - orbitals = arg.get_dyson_orbitals()[1] - left_i, right_i = util.unpack_vectors(orbitals) - left = np.concatenate([left, left_i], axis=1) - right = np.concatenate([right, right_i], axis=1) - - # Orthogonalise under the metric of the overlap matrix - overlap = sum([arg.get_overlap() for arg in args], np.zeros((nphys, nphys))) - orth, _ = util.matrix_power(overlap, -0.5, hermitian=hermitian, return_error=False) - left = util.rotate_subspace(left, orth) - right = util.rotate_subspace(right, orth.T.conj()) - - # Construct the eigenvectors - eigvecs = util.project_eigenvectors( - None, - left, - right if not hermitian else None, - ) + return combine_for_greens_function(*spectrals, overlap=overlap, chempot=chempot) - # Unorthogonalise the eigenvectors - unorth, _ = util.matrix_power(overlap, 0.5, hermitian=hermitian, return_error=False) - left = util.rotate_subspace(left, unorth.T.conj()) - right = util.rotate_subspace(right, unorth) - eigvecs = np.array([left, right]) + @classmethod + def combine_for_self_energy( + cls, + *spectrals: Spectral, + overlap: Array | None = None, + chempot: float | None = None, + ) -> Spectral: + """Combine multiple spectral representations to conserve the self-energies. - # Construct the result - result = Spectral(eigvals, eigvecs, nphys, chempot=chempot, sort=True) + Args: + *spectrals: Spectral representations to combine. + overlap: Overlap matrix for the physical space. If not passed, the overlap matrices from the + individual spectral representations will be summed. + chempot: Chemical potential to use for the combined result. If not provided, the + chemical potential from the individual results will be used, raising an exception + if they are inconsistent. - return result + Returns: + Combined spectral representation. + """ + return combine_for_self_energy(*spectrals, overlap=overlap, chempot=chempot) @cached_property def overlap(self) -> Array: @@ -351,3 +334,147 @@ def __eq__(self, other: object) -> bool: def __hash__(self) -> int: """Hash the object.""" return hash((tuple(self.eigvals), tuple(self.eigvecs.flatten()), self.nphys, self.chempot)) + + +def _get_physical_space_size(spectrals: tuple[Spectral, ...]) -> int: + """Get the size of the physical space from multiple spectral representations. + + Args: + spectrals: Spectral representations. + + Returns: + Size of the physical space. + + Raises: + ValueError: If the physical space sizes are inconsistent. + """ + nphys_set = {s.nphys for s in spectrals} + if len(nphys_set) != 1: + raise ValueError("Inconsistent physical space sizes in spectral representations.") + return nphys_set.pop() + + +def _get_chemical_potential(spectrals: tuple[Spectral, ...]) -> float | None: + """Get the chemical potential from multiple spectral representations. + + Args: + spectrals: Spectral representations. + + Returns: + Chemical potential. + + Raises: + ValueError: If the chemical potentials are inconsistent. + """ + chempot_set = {s.chempot for s in spectrals} + chempot_set.discard(None) + if not all(np.isclose(c1, c2) for c1 in chempot_set for c2 in chempot_set): + raise ValueError("Inconsistent chemical potentials in spectral representations.") + return chempot_set.pop() if chempot_set else None + + +def combine_for_greens_function( + *spectrals: Spectral, overlap: Array | None = None, chempot: float | None = None +) -> Spectral: + """Combine multiple spectral representations to conserve the Green's functions. + + Args: + *spectrals: Spectral representations to combine. + overlap: Overlap matrix for the physical space. If not passed, the overlap matrices from the + individual spectral representations will be summed. + chempot: Chemical potential to use for the combined result. If not provided, the chemical + potential from the individual results will be used, raising an exception if they are + inconsistent. + + Returns: + Combined spectral representation. + + Note: + This function combines the spectral representations in such a way that the resulting Green's + function is equivalent to the combination of the individual Green's functions. This is not + guaranteed to conserve the self-energy in the same way. + """ + if len(spectrals) == 1: + return spectrals[0] + nphys = _get_physical_space_size(spectrals) + hermitian = all(spectral.hermitian for spectral in spectrals) + if chempot is None: + chempot = _get_chemical_potential(spectrals) + if overlap is None: + overlap = sum(spectral.overlap for spectral in spectrals) + + # Combine the eigenvalues and Dyson orbitals + eigvals = np.concatenate([spectral.eigvals for spectral in spectrals], axis=0) + left = np.zeros((nphys, 0)) + right = np.zeros((nphys, 0)) + for spectral in spectrals: + orbitals = spectral.get_dyson_orbitals()[1] + l, r = util.unpack_vectors(orbitals) + left = np.concatenate((left, l), axis=1) + right = np.concatenate((right, r), axis=1) + + # Get the orthogonalisation matrices + orth, _ = util.matrix_power(overlap, -0.5, hermitian=hermitian, return_error=False) + unorth, _ = util.matrix_power(overlap, 0.5, hermitian=hermitian, return_error=False) + + # Find a basis for the null space + projector = right.T @ left.conj() + vectors = util.null_space_basis(projector, hermitian=hermitian) + + # Orthogonalise the vectors + left = left.T @ orth + right = right.T @ orth.T.conj() + + # Biorthonormalise full set of vectors + left = np.concatenate([left, vectors[0]], axis=1) + right = np.concatenate([right, vectors[1]], axis=1) + left, right = util.biorthonormalise(left, right) + + # Unorthogonalise the vectors + left[:, : nphys] = left[:, : nphys] @ unorth.T.conj() + right[:, : nphys] = right[:, : nphys] @ unorth + + # Get the combined result + eigvecs = np.array([left.T.conj(), right.T.conj()]) if not hermitian else left.T.conj() + spectral = Spectral(eigvals, eigvecs, nphys, chempot=chempot, sort=True) + + return spectral + + +def combine_for_self_energy( + *spectrals: Spectral, overlap: Array | None = None, chempot: float | None = None +) -> Spectral: + """Combine multiple spectral representations to conserve the self-energies. + + Args: + *spectrals: Spectral representations to combine. + overlap: Overlap matrix for the physical space. + chempot: Chemical potential to use for the combined result. If not provided, the chemical + potential from the individual results will be used, raising an exception if they are + inconsistent. + + Returns: + Combined spectral representation. + + Note: + This function combines the spectral representations in such a way that the resulting + self-energy is equivalent to the combination of the individual self-energies. This is not + guaranteed to conserve the Green's function in the same way. + """ + if len(spectrals) == 1: + return spectrals[0] + if chempot is None: + chempot = _get_chemical_potential(spectrals) + if overlap is None: + overlap = sum(spectral.overlap for spectral in spectrals) + + # Combine the self-energies + static = sum(spectral.get_static_self_energy() for spectral in spectrals) + self_energy = spectrals[0].get_self_energy(chempot=chempot).copy() + for spectral in spectrals[1:]: + self_energy = self_energy.concatenate(spectral.get_self_energy(chempot=chempot)) + + # Create the combined spectral representation + spectral = Spectral.from_self_energy(static, self_energy, overlap=overlap, sort=True) + + return spectral diff --git a/dyson/solvers/static/mblse.py b/dyson/solvers/static/mblse.py index 266b024..194a2ff 100644 --- a/dyson/solvers/static/mblse.py +++ b/dyson/solvers/static/mblse.py @@ -398,6 +398,7 @@ def _recurrence_iteration_non_hermitian( return error_sqrt, error_inv_sqrt, error_moments + #FIXME: being called twice? def solve(self, iteration: int | None = None) -> Spectral: """Solve the eigenvalue problem at a given iteration. diff --git a/tests/conftest.py b/tests/conftest.py index 147558c..68f9dbb 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -177,4 +177,4 @@ def _get_central_result( exact_p = exact_cache(mf, expression_method.p) assert exact_h.result is not None assert exact_p.result is not None - return Spectral.combine(exact_h.result, exact_p.result) + return Spectral.combine_for_greens_function(exact_h.result, exact_p.result) diff --git a/tests/test_davidson.py b/tests/test_davidson.py index 375d23f..b3cd0a0 100644 --- a/tests/test_davidson.py +++ b/tests/test_davidson.py @@ -9,7 +9,7 @@ from dyson.representations.lehmann import Lehmann from dyson.representations.spectral import Spectral -from dyson.solvers import Davidson +from dyson.solvers import Davidson, Exact if TYPE_CHECKING: from pyscf import scf @@ -149,8 +149,13 @@ def test_vs_exact_solver_central( assert helper.have_equal_moments(self_energy, self_energy_exact, 2) # Use the component-wise solvers - result_exact = Spectral.combine(exact_h.result, exact_p.result) - result_davidson = Spectral.combine(davidson_h.result, davidson_p.result) + overlap = expression_h.build_overlap() + expression_p.build_overlap() + result_exact = Spectral.combine_for_greens_function( + exact_h.result, exact_p.result, overlap=overlap + ) + result_davidson = Spectral.combine_for_greens_function( + exact_h.result, exact_p.result, overlap=overlap + ) # Get the self-energy and Green's function from the Davidson solver static = result_davidson.get_static_self_energy() diff --git a/tests/test_density.py b/tests/test_density.py index dafb6e7..1b1b805 100644 --- a/tests/test_density.py +++ b/tests/test_density.py @@ -8,7 +8,7 @@ import pytest from dyson.representations.spectral import Spectral -from dyson.solvers import DensityRelaxation +from dyson.solvers import DensityRelaxation, Exact from dyson.solvers.static.density import get_fock_matrix_function if TYPE_CHECKING: @@ -36,7 +36,7 @@ def test_vs_exact_solver( exact_p = exact_cache(mf, expression_method.p) assert exact_h.result is not None assert exact_p.result is not None - result_exact = Spectral.combine(exact_h.result, exact_p.result) + result_exact = Spectral.combine_for_greens_function(exact_h.result, exact_p.result) # Solve the Hamiltonian with DensityRelaxation get_fock = get_fock_matrix_function(mf) diff --git a/tests/test_downfolded.py b/tests/test_downfolded.py index 57c4de5..bc6c719 100644 --- a/tests/test_downfolded.py +++ b/tests/test_downfolded.py @@ -8,7 +8,7 @@ import pytest from dyson.representations.spectral import Spectral -from dyson.solvers import Downfolded +from dyson.solvers import Downfolded, Exact if TYPE_CHECKING: from pyscf import scf @@ -36,7 +36,7 @@ def test_vs_exact_solver( exact_p = exact_cache(mf, expression_method.p) assert exact_h.result is not None assert exact_p.result is not None - result_exact = Spectral.combine(exact_h.result, exact_p.result) + result_exact = Spectral.combine_for_greens_function(exact_h.result, exact_p.result) # Solve the Hamiltonian with Downfolded downfolded = Downfolded.from_self_energy( diff --git a/tests/test_exact.py b/tests/test_exact.py index dd514d6..ef74977 100644 --- a/tests/test_exact.py +++ b/tests/test_exact.py @@ -79,7 +79,7 @@ def test_vs_exact_solver_central( exact_p = exact_cache(mf, expression_method.p) assert exact_h.result is not None assert exact_p.result is not None - result_ph = Spectral.combine(exact_h.result, exact_p.result) + result_ph = Spectral.combine_for_greens_function(exact_h.result, exact_p.result) # Recover the hole self-energy and Green's function static = exact_h.result.get_static_self_energy() diff --git a/tests/test_expressions.py b/tests/test_expressions.py index f4c7989..67e8121 100644 --- a/tests/test_expressions.py +++ b/tests/test_expressions.py @@ -12,6 +12,7 @@ from dyson import util from dyson.expressions import ADC2, CCSD, FCI, HF, ADC2x from dyson.solvers import Davidson, Exact +from dyson.representations.spectral import Spectral if TYPE_CHECKING: from typing import Any @@ -159,7 +160,8 @@ def test_hf( assert helper.have_equal_moments(solver_p.result.get_greens_function(), gf_moments_p, 2) # Get the combined Green's function - result = solver_h.result.combine(solver_p.result) + overlap = hf_h.build_overlap() + hf_p.build_overlap() + result = Spectral.combine_for_greens_function(solver_h.result, solver_p.result, overlap=overlap) assert np.allclose( result.get_greens_function().as_perturbed_mo_energy()[nocc - 1], mf.mo_energy[nocc - 1] @@ -231,10 +233,11 @@ def test_ccsd( assert helper.have_equal_moments(solver_p.result.get_greens_function(), gf_moments_p, 2) # Get the combined Green's function - result = solver_h.result.combine(solver_p.result) + overlap = ccsd_h.build_overlap() + ccsd_p.build_overlap() + result = Spectral.combine_for_greens_function(solver_h.result, solver_p.result, overlap=overlap) - assert np.allclose(result.get_greens_function().physical().occupied().energies[-1], -ip_ref[0]) - assert np.allclose(result.get_greens_function().physical().virtual().energies[0], ea_ref[0]) + assert np.allclose(result.get_greens_function().occupied().energies[-1], -ip_ref[0]) + assert np.allclose(result.get_greens_function().virtual().energies[0], ea_ref[0]) if solver_cls is Exact: assert helper.have_equal_moments( result.get_greens_function(), gf_moments_h + gf_moments_p, 2 @@ -309,10 +312,11 @@ def test_fci( assert helper.have_equal_moments(solver_p.result.get_greens_function(), gf_moments_p, 2) # Get the combined Green's function - result = solver_h.result.combine(solver_p.result) + overlap = fci_h.build_overlap() + fci_p.build_overlap() + result = Spectral.combine_for_greens_function(solver_h.result, solver_p.result, overlap=overlap) - assert np.allclose(result.get_greens_function().physical().occupied().energies[-1], ip_ref) - assert np.allclose(result.get_greens_function().physical().virtual().energies[0], ea_ref) + assert np.allclose(result.get_greens_function().occupied().energies[-1], ip_ref) + assert np.allclose(result.get_greens_function().virtual().energies[0], ea_ref) if solver_cls is Exact: assert helper.have_equal_moments( result.get_greens_function(), gf_moments_h + gf_moments_p, 2 @@ -385,10 +389,11 @@ def test_adc2( assert helper.have_equal_moments(solver_p.result.get_greens_function(), gf_moments_p, 2) # Get the combined Green's function from the Davidson solver - result = solver_h.result.combine(solver_p.result) + overlap = adc_h.build_overlap() + adc_p.build_overlap() + result = Spectral.combine_for_greens_function(solver_h.result, solver_p.result, overlap=overlap) - assert np.allclose(result.get_greens_function().physical().occupied().energies[-1], -ip_ref[0]) - assert np.allclose(result.get_greens_function().physical().virtual().energies[0], ea_ref[0]) + assert np.allclose(result.get_greens_function().occupied().energies[-1], -ip_ref[0]) + assert np.allclose(result.get_greens_function().virtual().energies[0], ea_ref[0]) if solver_cls is Exact: assert helper.have_equal_moments( result.get_greens_function(), gf_moments_h + gf_moments_p, 2 @@ -455,10 +460,11 @@ def test_adc2x( assert helper.have_equal_moments(solver_p.result.get_greens_function(), gf_moments_p, 2) # Get the combined Green's function from the Davidson solver - result = solver_h.result.combine(solver_p.result) + overlap = adc_h.build_overlap() + adc_p.build_overlap() + result = Spectral.combine_for_greens_function(solver_h.result, solver_p.result, overlap=overlap) - assert np.allclose(result.get_greens_function().physical().occupied().energies[-1], -ip_ref[0]) - assert np.allclose(result.get_greens_function().physical().virtual().energies[0], ea_ref[0]) + assert np.allclose(result.get_greens_function().occupied().energies[-1], -ip_ref[0]) + assert np.allclose(result.get_greens_function().virtual().energies[0], ea_ref[0]) if solver_cls is Exact: assert helper.have_equal_moments( result.get_greens_function(), gf_moments_h + gf_moments_p, 2 diff --git a/tests/test_mblgf.py b/tests/test_mblgf.py index 611b906..4e60dbc 100644 --- a/tests/test_mblgf.py +++ b/tests/test_mblgf.py @@ -8,7 +8,7 @@ from dyson import util from dyson.representations.spectral import Spectral -from dyson.solvers import MBLGF +from dyson.solvers import Exact, MBLGF if TYPE_CHECKING: from pyscf import scf @@ -83,7 +83,7 @@ def test_vs_exact_solver_central( exact_p = exact_cache(mf, expression_method.p) assert exact_h.result is not None assert exact_p.result is not None - result_exact_ph = Spectral.combine(exact_h.result, exact_p.result) + result_exact_ph = Spectral.combine_for_greens_function(exact_h.result, exact_p.result) # Get the self-energy and Green's function from the exact solver static_exact = result_exact_ph.get_static_self_energy() @@ -99,7 +99,7 @@ def test_vs_exact_solver_central( mblgf_p.kernel() assert mblgf_h.result is not None assert mblgf_p.result is not None - result_ph = Spectral.combine(mblgf_h.result, mblgf_p.result) + result_ph = Spectral.combine_for_greens_function(mblgf_h.result, mblgf_p.result) assert helper.have_equal_moments( mblgf_h.result.get_self_energy(), exact_h.result.get_self_energy(), nmom_gf - 2 diff --git a/tests/test_mblse.py b/tests/test_mblse.py index f0f7468..085e80d 100644 --- a/tests/test_mblse.py +++ b/tests/test_mblse.py @@ -8,7 +8,7 @@ from dyson import util from dyson.representations.spectral import Spectral -from dyson.solvers import MBLSE +from dyson.solvers import Exact, MBLSE if TYPE_CHECKING: from pyscf import scf @@ -73,7 +73,8 @@ def test_vs_exact_solver_central( exact_p = exact_cache(mf, expression_method.p) assert exact_h.result is not None assert exact_p.result is not None - result_exact_ph = Spectral.combine(exact_h.result, exact_p.result) + overlap = expression_h.build_overlap() + expression_p.build_overlap() + result_exact_ph = Spectral.combine_for_greens_function(exact_h.result, exact_p.result, overlap=overlap) # Get the self-energy and Green's function from the exact solver static_exact = result_exact_ph.get_static_self_energy() @@ -101,15 +102,20 @@ def test_vs_exact_solver_central( hermitian=hermitian, ) result_p = mblse_p.kernel() - result_ph = Spectral.combine(result_h, result_p) - # Recover the self-energy and Green's function from the MBLSE solver + # Combine for Green's function + result_ph = Spectral.combine_for_greens_function(result_h, result_p, overlap=overlap) static = result_ph.get_static_self_energy() self_energy = result_ph.get_self_energy() greens_function = result_ph.get_greens_function() assert helper.are_equal_arrays(static, static_exact) - assert helper.have_equal_moments(self_energy, se_h_moments_exact + se_p_moments_exact, nmom_se) assert helper.have_equal_moments(self_energy, self_energy_exact, nmom_se) assert helper.recovers_greens_function(static, self_energy, greens_function, 4) assert helper.have_equal_moments(greens_function, greens_function_exact, nmom_se) + + # Combine for self-energy + result_ph = Spectral.combine_for_self_energy(result_h, result_p, overlap=overlap) + self_energy = result_ph.get_self_energy() + + assert helper.have_equal_moments(self_energy, se_h_moments_exact + se_p_moments_exact, nmom_se) From 55045f65422aab564911115f4f978cb1dfad2b0d Mon Sep 17 00:00:00 2001 From: Oliver Backhouse Date: Fri, 23 Jan 2026 17:59:31 +0000 Subject: [PATCH 07/17] linting --- dyson/representations/spectral.py | 24 +++++++++++++----------- dyson/solvers/static/mblse.py | 2 +- tests/test_davidson.py | 2 +- tests/test_density.py | 2 +- tests/test_downfolded.py | 2 +- tests/test_expressions.py | 2 +- tests/test_mblgf.py | 2 +- tests/test_mblse.py | 6 ++++-- 8 files changed, 23 insertions(+), 19 deletions(-) diff --git a/dyson/representations/spectral.py b/dyson/representations/spectral.py index 5d24e2c..c1fee77 100644 --- a/dyson/representations/spectral.py +++ b/dyson/representations/spectral.py @@ -252,8 +252,8 @@ def combine_for_greens_function( Args: *spectrals: Spectral representations to combine. - overlap: Overlap matrix for the physical space. If not passed, the overlap matrices from the - individual spectral representations will be summed. + overlap: Overlap matrix for the physical space. If not passed, the overlap matrices from + the individual spectral representations will be summed. chempot: Chemical potential to use for the combined result. If not provided, the chemical potential from the individual results will be used, raising an exception if they are inconsistent. @@ -274,8 +274,8 @@ def combine_for_self_energy( Args: *spectrals: Spectral representations to combine. - overlap: Overlap matrix for the physical space. If not passed, the overlap matrices from the - individual spectral representations will be summed. + overlap: Overlap matrix for the physical space. If not passed, the overlap matrices from + the individual spectral representations will be summed. chempot: Chemical potential to use for the combined result. If not provided, the chemical potential from the individual results will be used, raising an exception if they are inconsistent. @@ -366,8 +366,7 @@ def _get_chemical_potential(spectrals: tuple[Spectral, ...]) -> float | None: Raises: ValueError: If the chemical potentials are inconsistent. """ - chempot_set = {s.chempot for s in spectrals} - chempot_set.discard(None) + chempot_set = {s.chempot for s in spectrals if s.chempot is not None} if not all(np.isclose(c1, c2) for c1 in chempot_set for c2 in chempot_set): raise ValueError("Inconsistent chemical potentials in spectral representations.") return chempot_set.pop() if chempot_set else None @@ -401,7 +400,7 @@ def combine_for_greens_function( if chempot is None: chempot = _get_chemical_potential(spectrals) if overlap is None: - overlap = sum(spectral.overlap for spectral in spectrals) + overlap = sum((spectral.overlap for spectral in spectrals), np.zeros((nphys, nphys))) # Combine the eigenvalues and Dyson orbitals eigvals = np.concatenate([spectral.eigvals for spectral in spectrals], axis=0) @@ -431,8 +430,8 @@ def combine_for_greens_function( left, right = util.biorthonormalise(left, right) # Unorthogonalise the vectors - left[:, : nphys] = left[:, : nphys] @ unorth.T.conj() - right[:, : nphys] = right[:, : nphys] @ unorth + left[:, :nphys] = left[:, :nphys] @ unorth.T.conj() + right[:, :nphys] = right[:, :nphys] @ unorth # Get the combined result eigvecs = np.array([left.T.conj(), right.T.conj()]) if not hermitian else left.T.conj() @@ -463,13 +462,16 @@ def combine_for_self_energy( """ if len(spectrals) == 1: return spectrals[0] + nphys = _get_physical_space_size(spectrals) if chempot is None: chempot = _get_chemical_potential(spectrals) if overlap is None: - overlap = sum(spectral.overlap for spectral in spectrals) + overlap = sum((spectral.overlap for spectral in spectrals), np.zeros((nphys, nphys))) # Combine the self-energies - static = sum(spectral.get_static_self_energy() for spectral in spectrals) + static = sum( + (spectral.get_static_self_energy() for spectral in spectrals), np.zeros((nphys, nphys)) + ) self_energy = spectrals[0].get_self_energy(chempot=chempot).copy() for spectral in spectrals[1:]: self_energy = self_energy.concatenate(spectral.get_self_energy(chempot=chempot)) diff --git a/dyson/solvers/static/mblse.py b/dyson/solvers/static/mblse.py index 194a2ff..a1dadd4 100644 --- a/dyson/solvers/static/mblse.py +++ b/dyson/solvers/static/mblse.py @@ -398,7 +398,7 @@ def _recurrence_iteration_non_hermitian( return error_sqrt, error_inv_sqrt, error_moments - #FIXME: being called twice? + # FIXME: being called twice? def solve(self, iteration: int | None = None) -> Spectral: """Solve the eigenvalue problem at a given iteration. diff --git a/tests/test_davidson.py b/tests/test_davidson.py index b3cd0a0..725121e 100644 --- a/tests/test_davidson.py +++ b/tests/test_davidson.py @@ -9,7 +9,7 @@ from dyson.representations.lehmann import Lehmann from dyson.representations.spectral import Spectral -from dyson.solvers import Davidson, Exact +from dyson.solvers import Davidson if TYPE_CHECKING: from pyscf import scf diff --git a/tests/test_density.py b/tests/test_density.py index 1b1b805..9af7f4e 100644 --- a/tests/test_density.py +++ b/tests/test_density.py @@ -8,7 +8,7 @@ import pytest from dyson.representations.spectral import Spectral -from dyson.solvers import DensityRelaxation, Exact +from dyson.solvers import DensityRelaxation from dyson.solvers.static.density import get_fock_matrix_function if TYPE_CHECKING: diff --git a/tests/test_downfolded.py b/tests/test_downfolded.py index bc6c719..2998db0 100644 --- a/tests/test_downfolded.py +++ b/tests/test_downfolded.py @@ -8,7 +8,7 @@ import pytest from dyson.representations.spectral import Spectral -from dyson.solvers import Downfolded, Exact +from dyson.solvers import Downfolded if TYPE_CHECKING: from pyscf import scf diff --git a/tests/test_expressions.py b/tests/test_expressions.py index 67e8121..1ae996a 100644 --- a/tests/test_expressions.py +++ b/tests/test_expressions.py @@ -11,8 +11,8 @@ from dyson import util from dyson.expressions import ADC2, CCSD, FCI, HF, ADC2x -from dyson.solvers import Davidson, Exact from dyson.representations.spectral import Spectral +from dyson.solvers import Davidson, Exact if TYPE_CHECKING: from typing import Any diff --git a/tests/test_mblgf.py b/tests/test_mblgf.py index 4e60dbc..15b290a 100644 --- a/tests/test_mblgf.py +++ b/tests/test_mblgf.py @@ -8,7 +8,7 @@ from dyson import util from dyson.representations.spectral import Spectral -from dyson.solvers import Exact, MBLGF +from dyson.solvers import MBLGF if TYPE_CHECKING: from pyscf import scf diff --git a/tests/test_mblse.py b/tests/test_mblse.py index 085e80d..4eb6efe 100644 --- a/tests/test_mblse.py +++ b/tests/test_mblse.py @@ -8,7 +8,7 @@ from dyson import util from dyson.representations.spectral import Spectral -from dyson.solvers import Exact, MBLSE +from dyson.solvers import MBLSE if TYPE_CHECKING: from pyscf import scf @@ -74,7 +74,9 @@ def test_vs_exact_solver_central( assert exact_h.result is not None assert exact_p.result is not None overlap = expression_h.build_overlap() + expression_p.build_overlap() - result_exact_ph = Spectral.combine_for_greens_function(exact_h.result, exact_p.result, overlap=overlap) + result_exact_ph = Spectral.combine_for_greens_function( + exact_h.result, exact_p.result, overlap=overlap + ) # Get the self-energy and Green's function from the exact solver static_exact = result_exact_ph.get_static_self_energy() From 9a7d6b97b8a108b8fd3cf7161736040afad44e5a Mon Sep 17 00:00:00 2001 From: Oliver Backhouse Date: Tue, 27 Jan 2026 17:38:45 +0000 Subject: [PATCH 08/17] Fix examples for new combination functions --- examples/particle-hole-separation.py | 4 ++-- examples/spectra.py | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/examples/particle-hole-separation.py b/examples/particle-hole-separation.py index 41eef27..09220fe 100644 --- a/examples/particle-hole-separation.py +++ b/examples/particle-hole-separation.py @@ -23,10 +23,10 @@ solver_p = MBLGF.from_expression(exp_p, max_cycle=1) solver_p.kernel() -# Combine the results -- this function operators by projecting the result back into a self-energy +# Combine the results -- this function operates by projecting the result back into a self-energy # and combining the two self-energies, before diagonalising the combined self-energy to get a new # result spectrum. This may have unwanted consequences for some methodology, so use with care. -result = Spectral.combine(solver_h.result, solver_p.result) +result = Spectral.combine_for_greens_function(solver_h.result, solver_p.result) # Get the spectral functions grid = GridRF.from_uniform(-3.0, 3.0, 1024, eta=0.05) diff --git a/examples/spectra.py b/examples/spectra.py index c91c959..b244af5 100644 --- a/examples/spectra.py +++ b/examples/spectra.py @@ -24,7 +24,7 @@ exact_h.kernel() exact_p = Exact.from_expression(exp_p) exact_p.kernel() -result = exact_h.result.combine(exact_p.result) +result = exact_h.result.combine_for_greens_function(exact_p.result) static = result.get_static_self_energy() self_energy = result.get_self_energy() From 23c069439bcc17da524b560cbbdfcf89476c71c1 Mon Sep 17 00:00:00 2001 From: Oliver Backhouse Date: Tue, 27 Jan 2026 17:42:18 +0000 Subject: [PATCH 09/17] Fix for numpy promotion --- dyson/representations/lehmann.py | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/dyson/representations/lehmann.py b/dyson/representations/lehmann.py index 941173b..134b7fb 100644 --- a/dyson/representations/lehmann.py +++ b/dyson/representations/lehmann.py @@ -568,11 +568,13 @@ def matvec(self, physical: Array, vector: Array, chempot: bool | float = False) # Contract the supermatrix vector_phys, vector_aux = np.split(vector, [self.nphys]) - result_phys = util.einsum("pq,q...->p...", physical, vector_phys) - result_phys += util.einsum("pk,k...->p...", right, vector_aux) - result_aux = util.einsum("pk,p...->k...", left.conj(), vector_phys) - result_aux += util.einsum("k,k...->k...", energies, vector_aux) - result = np.concatenate((result_phys, result_aux), axis=0) + result_phys_phys = util.einsum("pq,q...->p...", physical, vector_phys) + result_phys_aux = util.einsum("pk,k...->p...", right, vector_aux) + result_aux_phys = util.einsum("pk,p...->k...", left.conj(), vector_phys) + result_aux_aux = util.einsum("k,k...->k...", energies, vector_aux) + result = np.concatenate( + (result_phys_phys + result_phys_aux, result_aux_phys + result_aux_aux), axis=0 + ) return result From 8071fa8980576f21093aa0fd82fe7156d913f9fd Mon Sep 17 00:00:00 2001 From: Oliver Backhouse Date: Tue, 27 Jan 2026 17:51:26 +0000 Subject: [PATCH 10/17] Remove GW from examples --- examples/solver-aufbau.py | 20 ++++++++++++-------- examples/solver-density.py | 20 ++++++++++++-------- examples/solver-shift.py | 20 ++++++++++++-------- 3 files changed, 36 insertions(+), 24 deletions(-) diff --git a/examples/solver-aufbau.py b/examples/solver-aufbau.py index b81896e..f403f0a 100644 --- a/examples/solver-aufbau.py +++ b/examples/solver-aufbau.py @@ -7,7 +7,7 @@ from pyscf import gto, scf -from dyson import MBLGF, TDAGW, AufbauPrinciple, Exact +from dyson import MBLGF, FCI, AufbauPrinciple, Exact, Spectral from dyson.solvers.static.chempot import search_aufbau_global # Get a molecule and mean-field from PySCF @@ -15,15 +15,19 @@ mf = scf.RHF(mol) mf.kernel() -# Use a TDA-GW Dyson expression for the Hamiltonian -exp = TDAGW.dyson.from_mf(mf) +# Use an FCI expression for the Hamiltonian +exp_h = FCI.h.from_mf(mf) +exp_p = FCI.p.from_mf(mf) # Use the exact solver to get the self-energy for demonstration purposes -exact = Exact.from_expression(exp) -exact.kernel() -static = exact.result.get_static_self_energy() -self_energy = exact.result.get_self_energy() -overlap = exact.result.get_overlap() +exact_h = Exact.from_expression(exp_h) +exact_h.kernel() +exact_p = Exact.from_expression(exp_p) +exact_p.kernel() +result = Spectral.combine_for_self_energy(exact_h.result, exact_p.result) +static = result.get_static_self_energy() +self_energy = result.get_self_energy() +overlap = result.get_overlap() # Solve the Hamiltonian using the Aufbau solver, initialisation via either: diff --git a/examples/solver-density.py b/examples/solver-density.py index f8a009e..6eb4112 100644 --- a/examples/solver-density.py +++ b/examples/solver-density.py @@ -15,7 +15,7 @@ from pyscf import gto, scf -from dyson import MBLSE, TDAGW, AufbauPrinciple, AuxiliaryShift, DensityRelaxation, Exact +from dyson import MBLSE, FCI, AufbauPrinciple, AuxiliaryShift, DensityRelaxation, Exact, Spectral from dyson.solvers.static.density import get_fock_matrix_function # Get a molecule and mean-field from PySCF @@ -23,15 +23,19 @@ mf = scf.RHF(mol) mf.kernel() -# Use a TDA-GW Dyson expression for the Hamiltonian -exp = TDAGW.dyson.from_mf(mf) +# Use an FCI expression for the Hamiltonian +exp_h = FCI.h.from_mf(mf) +exp_p = FCI.p.from_mf(mf) # Use the exact solver to get the self-energy for demonstration purposes -exact = Exact.from_expression(exp) -exact.kernel() -static = exact.result.get_static_self_energy() -self_energy = exact.result.get_self_energy() -overlap = exact.result.get_overlap() +exact_h = Exact.from_expression(exp_h) +exact_h.kernel() +exact_p = Exact.from_expression(exp_p) +exact_p.kernel() +result = Spectral.combine_for_self_energy(exact_h.result, exact_p.result) +static = result.get_static_self_energy() +self_energy = result.get_self_energy() +overlap = result.get_overlap() # Solve the Hamiltonian using the density relaxation solver, initialisation via either: diff --git a/examples/solver-shift.py b/examples/solver-shift.py index 0a93525..8b24fac 100644 --- a/examples/solver-shift.py +++ b/examples/solver-shift.py @@ -8,22 +8,26 @@ from pyscf import gto, scf -from dyson import MBLSE, TDAGW, AufbauPrinciple, AuxiliaryShift, Exact +from dyson import MBLSE, FCI, AufbauPrinciple, AuxiliaryShift, Exact, Spectral # Get a molecule and mean-field from PySCF mol = gto.M(atom="Li 0 0 0; H 0 0 1.64", basis="sto-3g", verbose=0) mf = scf.RHF(mol) mf.kernel() -# Use a TDA-GW Dyson expression for the Hamiltonian -exp = TDAGW.dyson.from_mf(mf) +# Use an FCI expression for the Hamiltonian +exp_h = FCI.h.from_mf(mf) +exp_p = FCI.p.from_mf(mf) # Use the exact solver to get the self-energy for demonstration purposes -exact = Exact.from_expression(exp) -exact.kernel() -static = exact.result.get_static_self_energy() -self_energy = exact.result.get_self_energy() -overlap = exact.result.get_overlap() +exact_h = Exact.from_expression(exp_h) +exact_h.kernel() +exact_p = Exact.from_expression(exp_p) +exact_p.kernel() +result = Spectral.combine_for_self_energy(exact_h.result, exact_p.result) +static = result.get_static_self_energy() +self_energy = result.get_self_energy() +overlap = result.get_overlap() # Solve the Hamiltonian using the auxiliary shift solver, initialisation via either: From 3264da948d337bb31cea34cba60a50b1e9c47f55 Mon Sep 17 00:00:00 2001 From: Oliver Backhouse Date: Tue, 27 Jan 2026 17:53:26 +0000 Subject: [PATCH 11/17] linting --- examples/solver-aufbau.py | 2 +- examples/solver-density.py | 2 +- examples/solver-shift.py | 2 +- 3 files changed, 3 insertions(+), 3 deletions(-) diff --git a/examples/solver-aufbau.py b/examples/solver-aufbau.py index f403f0a..d052c0c 100644 --- a/examples/solver-aufbau.py +++ b/examples/solver-aufbau.py @@ -7,7 +7,7 @@ from pyscf import gto, scf -from dyson import MBLGF, FCI, AufbauPrinciple, Exact, Spectral +from dyson import FCI, MBLGF, AufbauPrinciple, Exact, Spectral from dyson.solvers.static.chempot import search_aufbau_global # Get a molecule and mean-field from PySCF diff --git a/examples/solver-density.py b/examples/solver-density.py index 6eb4112..1dd7fe8 100644 --- a/examples/solver-density.py +++ b/examples/solver-density.py @@ -15,7 +15,7 @@ from pyscf import gto, scf -from dyson import MBLSE, FCI, AufbauPrinciple, AuxiliaryShift, DensityRelaxation, Exact, Spectral +from dyson import FCI, MBLSE, AufbauPrinciple, AuxiliaryShift, DensityRelaxation, Exact, Spectral from dyson.solvers.static.density import get_fock_matrix_function # Get a molecule and mean-field from PySCF diff --git a/examples/solver-shift.py b/examples/solver-shift.py index 8b24fac..cc814ed 100644 --- a/examples/solver-shift.py +++ b/examples/solver-shift.py @@ -8,7 +8,7 @@ from pyscf import gto, scf -from dyson import MBLSE, FCI, AufbauPrinciple, AuxiliaryShift, Exact, Spectral +from dyson import FCI, MBLSE, AufbauPrinciple, AuxiliaryShift, Exact, Spectral # Get a molecule and mean-field from PySCF mol = gto.M(atom="Li 0 0 0; H 0 0 1.64", basis="sto-3g", verbose=0) From cecca1cdb8ac6af2932d516f3b6ba3cdd3f31f56 Mon Sep 17 00:00:00 2001 From: Oliver Backhouse Date: Tue, 27 Jan 2026 18:03:12 +0000 Subject: [PATCH 12/17] Swap test case because He cc-pvdz sucks --- tests/conftest.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/tests/conftest.py b/tests/conftest.py index 68f9dbb..4cafe67 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -33,9 +33,9 @@ basis="sto-3g", verbose=0, ), - "he-ccpvdz": gto.M( - atom="He 0 0 0", - basis="cc-pvdz", + "lih-sto3g": gto.M( + atom="Li 0 0 0; H 0 0 1.64", + basis="sto-3g", verbose=0, ), } @@ -43,7 +43,7 @@ MF_CACHE = { "h2-631g": scf.RHF(MOL_CACHE["h2-631g"]).run(conv_tol=1e-12), "h2o-sto3g": scf.RHF(MOL_CACHE["h2o-sto3g"]).run(conv_tol=1e-12), - "he-ccpvdz": scf.RHF(MOL_CACHE["he-ccpvdz"]).run(conv_tol=1e-12), + "lih-sto3g": scf.RHF(MOL_CACHE["lih-sto3g"]).run(conv_tol=1e-12), } for key, mf in MF_CACHE.items(): From b623df93d90326d7a4bbb0d183c4ba1462aad4f3 Mon Sep 17 00:00:00 2001 From: Oliver Backhouse Date: Tue, 27 Jan 2026 18:19:08 +0000 Subject: [PATCH 13/17] Skip numerically unstable tests --- tests/test_mblgf.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/tests/test_mblgf.py b/tests/test_mblgf.py index 15b290a..ece76d9 100644 --- a/tests/test_mblgf.py +++ b/tests/test_mblgf.py @@ -71,7 +71,8 @@ def test_vs_exact_solver_central( if expression_h.nconfig > 1024 or expression_p.nconfig > 1024: pytest.skip("Skipping test for large Hamiltonian") if request.node.name in ( - "test_vs_exact_solver_central[lih-631g-CCSD-3]", + "test_vs_exact_solver_central[lih-sto3-CCSD-2]", + "test_vs_exact_solver_central[lih-sto3-CCSD-3]", "test_vs_exact_solver_central[h2o-sto3g-CCSD-2]", "test_vs_exact_solver_central[h2o-sto3g-CCSD-3]", ): From b52e6eeb673872b59649ab896ec44ac711efc69a Mon Sep 17 00:00:00 2001 From: Oliver Backhouse Date: Tue, 27 Jan 2026 18:23:21 +0000 Subject: [PATCH 14/17] Fix typo --- tests/test_mblgf.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/test_mblgf.py b/tests/test_mblgf.py index ece76d9..ab41f38 100644 --- a/tests/test_mblgf.py +++ b/tests/test_mblgf.py @@ -71,8 +71,8 @@ def test_vs_exact_solver_central( if expression_h.nconfig > 1024 or expression_p.nconfig > 1024: pytest.skip("Skipping test for large Hamiltonian") if request.node.name in ( - "test_vs_exact_solver_central[lih-sto3-CCSD-2]", - "test_vs_exact_solver_central[lih-sto3-CCSD-3]", + "test_vs_exact_solver_central[lih-sto3g-CCSD-2]", + "test_vs_exact_solver_central[lih-sto3g-CCSD-3]", "test_vs_exact_solver_central[h2o-sto3g-CCSD-2]", "test_vs_exact_solver_central[h2o-sto3g-CCSD-3]", ): From 40ba85f37added2e173a6d87c0fd03dc7e5eda32 Mon Sep 17 00:00:00 2001 From: Oliver Backhouse Date: Tue, 27 Jan 2026 18:36:27 +0000 Subject: [PATCH 15/17] Update tests/test_davidson.py Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- tests/test_davidson.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/test_davidson.py b/tests/test_davidson.py index 725121e..f9f6430 100644 --- a/tests/test_davidson.py +++ b/tests/test_davidson.py @@ -154,7 +154,7 @@ def test_vs_exact_solver_central( exact_h.result, exact_p.result, overlap=overlap ) result_davidson = Spectral.combine_for_greens_function( - exact_h.result, exact_p.result, overlap=overlap + davidson_h.result, davidson_p.result, overlap=overlap ) # Get the self-energy and Green's function from the Davidson solver From 2fe893a398670a5b72ff404acb4ca90c126a3a87 Mon Sep 17 00:00:00 2001 From: Oliver Backhouse Date: Tue, 27 Jan 2026 18:35:19 +0000 Subject: [PATCH 16/17] Don't return matrix power error --- dyson/util/linalg.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/dyson/util/linalg.py b/dyson/util/linalg.py index f9a214f..24d80f7 100644 --- a/dyson/util/linalg.py +++ b/dyson/util/linalg.py @@ -340,8 +340,8 @@ def project_eigenvectors( # If the system is not hermitian, we need to ensure biorthonormality overlap = ket.conj() @ bra.T - orth, orth_error = matrix_power(overlap, -0.5, hermitian=hermitian, return_error=True) - unorth, unorth_error = matrix_power(overlap, 0.5, hermitian=hermitian, return_error=True) + orth, _ = matrix_power(overlap, -0.5, hermitian=hermitian, return_error=False) + unorth, _ = matrix_power(overlap, 0.5, hermitian=hermitian, return_error=False) # Work in an orthonormal physical basis bra = bra.T @ orth From 6a717da0221c3f8402a0ba0c308a980f61453263 Mon Sep 17 00:00:00 2001 From: Oliver Backhouse Date: Tue, 27 Jan 2026 18:37:45 +0000 Subject: [PATCH 17/17] Fix use of classmethod --- examples/spectra.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/examples/spectra.py b/examples/spectra.py index b244af5..da2f554 100644 --- a/examples/spectra.py +++ b/examples/spectra.py @@ -7,6 +7,7 @@ from dyson.expressions import ADC2 from dyson.grids import GridRF from dyson.plotting import format_axes_spectral_function, plot_dynamic +from dyson.representations import Spectral from dyson.solvers import CPGF, MBLGF, MBLSE, CorrectionVector, Downfolded, Exact # Get a molecule and mean-field from PySCF @@ -24,7 +25,7 @@ exact_h.kernel() exact_p = Exact.from_expression(exp_p) exact_p.kernel() -result = exact_h.result.combine_for_greens_function(exact_p.result) +result = Spectral.combine_for_greens_function(exact_h.result, exact_p.result) static = result.get_static_self_energy() self_energy = result.get_self_energy()