Small fixes (pyscf#1456)

* MAINT: Consider lo.boys converged only when #iter < maxiter

* MAINT: More informative printout upon symmetrization fail

* MAINT: Handle case if no occ is available when writing molden file

* MAINT: Handle case of wrong energies assigned when writing to molden file

* MAINT: Ensure self._scf._eri is there before checking if it is None

* MAINT: Newline for more readable printout in AVAS

* MAINT: Prevent taking arccos of arguments > 1 and < -1 occuring due to finite precision

* MAINT: Automatically use mo_occ and mo_energy in chkfile.dump_mcscf

* MAINT: Allow setting pop_method upon init of PM localization

* BUG: fix integer division for python3

Co-authored-by: Xing Zhang <zhangxing.nju@gmail.com>

pyscf/lo/boys.py

@@ -78,7 +78,8 @@ def kernel(localizer, mo_coeff=None, callback=None, verbose=None):
                  imacro+1, e, de, norm_gorb, stat.tot_kf+1, stat.tot_hop)
         cput1 = log.timer('cycle= %d'%(imacro+1), *cput1)
 
-        if (norm_gorb < conv_tol_grad and abs(de) < localizer.conv_tol):
+        if (norm_gorb < conv_tol_grad and abs(de) < localizer.conv_tol
+                and stat.tot_hop < localizer.ah_max_cycle):
             conv = True
 
         if callable(callback):

pyscf/lo/ibo.py

@@ -287,17 +287,17 @@ def MakeAtomInfos():
     for At in "Li Be B C O N F Ne".split(): nCoreX[At] = 1
     for At in "Na Mg Al Si P S Cl Ar".split(): nCoreX[At] = 5
     for At in "Na Mg Al Si P S Cl Ar".split(): nCoreX[At] = 5
-    for At in "K Ca".split(): nCoreX[At] = 18/2
-    for At in "Sc Ti V Cr Mn Fe Co Ni Cu Zn".split(): nCoreX[At] = 18/2
-    for At in "Ga Ge As Se Br Kr".split(): nCoreX[At] = 18/2+5 # [Ar] and the 5 d orbitals.
+    for At in "K Ca".split(): nCoreX[At] = 18//2
+    for At in "Sc Ti V Cr Mn Fe Co Ni Cu Zn".split(): nCoreX[At] = 18//2
+    for At in "Ga Ge As Se Br Kr".split(): nCoreX[At] = 18//2 + 5  # [Ar] and the 5d orbitals.
     nAoX = {"H": 1, "He": 1}
     for At in "Li Be".split(): nAoX[At] = 2
     for At in "B C O N F Ne".split(): nAoX[At] = 5
     for At in "Na Mg".split(): nAoX[At] = 3*1 + 1*3
     for At in "Al Si P S Cl Ar".split(): nAoX[At] = 3*1 + 2*3
-    for At in "K Ca".split(): nAoX[At] = 18/2+1
-    for At in "Sc Ti V Cr Mn Fe Co Ni Cu Zn".split(): nAoX[At] = 18/2+1+5   # 4s, 3d
-    for At in "Ga Ge As Se Br Kr".split(): nAoX[At] = 18/2+1+5+3
+    for At in "K Ca".split(): nAoX[At] = 18//2 + 1
+    for At in "Sc Ti V Cr Mn Fe Co Ni Cu Zn".split(): nAoX[At] = 18//2 + 1 + 5  # 4s, 3d
+    for At in "Ga Ge As Se Br Kr".split(): nAoX[At] = 18//2 + 1 + 5 + 3
 
     AoLabels = {}
 

pyscf/lo/pipek.py

@@ -186,9 +186,11 @@ class PipekMezey(boys.Boys):
     conv_tol = getattr(__config__, 'lo_pipek_PM_conv_tol', 1e-6)
     exponent = getattr(__config__, 'lo_pipek_PM_exponent', 2)  # should be 2 or 4
 
-    def __init__(self, mol, mo_coeff=None, mf=None):
+    def __init__(self, mol, mo_coeff=None, mf=None, pop_method=None):
         boys.Boys.__init__(self, mol, mo_coeff)
         self._scf = mf
+        if pop_method is not None:
+            self.pop_method = pop_method
         self._keys = self._keys.union(['pop_method', 'exponent', '_scf'])
 
     def dump_flags(self, verbose=None):

pyscf/mcscf/avas.py

