Source code for moldesign.orbitals.orbitals

# Copyright 2016 Autodesk Inc.
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# See the License for the specific language governing permissions and
# limitations under the License.

""" Class definitions for atomic and molecular orbitals.

    In this documentation, we use the following conventions for labeling orbitals:
      - atomic orbitals using lower case greek labels and subscripts, e.g.,
          :math:`\left| \mu \right \rangle, F_{\nu \lambda}, etc.
      - molecular orbitals use lower case labels and subscripts, e.g.,
          :math:`\left| i \right \rangle, F_{kl}, etc.
      - adiabatic electronic states are indexed using capital letters, _N_, _L_, _M_, etc.

import numpy as np

from moldesign import units as u
from moldesign.utils import Alias

SHELLS = {0: 's', 1: 'p', 2: 'd', 3: 'f', 4: 'g', 5: 'h'}
ANGMOM = {v: k for k, v in SHELLS.iteritems()}

# See
SPHERICALNAMES = {(0, 0): 's', (1, -1): 'p(x)', (1, 0): 'p(z)', (1, 1): 'p(y)',
                  (2, 0): 'd(z^2)', (2, -2): 'd(xy)', (2, -1): 'd(yz)',
                  (2, 1): 'd(xz)', (2, 2): 'd(x^2-y^2)',
                  (3, -3): 'f(3yx^2-y^3)', (3, -2): 'f(xyz)', (3, -1): 'f(yz^2)',
                  (3, 0): 'f(z^3)', (3, 1): 'f(xz^2)', (3, 2): 'f(zx^2-zy^2)',
                  (3, 3): 'f(x^3-3xy^2)'}
ANGULAR_NAME_TO_COMPONENT = {'': (0, 0), 'x': (1, -1), 'y': (1, 1), 'z': (1, 0),
                             'z^2': (2, 0), 'xy': (2, -2), 'yz': (2, -1), 'xz': (2, 1),
                             'x^2-y^2': (2, 2),
                             'zx^2-zy^2': (3, 2), 'xyz': (3, -2), 'z^3': (3, 0),
                             '3yx^2-y^3': (3, -3), 'x^3 - 3xy^2': (3, 3), 'xz^2': (3, 1),
                             'yz^2': (3, -1)}

[docs]class Orbital(object): r""" Stores a single orbital and its meta-data Generally wants to be part of a set of MolecularOrbitals The orbital is defined as .. math:: \left| i \right \rangle = \sum_\mu c_{i \mu} \left| \mu \right \rangle where the coefficients :math:`c_{i \mu}` are stored in ``self.coeffs`` and the basis orbitals :math:`\left| \mu \right \rangle` are stored at ``self.basis`` """ def __init__(self, coeffs, basis=None, wfn=None, occupation=None, name='unnamed'): """ Initialization: Args: coeffs (numpy.array): orbital coefficients basis (moldesign.orbitals.basis.BasisSet): basis for this orbital wfn (moldesign.orbitals.wfn.ElectronicWfn): total electronic wavefunction that this orbital is a part of occupation (float): occupation of this orbital name (str): optional label for this orbital: """ self.coeffs = np.array(coeffs) = name self.molecule = None self.basis = basis self.occupation = occupation self.wfn = wfn self.index = None # will be set by the containing MolecularOrbitals object # Assign the basis functions if wfn is not None and wfn.aobasis is not None: if self.basis is not None: assert wfn.aobasis is self.basis else: self.basis = wfn.aobasis if self.basis is not None: assert len(self.coeffs) == len(self.basis)
[docs] def overlap(self, other): """ Calculate overlap with another orbital Args: other (Orbital): calculate the overlap with this orbital Returns: float: orbital overlap """ return
[docs] def fock_element(self, other): """ Calculate fock matrix element with another orbital Args: other (Orbital): calculate the fock element with this orbital Returns: u.Scalar[energy]: fock matrix element """ return # uses ldot to preserve units
@property def energy(self): """ u.Scalar[energy]: This orbital's energy Note: This is equivalent to self.fock(self) """ return self.fock_element(self) def __call__(self, coords): """ Calculate the orbital's value at a position (or list of positions) Args: coords (u.Vector[length]): Coordinates (shape ``(len(coords), 3)``) Returns: u.Scalar[length**(-3/2)]: value of the orbital """ return self.basis(coords, coeffs=self.coeffs) def __repr__(self): return '<%s %s>' % (self.__class__.__name__, __len__ = Alias('coeffs.__len__')
[docs]class MolecularOrbitals(object): """ Stores a wfn of molecular orbitals in an AO wfn Orbitals are accessed as orbs[orbital index, ao index] """ def __init__(self, orbitals, wfn=None, basis=None, canonical=False, orbtype=None): """ Initialization: Args: orbitals (List[Orbital] OR numpy.array): EITHER a list of orbitals OR an array of coefficients (indexed as coeffs[orbital_index, basis_index]) basis (moldesign.orbitals.basis.BasisSet): orbital basis set wfn (moldesign.orbitals.wfn.ElectronicWfn): total electronic wavefunction that this orbital is a part of canonical (bool): designate these as the canonical orbitals orbtype (str): description of these orbitals """ # Determine if these are the canonical orbitals if canonical: assert orbtype is None or orbtype == 'canonical' orbtype = 'canonical' if orbtype == 'canonical': canonical = True if not hasattr(orbitals[0], 'basis'): coeffs = orbitals orbitals = [] for c in coeffs: orbitals.append(Orbital(c, basis=basis, wfn=wfn)) self.orbitals = orbitals self.coeffs = np.array([orb.coeffs for orb in self.orbitals]) self.wfn = wfn self.basis = basis self.orbtype = orbtype if self.wfn is None: self.wfn = self.orbitals[0].wfn if self.basis is None: self.basis = self.orbitals[0].basis for iorb, orbital in enumerate(self.orbitals): orbital.index = iorb assert orbital.basis == self.basis assert orbital.wfn == self.wfn orbital.coeffs = self.coeffs[iorb, :] if canonical: self._set_cmo_names()
[docs] def align_phases(self, other, threshold=0.5, assert_same_type=True): """ Flip the signs of these orbitals to bring them into maximum coincidence with another set of orbitals Args: other (MolecularOrbitals): the "reference" set of orbitals to match phases with threshold (float): only flip orbital if the overlap is less than -1*threshold assert_same_type (bool): require that ``self.orbtype == other.orbtype`` Note: This function assumes that the overlap matrix is the same for both sets of orbitals - this is a reasonable assumption if the two sets of orbitals were calculated at very similar molecular geometries. """ if assert_same_type: assert self.orbtype == other.orbtype, "Orbital type mismatch: %s vs. %s" % ( self.orbtype, other.orbtype) for thisorb, otherorb in zip(self, other): if thisorb.overlap(otherorb) < -1.0 * threshold: thisorb.coeffs *= -1.0
# TODO: print a warning if overlap is small?
[docs] def overlap(self, other): """ Calculate overlaps between this and another set of orbitals Args: other (MolecularOrbitals): Returns: numpy.ndarray: overlaps between the two sets of orbitals Example: >>> canonical = mol.wfn.canonical >>> atomic = mol.wfn.basis >>> overlaps = canonical.overlap(atomic) >>> overlaps[i, j] == canonical.orbitals[i].overlap(atomic.orbitals[j]) True """ return
def __iter__(self): return iter(self.orbitals) def __len__(self): return len(self.orbitals) def __str__(self): return '%s orbitals' % self.orbtype def __repr__(self): return '<%d %s %s in %s>' % (len(self), self.orbtype, self.__class__.__name__, str(self.wfn)) def __getitem__(self, item): return self.orbitals[item] def _to_ao_density_matrix(self): c = self.coeffs * self.occupations[:, None]/2.0 return 2.0* @property def energies(self): """u.Vector[energy]: energies of the molecular orbitals This is just the diagonal of the fock matrix""" return self.fock.diagonal() @property def occupations(self): """ np.ndarray: orbital occupation numbers """ return np.array([orb.occupation for orb in self.orbitals]) @property def fock(self): """u.Array[energy]: Fock matrix for these orbitals""" return self.from_ao(self.wfn.fock_ao) @property def overlaps(self): """np.array: overlap matrix for these orbitals""" return self.from_ao(self.wfn.aobasis.overlaps) @property def h1e(self): """u.Array[energy]: 1-electron matrix elements for these orbitals""" return self.from_ao(self.wfn.aobasis.h1e) @property def h2e(self): """u.Array[energy]: 2-electron matrix elements for these orbitals""" return self.fock - self.h1e
[docs] def from_ao(self, ao_operator): """ Transform an operator into this orbital basis from the ao basis Given the matrix elements :math:`\hat O_{\mu \nu}` of an operator over AO basis indices :math:`\mu,\nu`, returns the operator's matrix elements :math:`\hat O_{ij}` over orbital indices :math:`i,j`: ..math:: \hat O_{ij} = \left \langle i \right| \hat O \left| j \right \rangle = \sum_{\mu \nu}C_{i \mu} O_{\mu \nu} C_{j \nu} where :math:`C_{i \mu}` is the expansion coefficient for AO basis function :math:`\mu` in molecular orbital _i_. Args: ao_operator (u.Array): matrix elements of the operator in the ao basis Returns: u.Array: matrix elements of the operator in this orbital basis Note: Assumes that this set of orbitals is orthogonal """ # Dot doesn't place nice with units, so we need to pass them explicitly ao_units = u.get_units(ao_operator) return * ao_units
[docs] def to_ao(self, mo_operator): """ Transform an operator from this orbital basis into the AO basis Given the matrix elements :math:`\hat O_{ij}` of an operator over orbital basis indices :math:`i,j`, returns the operator's matrix elements :math:`\hat O_{\mu \nu}` over orbital indices :math:`\mu, \nu`: ..math:: \hat O_{\mu \nu} = \left \langle \mu \right| \hat O \left| \nu \right \rangle = \sum_{i,j,\lambda,\kappa}S_{\mu \lambda} C_{i \lambda} O_{ij} C_{j \kappa} S_{\kappa \nu} where :math:`S_{\mu \nu} = \left \langle \mu | \nu \right \rangle` is the AO overlap matrix and :math:`C_{i \mu}` is the expansion coefficient for AO basis function :math:`\mu` in molecular orbital _i_. Args: mo_operator (u.Array): matrix elements of the operator in this orbital basis Returns: u.Array: matrix elements of the operator in the AO basis """ units = u.get_units(mo_operator) s = self.wfn.aobasis.overlaps o_ao = return o_ao * units
def _set_cmo_names(self): for i, orb in enumerate(self.orbitals): if != 'unnamed' and is not None: continue if i <= self.wfn.homo: if i < self.wfn.homo - 2: = 'cmo %d' % i elif i == self.wfn.homo: = 'HOMO' else: = 'HOMO-%d' % (self.wfn.homo - i) else: if i == self.wfn.lumo: = 'LUMO' elif i <= self.wfn.lumo + 2: = 'LUMO+%d' % (i - self.wfn.lumo) else: = 'virt cmo %d' % i