Dumping Cluster Hamiltonians

Some users may want to utilize Vayesta to easily define fragments within the system and obtain the corresponding cluster Hamiltonians, but solve the embedding problems externally and with their own solvers. To accomodate for this, the EWF class allows setting solver='Dump', which will dump orbitals and integrals of all fragments to an HDF5 file and exit.

The name of the HDF5 file is clusters.h5 by default, but can be adjusted via an additional solver option:

    mf, solver="DUMP", bath_options=dict(threshold=1e-6), solver_options=dict(dumpfile="clusters-rhf.h5")
)

Note

The full example can be found at examples/ewf/molecules/20-dump-clusters.py

The dump file contains a separate group for each fragment which was defined for the embedding. The content of each group can be best illustrated via this code snippet:

import h5py

with h5py.File("clusters-rhf.h5", "r") as f:
    # The HDF5 file contains a separate group for each fragment in the system:
    for key, frag in f.items():
        # The HDF5-group key for each fragment is constructed as 'fragment_%d' % id
        print("\nKey= %s" % key)
        # Name and ID of fragment:
        print("name= %s, id= %d" % (frag.attrs["name"], frag.attrs["id"]))
        # Number of all/occupied/virtual orbitals:
        print("Full cluster:")
        norb, nocc, nvir = frag.attrs["norb"], frag.attrs["nocc"], frag.attrs["nvir"]
        print("norb= %d, nocc= %d, nvir= %d" % (norb, nocc, nvir))
        # Orbital coefficients:
        # The first dimension corresponds to the atomic orbitals,
        # the second dimension to the cluster or fragment orbitals, respectively.
        print("c_cluster.shape= (%d, %d)" % frag["c_cluster"].shape)
        print("c_frag.shape=    (%d, %d)" % frag["c_frag"].shape)
        # Integral arrays:
        # fock, and eris are the Fock, and 2-electron Hamiltonian
        # matrix elements in the cluster space.
        # heff is the 1-electron Hamiltonian, plus an additional potential due to the Coulomb and exchange
        # interaction with the occupied environment orbitals outside the cluster.
        print("heff.shape=      (%d, %d)" % frag["heff"].shape)
        print("fock.shape=      (%d, %d)" % frag["fock"].shape)
        # The 2-electron integrals are in chemical ordering: eris[i,j,k,l] = (ij|kl)
        print("eris.shape=      (%d, %d, %d, %d)" % frag["eris"].shape)
        # DMET cluster:
        print("DMET cluster:")
        norb, nocc, nvir = [frag.attrs["%s_dmet_cluster" % x] for x in ("norb", "nocc", "nvir")]
        print("norb= %d, nocc= %d, nvir= %d" % (norb, nocc, nvir))

For a spin-unrestricted calculation, the shapes and dataset names are slighlty different:

with h5py.File("clusters-uhf.h5", "r") as f:
    for key, frag in f.items():
        print("\nKey= %s" % key)
        # Name and ID:
        print("name= %s, id= %d" % (frag.attrs["name"], frag.attrs["id"]))
        # The orbital sizes are now arrays of length 2, representing alpha and beta spin dimension
        norb, nocc, nvir = frag.attrs["norb"], frag.attrs["nocc"], frag.attrs["nvir"]
        print("norb= (%d, %d), nocc= (%d, %d), nvir= (%d, %d)" % (*norb, *nocc, *nvir))
        # The suffixes _a and _b correspond to alpha and beta spin:
        print("c_cluster_a.shape= (%d, %d)" % frag["c_cluster_a"].shape)
        print("c_cluster_b.shape= (%d, %d)" % frag["c_cluster_b"].shape)
        print("c_frag_a.shape=    (%d, %d)" % frag["c_frag_a"].shape)
        print("c_frag_b.shape=    (%d, %d)" % frag["c_frag_b"].shape)
        # Integral arrays:
        print("heff_a.shape=      (%d, %d)" % frag["heff_a"].shape)
        print("heff_b.shape=      (%d, %d)" % frag["heff_b"].shape)
        print("fock_a.shape=      (%d, %d)" % frag["fock_a"].shape)
        print("fock_b.shape=      (%d, %d)" % frag["fock_b"].shape)
        # The spin blocks correspond to (aa|aa), (aa|bb), and (bb|bb)
        # (Note that (bb|aa) = (aa|bb).transpose(2,3,0,1) for real orbitals):
        print("eris_aa.shape=     (%d, %d, %d, %d)" % frag["eris_aa"].shape)
        print("eris_ab.shape=     (%d, %d, %d, %d)" % frag["eris_ab"].shape)
        print("eris_bb.shape=     (%d, %d, %d, %d)" % frag["eris_bb"].shape)
        # DMET cluster:
        print("DMET cluster:")
        norb, nocc, nvir = [frag.attrs["%s_dmet_cluster" % x] for x in ("norb", "nocc", "nvir")]
        print("norb= (%d, %d), nocc= (%d, %d), nvir= (%d, %d)" % (*norb, *nocc, *nvir))
        print("c_dmet_cluster_a.shape=    (%d, %d)" % frag["c_dmet_cluster_a"].shape)
        print("c_dmet_cluster_b.shape=    (%d, %d)" % frag["c_dmet_cluster_b"].shape)