21177357553

Committed 20 Jan 2026 03:33PM UTC coverage: 74.598% (+0.3%) from 74.293%

Build # 21177357553

Build Type

Pull #200

github

Committed by

web-flow

Commit Message

Merge 3c14af03c 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 %.

555 of 621 new or added lines in 22 files covered. (89.37%)

3 existing lines in 1 file now uncovered.

7787 of 10141 relevant lines covered (76.79%)

168782.13 hits per line

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

73.74

/adcc/ReferenceState.py

#!/usr/bin/env python3
## vi: tabstop=4 shiftwidth=4 softtabstop=4 expandtab
## ---------------------------------------------------------------------
##
## Copyright (C) 2019 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/>.
##
## ---------------------------------------------------------------------
import numpy as np

from .misc import cached_property
from .Tensor import Tensor
from .MoSpaces import MoSpaces
from .backends import import_scf_results
from .OperatorIntegrals import OperatorIntegrals
from .OneParticleDensity import OneParticleDensity
from .TwoParticleDensity import TwoParticleDensity
from .NParticleOperator import product_trace, OperatorSymmetry

import libadcc

__all__ = ["ReferenceState"]


class ReferenceState(libadcc.ReferenceState):
    def __init__(self, hfdata, core_orbitals=None, frozen_core=None,
                 frozen_virtual=None, symmetry_check_on_import=False,
                 import_all_below_n_orbs=10):
        """Construct a ReferenceState holding information about the employed
        SCF reference.

        The constructed object is lazy and will at construction only setup
        orbital energies and coefficients. Fock matrix blocks and
        electron-repulsion integral blocks are imported as needed.

        Orbital subspace selection: In order to specify `frozen_core`,
        `core_orbitals` and `frozen_virtual`, adcc allows a range of
        specifications including

           a. A number: Just put this number of alpha orbitals and this
              number of beta orbitals into the respective space. For frozen
              core and core orbitals these are counted from below, for
              frozen virtual orbitals, these are counted from above. If both
              frozen core and core orbitals are specified like this, the
              lowest-energy, occupied orbitals will be put into frozen core.
           b. A range: The orbital indices given by this range will be put
              into the orbital subspace.
           c. An explicit list of orbital indices to be placed into the
              subspace.
           d. A pair of (a) to (c): If the orbital selection for alpha and
              beta orbitals should differ, a pair of ranges, or a pair of
              index lists or a pair of numbers can be specified.

        Parameters
        ----------
        hfdata
            Object with Hartree-Fock data (e.g. a molsturm scf state, a pyscf
            SCF object or any class implementing the
            :py:class:`adcc.HartreeFockProvider` interface or in fact any python
            object representing a pointer to a C++ object derived off
            the :cpp:class:`adcc::HartreeFockSolution_i`.

        core_orbitals : int or list or tuple, optional
            The orbitals to be put into the core-occupied space. For ways to
            define the core orbitals see the description above.

        frozen_core : int or list or tuple, optional
            The orbitals to be put into the frozen core space. For ways to
            define the core orbitals see the description above. For an automatic
            selection of the frozen core space one may also specify
            ``frozen_core=True``.

        frozen_virtuals : int or list or tuple, optional
            The orbitals to be put into the frozen virtual space. For ways to
            define the core orbitals see the description above.

        symmetry_check_on_import : bool, optional
            Should symmetry of the imported objects be checked explicitly during
            the import process. This massively slows down the import and has a
            dramatic impact on memory usage. Thus one should enable this only
            for debugging (e.g. for testing import routines from the host
            programs). Do not enable this unless you know what you are doing.

        import_all_below_n_orbs : int, optional
            For small problem sizes lazy make less sense, since the memory
            requirement for storing the ERI tensor is neglibile and thus the
            flexiblity gained by having the full tensor in memory is
            advantageous. Below the number of orbitals specified by this
            parameter, the class will thus automatically import all ERI tensor
            and Fock matrix blocks.

        Examples
        --------
        To start a calculation with the 2 lowest alpha and beta orbitals
        in the core occupied space, construct the class as

        >>> ReferenceState(hfdata, core_orbitals=2)

        or

        >>> ReferenceState(hfdata, core_orbitals=range(2))

        or

        >>> ReferenceState(hfdata, core_orbitals=[0, 1])

        or

        >>> ReferenceState(hfdata, core_orbitals=([0, 1], [0, 1]))

        There is no restriction to choose the core occupied orbitals
        from the bottom end of the occupied orbitals. For example
        to select the 2nd and 3rd orbital setup the class as

        >>> ReferenceState(hfdata, core_orbitals=range(1, 3))

        or

        >>> ReferenceState(hfdata, core_orbitals=[1, 2])

        If different orbitals should be placed in the alpha and
        beta orbitals, this can be achievd like so

        >>> ReferenceState(hfdata, core_orbitals=([1, 2], [0, 1]))

        which would place the 2nd and 3rd alpha and the 1st and second
        beta orbital into the core space.
        """
        if not isinstance(hfdata, libadcc.HartreeFockSolution_i):
            hfdata = import_scf_results(hfdata)

        self._mospaces = MoSpaces(hfdata, frozen_core=frozen_core,
                                  frozen_virtual=frozen_virtual,
                                  core_orbitals=core_orbitals)
        super().__init__(hfdata, self._mospaces, symmetry_check_on_import)

        if import_all_below_n_orbs is not None and \
           hfdata.n_orbs < import_all_below_n_orbs:
            super().import_all()

        self.operators = OperatorIntegrals(
            hfdata.operator_integral_provider, self._mospaces,
            self.orbital_coefficients, self.conv_tol
        )

        self.environment = None  # no environment attached by default
        for name in ["excitation_energy_corrections", "environment"]:
            if hasattr(hfdata, name):
                setattr(self, name, getattr(hfdata, name))

    def __getattr__(self, attr):
        from . import block as b

        if attr.startswith("f"):
            return self.fock(b.__getattr__(attr[1:]))
        else:
            return self.eri(b.__getattr__(attr))

    @property
    def mospaces(self):
        return self._mospaces

    @property
    def timer(self):
        ret = super().timer
        ret.attach(self.operators.timer)
        return ret

    @property
    def is_aufbau_occupation(self):
        """
        Returns whether the molecular orbital occupation in this reference
        is according to the Aufbau principle (lowest-energy orbitals are occupied)
        """
        eHOMO = max(max(self.orbital_energies(space).to_ndarray())
                    for space in self.mospaces.subspaces_occupied)
        eLUMO = min(min(self.orbital_energies(space).to_ndarray())
                    for space in self.mospaces.subspaces_virtual)
        return eHOMO < eLUMO

    def to_qcvars(self, properties=False, recurse=False):
        """
        Return a dictionary with property keys compatible to a Psi4 wavefunction
        or a QCEngine Atomicresults object.
        """
        qcvars = {"HF TOTAL ENERGY": self.energy_scf, }
        if properties:
            qcvars["HF DIPOLE"] = self.dipole_moment
        return qcvars

    @property
    def density(self):
        """
        Return the Hartree-Fock density in the MO basis
        """
        density = OneParticleDensity(self.mospaces,
                                     symmetry=OperatorSymmetry.HERMITIAN)
        for block in density.canonical_blocks:
            sym = libadcc.make_symmetry_operator(self.mospaces, block,
                                                 density.symmetry.to_str(), "1")
            density[block] = Tensor(sym)
        for ss in self.mospaces.subspaces_occupied:
            density[ss + ss].set_mask("ii", 1)
        density.reference_state = self
        return density

    @property
    def density_2p(self):
        """
        Return the two-particle Hartree-Fock density in the MO basis
        """
        density = TwoParticleDensity(self.mospaces,
                                     symmetry=OperatorSymmetry.HERMITIAN)
        for block in density.canonical_blocks:
            sym = libadcc.make_symmetry_operator(self.mospaces, block,
                                                 density.symmetry.to_str(), "1")
            density[block] = Tensor(sym)
        for ss in self.mospaces.subspaces_occupied:
            density[ss + ss + ss + ss].set_mask("ijij", 1)
            density[ss + ss + ss + ss].set_mask("ijji", -1)
            density[ss + ss + ss + ss].set_mask("iijj", 0)
        density.reference_state = self
        return density

    @cached_property
    def dipole_moment(self):
        """
        Return the HF dipole moment of the reference state (that is the sum of
        the electronic and the nuclear contribution.)
        """
        dipole_integrals = self.operators.electric_dipole
        return self.nuclear_dipole + np.array([product_trace(comp, self.density)
                                               for comp in dipole_integrals])

    def nuclear_quadrupole(self, gauge_origin="origin"):
        if isinstance(gauge_origin, str):
            gauge_origin = self.gauge_origin_to_xyz(gauge_origin)
        return super().nuclear_quadrupole(gauge_origin)

