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:

  1. Setting up the (molecular electronic) Hamiltonian
  2. Generating an initial guess of the wavefunction
  3. Defining the effective Hamiltonian
  4. Choose an orbital occupation scheme
  5. Optimize the wavefunction with a self-consistent field algorithm
  6. Optional (if needed for 7 or 8): Conversion of density and Fock matrix to orbitals
  7. Optional: Writing SCF results to a file
  8. 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:

data/examples/hf_dft/rhf_water_dense.py, lines 31–35
# 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:

data/examples/hf_dft/uhf_methyl_dense.py, lines 27–32

# Create alpha orbitals
exp_alpha = lf.create_expansion()
exp_beta = lf.create_expansion()

# Initial guess

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:

    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.

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:

    data/examples/hf_dft/rhf_water_dense.py, lines 37–45
    # 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:

    data/examples/hf_dft/uhf_methyl_dense.py, lines 35–42
    # 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'),
]

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

    data/examples/hf_dft/rks_water_lda.py, lines 37–47
    # 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'),
]

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

    data/examples/hf_dft/uks_methyl_gga.py, lines 38–48
    # 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'),
]

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

    data/examples/hf_dft/rks_water_hybgga.py, lines 37–46
    # 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'),
]

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

    data/examples/hf_dft/uks_methyl_mgga.py, lines 39–47
    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'),

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

    data/examples/hf_dft/rks_water_hybgga.py, lines 38–45
    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'),

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

    data/examples/hf_dft/uks_methyl_numlda.py, lines 38–48
    # 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'),
]

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)

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)

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

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:

    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)

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

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

  • 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):

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:

  1. 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)

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

\[\tilde{t}_{pq} = t_{pq} + \sum_{i \in \textrm{ncore}} ( 2 \langle pi \vert qi \rangle - \langle pi \vert iq \rangle),\]

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

\[e_\text{core} = e_\text{nn} + 2\sum_{i \in \textrm{ncore}} t_{ii} + \sum_{i, j \in \textrm{ncore}} (2 \langle ij \vert ij \rangle - \langle ij \vert ji \rangle)\]

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.

data/examples/hf_dft/rhf_n2_dense.py

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


# CODE BELOW IS FOR horton-regression-test.py ONLY. IT IS NOT PART OF THE EXAMPLE.
rt_results = {
    'energy': ham.cache['energy'],
    'exp_alpha': exp_alpha.energies,
    'nn': ham.cache["energy_nn"],
    'kin': ham.cache["energy_kin"],
    'ne': ham.cache["energy_ne"],
    'hartree': ham.cache["energy_hartree"],
    'x_hf': ham.cache["energy_x_hf"],
}
# BEGIN AUTOGENERATED CODE. DO NOT CHANGE MANUALLY.
rt_previous = {
    'energy': -108.95408660509983,
    'exp_alpha': np.array([
        -15.686532157922816, -15.683135679259996, -1.4709708458239139,
        -0.77419994645767287, -0.62619794227188197, -0.60802326936155027,
        -0.60802301978559303, 0.17548863834889036, 0.17548869520786309,
        0.59453400073264484, 0.82020385854776667, 0.87230499592573318,
        0.87230519090053693, 0.99220657829527714, 1.0509691768342646, 1.050969299052142,
        1.1439197756891639, 1.6428993715359781, 1.7574964669593922, 1.7574964669593935,
        1.8817291134001792, 1.8817291219136119, 2.2987680720779378, 2.2987680720779395,
        2.8741336841539828, 2.9961172837814947, 2.9961172850210782, 3.28312430564872
    ]),
    'hartree': 74.76419388482387,
    'kin': 108.71316284694453,
    'ne': -302.93789892354596,
    'nn': 23.615376326684874,
    'x_hf': -13.108920740007134,
}

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* ODASCFSolver Basic RKS DFT example with LDA exhange-correlation functional (Dirac+VWN)
uks_methyl_lda.py UKS/6-31G* ODASCFSolver Basic UKS DFT example with LDA exhange-correlation functional (Dirac+VWN)
rks_water_gga.py RKS/6-31G* CDIISSCFSolver Basic RKS DFT example with GGA exhange-correlation functional (PBE)
uks_methyl_gga.py UKS/6-31G* CDIISSCFSolver Basic UKS DFT example with GGA exhange-correlation functional (PBE)
rks_water_hybgga.py RKS/6-31G* EDIIS2SCFSolver Basic RKS DFT example with hyrbid GGA exhange-correlation functional (B3LYP)
uks_methyl_hybgga.py UKS/6-31G* EDIIS2SCFSolver Basic UKS DFT example with hyrbid GGA exhange-correlation functional (B3LYP)
rks_water_mgga.py RKS/6-31G* CDIISSCFSolver Basic RKS DFT example with MGGA exhange-correlation functional (TPSS)
uks_methyl_mgga.py UKS/6-31G* CDIISSCFSolver Basic UKS DFT example with MGGA exhange-correlation functional (TPSS)
rks_water_hybmgga.py RKS/6-31G* CDIISSCFSolver Basic RKS DFT example with hybrid MGGA exhange-correlation functional (TPSS)
uks_methyl_hybmgga.py UKS/6-31G* CDIISSCFSolver Basic UKS DFT example with hybrid MGGA exhange-correlation functional (M05)
rks_water_numlda.py RKS/6-31G* CDIISSCFSolver RKS DFT example with LDA and numerical Hartree
uks_methyl_numlda.py UKS/6-31G* CDIISSCFSolver UKS DFT example with LDA and numerical Hartree
rks_water_numgga.py RKS/6-31G* CDIISSCFSolver RKS DFT example with GGA and numerical Hartree
uks_methyl_numgga.py UKS/6-31G* 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.