Wrong RDM evaluations for very small systems

[yannra]

Hi all, I stumbled over a slightly weird behaviour pointing to a bug in either the interface or libgnme itself. Essentially, for very small systems/active spaces, the returned RDMs are all zero. The following is a slight modification to the example from the repository, where however the system size is decreased to 4 H atoms, and the same active space (spanning the complete orbital space in this instance) is considered for the bra and the ket state (allowing for a comparison to the pySCF RDM evaluation).

import numpy as np
from pyscf import gto, scf, mcscf, ao2mo
from pygnme import wick, utils

def owndata(x):
    # CARMA requires numpy arrays to have data ownership
    if not x.flags["OWNDATA"]:
        y = np.zeros(x.shape, order="C")
        y[:] = x
        x = y
    assert x.flags["OWNDATA"]
    return x

n_h_atoms = 4

mol = gto.Mole()
mol.atom = ";".join(["H %d 0 0" % i for i in range(n_h_atoms)])
mol.basis = "sto-6g"
mol.verbose = 0
mol.unit = "Bohr"
mol.build()

mf = scf.RHF(mol)
mf.kernel()

h1e = owndata(mf.get_hcore())
h2e = owndata(ao2mo.restore(1, mf._eri, mol.nao).reshape(mol.nao**2, mol.nao**2))
ovlp = owndata(mf.get_ovlp())
nmo, nocc = mf.mo_occ.size, np.sum(mf.mo_occ > 0)

# Choose an active space which is the complete orbital space
mc1 = mcscf.CASCI(mf, nmo, np.sum(mol.nelec))
e1, _, ci, mo1, _ = mc1.kernel()

ci = owndata(ci)
mo = owndata(mo1)
nact = mc1.ncas
ncore = mc1.ncore

# Compute coupling terms
rdm1_ = np.zeros((nmo, nmo))
rdm2_ = np.zeros((nmo, nmo, nmo, nmo))
# Setup biorthogonalised orbital pair
refx = wick.reference_state[float](nmo, nmo, nocc, nact, ncore, mo)
refw = wick.reference_state[float](nmo, nmo, nocc, nact, ncore, mo)

# Setup paired orbitals
orbs = wick.wick_orbitals[float, float](refx, refw, ovlp)

# Setup matrix builder object
mb = wick.wick_rscf[float, float, float](orbs, mol.energy_nuc())

vx = utils.fci_bitset_list(nocc - ncore, nact)
vw = utils.fci_bitset_list(nocc - ncore, nact)

# Loop over FCI occupation strings
for iwa in range(len(vw)):
    for iwb in range(len(vw)):
        for ixa in range(len(vx)):
            for ixb in range(len(vx)):
                # Compute RDM1 contribution for this pair of determinants
                rdm1_tmp = np.zeros((orbs.m_nmo, orbs.m_nmo))
                rdm2_tmp = np.zeros((orbs.m_nmo * orbs.m_nmo, orbs.m_nmo * orbs.m_nmo))
                mb.evaluate_rdm12(
                    vx[ixa],
                    vx[ixb],
                    vw[iwa],
                    vw[iwb],
                    np.array(1.0),
                    rdm1_tmp,
                    rdm2_tmp,
                )
                rdm1_ += rdm1_tmp * ci[iwa, iwb] * ci[ixa, ixb]
                rdm2_ += rdm2_tmp.reshape(rdm2_.shape) * ci[iwa, iwb] * ci[ixa, ixb]

rdm1_compare, rdm2_compare = mc1.fcisolver.make_rdm12(ci, nact, mol.nelec)

print("Max discrepancy 1-RDM ", abs(rdm1_ - rdm1_compare).max())

print("Max discrepancy 2-RDM ", abs(rdm2_ - rdm2_compare).max())

print("Computed 1-RDM\n", rdm1_)

For four Hydrogen atoms, the script gives an all zero 1RDM but the correct 2RDM, if I decrease the number of H atoms to two, also the 2RDM is all zero.

Maybe you guys can have a quick check what is going on here. I already had a brief chat with @obackhouse about this and he reckons the issue is in libgnme where some of the evaluations return early (incorrectly) due to break statements which are applied for small numbers of orbitals. Either way, I'm raising this here as I haven't tested the setup directly within a direct call to libgnme not going via this python wrapper.

Thanks a lot and best wishes, Yannic

[ghb24]

Tagging @hgaburton here, in case he has any ideas...?

[hgaburton]

Hmmm... It looks like P1 and P2 are correctly evaluated in the libgnme::wick_rscf::evaluate_rdm12 routine (you should be able to verify with print statements at the end of this routine).

e.g. for reference determinant in the H2 example...

P1 # evaluated in libgnme::wick_rscf::evaluate_rdm12
   2.0000        0
        0        0
P2 # evaluated in libgnme::wick_rscf::evaluate_rdm12
   2.0000        0        0        0
        0        0        0        0
        0        0        0        0
        0        0        0        0
rdm1_tmp
[[0. 0.]
 [0. 0.]]

This suggests that the error lies somewhere in the python wrapping?

[ghb24]

Would make sense - might be worth a look @obackhouse ....?

[obackhouse]

I spent some time looking into this, evaluate_rdm1 also has the problem, and therefore I think it's something internal in carma (probably my misuse rather than their bug). Their documentation isn't amazing so it might not be an easy fix but I will keep looking