21367620794

Committed 26 Jan 2026 05:37PM UTC coverage: 74.549% (+0.3%) from 74.293%

Build # 21367620794

Build Type

Pull #200

github

Committed by

web-flow

Commit Message

Merge 92d6b1ccc into 9ff4abd52

Pull Request Pull Request #200: Generalize OneParticleOperator to NParticleOperator and separate operators from densities

Run Details

1257 of 1986 branches covered (63.29%)

Branch coverage included in aggregate %.

559 of 630 new or added lines in 23 files covered. (88.73%)

11 existing lines in 1 file now uncovered.

7785 of 10143 relevant lines covered (76.75%)

168646.41 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

90.96

/adcc/NParticleOperator.py

#!/usr/bin/env python3
## vi: tabstop=4 shiftwidth=4 softtabstop=4 expandtab
## ---------------------------------------------------------------------
##
## Copyright (C) 2026 by the adcc authors
##
## This file is part of adcc.
##
## adcc is free software: you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published
## by the Free Software Foundation, either version 3 of the License, or
## (at your option) any later version.
##
## adcc is distributed in the hope that it will be useful,
## but WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
## GNU General Public License for more details.
##
## You should have received a copy of the GNU General Public License
## along with adcc. If not, see <http://www.gnu.org/licenses/>.
##
## ---------------------------------------------------------------------
from collections import Counter
from dataclasses import dataclass
from enum import Enum
from itertools import product
import math
import numpy as np

import libadcc

from .functions import evaluate
from .MoSpaces import split_spaces
from .Tensor import Tensor


class OperatorSymmetry(Enum):
    NOSYMMETRY = 0
    HERMITIAN = 1
    ANTIHERMITIAN = 2

    def to_str(self) -> str:
        return self.name.lower()


@dataclass(frozen=True)
class BlockInfo:
    canonical: str
    factor: int
    transpose: tuple[int, ...]


