14.1. Localization of molecular orbitals¶

In general, the localization algorithm optimizes some localization function by an orthogonal transformation of the orbitals. Given orbitals \(\vert i \rangle\), the localized orbital \(\vert \tilde{i} \rangle\) can be obtained by some transformation

\[\vert \tilde{i} \rangle = \sum_k \vert k \rangle \exp(-\mathbf{\kappa})_{ki},\]

where

\[\mathbf{\kappa} = \sum_{k > l} \kappa_{kl} (a^\dagger_k a_l - a^\dagger_l a_k)\]

and \(\kappa_{kl}\) is determined by the optimization of the localization function.

Many of the localization schemes, and thus the result of the localization, differ by the localization function. The localization function somehow measures the localization of the orbitals. So far, PyBEST only supports the Pipek-Mezey localization. [pipek1989]

14.1.1. Pipek-Mezey localization¶

In the Pipek-Mezey scheme, the Pipek-Mezey localization function, \(D\), is maximized.

\[D = \sum_{i} \sum_{A \in \textrm{atoms}} (Q_{ii}^A)^2,\]

where \(\mathbf{Q}^A\) is the atomic population matrix. The atomic population matrix can be obtained from the overlap of the atomic basis, the molecular orbitals from the atomic basis, the occupation of each molecular orbital, and some weighted projection of atomic basis function within each atom.

For example, if the Mulliken population analysis is used, the projectors are obtained through get_mulliken_operators(). The Pipek-Mezey localization is initialized through a series of PipekMezey function calls, __call__().

Note

The virtual and occupied blocks are localized separately. Thus, each block requires a separate function call.

The code snippet below summarizes all required steps (see also Orbital localization, where all arguments are explained)

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

14.1.2. Example Python scripts¶

14.1.2.1. Pipek-Mezey localization of restricted Hartree-Fock orbitals for the water molecule¶

This is a basic example on how to perform a Pipek-Mezey localization in PyBEST. This script performs a Pipek-Mezey localization for the water molecule using the cc-pVDZ basis set and Mulliken projectors. The localized orbitals are then dump to disk in the molden file format (see also Generating molden files).

Listing 14.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_ = hf(kin, na, er, external, olp, orb_a)

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

###############################################################################
## dump Molden file ###########################################################
###############################################################################
molden = IOData(gobasis=gobasis, orb_a=orb_a)
molden.to_file("water.molden")