Symmetry detection bug due to np.round tolerance (pyscf#2715)

* Fix point-group symmetry verficication due to the np.round issue

* indexing bug

* Fix symmetry detection tests

* Fix argsort_coords

* Fix bug in sort_coords

* Update comments

Author: sunqm

pyscf/symm/geom.py

@@ -56,18 +56,31 @@ def parallel_vectors(v1, v2, tol=TOLERANCE):
         cos = numpy.dot(_normalize(v1), _normalize(v2))
         return (abs(cos-1) < TOLERANCE) | (abs(cos+1) < TOLERANCE)
 
-def argsort_coords(coords, decimals=None):
+def argsort_coords(coords, decimals=None, tol=0.05):
+    # * np.round for decimal places can lead to more errors than the actual
+    # difference between two numbers. For example,
+    # np.round([0.1249999999,0.1250000001], 2) => [0.12, 0.13]
+    # np.round([0.1249999999,0.1250000001], 3) => [0.125, 0.125]
+    # When loosen tolerance is used, compared to the more strict tolerance,
+    # the coordinates might look more different.
+    # * Using the power of two as the factor can reduce such errors, although not
+    # faithfully rounding to the required decimals.
+    # * For normal molecules, tol~=0.1 in coordinates is enough to distinguish
+    # atoms in molecule. A very tight threshold is not appropriate here. With
+    # tight threshold, small differences in coordinates may lead to different
+    # arg orders.
     if decimals is None:
-        decimals = int(-numpy.log10(TOLERANCE)) - 1
-    coords = numpy.around(coords, decimals=decimals)
+        fac = 2**int(-numpy.log2(tol))
+    else:
+        fac = 2**int(3.3219281 * decimals)
+    # +.5 for rounding to the nearest integer
+    coords = (coords*fac + .5).astype(int)
     idx = numpy.lexsort((coords[:,2], coords[:,1], coords[:,0]))
     return idx
 
-def sort_coords(coords, decimals=None):
-    if decimals is None:
-        decimals = int(-numpy.log10(TOLERANCE)) - 1
+def sort_coords(coords, decimals=None, tol=0.05):
     coords = numpy.asarray(coords)
-    idx = argsort_coords(coords, decimals=decimals)
+    idx = argsort_coords(coords, tol=tol)
     return coords[idx]
 
 # ref. http://en.wikipedia.org/wiki/Rotation_matrix
@@ -482,7 +495,11 @@ def symm_identical_atoms(gpname, atoms):
         newc = numpy.dot(coords, op)
         idx = argsort_coords(newc)
         if not numpy.allclose(coords0, newc[idx], atol=TOLERANCE):
-            raise PointGroupSymmetryError('Symmetry identical atoms not found')
+            raise PointGroupSymmetryError(
+                'Symmetry identical atoms not found. This may be due to '
+                'the strict setting of the threshold symm.geom.TOLERANCE. '
+                'Consider adjusting the tolerance.')
+
         dup_atom_ids.append(idx)
 
     dup_atom_ids = numpy.sort(dup_atom_ids, axis=0).T
@@ -524,6 +541,14 @@ def check_symm(gpname, atoms, basis=None):
     opdic = symm_ops(gpname)
     ops = [opdic[op] for op in OPERATOR_TABLE[gpname]]
     rawsys = SymmSys(atoms, basis)
+
+    # A fast check using Casimir tensors
+    coords = rawsys.atoms[:,1:]
+    weights = rawsys.atoms[:,0]
+    for op in ops:
+        if not is_identical_geometry(coords, coords.dot(op), weights):
+            return False
+
     for lst in rawsys.atomtypes.values():
         coords = rawsys.atoms[lst,1:]
         idx = argsort_coords(coords)
@@ -543,6 +568,27 @@ def shift_atom(atoms, orig, axis):
     c = numpy.dot(c - orig, numpy.array(axis).T)
     return [[atoms[i][0], c[i]] for i in range(len(atoms))]
 
+def is_identical_geometry(coords1, coords2, weights):
+    '''A fast check to compare the geometry of two molecules using Casimir tensors'''
+    if coords1.shape != coords2.shape:
+        return False
+    for order in range(1, 4):
+        if abs(casimir_tensors(coords1, weights, order) -
+               casimir_tensors(coords2, weights, order)).max() > TOLERANCE:
+            return False
+    return True
+
+def casimir_tensors(r, q, order=1):
+    if order == 1:
+        return q.dot(r)
+    elif order == 2:
+        return numpy.einsum('i,ix,iy->xy', q, r, r)
+    elif order == 3:
+        return numpy.einsum('i,ix,iy,iz->xyz', q, r, r, r)
+    else:
+        raise NotImplementedError
+
+
 class RotationAxisNotFound(PointGroupSymmetryError):
     pass
 

pyscf/symm/test/test_basis.py

@@ -141,6 +141,11 @@ def test_symm_orb_c1(self):
         self.assertEqual(get_so(atoms,basis)[0], 220)
 
     def test_symm_orb_c3v_as_cs(self):
+        # This molecule has an approximate C3v symmetry for TOLERANCE>1e-5
+        # When using its subgroup cs, a mirror perpedicular to [-.5, -.866, 0]
+        # will be adopted as the mirror. This mirror leads to more errors in
+        # atomic coordinates than the yz-mirror. In this case, one can adjust
+        # TOLERANCE to pass the check_symm or symm_identical_atoms functions.
         atoms = [['Fe', ( 0.000000  , 0.000000  , 0.015198 )],
                  ['C',  ( 0.000000  , 0.000000  , -1.938396)],
                  ['C',  ( 0.000000  , -1.394127 , -1.614155)],
@@ -158,6 +163,10 @@ def test_symm_orb_c3v_as_cs(self):
                  ['O',  ( 0.000000  , 2.572496  , 1.441607 )],
                  ['O',  ( 2.227847  , -1.286248 , 1.441607 )],
                  ['O',  ( -2.227847 , -1.286248 , 1.441607 )],]
+        numpy.random.seed(2)
+        u = numpy.linalg.svd(numpy.random.rand(3,3))[0]
+        r = numpy.array([a[1] for a in atoms])
+        atoms = [[a[0], x] for a, x in zip(atoms, r.dot(u))]
         basis = {'Fe':gto.basis.load('def2svp', 'C'),
                  'C': gto.basis.load('def2svp', 'C'),
                  'H': gto.basis.load('def2svp', 'C'),

pyscf/symm/test/test_geom.py

@@ -886,6 +886,16 @@ def test_c2v_shifted(self):
         self.assertEqual(l, 'C2v')
         self.assertAlmostEqual(abs(axes - numpy.diag(axes.diagonal())).max(), 0, 9)
 
+    def test_geometry_small_discrepancy(self):
+        # issue 2713
+        mol = gto.M(
+            atom='''
+                O        0.000000    0.000000    0.900000
+                H        0.000000    0.000000    0.000000
+                H        0.914864    0.000000    1.249646
+                H       -0.498887    0.766869    1.249646''',
+            charge=1, symmetry=True)
+        self.assertEqual(mol.groupname, 'Cs')
 
 def ring(n, start=0):
     r = 1. / numpy.sin(numpy.pi/n)