3.2. The Occupation Model in PyBEST

In PyBEST, the number of occupied orbitals (for the restricted and unrestricted case) is specified using various OccupationModel classes. In the current version, the number of occupied orbitals in atoms and molecules (doubly and singly occupied orbitals) is assigned using a Basis instance. The charge and/or the number of unpaired electrons can be specified separately. By default PyBEST always assumes a closed-shell system with either doubly occupied or unoccupied orbitals (restricted orbitals; even number of electrons) or one unpaired electron (unrestricted orbitals; odd number of electrons), which will be placed in an alpha orbital. Specifically, there are five different occupation models implemented. Each occupation model is explained in detail below.

Alternatively, an occupation number instance can be created from a LinalgFactory instance. In this case, we need to specify the total number of electrons or occupation pattern explicitly. This feature has to be chosen for model Hamiltonians or for external Hamiltonians (that is, one- and two-electron integrals) not generated in PyBEST. See Example Python scripts for examples.

The input structure for all supported occupation models is similar. To create an instance, do as follows if a Basis instance gobasis is present,

# Assuming a gobasis instance ``gobasis``
# Initialize some occupation model ``OccModel`` using default settings
occ_model = OccModel(gobasis)

while, for a LinalgFactory instance lf, the total number of electrons needs to be specified using the keyword argument nel,

# Assuming an lf instance ``lf`` with 10 electrons
# Initialize some occupation model ``OccModel`` using default settings
occ_model = OccModel(lf, nel=10)

Note that OccModel is used as an example and encodes one of the supported occupation model classes of PyBEST mentioned below. It does not represent any class implementation in PyBEST.

All supported occupation models share some common features and functionality. For instance, they can be defined for restricted and unrestricted orbitals. The table below summarizes the most important arguments and keyword arguments supported in each class.

Occupation model class	Arguments		Keyword Arguments
	gobasis	lf	charge	alpha	unrestricted	ncore	nel
AufbauOccModel	✅	✅	✅	✅	✅	✅	✅
FixedOccModel	✅	✅	✅		✅	✅
FermiOccModel	✅	✅	✅	✅	✅	✅	✅
FractionalOccModel	✅	✅	✅		✅	✅
AufbauSpinOccModel	✅	✅	✅		fixed	✅	✅

The keyword argument nel has to be specified if an lf instance is passed to create an occupation number instance. It can only be omitted if the occupation model uses different keyword arguments to specify the occupation pattern of orbitals (more details on how to specify occupation patterns are given in each occupation model subsection). Below, each argument and keyword argument is explained.

gobasis:

(Basis) an instance of Basis. It contains the information on the total number of electrons

lf:

(Linalg) an instance of LinalgFactory. It does NOT contain any information on the total number of electrons. Thus, nel (or any other occupation pattern) needs to be specified.

charge:

(int) the total charge of the system

alpha:

(int) number of unpaired electrons

unrestricted:

(boolean) enforce unrestricted occupiations for both even and odd numbers of electrons

ncore:

(int) the number of frozen core orbitals. The same value for alpha and beta orbitals is assumed. They are not allowed to differ.

nel:

(int) the total number of electrons. Only mandatory if an lf instance is passed and for selected occupation models (if occupations are not specified otherwise).

Note

Currently, post-HF calculations only support the Aufbau occupation model. All other occupation models are only supported in SCF calculations (UHF and RHF).

The occupation module not only stores information on the occupation pattern of the orbitals but also all information concerning active orbital spaces. The current version of PyBEST only supports frozen core orbitals. By defining a frozen core using some occupation model, all attributes are set accordingly. This includes the number of active orbitals nact, the number of active occupied orbitals nacto, and the number of active virtual orbitals nactv. All attributes are stored in an instance of some occupation model class. To specify a frozen core, the ncore keyword argument has to be set during initialization,

# Assuming a gobasis instance ``gobasis``
# Initialize some occupation model ``OccModel`` with three frozen core orbitals
occ_model = OccModel(gobasis, ncore=3)

# Assuming an lf instance ``lf`` with 10 electrons
# Initialize some occupation model ``OccModel`` with three frozen core orbitals
occ_model = OccModel(lf, nel=10, ncore=3)

3.2.1. The Aufbau occupation model

