Background
My colleague and I have built an interface between PySCF and MRCC which seems to work all the time for RHF references, but only sometimes for UHF references. We use PySCF to calculate the integrals and print them in FCIDUMP format, then we convert the FCIDUMP file into MRCC's fort.55 format.
The format of the FCIDUMP follows the UHF standard designed by Knowles and Handy, in which after the header part, separate 4-index and 2-index integrals are written for alpha-alpha, beta-beta, and alph-beta orbitals, with a line of zeroes separating them. Unfortunately, despite this format being used in MOLPRO, and being compatible with the fort.55 file in some major quantum chemistry codes like MRCC or CFOUR, it is not natively supported in PySCF.
The Problem
We considered two open-shell systems, Li and B++. Both systems have 3 electrons occupying s-type orbitals. For B++, the coupled cluster energies obtained from our PySCF-MRCC interface are identical to the ones obtained using MRCC-MRCC (which is how we will label "standalone MRCC" in which the integrals, SCF and AO->MO transformation are calculated using MRCC, and the post-SCF calculations are also done by MRCC).
B++ (6-31G) (symmetry is on)
mrcc-mrcc
CCSD -23.372809509970
CCSD(T) -23.372810050030
CCSDT -23.372810049178
pyscf-mrcc
CCSD -23.372809509970
CCSD(T) -23.372810050030
CCSDT -23.372810049178
But for Li, this is what we get:
Li (6-31G) (symmetry is on)
mrcc-mrcc
UHF -7.431235814771
CCSD -7.431553874113
CCSD(T) -7.431554248616
CCSDT -7.431554228106
pysc-mrcc
UHF -7.431235814771
CCSD -7.431274408997
CCSD(T) -7.431274447587
CCSDT -7.431274443531
Upon some further debugging, we are finding that the generated fort.55 in the case B++ is similar to the one generated in MRCC in terms of having the same ORBSYM and the same indices for the integrals. However, for the case of Li, the generated fort.55 is having similar ORBSYM with that generated with MRCC but the indices and the number of integrals are not the same. Specifically, the irreducible representationsin of the orbitals:
Li (6-31G) mrcc-mrcc:
'Ag' 'Ag' 'B1u' 'B2u' 'B3u' 'Ag' 'B1u' 'B2u' 'B3u'
Li (6-31G) pyscf-mrcc:
'Ag' 'Ag' 'B1u' 'B2u' 'B3u' 'Ag' 'B1u' 'B2u' 'B3u'
About the integrals, there are 1163 integrals generated from mrcc-mrcc while 1166 integrals from pyscf-mrcc in the case of Li, thus I won't be able to show all of them here but rather the trend of differences. Interested readers are encouraged to refer to the fort.55 generated with pyscf-mrcc and the other fort.55 generated in mrcc-mrcc. By comparing those two files, you will find that the first section of integrals which corresponds to the alpha-alpha 2e integrals have similar indices, for instance:
Li fort.55 mrcc-mrcc
3.314991659E-01 2 2 1 1
2.189901842E-01 3 3 1 1
Li fort.55 pyscf-mrcc
3.314991659E-01 2 2 1 1
2.189901842E-01 3 3 1 1
B++ fort.55 mrcc-mrcc
7.534175484E-01 2 2 1 1
7.056882647E-01 3 3 1 1
B++ fort.55 pyscf-mrcc
7.534175484E-01 2 2 1 1
7.056882647E-01 3 3 1 1
However, in the beta-beta 2e integrals, we start to see differences in the indices and the corresponding coefficients, for instance:
Li fort.55 mrcc-mrcc
1.945413288E-01 5 5 1 1
1.120655588E-01 6 1 1 1
Li fort.55 pyscf-mrcc
1.945413288E-01 5 5 1 1
7.570801776E-02 6 3 1 1
B++ fort.55 mrcc-mrcc
6.73635998E-01 5 5 1 1
-2.84983078E-01 6 3 1 1
B++ fort.55 pyscf-mrcc
6.73635998E-01 5 5 1 1
-2.84983078E-01 6 3 1 1
Similar thing could be observed for alpha-beta integrals as well. In the case of Li, the indices or coefficients are sometimes different between mrcc-mrcc and pyscf-mrcc.
I will point out that the energies actually do match if we turn off symmetry, but then the RAM and CPU time are significantly greater, which means that we won't be able to do very-high-order coupled cluster calculations with big basis sets.
Why is the code described below working for B++ but not for Li?
The code
Suppose that we have defined a Mol object for an open-shell system, followed by running a UHF calculation. In the next step, we convert the 2e integrals using 4-fold permutation symmetry for alpha-alpha, beta-beta, and alpha-beta orbitals which is saved in corresponding eri_aaaa, eri_bbbb, and eri_aabb, respectively. Moreover, the transformed 1e core Hamiltonian for alpha and beta are also calculated and stored in h_aa and h_bb, respectively.
orbs = mf.mo_coeff
nmo = orbs[0].shape[0]
eri_aaaa = pyscf.ao2mo.restore(4,pyscf.ao2mo.incore.general(mf._eri, (orbs[0],orbs[0],orbs[0],orbs[0]), compact=False),nmo)
eri_bbbb = pyscf.ao2mo.restore(4,pyscf.ao2mo.incore.general(mf._eri, (orbs[1],orbs[1],orbs[1],orbs[1]), compact=False),nmo)
eri_aabb = pyscf.ao2mo.restore(4,pyscf.ao2mo.incore.general(mf._eri, (orbs[0],orbs[0],orbs[1],orbs[1]), compact=False),nmo)
h_core = mf.get_hcore(mol)
h_aa = reduce(numpy.dot, (orbs[0].T, h_core, orbs[0]))
h_bb = reduce(numpy.dot, (orbs[1].T, h_core, orbs[1]))
nuc = mol.energy_nuc()
In the next section of the code, we try to handle the symmetry of the orbitals in a way that will match the format of fort.55. We also define the range of values for integral indices:
if mol.symmetry:
groupname = mol.groupname
if groupname in ('SO3', 'Dooh'):
logger.info(mol, 'Lower symmetry from %s to D2h', groupname)
raise RuntimeError('Lower symmetry from %s to D2h' % groupname)
elif groupname == 'Coov':
logger.info(mol, 'Lower symmetry from Coov to C2v')
raise RuntimeError('''Lower symmetry from Coov to C2v''')
orbsym = pyscf.symm.label_orb_symm(mol,mol.irrep_name,mol.symm_orb,orbs[0])
orbsym = numpy.array(orbsym)
orbsym = [param.IRREP_ID_TABLE[groupname][i]+1 for i in orbsym]
a_inds = [i+1 for i in range(orbs[0].shape[0])]
b_inds = [i+1 for i in range(orbs[1].shape[1])]
nelec = mol.nelec
tol=1e-18
The last part is considered the main part of the script, which is supposed to handle the header of the generated file in a way that will match the fort.55 format. It is supposed to deal with the 4-index, and 2-index integrals in which we are assuming that we are using a 4-fold permutation symmetry. The sections after '4-fold symmetry' are in the sequence of alpha-alpha, beta-beta, alpha-beta 2e integrals, followed by alpha, beta 1e integrals, and finally the nuclear repulsion energy.
with open('fort.55', 'w') as fout:
if not isinstance(nelec, (int, numpy.number)):
ms = abs(nelec[0] - nelec[1])
nelec = nelec[0] + nelec[1]
else: ms=0
fout.write(f"{nmo:1d} {nelec:1d}\n")
if orbsym is not None and len(orbsym) > 0:
fout.write(f"{' '.join([str(x) for x in orbsym])}\n")
else:
fout.write(f"{' 1' * nmo}\n")
fout.write(' 150000\n')
output_format = float_format + ' %5d %5d %5d %5d\n'
#4-fold symmetry
kl = 0
for l in range(nmo):
for k in range(0, l+1):
ij = 0
for i in range(0, nmo):
for j in range(0, i+1):
if i >= k:
if abs(eri_aaaa[ij,kl]) > tol:
fout.write(output_format % (eri_aaaa[ij,kl], a_inds[i], a_inds[j], a_inds[k], a_inds[l]))
ij += 1
kl += 1
fout.write(' 0.00000000000000000000E+00' + ' 0 0 0 0\n')
kl = 0
for l in range(nmo):
for k in range(0, l+1):
ij = 0
for i in range(0, nmo):
for j in range(0, i+1):
if i >= k:
if abs(eri_bbbb[ij,kl]) > tol:
fout.write(output_format % (eri_bbbb[ij,kl], b_inds[i], b_inds[j], b_inds[k], b_inds[l]))
ij += 1
kl += 1
fout.write(' 0.00000000000000000000E+00' + ' 0 0 0 0\n')
ij = 0
for j in range(nmo):
for i in range(0, j+1):
kl = 0
for k in range(nmo):
for l in range(0, k+1):
if abs(eri_aabb[ij,kl]) > tol:
fout.write(output_format % (eri_aabb[ij,kl], a_inds[i], a_inds[j], b_inds[k], b_inds[l]))
kl += 1
ij +=1
fout.write(' 0.00000000000000000000E+00' + ' 0 0 0 0\n')
h_aa = h_aa.reshape(nmo,nmo)
h_bb = h_bb.reshape(nmo,nmo)
output_format = float_format + ' %5d %5d 0 0\n'
for i in range(nmo):
for j in range(nmo):
fout.write(output_format % (h_aa[i,j], a_inds[i], a_inds[j]))
fout.write(' 0.00000000000000000000E+00' + ' 0 0 0 0\n')
for i in range(nmo):
for j in range(nmo):
fout.write(output_format % (h_bb[i,j], b_inds[i], b_inds[j]))
fout.write(' 0.00000000000000000000E+00' + ' 0 0 0 0\n')
output_format = float_format + ' 0 0 0 0\n'
fout.write(output_format % nuc)
