4.3. The perturbation theory module

4.3.1. Supported features

The perturbation theory module supports spin-restricted orbitals and the DenseLinalgFactory. HORTON offers the following flavors of perturbation theory:

1. Moller-Plesset Perturbation theory of second order with a restricted, closed-shell Hartree-Fock reference function (see Moller-Plesset perturbation theory of second order)
2. PTa with an AP1roG reference function (see The PTa module)
3. PTb with an AP1roG reference function (see The PTb module)

4.3.2. Moller-Plesset perturbation theory of second order

4.3.2.1. Getting started

The RMP2 module requires a restricted Hartree-Fock reference wavefunction (see How to use HORTON as a Hartree-Fock/DFT program), its energy, and the one- and two-electron integrals as input arguments.

The Hartree-Fock energy can be calculated after SCF convergence as follows

ehf = ham.compute_energy()

where ham is an instance of the restricted effective Hamiltonian class REffHam (see FIXME), which contains all one- and two-electron terms. Note that all one-electron terms have to be combined into one single operator term. This can be done in the following way if the Hamiltonian contains the kinetic energy of the electrons and the electron-nuclear attraction,

###############################################################################
## Calculate kinetic energy (kin) and nuclear attraction (na) term ############
###############################################################################
kin = obasis.compute_kinetic(lf)
na = obasis.compute_nuclear_attraction(mol.coordinates, mol.pseudo_numbers, lf)
###############################################################################
## Combine one-electron integrals to single Hamiltonian #######################
###############################################################################
one = kin.copy()
one.iadd(na)

4.3.2.2. How to set-up an MP2 calculation

First, create an instance of the RMP2 class

mp2 = RMP2(lf, occ_model)

with arguments

lf:A LinalgFactory instance (see FIXME). Only DenseLinalgFactory is supported occ_model:(AufbauOccModel instance) an Aufbau occupation model

A function call initiates an MP2 calculation,

emp2, tmp2 = mp2(one, two, orb, **{'eref': ehf, 'FIXME': 'tensordot'})

with arguments

one:(TwoIndex instance) the one-electron integrals two:(FourIndex instance) the two-electron integrals orb:(Expansion instance) the AO/MO coefficient matrix

and keyword arguments

eref:(float) the Hartree-Fock reference energy (default float('nan')) (see Getting started) indextrans:(str, optional) the 4-index transformation. Choice between tensordot (default) and einsum. tensordot is faster than einsum, requires, however, more memory. If DenseLinalgFactory is used, the memory requirement scales as \(2N^4\) for einsum and \(3N^4\) for tensordot, respectively. Due to the storage of the two-electron integrals, the total amount of memory increases to \(3N^4\) for einsum and \(4N^4\) for tensordot, respectively.

The function call gives 2 return values:

emp2:(list of float) the MP2 energy correction (first element) tmp2:(list of FourIndex instances) the MP2 amplitudes (first element). The double excitation amplitudes \(t_{ij}^{ab}\) are stored as t[i,a,j,b]

4.3.3. The PTa module

4.3.3.1. Getting started

PTa [Limacher2014] adds dynamic electron correlation effects on top of an AP1roG wavefunction (see The AP1roG module). You have to optimize an AP1roG wavefunction (see AP1roG with orbital optimization or AP1roG without orbital optimization), before the PTa energy correction can be determined.

4.3.3.2. How to set-up a calculation

First, create an instance of the PTa class

pta = PTa(lf, occ_model)

with arguments

lf:A LinalgFactory instance (see FIXME). Only DenseLinalgFactory is supported occ_model:(AufbauOccModel instance) an Aufbau occupation model

A function call initiates an PTa calculation,

epta, tpta = pta(one, two, orb, c, **{'eref': energy, 'ecore': ecore, 'indextrans': 'tensordot'})

with arguments

one:(TwoIndex instance) the one-electron integrals (the same integrals as used in the AP1roG module The AP1roG module) two:(FourIndex instance) the two-electron integrals (the same integrals as used in the AP1roG module The AP1roG module) orb:(Expansion instance) the optimized AP1roG MO coefficient matrix c:(TwoIndex instance) the geminal coefficients (see AP1roG with orbital optimization)

and keyword arguments

eref:(float) the AP1roG reference energy (default float('nan')) (see AP1roG with orbital optimization how to get the AP1roG reference energy) ecore:(float) the core energy (default float('nan')). Usually, the nuclear repulsion term indextrans:(str, optional) the 4-index transformation. Choice between tensordot (default) and einsum. tensordot is faster than einsum, requires, however, more memory. If DenseLinalgFactory is used, the memory requirement scales as \(2N^4\) for einsum and \(3N^4\) for tensordot, respectively. Due to the storage of the two-electron integrals, the total amount of memory increases to \(3N^4\) for einsum and \(4N^4\) for tensordot, respectively.

The function call gives 2 return values:

epta:(list of float) the PTa energy corrections. Contains the total PTa energy correction (first element), its seniority-zero contribution (second element), its seniority-two contribution (third element), and its seniority-four contribution (fourth element) tpta:(list of FourIndex instances) the PTa amplitudes (first element). The double excitation amplitudes \(t_{ij}^{ab}\) are stored as t[i,a,j,b]

4.3.4. The PTb module

4.3.4.1. Getting started

PTb [Limacher2014] represents a different flavor to add dynamic electron correlation effects on top of an AP1roG wavefunction (see The AP1roG module). Similarly to PTa (The PTa module), you have to optimize an AP1roG wavefunction (see AP1roG with orbital optimization or AP1roG without orbital optimization), before the PTb energy correction can be determined.

4.3.4.2. How to set-up a calculation

First, create an instance of the PTb class

ptb = PTb(lf, occ_model)

with arguments

lf:A LinalgFactory instance (see FIXME). Only DenseLinalgFactory is supported occ_model:(AufbauOccModel instance) an Aufbau occupation model

A function call initiates an PTb calculation,

eptb, tptb = ptb(one, two, orb, c, **{'eref': energy, 'ecore': ecore})

with arguments

one:(TwoIndex instance) the one-electron integrals (the same integrals as used in the AP1roG module The AP1roG model) two:(FourIndex instance) the two-electron integrals (the same integrals as used in the AP1roG module The AP1roG model) orb:(Expansion instance) the optimized AP1roG MO coefficient matrix c:(TwoIndex instance) the geminal coefficients (see AP1roG with orbital optimization)

Note that optional keyword arguments have been omitted. All keyword arguments are summarized in Summary of keyword arguments.

The function call gives 2 return values:

eptb:(list of float) the PTb energy corrections. Contains the total PTb energy correction (first element), its seniority-zero contribution (second element), its seniority-two contribution (third element), and its seniority-four contribution (fourth element) tptb:(list of FourIndex instances) the PTb amplitudes (first element). The double excitation amplitudes \(t_{ij}^{ab}\) are stored as t[i,a,j,b]

4.3.4.3. Summary of keyword arguments

indextrans:(str) 4-index Transformation. Choice between tensordot (default) and einsum. tensordot is faster than einsum, requires, however, more memory. If DenseLinalgFactory is used, the memory requirement scales as \(2N^4\) for einsum and \(3N^4\) for tensordot, respectively. Due to the storage of the two-electron integrals, the total amount of memory increases to \(3N^4\) for einsum and \(4N^4\) for tensordot, respectively. eref:(float) AP1roG reference energy (default float('nan')) ecore:(float) core energy (default float('nan')) threshold:(float) optimization threshold for amplitudes (default 1e-6) maxiter:(int) maximum number of iterations (default 200) guess:(1-dim np.array) initial guess (default None). In not provided, an initial guess containing random numbers in the interval \((0,0.01]\) is generated. The random guess preserves the symmetry of the PTb amplitudes, that is, \(t_{ij}^{ab}=t_{ji}^{ba}\). If a user-defined guess is provided, the elements of \(t_{ij}^{ab}\) have to be indexed in C-like order

4.3.5. Example Python scripts

Several complete examples can be found in the directory data/examples/perturbation_theory. Three of these are discussed in the following subsections.

4.3.5.1. MP2 calculation on the water molecule

This is a basic example on how to perform a RMP2 calculation in HORTON. This script performs a RMP2 calculation on the water molecule using the cc-pVDZ basis set.

data/examples/perturbation_theory/mp2_water_cc-pvdz.py

from horton import *
from horton.test.common import numpy_seed

###############################################################################
## Set up molecule, define basis set ##########################################
###############################################################################
# get the XYZ file from HORTON's test data directory
fn_xyz = context.get_fn('test/water.xyz')
mol = IOData.from_file(fn_xyz)
obasis = get_gobasis(mol.coordinates, mol.numbers, 'cc-pvdz')
###############################################################################
## Define Occupation model, expansion coefficients and overlap ################
###############################################################################
lf = DenseLinalgFactory(obasis.nbasis)
occ_model = AufbauOccModel(5)
orb = lf.create_expansion(obasis.nbasis)
olp = obasis.compute_overlap(lf)
###############################################################################
## Construct Hamiltonian ######################################################
###############################################################################
kin = obasis.compute_kinetic(lf)
na = obasis.compute_nuclear_attraction(mol.coordinates, mol.pseudo_numbers, lf)
er = obasis.compute_electron_repulsion(lf)
external = {'nn': compute_nucnuc(mol.coordinates, mol.pseudo_numbers)}
terms = [
    RTwoIndexTerm(kin, 'kin'),
    RDirectTerm(er, 'hartree'),
    RExchangeTerm(er, 'x_hf'),
    RTwoIndexTerm(na, 'ne'),
]
ham = REffHam(terms, external)
###############################################################################
## Perform initial guess ######################################################
###############################################################################
guess_core_hamiltonian(olp, kin, na, orb)
###############################################################################
## Do a Hartree-Fock calculation ##############################################
###############################################################################
scf_solver = PlainSCFSolver(1e-6)
scf_solver(ham, lf, olp, occ_model, orb)
###############################################################################
## Get Hartree-Fock energy ####################################################
###############################################################################
ehf = ham.compute_energy()
###############################################################################
## Combine one-electron integrals to single Hamiltonian #######################
###############################################################################
one = kin.copy()
one.iadd(na)

###############################################################################
## Do RMP2 calculation ########################################################
###############################################################################
mp2 = RMP2(lf, occ_model)
with numpy_seed():  # reproducible 'random' numbers to make sure it always works
    emp2, tmp2 = mp2(one, er, orb, **{'eref': ehf})

4.3.5.2. PTa calculation on the water molecule

This is a basic example on how to perform a PTa calculation in HORTON. This script performs a PTa calculation on the water molecule using the cc-pVDZ basis set.

data/examples/perturbation_theory/pta_water_cc-pvdz.py

from horton import *
from horton.test.common import numpy_seed

###############################################################################
## Set up molecule, define basis set ##########################################
###############################################################################
# get the XYZ file from HORTON's test data directory
fn_xyz = context.get_fn('test/water.xyz')
mol = IOData.from_file(fn_xyz)
obasis = get_gobasis(mol.coordinates, mol.numbers, 'cc-pvdz')
###############################################################################
## Define Occupation model, expansion coefficients and overlap ################
###############################################################################
lf = DenseLinalgFactory(obasis.nbasis)
occ_model = AufbauOccModel(5)
orb = lf.create_expansion(obasis.nbasis)
olp = obasis.compute_overlap(lf)
###############################################################################
## Construct Hamiltonian ######################################################
###############################################################################
kin = obasis.compute_kinetic(lf)
na = obasis.compute_nuclear_attraction(mol.coordinates, mol.pseudo_numbers, lf)
er = obasis.compute_electron_repulsion(lf)
external = {'nn': compute_nucnuc(mol.coordinates, mol.pseudo_numbers)}
terms = [
    RTwoIndexTerm(kin, 'kin'),
    RDirectTerm(er, 'hartree'),
    RExchangeTerm(er, 'x_hf'),
    RTwoIndexTerm(na, 'ne'),
]
ham = REffHam(terms, external)
###############################################################################
## Perform initial guess ######################################################
###############################################################################
guess_core_hamiltonian(olp, kin, na, orb)
###############################################################################
## Do a Hartree-Fock calculation ##############################################
###############################################################################
scf_solver = PlainSCFSolver(1e-6)
scf_solver(ham, lf, olp, occ_model, orb)
###############################################################################
## Combine one-electron integrals to single Hamiltonian #######################
###############################################################################
one = kin.copy()
one.iadd(na)

###############################################################################
## Do OO-AP1roG optimization ##################################################
###############################################################################
ap1rog = RAp1rog(lf, occ_model)
with numpy_seed():  # reproducible 'random' numbers to make sure it always works
    energy, g, l = ap1rog(one, er, external['nn'], orb, olp, True)

###############################################################################
## Do PTa calculation #########################################################
###############################################################################
pta = PTa(lf, occ_model)
with numpy_seed():  # reproducible 'random' numbers to make sure it always works
    energypta, amplitudes = pta(one, er, orb, g, **{'eref': energy, 'ecore': external['nn']})

4.3.5.3. PTb calculation on the water molecule

This is a basic example on how to perform a PTb calculation in HORTON. This script performs a PTb calculation on the water molecule using the cc-pVDZ basis set.

data/examples/perturbation_theory/ptb_water_cc-pvdz.py

from horton import *
from horton.test.common import numpy_seed

###############################################################################
## Set up molecule, define basis set ##########################################
###############################################################################
# get the XYZ file from HORTON's test data directory
fn_xyz = context.get_fn('test/water.xyz')
mol = IOData.from_file(fn_xyz)
obasis = get_gobasis(mol.coordinates, mol.numbers, 'cc-pvdz')
###############################################################################
## Define Occupation model, expansion coefficients and overlap ################
###############################################################################
lf = DenseLinalgFactory(obasis.nbasis)
occ_model = AufbauOccModel(5)
orb = lf.create_expansion(obasis.nbasis)
olp = obasis.compute_overlap(lf)
###############################################################################
## Construct Hamiltonian ######################################################
###############################################################################
kin = obasis.compute_kinetic(lf)
na = obasis.compute_nuclear_attraction(mol.coordinates, mol.pseudo_numbers, lf)
er = obasis.compute_electron_repulsion(lf)
external = {'nn': compute_nucnuc(mol.coordinates, mol.pseudo_numbers)}
terms = [
    RTwoIndexTerm(kin, 'kin'),
    RDirectTerm(er, 'hartree'),
    RExchangeTerm(er, 'x_hf'),
    RTwoIndexTerm(na, 'ne'),
]
ham = REffHam(terms, external)
###############################################################################
## Perform initial guess ######################################################
###############################################################################
guess_core_hamiltonian(olp, kin, na, orb)
###############################################################################
## Do a Hartree-Fock calculation ##############################################
###############################################################################
scf_solver = PlainSCFSolver(1e-6)
scf_solver(ham, lf, olp, occ_model, orb)
###############################################################################
## Combine one-electron integrals to single Hamiltonian #######################
###############################################################################
one = kin.copy()
one.iadd(na)

###############################################################################
## Do OO-AP1roG optimization ##################################################
###############################################################################
ap1rog = RAp1rog(lf, occ_model)
with numpy_seed():  # reproducible 'random' numbers to make sure it always works
    energy, g, l = ap1rog(one, er, external['nn'], orb, olp, True)

###############################################################################
## Do PTb calculation #########################################################
###############################################################################
ptb = PTb(lf, occ_model)
with numpy_seed():  # reproducible 'random' numbers to make sure it always works
    energyptb, amplitudes = ptb(one, er, orb, g, **{'eref': energy, 'ecore': external['nn'], 'threshold': 1e-6})