In most electronic structure calculations, this occupation model will be first and most decent choice. The orbitals are occupied according to the Aufbau principle, that is, the energetically lowest-lying orbitals become occupied. To initialize the occupation model, create an instance of the AufbauOccModel class. The number of doubly- or singly-occupied orbitals is determined from a Basis instance as follows

• Even number of electrons: restricted orbitals

• Odd number of electrons: unrestricted orbitals with 1 electron in an alpha orbital

# Assuming a gobasis instance ``gobasis``
# Initialize an instance of the AufbauOccModel using default settings
occ_model = AufbauOccModel(gobasis)

The default assignment of doubly occupied and singly occupied orbitals can be overwritten using the following keyword arguments

charge:

(int) the total charge of the system

alpha:

(int) number of unpaired electrons

unrestricted:

(boolean) enforce unrestricted occupations even for an even number of electrons

The examples listed below (see Example Python scripts) demonstrate how each keyword argument is used to specify the occupation pattern of a molecule or a model Hamiltonian.

Alternatively, an instance of the AufbauOccModel can be created from a LinalgFactory instance. To do so, the Basis instance gobasis has to be replaced by a LinalgFactory instance lf. In addition, we have to specify the total number of electrons otherwise, the initialization of the AufbauOccModel will raise an error.

# Defining some LinalgFactory instance with ten orbitals
lf = DenseLinalgFactory(10)
# Initialize an instance of the AufbauOccModel using an lf instance
# The total number of electrons is specified using nel
# We assume four electrons
occ_model = AufbauOccModel(lf, nel=4)

Each OccupationModel offers the feature to specify the number of frozen core orbitals which are used in post-HF modules.

ncore:

(int) the number of frozen core orbitals. The same value for alpha and beta orbitals is assumed. They are not allowed to differ.

# Assuming a gobasis instance ``gobasis``
# Initialize an instance of the AufbauOccModel using default settings
# Enforce four frozen core orbitals
occ_model = AufbauOccModel(gobasis, ncore=4)

Note

A frozen core is only supported in post-HF methods. It has no effect in HF calculations and will be simply ignored.

3.2.2. The fixed occupation model

The FixedOccModel freezes the occupation numbers to some pre-defined value. It features similar functionality as the AufbauOccModel using similar keyword arguments to steer the occupation of the orbitals.

charge:

(int) the total charge of the system

unrestricted:

(boolean) enforce unrestricted occupations even for an even number of electrons

ncore:

(int) the number of frozen core orbitals. The same number of alpha and beta orbitals is frozen. They are not allowed to differ.

The occupations have to be passed as a one-dimensional numpy array, where occ_a assigns an occupation vector to alpha orbitals, while occ_b (optional) specifies the occupation vector of beta orbitals. If only occ_a is given, PyBEST assumes occ_b = occ_a and restricted orbitals.

# Assume a gobasis instance ``gobasis`` with eight electrons
# Restricted case
occ_model = FixedOccModel(gobasis, occ_a=np.array([1.0, 1.0, 1.0, 0.5, 0.5, 0.0]))

# Unrestricted cases
# Same occupation number vector for alpha and beta
occ_model = FixedOccModel(
    gobasis, unrestricted=True, occ_a=np.array([1.0, 1.0, 1.0, 0.5, 0.5, 0.0])
)
# Different occupation number vector for alpha and beta
# Defaults to unrestricted orbitals
occ_model = FixedOccModel(
    gobasis,
    occ_a=np.array([1.0, 1.0, 1.0, 0.5, 0.5, 0.0]),
    occ_b=np.array([1.0, 0.7, 1.0, 1.0, 0.0, 0.3])
)

More examples are listed below (see Example Python scripts).

Note

The number of electrons stored in gobasis has to agree with the number of electrons in the occupation number vectors occ_a and occ_b. Otherwise, an error will be raised.

3.2.3. The Fermi occupation model

In case of convergence difficulties in HF calculations, the Fermi-smearing method can be applied to fill up the orbitals [rabuck1999].

# Evoke an instance of the FermiOccModel using default parameters
# (default occupations, T=250K, dT=50K, method=pfon - pseudo fractional occupation numbers)
occ_model = FermiOccModel(gobasis)

In general, the FermiOccModel is a child class of the AufbauOccModel. Thus, both use the same set of keyword arguments charge, unrestricted, alpha, and ncore. The convergence of the Fermi-smearing method can be steered using the following keyword arguments,

