PySCF: How to do CCSD calculations from an FCIDUMP file in HDF5 format?

Question

Recently my question: Why does CCSD work, but not CCSD(T), after reading an FCIDUMP in PySCF? lead to a bug fix that allows users to do CCSD(T) calculations from an FCIDUMP file in PySCF. That was with the FCIDUMP file being in ASCII format though!

When the FCIDUMP file is in ASCII format, it's easy to create an "SCF object" from the dump file:

mf = pyscf.tools.fcidump.to_scf(name+'.fcidump', molpro_orbsym=False, mf=None)
mf.run()

That object is what's required to do CCSD calculations:

mcc = cc.CCSD(mf)
e = mcc.kernel()

However the function that "loads" an FCIDUMP file in HDF5 format gives an eri object instead of an scf object:

eri = ao2mo.load(name+'.fcidump.h5',dataname='dummy_label')

Attempting to do CCSD from only the ERI object

If I were just to replace the above mf with eri as in the line below:

mcc = cc.CCSD(eri)

then it will trigger the convert_to_rhf function in pyscf/scf/addons.py which will crash at assert (isinstance(mf, hf.SCF)) presumably because eri is not an instance of an hf.SCF object.

By simply doing:

mcc = cc.CCSD()

there will be a complaint about the required mf input being missing. If we do:

mcc = pyscf.cc.ccsd.kernel(eris=eri, callback=None)

there will be an error because cc.ccsd.kernel() requries a mycc object as an input argument. The following gets us closer:

mcc = pyscf.cc.ccsd.kernel(0,eris=eri, verbose=9,callback=None)

but ends in pyscf/cc/ccsd.py at mycc.get_init_guess(eris) presumably because the integer 0 doesn't contain information about the initial guess for the CCSD calculation:

AttributeError: 'int' object has no attribute 'get_init_guess'

In most programs, asking for a CCSD calculation automatically triggers an MP2 calculation so that there's a decent initial guess for the CCSD, and I'm sure that's also happening in PySCF if I do cc.CCSD(mf), but it seems that if we run pyscf.cc.ccsd.kernel() instead (due to not having an mf object), then we need to do explicitly tell the program to do an MP2 calcualtion first. The brings us to the question:

Attempting to do MP2 from only the ERI object

Running mf.MP2().run() requires an mf object. The following:

eri.MP2().run(verbose=9)

won't work because eri is a load object and doesn't have an mp2 attribute:

AttributeError: 'load' object has no attribute 'MP2'

In fact pyscf.ci.cisd.CISD() would also complain about not having an mf object.

Attempting to get an SCF object from an ERI object

As long as we can get an mf object from the eri, we can do CISD, MP2, CCSD and much more. Attempting to create the mf object from the FCIDUMP in the way that we did when the FCIDUMP was in ASCII format:

mf = pyscf.tools.fcidump.to_scf(name+'.fcidump.h5', molpro_orbsym=False, mf=None)

just leads to a complaint which seems to be about the FCIDUMP file not being in ASCII format:

python3.8/codecs.py", line 322, in decode
    (result, consumed) = self._buffer_decode(data, self.errors, final)
UnicodeDecodeError: 'utf-8' codec can't decode byte 0x89 in position 0: invalid start byte

There's not much documentation for pyscf.tools.fcidump.to_scf() but the source code suggests that it will only work if the FCIDUMP file is in ASCII format. Perhaps it can be hacked to go into the HDF5 file and fetch the integrals, but I would think that since PySCF has the capability to load the integrals with ao2mo.load(), it would have the ability to use loaded integrals, meaning that there's already a way to get the mf object from an eri object.

I've also looked at pyscf.tools.fcidump.scf_from_fcidump() but it seems to require the user to already have an mf object. Also in pyscf/tools/fcidump.py there are the from_mo() and from_integrals(), but they are for writing FCIDUMP files, not loading them.

Question

My question was how to do CCSD from an FCIDUMP file in HDF5 format. Based on the above troubleshooting, I've discovered that it boils down to converting the eri object obtained from the FCIDUMP.h5 file into an scf object, which I don't know how to do!

The following script will create the FCIDUMP.h5 file (only about 11 kB since I used a STO-3G basis set) and then will load the FCIDUMP.h5 file into an eri object which is classified as a "load" object:

#!/usr/bin/env python
# Author: Nike Dattani, [email protected]

import pyscf
from pyscf import gto, scf, ao2mo, cc

mol = gto.Mole()
mol.atom = ''' O
               H  1  0.9576257
               H  1  0.9576257  2 104.51 '''
mol.unit = 'angstrom'
mol.basis = 'STO-3G'
mol.charge = 0
mol.spin = 0
mol.verbose = 9
mol.symmetry = True
mol.symmetry_subgroup = 'C2v'
name = 'out'
#mol.output = name+'.txt' # only worked up to end of SCF (CC output went to stdout)
mol.build()

#####################
## Hartree-Fock:
#####################

mf = mol.RHF().set(chkfile = name + '.chk')
mf.kernel()

pyscf.ao2mo.kernel(mol, mf.mo_coeff, name+'.fcidump.h5', dataname='hdf5_group_name')
eri = ao2mo.load(name+'.fcidump.h5',dataname='dummy_label')

Running the above script will also give you an mf object, but the question is how to continue to post-HF calculations without using it at all.

Answer 1

Thanks for all these troubleshooting. As you noticed, ao2mo.kernel does not give "an FCIDUMP file in HDF5 format", but only provides 2-electron integrals in MO basis. Similarly ao2mo.load only loads 2-e integrals. You will need more information, such as 1-e integrals, to re-create a mf object.

