Wrong RDM evaluations for very small systems #3
Description
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