temperature:

(float) controls the width of the distribution (derivative). Default value: 250 K.

eps:

(float) the error on the sum of the occupation number when searching for the right Fermi level. Only required for the FON method. Default value: 1e-8.

method:

(str) the method to define fractional occupation numbers FON as presented in [rabuck1999]. FON ("fon") and pFON ("pfon") are supported. Default value: "pfon"

delta_t:

(float) the amount in K by which the temperature is reduced. Default value: 50 K.

More examples of how to specify different occupation models using Fermi-smearing are listed below (see Example Python scripts). Note that at the end of an SCF calculation that exploited Fermi smearing, a conventional calculation using the Aufbau occupation model should be performed.

3.2.4. The fractional occupation model

Fractional occupation numbers can be defined using the FractionalOccModel class. This occupation model is a child class of the AufbauOccModel class. Thus, the orbitals are occupied in accordance with the Aufbau principle, where the highest- lying orbitals might have fractional occupation numbers. In contrast, the AufbauOccModel allows for integer occupations only. These fractional occupation numbers are set during the initialization of the class and are not updated afterward.

Note

The sum of all occupation numbers has to equal the total number of electrons, otherwise an error will be raised.

The FractionalOccModel should NOT be used in combination with post-Hartree-Fock wave functions as the number of occupied orbitals might be assigned incorrectly. This occupation model might be useful for single Slater determinant wave functions in case of convergence difficulties.

Since the occupation model is a child class of AufbauOccModel, it contains similar functionality and supports similar keyword arguments to steer the occupation of the orbitals.

charge:

(int) the total charge of the system

unrestricted:

(boolean) enforce unrestricted occupations even for an even number of electrons

ncore:

(int) the number of frozen core orbitals. The same number of alpha and beta orbitals is frozen. They are not allowed to differ.

The fractional occupation numbers are set using the keyword arguments nocc_a and nocc_b, which specify the total number of (occupied) alpha and beta orbitals. Thus, their sum has to equal the total number of electrons.

nocc_a:

(float) the number of alpha electrons

nocc_b:

(float) the number of beta electrons (optional). If not provided, we assume nocc_b = nocc_a.

# Take gobasis information from ``gobasis`` (assuming 10 electrons)
# Enforce 5.5 alpha and 4.5 beta electrons
# Defaults to unrestricted orbitals
occ_model = FractionalOccModel(gobasis, nocc_a=5.5, nocc_b=4.5)

More examples of how to specify different fractional occupation patterns are listed below (see Example Python scripts).

3.2.5. The spin Aufbau occupation model

An Aufbau occupation model for spin orbitals can be defined using the AufbauSpinOccModel class. This occupation model is a child class of the AufbauOccModel class. Thus, the orbitals are occupied in accordance with the Aufbau principle. The AufbauSpinOccModel always enforces unrestricted orbitals. In contrast to the occupation models mentioned above, which assign “fixed” occupations to orbitals resulting in a well-defined number of occupied alpha and beta orbitals, this occupation model assigns the occupations of alpha and beta orbitals on-the-fly. Meaning that during an SCF optimization step, the energetically lowest-lying orbitals are occupied, independent of their spin degree of freedom. Therefore, the number of alpha and beta electrons may change during the optimization.

In contrast to the AufbauOccModel, only a subset of keyword arguments is supported.

charge:

(int) the total charge of the system

ncore:

(int) the number of frozen core orbitals. The same number of alpha and beta orbitals is frozen. They are not allowed to differ.

To create an instance of the AufbauSpinOccModel class, a similar syntax as for the AufbauOccModel is used,

# Take gobasis information from ``gobasis``
# (contains information on the number of electrons)
# Defaults to unrestricted orbitals
occ_model = AufbauSpinOccModel(gobasis)

More examples of how to specify different Aufbau occupation patterns for spin orbitals are listed below (see Example Python scripts).

3.2.6. Example Python scripts

Several complete examples can be found in the directory data/examples/occ_model.

3.2.6.1. AufbauOccModel

This is a basic example of how to create an instance of the AufbauOccModel using various settings.

Listing 3.2 data/data/examples/occ_model/aufbau_model.py

from pybest import context
from pybest.gbasis import get_gobasis
from pybest.linalg import DenseLinalgFactory
from pybest.occ_model import AufbauOccModel