# TODO some nice describe method

1	#!/usr/bin/env python3
2	## vi: tabstop=4 shiftwidth=4 softtabstop=4 expandtab
3	## ---------------------------------------------------------------------
4	##
5	## Copyright (C) 2019 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	import numpy as np	2✔
24
25	from .misc import cached_property	2✔
26	from .Tensor import Tensor	2✔
27	from .MoSpaces import MoSpaces	2✔
28	from .backends import import_scf_results	2✔
29	from .OperatorIntegrals import OperatorIntegrals	2✔
30	from .OneParticleDensity import OneParticleDensity	2✔
31	from .TwoParticleDensity import TwoParticleDensity	2✔
32	from .NParticleOperator import product_trace, OperatorSymmetry	2✔
33
34	import libadcc	2✔
35
36	__all__ = ["ReferenceState"]	2✔
37
38
39	class ReferenceState(libadcc.ReferenceState):	2✔
40	def __init__(self, hfdata, core_orbitals=None, frozen_core=None,	2✔
41	frozen_virtual=None, symmetry_check_on_import=False,
42	import_all_below_n_orbs=10):
43	"""Construct a ReferenceState holding information about the employed
44	SCF reference.
45
46	The constructed object is lazy and will at construction only setup
47	orbital energies and coefficients. Fock matrix blocks and
48	electron-repulsion integral blocks are imported as needed.
49
50	Orbital subspace selection: In order to specify `frozen_core`,
51	`core_orbitals` and `frozen_virtual`, adcc allows a range of
52	specifications including
53
54	a. A number: Just put this number of alpha orbitals and this
55	number of beta orbitals into the respective space. For frozen
56	core and core orbitals these are counted from below, for
57	frozen virtual orbitals, these are counted from above. If both
58	frozen core and core orbitals are specified like this, the
59	lowest-energy, occupied orbitals will be put into frozen core.
60	b. A range: The orbital indices given by this range will be put
61	into the orbital subspace.
62	c. An explicit list of orbital indices to be placed into the
63	subspace.
64	d. A pair of (a) to (c): If the orbital selection for alpha and
65	beta orbitals should differ, a pair of ranges, or a pair of
66	index lists or a pair of numbers can be specified.
67
68	Parameters
69	----------
70	hfdata
71	Object with Hartree-Fock data (e.g. a molsturm scf state, a pyscf
72	SCF object or any class implementing the
73	:py:class:`adcc.HartreeFockProvider` interface or in fact any python
74	object representing a pointer to a C++ object derived off
75	the :cpp:class:`adcc::HartreeFockSolution_i`.
76
77	core_orbitals : int or list or tuple, optional
78	The orbitals to be put into the core-occupied space. For ways to
79	define the core orbitals see the description above.
80
81	frozen_core : int or list or tuple, optional
82	The orbitals to be put into the frozen core space. For ways to
83	define the core orbitals see the description above. For an automatic
84	selection of the frozen core space one may also specify
85	``frozen_core=True``.
86
87	frozen_virtuals : int or list or tuple, optional
88	The orbitals to be put into the frozen virtual space. For ways to
89	define the core orbitals see the description above.
90
91	symmetry_check_on_import : bool, optional
92	Should symmetry of the imported objects be checked explicitly during
93	the import process. This massively slows down the import and has a
94	dramatic impact on memory usage. Thus one should enable this only
95	for debugging (e.g. for testing import routines from the host
96	programs). Do not enable this unless you know what you are doing.
97
98	import_all_below_n_orbs : int, optional
99	For small problem sizes lazy make less sense, since the memory
100	requirement for storing the ERI tensor is neglibile and thus the
101	flexiblity gained by having the full tensor in memory is
102	advantageous. Below the number of orbitals specified by this
103	parameter, the class will thus automatically import all ERI tensor
104	and Fock matrix blocks.
105
106	Examples
107	--------
108	To start a calculation with the 2 lowest alpha and beta orbitals
109	in the core occupied space, construct the class as
110
111	>>> ReferenceState(hfdata, core_orbitals=2)
112
113	or
114
115	>>> ReferenceState(hfdata, core_orbitals=range(2))
116
117	or
118
119	>>> ReferenceState(hfdata, core_orbitals=[0, 1])
120
121	or
122
123	>>> ReferenceState(hfdata, core_orbitals=([0, 1], [0, 1]))
124
125	There is no restriction to choose the core occupied orbitals
126	from the bottom end of the occupied orbitals. For example
127	to select the 2nd and 3rd orbital setup the class as
128
129	>>> ReferenceState(hfdata, core_orbitals=range(1, 3))
130
131	or
132
133	>>> ReferenceState(hfdata, core_orbitals=[1, 2])
134
135	If different orbitals should be placed in the alpha and
136	beta orbitals, this can be achievd like so
137
138	>>> ReferenceState(hfdata, core_orbitals=([1, 2], [0, 1]))
139
140	which would place the 2nd and 3rd alpha and the 1st and second
141	beta orbital into the core space.
142	"""
143	if not isinstance(hfdata, libadcc.HartreeFockSolution_i):	2✔
144	hfdata = import_scf_results(hfdata)	2✔
145
146	self._mospaces = MoSpaces(hfdata, frozen_core=frozen_core,	2✔
147	frozen_virtual=frozen_virtual,
148	core_orbitals=core_orbitals)
149	super().__init__(hfdata, self._mospaces, symmetry_check_on_import)	2✔
150
151	if import_all_below_n_orbs is not None and \	2!
152	hfdata.n_orbs < import_all_below_n_orbs:
153	super().import_all()	×
154
155	self.operators = OperatorIntegrals(	2✔
156	hfdata.operator_integral_provider, self._mospaces,
157	self.orbital_coefficients, self.conv_tol
158	)
159
160	self.environment = None # no environment attached by default	2✔
161	for name in ["excitation_energy_corrections", "environment"]:	2✔
162	if hasattr(hfdata, name):	2✔
163	setattr(self, name, getattr(hfdata, name))	2✔
164
165	def __getattr__(self, attr):	2✔
166	from . import block as b	2✔
167
168	if attr.startswith("f"):	2✔
169	return self.fock(b.__getattr__(attr[1:]))	2✔
170	else:
171	return self.eri(b.__getattr__(attr))	2✔
172
173	@property	2✔
174	def mospaces(self):	2✔
175	return self._mospaces	2✔
176
177	@property	2✔
178	def timer(self):	2✔
179	ret = super().timer	2✔
180	ret.attach(self.operators.timer)	2✔
181	return ret	2✔
182
183	@property	2✔
184	def is_aufbau_occupation(self):	2✔
185	"""
186	Returns whether the molecular orbital occupation in this reference
187	is according to the Aufbau principle (lowest-energy orbitals are occupied)
188	"""
189	eHOMO = max(max(self.orbital_energies(space).to_ndarray())	×
190	for space in self.mospaces.subspaces_occupied)
191	eLUMO = min(min(self.orbital_energies(space).to_ndarray())	×
192	for space in self.mospaces.subspaces_virtual)
193	return eHOMO < eLUMO	×
194
195	def to_qcvars(self, properties=False, recurse=False):	2✔
196	"""
197	Return a dictionary with property keys compatible to a Psi4 wavefunction
198	or a QCEngine Atomicresults object.
199	"""
200	qcvars = {"HF TOTAL ENERGY": self.energy_scf, }	×
201	if properties:	×
202	qcvars["HF DIPOLE"] = self.dipole_moment	×
203	return qcvars	×
204
205	@property	2✔
206	def density(self):	2✔
207	"""
208	Return the Hartree-Fock density in the MO basis
209	"""
210	density = OneParticleDensity(self.mospaces,	2✔
211	symmetry=OperatorSymmetry.HERMITIAN)
212	for block in density.canonical_blocks:	2✔
213	sym = libadcc.make_symmetry_operator(self.mospaces, block,	2✔
214	density.symmetry.to_str(), "1")
215	density[block] = Tensor(sym)	2✔
216	for ss in self.mospaces.subspaces_occupied:	2✔
217	density[ss + ss].set_mask("ii", 1)	2✔
218	density.reference_state = self	2✔
219	return density	2✔
220
221	@property	2✔
222	def density_2p(self):	2✔
223	"""
224	Return the two-particle Hartree-Fock density in the MO basis
225	"""
NEW 226	density = TwoParticleDensity(self.mospaces,	×
227	symmetry=OperatorSymmetry.HERMITIAN)
NEW 228	for block in density.canonical_blocks:	×
NEW 229	sym = libadcc.make_symmetry_operator(self.mospaces, block,	×
230	density.symmetry.to_str(), "1")
NEW 231	density[block] = Tensor(sym)	×
NEW 232	for ss in self.mospaces.subspaces_occupied:	×
NEW 233	density[ss + ss + ss + ss].set_mask("ijij", 1)	×
NEW 234	density[ss + ss + ss + ss].set_mask("ijji", -1)	×
NEW 235	density[ss + ss + ss + ss].set_mask("iijj", 0)	×
NEW 236	density.reference_state = self	×
NEW 237	return density	×
238
239	@cached_property	2✔
240	def dipole_moment(self):	2✔
241	"""
242	Return the HF dipole moment of the reference state (that is the sum of
243	the electronic and the nuclear contribution.)
244	"""
245	dipole_integrals = self.operators.electric_dipole	2✔
246	return self.nuclear_dipole + np.array([product_trace(comp, self.density)	2✔
247	for comp in dipole_integrals])
248
249	def nuclear_quadrupole(self, gauge_origin="origin"):	2✔
250	if isinstance(gauge_origin, str):	2!
251	gauge_origin = self.gauge_origin_to_xyz(gauge_origin)	2✔
252	return super().nuclear_quadrupole(gauge_origin)	2✔
253
254	# TODO some nice describe method

adc-connect / adcc / 21177357553

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