class NParticleOperator:
    def __init__(self, spaces, n_particle_op, symmetry=OperatorSymmetry.HERMITIAN):
        """
        Construct an NParticleOperator object. All blocks are initialised
        as zero blocks.

        Parameters
        ----------
        spaces : adcc.MoSpaces or adcc.ReferenceState or adcc.LazyMp
            MoSpaces object

        n_particle_operator: int
            Particle rank of the operator, i.e. the number of creation and
            annihilation operators involved.

        symmetry : OperatorSymmetry, optional
            Symmetry type of the operator. Can be:
            - `OperatorSymmetry.NOSYMMETRY` : No symmetry is enforced
            - `OperatorSymmetry.HERMITIAN` : Operator is Hermitian (O^\\dagger = O)
            - `OperatorSymmetry.ANTIHERMITIAN` : Operator is Antihermitian
                                                 (O^\\dagger = -O)
            Default is `OperatorSymmetry.HERMITIAN`.
        """
        self._n_particle_op = n_particle_op
        if hasattr(spaces, "mospaces"):
            self.mospaces = spaces.mospaces
        else:
            self.mospaces = spaces
        self.symmetry = symmetry

        # Set reference_state if possible
        if isinstance(spaces, libadcc.ReferenceState):
            self.reference_state = spaces
        elif hasattr(spaces, "reference_state"):
            assert isinstance(spaces.reference_state, libadcc.ReferenceState)
            self.reference_state = spaces.reference_state

        occs = sorted(self.mospaces.subspaces_occupied)
        virts = sorted(self.mospaces.subspaces_virtual)
        self.orbital_subspaces = occs + virts

        # check that orbital subspaces are correct
        assert sum(self.mospaces.n_orbs(ss) for ss in self.orbital_subspaces) \
            == self.mospaces.n_orbs("f")
        self._tensors = {}

        # Initialize all blocks; symmetry rules are applied lazily upon access.
        combs = list(product(self.orbital_subspaces,
                             repeat=2 * self._n_particle_op))
        self.blocks = ["".join((comb)) for comb in combs]
        # Mapping of block onto canonical blcok
        self._block_info = {}           # block -> BlockInfo
        self._canonical_factors = Counter()  # canonical_block -> counter (factor)

        for block in self.blocks:
            from .block import get_canonical_block
            bra = block[:2 * self.n_particle_op]
            ket = block[2 * self.n_particle_op:]
            canonical_block, factor, transpose = get_canonical_block(
                bra, ket, self.symmetry
            )
            self._block_info[block] = BlockInfo(
                canonical=canonical_block,
                factor=factor,
                transpose=transpose,
            )
            self._canonical_factors[canonical_block] += 1

    @property
    def canonical_factors(self) -> dict[str, int]:
        """
        Returns canonical block multiplicity factors.
        """
        return dict(self._canonical_factors)

    @property
    def canonical_blocks(self) -> tuple[str, ...]:
        """
        Returns tuple of canonical block labels.
        """
        return tuple(self._canonical_factors.keys())

    @property
    def n_particle_op(self) -> int:
        """
        Returns the particle rank of the operator.
        """
        return self._n_particle_op

    @property
    def shape(self) -> tuple[int, ...]:
        """
        Returns the shape tuple of the NParticleOperator
        """
        size = self.mospaces.n_orbs("f")
        return (size, size,) * self.n_particle_op

    @property
    def size(self) -> int:
        """
        Returns the number of elements of the NParticleOperator
        """
        return np.prod(self.shape)

    @property
    def blocks_nonzero(self) -> list[str, ...]:
        """
        Returns a list of the non-zero block labels
        """
        return [b for b in self.canonical_blocks if b in self._tensors]

    def is_zero_block(self, block):
        """
        Checks if block is explicitly marked as zero block.
        Returns False if the block does not exist.
        """
        block_info = self._block_info.get(block, None)
        if block_info is None:
            return False
        return block_info.canonical not in self._tensors

    def block(self, block):
        """
        Returns tensor of the given block.
        Does not create a block in case it is marked as a zero block.
        Use __getitem__ for that purpose.
        """
        block_info = self._block_info.get(block)
        if block_info is None:
            raise KeyError(f"Invalid block {block} requested. "
                           f"Available blocks are: {self.blocks}.")

        c = block_info.canonical
        if c not in self._tensors:
            raise KeyError("The block function does not support "
                           "access to zero-blocks. Available non-zero "
                           f"blocks are: {self.blocks_nonzero}.")

        if block == c:
            return self._tensors[c]
        return block_info.factor * self._tensors[c].transpose(block_info.transpose)

    def __getitem__(self, blk):
        if blk not in self.blocks:
            raise KeyError(f"Invalid block {blk} requested. "
                           f"Available blocks are: {self.blocks}.")

        block_info = self._block_info[blk]
        c = block_info.canonical

        if c not in self._tensors:
            sym = libadcc.make_symmetry_operator(
                self.mospaces, c, self.symmetry.to_str(), "1"
            )
            self._tensors[c] = Tensor(sym)

        if blk == c:
            return self._tensors[c]
        return block_info.factor * self._tensors[c].transpose(block_info.transpose)

    def __getattr__(self, attr) -> Tensor:
        from . import block as b
        return self.__getitem__(b.__getattr__(attr))

    def __setitem__(self, block, tensor):
        block_info = self._block_info.get(block)
        if block_info is None or block_info.canonical != block:
            raise KeyError(f"Invalid block {block} assigned. "
                           f"Available blocks are: {self.blocks}.")

        spaces = split_spaces(block)
        expected_shape = tuple(self.mospaces.n_orbs(space) for space in spaces)
        if expected_shape != tensor.shape:
            raise ValueError("Invalid shape of incoming tensor. "
                             f"Expected shape {expected_shape}, but "
                             f"got shape {tensor.shape} instead.")

        self._tensors[block] = tensor

    def __setattr__(self, attr, value):
        try:
            from . import block as b
            self.__setitem__(b.__getattr__(attr), value)
        except AttributeError:
            super().__setattr__(attr, value)

    def set_zero_block(self, block):
        """
        Set a given block as zero block
        """
        if block not in self.canonical_blocks:
            raise KeyError(f"Invalid block {block} set as zero block. "
                           f"Available blocks are: {self.blocks}.")
        self._tensors.pop(block)

    def to_ndarray(self) -> np.ndarray:
        """
        Returns the NParticleOperator as a contiguous
        np.ndarray instance including all blocks
        """
        # offsets to start index of spaces
        offsets = {
            sp: sum(
                self.mospaces.n_orbs(ss)
                for ss in self.orbital_subspaces[:self.orbital_subspaces.index(sp)]
            )
            for sp in self.orbital_subspaces
        }
        # slices for each space
        slices = {
            sp: slice(offsets[sp], offsets[sp] + self.mospaces.n_orbs(sp))
            for sp in self.orbital_subspaces
        }
        # build complete np.ndarray -> go through all blocks
        ret = np.zeros((self.shape))
        for block in self.blocks:
            if self.is_zero_block(block):
                continue
            spaces = split_spaces(block)
            slice_list = tuple(slices[sp] for sp in spaces)
            dm_block = self[block].to_ndarray()
            ret[slice_list] = dm_block
        return ret

    def _construct_empty(self):
        """
        Create an empty instance of an NParticleOperator
        """
        return self.__class__(
            self.mospaces,
            self.n_particle_op,
            symmetry=self.symmetry,
        )

    def copy(self):
        """
        Return a deep copy of the NParticleOperator
        """
        ret = self._construct_empty()
        for b in self.blocks_nonzero:
            ret[b] = self.block(b).copy()
        if hasattr(self, "reference_state"):
            ret.reference_state = self.reference_state
        return ret

    def _transform_to_ao(self, refstate):
        raise NotImplementedError("Needs to be implemented on the "
                                  f"{self.__class__.__name__} class.")

    def to_ao_basis(self, refstate_or_coefficients=None):
        """
        Transforms the NParticleOperator to the atomic orbital
        basis using a ReferenceState or a coefficient map. If no
        ReferenceState or coefficient map is given, the ReferenceState
        used to construct the NParticleOperator is taken instead.
        """
        if isinstance(refstate_or_coefficients, (dict, libadcc.ReferenceState)):
            return self._transform_to_ao(refstate_or_coefficients)
        elif refstate_or_coefficients is None:
            if not hasattr(self, "reference_state"):
                raise ValueError("Argument reference_state is required if no "
                                 "reference_state is stored in the "
                                 "NParticleOperator")
            return self._transform_to_ao(self.reference_state)
        else:
            raise TypeError("Argument type not supported.")

    def __iadd__(self, other):
        from . import block as b
        if self.mospaces != other.mospaces:
            raise ValueError(f"Cannot add {self.__class__.__name__}s with "
                             "differing mospaces.")
        if self.symmetry is not OperatorSymmetry.NOSYMMETRY and \
                self.symmetry is not other.symmetry:
            raise ValueError("Cannot add non-symmetric matrix "
                             "in-place to symmetric one.")

        for block in other.blocks_nonzero:
            if self.is_zero_block(block):
                self[block] = other.block(block).copy()
            else:
                self[block] = self.block(block) + other.block(block)

        if self.symmetry is OperatorSymmetry.NOSYMMETRY \
                and other.symmetry is not OperatorSymmetry.NOSYMMETRY:
            for block in self.blocks_nonzero:
                bra = block[:2 * self.n_particle_op]
                ket = block[2 * self.n_particle_op:]
                c_block, factor, transpose = b.get_canonical_block(
                    bra, ket, other.symmetry
                )
                if block == c_block:
                    continue  # Done already
                obT = 0
                if not other.is_zero_block(c_block):
                    obT = factor * other.block(c_block).transpose(transpose)
                if not self.is_zero_block(block):
                    obT += self.block(block)
                self[block] = evaluate(obT)

        # Update ReferenceState pointer
        if hasattr(self, "reference_state"):
            if hasattr(other, "reference_state") \
                    and self.reference_state != other.reference_state:
                delattr(self, "reference_state")
        return self

    def __isub__(self, other):
        from . import block as b
        if self.mospaces != other.mospaces:
            raise ValueError(f"Cannot add {self.__class__.__name__}s with "
                             "differing mospaces.")
        if self.symmetry is not OperatorSymmetry.NOSYMMETRY and \
                self.symmetry is not other.symmetry:
            raise ValueError("Cannot add non-symmetric matrix "
                             "in-place to symmetric one.")

        for block in other.blocks_nonzero:
            if self.is_zero_block(block):
                self[block] = -1.0 * other.block(block)  # The copy is implicit
            else:
                self[block] = self.block(block) - other.block(block)

        if self.symmetry is OperatorSymmetry.NOSYMMETRY \
                and other.symmetry is not OperatorSymmetry.NOSYMMETRY:
            for block in self.blocks_nonzero:
                bra = block[:2 * self.n_particle_op]
                ket = block[2 * self.n_particle_op:]
                c_block, factor, transpose = b.get_canonical_block(
                    bra, ket, other.symmetry
                )
                if block == c_block:
                    continue  # Done already
                obT = 0
                if not other.is_zero_block(c_block):
                    obT = -1 * factor * other.block(c_block).transpose(transpose)
                if not self.is_zero_block(block):
                    obT += self.block(block)
                self[block] = evaluate(obT)

        # Update ReferenceState pointer
        if hasattr(self, "reference_state"):
            if hasattr(other, "reference_state") \
                    and self.reference_state != other.reference_state:
                delattr(self, "reference_state")
        return self

    def __imul__(self, other):
        if not isinstance(other, (float, int)):
            return NotImplemented
        for b in self.blocks_nonzero:
            self[b] = self.block(b) * other
        return self

    def __add__(self, other):
        if (
            self.symmetry is not other.symmetry
            and self.symmetry is not OperatorSymmetry.NOSYMMETRY
            and other.symmetry is not OperatorSymmetry.NOSYMMETRY
        ):
            raise ValueError("Addition of Hermitian and Antihermitian "
                             "operators is not implemented.")

        if self.symmetry is OperatorSymmetry.NOSYMMETRY \
                or other.symmetry is not OperatorSymmetry.NOSYMMETRY:
            return self.copy().__iadd__(other)
        else:
            return other.copy().__iadd__(self)

    def __sub__(self, other):
        if (
            self.symmetry is not other.symmetry
            and self.symmetry is not OperatorSymmetry.NOSYMMETRY
            and other.symmetry is not OperatorSymmetry.NOSYMMETRY
        ):
            raise ValueError("Substraction of Hermitian and Antihermitian "
                             "operators is not implemented.")

        if self.symmetry is OperatorSymmetry.NOSYMMETRY \
                or other.symmetry is not OperatorSymmetry.NOSYMMETRY:
            return self.copy().__isub__(other)
        else:
            return (-1.0 * other).__iadd__(self)

    def __mul__(self, other):
        return self.copy().__imul__(other)

    def __rmul__(self, other):
        return self.copy().__imul__(other)

    def evaluate(self):
        for b in self.blocks_nonzero:
            self.block(b).evaluate()
        return self

    def set_random(self):
        for b in self.canonical_blocks:
            self[b].set_random()
        return self


