3.1. Working with Orbitals in PyBEST

The (molecular) orbitals in PyBEST contain information on the AO/MO coefficient matrix, the orbital occupation numbers, and the orbital energies. Internally, the orbitals are an instance of the DenseOrbital class. Specifically, the AO/MO coefficients, occupation numbers, and energies are attributes of the DenseOrbital class. In the following, some simple orbital manipulations are explained that might be useful for users.

3.1.1. Creating an empty set of orbitals

Prior to any wave function calculation, the orbitals need to be created as an instance of the DenseOrbital class

# Create an empty set of orbitals
orb_a = lf.create_orbital()

where lf is an instance of the LinalgFactory class, where all linear algebra and tensor contraction operations are defined.

3.1.2. Printing orbital information

The AO/MO coefficient matrix, the orbital occupation numbers, and the orbital energies can be accessed as follows (assuming orb to be an instance of DenseOrbital)

# Print AO/MO coefficients
print(orb.coeffs)

# Print occupations
print(orb.occupations)

# Print orbital energies
print(orb.energies)

The AO/MO coefficients can also be printed in a more readable format using the print_ao_mo_coeffs() function

# Print all AO/MO coefficients in a more readable format
print_ao_mo_coeffs(gobasis, orb)

where gobasis is an instance of the GOBasis class. If only a subset of orbitals is required, this can be achieved using optional arguments

# Print AO/MO coefficients in a more readable format for orbitals
# no. 4 to no. 31 (including)
print_ao_mo_coeffs(gobasis, orb, begin=4, end=31)

where begin and end are the first and last index, respectively, of the orbital to be printed. Note that the orbitals are indexed using Fortran indexing, that is, beginning with 1. The order of the AOs is according to the list of atoms in the geometry input. In the output, each column corresponds to one molecular orbital.

Note

In all electronic structure methods that optimize the orbitals, the current orbitals are always overwritten with the most recent ones.

To save a copy of the current orbitals during the execution of the PyBEST code, use the copy() method of DenseOrbital,

# Create a true copy of the current orbitals
orb_output = orb.copy()

3.1.3. Writing orbitals to a checkpoint file

In PyBEST, all electronic structure methods save their own copy of the molecular orbitals and dump them in their method-dependent checkpoint files. In case, you want to write the molecular orbitals to a separate file at any stage of the program (represented in PyBEST’s internal format), you can use the IOData container (see also Input/Output Operations: the IOData Container for more details on how to use the container),

# Dump orbitals for restart
mol = IOData(olp=olp, orb_a=orb)
mol.to_file('checkpoint.h5')

where the overlap matrix is stored in olp and the orbitals in orb_a. The molecular orbitals (for restricted orbitals or alpha spin) are always associated with orb_a, while the overlap matrix is associated with olp. In case of unrestricted orbitals, PyBEST uses the attribute name orb_b. If the user chooses a different naming convention, PyBEST might run into unexpected errors (for instance, during restarts). See also Input/Output Operations: the IOData Container for more information on how to use the IOData container and all default variable names.

Note

To use previously dumped orbitals for restart purposes, the overlap matrix has to be written to file as well.

3.1.4. Reading orbitals from a checkpoint file

Most electronic structure methods in PyBEST feature a restart option, which reads and projects some previously stored orbitals automatically. The documentation below is mentioned for reasons of completeness to allow the user to manually read in orbitals from disk.

3.1.4.1. PyBEST’s internal format

Similar to the recipe of dumping molecular orbitals to file, they can be read from disk using the IOData container functionality. First, some previously stored PyBEST checkpoint file that contains the molecular orbitals and overlap integrals has to be read in,

# Read some PyBEST checkpoint file and store the corresponding data  in 'old'
old = IOData.from_file('checkpoint.h5')

The molecular orbitals and overlap matrix can be accessed from old as old.orb_a and old.olp, respectively. If the molecular geometry of the restart orbitals differs from the current molecular geometry, the molecular orbitals have to be projected against the current molecular orbitals orb and current overlap matrix olp,

# Project orbitals read from file against current orbitals
project_orbitals(old.olp, olp, old.orb_a, orb)

The project_orbitals() function overwrites orb and olp with the projected restart orbitals and overlap matrix, respectively.

Note

If the restart orbitals are obtained for a different molecular geometry, the project_orbitals() function has to be used, otherwise all electronic structure calculations will be corrupted.

