11687183074

Committed 05 Nov 2024 03:20PM UTC coverage: 73.291% (-0.5%) from 73.743%

Build # 11687183074

Build Type

Pull #169

github

Committed by

web-flow

Commit Message

Merge e14f10df3 into e2e1bce49

Pull Request Pull Request #169: Additional operator integrals and selection of gauge origin

Run Details

1166 of 1868 branches covered (62.42%)

Branch coverage included in aggregate %.

132 of 233 new or added lines in 12 files covered. (56.65%)

5 existing lines in 3 files now uncovered.

7313 of 9701 relevant lines covered (75.38%)

271899.11 hits per line

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

80.25

/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 .OneParticleOperator import OneParticleOperator, product_trace

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 = OneParticleOperator(self.mospaces, is_symmetric=True)
        for block in density.blocks:
            sym = libadcc.make_symmetry_operator(self.mospaces, block, True, "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

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

adc-connect / adcc / 11687183074

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