# Aufbau occupation model
# -----------------------

# The H2O molecule
# ----------------

# Load the coordinates from file for H2O.
# Use the XYZ file from PyBEST's test data directory.
fn_xyz = context.get_fn("test/water.xyz")

# Create a Gaussian basis set
factory = get_gobasis("cc-pvdz", fn_xyz)

# Use Aufbau occupation model to occupy orbitals

# If we do not specify the number of frozen core orbitals (ncore),
# then ncore will be calculated automatically

# Example 1: H2O with only doubly occupied orbitals
# defaults to restricted orbitals
occ_model = AufbauOccModel(factory, ncore=0)

# Example 2: H2O+ with one singly occupied (alpha) orbital
# defaults to unrestricted orbitals with one unpaired alpha electron
occ_model = AufbauOccModel(factory, charge=1, ncore=0)

# Example 3: H2O with five alpha and five beta electrons
# enforce unrestricted orbitals
occ_model = AufbauOccModel(factory, unrestricted=True, ncore=0)

# Example 4: H2O+ with three singly occupied (alpha) orbitals
# defaults to unrestricted orbitals
occ_model = AufbauOccModel(factory, charge=1, alpha=3, ncore=0)


# The NO molecule
# ---------------

# Load the coordinates from file for NO.
# Use the XYZ file from PyBEST's test data directory.
fn_xyz = context.get_fn("test/no.xyz")

# Create a Gaussian basis set
factory = get_gobasis("cc-pvdz", fn_xyz)

# Use Aufbau occupation model to occupy orbitals

# Example 1: NO molecule
# defaults to unrestricted orbitals with one unpaired electron
occ_model = AufbauOccModel(factory, ncore=0)


# Using an LF instance: here, we have to pass the number of electrons
# -------------------------------------------------------------------

# Create an LF instance for 24 orbitals and 10 electrons
lf = DenseLinalgFactory(24)

# Use Aufbau occupation model to occupy orbitals

# Select default parameters
# restricted orbitals with 10 electrons
occ_model = AufbauOccModel(lf, nel=10, ncore=0)

# enforce unrestricted representation
occ_model = AufbauOccModel(lf, nel=10, unrestricted=True, ncore=0)

# unrestricted orbitals with 9 electrons (with 1 unpaired electron)
# defaults to unrestricted orbitals
occ_model = AufbauOccModel(lf, nel=9, ncore=0)

# unrestricted orbitals with 9 electrons (with 3 unpaired electrons)
# defaults to unrestricted orbitals
occ_model = AufbauOccModel(lf, nel=9, alpha=3, ncore=0)

3.2.6.2. FixedOccModel

This is a basic example on how to create an instance of the FixedOccModel using various settings.

Listing 3.3 data/data/examples/occ_model/fixed_model.py

import numpy as np

from pybest import context
from pybest.gbasis import get_gobasis
from pybest.linalg import DenseLinalgFactory
from pybest.occ_model import FixedOccModel

# Occupation model with fixed (user-defined) occupations
# ------------------------------------------------------

# In the FixedOccModel occupation model, the alpha and beta orbitals are
# occupied with respect to user-defined occupation number vectors.
# The occupation numbers have to be any real number in the interval [0,1].
# Only the occupation numbers for all occupied orbitals have to be specified.
# All remaining (virtual) orbitals will be assigned occupation numbers of 0.

# The H2O molecule
# ----------------

# Load the coordinates from file for H2O.
# Use the XYZ file from PyBEST's test data directory.
fn_xyz = context.get_fn("test/water.xyz")

# Create a Gaussian basis set
factory = get_gobasis("cc-pvdz", fn_xyz)

# Use FixedOccModel occupation model

# Example 1: H2O with a total number of 10 electrons and occupying the 5
# energetically lowest orbitals (this is equivalent to the AufbauOccModel)
# defaults to restricted orbitals
occ_model = FixedOccModel(factory, occ_a=np.array([1, 1, 1, 1, 1]))

# Example 2: H2O+ with a 5 alpha and 4 beta electrons (Aufbau occupation)
# defaults to unrestricted orbitals
occ_model = FixedOccModel(
    factory,
    charge=1,
    occ_a=np.array([1, 1, 1, 1, 1]),
    occ_b=np.array([1, 1, 1, 1]),
)

