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
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')
# CODE BELOW IS FOR horton-regression-test.py ONLY. IT IS NOT PART OF THE EXAMPLE.
rt_results = {'coordinates': mol.coordinates}
# BEGIN AUTOGENERATED CODE. DO NOT CHANGE MANUALLY.
rt_previous = {
'coordinates': np.array([
2.3071525963747455, 0.0, 0.0, 1.9409002724808868, 1.2473408656988467, 0.0,
0.95842582582035096, 2.0986598198277173, 0.0, -0.32834204862486149,
2.2836691095829882, 0.0, -1.5108636425857356, 1.7436295926805356, 0.0,
-2.2136967052780125, 0.65000000000000002, 0.0, -2.213696705278013,
-0.64999999999999936, 0.0, -1.5108636425857362, -1.7436295926805354, 0.0,
-0.32834204862486199, -2.2836691095829882, 0.0, 0.95842582582035007,
-2.0986598198277178, 0.0, 1.9409002724808868, -1.2473408656988465, 0.0
]).reshape(11, 3),
}
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):
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)
# CODE BELOW IS FOR horton-regression-test.py ONLY. IT IS NOT PART OF THE EXAMPLE.
rt_results = {
'scales': obasis.get_scales()
}
# BEGIN AUTOGENERATED CODE. DO NOT CHANGE MANUALLY.
rt_previous = {
'scales': np.array([
0.013401011981382846, 0.018048783700127666, 0.024308506962500243,
0.032739242741423244, 0.044093946902429959, 0.059386717304072578,
0.079983363702002319, 0.10772338932847644, 0.14508427842131391,
0.1954027623550548, 0.26317282583235962, 0.35444706831084288, 0.47737726658062746,
0.64294241657582896, 0.86592927642596351, 1.1662529838442486, 1.5707356932652521,
2.1155020842604459, 2.8492056860355119, 3.8373741636728487, 5.168261647165739,
6.9607307795074735, 9.3748684359572785, 12.626283213000068, 17.005361607360911,
22.903202670074823, 30.846547380648691, 41.544822312114547, 55.953499094940575,
75.359428365025664
]),
}
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:
# 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.