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

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¶

# Create alpha orbitals
exp_alpha = lf.create_expansion()

# Initial guess
guess_core_hamiltonian(olp, kin, na, exp_alpha)

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

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

# Load fchk file
mol = IOData.from_file('water.fchk')

# Print the number of alpha orbitals (occupied and virtual)
print mol.exp_alpha.nfn

# 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.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 relativey 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:

  data/examples/hf_dft/rhf_water_dense.py, lines 50–52¶

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

• Usage in the unrestricted case:

  data/examples/hf_dft/uhf_methyl_dense.py, lines 47–49¶

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

  data/examples/hf_dft/rks_water_lda.py, lines 52–58¶

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

  data/examples/hf_dft/uks_methyl_lda.py, lines 53–60¶

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

• Usage in the restricted case:

  data/examples/hf_dft/rks_water_lda.py, lines 60–67¶

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

  data/examples/hf_dft/uks_methyl_lda.py, lines 62–71¶

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

  data/examples/hf_dft/rks_water_lda.py, lines 68–80¶

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

  data/examples/hf_dft/uks_methyl_lda.py, lines 72–86¶

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

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.

from horton import *
import numpy as np


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

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.