Handle smearing-ROHF

Author: sunqm

pyscf/scf/smearing.py

@@ -31,12 +31,28 @@ def smearing(mf, sigma=None, method=SMEARING_METHOD, mu0=None, fix_spin=False):
         mf.fix_spin = fix_spin
         return mf
 
-    assert not mf.istype('KSCF')
-    if mf.istype('ROHF'):
-        # ROHF leads to two Fock matrices. It's not clear how to define the
-        # Roothaan effective Fock matrix from the two.
-        raise NotImplementedError('Smearing-ROHF')
+    if mf.istype('_CIAH_SOSCF'):
+        raise NotImplementedError('Smearing with second order SCF is not supported')
+
+    if mf.istype('KSCF'):
+        from pyscf.pbc.scf.smearing import smearing
+        return smearing(mf, sigma, method, mu0, fix_spin)
 
+    if mf.istype('ROHF'):
+        if fix_spin:
+            raise RuntimeError('Smearing-ROHF with fix_spin not supported. Use UHF instead.')
+        # Roothaan Fock matrix for ROHF is not supported by the smearing method.
+        # The single occupancy can be handled using the regular RHF class
+        from pyscf.scf.addons import _object_without_soscf
+        from pyscf import scf
+        from pyscf import dft
+        known_class = {
+            dft.rks_symm.ROKS: dft.rks_symm.RKS,
+            dft.roks.ROKS    : dft.rks.RKS     ,
+            scf.hf_symm.ROHF : scf.hf_symm.RHF ,
+            scf.rohf.ROHF    : scf.hf.RHF      ,
+        }
+        mf = _object_without_soscf(mf, known_class)
     return lib.set_class(_SmearingSCF(mf, sigma, method, mu0, fix_spin),
                          (_SmearingSCF, mf.__class__))
 
@@ -136,8 +152,12 @@ def get_occ(self, mo_energy=None, mo_coeff=None):
             mo_occ = super().get_occ(mo_energy, mo_coeff)
             return mo_occ
 
+        if self.istype('ROHF'):
+            # ROHF leads to two Fock matrices. It's not clear how to define the
+            # Roothaan effective Fock matrix from the two.
+            raise NotImplementedError('Smearing-ROHF')
         is_uhf = self.istype('UHF')
-        is_rhf = self.istype('RHF')
+        is_rhf = not is_uhf
 
         sigma = self.sigma
         if self.smearing_method.lower() == 'fermi':

pyscf/scf/test/test_addons.py

@@ -413,6 +413,44 @@ def test_remove_lindep(self):
                                        lindep=1e-9).run()
         self.assertAlmostEqual(mf.e_tot, -1.6291001503057689, 7)
 
+    @unittest.skip('Smearing ROHF with fix_spin does not work')
+    def test_rohf_smearing(self):
+        # Fe2 from https://doi.org/10.1021/acs.jpca.1c05769
+        mol = gto.M(
+            atom='''
+        Fe       0. 0. 0.
+        Fe       2.01 0. 0.
+        ''', basis="lanl2dz", ecp="lanl2dz", symmetry=False, unit='Angstrom',
+            spin=6, charge=0, verbose=0)
+        myhf_s = scf.ROHF(mol)
+        myhf_s = addons.smearing(myhf_s, sigma=0.01, method='gaussian', fix_spin=True)
+        myhf_s.kernel()
+        # macos py3.7 CI -242.4828982467762
+        # linux py3.11 CI -242.48289824670388
+        self.assertAlmostEqual(myhf_s.e_tot, -242.482898246, 6)
+        self.assertAlmostEqual(myhf_s.entropy, 0.45197, 4)
+        myhf2 = myhf_s.undo_smearing().newton()
+        myhf2.kernel(myhf_s.make_rdm1())
+        # FIXME: note 16mHa lower energy than myhf
+        self.assertAlmostEqual(myhf2.e_tot, -244.9808750, 6)
+
+        myhf_s.smearing_method = 'fermi'
+        myhf_s.kernel()
+        self.assertAlmostEqual(myhf_s.e_tot, -244.200255453, 6)
+        self.assertAlmostEqual(myhf_s.entropy, 3.585155, 4)
+
+    def test_rohf_smearing1(self):
+        mol = gto.M(atom = '''
+            7      0.   0  -0.7
+            7      0.   0   0.7''',
+            charge = -1,
+            spin = 1)
+        mf = mol.RHF()
+        mf = addons.smearing(mf, sigma=0.1)
+        mf.kernel()
+        self.assertAlmostEqual(mf.mo_occ.sum(), 15, 8)
+        self.assertAlmostEqual(mf.e_tot, -106.9310800402, 8)
+
     def test_uhf_smearing(self):
         mol = gto.M(
             atom='''