3.1. Molecular Hamiltonians¶

A molecular Hamiltonian is typically set up in three steps. First, you load or construct a molecular geometry. Then, you generate a Gaussian basis set for this molecule. Finally, you can calculate all kinds of matrix elements with that basis to define a Hamiltonian.

3.1.1. Specifying the molecular geometry¶

HORTON can load/dump the molecular geometry from/to different file formats. The file format is determined automatically using the extension or prefix of the filename. For more information on supported file formats, refer to Data file formats (input and output).

3.1.1.1. Reading the molecular geometry from file¶

The molecular geometry can be read from file using the method from_file() of the IOData class,

mol = IOData.from_file(filename)

3.1.1.2. Constructing a molecular geometry from scratch¶

You can also generate molecular geometries with Python code. The following example constructs a ring of Hydrogen atoms and writes it to an XYZ file

data/examples/hamiltonian/hydrogen_ring.py¶

import numpy as np

from horton import *  # pylint: disable=wildcard-import,unused-wildcard-import


# Define the ring.
natom = 11
spacing = 1.3  # distance between two neighboring atoms in bohr
radius = spacing/(2*np.sin(np.pi/natom))

# Define the coordinates and elements.
coordinates = np.zeros((natom, 3))
numbers = np.ones(natom, dtype=int)  # must be integers
for iatom in xrange(natom):
    angle = (2*np.pi/natom)*iatom
    coordinates[iatom, 0] = radius*np.cos(angle)
    coordinates[iatom, 1] = radius*np.sin(angle)

# Write the molecule to an XYZ file (optional).
mol = IOData(coordinates=coordinates, numbers=numbers, title='H Ring')
mol.to_file('ring.xyz')

3.1.2. Specifying the basis set¶

HORTON supports basis sets consisting of generally contracted Cartesian Gaussian functions up to a maximum angular momentum of seven (I functions). HORTON is using the same basis set format as NWChem, and the basis sets can be downloaded from the EMSL webpage (https://bse.pnl.gov/bse/portal).

HORTON is distributed with most of the popular basis sets. A list of currently supported built-in basis sets can be found here: Standard basis sets. The basis set for a given molecule is constructed with the function get_gobasis()

3.1.2.1. Unique basis set for all atoms¶

Usually, you want to define a unique (the same) basis set for the whole system. This can be done by a function call.

obasis = get_gobasis(mol.coordinates, mol.numbers, 'cc-pvdz')

where mol.coordinates and mol.numbers are read from file (see Reading the molecular geometry from file), and cc-pvdz is the cc-pVDZ basis set.

3.1.2.2. Specifying different basis sets for different atoms¶

In some cases, you may want to specify different basis sets for different atoms. For example, you might like to use the 3-21G basis set for the hydrogen atom, the 6-31G basis set for the carbon atom, and STO-3G for all remaining atoms:

obasis = get_gobasis(mol.coordinates, mol.numbers, 'sto-3g',
                     element_map={'H':'3-21g', 'C':'6-31g'})

where mol.coordinates and mol.numbers are read from file (see Reading the molecular geometry from file), and sto-3g, 3-21g and 6-31g are the basis set names (see Standard basis sets)

Alternatively, the same result can be obtained by substituting the H and C symbols with their atomic numbers:

obasis = get_gobasis(mol.coordinates, mol.numbers, 'sto-3g',
                     element_map={1:'3-21g', 6:'6-31g'})

You can also override the default basis for selected atoms based on their index, i.e. position in the list of atoms that specify the molecule:

obasis = get_gobasis(mol.coordinates, mol.numbers, 'sto-3g',
                     index_map={0:'3-21g', 2:'6-31g'})

The above example uses the 3-21g basis for the first atom, the 6-31g basis for the third atom and the sto-3g basis for all other atoms.

3.1.2.3. Loading custom basis sets from file¶

You can also use other basis sets besides the ones that are shipped with HORTON. It is assumed that the basis is available in NWChem format:

mybasis = GOBasisFamily('myname', filename='mybasis.nwchem'),
obasis = get_gobasis(mol.coordinates, mol.numbers, mybasis)

Anywhere you can specify a built-in basis set with a string, you can also use instance of the GOBasisFamily class (mybasis in the example above), e.g. in the arguments default, element_map and index_map of get_gobasis.

3.1.2.4. Defining basis sets with Python code¶

In some circumstances, it may be useful to generate the basis set with some Python code. For example, the following code generates an even tempered basis for Lithium (without polarization functions):

data/examples/hamiltonian/even_tempered_li.py¶

import numpy as np

from horton import *  # pylint: disable=wildcard-import,unused-wildcard-import


# specify the even tempered basis set
alpha_low = 5e-3
alpha_high = 5e2
nbasis = 30
lnratio = (np.log(alpha_high) - np.log(alpha_low))/(nbasis-1)

# build a list of "contractions". These aren't real contractions as every
# contraction only contains one basis function.
bcs = []
for ibasis in xrange(nbasis):
    alpha = alpha_low * np.exp(lnratio * ibasis)
    # arguments of GOBasisContraction:
    #     shell_type, list of exponents, list of contraction coefficients
    bcs.append(GOBasisContraction(0, np.array([alpha]), np.array([1.0])))

# Finish setting up the basis set:
ba = GOBasisAtom(bcs)
obasis = get_gobasis(np.array([[0.0, 0.0, 0.0]]), np.array([3]), default=ba)

All basis functions in this example are just single s-type primitives, i.e. no contactions are used. At the end of the example, the basis set is constructed for a single Li atom in the origin.

Note that get_gobasis also accepts instances of GOBasisAtom for the arguments default, element_map and index_map.

3.1.2.5. Loading geometry and basis set info from one file¶

When post-processing results from other programs, it may be desirable to use exactly the same basis set that was used in the other program (Gaussian, ORCA, PSI4, etc). This can be achieved by loading the geometry, basis set and wavefunction from one of the following formats: .mkl, .molden, .fchk, or HORTON’s internal .h5 format. In principle, it is also possible to load the basis set from a .wfn file, but keep in mind that this format does not support contracted basis functions, so HORTON will then use a decontracted basis set, which is less efficient.

You simply use the IOData.from_file method to load the file. The orbital basis object is then available in the obasis attribute if the return value. For example:

# Load the geometry, orbital basis and wavefunction from a Gaussian
# formatted checkpoint file:
mol = IOData.from_file('water.fchk')

# Print the number of basis functions
print mol.obasis.nbasis

3.1.3. Computing (Hamiltonian) matrix elements¶

Given a GOBasis instance (the obasis object from the examples in the previous section), you can generate the two-index and four-index objects for the molecular electronic Hamiltonian. It may also be useful to construct the overlap operator as the Gaussian basis sets are not orthonormal.

Before computing the matrix elements, you first have to specify how the two- and four-index objects will be represented. By default, HORTON uses a dense matrix representation, which is implemented in the DenseLinalgFactory class. An instance of this class must be passed to the methods that compute the matrix elements. Alternatively, you may also use the CholeskyLinalgFactory, which represents all four-index objects with a Cholesky decomposition. Note that the four-center electron-repulsion integrals are computed with LibInt [valeev2014].

This is a basic example, assuming some of the preceding code has created the obasis object:

data/examples/hf_dft/rhf_water_dense.py, lines 22–29¶

# 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)

For the nuclear attraction integrals, you also have to specify arrays with atomic coordinates and nuclear charges.