# Example 3: H2O+ with a 5 alpha and 4 beta electrons (non Aufbau occupation)
# defaults to unrestricted orbitals
occ_model = FixedOccModel(
    factory,
    charge=1,
    occ_a=np.array([1, 1, 0, 0, 1, 1, 1]),
    occ_b=np.array([1, 1, 1, 0, 1]),
)

# Example 4: H2O with fractional occupation numbers in some orbitals (non Aufbau)
# defaults to restricted orbitals (enforces fractional occupations on alpha and
# beta orbitals)
occ_model = FixedOccModel(factory, occ_a=np.array([1, 1, 0.3, 0, 0.7, 1, 1]))

# Example 5: H2O with fractional occupation numbers in some orbitals (non Aufbau)
# enforces unrestricted orbitals (enforces fractional occupations on alpha and
# beta orbitals)
occ_model = FixedOccModel(
    factory, unrestricted=True, occ_a=np.array([1, 1, 0.3, 0, 0.7, 1, 1])
)

# The NO molecule
# ---------------

# Load the coordinates from file for NO.
# Use the XYZ file from PyBEST's test data directory.
fn_xyz = context.get_fn("test/no.xyz")

# Create a Gaussian basis set
factory = get_gobasis("cc-pvdz", fn_xyz)

# Use FixedOccModel occupation model

# Example 1: NO molecule with a total number of 15 electrons (since we have an
# odd number of electrons, both occ_a and occ_b need to be specified).
# defaults to unrestricted orbitals
occ_model = FixedOccModel(
    factory,
    occ_a=np.array([1, 1, 1, 1, 1, 1, 1, 0, 1]),
    occ_b=np.array([1, 1, 1, 1, 1, 1, 1]),
)


# Using an LF instance: since we have to pass occ_a, we do not need to pass nel
# -----------------------------------------------------------------------------

# Create an LF instance for 24 orbitals and 10 electrons
lf = DenseLinalgFactory(24)

# Use FixedOccModel occupation model

# restricted orbitals with 10 electrons and fractional occupation numbers
# defaults to restricted orbitals
occ_model = FixedOccModel(lf, occ_a=np.array([1, 1, 0.3, 0, 0.7, 1, 1]))

# unrestricted orbitals with 15 electrons (8 alpha and 7 beta)
# defaults to unrestricted orbitals
occ_model = FixedOccModel(
    lf,
    occ_a=np.array([1, 1, 1, 1, 1, 1, 1, 0, 1]),
    occ_b=np.array([1, 1, 1, 1, 1, 1, 1]),
)

3.2.6.3. FermiOccModel

This is a basic example on how to create an instance of the FermiOccModel using various settings.

Listing 3.4 data/data/examples/occ_model/aufbau_fermi.py

from pybest import context
from pybest.gbasis import get_gobasis
from pybest.linalg import DenseLinalgFactory
from pybest.occ_model import FermiOccModel

# Aufbau occupation model with Fermi smearing
# -------------------------------------------

# In the FermiOccModel occupation model, the alpha and beta orbitals are
# occupied with respect to their energy. That is, the energetically lowest
# (alpha and beta) orbitals are occupied first, while a Fermi distribution is
# applied to the occupation numbers. Thus, the occupation numbers are real
# numbers in the range [0, 1] for each spin channel. During the SCF
# optimization the temperature of the Fermi distribution is (by default)
# gradually lowered.


# The H2O molecule
# ----------------

# Load the coordinates from file for H2O.
# Use the XYZ file from PyBEST's test data directory.
fn_xyz = context.get_fn("test/water.xyz")

# Create a Gaussian basis set
factory = get_gobasis("cc-pvdz", fn_xyz)

# Use Aufbau occupation model to occupy orbitals, while occupation numbers are
# assigned according to a Fermi distribution

# Example 1: H2O with default parameters (T=250K, dT=50K, method=pfon - pseudo
# fractional occupation numbers)
# defaults to restricted orbitals
occ_model = FermiOccModel(factory)

# Example 2: H2O+ with one (in total) alpha electron and default parameters
# (T=250K, dT=50K, method=pfon - pseudo fractional occupation numbers)
# defaults to unrestricted orbitals
occ_model = FermiOccModel(factory, charge=1)