@@ -119,7 +119,7 @@ def _kernel(avas_obj):
     pmol.build(False, False)
 
     baslst = pmol.search_ao_label(avas_obj.aolabels)
-    log.info('reference AO indices for %s %s: %s',
+    log.info('reference AO indices for %s %s:\n %s',
              avas_obj.minao, avas_obj.aolabels, baslst)
 
     if avas_obj.with_iao:

pyscf/mcscf/casci.py

@@ -918,7 +918,7 @@ def ao2mo(self, mo_coeff=None):
         elif mo_coeff.shape[1] != ncas:
             mo_coeff = mo_coeff[:,ncore:nocc]
 
-        if self._scf._eri is not None:
+        if hasattr(self._scf, '_eri') and self._scf._eri is not None:
             eri = ao2mo.full(self._scf._eri, mo_coeff,
                              max_memory=self.max_memory)
         else:

pyscf/mcscf/chkfile.py

@@ -40,6 +40,8 @@ def dump_mcscf(mc, chkfile=None, key='mcscf',
     if e_tot is None: e_tot = mc.e_tot
     if e_cas is None: e_cas = mc.e_cas
     if mo_coeff is None: mo_coeff = mc.mo_coeff
+    if mo_occ is None: mo_occ = mc.mo_occ
+    if mo_energy is None: mo_energy = mc.mo_energy
     #if ci_vector is None: ci_vector = mc.ci
 
     if h5py.is_hdf5(chkfile):

pyscf/symm/Dmatrix.py

@@ -153,7 +153,13 @@ def get_euler_angles(c1, c2):
         if numpy.cross(c1[1], yp).dot(c1[2]) < 0:
             alpha = -alpha
 
-        gamma = numpy.arccos(yp.dot(c2[1]))
+        tmp = yp.dot(c2[1])
+        if abs(tmp - 1.0) < 1e-12:
+            gamma = numpy.arccos(1.0)
+        elif abs(tmp + 1.0) < 1e-12:
+            gamma = numpy.arccos(-1.0)
+        else:
+            gamma = numpy.arccos(tmp)
         if numpy.cross(yp, c2[1]).dot(c2[2]) < 0:
             gamma = -gamma
 

pyscf/symm/addons.py

@@ -221,9 +221,9 @@ def symmetrize_space(mol, mo, s=None,
             mo1.append(numpy.dot(csym, u[:,abs(1-e) < 1e-6]))
     mo1 = numpy.hstack(mo1)
     if mo1.shape[1] != nmo:
-        raise ValueError('The input orbital space is not symmetrized.\n One '
+        raise ValueError('mo1.shape[1] != nmo: %d != %d The input orbital space is not symmetrized.\n One '
                          'possible reason is that the input mol and orbitals '
-                         'are of different orientation.')
+                         'are of different orientation.' % (mo1.shape[1], nmo))
     if (check and
         abs(reduce(numpy.dot, (mo1.conj().T, s, mo1)) - numpy.eye(nmo)).max() > tol):
         raise ValueError('Orbitals are not orthogonalized')

pyscf/tools/molden.py

@@ -53,7 +53,7 @@ def orbital_coeff(mol, fout, mo_coeff, spin='Alpha', symm=None, ene=None,
                                       mo_coeff, tol=1e-5)
             except ValueError as e:
                 logger.warn(mol, str(e))
-    if ene is None:
+    if ene is None or len(ene) != nmo:
         ene = numpy.arange(nmo)
     assert (spin == 'Alpha' or spin == 'Beta')
     if occ is None:
@@ -135,14 +135,17 @@ def from_chkfile(filename, chkfile, key='scf/mo_coeff', ignore_h=IGNORE_H):
                 ene = dat['mo_energy']
             else:
                 ene = None
-            occ = dat['mo_occ']
+            if 'mo_occ' in dat:
+                occ = dat['mo_occ']
+            else:
+                occ = None
             mo = dat['mo_coeff']
 
         if isinstance(ene, str) and ene == 'None':
             ene = None
         if isinstance(ene, str) and occ == 'None':
             occ = None
-        if occ.ndim == 2:
+        if occ is not None and occ.ndim == 2:
             orbital_coeff(mol, f, mo[0], spin='Alpha', ene=ene[0], occ=occ[0],
                           ignore_h=ignore_h)
             orbital_coeff(mol, f, mo[1], spin='Beta', ene=ene[1], occ=occ[1],