21211666531

Committed 21 Jan 2026 01:38PM UTC coverage: 74.578% (+0.3%) from 74.293%

Build # 21211666531

Build Type

Pull #200

github

Committed by

web-flow

Commit Message

Merge cb6bc1db0 into 9ff4abd52

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

Run Details

1270 of 2000 branches covered (63.5%)

Branch coverage included in aggregate %.

552 of 621 new or added lines in 22 files covered. (88.89%)

2 existing lines in 1 file now uncovered.

7783 of 10139 relevant lines covered (76.76%)

169253.42 hits per line

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

90.91

/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 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()


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, reverse=True)
        virts = sorted(self.mospaces.subspaces_virtual, reverse=True)
        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]
        self.canonical_blocks = []
        self.canonical_factors = {}

        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, _, _ = get_canonical_block(bra, ket, self.symmetry)
            if canonical_block not in self.canonical_blocks:
                self.canonical_blocks.append(canonical_block)
                self.canonical_factors[canonical_block] = 1
            else:
                self.canonical_factors[canonical_block] += 1

    @property
    def n_particle_op(self) -> int:
        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.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.
        """
        if block not in self.canonical_blocks:
            return False
        return block not in self.blocks_nonzero

    def block(self, block) -> Tensor:
        """
        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.
        """
        from . import block as b
        if block in self.blocks_nonzero:
            return self._tensors[block]
        bra = block[:2 * self.n_particle_op]
        ket = block[2 * self.n_particle_op:]
        canonical_block, factor, transpose = b.get_canonical_block(
            bra, ket, self.symmetry
        )
        if canonical_block in self.blocks_nonzero:
            return factor * self._tensors[canonical_block].transpose(transpose)
        else:
            raise KeyError("The block function does not support "
                           "access to zero-blocks. Available non-zero "
                           f"blocks are: {self.blocks_nonzero}.")
        pass

    def __getitem__(self, blk) -> Tensor:
        if blk not in self.blocks:
            raise KeyError(f"Invalid block {blk} requested. "
                           f"Available blocks are: {self.blocks}.")
        if blk not in self._tensors:
            if blk in self.canonical_blocks:
                sym = libadcc.make_symmetry_operator(
                    self.mospaces, blk, self.symmetry.to_str(), "1"
                )
                self._tensors[blk] = Tensor(sym)
                return self._tensors[blk]
            else:
                from . import block as b
                bra = blk[:2 * self.n_particle_op]
                ket = blk[2 * self.n_particle_op:]
                canonical_block, factor, transpose = b.get_canonical_block(
                    bra, ket, self.symmetry
                )
                if canonical_block not in self._tensors:
                    sym = libadcc.make_symmetry_operator(
                        self.mospaces, canonical_block, self.symmetry.to_str(), "1"
                    )
                    self._tensors[canonical_block] = Tensor(sym)
                return factor * self._tensors[canonical_block].transpose(transpose)
        return self._tensors[blk]

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

    def __setitem__(self, block, tensor):
        """
        Assigns a tensor to the specified block
        """
        if block not in self.blocks:
            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.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) -> tuple[Tensor, Tensor]:
        raise NotImplementedError("Needs to be implemented on the Child 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("Cannot add OneParticleOperators with "
                             "differing mospaces.")
        if self.symmetry is not OperatorSymmetry.NOSYMMETRY \
                and other.symmetry == OperatorSymmetry.NOSYMMETRY:
            raise ValueError("Cannot add non-symmetric matrix "
                             "in-place to symmetric one.")
        if (
            self.symmetry != other.symmetry
            and self.symmetry != OperatorSymmetry.NOSYMMETRY
            and other.symmetry != OperatorSymmetry.NOSYMMETRY
        ):
            raise ValueError("Cannot add Hermitian and "
                             "anti-Hermitian matrices in-place.")

        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 == OperatorSymmetry.NOSYMMETRY \
                and other.symmetry != 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("Cannot add OneParticleOperators with "
                             "differing mospaces.")
        if self.symmetry is not OperatorSymmetry.NOSYMMETRY \
                and other.symmetry == OperatorSymmetry.NOSYMMETRY:
            raise ValueError("Cannot add non-symmetric matrix "
                             "in-place to symmetric one.")
        if (
            self.symmetry != other.symmetry
            and self.symmetry != OperatorSymmetry.NOSYMMETRY
            and other.symmetry != OperatorSymmetry.NOSYMMETRY
        ):
            raise ValueError("Cannot add Hermitian and "
                             "anti-Hermitian matrices in-place.")

        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 == 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
        # raise NotImplementedError

    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 != other.symmetry
            and self.symmetry != OperatorSymmetry.NOSYMMETRY
            and other.symmetry != OperatorSymmetry.NOSYMMETRY
        ):
            other = other._expand_to_nosymmetry()

        if self.symmetry == 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 != other.symmetry
            and self.symmetry != OperatorSymmetry.NOSYMMETRY
            and other.symmetry != OperatorSymmetry.NOSYMMETRY
        ):
            other = other._expand_to_nosymmetry()

        if self.symmetry == 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 _expand_to_nosymmetry(self):
        from . import block as b

        if self.symmetry == OperatorSymmetry.NOSYMMETRY:
            return self

        sym = self.symmetry

        for block in self.blocks:
            if self.is_zero_block(block):
                continue

            bra = block[:2 * self.n_particle_op]
            ket = block[2 * self.n_particle_op:]

            c_block, factor, transpose = b.get_canonical_block(
                bra, ket, sym
            )

            if block == c_block:
                continue

            if self.is_zero_block(c_block):
                continue

            value = factor * self.block(c_block).transpose(transpose)
            self[block] = evaluate(value)

        self.symmetry = OperatorSymmetry.NOSYMMETRY

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

adc-connect / adcc / 21211666531

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