22406590854

Committed 25 Feb 2026 04:44PM UTC coverage: 74.749% (+0.5%) from 74.248%

Build # 22406590854

Build Type

Pull #202

github

Committed by

web-flow

Commit Message

Merge a12231cd6 into b49d31427

Pull Request Pull Request #202: <S^2> calculation for UHF

Coverage Stats

1299 of 2028 branches covered (64.05%)

Branch coverage included in aggregate %.

189 of 207 new or added lines in 9 files covered. (91.3%)

1 existing line in 1 file now uncovered.

8100 of 10546 relevant lines covered (76.81%)

169226.74 hits per line

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

87.27

/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.orbital_coefficients_alpha,
            self.orbital_coefficients_beta, 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)

    @cached_property
    def ssq(self):
        """
        Return <S^2> of the HF reference state.
        """
        if self.restricted:
            raise NotImplementedError(
                "<S^2> is not implemented for restricted HF references."
            )
        ssq_1p_op = self.operators.ssq_1p
        ssq_2p_op = self.operators.ssq_2p
        ssq_1p = product_trace(ssq_1p_op, self.density)
        ssq_2p = product_trace(ssq_2p_op, self.density_2p)
        return (ssq_1p + ssq_2p)

# 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.orbital_coefficients_alpha,
158	self.orbital_coefficients_beta, self.conv_tol
159	)
160
161	self.environment = None # no environment attached by default	2✔
162	for name in ["excitation_energy_corrections", "environment"]:	2✔
163	if hasattr(hfdata, name):	2✔
164	setattr(self, name, getattr(hfdata, name))	2✔
165
166	def __getattr__(self, attr):	2✔
167	from . import block as b	2✔
168
169	if attr.startswith("f"):	2✔
170	return self.fock(b.__getattr__(attr[1:]))	2✔
171	else:
172	return self.eri(b.__getattr__(attr))	2✔
173
174	@property	2✔
175	def mospaces(self):	2✔
176	return self._mospaces	2✔
177
178	@property	2✔
179	def timer(self):	2✔
180	ret = super().timer	2✔
181	ret.attach(self.operators.timer)	2✔
182	return ret	2✔
183
184	@property	2✔
185	def is_aufbau_occupation(self):	2✔
186	"""
187	Returns whether the molecular orbital occupation in this reference
188	is according to the Aufbau principle (lowest-energy orbitals are occupied)
189	"""
190	eHOMO = max(max(self.orbital_energies(space).to_ndarray())	×
191	for space in self.mospaces.subspaces_occupied)
192	eLUMO = min(min(self.orbital_energies(space).to_ndarray())	×
193	for space in self.mospaces.subspaces_virtual)
194	return eHOMO < eLUMO	×
195
196	def to_qcvars(self, properties=False, recurse=False):	2✔
197	"""
198	Return a dictionary with property keys compatible to a Psi4 wavefunction
199	or a QCEngine Atomicresults object.
200	"""
201	qcvars = {"HF TOTAL ENERGY": self.energy_scf, }	×
202	if properties:	×
203	qcvars["HF DIPOLE"] = self.dipole_moment	×
204	return qcvars	×
205
206	@property	2✔
207	def density(self):	2✔
208	"""
209	Return the Hartree-Fock density in the MO basis
210	"""
211	density = OneParticleDensity(self.mospaces,	2✔
212	symmetry=OperatorSymmetry.HERMITIAN)
213	for block in density.canonical_blocks:	2✔
214	sym = libadcc.make_symmetry_operator(self.mospaces, block,	2✔
215	density.symmetry.to_str(), "1")
216	density[block] = Tensor(sym)	2✔
217	for ss in self.mospaces.subspaces_occupied:	2✔
218	density[ss + ss].set_mask("ii", 1)	2✔
219	density.reference_state = self	2✔
220	return density	2✔
221
222	@property	2✔
223	def density_2p(self):	2✔
224	"""
225	Return the two-particle Hartree-Fock density in the MO basis
226	"""
227	density = TwoParticleDensity(self.mospaces,	2✔
228	symmetry=OperatorSymmetry.HERMITIAN)
229	for block in density.canonical_blocks:	2✔
230	sym = libadcc.make_symmetry_operator(self.mospaces, block,	2✔
231	density.symmetry.to_str(), "1")
232	density[block] = Tensor(sym)	2✔
233	for ss in self.mospaces.subspaces_occupied:	2✔
234	density[ss + ss + ss + ss].set_mask("ijij", 1)	2✔
235	density[ss + ss + ss + ss].set_mask("ijji", -1)	2✔
236	density[ss + ss + ss + ss].set_mask("iijj", 0)	2✔
237	density.reference_state = self	2✔
238	return density	2✔
239
240	@cached_property	2✔
241	def dipole_moment(self):	2✔
242	"""
243	Return the HF dipole moment of the reference state (that is the sum of
244	the electronic and the nuclear contribution.)
245	"""
246	dipole_integrals = self.operators.electric_dipole	2✔
247	return self.nuclear_dipole + np.array([product_trace(comp, self.density)	2✔
248	for comp in dipole_integrals])
249
250	def nuclear_quadrupole(self, gauge_origin="origin"):	2✔
251	if isinstance(gauge_origin, str):	2!
252	gauge_origin = self.gauge_origin_to_xyz(gauge_origin)	2✔
253	return super().nuclear_quadrupole(gauge_origin)	2✔
254
255	@cached_property	2✔
256	def ssq(self):	2✔
257	"""
258	Return <S^2> of the HF reference state.
259	"""
260	if self.restricted:	2!
NEW 261	raise NotImplementedError(	×
262	"<S^2> is not implemented for restricted HF references."
263	)
264	ssq_1p_op = self.operators.ssq_1p	2✔
265	ssq_2p_op = self.operators.ssq_2p	2✔
266	ssq_1p = product_trace(ssq_1p_op, self.density)	2✔
267	ssq_2p = product_trace(ssq_2p_op, self.density_2p)	2✔
268	return (ssq_1p + ssq_2p)	2✔
269
270	# TODO some nice describe method

adc-connect / adcc / 22406590854

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