def product_trace(op1, op2) -> float:
    """
    Compute the expectation value (inner product) of two N-particle operators.

    For the special case of a density operator and a general operator, the
    inner product is identical in the AO and MO basis.

    Parameters
    ----------
    op1, op2 : NParticleOperator or NParticleDensity
        Operators or densities to compute the inner product of.

    Returns
    -------
    float
        The trace / expectation value <op1, op2>.
    """
    # TODO use blocks_nonzero and build the set intersection
    #      to avoid the is_zero_block( ) checks below.
    #      I'm a bit hesitant to do this right now, because I'm lacking
    #      the time at the moment to build a more sophisticated test,
    #      which could potentially catch an arising error.
    if op1.n_particle_op != op2.n_particle_op:
        raise TypeError(
            "Both operators must have the same number of N-particle operators "
            f"(got {op1.n_particle_op} and {op2.n_particle_op})."
        )
    assert op1.blocks == op2.blocks
    if op1.symmetry == OperatorSymmetry.NOSYMMETRY:
        factors = op1.canonical_factors.copy()
    elif op2.symmetry == OperatorSymmetry.NOSYMMETRY:
        factors = op2.canonical_factors.copy()
    else:
        assert op1.canonical_factors == op2.canonical_factors
        factors = op1.canonical_factors.copy()
        # Off-diagonal terms can be removed since a Hermitian times
        # an anti-Hermitian operator yields for off diagonal terms
        # (e.g. for a one-particle operator) ov * ov and vo * vo = (ov).T * (−ov.T),
        # which cancel exactly to zero.
        if op1.symmetry is not op2.symmetry:
            to_remove = []
            for b in list(factors.keys()):
                spaces = split_spaces(b)
                n = op1.n_particle_op
                if spaces[:n] != spaces[n:]:
                    to_remove.append(b)
            for b in to_remove:
                factors.pop(b)
    ret = 0
    for b, factor in factors.items():
        if op1.is_zero_block(b) or op2.is_zero_block(b):
            continue

        ret += (1 / math.factorial(op1.n_particle_op) ** 2 * factor
                * op1.block(b).dot(op2.block(b)))
    return ret

