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

Author: sunqm

pyscf/symm/geom.py

@@ -57,9 +57,19 @@ def parallel_vectors(v1, v2, tol=TOLERANCE):
         return (abs(cos-1) < TOLERANCE) | (abs(cos+1) < TOLERANCE)
 
 def argsort_coords(coords, decimals=None):
+    # * 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.124, 0.125]
+    # 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.
     if decimals is None:
-        decimals = int(-numpy.log10(TOLERANCE)) - 1
-    coords = numpy.around(coords, decimals=decimals)
+        fac = 2**int(-numpy.log2(TOLERANCE)+.5)
+    else:
+        fac = 2**int(3.3219281 * decimals + .5)
+    coords = numpy.around(coords*fac)
     idx = numpy.lexsort((coords[:,2], coords[:,1], coords[:,0]))
     return idx
 
@@ -482,7 +492,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
@@ -521,6 +535,13 @@ 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:]
+    for op in ops:
+        if not is_identical_geometry(coords, coords.dot(op), rawsys.weights):
+            return False
+
     for lst in rawsys.atomtypes.values():
         coords = rawsys.atoms[lst,1:]
         idx = argsort_coords(coords)
@@ -540,6 +561,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[lst], weights, order) -
+               casimir_tensors(coords2[lst], 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
 
@@ -572,6 +614,7 @@ def __init__(self, atoms, basis=None):
                 fake_chgs.append([mole.charge(ksymb)] * len(lst))
         coords = numpy.array(numpy.vstack(coords), dtype=float)
         fake_chgs = numpy.hstack(fake_chgs)
+        self.weights = fake_chgs
         self.charge_center = numpy.einsum('i,ij->j', fake_chgs, coords)/fake_chgs.sum()
         coords = coords - self.charge_center
 

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)