# Example 3: H2O with default parameters (T=250K, dT=50K, method=pfon - pseudo
# fractional occupation numbers)
# enforce unrestricted orbitals
occ_model = FermiOccModel(factory, unrestricted=True)

# Example 4: H2O+ with (in total) three unpaired alpha electrons and default
# parameters (T=250K, dT=50K, method=pfon - pseudo fractional occupation numbers)
# defaults to unrestricted orbitals
occ_model = FermiOccModel(factory, charge=1, alpha=3)

# Example 5: H2O+ with (in total) three unpaired alpha electrons and default
# parameters (T=500K, dT=25K, method=fon - fractional occupation numbers)
# defaults to unrestricted orbitals
occ_model = FermiOccModel(
    factory, charge=1, alpha=3, temperature=500, delta_t=25, method="fon"
)

# The NO molecule
# ---------------

# Load the coordinates from file for NO.
# Use the XYZ file from PyBEST's test data directory.
fn_xyz = context.get_fn("test/no.xyz")

# Create a Gaussian basis set
factory = get_gobasis("cc-pvdz", fn_xyz)

# Use Aufbau occupation model to occupy orbitals, while occupation numbers are
# assigned according to a Fermi distribution

# Example 1: NO molecule
# defaults to unrestricted orbitals with one unpaired electron (in total)
occ_model = FermiOccModel(factory)

# Example 2: NO molecule
# parameters (T=500K, dT=25K, method=fon - fractional occupation numbers)
# defaults to unrestricted orbitals with one unpaired electron (in total)
occ_model = FermiOccModel(factory, temperature=500, delta_t=25, method="fon")

# Example 3: NO molecule with 3 unpaired electrons
# parameters (T=500K, dT=25K, method=fon - fractional occupation numbers)
# defaults to unrestricted orbitals with one unpaired electron (in total)
occ_model = FermiOccModel(
    factory, alpha=3, temperature=500, delta_t=25, method="fon"
)

# Using an LF instance: here, we have to pass the number of electrons
# -------------------------------------------------------------------

# Create an LF instance for 24 orbitals and 10 electrons
lf = DenseLinalgFactory(24)

# Use Aufbau occupation model to occupy orbitals

# Select default parameters
# restricted orbitals with 10 electrons
occ_model = FermiOccModel(lf, nel=10)

# enforce unrestricted representation
occ_model = FermiOccModel(lf, nel=10, unrestricted=True)

# unrestricted orbitals with 9 electrons (with 1 unpaired electron)
# defaults to unrestricted orbitals
occ_model = FermiOccModel(lf, nel=9)

# unrestricted orbitals with 9 electrons (with 3 unpaired electrons in total)
# defaults to unrestricted orbitals
occ_model = FermiOccModel(lf, nel=9, alpha=3)

3.2.6.4. FractionalOccModel

This is a basic example on how to create an instance of the FractionalOccModel using various settings.

Listing 3.5 data/data/examples/occ_model/aufbau_fractional_model.py

from pybest import context
from pybest.gbasis import get_gobasis
from pybest.linalg import DenseLinalgFactory
from pybest.occ_model import FractionalOccModel

# Aufbau occupation model allowing for fractional occupation numbers in the
# highest occupied molecular orbital
# -------------------------------------------------------------------------

# In the FractionalOccModel occupation model, the HOMO of the alpha and beta
# electrons can be occupied with a fraction of an electron. During the
# SCF optimization, this fractional occupation numbers are kept fixed.

# The H2O molecule
# ----------------

# Load the coordinates from file for H2O.
# Use the XYZ file from PyBEST's test data directory.
fn_xyz = context.get_fn("test/water.xyz")

# Create a Gaussian basis set
factory = get_gobasis("cc-pvdz", fn_xyz)

# Use Aufbau occupation model for fractional occupation numbers

# Example 1: H2O with 5.5 alpha and 4.5 beta electrons (sums up to 10 electrons)
# defaults to unrestricted orbitals
occ_model = FractionalOccModel(factory, nocc_a=5.5, nocc_b=4.5)

# Example 2: H2O+ with a 4.5 alpha and beta electrons (sums up 9 electrons)
# defaults to restricted orbitals
occ_model = FractionalOccModel(factory, charge=1, nocc_a=4.5)

# Example 2: H2O+ with a 4.5 alpha and beta electrons (sums up 9 electrons)
# enforce unrestricted orbitals
occ_model = FractionalOccModel(
    factory, charge=1, nocc_a=4.5, unrestricted=True
)


# The NO molecule
# ---------------

# Load the coordinates from file for NO.
# Use the XYZ file from PyBEST's test data directory.
fn_xyz = context.get_fn("test/no.xyz")

# Create a Gaussian basis set
factory = get_gobasis("cc-pvdz", fn_xyz)

# Use Aufbau occupation model for fractional occupation numbers

# Example 1: NO with a 7.5 alpha and beta electrons (sums up 15 electrons)
# defaults to restricted orbitals
occ_model = FractionalOccModel(factory, nocc_a=7.5)

# Example 1: NO with a 7.5 alpha and beta electrons (sums up 15 electrons)
# enforce unrestricted orbitals
occ_model = FractionalOccModel(factory, nocc_a=7.5, unrestricted=True)


# Using an LF instance: since we have to pass nocc_a, we do not need to pass nel
# ------------------------------------------------------------------------------

# Create an LF instance for 24 orbitals and 10 electrons
lf = DenseLinalgFactory(24)

# Use Aufbau occupation model for fractional occupation numbers

# Select default parameters
# restricted orbitals with 15 electrons (7.5 alpha and 7.5 beta)
# defaults to restricted orbitals
occ_model = FractionalOccModel(lf, nocc_a=7.5)

# enforce unrestricted representation
occ_model = FractionalOccModel(lf, nocc_a=7.5, unrestricted=True)

# unrestricted orbitals with 15 electrons (9.1 alpha and 5.9 beta)
# defaults to unrestricted orbitals
occ_model = FractionalOccModel(lf, nocc_a=9.1, nocc_b=5.9)

# For an LF instance, we can have a fractional number of electrons
# unrestricted orbitals with 15.1 electrons (9.2 alpha and 5.9 beta)
# defaults to unrestricted orbitals
occ_model = FractionalOccModel(lf, nocc_a=9.2, nocc_b=5.9)

3.2.6.5. AufbauSpinOccModel

This is a basic example on how to create an instance of the AufbauSpinOccModel using various settings.

Listing 3.6 data/data/examples/occ_model/aufbau_spin_model.py

from pybest import context
from pybest.gbasis import get_gobasis
from pybest.linalg import DenseLinalgFactory
from pybest.occ_model import AufbauSpinOccModel

# Aufbau occupation model for unrestricted wave functions
# -------------------------------------------------------

# In the AufbauSpinOccModel occupation model, the alpha and beta orbitals are
# occupied with respect to their energy. That is, the energetically lowest
# (alpha and beta) orbitals are occupied first. Only the total number of
# electrons is fixed, while the number of electrons in each spin channel can
# change during the SCF optimization.

# The H2O molecule
# ----------------

# Load the coordinates from file for H2O.
# Use the XYZ file from PyBEST's test data directory.
fn_xyz = context.get_fn("test/water.xyz")

# Create a Gaussian basis set
factory = get_gobasis("cc-pvdz", fn_xyz)

# Use Aufbau occupation model for unrestricted wave functions

# Example 1: H2O with a total number of 10 electrons
# defaults to unrestricted orbitals
occ_model = AufbauSpinOccModel(factory)

# Example 2: H2O+ with a total number of 9 electrons
# defaults to unrestricted orbitals
occ_model = AufbauSpinOccModel(factory, charge=1)


# The NO molecule
# ---------------

# Load the coordinates from file for NO.
# Use the XYZ file from PyBEST's test data directory.
fn_xyz = context.get_fn("test/no.xyz")

# Create a Gaussian basis set
factory = get_gobasis("cc-pvdz", fn_xyz)

# Use Aufbau occupation model for unrestricted wave functions

# Example 1: NO molecule with a total number of 15 electrons
# defaults to unrestricted orbitals
occ_model = AufbauSpinOccModel(factory)


# Using an LF instance: we need to pass the number of electrons nel
# -----------------------------------------------------------------

# Create an LF instance for 24 orbitals and 10 electrons
lf = DenseLinalgFactory(24)

# Use Aufbau occupation model for unrestricted wave functions

# Select default parameters
# 15 electrons
# defaults to unrestricted orbitals
occ_model = AufbauSpinOccModel(lf, nel=15)