3.1.4.2. External file format

Orbital information can also be imported from external programs. Supported file formats are .mkl and .molden. Similar to the procedure mentioned above (PyBEST’s internal format), the IOData container is used to read in some external file. The file extension is handled automatically.

# Load molden file file
mol = IOData.from_file('water.molden')

# Print the total number of molecular orbitals
print(mol.orb_a.nfn)

Similarly, if the orbitals are not projected onto a new basis set, the Hamiltonian has to be recomputed for the external basis set stored as mol.gobasis.

The mol instance has the following attributes

mol.gobasis:

A Basis instance. It contains all basis set information.

mol.orb_a:

A DenseOrbital instance containing the molecular orbitals for alpha spin

mol.orb_b:

A DenseOrbital instance containing the molecular orbitals for beta spin (if present)

mol.coordinates:

A np.array containing the xyz coordinates

mol.atom:

A list of strings with the element names

3.1.5. Rotating orbitals

The molecular orbitals can be perturbed by various rotations. The user has the choice between (a series) of 2-by-2 rotations or a random rotation of all molecular orbitals.

3.1.5.1. Givens rotation (2-by-2 rotations)

Two molecular orbitals can be rotated using the rotate_2orbitals() method. Depending on the angle of rotation, the selected orbitals are either mixed or swapped. By default, the function rotates the HOMO and LUMO orbitals by 45 degrees.

# Mix HOMO and LUMO orbitals
orb.rotate_2orbitals()

If a different pair of orbitals is to be rotated, the corresponding orbital indices can be specified as optional arguments,

# Rotate 7th and 12th orbital by 60 deg; note zero-based indexing
orb.rotate_2orbitals(60*deg, 6, 11)

Keep in mind that PyBEST uses radians as unit for angles, that is, 60*deg == 60*np.pi/180. Since rotate_2orbitals() is a class method (defined in DenseOrbital), all quantities use Python-based indexing, starting with 0. The rotate_2orbitals() method can be applied in sequence if several orbitals are to be rotated with each other,

# Perform series of rotations: 7/12 by 60 degree, 8/14 by 20 degree, 9/26
# by 80 degree
orb.rotate_2orbitals(60*deg, 6, 11)
orb.rotate_2orbitals(20*deg, 7, 13)
orb.rotate_2orbitals(80*deg, 8, 25)

Note that an angle of 90 degrees is equivalent to swapping the corresponding orbitals. This functionality can also be achieved using the swap_orbitals() method (see Swapping orbitals below).

3.1.5.2. Random orbital rotations

All molecular orbitals can be perturbed simultaneously using the rotate_random() method. This function method performs a random unitary rotation to all molecular orbitals.

# perform random unitary rotation of orbitals orb
orb.rotate_random()

Note

All orbital rotation operations leave the occupation numbers and orbital energies unchanged. Only the MO coefficients are perturbed. The AOs remain always untouched.

3.1.6. Swapping orbitals

Two molecular orbitals can be swapped using the swap_orbitals() method. The list of orbital swaps is passed as a 2-dimensional numpy array, where each row corresponds to one transposition of orbitals. Since swap_orbitals() is again a class method, the indexing starts with 0.

# Swap some orbitals; note zero-based indexing
swaps = np.array([[0, 6], [5, 9]])
orb.swap_orbitals(swaps)

By default, the orbital occupations and orbital energies are swapped as well. In some cases (like in SCF calculations), the swapping of orbital occupations and/or orbital energies should be suppressed to provide a perturbed guess for the molecular orbitals. These features can be switched on using optional arguments,

# Swap only AO/MO coeffients and leave occupations and energies unchanged;
# note zero-based indexing
swaps = np.array([[0, 6], [5, 9]])
orb.swap_orbitals(swaps, skip_occs=True, skip_energies=True)

3.1.7. Orbital localization

The current version of PyBEST supports only Pipek-Mezey localization. For more details on Pipek-Mezey localization please refer to Localization of molecular orbitals. Due to reasons of completeness, we will only briefly explain the syntax on how to perform a localization of molecular orbitals. In the Pipek-Mezey localization algorithm, the occupied and virtual block of the (canoncial Hartree-Fock) orbitals are localized using Mulliken projectors. In order to perform a localization, the Mulliken projectors need to be determined. To do so, an instance of the Basis class (here gobasis) and of the LinalgFactory class (here lf) are required.

