LRC hybrid functional energy difference between PySCFAD and PySCF

[jstoppelman]

Hello, I was trying to run a calculation with a long-range corrected hybrid density functional in PySCFAD. I was able to run a single-point energy SCF calculation with this functional, but I am getting a significant difference in the output compared to a standalone PySCF calculation. I have displayed the code below for a fairly simple example case:

 from pyscfad import gto, dft

mol = gto.Mole(atom='H 0 0 0; H 0 0 1.4', basis='ccpvdz')
mol.build()
mf_dft = dft.RKS(mol)
mf_dft.xc = 'HYB_GGA_XC_LRC_WPBE'
mf_dft.kernel()
print(mf_dft.e_tot)

from pyscf import gto, dft
mol = gto.Mole(atom='H 0 0 0; H 0 0 1.4', basis='ccpvdz')
mol.build()
mf_dft = dft.RKS(mol)
mf_dft.xc = 'HYB_GGA_XC_LRC_WPBE'
mf_dft.kernel()
print(mf_dft.e_tot)

I also tried defining a trial LRC functional using a string for debugging purposes and see a similar issue.

from pyscfad import gto, dft

mol = gto.Mole(atom='H 0 0 0; H 0 0 1.4', basis='ccpvdz')
mol.build()
mf_dft = dft.RKS(mol)
mf_dft.xc = 'PBE + 0.1 * SR_HF(0.26) + 0.9 * LR_HF(0.26)'
mf_dft.kernel()
print(mf_dft.e_tot)

from pyscf import gto, dft

mol = gto.Mole(atom='H 0 0 0; H 0 0 1.4', basis='ccpvdz')
mol.build()
mf_dft = dft.RKS(mol)
mf_dft.xc = 'PBE + 0.1 * SR_HF(0.26) + 0.9 * LR_HF(0.26)'
mf_dft.kernel()
print(mf_dft.e_tot)
``
Are LRC functionals designed to be available with PySCFAD? Thanks in advance for any information you can provide.

[jstoppelman]

I see what the problem is now. PySCFAD has the code needed for LRC functionals in pyscfad/dft/rks.py, but there's a missing section in pyscfad/scf/hf.py. I modified my local version of pyscfad/scf/hf.py with the needed code, although you may have a cleaner fix. Using this code, I obtain the same energy with a LRC functional in PySCFAD as I obtain using PySCF for the above test. I also get the same gradient using PySCFAD as I do with PySCF without any further modifications. Very cool code! Basically all I added to pyscfad/scf/hf.py was:

with mol.with_range_coulomb(omega):
        eri_lr = mol.intor('int2e', aosym='s1')

which is used in pyscf.scf.hf.get_jk for computing the long-range ERI. The full modified pyscfad.scf.hf.SCF.get_jk method is shown below if it is useful.

    def get_jk(self, mol=None, dm=None, hermi=1, with_j=True, with_k=True,
               omega=None):
        if mol is None:
            mol = self.mol
        if dm is None:
            dm = self.make_rdm1()
        if self._eri is None and omega is None:
            if config.moleintor_opt:
                self._eri = self.mol.intor('int2e', aosym='s4')
            else:
                self._eri = self.mol.intor('int2e', aosym='s1')
        elif omega:
            #Use the mol.with_range_coulomb code to get the "screened" ERI
            if config.moleintor_opt:
                with mol.with_range_coulomb(omega):
                    eri_lr = mol.intor('int2e', aosym='s4')
            else:
                with mol.with_range_coulomb(omega):
                    eri_lr = mol.intor('int2e', aosym='s1')
            #The function call with omega happens after the code has already been called once
            #self._eri currently contains the "unscreened" ERI, so we want to set that back
            #once we have computed the LR vj and vk
            _eri = self._eri
            self._eri = eri_lr
        vj, vk = dot_eri_dm(self._eri, dm, hermi, with_j, with_k)
        if omega:
            #Reset self._eri to the "unscreened" ERI
            self._eri = _eri
        return vj, vk

[fishjojo]

Hi @jstoppelman, thanks for looking into this problem. I have adopted your fix in #24. The dft module is not well tested and only provides basic functionalities. It may be redesigned in the future based on feedback from the community about what kinds of tasks AD would be useful for.

[jstoppelman]

Ok thanks!