1	#!/usr/bin/env python3
2	## vi: tabstop=4 shiftwidth=4 softtabstop=4 expandtab
3	## ---------------------------------------------------------------------
4	##
5	## Copyright (C) 2026 by the adcc authors
6	##
7	## This file is part of adcc.
8	##
9	## adcc is free software: you can redistribute it and/or modify
10	## it under the terms of the GNU General Public License as published
11	## by the Free Software Foundation, either version 3 of the License, or
12	## (at your option) any later version.
13	##
14	## adcc is distributed in the hope that it will be useful,
15	## but WITHOUT ANY WARRANTY; without even the implied warranty of
16	## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17	## GNU General Public License for more details.
18	##
19	## You should have received a copy of the GNU General Public License
20	## along with adcc. If not, see <http://www.gnu.org/licenses/>.
21	##
22	## ---------------------------------------------------------------------
23	from collections import Counter	2✔
24	from dataclasses import dataclass	2✔
25	from enum import Enum	2✔
26	from itertools import product	2✔
27	import math	2✔
28	import numpy as np	2✔
29
30	import libadcc	2✔
31
32	from .functions import evaluate	2✔
33	from .MoSpaces import split_spaces	2✔
34	from .Tensor import Tensor	2✔
35
36
37	class OperatorSymmetry(Enum):	2✔
38	NOSYMMETRY = 0	2✔
39	HERMITIAN = 1	2✔
40	ANTIHERMITIAN = 2	2✔
41
42	def to_str(self) -> str:	2✔
43	return self.name.lower()	2✔
44
45
46	@dataclass(frozen=True)	2✔
47	class BlockInfo:	2✔
48	canonical: str	2✔
49	factor: int	2✔
50	transpose: tuple[int, ...]	2✔
51
52
53	class NParticleOperator:	2✔
54	def __init__(self, spaces, n_particle_op, symmetry=OperatorSymmetry.HERMITIAN):	2✔
55	"""
56	Construct an NParticleOperator object. All blocks are initialised
57	as zero blocks.
58
59	Parameters
60	----------
61	spaces : adcc.MoSpaces or adcc.ReferenceState or adcc.LazyMp
62	MoSpaces object
63
64	n_particle_operator: int
65	Particle rank of the operator, i.e. the number of creation and
66	annihilation operators involved.
67
68	symmetry : OperatorSymmetry, optional
69	Symmetry type of the operator. Can be:
70	- `OperatorSymmetry.NOSYMMETRY` : No symmetry is enforced
71	- `OperatorSymmetry.HERMITIAN` : Operator is Hermitian (O^\\dagger = O)
72	- `OperatorSymmetry.ANTIHERMITIAN` : Operator is Antihermitian
73	(O^\\dagger = -O)
74	Default is `OperatorSymmetry.HERMITIAN`.
75	"""
76	self._n_particle_op = n_particle_op	2✔
77	if hasattr(spaces, "mospaces"):	2✔
78	self.mospaces = spaces.mospaces	2✔
79	else:
80	self.mospaces = spaces	2✔
81	self.symmetry = symmetry	2✔
82
83	# Set reference_state if possible
84	if isinstance(spaces, libadcc.ReferenceState):	2✔
85	self.reference_state = spaces	2✔
86	elif hasattr(spaces, "reference_state"):	2✔
87	assert isinstance(spaces.reference_state, libadcc.ReferenceState)	2✔
88	self.reference_state = spaces.reference_state	2✔
89
90	occs = sorted(self.mospaces.subspaces_occupied)	2✔
91	virts = sorted(self.mospaces.subspaces_virtual)	2✔
92	self.orbital_subspaces = occs + virts	2✔
93
94	# check that orbital subspaces are correct
95	assert sum(self.mospaces.n_orbs(ss) for ss in self.orbital_subspaces) \	2✔
96	== self.mospaces.n_orbs("f")
97	self._tensors = {}	2✔
98
99	# Initialize all blocks; symmetry rules are applied lazily upon access.
100	combs = list(product(self.orbital_subspaces,	2✔
101	repeat=2 * self._n_particle_op))
102	self.blocks = ["".join((comb)) for comb in combs]	2✔
103	# Mapping of block onto canonical blcok
104	self._block_info = {} # block -> BlockInfo	2✔
105	self._canonical_factors = Counter() # canonical_block -> counter (factor)	2✔
106
107	for block in self.blocks:	2✔
108	from .block import get_canonical_block	2✔
109	bra = block[:2 * self.n_particle_op]	2✔
110	ket = block[2 * self.n_particle_op:]	2✔
111	canonical_block, factor, transpose = get_canonical_block(	2✔
112	bra, ket, self.symmetry
113	)
114	self._block_info[block] = BlockInfo(	2✔
115	canonical=canonical_block,
116	factor=factor,
117	transpose=transpose,
118	)
119	self._canonical_factors[canonical_block] += 1	2✔
120
121	@property	2✔
122	def canonical_factors(self) -> dict[str, int]:	2✔
123	"""
124	Returns canonical block multiplicity factors.
125	"""
126	return dict(self._canonical_factors)	2✔
127
128	@property	2✔
129	def canonical_blocks(self) -> tuple[str, ...]:	2✔
130	"""
131	Returns tuple of canonical block labels.
132	"""
133	return tuple(self._canonical_factors.keys())	2✔
134
135	@property	2✔
136	def n_particle_op(self) -> int:	2✔
137	"""
138	Returns the particle rank of the operator.
139	"""
140	return self._n_particle_op	2✔
141
142	@property	2✔
143	def shape(self) -> tuple[int, ...]:	2✔
144	"""
145	Returns the shape tuple of the NParticleOperator
146	"""
147	size = self.mospaces.n_orbs("f")	2✔
148	return (size, size,) * self.n_particle_op	2✔
149
150	@property	2✔
151	def size(self) -> int:	2✔
152	"""
153	Returns the number of elements of the NParticleOperator
154	"""
155	return np.prod(self.shape)	2✔
156
157	@property	2✔
158	def blocks_nonzero(self) -> list[str, ...]:	2✔
159	"""
160	Returns a list of the non-zero block labels
161	"""
162	return [b for b in self.canonical_blocks if b in self._tensors]	2✔
163
164	def is_zero_block(self, block):	2✔
165	"""
166	Checks if block is explicitly marked as zero block.
167	Returns False if the block does not exist.
168	"""
169	block_info = self._block_info.get(block, None)	2✔
170	if block_info is None:	2!
NEW 171	return False	×
172	return block_info.canonical not in self._tensors	2✔
173
174	def block(self, block):	2✔
175	"""
176	Returns tensor of the given block.
177	Does not create a block in case it is marked as a zero block.
178	Use __getitem__ for that purpose.
179	"""
180	block_info = self._block_info.get(block)	2✔
181	if block_info is None:	2!
NEW 182	raise KeyError(f"Invalid block {block} requested. "	×
183	f"Available blocks are: {self.blocks}.")
184
185	c = block_info.canonical	2✔
186	if c not in self._tensors:	2✔
187	raise KeyError("The block function does not support "	2✔
188	"access to zero-blocks. Available non-zero "
189	f"blocks are: {self.blocks_nonzero}.")
190
191	if block == c:	2✔
192	return self._tensors[c]	2✔
193	return block_info.factor * self._tensors[c].transpose(block_info.transpose)	2✔
194
195	def __getitem__(self, blk):	2✔
196	if blk not in self.blocks:	2✔
197	raise KeyError(f"Invalid block {blk} requested. "	2✔
198	f"Available blocks are: {self.blocks}.")
199
200	block_info = self._block_info[blk]	2✔
201	c = block_info.canonical	2✔
202
203	if c not in self._tensors:	2✔
204	sym = libadcc.make_symmetry_operator(	2✔
205	self.mospaces, c, self.symmetry.to_str(), "1"
206	)
207	self._tensors[c] = Tensor(sym)	2✔
208
209	if blk == c:	2✔
210	return self._tensors[c]	2✔
211	return block_info.factor * self._tensors[c].transpose(block_info.transpose)	2✔
212
213	def __getattr__(self, attr) -> Tensor:	2✔
214	from . import block as b	2✔
215	return self.__getitem__(b.__getattr__(attr))	2✔
216
217	def __setitem__(self, block, tensor):	2✔
218	block_info = self._block_info.get(block)	2✔
219	if block_info is None or block_info.canonical != block:	2✔
220	raise KeyError(f"Invalid block {block} assigned. "	2✔
221	f"Available blocks are: {self.blocks}.")
222
223	spaces = split_spaces(block)	2✔
224	expected_shape = tuple(self.mospaces.n_orbs(space) for space in spaces)	2✔
225	if expected_shape != tensor.shape:	2✔
226	raise ValueError("Invalid shape of incoming tensor. "	2✔
227	f"Expected shape {expected_shape}, but "
228	f"got shape {tensor.shape} instead.")
229
230	self._tensors[block] = tensor	2✔
231
232	def __setattr__(self, attr, value):	2✔
233	try:	2✔
234	from . import block as b	2✔
235	self.__setitem__(b.__getattr__(attr), value)	2✔
236	except AttributeError:	2✔
237	super().__setattr__(attr, value)	2✔
238
239	def set_zero_block(self, block):	2✔
240	"""
241	Set a given block as zero block
242	"""
243	if block not in self.canonical_blocks:	2✔
244	raise KeyError(f"Invalid block {block} set as zero block. "	2✔
245	f"Available blocks are: {self.blocks}.")
246	self._tensors.pop(block)	2✔
247
248	def to_ndarray(self) -> np.ndarray:	2✔
249	"""
250	Returns the NParticleOperator as a contiguous
251	np.ndarray instance including all blocks
252	"""
253	# offsets to start index of spaces
254	offsets = {	2✔
255	sp: sum(
256	self.mospaces.n_orbs(ss)
257	for ss in self.orbital_subspaces[:self.orbital_subspaces.index(sp)]
258	)
259	for sp in self.orbital_subspaces
260	}
261	# slices for each space
262	slices = {	2✔
263	sp: slice(offsets[sp], offsets[sp] + self.mospaces.n_orbs(sp))
264	for sp in self.orbital_subspaces
265	}
266	# build complete np.ndarray -> go through all blocks
267	ret = np.zeros((self.shape))	2✔
268	for block in self.blocks:	2✔
269	if self.is_zero_block(block):	2!
NEW 270	continue	×
271	spaces = split_spaces(block)	2✔
272	slice_list = tuple(slices[sp] for sp in spaces)	2✔
273	dm_block = self[block].to_ndarray()	2✔
274	ret[slice_list] = dm_block	2✔
275	return ret	2✔
276
277	def _construct_empty(self):	2✔
278	"""
279	Create an empty instance of an NParticleOperator
280	"""
281	return self.__class__(	2✔
282	self.mospaces,
283	self.n_particle_op,
284	symmetry=self.symmetry,
285	)
286
287	def copy(self):	2✔
288	"""
289	Return a deep copy of the NParticleOperator
290	"""
291	ret = self._construct_empty()	2✔
292	for b in self.blocks_nonzero:	2✔
293	ret[b] = self.block(b).copy()	2✔
294	if hasattr(self, "reference_state"):	2!
295	ret.reference_state = self.reference_state	2✔
296	return ret	2✔
297
298	def _transform_to_ao(self, refstate):	2✔
NEW 299	raise NotImplementedError("Needs to be implemented on the "	×
300	f"{self.__class__.__name__} class.")
301
302	def to_ao_basis(self, refstate_or_coefficients=None):	2✔
303	"""
304	Transforms the NParticleOperator to the atomic orbital
305	basis using a ReferenceState or a coefficient map. If no
306	ReferenceState or coefficient map is given, the ReferenceState
307	used to construct the NParticleOperator is taken instead.
308	"""
309	if isinstance(refstate_or_coefficients, (dict, libadcc.ReferenceState)):	2✔
310	return self._transform_to_ao(refstate_or_coefficients)	2✔
311	elif refstate_or_coefficients is None:	2!
312	if not hasattr(self, "reference_state"):	2!
NEW 313	raise ValueError("Argument reference_state is required if no "	×
314	"reference_state is stored in the "
315	"NParticleOperator")
316	return self._transform_to_ao(self.reference_state)	2✔
317	else:
NEW 318	raise TypeError("Argument type not supported.")	×
319
320	def __iadd__(self, other):	2✔
321	from . import block as b	2✔
322	if self.mospaces != other.mospaces:	2!
NEW 323	raise ValueError(f"Cannot add {self.__class__.__name__}s with "	×
324	"differing mospaces.")
325	if self.symmetry is not OperatorSymmetry.NOSYMMETRY and \	2✔
326	self.symmetry is not other.symmetry:
327	raise ValueError("Cannot add non-symmetric matrix "	2✔
328	"in-place to symmetric one.")
329
330	for block in other.blocks_nonzero:	2✔
331	if self.is_zero_block(block):	2!
NEW 332	self[block] = other.block(block).copy()	×
333	else:
334	self[block] = self.block(block) + other.block(block)	2✔
335
336	if self.symmetry is OperatorSymmetry.NOSYMMETRY \	2✔
337	and other.symmetry is not OperatorSymmetry.NOSYMMETRY:
338	for block in self.blocks_nonzero:	2✔
339	bra = block[:2 * self.n_particle_op]	2✔
340	ket = block[2 * self.n_particle_op:]	2✔
341	c_block, factor, transpose = b.get_canonical_block(	2✔
342	bra, ket, other.symmetry
343	)
344	if block == c_block:	2✔
345	continue # Done already	2✔
346	obT = 0	2✔
347	if not other.is_zero_block(c_block):	2!
348	obT = factor * other.block(c_block).transpose(transpose)	2✔
349	if not self.is_zero_block(block):	2!
350	obT += self.block(block)	2✔
351	self[block] = evaluate(obT)	2✔
352
353	# Update ReferenceState pointer
354	if hasattr(self, "reference_state"):	2!
355	if hasattr(other, "reference_state") \	2!
356	and self.reference_state != other.reference_state:
NEW 357	delattr(self, "reference_state")	×
358	return self	2✔
359
360	def __isub__(self, other):	2✔
361	from . import block as b	2✔
362	if self.mospaces != other.mospaces:	2!
NEW 363	raise ValueError(f"Cannot add {self.__class__.__name__}s with "	×
364	"differing mospaces.")
365	if self.symmetry is not OperatorSymmetry.NOSYMMETRY and \	2✔
366	self.symmetry is not other.symmetry:
367	raise ValueError("Cannot add non-symmetric matrix "	2✔
368	"in-place to symmetric one.")
369
370	for block in other.blocks_nonzero:	2✔
371	if self.is_zero_block(block):	2!
NEW 372	self[block] = -1.0 * other.block(block) # The copy is implicit	×
373	else:
374	self[block] = self.block(block) - other.block(block)	2✔
375
376	if self.symmetry is OperatorSymmetry.NOSYMMETRY \	2✔
377	and other.symmetry is not OperatorSymmetry.NOSYMMETRY:
378	for block in self.blocks_nonzero:	2✔
379	bra = block[:2 * self.n_particle_op]	2✔
380	ket = block[2 * self.n_particle_op:]	2✔
381	c_block, factor, transpose = b.get_canonical_block(	2✔
382	bra, ket, other.symmetry
383	)
384	if block == c_block:	2✔
385	continue # Done already	2✔
386	obT = 0	2✔
387	if not other.is_zero_block(c_block):	2!
388	obT = -1 * factor * other.block(c_block).transpose(transpose)	2✔
389	if not self.is_zero_block(block):	2!
390	obT += self.block(block)	2✔
391	self[block] = evaluate(obT)	2✔
392
393	# Update ReferenceState pointer
394	if hasattr(self, "reference_state"):	2!
395	if hasattr(other, "reference_state") \	2!
396	and self.reference_state != other.reference_state:
NEW 397	delattr(self, "reference_state")	×
398	return self	2✔
399
400	def __imul__(self, other):	2✔
401	if not isinstance(other, (float, int)):	2!
NEW 402	return NotImplemented	×
403	for b in self.blocks_nonzero:	2✔
404	self[b] = self.block(b) * other	2✔
405	return self	2✔
406
407	def __add__(self, other):	2✔
408	if (	2✔
409	self.symmetry is not other.symmetry
410	and self.symmetry is not OperatorSymmetry.NOSYMMETRY
411	and other.symmetry is not OperatorSymmetry.NOSYMMETRY
412	):
413	raise ValueError("Addition of Hermitian and Antihermitian "	2✔
414	"operators is not implemented.")
415
416	if self.symmetry is OperatorSymmetry.NOSYMMETRY \	2✔
417	or other.symmetry is not OperatorSymmetry.NOSYMMETRY:
418	return self.copy().__iadd__(other)	2✔
419	else:
420	return other.copy().__iadd__(self)	2✔
421
422	def __sub__(self, other):	2✔
423	if (	2✔
424	self.symmetry is not other.symmetry
425	and self.symmetry is not OperatorSymmetry.NOSYMMETRY
426	and other.symmetry is not OperatorSymmetry.NOSYMMETRY
427	):
428	raise ValueError("Substraction of Hermitian and Antihermitian "	2✔
429	"operators is not implemented.")
430
431	if self.symmetry is OperatorSymmetry.NOSYMMETRY \	2✔
432	or other.symmetry is not OperatorSymmetry.NOSYMMETRY:
433	return self.copy().__isub__(other)	2✔
434	else:
435	return (-1.0 * other).__iadd__(self)	2✔
436
437	def __mul__(self, other):	2✔
438	return self.copy().__imul__(other)	2✔
439
440	def __rmul__(self, other):	2✔
441	return self.copy().__imul__(other)	2✔
442
443	def evaluate(self):	2✔
444	for b in self.blocks_nonzero:	2✔
445	self.block(b).evaluate()	2✔
446	return self	2✔
447
448	def set_random(self):	2✔
449	for b in self.canonical_blocks:	2✔
450	self[b].set_random()	2✔
451	return self	2✔
452
453
454	def product_trace(op1, op2) -> float:	2✔
455	"""
456	Compute the expectation value (inner product) of two N-particle operators.
457
458	For the special case of a density operator and a general operator, the
459	inner product is identical in the AO and MO basis.
460
461	Parameters
462	----------
463	op1, op2 : NParticleOperator or NParticleDensity
464	Operators or densities to compute the inner product of.
465
466	Returns
467	-------
468	float
469	The trace / expectation value <op1, op2>.
470	"""
471	# TODO use blocks_nonzero and build the set intersection
472	# to avoid the is_zero_block( ) checks below.
473	# I'm a bit hesitant to do this right now, because I'm lacking
474	# the time at the moment to build a more sophisticated test,
475	# which could potentially catch an arising error.
476	if op1.n_particle_op != op2.n_particle_op:	2!
NEW 477	raise TypeError(	×
478	"Both operators must have the same number of N-particle operators "
479	f"(got {op1.n_particle_op} and {op2.n_particle_op})."
480	)
481	assert op1.blocks == op2.blocks	2✔
482	if op1.symmetry == OperatorSymmetry.NOSYMMETRY:	2✔
483	factors = op1.canonical_factors.copy()	2✔
484	elif op2.symmetry == OperatorSymmetry.NOSYMMETRY:	2✔
485	factors = op2.canonical_factors.copy()	2✔
486	else:
487	assert op1.canonical_factors == op2.canonical_factors	2✔
488	factors = op1.canonical_factors.copy()	2✔
489	# Off-diagonal terms can be removed since a Hermitian times
490	# an anti-Hermitian operator yields for off diagonal terms
491	# (e.g. for a one-particle operator) ov * ov and vo * vo = (ov).T * (−ov.T),
492	# which cancel exactly to zero.
493	if op1.symmetry is not op2.symmetry:	2✔
494	to_remove = []	2✔
495	for b in list(factors.keys()):	2✔
496	spaces = split_spaces(b)	2✔
497	n = op1.n_particle_op	2✔
498	if spaces[:n] != spaces[n:]:	2✔
499	to_remove.append(b)	2✔
500	for b in to_remove:	2✔
501	factors.pop(b)	2✔
502	ret = 0	2✔
503	for b, factor in factors.items():	2✔
504	if op1.is_zero_block(b) or op2.is_zero_block(b):	2✔
505	continue	2✔
506
507	ret += (1 / math.factorial(op1.n_particle_op) ** 2 * factor	2✔
508	* op1.block(b).dot(op2.block(b)))
509	return ret	2✔

adc-connect / adcc / 21367620794

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous