# 4.1. How to use HORTON as a Hartree-Fock/DFT program¶

This section explains how the relevant Hartree-Fock and DFT features in HORTON are used with Python scripts. This is different from most regular quantum chemistry codes, where you just prepare an input file with a molecular geometry and some control options. Such an input-file interface for HORTON is under construction. The script interface, explained here, is more expressive and laborsome than a regular input file, because you basically have to write the main program. However, it has the advantage of flexibility, since you can tinker more with the individual steps in the calculation.

HORTON can do restricted and unrestricted HF/DFT SCF calculations in various different ways. The following sections cover the typical steps of a calculation:

- Setting up the (molecular electronic) Hamiltonian
- Generating an initial guess of the wavefunction
- Defining the effective Hamiltonian
- Choose an orbital occupation scheme
- Optimize the wavefunction with a self-consistent field algorithm
- Optional (if needed for 7 or 8): Conversion of density and Fock matrix to orbitals
- Optional: Writing SCF results to a file
- Optional: Preparing for a Post-Hartree-Fock calculation

The last section contains an overview of Complete examples that are shipped (and tested) with HORTON. These may be convenient as a starting point when preparing your own scripts. You may also just dive straight into these examples and consult the first five sections to better understand how each example works.

The code snippets below use the following conventions for variable names. The following are defined by setting up the Hamiltonian:

`olp`

: two-index object with the overlap integrals.`kin`

: two-index object with the kinetic energy integrals.`na`

: two-index object with the nuclear attraction integrals.`er`

: four-index object with the electron repulsion integrals.

The following names are also used systematically:

`exp_alpha`

,`exp_beta`

: the expansion coefficient of the alpha and beta orbitals in a (local) basis.

## 4.1.1. Setting up the (molecular electronic) Hamiltonian¶

See Hamiltonians.

You first have to compute/load the two- and four-index objects for the one- and two-body terms in the Hamiltonian. Some DFT implementations in HORTON do not require pre-computed four-center integrals.

## 4.1.2. Generating an initial guess of the wavefunction¶

### 4.1.2.1. Core Hamiltonian guess¶

This guess involves only the one-body terms of the electronic Hamiltonian and completely neglects electron-electron interactions. The orbitals for such a one-body Hamiltonian can be computed without any prior guess.

The function `guess_core_hamiltonian()`

can be
used to compute the core hamiltonian guess. This type of guess only relies on a
simple one-body Hamiltonian, usually consisting of the kinetic energy and the
interaction with the external field. In case of an unrestricted calculation,
this means that the same initial orbital guess is made for the alpha and beta
spins.

The guess for a restricted wavefunction is done as follows:

```
# Create alpha orbitals
exp_alpha = lf.create_expansion()
# Initial guess
guess_core_hamiltonian(olp, kin, na, exp_alpha)
```

For an unrestricted wavefunction, the procedure is very similar:

```
# Create alpha orbitals
exp_alpha = lf.create_expansion()
exp_beta = lf.create_expansion()
# Initial guess
guess_core_hamiltonian(olp, kin, na, exp_alpha, exp_beta)
```

The arguments `exp_alpha`

and `exp_beta`

are treated as output arguments.
Instead of `kin`

and `na`

, you may provide just any set of one-body
operators to construct different types of guesses.

### 4.1.2.2. Randomizing an initial guess¶

The method `rotate_random()`

can be
used to apply a random unitary rotation to all the orbitals (occupied and virtual).
The orbital energies and occupation numbers are not affected. For example, after
constructing a Core Hamiltonian guess, the following line randomizes the
orbitals:

```
# randomly rotate the orbitals (irrespective of occupied or virtual)
exp_alpha.rotate_random()
```

### 4.1.2.3. Modifying the initial guess¶

If needed, you may fine-tune the initial guess by making fine-grained modifications to the orbitals. (These may also be useful for fixing the orbitals that come out of a failed SCF.)

The method

`rotate_2orbitals()`

allows you to mix two orbitals. By default, it rotates the HOMO and LUMO orbitals by 45 degrees:# Mix HOMO and LUMO orbitals exp_alpha.rotate_2orbitals() # Rotate 1st and 6th orbital by 30 deg exp._alpha.rotate_2orbitals(30*deg, 0, 5)

Note that HORTON uses radians as unit for angles, i.e.

`30*deg == np.pi/6`

. Also, HORTON uses zero-based indices, so the arguments`0, 5`

refer to the first and the sixth orbital.The method

`swap_orbitals()`

allows you to swap several orbitals. It takes as an argument an array where each row is a pair of orbitals to swap. For example, the following swaps 1st and 3rd, followed by a swap of 2nd and 4th:# Swap some orbitals swaps = np.array([[0, 2], [1, 3]]) exp_alpha.swap_orbitals(swaps)

### 4.1.2.4. Reading a guess from a file¶

You may also load orbitals from an external file. The file formats `.mkl`

,
`.molden`

, `.fchk`

, or HORTON’s internal `.h5`

are all viable
sources of orbitals. For example, the orbitals from a Gaussian formatted
checkpoint file, `*.fchk`

, may be loaded as follows:

```
# Load fchk file
mol = IOData.from_file('water.fchk')
# Print the number of alpha orbitals (occupied and virtual)
print mol.exp_alpha.nfn
```

Obviously, if you would like to use these orbitals without projecting them onto
a new basis set (as explained in Projecting orbitals from a smaller basis onto a larger one), you are
forced to continue working in the same basis set, which can be accessed
in this example as `mol.obasis`

. See Loading geometry and basis set info from one file
for more details.

### 4.1.2.5. Projecting orbitals from a smaller basis onto a larger one¶

Assuming you have obtained (converged) orbitals in a smaller basis set, you can
try to use these as initial guess after projecting the orbitals onto the
larger basis set. This is exactly what the function
`project_orbitals_mgs()`

does. The following
snippet assumes that the `obasis0`

and `exp_alpha0`

are the small basis set
and a set of orbitals in that basis for the `IOData`

instance `mol`

.

```
# Definition of the bigger basis set
obasis1 = get_gobasis(mol.coordinates, mol.numbers, 'aug-cc-pvtz'):
# Linalg factory for the bigger basis set
lf1 = DenseLinalgFactory(obasis1.nbasis)
# Create a expansion object for the alpha orbitals in the large basis
exp_alpha1 = lf1.create_expansion()
# The actual projection
project_orbitals_msg(obasis0, obasis1, exp_alpha0, exp_alpha1)
```

## 4.1.3. Defining the effective Hamiltonian¶

HORTON implements spin-restricted and spin-unrestricted effective Hamiltonians. Mathematically speaking, these are models for the energy as function of a set of density matrices. The implementation also provides an API to compute the corresponding Fock matrix for every density matrix, i.e. the derivative of the energy toward the density matrix elements.

For the restricted case, the alpha and beta density matrices are assumed to be identical. Hence the energy is only a function of the alpha density matrix. When constructing the Fock matrix, the derivative is divided by two to obtain such that the Fock matrix has the conventional orbital energies as eigenvalues.

\[\begin{split}D^\alpha &\rightarrow E(D^\alpha) \\ &\rightarrow F^\alpha_{\mu\nu} = \frac{1}{2}\frac{\partial E}{\partial D^\alpha_{\nu\mu}}\end{split}\]For the unrestricted case, the alpha and beta density matrices are allowed to be different. Hence, there are also alpha and beta Fock matrices.

\[\begin{split}D^\alpha, D^\beta &\rightarrow E(D^\alpha, D^\beta) \\ &\rightarrow F^\alpha_{\mu\nu} = \frac{\partial E}{\partial D^\alpha_{\nu\mu}} \\ &\rightarrow F^\beta_{\mu\nu} = \frac{\partial E}{\partial D^\beta_{\nu\mu}}\end{split}\]

A generic API for such
density matrix functionals is implemented in the classes
`REffHam`

and
`UEffHam`

. The prefixes `R`

and `U`

are used (also below) to differentiate between restricted and unrestricted
implementations. A Hartree-Fock or DFT effective Hamiltonian is defined by
constructing an instance of the `REffHam`

or `UEffHam`

classes and providing
the necessary energy terms to the constructor.

### 4.1.3.1. Supported energy terms for the effective Hamiltonians¶

All classes below take a `label`

argument to give each term in the Hamiltonian
a name, e.g. used for storing/displaying results. For each class listed below,
follow the hyperlinks to the corresponding documentation for a description of
the constructor arguments.

Simple one-body terms are specified with

`RTwoIndexTerm`

, or`UTwoIndexTerm`

.The direct term of a two-body interaction is specified with

`RDirectTerm`

, or`UDirectTerm`

.The exchange term of a two-body interaction is specified with

`RExchangeTerm`

, or`UExchangeTerm`

.Functionals of the density (or its derivatives) that require numerical integration are all grouped into on term using

`RGridGroup`

, or`UGridGroup`

. This makes it possible to compute at every SCF iteration the density (and its gradients) only once for all terms that depend on the density. This also allows for a similar gain in efficiency when building the Fock matrices. The constructor of a`GridGroup`

class takes a numerical integration grid and a list of instances of the following classes as arguments:- An LDA functional from LibXC can be specified with
`RLibXCLDA`

or`ULibXCLDA`

. - A GGA functional from LibXC can be specified with
`RLibXCGGA`

or`ULibXCGGA`

. - A Hybrid GGA functional from LibXC can be specified with
`RLibXCHybridGGA`

or`ULibXCHybridGGA`

. - An MGGA functional from LibXC can be specified with
`RLibXCMGGA`

or`ULibXCMGGA`

. - A Hybrid MGGA functional from LibXC can be specified with
`RLibXCHybridMGGA`

or`ULibXCHybridMGGA`

. - A numerical implementation of the Hartree term (using an improved version
of Becke’s Poisson solver) can be used instead of the
`RDirectTerm`

or`UDirectTerm`

classes, which require four-center integrals. The relevant classes are`RBeckeHartree`

or`UBeckeHartree`

.

Integration grids are discussed in more detail in the section Building an integration grid. A list of the supported LibXC functionals can be found in LibXC Functionals.

- An LDA functional from LibXC can be specified with

Using these classes, you can construct the Hartree-Fock or a DFT effective Hamiltonian.

### 4.1.3.2. A few typical examples¶

The examples below assume that some or all of the following variables are already defined:

`obasis`

: An orbital basis set.`olp`

: two-index object with the overlap integrals.`kin`

: two-index object with the kinetic energy integrals.`na`

: two-index object with the nuclear attraction integrals.`er`

: four-index object with the electron repulsion integrals.`grid`

: a numerical integration grid.

If you are unfamiliar with any of these, please read the sections
Hamiltonians and How to use HORTON for numerical integration?. The examples below
also make use of the external argument of
`REffHam`

or
`UEffHam`

to add the nuclear-nuclear
repulsion energy to the total energy reported by the effective Hamiltonian.

Restricted Hartree-Fock:

# Construct the restricted HF effective Hamiltonian 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)

Unrestricted Hartree-Fock:

# Construct the restricted HF effective Hamiltonian external = {'nn': compute_nucnuc(mol.coordinates, mol.pseudo_numbers)} terms = [ UTwoIndexTerm(kin, 'kin'), UDirectTerm(er, 'hartree'), UExchangeTerm(er, 'x_hf'), UTwoIndexTerm(na, 'ne'), ] ham = UEffHam(terms, external)

Restricted Kohn-Sham DFT with the Dirac exchange and the VWN correlation functionals:

# Construct the restricted HF effective Hamiltonian external = {'nn': compute_nucnuc(mol.coordinates, mol.pseudo_numbers)} terms = [ RTwoIndexTerm(kin, 'kin'), RDirectTerm(er, 'hartree'), RGridGroup(obasis, grid, [ RLibXCLDA('x'), RLibXCLDA('c_vwn'), ]), RTwoIndexTerm(na, 'ne'), ] ham = REffHam(terms, external)

Unrestricted Kohn-Sham DFT with the PBE GGA exchange and correlation functionals:

# Construct the restricted HF effective Hamiltonian external = {'nn': compute_nucnuc(mol.coordinates, mol.pseudo_numbers)} terms = [ UTwoIndexTerm(kin, 'kin'), UDirectTerm(er, 'hartree'), UGridGroup(obasis, grid, [ ULibXCGGA('x_pbe'), ULibXCGGA('c_pbe'), ]), UTwoIndexTerm(na, 'ne'), ] ham = UEffHam(terms, external)

Restricted Kohn-Sham DFT with the Hybrid GGA functional B3LYP:

# Construct the restricted HF effective Hamiltonian external = {'nn': compute_nucnuc(mol.coordinates, mol.pseudo_numbers)} libxc_term = RLibXCHybridGGA('xc_b3lyp') terms = [ RTwoIndexTerm(kin, 'kin'), RDirectTerm(er, 'hartree'), RGridGroup(obasis, grid, [libxc_term]), RExchangeTerm(er, 'x_hf', libxc_term.get_exx_fraction()), RTwoIndexTerm(na, 'ne'), ] ham = REffHam(terms, external)

Unrestricted Kohn-Sham DFT with the TPSS MGGA exchange and correlation functionals:

# Construct the restricted HF effective Hamiltonian external = {'nn': compute_nucnuc(mol.coordinates, mol.pseudo_numbers)} terms = [ UTwoIndexTerm(kin, 'kin'), UDirectTerm(er, 'hartree'), UGridGroup(obasis, grid, [ ULibXCMGGA('x_tpss'), ULibXCMGGA('c_tpss'), ]), UTwoIndexTerm(na, 'ne'), ] ham = UEffHam(terms, external)

Restricted Kohn-Sham DFT with the Hybrid MGGA functional M05:

# Construct the restricted HF effective Hamiltonian external = {'nn': compute_nucnuc(mol.coordinates, mol.pseudo_numbers)} libxc_term = RLibXCHybridMGGA('xc_m05') terms = [ RTwoIndexTerm(kin, 'kin'), RDirectTerm(er, 'hartree'), RGridGroup(obasis, grid, [libxc_term]), RExchangeTerm(er, 'x_hf', libxc_term.get_exx_fraction()), RTwoIndexTerm(na, 'ne'), ] ham = REffHam(terms, external)

Unrestricted Kohn-Sham DFT with LDA exchange and correlation functionals and with a numerical integration of the Hartree term:

# Construct the restricted HF effective Hamiltonian external = {'nn': compute_nucnuc(mol.coordinates, mol.pseudo_numbers)} terms = [ UTwoIndexTerm(kin, 'kin'), UGridGroup(obasis, grid, [ UBeckeHartree(lmax=8), ULibXCLDA('x'), ULibXCLDA('c_vwn'), ]), UTwoIndexTerm(na, 'ne'), ] ham = UEffHam(terms, external)

## 4.1.4. Choose an orbital occupation scheme¶

Before calling an SCF solver, you have to select a scheme to set the orbital occupations after each SCF iteration, even when the occupation numbers are to remain fixed throughout the calculation. You can use any of the following three options:

`FixedOccModel`

. Keep all occupation numbers fixed at preset values. For example,# Restricted case occ_model = FixedOccModel(np.array([1.0, 1.0, 0.5, 0.5, 0.0])) # Unrestricted case occ_model = FixedOccModel(np.array([1.0, 1.0, 0.5, 0.5, 0.0]), np.array([1.0, 0.7, 1.0, 0.0, 0.3]))

`AufbauOccModel`

. Fill all orbitals according to the Aufbau principle. For example,# Restricted case (three alpha and three beta electrons) occ_model = AufbauOccModel(3.0) # Unrestricted case (two alpha and three beta electrons) occ_model = AufbauOccModel(2.0, 3.0)

`FermiOccModel`

. Use the Fermi-smearing method to fill up the orbitals. [rabuck1999] Only part of the methodology presented Rabuck is implemented. See`FermiOccModel`

for details. For example,# Restricted case (three alpha and three beta electrons, 300K) occ_model = FermiOccModel(3.0, temperature=300) # Unrestricted case (two alpha and three beta electrons, 500K) occ_model = AufbauOccModel(2.0, 3.0, temperature=500)

## 4.1.5. Optimize the wavefunction with a self-consistent field algorithm¶

HORTON supports the following SCF algorithms:

`PlainSCFSolver`

: the ordinary SCF solver. This method just builds and diagonalizes the Fock matrices at every iteration.`ODASCFSolver`

: the optimal damping SCF solver. [cances2001] It uses a cubic interpolation to estimate the optimal mixing between the old and the new density matrices. This is relatively robust but slow.`CDIISSCFSolver`

: the (traditional) commutator direct inversion of the iterative subspace (CDIIS) algorithm, also know as Pulay mixing. [pulay1980] This is usually very efficient but sensitive to the initial guess.`EDIISSCFSolver`

: the energy direct inversion of the iterative subspace (EDIIS) method. [kudin2002] This method works well for the initial iterations but becomes numerically unstable close to the solution. It typically works better with a relatively poor initial guess.`EDIIS2SCFSolver`

: a combination of CDIIS and EDIIS. [kudin2002] This method tries to combine the benefits of both approaches.

The plain SCF solver starts from an initial guess of the orbitals and updates these in-place.

Usage in the restricted case:

# Converge WFN with plain SCF scf_solver = PlainSCFSolver(1e-6) scf_solver(ham, lf, olp, occ_model, exp_alpha)

Usage in the unrestricted case:

# Converge WFN with plain SCF scf_solver = PlainSCFSolver(1e-6) scf_solver(ham, lf, olp, occ_model, exp_alpha, exp_beta)

All other solvers start from an initial guess of the density matrix and update that quantity in-place. The usage pattern is as follow:

Usage in the restricted case:

# Converge WFN with Optimal damping algorithm (ODA) SCF # - Construct the initial density matrix (needed for ODA). occ_model.assign(exp_alpha) dm_alpha = exp_alpha.to_dm() # - SCF solver scf_solver = ODASCFSolver(1e-6) scf_solver(ham, lf, olp, occ_model, dm_alpha)

Usage in the unrestricted case:

# Converge WFN with Optimal damping algorithm (ODA) SCF # - Construct the initial density matrix (needed for ODA). occ_model.assign(exp_alpha, exp_beta) dm_alpha = exp_alpha.to_dm() dm_beta = exp_beta.to_dm() # - SCF solver scf_solver = ODASCFSolver(1e-6) scf_solver(ham, lf, olp, occ_model, dm_alpha, dm_beta)

All SCF solvers support the following two options:

threshold: The convergence threshold for the SCF solver. The exact meaning of this threshold is algorithm-dependent. In case of the `PlainSCFSolver`

, convergence is tested with the RSMD error on the Roothaan equations. In all other cases, the square root of the commutator of the density and the Fock matrix is tested.maxiter: The maximum allowed number of SCF iterations. If convergence is not reached within this number of iterations, a `NoSCFConvergence`

error is raised.

Click on the links of the SCF solver classes above for more optional arguments and their default values.

## 4.1.6. Conversion of density and Fock matrix to orbitals¶

The `PlainSCFSolver`

iteratively updates the
orbitals. Hence, for this SCF solver, there is no need to reconstruct the orbitals
after SCF convergence. For all other SCF algorithms, however, you need to iteratively
update the density matrix and explicitly reconstruct the orbitals after SCF
convergence.

Note

SCF algorithms that update the density matrix may produce a converged
density matrix with fractional occupations for degenerate orbitals at the
Fermi level. The code snippets below properly handle such cases as well. For
more details, refer to
`from_fock_and_dm()`

.

Usage in the restricted case:

# Derive orbitals (coeffs, energies and occupations) from the Fock and density # matrices. The energy is also computed to store it in the output file below. fock_alpha = lf.create_two_index() ham.reset(dm_alpha) ham.compute_energy() ham.compute_fock(fock_alpha) exp_alpha.from_fock_and_dm(fock_alpha, dm_alpha, olp)

Usage in the unrestricted case:

# Derive orbitals (coeffs, energies and occupations) from the Fock and density # matrices. The energy is also computed to store it in the output file below. fock_alpha = lf.create_two_index() fock_beta = lf.create_two_index() ham.reset(dm_alpha, dm_beta) ham.compute_energy() ham.compute_fock(fock_alpha, fock_beta) exp_alpha.from_fock_and_dm(fock_alpha, dm_alpha, olp) exp_beta.from_fock_and_dm(fock_beta, dm_beta, olp)

## 4.1.7. Writing SCF results to a file¶

Once you have obtained a converged set of molecular orbitals, they can be stored in a file. Two file formats can be used: (i) HORTON’s internal format based on HDF5 and (ii) the Molden format. The internal format has the advantage that all data is stored in full precision and that you can include as much additional information as you like. The Molden format is useful for visualization purposes or when you would like to use the orbitals in another quantum chemistry code. Both formats can be used to store the wavefunction and to use it as an initial guess afterwards. For more background information on data formats in HORTON, refer to Data file formats (input and output).

The two code snippets below show how a wavefunction can be written to a file.
In both examples, the object `mol`

is an instance of `IOData`

that was
created earlier in the script by loading an `.xyz`

file as follows:

```
mol = IOData.from_file('some_molecule.xyz')
```

When the `.xyz`

file is loaded, the attributes `numbers`

and `coordinates`

of the object `mol`

are already set as required for the `.molden`

format.

Usage in the restricted case:

# Assign results to the molecule object and write it to a file, e.g. for # later analysis. Note that the ODA algorithm can only really construct an # optimized density matrix and no orbitals. mol.title = 'RKS computation on water' mol.energy = ham.cache['energy'] mol.obasis = obasis mol.exp_alpha = exp_alpha mol.dm_alpha = dm_alpha # useful for visualization: mol.to_file('water.molden') # useful for post-processing (results stored in double precision): mol.to_file('water.h5')

Usage in the unrestricted case:

# Assign results to the molecule object and write it to a file, e.g. for # later analysis. Note that the ODA algorithm can only really construct an # optimized density matrix and no orbitals. mol.title = 'UKS computation on methyl' mol.energy = ham.cache['energy'] mol.obasis = obasis mol.exp_alpha = exp_alpha mol.exp_beta = exp_beta mol.dm_alpha = dm_alpha mol.dm_beta = dm_beta # useful for visualization: mol.to_file('methyl.molden') # useful for post-processing (results stored in double precision): mol.to_file('methyl.h5')

## 4.1.8. Preparing for a Post-Hartree-Fock calculation¶

Once the SCF has converged and you have obtained a set of orbitals, you can use
these orbitals to convert the integrals from the atomic-orbital (AO) basis to
the molecular-orbital (MO) basis. Two implementations are available in HORTON:
(i) using all molecular orbitals or (ii) by specifying a frozen core and active
set of orbitals. A full example, `rhf_n2_dense.py`

, which covers both options
and which includes dumping the transformed integrals to a file, is given in the
section Complete examples below.

The conversion to an MO basis is useful for post-HF calculations. For such purposes, it is also of interest to sum all one-body operators into a single term:

```
lf = DenseLinalgFactory(obasis.nbasis)
kin = obasis.compute_kinetic(lf)
na = obasis.compute_nuclear_attraction(mol.coordinates, mol.pseudo_numbers, lf)
one = kin.copy()
one.iadd(na)
```

### 4.1.8.1. Transforming the Hamiltonian to the molecular-orbital (MO) basis¶

The function `transform_integrals()`

can be used for
this purpose. It works differently for restricted and unrestricted orbitals:

In the case of restricted (Hartree-Fock) orbitals, there is just one transformed one-body and one transformed two-body operator:

(one_mo,), (two_mo,) = transform_integrals(one, er, 'tensordot', exp_alpha)

In the case of unrestricted (Hartree-Fock) orbitals, there are two transformed one-body and three transformed two-body operators:

one_mo_ops, two_mo_ops = transform_integrals(one, er, 'tensordot', exp_alpha, exp_beta) one_mo_alpha, one_mo_beta = one_mo_ops two_mo_alpha_alpha, two_mo_alpha_beta, two_mo_beta_beta = two_mo_ops

### 4.1.8.2. Reducing the Hamiltonian to an active space¶

The active-space feature only works for restricted orbitals. The new one-electron integrals \(\tilde{t}_{pq}\) become

where \(t_{pq}\) is the element \(pq\) of the old one-electron integrals and \(\langle pi \vert qi \rangle\) is the appropriate two-electron integral in physicist’s notation. The core energy of the active space is calculated as

where the two-electron integrals \(\langle pq \vert rs \rangle\) contain only
the elements with active orbital indices \(p,q,r,s\). This type of conversion
is implemented in the function `split_core_active()`

.
It is used as follows:

```
one_small, two_small, core_energy = split_core_active(one, er,
external['nn'], exp_alpha, ncore, nactive)
```

## 4.1.9. Complete examples¶

The following is a basic example of a restricted Hartree-Fock calculation of water. It contains all the steps discussed in the previous sections.

```
import numpy as np
from horton import * # pylint: disable=wildcard-import,unused-wildcard-import
# Hartree-Fock calculation
# ------------------------
# Construct a molecule from scratch
bond_length = 1.098*angstrom
mol = IOData(title='dinitrogen')
mol.coordinates = np.array([[0.0, 0.0, 0.0], [0.0, 0.0, bond_length]])
mol.numbers = np.array([7, 7])
# Create a Gaussian basis set
obasis = get_gobasis(mol.coordinates, mol.numbers, 'cc-pvdz')
# Create a linalg factory
lf = DenseLinalgFactory(obasis.nbasis)
# Compute Gaussian integrals
olp = obasis.compute_overlap(lf)
kin = obasis.compute_kinetic(lf)
na = obasis.compute_nuclear_attraction(mol.coordinates, mol.pseudo_numbers, lf)
er = obasis.compute_electron_repulsion(lf)
# Create alpha orbitals
exp_alpha = lf.create_expansion()
# Initial guess
guess_core_hamiltonian(olp, kin, na, exp_alpha)
# Construct the restricted HF effective Hamiltonian
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)
# Decide how to occupy the orbitals (7 alpha electrons)
occ_model = AufbauOccModel(7)
# Converge WFN with plain SCF
scf_solver = PlainSCFSolver(1e-6)
scf_solver(ham, lf, olp, occ_model, exp_alpha)
# Write SCF results to a file
# ---------------------------
# Assign results to the molecule object and write it to a file, e.g. for
# later analysis
mol.title = 'RHF computation on dinitrogen'
mol.energy = ham.cache['energy']
mol.obasis = obasis
mol.exp_alpha = exp_alpha
# useful for visualization:
mol.to_file('n2-scf.molden')
# useful for post-processing (results stored in double precision)
mol.to_file('n2-scf.h5')
# Export Hamiltonian in Hartree-Fock molecular orbital basis (all orbitals active)
# --------------------------------------------------------------------------------
# Transform orbitals
one = kin.copy()
one.iadd(na)
two = er
(one_mo,), (two_mo,) = transform_integrals(one, two, 'tensordot', mol.exp_alpha)
# Write files
mol_all_active = IOData(core_energy=external['nn'], one_mo=one_mo, two_mo=two_mo)
# useful for exchange with other codes
mol_all_active.to_file('n2.FCIDUMP')
# useful for exchange with other HORTON scripts
mol_all_active.to_file('n2-hamiltonian.h5')
# Export Hamiltonian in Hartree-Fock molecular orbital basis for CAS(8,8)
# -----------------------------------------------------------------------
# Transform orbitals
one_small, two_small, core_energy = split_core_active(one, er,
external['nn'], exp_alpha, ncore=2, nactive=8)
# Write files
mol_cas88 = IOData(core_energy=core_energy, one_mo=one_mo, two_mo=two_mo, nelec=8, ms2=0, lf=lf)
# useful for exchange with other codes
mol_cas88.to_file('n2-cas8-8.FCIDUMP')
# useful for exchange with other HORTON scripts
mol_cas88.to_file('n2-hamiltonian-cas8-8.h5')
```

The directory `data/examples/hf_dft`

contains many more examples that use the
different options discussed above. The following table shows which features are
used in which example.

File name | LOT | SCF | Description |
---|---|---|---|

rhf_h2_cholesky.py | RHF/6-31G | PlainSCFSolver | Basic RHF example with Cholesky matrices, includes export of Hamiltonian |

rhf_n2_dense.py | RHF/cc-pvtz | PlainSCFSolver | Basic RHF example with dense matrices, includes export of Hamiltonian |

rhf_water_dense.py | RHF/3-21G | PlainSCFSolver | Basic RHF example with dense matrices |

uhf_methyl_dense.py | UHF/3-21G | PlainSCFSolver | Basic UHF example with dense matrices |

rhf_water_cholesky.py | RHF/cc-pVDZ | PlainSCFSolver | Basic RHF example with Cholesky decomposition of the ERI |

uhf_methyl_cholesky.py | UHF/cc-pVDZ | PlainSCFSolver | Basic UHF example with Cholesky decomposition of the ERI |

rks_water_lda.py | RKS/6-31G(d) | ODASCFSolver | Basic RKS DFT example with LDA exhange-correlation functional (Dirac+VWN) |

uks_methyl_lda.py | UKS/6-31G(d) | ODASCFSolver | Basic UKS DFT example with LDA exhange-correlation functional (Dirac+VWN) |

rks_water_gga.py | RKS/6-31G(d) | CDIISSCFSolver | Basic RKS DFT example with GGA exhange-correlation functional (PBE) |

uks_methyl_gga.py | UKS/6-31G(d) | CDIISSCFSolver | Basic UKS DFT example with GGA exhange-correlation functional (PBE) |

rks_water_hybgga.py | RKS/6-31G(d) | EDIIS2SCFSolver | Basic RKS DFT example with hyrbid GGA exhange-correlation functional (B3LYP) |

uks_methyl_hybgga.py | UKS/6-31G(d) | EDIIS2SCFSolver | Basic UKS DFT example with hyrbid GGA exhange-correlation functional (B3LYP) |

rks_water_mgga.py | RKS/6-31G(d) | CDIISSCFSolver | Basic RKS DFT example with MGGA exhange-correlation functional (TPSS) |

uks_methyl_mgga.py | UKS/6-31G(d) | CDIISSCFSolver | Basic UKS DFT example with MGGA exhange-correlation functional (TPSS) |

rks_water_hybmgga.py | RKS/6-31G(d) | CDIISSCFSolver | Basic RKS DFT example with hybrid MGGA exhange-correlation functional (TPSS) |

uks_methyl_hybmgga.py | UKS/6-31G(d) | CDIISSCFSolver | Basic UKS DFT example with hybrid MGGA exhange-correlation functional (M05) |

rks_water_numlda.py | RKS/6-31G(d) | CDIISSCFSolver | RKS DFT example with LDA and numerical Hartree |

uks_methyl_numlda.py | UKS/6-31G(d) | CDIISSCFSolver | UKS DFT example with LDA and numerical Hartree |

rks_water_numgga.py | RKS/6-31G(d) | CDIISSCFSolver | RKS DFT example with GGA and numerical Hartree |

uks_methyl_numgga.py | UKS/6-31G(d) | CDIISSCFSolver | UKS DFT example with GGA and numerical Hartree |

A more elaborate example can be found in `data/examples/hf_compare`

. It
contains a script that systematically computes all elements in the periodic table
(for different charges and multiplicities), and compares the results with
outputs obtained with Gaussian. See the `README`

file for instructions how to
run this example.