Molecular Hamiltonian as LinearMap?

[msmerlak]

Hi,

I would have a feature request. (Would love to code it myself and PR, but, not being an expert, this would take me several days. Hoping you can help!)

Fermi runs HF with a simple interface. I would like to use the molecular orbitals and ERIs to form the molecular Hamiltonian as a LinearMap, for use in FCI calculations. Right now I’m using a workaround with ‘PyCall’ and ‘PySCF’. A pure Julia version would be more efficient and more elegant.

Thanks!

[gustavojra]

Hi @msmerlak,

I have worked in the past with an approach where I'd simply build the whole Hamiltonian matrix... It is expensive, but it's fairly simple to implement, I must still have the code somewhere.

Would that work for you, or do you need something more efficient like Olsen's fci (where you don't compute H, but rather just the effect of it onto the CI vector).

[msmerlak]

Hi @gustavojra, and thanks for your trouble. I would indeed need the matrix-free version, as building the whole matrix would be too expensive for all but the smallest molecules. (My goal is to show that an eigenvalue algorithm I've designed can outperform Davidson. But to show this I need an example where H * CI is the computational bottleneck, i.e. big enough).

Using PySCF this is how I do it:

from pyscf import gto, scf, fci, ao2mo
from functools import reduce
from numpy import dot

def hop(coordinates, basis):

    mol = gto.M(
        atom = coordinates,
        basis = basis, 
        symmetry = True)

    hf = scf.HF(mol)
    hf.kernel()
    nelec = mol.nelectron
    norb = hf.mo_coeff.shape[0]

    neleca, nelecb = fci.addons._unpack_nelec(nelec)

    na = fci.cistring.num_strings(norb, neleca)
    nb = fci.cistring.num_strings(norb, nelecb)


    h1e = reduce(dot, (hf.mo_coeff.T, hf.get_hcore(), hf.mo_coeff))
    eri = ao2mo.incore.general(hf._eri, (hf.mo_coeff,)*4, compact=False)
    h2e = fci.direct_spin1.absorb_h1e(h1e, eri, norb, nelec, .5)

    fci_solver = fci.FCI(hf)

    return lambda c: fci.direct_spin1.contract_2e(h2e, c.reshape(na,nb), norb, nelec).ravel(), na, nb, fci.direct_spin1.make_hdiag(h1e, eri, norb, nelec), hf, fci_solver

I've tried to do this myself in Julia using Fermi.jl, but got confused with the construction of the 2e Hamiltonian and the contractions that go into it.