GHF with fractional occupancy (pyscf#1826)

* Fix bug in GHF with fractional occupancy

* Remove unused code in X2C

* Fix setup.py

pyscf/scf/addons.py

@@ -48,7 +48,7 @@ def frac_occ_(mf, tol=1e-3):
         >>> mf = scf.addons.frac_occ(mf)
         >>> mf.run()
     '''
-    from pyscf.scf import uhf, rohf
+    from pyscf.scf import hf, uhf, rohf
     old_get_occ = mf.get_occ
     mol = mf.mol
 
@@ -145,7 +145,7 @@ def get_grad(mo_coeff, mo_occ, fock):
             g[uniq_var_b] += fockb[uniq_var_b]
             return g[uniq_var_a | uniq_var_b]
 
-    else:  # RHF
+    elif isinstance(mf, hf.RHF):
         def get_occ(mo_energy, mo_coeff=None):
             nocc = (mol.nelectron+1) // 2  # n_docc + n_socc
             mo_occ, frac_lst, homo, lumo = guess_occ(mo_energy, nocc)
@@ -162,6 +162,18 @@ def get_occ(mo_energy, mo_coeff=None):
             else:
                 mo_occ = old_get_occ(mo_energy, mo_coeff)
             return mo_occ
+    else:
+        def get_occ(mo_energy, mo_coeff=None):
+            nocc = mol.nelectron
+            mo_occ, frac_lst, homo, lumo = guess_occ(mo_energy, nocc)
+            if abs(homo - lumo) < tol:
+                logger.warn(mf, 'fraction occ = %6g  for orbitals %s',
+                            mo_occ[frac_lst[0]], frac_lst)
+                logger.info(mf, 'HOMO = %.12g  LUMO = %.12g', homo, lumo)
+                logger.debug(mf, '  mo_energy = %s', mo_energy)
+            else:
+                mo_occ = old_get_occ(mo_energy, mo_coeff)
+            return mo_occ
 
     mf.get_occ = get_occ
     if get_grad is not None:

pyscf/x2c/x2c.py

@@ -758,13 +758,12 @@ def dip_moment(self, mol=None, dm=None, unit='Debye', verbose=logger.NOTE,
             log = logger.new_logger(mol, verbose)
 
             nao = mol.nao
-            dm = dm[:nao,:nao] + dm[nao:,nao:]
-
             charges = mol.atom_charges()
             coords  = mol.atom_coords()
             nucl_dip = numpy.einsum('i,ix->x', charges, coords)
             with mol.with_common_orig(nucl_dip):
-                ao_dip = _block_diag(mol.intor_symmetric('int1e_r'))
+                r = mol.intor_symmetric('int1e_r')
+                ao_dip = numpy.array([_block_diag(x) for x in r])
                 if picture_change:
                     xmol = self.with_x2c.get_xmol()[0]
                     nao = xmol.nao

setup.py

@@ -96,7 +96,7 @@ def build_cmake(self, extension):
         self.announce('Building binaries', level=3)
         # Do not use high level parallel compilation. OOM may be triggered
         # when compiling certain functionals in libxc.
-        cmd = ['cmake', '--build', self.build_temp, '-j2']
+        cmd = ['cmake', '--build', self.build_temp, '-j', '2']
         build_args = os.getenv('CMAKE_BUILD_ARGS')
         if build_args:
             cmd.extend(build_args.split(' '))