@@ -71,7 +71,7 @@ def get_molint_from_C(cell, C_ks, kpts, mo_slices=None, exxdiv=None,
         eri_mesh = cell.cutoff_to_mesh(ecut_eri)
         if (np.array(mesh) == eri_mesh).all():
             eri_mesh = None
-        if (eri_mesh > mesh).any():
+        elif (eri_mesh > mesh).any():
             logger.warn(
                 cell,
                 "eri_mesh larger than mesh, so eri_mesh is ignored."
 
@@ -58,6 +58,10 @@ class PWKRCCSD(cc.kccsd_rhf.RCCSD):
     Other CCSD results can be accessed from the underlaying
     `mcc` (Periodic RCCSD, which contains e_tot, t-amplitudes, etc.)
     after `kernel` is called.
+
+    `ecut_eri` can be set to use a sparser grid for ERI evaluation,
+    reducing computational cost. NOTE that setting `ecut_eri` too
+    small can cause numerical issues.
     """
 
     _keys = {'ecut_eri'}
@@ -101,15 +105,28 @@ def __init__(self, cc):
         self.e_hf = mf.e_tot
         self.mo_energy = np.asarray(mo_energy)
 
-        if False:#cc.keep_exxdiv:
-            moe_noewald = self.mo_energy
-        else:
-            # remove ewald correction
+        if cc._scf.exxdiv == "ewald" and not cc.keep_exxdiv:
+            # remove ewald correction form mo_energy to make self.fock
             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
+        else:
+            moe_noewald = self.mo_energy
+            if not cc.keep_exxdiv:
+                # To match the GTO implementation, we add Madelung correction
+                # to the mo_energy even if exxdiv=None.
+                # TODO is this actually the best practice? One might expect
+                # exxdiv to not be used anywhere if it is explicitly
+                # set to None by the user. For now, we match the
+                # behavior of the equivalent GTO CCSD calculation.
+                # This choice affects initial amplitude guess as well
+                # as the perturbative triples energy.
+                from pyscf.pbc.cc.ccsd import _adjust_occ
+                mf._set_madelung()
+                self.mo_energy = [_adjust_occ(mo_e, cc.nocc, -mf._madelung)
+                                  for k, mo_e in enumerate(self.mo_energy)]
         self.fock = np.asarray(
             [np.diag(moe.astype(np.complex128)) for moe in moe_noewald]
         )
 
@@ -44,7 +44,7 @@
 from pyscf.pbc.pwscf import chkfile
 from pyscf.pbc import gto, scf, tools
 from pyscf.pbc.pwscf import pw_helper
-from pyscf.pbc.pwscf.pw_helper import get_kcomp, set_kcomp
+from pyscf.pbc.pwscf.pw_helper import get_kcomp, set_kcomp, ewald_correction
 from pyscf.pbc.pwscf import pseudo as pw_pseudo
 from pyscf.pbc.pwscf import jk as pw_jk
 from pyscf.lib import logger
@@ -1159,22 +1159,6 @@ def apply_Fock_kpt(mf, C_k, kpt, mocc_ks, mesh, Gv, vj_R, exxdiv,
         return Cbar_k
 
 
-def ewald_correction(moe_ks, mocc_ks, madelung):
-    if isinstance(moe_ks[0][0], float): # RHF
-        nkpts = len(moe_ks)
-        moe_ks_new = [None] * nkpts
-        for k in range(nkpts):
-            moe_ks_new[k] = moe_ks[k].copy()
-            moe_ks_new[k][:] -= 0.5 * mocc_ks[k] * madelung
-    else:                               # UHF
-        ncomp = len(moe_ks)
-        moe_ks_new = [None] * ncomp
-        for comp in range(ncomp):
-            moe_ks_new[comp] = ewald_correction(moe_ks[comp], mocc_ks[comp],
-                                                ncomp * madelung)
-    return moe_ks_new
-
-
 def get_mo_energy(mf, C_ks, mocc_ks, mesh=None, Gv=None, exxdiv=None,
                   vj_R=None, comp=None, ret_mocc=True, full_ham=False):
     """
 
@@ -426,6 +426,10 @@ def PWKRMP2_from_gtomf(mf, chkfile=None):
 class PWKRMP2:
     """
     Restriced MP2 perturbation theory in a plane-wave basis.
+
+    `ecut_eri` can be set to use a sparser grid for ERI evaluation,
+    reducing computational cost. NOTE that setting `ecut_eri` too
+    small can cause numerical issues.
     """
     def __init__(self, mf, nvir=None, frozen=None):
         self.cell = self.mol = mf.cell
 
@@ -387,6 +387,10 @@ def kernel_dx_(cell, kpts, chkfile_name, summary, nvir=None, nvir_lst=None,
 class PWKUMP2(kmp2.PWKRMP2):
     """
     Spin-unrestriced MP2 in a plane-wave basis.
+
+    `ecut_eri` can be set to use a sparser grid for ERI evaluation,
+    reducing computational cost. NOTE that setting `ecut_eri` too
+    small can cause numerical issues.
     """
     def __init__(self, mf, nvir=None):
         kmp2.PWKRMP2.__init__(self, mf, nvir=nvir)
 
@@ -105,6 +105,22 @@ def get_basis_data(cell, kpts, ecut_wf, ecut_rho=None, wf_mesh=None,
     return np.array(wf_mesh), np.array(xc_mesh), wf2xc, basis_data
 
 
+def ewald_correction(moe_ks, mocc_ks, madelung):
+    if isinstance(moe_ks[0][0], float): # RHF
+        nkpts = len(moe_ks)
+        moe_ks_new = [None] * nkpts
+        for k in range(nkpts):
+            moe_ks_new[k] = moe_ks[k].copy()
+            moe_ks_new[k][:] -= 0.5 * mocc_ks[k] * madelung
+    else:                               # UHF
+        ncomp = len(moe_ks)
+        moe_ks_new = [None] * ncomp
+        for comp in range(ncomp):
+            moe_ks_new[comp] = ewald_correction(moe_ks[comp], mocc_ks[comp],
+                                                ncomp * madelung)
+    return moe_ks_new
+
+
 def wf_fft(C_k_R, mesh, basis=None):
     """
     Fourier transform a wave function from real to reciprocal space.
@@ -468,9 +484,20 @@ def gto2cpw(cell, basis, kpts, amin=None, amax=None, ke_or_mesh=None, out=None):
     return out
 
 
+def _reduce_fock_to_moe(fockao, mo_coeff):
+    return [np.einsum("ai,ai->i", mo.conj(), fockao[k].dot(mo)).real
+            for k, mo in enumerate(mo_coeff)]
+
+
 def gtomf2pwmf(mf, chkfile=None):
     """
+    Create a PW HF calculation from a GTO HF calculation.
+    If mf.exxdiv is "ewald", the returned PW HF object
+    has exxdiv="ewald", otherwise exxdiv=None.
+
     Args:
+        mf:
+            GTO HF calculation
         chkfile (str):
             A hdf5 file to store chk variables (mo_energy, mo_occ, etc.).
             If not provided, a temporary file is generated.
@@ -484,44 +511,79 @@ def gtomf2pwmf(mf, chkfile=None):
     nkpts = len(kpts)
 # transform GTO MO coeff to PW MO coeff
     Cgto_ks = mf.mo_coeff
+
+    exxdiv = "ewald" if mf.exxdiv == "ewald" else None
+    if mf.exxdiv is not None:
+        dm = mf.make_rdm1(mf.mo_coeff, mf.mo_occ)
+        with lib.temporary_env(mf, exxdiv=None):
+            # Calculate veff without the exxdiv
+            vhf = mf.get_veff(cell, dm)
+        fockao = mf.get_hcore() + vhf
+
     if isinstance(mf, scf.khf.KRHF):
-        pwmf = pwscf.KRHF(cell, kpts)
+        pwmf = pwscf.KRHF(cell, kpts, exxdiv=exxdiv)
         nmo_ks = [Cgto_ks[k].shape[1] for k in range(nkpts)]
         pwmf.mo_coeff = C_ks = get_C_ks_G(cell, kpts, Cgto_ks, nmo_ks)
-        pwmf.mo_energy = moe_ks = mf.mo_energy
+        if mf.exxdiv is None:
+            pwmf.mo_energy = moe_ks = mf.mo_energy
+        else:
+            pwmf.mo_energy = moe_ks = _reduce_fock_to_moe(fockao, mf.mo_coeff)
         pwmf.mo_occ = mocc_ks = mf.mo_occ
         pwmf.e_tot = mf.e_tot
     elif isinstance(mf, scf.kuhf.KUHF):
-        pwmf = pwscf.KUHF(cell, kpts)
+        pwmf = pwscf.KUHF(cell, kpts, exxdiv=exxdiv)
         C_ks = [None] * 2
         for s in [0,1]:
             nmo_ks = [Cgto_ks[s][k].shape[1] for k in range(nkpts)]
             C_ks[s] = get_C_ks_G(cell, kpts, Cgto_ks[s], nmo_ks)
         pwmf.mo_coeff = C_ks
-        pwmf.mo_energy = moe_ks = mf.mo_energy
+        if mf.exxdiv is None:
+            pwmf.mo_energy = moe_ks = mf.mo_energy
+        else:
+            pwmf.mo_energy = moe_ks = [
+                _reduce_fock_to_moe(fockao[0], mf.mo_coeff[0]),
+                _reduce_fock_to_moe(fockao[1], mf.mo_coeff[1]),
+            ]
         pwmf.mo_occ = mocc_ks = mf.mo_occ
         pwmf.e_tot = mf.e_tot
     elif isinstance(mf, scf.hf.RHF):
-        pwmf = pwscf.KRHF(cell, kpts)
+        pwmf = pwscf.KRHF(cell, kpts, exxdiv=exxdiv)
         nmo_ks = [Cgto_ks.shape[1]]
         C_ks = get_C_ks_G(cell, kpts, [Cgto_ks], nmo_ks)
         pwmf.mo_coeff = C_ks
-        pwmf.mo_energy = moe_ks = [mf.mo_energy]
+        if mf.exxdiv is None:
+            pwmf.mo_energy = moe_ks = [mf.mo_energy]
+        else:
+            pwmf.mo_energy = moe_ks = _reduce_fock_to_moe(
+                [fockao], [mf.mo_coeff]
+            )
         pwmf.mo_occ = mocc_ks = [mf.mo_occ]
         pwmf.e_tot = mf.e_tot
     elif isinstance(mf, scf.uhf.UHF):
-        pwmf = pwscf.KUHF(cell, kpts)
+        pwmf = pwscf.KUHF(cell, kpts, exxdiv=exxdiv)
         C_ks = [None] * 2
         for s in [0,1]:
             nmo_ks = [Cgto_ks[s].shape[1]]
             C_ks[s] = get_C_ks_G(cell, kpts, [Cgto_ks[s]], nmo_ks)
         pwmf.mo_coeff = C_ks
-        pwmf.mo_energy = moe_ks = [[mf.mo_energy[s]] for s in [0,1]]
+        if mf.exxdiv is None:
+            pwmf.mo_energy = moe_ks = [[mf.mo_energy[s]] for s in [0,1]]
+        else:
+            pwmf.mo_energy = moe_ks = [
+                _reduce_fock_to_moe([fockao[0]], [mf.mo_coeff[0]]),
+                _reduce_fock_to_moe([fockao[1]], [mf.mo_coeff[1]]),
+            ]
         pwmf.mo_occ = mocc_ks = [[mf.mo_occ[s]] for s in [0,1]]
         pwmf.e_tot = mf.e_tot
     else:
         raise TypeError
-# update chkfile
+
+    if exxdiv == "ewald":
+        pwmf.mo_energy = moe_ks = ewald_correction(
+            pwmf.mo_energy, pwmf.mo_occ, pwmf.madelung
+        )
+
+    # update chkfile
     if chkfile is None:
         swapfile = tempfile.NamedTemporaryFile(dir=lib.param.TMPDIR)
         chkfile = swapfile.name
 
@@ -20,7 +20,7 @@
 
 class KnownValues(unittest.TestCase):
     def _run_test(self, atom, a, basis, pseudo, ke_cutoff, kmesh, exxdiv,
-                  test_scf=True):
+                  test_scf=True, keep_exxdiv=False):
         # cell
         cell = gto.Cell(
             atom=atom,
@@ -37,33 +37,39 @@ def _run_test(self, atom, a, basis, pseudo, ke_cutoff, kmesh, exxdiv,
         gmf = scf.KRHF(cell, kpts)
         gmf.exxdiv = exxdiv
         gmf.kernel()
+        assert gmf.converged
         gcc = cc.KCCSD(gmf)
+        gcc.keep_exxdiv = keep_exxdiv
         gcc.kernel()
-        if test_scf:
-            etrip = gcc.ccsd_t()
-            ips_ref = gcc.ipccsd()
+        assert gcc.converged
+        etrip = gcc.ccsd_t()
+        ips_ref = gcc.ipccsd()
 
         # PW
         pmf = pw_helper.gtomf2pwmf(gmf)
         pmf.build()
         pmf.update_pp(pmf.mo_coeff)
         pmf.update_k(pmf.mo_coeff, pmf.mo_occ)
+        assert pmf.converged
         pcc = pwscf.PWKRCCSD(pmf)
+        pcc.keep_exxdiv = keep_exxdiv
         pcc.kernel()
-        if test_scf:
-            # Check CCSD(T)
-            pcc.ccsd_t()
-            assert(abs(gcc.e_corr - pcc.e_corr) < 1.e-6)
-        pcc.ecut_eri = 15
+        assert pcc.converged
+        etrip_test = pcc.ccsd_t()
+        # Check CCSD(T)
+        assert_allclose(gcc.e_corr, pcc.e_corr, atol=1.e-6, rtol=0)
+        assert_allclose(etrip, etrip_test, atol=1.e-5, rtol=0)
+
+        pcc.ecut_eri = 25
         pcc.kernel()
-        if test_scf:
-            # Check EOM-CCSD
-            etrip_test = pcc.ccsd_t()
-            ips_test = pcc.ipccsd()
-            assert(abs(gcc.e_corr - pcc.e_corr) < 1.e-4)
-            assert(abs(etrip - etrip_test) < 1.e-5)
-            assert(np.max(np.abs(ips_ref[0] - ips_test[0])) < 3.e-5)
-            assert(np.max(np.abs(ips_ref[1] - ips_test[1])) < 1.e-3)
+        assert pcc.converged
+        # Check EOM-CCSD
+        etrip_test = pcc.ccsd_t()
+        ips_test = pcc.ipccsd()
+        assert_allclose(gcc.e_corr, pcc.e_corr, atol=1.e-4, rtol=0)
+        assert_allclose(etrip, etrip_test, atol=1.e-5, rtol=0)
+        assert_allclose(ips_ref[0], ips_test[0], atol=3.e-5, rtol=0)
+        assert_allclose(ips_ref[1], ips_test[1], atol=1.e-3, rtol=0)
 
         if test_scf:
             pwmf = pwscf.KRHF(cell, kpts, ecut_wf=20)
@@ -75,8 +81,10 @@ def _run_test(self, atom, a, basis, pseudo, ke_cutoff, kmesh, exxdiv,
             # pwmf.damp_type = "simple"
             pwmf.damp_factor = 0.0
             pwmf.kernel()
+            assert pwmf.converged
             pwcc = pwscf.PWKRCCSD(pwmf)
             pwcc.kernel()
+            assert pwcc.converged
             # Just to make sure the code stays consistent
             assert_allclose(pwcc.e_corr, -0.03234330656841895, atol=1.e-4, rtol=0)
             pwcc.ecut_eri = 15
@@ -87,16 +95,23 @@ def test_krccsd(self):
         ke_cutoff = 50
         basis = "gth-szv"
         pseudo = "gth-pade"
-        exxdiv = "ewald"
         atom = "Li 0 0 0; Li 1.75 1.75 1.75"
         a = np.eye(3) * 3.5
 
         # same occ per kpt
         kmesh = [2,1,1]
-        self._run_test(atom, a, basis, pseudo, ke_cutoff, kmesh, exxdiv)
+        self._run_test(atom, a, basis, pseudo, ke_cutoff, kmesh, "ewald")
+        kmesh = [3,1,1]
+        self._run_test(atom, a, basis, pseudo, ke_cutoff, kmesh, "ewald",
+                       test_scf=False, keep_exxdiv=True)
+        self._run_test(atom, a, basis, pseudo, ke_cutoff, kmesh, None,
+                       test_scf=False, keep_exxdiv=True)
+        self._run_test(atom, a, basis, pseudo, ke_cutoff, kmesh, None,
+                       test_scf=False)
+
         # diff occ per kpt (i.e., needs padding)
         kmesh = [2,2,1]
-        self._run_test(atom, a, basis, pseudo, ke_cutoff, kmesh, exxdiv,
+        self._run_test(atom, a, basis, pseudo, ke_cutoff, kmesh, "ewald",
                        test_scf=False)
 
 
 
@@ -19,8 +19,7 @@
 
 
 class KnownValues(unittest.TestCase):
-    def test_krmp2(self):
-        kmesh = [2,1,1]
+    def _check_krmp2(self, kmesh):
         ke_cutoff = 100
         pseudo = "gth-pade"
         exxdiv = "ewald"
@@ -40,18 +39,29 @@ def test_krmp2(self):
         kpts = cell.make_kpts(kmesh)
 
         # GTO
-        gmf = scf.KRHF(cell, kpts)
-        gmf.exxdiv = exxdiv
-        gmf.kernel()
-        gmp = mp.KMP2(gmf)
-        gmp.kernel()
+        if (np.array(kmesh) == 1).all():
+            gmf = scf.RHF(cell)
+            gmf.exxdiv = exxdiv
+            gmf.kernel()
+            gmp = mp.MP2(gmf)
+            gmp.kernel()
+        else:
+            gmf = scf.KRHF(cell, kpts)
+            gmf.exxdiv = exxdiv
+            gmf.kernel()
+            gmp = mp.KMP2(gmf)
+            gmp.kernel()
 
         # PW
         pmf = pw_helper.gtomf2pwmf(gmf)
         pmp = kmp2.PWKRMP2(pmf)
         pmp.kernel()
         assert(abs(gmp.e_corr - pmp.e_corr) < 1.e-6)
 
+    def test_kump2(self):
+        self._check_krmp2([1, 1, 1])
+        self._check_krmp2([2, 1, 1])
+
 
 if __name__ == "__main__":
     unittest.main()