2.4. Dumping/Loading a Hamiltonian to/from a file¶
Before reading this section, please refer to Input/Output Operations: the IOData Container
for more informations on how to load and dump data in different data file
formats with PyBEST. For reasons of completeness, the instruction on how to
dump/load a Hamiltonian in PyBEST are summarized here. Nonetheless, we recommend
that the user should know the behavior of and basic operations on the
IOData class before continuing reading.
PyBEST supports two formats for Hamiltonians: (i) an internal binary format based
on HDF5 (extension .h5) and (ii) Molpro’s FCIDUMP ASCII format (containing FCIDUMP
somewhere in the file name). The internal format is more flexible and can store
Hamiltonians in various ways, for instance, all one-electron integrals can be
stored separately. The FCIDUMP format is more restricted
but can be used to interface PyBEST with different codes, e.g., Molpro or
Dalton (both can generate FCIDUMP files).
2.4.1. PyBEST’s internal format¶
2.4.1.1. Dumping¶
You can store all one- and two-electron integrals in the atomic-orbital (AO) basis. For example, to store the molecular Hamiltonian of the water molecule, follow the code snippet below
from pybest import context
from pybest.gbasis import get_gobasis
from pybest.gbasis.gobasis import (
compute_kinetic,
compute_nuclear,
compute_cholesky_eri,
compute_nuclear_repulsion,
)
from pybest.io.iodata import IOData
from pybest.linalg.cholesky import CholeskyLinalgFactory
# Set up molecule, define basis set
# ---------------------------------
# get the XYZ file from PyBEST's test data directory
fn_xyz = context.get_fn("test/water.xyz")
obasis = get_gobasis("cc-pvdz", fn_xyz)
lf = CholeskyLinalgFactory(obasis.nbasis)
# Construct Hamiltonian
# ---------------------
mol = IOData.from_file()
mol.lf = lf
mol.kin = compute_kinetic(obasis)
mol.ne = compute_nuclear(obasis)
mol.eri = compute_cholesky_eri(obasis)
mol.e_core = compute_nuclear_repulsion(obasis)
# Write to a HDF5 file
# --------------------
mol.to_file("hamiltonian_ao.h5")
The internal format will store all Hamiltonian terms (that is, the integrals and core
energy) as attributes of the IOData instance.
The IOData container class is rather flexible
and allows you to update the container on the fly. To delete or add attributes,
use e.g. del mol.title or mol.egg = "spam" (see also
Filling the IOData container for how to
manipulate the container).
In the HDF5 file, all data is stored in binary form with full precision, which means that all significant digits are written to file. In the example above, the dumped HDF5 file will have the following layout:
$ h5dump -n hamiltonian_ao.h5
HDF5 "hamiltonian_ao.h5" {
FILE_CONTENTS {
group /
dataset /e_core
group /eri
dataset /eri/array
group /kin
dataset /kin/array
group /lf
group /ne
dataset /ne/array
}
}
2.4.1.2. Loading¶
You can load one- and two-electron integrals stored in PyBEST’s internal format, as follows:
from pybest.io.iodata import IOData
# Load the integrals from the file
# Assuming you have run `dump_internal_ao.py` first.
data = IOData.from_file("hamiltonian_ao.h5")
# Print all attributes (stored as dictionary)
print(data.__dict__)
2.4.2. FCIDUMP format¶
2.4.2.1. Dumping¶
While the internal (binary) format can only be handled by PyBEST, the FCIDUMP format stores the Hamiltonian in ASCII format and can be easily read in by other codes. In the current version of PyBEST, the FCIDUMP format has the following restrictions:
One-electron integrals (e.g., the kinetic energy and nuclear attraction) are all stored as one single term (FCIDUMP standard).
Integrals can only be stored in a restricted (MO) basis set, that is, unrestricted orbitals (alpha and beta orbitals) are not supported, yet.
The two-electron integrals must be stored as a
DenseFourIndexobject, that is, the Cholesky decomposition of the ERI is not supported.
The FCIDUMP format is generally used for storing integrals in the MO basis.
Thus, the one- and two-electron integrals need to be first transformed from the
AO to some MO basis. In the following example, we will assume that a set of
one- (one) and two-electron integrals (two) are already available.
See section The AO/MO Transformation of the Hamiltonian for a detailed instruction on
how to transform (and store) integrals into the MO basis.
Furthermore, we have an additional energy associated with some external potential
labeled as e_core.
To store all (pre-transformed) integrals using the FCIDUMP format, proceed as
follows
# Write to a FCIDUMP file
# -----------------------
data = IOData(one=one, two=two, e_core=e_core, nelec=20, ms2=0)
data.to_file("hamiltonian_mo.FCIDUMP")
In this example, we set the IOData attributes by
using keyword arguments in the constructor. The complete example is also shown
below in Example Python scripts.
Note
Although the FCIDUMP format stores all one-electron integrals as one term, PyBEST can handle them separatly. The FCIDUMP dump function will combine them into one single term automatically.
The file hamiltonian_mo.FCIDUMP will contain the following:
&FCI NORB=28,NELEC=20,MS2=0,
ORBSYM= 1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,
ISYM=1
&END
3.1025487346395235e+00 1 1 1 1
-2.0955127330533898e-06 2 1 1 1
2.9025424092728871e+00 2 1 2 1
3.1025489371112021e+00 2 2 1 1
2.0955125859036790e-06 2 2 2 1
3.1025491395905633e+00 2 2 2 2
-3.4893251387724522e-01 3 1 1 1
.
.
-2.4806545706468341e-16 28 26 0 0
-4.1199682554449168e-16 28 27 0 0
-6.9042050091456133e+00 28 28 0 0
2.0000000009826579e+01 0 0 0 0
This file is divided into two blocks. The first block (between &FCI and
&END) contains system-specific information:
- NORB
number of orbitals/basis functions
- NELEC
number of electrons
- MS2
spin projection
- ORBSYM
irreducible representation of each orbital
- ISYM
total symmetry of the wave function
The second block (after &END) contains the one- and two-electron integrals
and the core energy in the following order
All symmetry-unique elements of the two-electron integrals are listed, where the first column is the value of the integral, followed by the orbital indices. Note that the orbital indices start from 1 (Fortran indexing). The last four columns encode the orbital indices (
i j k l) in chemists’ notation,
\[(ij\vert kl) = \langle ik \vert jl \rangle = \int \phi_i^*(\mathbf{x}_1) \phi_k^*(\mathbf{x}_2) \frac{1}{r_{12}} \phi_j(\mathbf{x}_1) \phi_l(\mathbf{x}_2) d\mathbf{x}_1 d\mathbf{x}_2\]
Then, all symmetry-unique elements of the one-electron integrals are listed, where the first column is the value of the integral, followed by the orbital indices. Note that the last two columns contain zeros.
The core energy (for instance, the nuclear repulsion term, etc.) is written on the last line with all orbital indices equal to 0.
PyBEST does not (yet) support point group symmetries (only \(C_1\)). Thus, all orbitals will be assigned the symmetry label 1, while the total symmetry of the wave function will be set to 1.
2.4.2.2. Loading¶
You can load integrals, stored in an FCIDUMP file, as follows
data = IOData.from_file("hamiltonian_mo.FCIDUMP")
# Access some attributes. In more realistic cases, some code follows that does
# a useful calculation.
print("Core energy: ", data.e_core)
print("Element [0, 0] of one-body integrals: ", data.one.get_element(0, 0))
The resulting IOData container data has the
following attributes
- one
All one-electron integrals combined in one single term
- two
The two-electron repulsion integrals
- e_core
The core energy
- nelec
The number of electrons comprised in the active space
- ms2
The spint projection
- lf
An instance of the
DenseLinalgFactory- orb_a
An instance of
DenseOrbital(created under the assumption that the orbitals are orthonormal)- olp
An overlap matrix (created under the assumption that the orbitals are orthonormal)
These attributes can be then passed individually to some wave function method.
2.4.3. Example Python scripts¶
Several complete examples can be found in the directory data/examples/hamiltonian.
2.4.3.1. Dumping in the FCIDUMP format: the Ne dimer¶
This is a basic example on how to dump all one- and two-electron integrals in the MO basis for the Neon dimer using the FCIDUMP format in PyBEST. This script performs a RHF calculation followed by a transfomration of the Hamiltonian into the MO basis for the Neon dimer using the cc-pVDZ basis set.
import numpy as np
from pybest.gbasis import get_gobasis
from pybest.gbasis.gobasis import (
compute_eri,
compute_kinetic,
compute_nuclear,
compute_nuclear_repulsion,
compute_overlap,
)
from pybest.io.iodata import IOData
from pybest.linalg.dense import DenseLinalgFactory
from pybest.orbital_utils import transform_integrals
from pybest.scf import AufbauOccModel
from pybest.wrappers.hf import RHF
# Set up Neon dimer, define basis set
# -----------------------------------
coordinates = np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 5.0]])
atom = np.array(["Ne", "Ne"])
mol = IOData(coordinates=coordinates, atom=atom)
obasis = get_gobasis("cc-pvdz", mol)
lf = DenseLinalgFactory(obasis.nbasis)
occ_model = AufbauOccModel(2)
orb = lf.create_orbital()
# Construct Hamiltonian
# ---------------------
kin_ao = compute_kinetic(obasis)
ne_ao = compute_nuclear(obasis)
eri_ao = compute_eri(obasis)
e_core = compute_nuclear_repulsion(obasis)
olp = compute_overlap(obasis)
# Get Hartree-Fock MOs
# --------------------
hf = RHF(lf, occ_model)
hf_output = hf(kin_ao, ne_ao, eri_ao, e_core, orb, olp)
# Transform Hamiltinian to MO basis
# ---------------------------------
# transform integrals for restricted orbitals orb
t_ints = transform_integrals(kin_ao, ne_ao, eri_ao, orb)
# transformed one-electron integrals: attribute 'one' (list)
(one,) = t_ints.one # or: one = ti_.one[0]
# transformed two-electron integrals: attribute 'two' (list)
(two,) = t_ints.two # or: two = ti_.two[0]
# Write to a FCIDUMP file
# -----------------------
data = IOData(one=one, two=two, e_core=e_core, nelec=20, ms2=0)
data.to_file("hamiltonian_mo.FCIDUMP")