If the FCIDUMP file in HDF5 format, instead of ASCII, is desired, basically you can change the relevant functions in tools/fcidump.py to enable it. PySCF currently doesn't have a HDF5 counterpart for the FCIDUMP read/write functions. One reason is that there are no well-adopted standards to store the information in HDF5. Some programs may have their internal format, but this limits the transferability between programs. The text-based FCIDUMP format, on the other hand, has been adopted in the community for many years.

As a temporary solution, I wrote a script to allow write and read HDF5 in a similar fashion to how you treat FCIDUMP. Please feel free to use/modify/optimize. Eventually, such codes can be merged into PySCF.

#!/usr/bin/env python

import h5py, numpy
import pyscf
from pyscf import gto, scf, ao2mo, cc, tools
from pyscf.tools.fcidump import TOL, DEFAULT_FLOAT_FORMAT, MOLPRO_ORBSYM

#####################
## Redefine functions to enable HDF5 format (these functions mirror those in pyscf/tools/fcidump.py)
#####################

# Write some basic information to HDF5
def write_head(fout, nmo, nelec, ms=0, orbsym=None):
    if not isinstance(nelec, (int, numpy.number)):
        ms = abs(nelec[0] - nelec[1])
        nelec = nelec[0] + nelec[1]
    fout['norb'] = nmo
    fout['nelec'] = nelec
    fout['ms2'] = ms
    if orbsym is not None and len(orbsym) > 0:
        fout['orbsym'] = numpy.asarray(orbsym) 
    else:
        fout['orbsym'] = numpy.full(1, nmo) 
    fout['isym'] = 1

# Write 2-e integrals to HDF5
def write_eri(fout, eri, nmo, tol=TOL, float_format=DEFAULT_FLOAT_FORMAT):
    if eri.size == nmo**4:
        eri = ao2mo.restore(8, eri, nmo)
    fout.create_dataset('h2', data=eri)

# Write 1-e integrals to HDF5
def write_hcore(fout, h, nmo, tol=TOL, float_format=DEFAULT_FLOAT_FORMAT):
    h = h.reshape(nmo,nmo)
    fout.create_dataset('h1', data=h)

# Convert the given 1-electron and 2-electron integrals to HDF5 format
def hdf5_from_integrals(filename, h1e, h2e, nmo, nelec, nuc=0, ms=0, orbsym=None,
                        tol=TOL, float_format=DEFAULT_FLOAT_FORMAT):
    with h5py.File(filename, 'w') as fout:
        write_head(fout, nmo, nelec, ms, orbsym)
        write_eri(fout, h2e, nmo, tol=tol, float_format=float_format)
        write_hcore(fout, h1e, nmo, tol=tol, float_format=float_format)
        fout['nuc'] = nuc

# Use the Hamiltonians stored in HDF5 format to build an SCF object
def hdf5_to_scf(filename, molpro_orbsym=MOLPRO_ORBSYM, mf=None, **kwargs):
    with h5py.File(filename, 'r') as f:
        mol = gto.M()
        mol.nelectron = int(f['nelec'][()])
        mol.spin = int(f['ms2'][()])
        norb = mol.nao = int(f['norb'][()])
        nuc = f['nuc'][()]
        mol.energy_nuc = lambda *args: nuc 
        mol.incore_anyway = True
        if 'orbsym' in f.keys():
            mol.symmetry = True
            mol.groupname = 'N/A'
            orbsym = f['orbsym'][()]
            mol.irrep_id = list(set(orbsym))
            mol.irrep_name = [('IR%d' % ir) for ir in mol.irrep_id]
            so = numpy.eye(norb)
            mol.symm_orb = []
            for ir in mol.irrep_id:
                mol.symm_orb.append(so[:,orbsym==ir])

        if mf is None:
            mf = mol.RHF(**kwargs)
        else:
            mf.mol = mol
        h1 = f['h1'][()]
        mf.get_hcore = lambda *args: h1
        mf.get_ovlp = lambda *args: numpy.eye(norb)
        
        norb_pair = norb * (norb+1) // 2
        h2e = f['h2'][()]
        mf._eri = h2e[numpy.tril_indices(norb_pair)]

    return mf


#####################
## Input starts here
#####################

mol = gto.Mole()
mol.atom = ''' O
               H  1  0.9576257
               H  1  0.9576257  2 104.51 '''
mol.unit = 'angstrom'
mol.basis = 'CC-PVTZ'
mol.charge = 0
mol.spin = 0
mol.verbose = 9
mol.symmetry = True
mol.symmetry_subgroup = 'C2v'
name = 'out'
mol.build()

#####################
## Hartree-Fock:
#####################

mf = mol.RHF().set(chkfile = name + '.chk')
mf.kernel()

# Replace two functions in tools/fcidump.py if you want to store information as HDF5
store_as_h5 = True
if store_as_h5: 
    from pyscf.tools import fcidump
    fcidump.from_integrals = hdf5_from_integrals
    fcidump.to_scf = hdf5_to_scf
    ext = 'h5'
else:
    ext = 'fcidump'

tools.fcidump.from_chkfile(f'{name}.{ext}',name+'.chk', tol=1e-18, float_format=' %.16g', molpro_orbsym=False, orbsym=None)
mf = tools.fcidump.to_scf(f'{name}.{ext}', molpro_orbsym=False, mf=None)
mf.kernel()

#####################
## post-Hartree-Fock:
#####################

mcc = cc.CCSD(mf).set(frozen=1,verbose=9)
e = mcc.kernel()
et = mcc.ccsd_t()