From 35835bcb950b525aa89df6d8eb6e0f0366627549 Mon Sep 17 00:00:00 2001 From: Qiming Sun Date: Thu, 25 Jan 2024 00:49:09 -0800 Subject: [PATCH] Chkfile initial guess problem in ROHF and ROKS (issue #1986) (#2010) * Fix chkfile initial guess problem for ROHF and ROKS (issue #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 | 2 +- pyscf/dft/test/test_he.py | 16 ++++++ pyscf/dft/test/test_xc_deriv.py | 10 ++-- pyscf/scf/hf.py | 66 ++++++++++++++++++++--- pyscf/scf/rohf.py | 23 +++++++- 5 files changed, 104 insertions(+), 13 deletions(-) diff --git a/.github/workflows/ci_linux/build_pyscf.sh b/.github/workflows/ci_linux/build_pyscf.sh index 04e2aac97b..d9bd1b39e7 100755 --- a/.github/workflows/ci_linux/build_pyscf.sh +++ b/.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 diff --git a/pyscf/dft/test/test_he.py b/pyscf/dft/test/test_he.py index 791593616d..27f6f0a08f 100644 --- a/pyscf/dft/test/test_he.py +++ b/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() diff --git a/pyscf/dft/test/test_xc_deriv.py b/pyscf/dft/test/test_xc_deriv.py index 03655f45e1..b6e430ffc2 100644 --- a/pyscf/dft/test/test_xc_deriv.py +++ b/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) diff --git a/pyscf/scf/hf.py b/pyscf/scf/hf.py index 8536176c47..37982c581f 100644 --- a/pyscf/scf/hf.py +++ b/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 diff --git a/pyscf/scf/rohf.py b/pyscf/scf/rohf.py index 777a9a5b17..f1243e9935 100644 --- a/pyscf/scf/rohf.py +++ b/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):