Chkfile initial guess problem in ROHF and ROKS (issue pyscf#1986) (py…

…scf#2010)

* Fix chkfile initial guess problem for ROHF and ROKS (issue pyscf#1986)

* Use UHF natural orbital for RHF chkfile initial guess

* Enable some xc_deriv tests

* Enable xcfun compilation

.github/workflows/ci_linux/build_pyscf.sh

@@ -5,7 +5,7 @@ set -e
 cd ./pyscf/lib
 curl -L "https://github.com/pyscf/pyscf-build-deps/blob/master/pyscf-2.4a-deps.tar.gz?raw=true" | tar xzf -
 mkdir build; cd build
-cmake -DBUILD_LIBXC=OFF -DBUILD_XCFUN=OFF -DBUILD_LIBCINT=OFF -DXCFUN_MAX_ORDER=4 ..
+cmake -DBUILD_LIBXC=OFF -DBUILD_XCFUN=ON -DBUILD_LIBCINT=OFF -DXCFUN_MAX_ORDER=4 ..
 make -j4
 cd ..
 rm -Rf build

pyscf/dft/test/test_he.py

@@ -14,6 +14,7 @@
 # limitations under the License.
 
 import unittest
+import tempfile
 import numpy
 from pyscf import gto
 from pyscf import lib
@@ -191,6 +192,21 @@ def test_convert(self):
         self.assertTrue(isinstance(udks.to_dhf(), scf.dhf.DHF))
         self.assertTrue(isinstance(udks.to_dks('pbe'), dft.dks.DKS))
 
+    # issue 1986
+    def test_init_guess_chkfile(self):
+        with tempfile.NamedTemporaryFile() as tmpf:
+            mol = gto.M(atom='He 0 0 0', basis='631g', charge=1, spin=1)
+            mf = dft.RKS(mol)
+            mf.chkfile = tmpf.name
+            e1 = mf.kernel()
+            mf = dft.RKS(mol)
+            mf.init_guess = 'chkfile'
+            mf.chkfile = tmpf.name
+            mf.max_cycle = 1
+            e2 = mf.kernel()
+            self.assertAlmostEqual(e1, e2, 9)
+
+
 if __name__ == "__main__":
     print("Full Tests for He")
     unittest.main()

pyscf/dft/test/test_xc_deriv.py

@@ -415,7 +415,7 @@ def test_libxc_mgga_deriv3(self):
         self.assertAlmostEqual(xc1.sum(), 19.977512805950784, 7)
         self.assertAlmostEqual(abs(ref[1] - xc1).max(), 0, 9)
 
-    @unittest.skip(dft.libxc.max_deriv_order('pbe,') <= 3)
+    @unittest.skipIf(dft.libxc.max_deriv_order('pbe,') <= 3, 'libxc order')
     def test_libxc_gga_deriv4(self):
         rho1 = rho[:,:4].copy()
         xc1 = dft.libxc.eval_xc_eff('PBE', rho1, deriv=4)
@@ -425,7 +425,7 @@ def test_libxc_gga_deriv4(self):
         xc1 = dft.libxc.eval_xc_eff('PBE', rho1, deriv=4)
         self.assertAlmostEqual(xc1.sum(), -869.6617638095072, 4)
 
-    @unittest.skip(not hasattr(dft, 'xcfun'))
+    @unittest.skipIf(not hasattr(dft, 'xcfun'), 'xcfun order')
     def test_xcfun_lda_deriv3(self):
         rho1 = rho[:,0].copy()
         ref = eval_xc_eff('LDA,', rho1, 3, dft.xcfun)
@@ -457,7 +457,7 @@ def test_xcfun_lda_deriv3(self):
         self.assertAlmostEqual(xc1.sum(), -2.7477984980958627, 9)
         self.assertAlmostEqual(abs(ref[1] - xc1).max(), 0, 9)
 
-    @unittest.skip(not hasattr(dft, 'xcfun'))
+    @unittest.skipIf(not hasattr(dft, 'xcfun'), 'xcfun order')
     def test_xcfun_gga_deriv3(self):
         rho1 = rho[:,:4].copy()
         ref = eval_xc_eff('PBE', rho1, 3, dft.xcfun)
@@ -489,7 +489,7 @@ def test_xcfun_gga_deriv3(self):
         self.assertAlmostEqual(xc1.sum(), -3.0715856471099032, 9)
         self.assertAlmostEqual(abs(ref[1] - xc1).max(), 0, 9)
 
-    @unittest.skip(not hasattr(dft, 'xcfun'))
+    @unittest.skipIf(not hasattr(dft, 'xcfun'), 'xcfun order')
     def test_xcfun_mgga_deriv3(self):
         rho1 = rho
         ref = eval_xc_eff('M06', rho1, 3, dft.xcfun)
@@ -521,7 +521,7 @@ def test_xcfun_mgga_deriv3(self):
         self.assertAlmostEqual(xc1.sum(), 19.977512805950784, 9)
         self.assertAlmostEqual(abs(ref[1] - xc1).max(), 0, 9)
 
-    @unittest.skip(not (hasattr(dft, 'xcfun') and dft.xcfun.MAX_DERIV_ORDER > 3))
+    @unittest.skipIf(not (hasattr(dft, 'xcfun') and dft.xcfun.MAX_DERIV_ORDER > 3), 'xcfun order')
     def test_xcfun_gga_deriv4(self):
         rho1 = rho[:,:4].copy()
         xc1 = dft.xcfun.eval_xc_eff('PBE', rho1, deriv=4)

pyscf/scf/hf.py

@@ -656,12 +656,66 @@ def init_guess_by_chkfile(mol, chkfile_name, project=None):
     '''Read the HF results from checkpoint file, then project it to the
     basis defined by ``mol``
 
+    Kwargs:
+        project : None or bool
+            Whether to project chkfile's orbitals to the new basis.  Note when
+            the geometry of the chkfile and the given molecule are very
+            different, this projection can produce very poor initial guess.
+            In PES scanning, it is recommended to switch off project.
+
+            If project is set to None, the projection is only applied when the
+            basis sets of the chkfile's molecule are different to the basis
+            sets of the given molecule (regardless whether the geometry of
+            the two molecules are different).  Note the basis sets are
+            considered to be different if the two molecules are derived from
+            the same molecule with different ordering of atoms.
+
     Returns:
         Density matrix, 2D ndarray
     '''
-    from pyscf.scf import uhf
-    dm = uhf.init_guess_by_chkfile(mol, chkfile_name, project)
-    return dm[0] + dm[1]
+    from pyscf.scf import addons
+    chk_mol, scf_rec = chkfile.load_scf(chkfile_name)
+    if project is None:
+        project = not gto.same_basis_set(chk_mol, mol)
+
+    # Check whether the two molecules are similar
+    im1 = scipy.linalg.eigvalsh(mol.inertia_moment())
+    im2 = scipy.linalg.eigvalsh(chk_mol.inertia_moment())
+    # im1+1e-7 to avoid 'divide by zero' error
+    if abs((im1-im2)/(im1+1e-7)).max() > 0.01:
+        logger.warn(mol, "Large deviations found between the input "
+                    "molecule and the molecule from chkfile\n"
+                    "Initial guess density matrix may have large error.")
+
+    if project:
+        s = get_ovlp(mol)
+
+    def fproj(mo):
+        if project:
+            mo = addons.project_mo_nr2nr(chk_mol, mo, mol)
+            norm = numpy.einsum('pi,pi->i', mo.conj(), s.dot(mo))
+            mo /= numpy.sqrt(norm)
+        return mo
+
+    mo = scf_rec['mo_coeff']
+    mo_occ = scf_rec['mo_occ']
+    if getattr(mo[0], 'ndim', None) == 1:  # RHF
+        if numpy.iscomplexobj(mo):
+            raise NotImplementedError('TODO: project DHF orbital to UHF orbital')
+        mo_coeff = fproj(mo)
+        dm = make_rdm1(mo_coeff, mo_occ)
+    else:  #UHF
+        if getattr(mo[0][0], 'ndim', None) == 2:  # KUHF
+            logger.warn(mol, 'k-point UHF results are found.  Density matrix '
+                        'at Gamma point is used for the molecular SCF initial guess')
+            mo = mo[0]
+        dma = make_rdm1(fproj(mo[0]), mo_occ[0])
+        dmb = make_rdm1(fproj(mo[1]), mo_occ[1])
+        dm = dma + dmb
+        s = get_ovlp(mol)
+        _, mo_coeff = scipy.linalg.eigh(dm, s, type=2)
+        dm = lib.tag_array(dm, mo_coeff=mo_coeff[:,::-1], mo_occ=mo_occ)
+    return dm
 
 
 def get_init_guess(mol, key='minao'):
@@ -1358,11 +1412,9 @@ def __call__(self, mol_or_geom, **kwargs):
             dm0 = kwargs.pop('dm0')
         elif self.mo_coeff is None:
             dm0 = None
-        elif self.chkfile and h5py.is_hdf5(self.chkfile):
-            dm0 = self.from_chk(self.chkfile)
         else:
             dm0 = None
-            # dm0 form last calculation cannot be used in the current
+            # dm0 form last calculation may not be used in the current
             # calculation if a completely different system is given.
             # Obviously, the systems are very different if the number of
             # basis functions are different.
@@ -1371,6 +1423,8 @@ def __call__(self, mol_or_geom, **kwargs):
             # last calculation.
             if numpy.array_equal(self._last_mol_fp, mol.ao_loc):
                 dm0 = self.make_rdm1()
+            elif self.chkfile and h5py.is_hdf5(self.chkfile):
+                dm0 = self.from_chk(self.chkfile)
         self.mo_coeff = None  # To avoid last mo_coeff being used by SOSCF
         e_tot = self.kernel(dm0=dm0, **kwargs)
         self._last_mol_fp = mol.ao_loc

pyscf/scf/rohf.py

@@ -49,7 +49,28 @@ def init_guess_by_atom(mol):
 
 init_guess_by_huckel = uhf.init_guess_by_huckel
 init_guess_by_mod_huckel = uhf.init_guess_by_mod_huckel
-init_guess_by_chkfile = uhf.init_guess_by_chkfile
+
+def init_guess_by_chkfile(mol, chkfile_name, project=None):
+    '''Read SCF chkfile and make the density matrix for ROHF initial guess.
+
+    Kwargs:
+        project : None or bool
+            Whether to project chkfile's orbitals to the new basis.  Note when
+            the geometry of the chkfile and the given molecule are very
+            different, this projection can produce very poor initial guess.
+            In PES scanning, it is recommended to switch off project.
+
+            If project is set to None, the projection is only applied when the
+            basis sets of the chkfile's molecule are different to the basis
+            sets of the given molecule (regardless whether the geometry of
+            the two molecules are different).  Note the basis sets are
+            considered to be different if the two molecules are derived from
+            the same molecule with different ordering of atoms.
+    '''
+    dm = uhf.init_guess_by_chkfile(mol, chkfile_name, project)
+    mo_coeff = dm.mo_coeff[0]
+    mo_occ = dm.mo_occ[0] + dm.mo_occ[1]
+    return lib.tag_array(dm, mo_coeff=mo_coeff, mo_occ=mo_occ)
 
 def get_fock(mf, h1e=None, s1e=None, vhf=None, dm=None, cycle=-1, diis=None,
              diis_start_cycle=None, level_shift_factor=None, damp_factor=None):