Source code for iodata.orbitals

# IODATA is an input and output module for quantum chemistry.
# Copyright (C) 2011-2019 The IODATA Development Team
#
# This file is part of IODATA.
#
# IODATA 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.
#
# IODATA 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 this program; if not, see <http://www.gnu.org/licenses/>
# --
"""Data structure for molecular orbitals."""

from typing import Optional

import attrs
import numpy as np
from numpy.typing import NDArray

from .attrutils import convert_array_to, validate_shape

__all__ = ["MolecularOrbitals"]


[docs] def validate_norbab(mo, attribute, value): """Validate the norba or norbb value assigned to a MolecularOrbitals object. Parameters ---------- mo The MolecularOrbitals instance. attribute Attribute instancce being changed. value The new value. """ if mo.kind == "generalized": if value is not None: raise ValueError( f"Attribute {attribute.name} must be None in case of generalized orbitals." ) return if value is None: raise ValueError( f"Attribute {attribute.name} cannot be None in case of (un)restricted orbitals." ) if mo.kind == "restricted": norb_other = mo.norbb if (attribute.name == "norba") else mo.norba if value != norb_other: raise ValueError("In case of restricted orbitals, norba must be equal to norbb.")
[docs] def validate_occs_aminusb(mo, _attribtue, value): """Validate the occs_aminusb attribute.""" if mo.kind != "restricted" and value is not None: raise ValueError("Attribute occs_aminusb can only be set for restricted wavefunctions.")
[docs] @attrs.define class MolecularOrbitals: """Class of Orthonormal Molecular Orbitals. Notes ----- For restricted wavefunctions, the occupation numbers are spin-summed values and several rules are used to deduce the alpha and beta occupation numbers: - When ``occs_aminusb`` is set, alpha and beta occupation numbers are derived trivially as ``(occs + occs_aminusb) / 2`` and ``(occs - occs_aminusb) / 2``, respectively. - When ``occs_aminusb`` is not set, there are two possibilities. When the occupation numbers are integers, it is assumed that the orbitals represent a restricted open-shell HF or KS wavefunction. An occupation number of 1 is then interpreted as an occupied alpha orbital and a virtual beta orbital. When the occupation numbers are fractional, it is assumed that the orbitals are closed-shell natural orbitals. One can always describe all cases by setting ``occs_aminusb``. While this seems appealing, keep in mind that most wavefunction file formats (FCHK, Molden, Molekel, WFN and WFX) do not support it. """ kind: str = attrs.field( validator=attrs.validators.in_(["restricted", "unrestricted", "generalized"]) ) """Type of molecular orbitals, which can be 'restricted', 'unrestricted', or 'generalized'.""" norba: int = attrs.field(validator=validate_norbab) """ Number of (occupied and virtual) alpha molecular orbitals. Set to `None` in case oftype=='generalized'. """ norbb: int = attrs.field(validator=validate_norbab) """ Number of (occupied and virtual) beta molecular orbitals. Set to `None` in case of type=='generalized'. This is expected to be equal to `norba` for the `restricted` kind. """ occs: Optional[NDArray[float]] = attrs.field( default=None, converter=convert_array_to(float), validator=attrs.validators.optional(validate_shape("norb")), ) """ Molecular orbital occupation numbers. The length equals the number of columns of coeffs. (optional) """ coeffs: Optional[NDArray[float]] = attrs.field( default=None, converter=convert_array_to(float), validator=attrs.validators.optional(validate_shape(None, "norb")), ) """ Molecular orbital coefficients. In case of restricted: shape = (nbasis, norba) = (nbasis, norbb). In case of unrestricted: shape = (nbasis, norba + norbb). In case of generalized: shape = (2 * nbasis, norb), where norb is the total number of orbitals. (optional) """ energies: Optional[NDArray[float]] = attrs.field( default=None, converter=convert_array_to(float), validator=attrs.validators.optional(validate_shape("norb")), ) """Molecular orbital energies. The length equals the number of columns of coeffs. (optional)""" irreps: Optional[NDArray] = attrs.field( default=None, validator=attrs.validators.optional(validate_shape("norb")) ) """Irreducible representation. The length equals the number of columns of coeffs. (optional)""" occs_aminusb: Optional[NDArray[float]] = attrs.field( default=None, converter=convert_array_to(float), validator=attrs.validators.and_( attrs.validators.optional(validate_shape("norb")), validate_occs_aminusb ), ) """ The difference between alpha and beta occupation numbers. The length equals the number of columns of coeffs. (optional and only allowed to be not None for restricted wavefunctions) """ @property def nelec(self) -> float: """Return the total number of electrons.""" if self.occs is None: return None return self.occs.sum() @property def nbasis(self): """Return the number of spatial basis functions.""" if self.coeffs is None: return None if self.kind == "generalized": return self.coeffs.shape[0] // 2 return self.coeffs.shape[0] @property def norb(self): """Return the number of spatially distinct orbitals. Notes ----- In case of restricted wavefunctions, this may be less than just the sum of ``norba`` and ``norbb``, because alpha and beta orbitals share the same spatical dependence. """ if self.kind == "restricted": return self.norba if self.kind == "unrestricted": return self.norba + self.norbb if self.coeffs is not None: return self.coeffs.shape[1] if self.occs is not None: return self.occs.shape[0] if self.energies is not None: return self.energies.shape[0] if self.irreps is not None: return len(self.irreps) return None @property def spinpol(self) -> float: """Return the spin polarization of the Slater determinant.""" if self.kind == "generalized": raise NotImplementedError if self.occs is None: return None if self.kind == "restricted": if self.occs_aminusb is None: # heuristics ... if (self.occs == self.occs.astype(int)).all(): # restricted open-shell HF/KS nbeta = np.clip(self.occs, 0, 1).sum() return abs(self.nelec - 2 * nbeta) # restricted closed-shell natural orbitals return 0.0 return self.occs_aminusb.sum() return abs(self.occsa.sum() - self.occsb.sum()) @property def occsa(self): """Return alpha occupation numbers. Notes ----- For restricted wavefunctions, in-place assignment to occsa will not work. In this case, the array is derived from ``mo.occs`` and optionally ``mo.occs_aminusb``. To avoid that in-place assignment of occsa is silently ignored, it is returned as a non-writeable array. To change occsa, one can assign a whole new array, e.g. ``mo.occsa = new_occsa`` will work, while ``mo.occsa[1] = 0.3`` will not. """ if self.kind == "generalized": raise NotImplementedError if self.occs is None: return None if self.kind == "restricted": if self.occs_aminusb is None: # heuristics ... if (self.occs == self.occs.astype(int)).all(): # restricted open-shell HF/KS result = np.clip(self.occs, 0, 1) else: # restricted closed-shell natural orbitals result = self.occs / 2 else: result = (self.occs + self.occs_aminusb) / 2 result.flags.writeable = False return result return self.occs[: self.norba] @occsa.setter def occsa(self, occsa): if self.kind == "generalized": raise NotImplementedError if self.kind == "restricted": occsa = np.array(occsa) if self.occs is None: self.occs = occsa self.occs_aminusb = occsa.copy() else: occsb = np.array(self.occsb) self.occs = occsa + occsb self.occs_aminusb = occsa - occsb else: self.occs[: self.norba] = occsa @property def occsb(self): """Return beta occupation numbers. Notes ----- For restricted wavefunctions, in-place assignment to occsb will not work. In this case, the array is derived from ``mo.occs`` and optionally ``mo.occs_aminusb``. To avoid that in-place assignment of occsb is silently ignored, it is returned as a non-writeable array. To change occsb, one can assign a whole new array, e.g. ``mo.occsb = new_occsb`` will work, while ``mo.occsb[1] = 0.3`` will not. """ if self.kind == "generalized": raise NotImplementedError if self.occs is None: return None if self.kind == "restricted": if self.occs_aminusb is None: # heuristics ... if (self.occs == self.occs.astype(int)).all(): # restricted open-shell HF/KS result = self.occs - np.clip(self.occs, 0, 1) else: # restricted closed-shell natural orbitals result = self.occs / 2 else: result = (self.occs - self.occs_aminusb) / 2 result.flags.writeable = False return result return self.occs[self.norba :] @occsb.setter def occsb(self, occsb): if self.kind == "generalized": raise NotImplementedError if self.kind == "restricted": occsb = np.array(occsb) if self.occs is None: self.occs = occsb self.occs_aminusb = -occsb else: occsa = np.array(self.occsa) self.occs = occsa + occsb self.occs_aminusb = occsa - occsb else: self.occs[self.norba :] = occsb @property def coeffsa(self): """Return alpha orbital coefficients.""" if self.kind == "generalized": raise NotImplementedError if self.coeffs is None: return None if self.kind == "restricted": return self.coeffs return self.coeffs[:, : self.norba] @property def coeffsb(self): """Return beta orbital coefficients.""" if self.kind == "generalized": raise NotImplementedError if self.coeffs is None: return None if self.kind == "restricted": return self.coeffs return self.coeffs[:, self.norba :] @property def energiesa(self): """Return alpha orbital energies.""" if self.kind == "generalized": raise NotImplementedError if self.energies is None: return None if self.kind == "restricted": return self.energies return self.energies[: self.norba] @property def energiesb(self): """Return beta orbital energies.""" if self.kind == "generalized": raise NotImplementedError if self.energies is None: return None if self.kind == "restricted": return self.energies return self.energies[self.norba :] @property def irrepsa(self): """Return alpha irreps.""" if self.kind == "generalized": raise NotImplementedError if self.irreps is None: return None if self.kind == "restricted": return self.irreps return self.irreps[: self.norba] @property def irrepsb(self): """Return beta irreps.""" if self.kind == "generalized": raise NotImplementedError if self.irreps is None: return None if self.kind == "restricted": return self.irreps return self.irreps[self.norba :]