###############################################################################
## Define Mulliken projectors #################################################
###############################################################################
mulliken = get_mulliken_operators(gobasis)
###############################################################################
## Pipek-Mezey localizaton ####################################################
###############################################################################
loc = PipekMezey(lf, occ_model, mulliken)

As soon as the Mulliken projectors are calculated, they can be passed to the PipekMezey class. The occupied and virtual blocks can be localized individually by two separate function calls,

###############################################################################
## occupied block #############################################################
###############################################################################
loc(hf_output, "occ")
###############################################################################
## virtual block ##############################################################
###############################################################################
loc(hf_output, "virt")

Note

The orbitals have to passed as an IOData container, here hf_results. The corresponding attribute orb_a will be overwritten in the localization procedure.

Frozen core orbitals are also supported during the localization of the occupied orbitals. To define a frozen core, the ncore argument has to be defined during the initialization of the PipekMezey class

###############################################################################
## Pipek-Mezey localizaton with one frozen core orbital #######################
###############################################################################
loc = PipekMezey(lf, occ_model, mulliken, ncore=1)

The actual localization is invoked as in the case of no frozen core orbitals.

3.1.7.1. Example Python scripts

The following, we list a complete example script for a restricted Hartree-Fock calculation of the water molecule followed by orbital localization.

Listing 3.1 data/examples/localization/water_pm.py

from pybest import context
from pybest.gbasis import get_gobasis
from pybest.gbasis.gobasis import (
    compute_eri,
    compute_kinetic,
    compute_nuclear,
    compute_nuclear_repulsion,
    compute_overlap,
)
from pybest.io.iodata import IOData
from pybest.linalg.dense import DenseLinalgFactory
from pybest.localization import PipekMezey
from pybest.part.mulliken import get_mulliken_operators
from pybest.scf import AufbauOccModel
from pybest.wrappers.hf import RHF


###############################################################################
## Set up molecule, define basis set ##########################################
###############################################################################
# Use the XYZ file from PyBEST's test data directory.
fn_xyz = context.get_fn("test/water.xyz")
gobasis = get_gobasis("cc-pvdz", fn_xyz)
###############################################################################
## Define Occupation model, expansion coefficients and overlap ################
###############################################################################
lf = DenseLinalgFactory(gobasis.nbasis)
occ_model = AufbauOccModel(5)
orb_a = lf.create_orbital(gobasis.nbasis)
olp = compute_overlap(gobasis)
###############################################################################
## Construct Hamiltonian ######################################################
###############################################################################
kin = compute_kinetic(gobasis)
na = compute_nuclear(gobasis)
er = compute_eri(gobasis)
external = compute_nuclear_repulsion(gobasis)
###############################################################################
## Do a Hartree-Fock calculation ##############################################
###############################################################################
hf = RHF(lf, occ_model)
hf_output = hf(kin, na, er, external, olp, orb_a)

3.1.8. Generating molden files

PyBEST uses libint2 to dump molden files. As all I/O operations, the generation of molden files is handled by the IOData container class. The corresponding class instance has to contain an instance of the orbitals, for which the molden file has to be generated, as well as an instance of the Basis class. The latter contains all necessary information to dump the molden file.

###############################################################################
## Define Mulliken projectors #################################################
###############################################################################
mulliken = get_mulliken_operators(gobasis)
###############################################################################
## Pipek-Mezey localizaton ####################################################
###############################################################################
loc = PipekMezey(lf, occ_model, mulliken)
###############################################################################
## occupied block #############################################################
###############################################################################
loc(hf_output, "occ")
###############################################################################
## virtual block ##############################################################
###############################################################################
loc(hf_output, "virt")

###############################################################################
## dump Molden file; hf_ already contains the orb_a attribute #################
###############################################################################
hf_output.gobasis = gobasis
hf_output.to_file("water.molden")

Specifying the file ending molden tells the to_file() method to dump a molden file to disk.

In the case of unrestricted orbitals, both alpha and beta orbitals have to be passed to the IOData container. We can also construct an IOData container from scratch that contains all required information,

molden = IOData(gobasis=gobasis,
                orb_a=orb_a,
                orb_b=orb_b
                )
molden.to_file("water.molden")