Unit tests failing again, most likely due to new PySCF version

[MatthewRHermes]

Several occri and trexio unit tests are now failing since (I'm guessing) the new PySCF version dropped.

[MatthewRHermes]

PR #192 is my debug attempt

[MatthewRHermes]

There is an actual discrepancy in the behavior of the PWKRCCSD amplitude equations between PySCF 2.12.0 and 2.13.0. Even if I strictly tighten convergence criteria, revert lindep behavior, and load all wave function data other than the amplitudes themselves from checkpoint files, the discrepancy driving the failure of the last unfixed unit test persists. @kylebystrom or @sunqm please help!

[kylebystrom]

I will try to take a look at this today. Thanks for bringing this to my attention!

[MatthewRHermes]

OK so I can make the test pass on my laptop if I absolutely clobber the plane wave basis, setting the cutoffs in lines 57 and 87 to 200, without changing anything else, at the price of it taking 3 times as long to complete. Is that reasonable? I don't have any intuition about the numerical precision of these methods and these thresholds are all over the place.

[kylebystrom]

That's great that increasing the PW cutoff fixes the issue, though it's still very strange that the PySCF updates break the test. The point of the second part of the test (after setting ecut_eri on line 57) is to show that precision is mostly preserved in a small PW basis, which drastically speeds up the calculations. So I will look into it a bit more before we settle on that solution.

The commit that breaks the unit test is pyscf/pyscf@8ae29f6. I suspect it has to do with the change in how the Ewald correction is handled, but I haven't pinpointed exactly where that comes in or why it only affects CCSD.

[kylebystrom]

As an additional debugging note, the correlation energy for GTO CCSD is identical in 2.12 and 2.13. It's the PWKRCCSD result that changes. The SCF energies are the same.

Edit: The SCF energies the same for GTO and PW because the PW is built from the GTO calculation in this case. Between PySCF 2.12 and 2.13, the total SCF energies differ by about 2e-6 Ha. The occupied orbital energies differ by ~3e-6 to 7e-5 Ha for the PW SCF calculation between 2.12 and 2.13, which might be where the discrepancy arises.

[MatthewRHermes]

Yeah I went back and checkpointed the mo_energy from v2.12.0 as well as the orbitals, and the orbital energies in the GTO calculation are different even from the values you get by recomputing them in v2.13.0 from the v2.12.0 mo_coeff array. If you slot in the v2.12.0 mo_energy, the test passes.

[MatthewRHermes]

As you suspected, this and only this change to PySCF core is responsible for the unittest failing:

diff --git a/pyscf/pbc/df/fft_jk.py b/pyscf/pbc/df/fft_jk.py
index 3e7601332..72f591322 100644
--- a/pyscf/pbc/df/fft_jk.py
+++ b/pyscf/pbc/df/fft_jk.py
@@ -26,8 +26,8 @@ from pyscf import lib
 from pyscf.lib import logger
 from pyscf.pbc import tools
 from pyscf.pbc.df.df_jk import _format_dms, _format_kpts_band, _format_jks
+from pyscf.pbc.df.df_jk import _ewald_exxdiv_for_G0
 from pyscf.pbc.lib.kpts_helper import is_zero
-from pyscf.pbc.tools.k2gamma import kpts_to_kmesh
 
 
 def get_j_kpts(mydf, dm_kpts, hermi=1, kpts=np.zeros((1,3)), kpts_band=None):
@@ -267,12 +267,13 @@ def get_k_kpts(mydf, dm_kpts, hermi=1, kpts=np.zeros((1,3)), kpts_band=None,
         for k1, ao1T in enumerate(ao1_kpts):
             kpt1 = kpts_band[k1]
 
-            # In PySCF v1.5 - v2.12, the G=0 term is evaluated analytically
-            # using _ewald_exxdiv_for_G0. The G=0 component obtained here may
-            # differ from _ewald_exxdiv_for_G0 due to discretization errors in
-            # the FFT-based density (especially when the mesh is not
-            # sufficiently dense), which can lead to small discrepancies.
-            coulG = tools.get_coulG(cell, kpt2-kpt1, exxdiv, mydf, mesh)
+            # If we have an ewald exxdiv, we add the G=0 correction near the
+            # end of the function to bypass any discretization errors
+            # that arise from the FFT.
+            if exxdiv == 'ewald':
+                coulG = tools.get_coulG(cell, kpt2-kpt1, False, mydf, mesh)
+            else:
+                coulG = tools.get_coulG(cell, kpt2-kpt1, exxdiv, mydf, mesh)
             if is_zero(kpt1-kpt2):
                 expmikr = np.array(1.)
             else:
@@ -294,7 +295,15 @@ def get_k_kpts(mydf, dm_kpts, hermi=1, kpts=np.zeros((1,3)), kpts_band=None,
 
             for i in range(nset):
                 vk_kpts[i,k1] += weight * lib.dot(vR_dm[i], ao1T.T)
-        t1 = log.timer_debug1('get_k_kpts: make_kpt (%d,*)'%k2, *t1)
+        t1 = logger.timer_debug1(mydf, 'get_k_kpts: make_kpt (%d,*)'%k2, *t1)
+
+    # Function _ewald_exxdiv_for_G0 to add back in the G=0 component to vk_kpts
+    # Note in the _ewald_exxdiv_for_G0 implementation, the G=0 treatments are
+    # different for 1D/2D and 3D systems.  The special treatments for 1D and 2D
+    # can only be used with AFTDF/GDF/MDF method.  In the FFTDF method, 1D, 2D
+    # and 3D should use the ewald probe charge correction.
+    if exxdiv == 'ewald' and cell.dimension != 0:
+        _ewald_exxdiv_for_G0(cell, kpts, dms, vk_kpts, kpts_band=kpts_band)
 
     return _format_jks(vk_kpts, dm_kpts, input_band, kpts)

[MatthewRHermes]

OK, so in v2.12.0 we do this:

def _ewald_exxdiv_for_G0(cell, kpts, dms, vk, kpts_band=None):
    s = cell.pbc_intor('int1e_ovlp', hermi=1, kpts=kpts)
    madelung = tools.pbc.madelung(cell, kpts)
    if kpts is None:
        for i,dm in enumerate(dms):
            vk[i] += madelung * reduce(numpy.dot, (s, dm, s))
    elif numpy.shape(kpts) == (3,):
        if kpts_band is None or is_zero(kpts_band-kpts):
            for i,dm in enumerate(dms):
                vk[i] += madelung * reduce(numpy.dot, (s, dm, s))

    elif kpts_band is None or numpy.array_equal(kpts, kpts_band):
        for k in range(len(kpts)):
            for i,dm in enumerate(dms):
                vk[i,k] += madelung * reduce(numpy.dot, (s[k], dm[k], s[k]))
    else:
        for k, kpt in enumerate(kpts):
            for kp in member(kpt, kpts_band.reshape(-1,3)):
                for i,dm in enumerate(dms):
                    vk[i,kp] += madelung * reduce(numpy.dot, (s[k], dm[k], s[k]))

and in v2.13.0 we instead do this:

    # For full-range Coulomb and long-range Coulomb,
    # the divergent part of periodic summation of (ii|ii) integrals in
    # Coulomb integrals were cancelled out by electron-nucleus
    # interaction. The periodic part of (ii|ii) in exchange cannot be
    # cancelled out by Coulomb integrals. Its leading term is calculated
    # using Ewald probe charge (the function madelung below)
    if cell.dimension > 0 and exxdiv == 'ewald' and len(G0_idx) > 0:
        if omega is None: # Affects DFT-RSH
            coulG[G0_idx] += Nk*cell.vol*madelung(cell, kpts)
        else: # for RangeSeparatedJKBuilder
            coulG[G0_idx] += Nk*cell.vol*madelung(cell, kpts, omega=0)

and PWKRCCSD has a block of code that does this:

pyscf-forge/pyscf/pbc/pwscf/kccsd_rhf.py

Lines 104 to 115 in 0286eaa

	if False:#cc.keep_exxdiv:
	moe_noewald = self.mo_energy
	else:
	# remove ewald correction
	moe_noewald = np.zeros_like(self.mo_energy)
	for k in range(nkpts):
	moe = self.mo_energy[k].copy()
	moe[mo_occ[k]>THR_OCC] += mf._madelung
	moe_noewald[k] = moe
	self.fock = np.asarray(
	[np.diag(moe.astype(np.complex128)) for moe in moe_noewald]
	)

in which the comment "remove ewald correction" is suggesting to me that it's removing the old Ewald correction, and needs to be updated for the new, v2.13.0, Ewald correction. Does that make any sense?

[MatthewRHermes]

If there were a straightforward way to just recompute mo_energy for the PW SCF classes I think this would be trivial; just set exxdiv to None temporarily. But canonicalize or get_bands etc. don't seem to be implemented for PW SCF.

[kylebystrom]

Yes this all looks correct, and thanks of much for identifying the specific chain of issues! I think I know how to fix this.

First off, I think that last chunk of code I wrote is a bad idea because it assumes exxdiv=ewald when it might be None, so I need to rewrite it anyway. Second, gtomf2pwmf copies the mo energies from the GTO calculation rather than re-computing them. This is important because the mo energies can be a bit different in the PW basis than the GTO basis due to finite basis error.

I think what should happen is that when gtomf2pwmf is called, the GTO mo energies are re-computed with exxdiv=None. Then the appropriate exxdiv in the PW mode is added back in. In the CC module, the _scf.exxdiv setting is then checked before determining if the Madelung term needs to be removed.

To your last comment, there is a get_mo_energy implemented for PW mode that I tried using yesterday, but I was having trouble with it. I will play around with that too to make sure it works properly.

I will take care of this today.

[kylebystrom]

@MatthewRHermes I made PR #193 to address the PW test issues. Thanks again for identifying the problem! Please let me know if you want any additional changes.

Based on this discussion and exploration of the code, there are definitely some additional features and utilities that PW mode is lacking that could be useful later, but in the meantime these fixes should address the immediate issues.

[sunqm]

A short explanation for the small changes in orbital energies or total energies you might observe in PySCF 2.13.

Two changes introduced in PySCF 2.13 can affect these values.

First, linear dependency removal for the orbital basis is enabled by default. This uses a threshold of 1e-6 for overlap eigenvalues to remove linearly dependent basis functions.

Second, the exxdiv=ewald correction at G=0 is changed. Previously, this correction was handled analytically, as in the additional code block mentioned by @MatthewRHermes. In PySCF 2.13, it is computed together with the FFTDF integrals. When the energy cutoff is insufficient, the numerical result can differ slightly from the analytical value. Increasing the energy cutoff (or cell.precision) can reduce this difference.

[MatthewRHermes]

Hey @sunqm, thanks. I think I've accounted for all of the changes in #192. Could you take a look at it?