diff --git a/HISTORY.rst b/HISTORY.rst
index 6b181c2..f20b91f 100644
--- a/HISTORY.rst
+++ b/HISTORY.rst
@@ -1,6 +1,10 @@
 =======
 History
 =======
+2023.4.6 -- Added gradients
+    * Added gradient on atoms as a separate table alongside atoms, so they take no space
+      unless actually used.
+      
 2024.3.13 -- Handle uppercase X, Y, Z in strings for symmetry operators
     * the Crystallographic Open Database CIF files seems to use upper case X, Y, Z in
       explicit symmetry operators. These need to be lowercased in the code.
diff --git a/molsystem/atoms.py b/molsystem/atoms.py
index 52f888d..8ee0f80 100644
--- a/molsystem/atoms.py
+++ b/molsystem/atoms.py
@@ -46,6 +46,7 @@ def __init__(self, configuration) -> None:
         self._atom_table = _Table(self.system_db, "atom")
         self._coordinates_table = _Table(self.system_db, "coordinates")
         self._velocities_table = _Table(self.system_db, "velocities")
+        self._gradients_table = _Table(self.system_db, "gradients")
 
         super().__init__(self._system_db, "atom")
 
@@ -70,6 +71,8 @@ def __delitem__(self, key) -> None:
             del self._coordinates_table[key]
         elif key in self._velocities_table.attributes:
             del self._velocities_table[key]
+        elif key in self._gradients_table.attributes:
+            del self._gradients_table[key]
 
     def __iter__(self) -> iter:
         """Allow iteration over the object"""
@@ -184,8 +187,8 @@ def atom_generators(self):
     @property
     def attributes(self) -> Dict[str, Any]:
         """The definitions of the attributes.
-        Combine the attributes of the atom, coordinates, and velocities tables to
-        make it look like a single larger table.
+        Combine the attributes of the atom, coordinates, velocities, and gradients
+        tables to make it look like a single larger table.
         """
 
         result = self._atom_table.attributes
@@ -198,6 +201,10 @@ def attributes(self) -> Dict[str, Any]:
             if key != "atom":  # atom key links the tables together, so ignore
                 result[key] = value
 
+        for key, value in self._gradients_table.attributes.items():
+            if key != "atom":  # atom key links the tables together, so ignore
+                result[key] = value
+
         return result
 
     @property
@@ -225,6 +232,16 @@ def db(self):
         """The database connection."""
         return self.system_db.db
 
+    @property
+    def gradients(self):
+        """The gradients as list of lists."""
+        return self.get_gradients()
+
+    @gradients.setter
+    def gradients(self, xyz):
+        """The gradients as list of lists."""
+        return self.set_gradients(xyz)
+
     @property
     def group(self):
         """The space or point group of the configuration."""
@@ -234,6 +251,19 @@ def group(self):
     def group(self, value):
         self.configuration.symmetry.group = value
 
+    @property
+    def have_gradients(self):
+        """Whether there are any gradients for this configuration."""
+        sql = (
+            "SELECT COUNT(*)"
+            " FROM atomset_atom AS aa, gradients AS ve "
+            "WHERE aa.atomset = ?"
+            "  AND ve.atom = aa.atom AND ve.configuration = ?"
+        )
+        parameters = [self.atomset, self.configuration.id]
+        self.cursor.execute(sql, parameters)
+        return self.cursor.fetchone()[0] > 0
+
     @property
     def have_velocities(self):
         """Whether there are any velocities for this configuration."""
@@ -476,6 +506,16 @@ def append(self, **kwargs: Dict[str, Any]) -> None:
         if have_velocities:
             self._velocities_table.append(n=n_rows, **data)
 
+        # And gradients table
+        data = {"configuration": configuration, "atom": ids}
+        have_gradients = False
+        for column in self._gradients_table.attributes:
+            if column != "id" and column in kwargs:
+                data[column] = kwargs.pop(column)
+                have_gradients = True
+        if have_gradients:
+            self._gradients_table.append(n=n_rows, **data)
+
         # And to the atomset
         table = _Table(self.system_db, "atomset_atom")
         table.append(atomset=self.atomset, atom=ids)
@@ -497,88 +537,79 @@ def atoms(self, *args):
         sqlite3.Cursor
             A cursor that returns sqlite3.Row objects for the atoms.
         """
-        if self.have_velocities:
-            columns = self._columns(velocities=True)
-            column_defs = ", ".join(columns)
-
-            # What tables are requested in the extra arguments?
-            tables = set()
-            if len(args) > 0:
-                atom_columns = [*self._atom_table.attributes]
-                coordinates_columns = [*self._coordinates_table.attributes]
-                velocities_columns = [*self._velocities_table.attributes]
-                for col, op, value in grouped(args, 3):
-                    if "." in col:
-                        tables.add(col.split(".")[0])
-                    elif col in atom_columns:
-                        tables.add("at")
-                    elif col in coordinates_columns:
-                        tables.add("co")
-                    elif col in velocities_columns:
-                        tables.add("ve")
-                    else:
-                        raise ValueError(f"Column '{col}' is not available")
+        have_velocities = self.have_velocities
+        have_gradients = self.have_gradients
+        columns = self._columns(velocities=have_velocities, gradients=have_gradients)
+        column_defs = ", ".join(columns)
 
-            # Build the query based on the tables needed
-            sql = (
-                f"SELECT {column_defs}"
-                " FROM atomset_atom AS aa, atom AS at, coordinates AS co, "
-                "      velocities AS ve "
-                "WHERE aa.atomset = ?"
-                "  AND at.id = aa.atom"
-                "  AND co.atom = at.id AND co.configuration = ?"
-                "  AND ve.atom = at.id AND ve.configuration = ?"
-            )
-            parameters = [self.atomset, self.configuration.id, self.configuration.id]
-
-            # And any extra selection criteria
-            if len(args) > 0:
-                for col, op, value in grouped(args, 3):
-                    if op == "==":
-                        op = "="
-                    sql += f' AND "{col}" {op} ?'
-                    parameters.append(value)
+        # What tables are requested in the extra arguments?
+        tables = set()
+        if len(args) == 0:
+            tables.add("at")
+            tables.add("co")
+            if have_velocities:
+                tables.add("ve")
+            if have_gradients:
+                tables.add("gr")
         else:
-            columns = self._columns(velocities=False)
-            column_defs = ", ".join(columns)
-
-            # What tables are requested in the extra arguments?
-            tables = set()
-            if len(args) > 0:
-                atom_columns = [*self._atom_table.attributes]
-                coordinates_columns = [*self._coordinates_table.attributes]
-                velocities_columns = [*self._velocities_table.attributes]
-                for col, op, value in grouped(args, 3):
-                    if "." in col:
-                        tables.add(col.split(".")[0])
-                    elif col in atom_columns:
-                        tables.add("at")
-                    elif col in coordinates_columns:
-                        tables.add("co")
-                    elif col in velocities_columns:
+            atom_columns = [*self._atom_table.attributes]
+            coordinates_columns = [*self._coordinates_table.attributes]
+            velocities_columns = [*self._velocities_table.attributes]
+            gradients_columns = [*self._gradients_table.attributes]
+            for col, op, value in grouped(args, 3):
+                if "." in col:
+                    tables.add(col.split(".")[0])
+                elif col in atom_columns:
+                    tables.add("at")
+                elif col in coordinates_columns:
+                    tables.add("co")
+                elif col in velocities_columns:
+                    if have_velocities:
+                        tables.add("ve")
+                    else:
                         raise ValueError(
                             "Query for atom has velocities, but the atoms don't"
                         )
+                elif col in gradients_columns:
+                    if have_gradients:
+                        tables.add("gr")
                     else:
-                        raise ValueError(f"Column '{col}' is not available")
+                        raise ValueError(
+                            "Query for atom has gradients, but the atoms don't"
+                        )
+                else:
+                    raise ValueError(f"Column '{col}' is not available")
 
-            # Build the query based on the tables needed
-            sql = (
-                f"SELECT {column_defs}"
-                " FROM atomset_atom AS aa, atom AS at, coordinates AS co "
-                "WHERE aa.atomset = ?"
-                "  AND at.id = aa.atom"
-                "  AND co.atom = at.id AND co.configuration = ?"
-            )
-            parameters = [self.atomset, self.configuration.id]
+        # Build the query based on the tables needed
+        from_string = ["atomset_atom AS aa", "atom AS at", "coordinates AS co"]
+        if "ve" in tables:
+            from_string.append("velocities AS ve")
+        if "gr" in tables:
+            from_string.append("gradients AS gr")
+        from_string = ", ".join(from_string)
 
-            # And any extra selection criteria
-            if len(args) > 0:
-                for col, op, value in grouped(args, 3):
-                    if op == "==":
-                        op = "="
-                    sql += f' AND "{col}" {op} ?'
-                    parameters.append(value)
+        sql = (
+            f"SELECT {column_defs}\n"
+            f" FROM {from_string}\n"
+            "WHERE aa.atomset = ?\n"
+            "  AND at.id = aa.atom\n"
+            "  AND co.atom = at.id AND co.configuration = ?"
+        )
+        parameters = [self.atomset, self.configuration.id]
+        if "ve" in tables:
+            sql += "\n  AND ve.atom = at.id AND ve.configuration = ?"
+            parameters.append(self.configuration.id)
+        if "gr" in tables:
+            sql += "\n  AND gr.atom = at.id AND gr.configuration = ?"
+            parameters.append(self.configuration.id)
+
+        # And any extra selection criteria
+        if len(args) > 0:
+            for col, op, value in grouped(args, 3):
+                if op == "==":
+                    op = "="
+                sql += f'\n AND "{col}" {op} ?'
+                parameters.append(value)
 
         logger.debug("atoms query:")
         logger.debug(sql)
@@ -590,7 +621,7 @@ def atoms(self, *args):
     def diff(self, other):
         """Difference between these atoms and another
 
-        Currently ignores velocities. Not sure what we want to do....
+        Currently ignores velocities and gradients. Not sure what we want to do....
 
         Parameters
         ----------
@@ -605,8 +636,8 @@ def diff(self, other):
         result = {}
 
         # Check the columns
-        columns = self._columns(velocities=False)
-        other_columns = other._columns(velocities=False)
+        columns = self._columns(velocities=False, gradients=False)
+        other_columns = other._columns(velocities=False, gradients=False)
 
         column_defs = ", ".join(columns)
         other_column_defs = ", ".join(other_columns)
@@ -1152,6 +1183,58 @@ def get_property(self, name, asymmetric=False):
         # Expand to the asymmetric atoms
         return [data[i] for i in symmetry.atom_to_asymmetric_atom]
 
+    def get_gradients(
+        self,
+        fractionals=True,
+        as_array=False,
+    ):
+        """Return the gradients.
+
+        Symmetry is not supported, because it makes no (little?) sense for gradients.
+
+        Parameters
+        ----------
+        fractionals : bool = True
+            Return the gradients as fractional gradients for periodic
+            systems. Non-periodic systems always use Cartesian gradients.
+        as_array : bool = False
+            Whether to return the results as a np array or as a list of
+            lists (the default).
+
+        Returns
+        -------
+        abc : [N][float*3]
+            The gradients, either Cartesian or fractional
+        """
+        gxs = self.get_column_data("gx")
+        gys = self.get_column_data("gy")
+        gzs = self.get_column_data("gz")
+        xyz = [[gx, gy, gz] for gx, gy, gz in zip(gxs, gys, gzs)]
+
+        periodicity = self.configuration.periodicity
+        if periodicity == 0:
+            if as_array:
+                return np.array(xyz)
+            else:
+                return xyz
+
+        cell = self.configuration.cell
+
+        if fractionals:
+            if self.configuration.coordinate_system == "Cartesian":
+                return cell.to_fractionals(xyz, as_array=as_array)
+            elif as_array:
+                return np.array(xyz)
+            else:
+                return xyz
+        else:
+            if self.configuration.coordinate_system == "fractional":
+                return cell.to_cartesians(xyz, as_array=as_array)
+            elif as_array:
+                return np.array(xyz)
+            else:
+                return xyz
+
     def get_velocities(
         self,
         fractionals=True,
@@ -1276,6 +1359,77 @@ def set_coordinates(self, xyz, fractionals=True):
         y_column[0:] = ys
         z_column[0:] = zs
 
+    def set_gradients(self, gxyz, fractionals=False):
+        """Set the gradients to new values.
+
+        Parameters
+        ----------
+        fractionals : bool = False
+            The gradients are fractional gradients for periodic
+            systems. Ignored for non-periodic systems.
+
+        Returns
+        -------
+        None
+        """
+        if self.n_symops > 1:
+            raise RuntimeError("Can't handle gradients with symmetry.")
+
+        as_array = isinstance(gxyz, np.ndarray)
+
+        gxs = []
+        gys = []
+        gzs = []
+
+        periodicity = self.configuration.periodicity
+        coordinate_system = self.configuration.coordinate_system
+        if (
+            periodicity == 0
+            or (coordinate_system == "Cartesian" and not fractionals)
+            or (coordinate_system == "fractional" and fractionals)
+        ):
+            if as_array:
+                for gx, gy, gz in gxyz.tolist():
+                    gxs.append(gx)
+                    gys.append(gy)
+                    gzs.append(gz)
+            else:
+                for gx, gy, gz in gxyz:
+                    gxs.append(gx)
+                    gys.append(gy)
+                    gzs.append(gz)
+        else:
+            cell = self.configuration.cell
+            if coordinate_system == "fractional":
+                # Convert gradients to fractionals
+                for gx, gy, gz in cell.to_fractionals(gxyz):
+                    gxs.append(gx)
+                    gys.append(gy)
+                    gzs.append(gz)
+            else:
+                for gx, gy, gz in cell.to_cartesians(gxyz):
+                    gxs.append(gx)
+                    gys.append(gy)
+                    gzs.append(gz)
+
+        gx_column = self.get_column("gx")
+        if len(gx_column) == 0:
+            # No gradients in the database, so need to add rather than setFormatter
+            self._gradients_table.append(
+                n=len(gxs),
+                gx=gxs,
+                gy=gys,
+                gz=gzs,
+                atom=self.ids,
+                configuration=self.configuration.id,
+            )
+        else:
+            gy_column = self.get_column("gy")
+            gz_column = self.get_column("gz")
+            gx_column[0:] = gxs
+            gy_column[0:] = gys
+            gz_column[0:] = gzs
+
     def set_velocities(self, vxyz, fractionals=False):
         """Set the velocities to new values.
 
@@ -1391,6 +1545,18 @@ def get_column(self, key: str) -> Any:
                 f"   AND aa.atomset = {self.atomset}"
             )
             return _Column(self._velocities_table, key, sql=sql)
+        elif key in self._gradients_table.attributes:
+            sql = (
+                f'SELECT gr.rowid, gr."{key}"'
+                "   FROM atom as at,"
+                "        gradients as gr,"
+                "        atomset_atom as aa"
+                "  WHERE gr.atom = at.id"
+                f"   AND gr.configuration = {self.configuration.id}"
+                "    AND at.id = aa.atom"
+                f"   AND aa.atomset = {self.atomset}"
+            )
+            return _Column(self._gradients_table, key, sql=sql)
         else:
             raise KeyError(f"'{key}' not in atoms")
 
@@ -1438,6 +1604,18 @@ def get_column_data(self, key: str) -> Any:
                 f"   AND aa.atomset = {self.atomset}"
             )
             return [row[0] for row in self.db.execute(sql)]
+        elif key in self._gradients_table.attributes:
+            sql = (
+                f'SELECT gr."{key}"'
+                "   FROM atom as at,"
+                "        gradients as gr,"
+                "        atomset_atom as aa"
+                "  WHERE gr.atom = at.id"
+                f"   AND gr.configuration = {self.configuration.id}"
+                "    AND at.id = aa.atom"
+                f"   AND aa.atomset = {self.atomset}"
+            )
+            return [row[0] for row in self.db.execute(sql)]
         else:
             raise KeyError(f"'{key}' not in atoms")
 
@@ -1501,8 +1679,8 @@ def delete(self, atoms) -> int:
         -------
         None
         """
-        # Delete the listed atoms, which will cascade to delete coordinates and
-        # velocities
+        # Delete the listed atoms, which will cascade to delete coordinates,
+        # velocities, and gradients
         if atoms == "all" or atoms == "*":
             sql = """
             DELETE FROM atom
@@ -1534,10 +1712,11 @@ def to_dataframe(self):
 
         return df
 
-    def _columns(self, velocities=True):
-        """The list of columns across the atom, coordinates and velocities tables.
+    def _columns(self, velocities=True, gradients=True):
+        """The list of columns across the atom, coordinates, velocities, and gradients
+        tables.
 
-        Uses 'at', 'co', and 've' as the shorthand for the full table names.
+        Uses 'at', 'co', 've', and 'gr' as the shorthand for the full table names.
         """
         atom_columns = [*self._atom_table.attributes]
         coordinates_columns = [*self._coordinates_table.attributes]
@@ -1548,6 +1727,10 @@ def _columns(self, velocities=True):
             velocities_columns = [*self._velocities_table.attributes]
             velocities_columns.remove("atom")
             columns += [f've."{x}"' for x in velocities_columns]
+        if gradients:
+            gradients_columns = [*self._gradients_table.attributes]
+            gradients_columns.remove("atom")
+            columns += [f'gr."{x}"' for x in gradients_columns]
 
         return columns
 
@@ -1725,51 +1908,43 @@ def atoms(self, *args):
         sqlite3.Cursor
             A cursor that returns sqlite3.Row objects for the atoms.
         """
-        if self.have_velocities:
-            columns = self._columns()
-            column_defs = ", ".join(columns)
-
-            sql = f"""
-            SELECT {column_defs}
-              FROM atom as at, coordinates as co, velocities as ve, subset_atom as sa
-             WHERE co.atom = at.id
-               AND co.configuration = ?
-               AND ve.atom = at.id
-               AND ve.configuration = ?
-               AND at.id = sa.atom
-               AND sa.subset = ?
-            """
+        have_velocities = self.have_velocities
+        have_gradients = self.have_gradients
+        columns = self._columns(velocities=have_velocities, gradients=have_gradients)
+        column_defs = ", ".join(columns)
 
-            parameters = [self.configuration.id, self.configuration.id, self.subset_id]
-            if len(args) > 0:
-                for col, op, value in grouped(args, 3):
-                    if op == "==":
-                        op = "="
-                    sql += f' AND "{col}" {op} ?'
-                    parameters.append(value)
-        else:
-            columns = self._columns(velocities=False)
-            column_defs = ", ".join(columns)
+        from_string = ["atom as at", "coordinates as co", "subset_atom as sa"]
+        if have_velocities:
+            from_string.append("velocities AS ve")
+        if have_gradients:
+            from_string.append("gradients AS gr")
+        from_string = ", ".join(from_string)
 
-            sql = f"""
-            SELECT {column_defs}
-              FROM atom as at, coordinates as co, subset_atom as sa
-             WHERE co.atom = at.id
-               AND co.configuration = ?
-               AND at.id = sa.atom
-               AND sa.subset = ?
-            """
+        sql = f"""
+        SELECT {column_defs}
+          FROM {from_string}
+         WHERE co.atom = at.id
+           AND co.configuration = ?
+           AND at.id = sa.atom
+           AND sa.subset = ?
+        """
+        parameters = [self.configuration.id, self.subset_id]
+        if have_velocities:
+            sql += "\n  AND ve.atom = at.id AND ve.configuration = ?"
+            parameters.append(self.configuration.id)
+        if have_gradients:
+            sql += "\n  AND gr.atom = at.id AND gr.configuration = ?"
+            parameters.append(self.configuration.id)
 
-            parameters = [self.configuration.id, self.subset_id]
-            if len(args) > 0:
-                for col, op, value in grouped(args, 3):
-                    if op == "==":
-                        op = "="
-                    sql += f' AND "{col}" {op} ?'
-                    parameters.append(value)
+        if len(args) > 0:
+            for col, op, value in grouped(args, 3):
+                if op == "==":
+                    op = "="
+                sql += f'\n AND "{col}" {op} ?'
+                parameters.append(value)
 
-            if self.template.is_full and self.template_order:
-                sql += "ORDER BY sa.templateatom"
+        if self.template.is_full and self.template_order:
+            sql += "\nORDER BY sa.templateatom"
 
         return self.db.execute(sql, parameters)
 
diff --git a/molsystem/cif.py b/molsystem/cif.py
index 2b0c0c8..3cad22f 100644
--- a/molsystem/cif.py
+++ b/molsystem/cif.py
@@ -139,7 +139,7 @@ def to_cif_text(self):
             elif symmetry.n_symops == 1:
                 lines.append("_space_group_name_H-M_full   'P 1'")
             else:
-                spgname_system = symmetry.spacegroup_names_to_system[spgname]
+                spgname_system = symmetry.spacegroup_names_to_system()[spgname]
                 lines.append(f"_space_group_{spgname_system}       '{spgname}'")
                 hall = symmetry.hall_symbol
                 lines.append(f"_symmetry_space_group_name_Hall '{hall}'")
@@ -285,6 +285,7 @@ def from_cif_text(self, text):
                 self.symmetry.symops = data_block[symdata]
                 used_symops = True
             else:
+                found = False
                 for section in (
                     "_space_group" + dot + "name_Hall",
                     "_space_group" + dot + "name_H-M_full",
@@ -293,9 +294,9 @@ def from_cif_text(self, text):
                     "_symmetry" + dot + "space_group_name_H-M",
                 ):
                     if section in data_block:
+                        found = True
                         self.symmetry.group = data_block[section]
-                        break
-                else:
+                if not found:
                     raise RuntimeError(
                         "CIF file does not contain required symmetry information. "
                         "Neither "
diff --git a/molsystem/symmetry.py b/molsystem/symmetry.py
index 311756b..04287bd 100644
--- a/molsystem/symmetry.py
+++ b/molsystem/symmetry.py
@@ -10,6 +10,7 @@
 
 logger = logging.getLogger(__name__)
 # logger.setLevel("INFO")
+logger.setLevel("WARNING")
 
 
 class _Symmetry(object):
@@ -47,6 +48,337 @@ def __init__(self, configuration):
 
         super().__init__()
 
+    @staticmethod
+    def filter_atoms(atoms, coordinates, group=None, sym_ops=None, operators=None):
+        """Check and filter the atoms in a periodic system.
+
+        There is an issue with some CIF files, notably those from the Cambridge
+        structural database, where there are atoms in the asymmetric unit that
+        are duplicates. It appears that the CIF file is written with entire molecules
+        in some cases. So if the molecule is on a symmetry element, the atoms that are
+        related by symmetry are explicitly in the file.
+
+        To handle this we need to check for duplicates and remove them. Not a bad check
+        to have, anyway.
+
+        Parameters
+        ----------
+        atoms : str or int
+            The atomic symbols or numbers of the atoms.
+        coordinates : array_like
+            The coordinates of the atoms.
+        group : str
+            The space group symbol.
+        sym_ops : array_like
+            The symmetry operations as strings.
+        operators : array_like
+            The symmetry operators.
+        """
+        if group is not None:
+            # Get the symmetry operators
+            if sym_ops is not None or operators is not None:
+                raise RuntimeError(
+                    "Cannot specify both group and sym_ops or operators."
+                )
+            operators = _Symmetry.get_operators(group)
+        elif sym_ops is not None:
+            if operators is not None or group is not None:
+                raise RuntimeError(
+                    "Cannot specify both sym_ops and group or operators."
+                )
+            operators = _Symmetry.get_operators(sym_ops)
+        elif operators is not None:
+            if sym_ops is not None or group is not None:
+                raise RuntimeError(
+                    "Cannot specify both operators and group or sym_ops."
+                )
+        else:
+            raise RuntimeError("Must specify either group or sym_ops or operators.")
+
+        self = _Symmetry()
+
+        symbols = self.configuration.atoms.asymmetric_symbols
+        n_atoms = self.configuration.n_asymmetric_atoms
+        logger.info(f"Expanding {n_atoms} asymmetric atoms to full cell.")
+        if n_atoms == 0:
+            self._atom_generators = None
+            self._symop_to_atom = None
+        self._atom_generators = []  # Symmetry ops that create the symmetric atoms
+        self._symop_to_atom = []  # The symmetry atom resulting from symmetry ops
+        if self.configuration.periodicity == 3:
+            if len(operators) == 1:
+                # P1 is a special case. Nothing to do.
+                self._atom_generators = [[0] for i in range(n_atoms)]
+                self._symop_to_atom = [[i] for i in range(n_atoms)]
+                self._atom_to_asymmetric_atom = [*range(n_atoms)]
+            else:
+                uvw0 = self.configuration.atoms.get_coordinates(
+                    as_array=True, asymmetric=True
+                )
+                if uvw0.shape[0] != n_atoms:
+                    raise RuntimeError(
+                        f"Mismatch of number of atoms in symmetry: {uvw0.shape[0]} != "
+                        f"{n_atoms}"
+                    )
+                logger.debug(
+                    f"""
+{n_atoms=}
+{uvw0.shape=}"
+Original coordinates
+{uvw0}
+                    """
+                )
+
+                uvw = np.ndarray((n_atoms, 4))
+                uvw[:, 0:3] = uvw0[:, :]
+                logger.debug(f"\nExpanded coordinates\n{uvw}")
+                uvw[:, 3] = 1
+                # logger.debug(self.symops)
+                logger.debug(f"\n{operators.shape=}")
+                # logger.debug(operators)
+                # logger.debug(f"{uvw.shape=}")
+                # logger.debug("Expanded coordinates")
+                # logger.debug(uvw)
+                xformed = np.einsum("ijk,lk", operators, uvw)
+                logger.debug(f"\n{xformed.shape=}\n{xformed}")
+
+                # For comparison, bring all atoms into the cell [0..1)
+                tmp = xformed - np.floor(xformed)
+
+                logger.debug(f"\ntmp = Transformed and floored coordinates\n{tmp}")
+
+                # Keep track of atoms and coordinates found so can check for duplicates
+                found = {}
+
+                n_sym_atoms = 0
+                n_asymmetric_atoms = 0
+                for i in range(n_atoms):
+                    values, I1, I2 = np.unique(
+                        np.round(tmp[:, :3, i], 4),
+                        axis=0,
+                        return_index=True,
+                        return_inverse=True,
+                    )
+                    # pprint.pprint(np.round(tmp[:, :3, i], 4).tolist())
+                    logger.debug(
+                        f"""
+{i=}: {symbols[i]}
+{np.round(tmp[:, :3, i], 4)}
+nvalues
+{values}
+I1
+{I1}
+I2
+{I2}
+                        """
+                    )
+                    # Check for duplicates
+                    if tuple(values[0]) in found:
+                        if found[tuple(values[0])] != symbols[i]:
+                            raise RuntimeError(
+                                "Duplicate atoms with different symbols: "
+                                f"{symbols[i]} != {found[tuple(values[0])]} at "
+                                f"{values[0]}"
+                            )
+                        else:
+                            # Same element, so ignore
+                            logger.debug("\nDuplicate atom, ignoring")
+                            continue
+
+                    for value in values:
+                        found[tuple(value)] = symbols[i]
+
+                    n_asymmetric_atoms += 1
+                    I2 += n_sym_atoms
+                    self._atom_generators.append(I1.tolist())
+                    self._symop_to_atom.append(I2.tolist())
+                    n_sym_atoms += I1.shape[0]
+                tmp = []
+                for i, generators in enumerate(self._atom_generators):
+                    tmp.extend([i] * len(generators))
+                self._atom_to_asymmetric_atom = tmp
+                logger.info(f"{n_asymmetric_atoms=} {n_sym_atoms=}")
+
+    @staticmethod
+    def get_operators(self, value, periodicity=3):
+        """Get the symmetry operators as 4x4 matrices.
+
+        Parameters
+        ----------
+        value : str or [str]
+            The symmetry group or space group, or symmetry operators as strings.
+        periodicity : int
+            The periodicity of the structure. 3 for 3D, 2 for 2D, 1 for 1D, or 0.
+
+        Returns
+        -------
+        operators : np.ndarray
+            The symmetry operators as 4x4 matrices.
+        """
+        if isinstance(value, str):
+            if value == "":
+                raise RuntimeError("No group specified.")
+            if self.periodicity == 0:
+                if value != "C1":
+                    raise NotImplementedError("Point groups not implemented yet!")
+                W4 = [[1, 0, 0], [0, 1, 0], [0, 0, 1]]
+
+                W4s = []
+                W4s.append(W4)
+                operators = np.array(W4s, dtype=float)
+            else:
+                # Get the 4x4 augmented matrices
+                hall = _Symmetry.to_hall(value)
+                value = _Symmetry.hall_to_spacegroup_name(hall)
+                data = spglib.get_symmetry_from_database(hall)
+                W4s = []
+                for W, w in zip(
+                    data["rotations"].tolist(), data["translations"].tolist()
+                ):
+                    W4 = []
+                    for Wrow, w_i in zip(W, w):
+                        W4.append([*Wrow, w_i])
+                    W4.append([0, 0, 0, 1])
+                    W4s.append(W4)
+                operators = np.array(W4s, dtype=float)
+        elif isinstance(value, list):
+            W4s = []
+            operators = np.array(W4s, dtype=float)
+        return operators
+
+    @staticmethod
+    def hall_to_spacegroup_name(hall):
+        """Return the International name including setting for a hall number."""
+        if _Symmetry.hall_to_spgname is None:
+            _Symmetry.spgname_to_hall = None
+            _Symmetry.spacegroup_names_to_hall()
+        return _Symmetry.hall_to_spgname[hall]
+
+    @staticmethod
+    def spacegroup_names_to_hall():
+        """Dictionary of Hall number for spacegroup names."""
+        if _Symmetry.spgname_to_hall is None:
+            system_name = {
+                "international_full": "name_H-M_full",
+                "international": "name_H-M_alt",
+                "international_short": "name_H-M_short",
+                "hall_symbol": "name_Hall",
+            }
+            # Initialize the symmetry data
+            _Symmetry.spgname_to_hall = {}
+            _Symmetry.hall_to_spgname = {}
+            _Symmetry.hall_to_hall_symbol = {}
+            _Symmetry.hall_to_IT_number = {}
+            if _Symmetry.spgname_to_system is None:
+                _Symmetry.spgname_to_system = {}
+            for hall in range(1, 530):
+                data = spglib.get_spacegroup_type(hall)
+                # pws pprint.pprint(data)
+                choice = data["choice"]
+                _Symmetry.hall_to_hall_symbol[hall] = data["hall_symbol"]
+                _Symmetry.hall_to_IT_number[hall] = data["number"]
+
+                # Handle Hall to spacegroup using the full H-M name
+                key = "international_full"
+                name = data[key]
+
+                # Default setting if there are multiple, leave unadorned
+                if choice in ("2", "b", "b1", "H"):
+                    if (
+                        hall in _Symmetry.hall_to_spgname
+                        and name != _Symmetry.hall_to_spgname[hall]
+                    ):
+                        raise RuntimeError(
+                            f"{hall=} {key} --> {name} exists: "
+                            f"{_Symmetry.hall_to_spgname[hall]}"
+                        )
+                    else:
+                        _Symmetry.hall_to_spgname[hall] = name
+                else:
+                    # Add other settings to the spacegroup name
+                    if choice != "":
+                        name += f" :{choice}"
+                    if hall not in _Symmetry.hall_to_spgname:
+                        # pws if choice != "":
+                        # pws print(f"{hall} = {name}   QQQQQQQQQQQQQQQQQQQQQQ")
+                        _Symmetry.hall_to_spgname[hall] = name
+
+                # Now handle all the rest
+                for key in (
+                    "international_full",
+                    "international",
+                    "international_short",
+                    "hall_symbol",
+                ):
+                    name = data[key]
+
+                    # SPGLib encodes the international name: 'C 2/c = B 2/n 1 1',
+                    if key == "international":
+                        name = name.split("=")[0].strip()
+
+                    # Default setting if there are multiple, leave unadorned
+                    if choice in ("2", "b", "b1", "H"):
+                        _Symmetry.spgname_to_hall[name] = hall
+                        _Symmetry.spgname_to_system[name] = system_name[key]
+                        for txt in ("_", " "):
+                            tmp = name.replace(txt, "")
+                            _Symmetry.spgname_to_hall[tmp] = hall
+                            _Symmetry.spgname_to_system[tmp] = system_name[key]
+                            if tmp[-2:] == ":H":
+                                tmp = tmp[:-2].strip()
+                                _Symmetry.spgname_to_hall[tmp] = hall
+                                _Symmetry.spgname_to_system[tmp] = system_name[key]
+
+                    if key in ("international", "international_short") and choice != "":
+                        name += f" :{choice}"
+
+                    _Symmetry.spgname_to_hall[name] = hall
+                    _Symmetry.spgname_to_system[name] = system_name[key]
+                    for txt in ("_", " "):
+                        tmp = name.replace(txt, "")
+                        _Symmetry.spgname_to_hall[tmp] = hall
+                        _Symmetry.spgname_to_system[tmp] = system_name[key]
+                        if tmp[-2:] == ":H":
+                            tmp = tmp[:-2].strip()
+                            _Symmetry.spgname_to_hall[tmp] = hall
+                            _Symmetry.spgname_to_system[tmp] = system_name[key]
+        return _Symmetry.spgname_to_hall
+
+    @staticmethod
+    def spacegroup_numbers_to_hall():
+        """List of the Hall spacegroup names for the IT number."""
+        if _Symmetry.spgno_to_hall is None:
+            # Initialize the symmetry data
+            _Symmetry.spgno_to_hall = {}
+            if _Symmetry.spgname_to_system is None:
+                _Symmetry.spgname_to_system = {}
+            for hall in range(1, 530):
+                data = spglib.get_spacegroup_type(hall)
+                spgno = int(data["number"])
+                if spgno not in _Symmetry.spgno_to_hall:
+                    _Symmetry.spgno_to_hall[spgno] = hall
+                    _Symmetry.spgname_to_system[spgno] = "Int_Tables_number"
+                    _Symmetry.spgname_to_system[str(spgno)] = "Int_Tables_number"
+
+        return _Symmetry.spgno_to_hall
+
+    @staticmethod
+    def to_hall(name):
+        """Hall number given full spacegroup name or number."""
+        if isinstance(name, int):
+            return _Symmetry.spacegroup_numbers_to_hall()[name]
+        return _Symmetry.spacegroup_names_to_hall()[name]
+
+    @staticmethod
+    def spacegroup_names_to_system():
+        """Dictionary of system (name_Hall, name_H-M_full,...) for spacegroup names."""
+        if _Symmetry.spgname_to_hall is None:
+            _Symmetry.spacegroup_names_to_hall()
+        if _Symmetry.spgno_to_hall is None:
+            _Symmetry.spacegroup_numbers_to_hall()
+
+        return _Symmetry.spgname_to_system
+
     @property
     def configuration(self):
         """Return the configuration."""
@@ -330,193 +662,6 @@ def system_db(self):
         """Return the SystemDB object that contains this cell."""
         return self._system_db
 
-    @property
-    def spacegroup_numbers_to_hall(self):
-        """List of the Hall spacegroup names for the IT number."""
-        if _Symmetry.spgno_to_hall is None:
-            # Initialize the symmetry data
-            _Symmetry.spgno_to_hall = {}
-            if _Symmetry.spgname_to_system is None:
-                _Symmetry.spgname_to_system = {}
-            for hall in range(1, 530):
-                data = spglib.get_spacegroup_type(hall)
-                spgno = int(data["number"])
-                if spgno not in _Symmetry.spgno_to_hall:
-                    _Symmetry.spgno_to_hall[spgno] = hall
-                    _Symmetry.spgname_to_system[spgno] = "Int_Tables_number"
-                    _Symmetry.spgname_to_system[str(spgno)] = "Int_Tables_number"
-
-        return _Symmetry.spgno_to_hall
-
-    @property
-    def spacegroup_names_to_hall(self):
-        """Dictionary of Hall number for spacegroup names."""
-        if _Symmetry.spgname_to_hall is None:
-            system_name = {
-                "international_full": "name_H-M_full",
-                "international": "name_H-M_alt",
-                "international_short": "name_H-M_short",
-                "hall_symbol": "name_Hall",
-            }
-            # Initialize the symmetry data
-            _Symmetry.spgname_to_hall = {}
-            _Symmetry.hall_to_spgname = {}
-            _Symmetry.hall_to_hall_symbol = {}
-            _Symmetry.hall_to_IT_number = {}
-            if _Symmetry.spgname_to_system is None:
-                _Symmetry.spgname_to_system = {}
-            for hall in range(1, 530):
-                data = spglib.get_spacegroup_type(hall)
-                # pws pprint.pprint(data)
-                choice = data["choice"]
-                _Symmetry.hall_to_hall_symbol[hall] = data["hall_symbol"]
-                _Symmetry.hall_to_IT_number[hall] = data["number"]
-
-                # Handle Hall to spacegroup using the full H-M name
-                key = "international_full"
-                name = data[key]
-
-                # Default setting if there are multiple, leave unadorned
-                if choice in ("2", "b", "b1", "H"):
-                    if (
-                        hall in _Symmetry.hall_to_spgname
-                        and name != _Symmetry.hall_to_spgname[hall]
-                    ):
-                        raise RuntimeError(
-                            f"{hall=} {key} --> {name} exists: "
-                            f"{_Symmetry.hall_to_spgname[hall]}"
-                        )
-                    else:
-                        _Symmetry.hall_to_spgname[hall] = name
-                else:
-                    # Add other settings to the spacegroup name
-                    if choice != "":
-                        name += f" :{choice}"
-                    if hall not in _Symmetry.hall_to_spgname:
-                        # pws if choice != "":
-                        # pws print(f"{hall} = {name}   QQQQQQQQQQQQQQQQQQQQQQ")
-                        _Symmetry.hall_to_spgname[hall] = name
-
-                # Now handle all the rest
-                for key in (
-                    "international_full",
-                    "international",
-                    "international_short",
-                    "hall_symbol",
-                ):
-                    name = data[key]
-
-                    # SPGLib encodes the international name: 'C 2/c = B 2/n 1 1',
-                    if key == "international":
-                        name = name.split("=")[0].strip()
-
-                    # Default setting if there are multiple, leave unadorned
-                    if choice in ("2", "b", "b1", "H"):
-                        _Symmetry.spgname_to_hall[name] = hall
-                        _Symmetry.spgname_to_system[name] = system_name[key]
-                        for txt in ("_", " "):
-                            tmp = name.replace(txt, "")
-                            _Symmetry.spgname_to_hall[tmp] = hall
-                            _Symmetry.spgname_to_system[tmp] = system_name[key]
-                            if tmp[-2:] == ":H":
-                                tmp = tmp[:-2].strip()
-                                _Symmetry.spgname_to_hall[tmp] = hall
-                                _Symmetry.spgname_to_system[tmp] = system_name[key]
-
-                    if key in ("international", "international_short") and choice != "":
-                        name += f" :{choice}"
-
-                    _Symmetry.spgname_to_hall[name] = hall
-                    _Symmetry.spgname_to_system[name] = system_name[key]
-                    for txt in ("_", " "):
-                        tmp = name.replace(txt, "")
-                        _Symmetry.spgname_to_hall[tmp] = hall
-                        _Symmetry.spgname_to_system[tmp] = system_name[key]
-                        if tmp[-2:] == ":H":
-                            tmp = tmp[:-2].strip()
-                            _Symmetry.spgname_to_hall[tmp] = hall
-                            _Symmetry.spgname_to_system[tmp] = system_name[key]
-        return _Symmetry.spgname_to_hall
-
-    def old_spacegroup_names_to_hall(self):
-        """Dictionary of Hall number for spacegroup names."""
-        if _Symmetry.spgname_to_hall is None:
-            system_name = {
-                "international_full": "name_H-M_full",
-                "international": "name_H-M_alt",
-                "international_short": "name_H-M_short",
-                "hall_symbol": "name_Hall",
-            }
-            # Initialize the symmetry data
-            _Symmetry.spgname_to_hall = {}
-            _Symmetry.hall_to_spgname = {}
-            _Symmetry.hall_to_hall_symbol = {}
-            _Symmetry.hall_to_IT_number = {}
-            if _Symmetry.spgname_to_system is None:
-                _Symmetry.spgname_to_system = {}
-            for hall in range(1, 530):
-                data = spglib.get_spacegroup_type(hall)
-                choice = data["choice"]
-                _Symmetry.hall_to_hall_symbol[hall] = data["hall_symbol"]
-                _Symmetry.hall_to_IT_number[hall] = data["number"]
-                for key in (
-                    "international_full",
-                    "international",
-                    "international_short",
-                    "hall_symbol",
-                ):
-                    name = data[key]
-                    if "international" in key and choice in ("2",):
-                        if (
-                            name in _Symmetry.spgname_to_hall
-                            and hall != _Symmetry.spgname_to_hall[name]
-                        ):
-                            raise RuntimeError(
-                                f"{hall=} {key} --> {name} exists: "
-                                f"{_Symmetry.spgname_to_hall[name]}"
-                            )
-                        if key == "international":
-                            _Symmetry.hall_to_spgname[hall] = name
-                        _Symmetry.spgname_to_hall[name] = hall
-                        _Symmetry.spgname_to_system[name] = system_name[key]
-                        for txt in ("_", " "):
-                            tmp = name.replace(txt, "")
-                            _Symmetry.spgname_to_hall[tmp] = hall
-                            _Symmetry.spgname_to_system[tmp] = system_name[key]
-                    if choice != "":
-                        name += f" :{choice}"
-                    if (
-                        name in _Symmetry.spgname_to_hall
-                        and hall != _Symmetry.spgname_to_hall[name]
-                    ):
-                        raise RuntimeError(
-                            f"{hall=} {key} --> {name} exists: "
-                            f"{_Symmetry.spgname_to_hall[name]}"
-                        )
-                    if key == "international" and hall not in _Symmetry.spgname_to_hall:
-                        _Symmetry.hall_to_spgname[hall] = name
-                    _Symmetry.spgname_to_hall[name] = hall
-                    _Symmetry.spgname_to_system[name] = system_name[key]
-                    for txt in ("_", " "):
-                        tmp = name.replace(txt, "")
-                        _Symmetry.spgname_to_hall[tmp] = hall
-                        _Symmetry.spgname_to_system[tmp] = system_name[key]
-                        if tmp[-2:] == ":H":
-                            tmp = tmp[:-2].strip()
-                            _Symmetry.spgname_to_hall[tmp] = hall
-                            _Symmetry.spgname_to_system[tmp] = system_name[key]
-        return _Symmetry.spgname_to_hall
-
-    @property
-    def spacegroup_names_to_system(self):
-        """Dictionary of system (name_Hall, name_H-M_full,...) for spacegroup names."""
-        if _Symmetry.spgname_to_hall is None:
-            self.spacegroup_names_to_hall
-        if _Symmetry.spgno_to_hall is None:
-            self.spacegroup_numbers_to_hall
-
-        return _Symmetry.spgname_to_system
-
     def find_spacegroup_from_operators(self):
         """Find the spacegroup from the symmetry operators."""
         if self.configuration.periodicity > 0:
@@ -548,13 +693,6 @@ def find_spacegroup_from_operators(self):
                     hall_number += 1
             return self.hall_to_spacegroup_name(hall_number)
 
-    def hall_to_spacegroup_name(self, hall):
-        """Return the International name including setting for a hall number."""
-        if _Symmetry.hall_to_spgname is None:
-            _Symmetry.spgname_to_hall = None
-            self.spacegroup_names_to_hall
-        return _Symmetry.hall_to_spgname[hall]
-
     def reset_atoms(self):
         """The atoms have changed, so need to recalculate the symmetric atoms."""
         self._operators = None
@@ -592,15 +730,15 @@ def symmetrize_atomic_scalar(self, v_sym):
             return v_sym, None
 
         v_in = np.array(v_sym)
-        print(f"{v_in=}")
+        # print(f"{v_in=}")
         generators = self.atom_generators
         v = np.ndarray((len(generators),), dtype=float)
         delta = np.zeros_like(v_in)
         start = 0
         for asym_atom, ops in enumerate(generators):
-            print(f"{asym_atom=} {ops}")
+            # print(f"{asym_atom=} {ops}")
             n = len(ops)
-            print(f"{n=}")
+            # print(f"{n=}")
             v[asym_atom] = np.average(v_in[start : start + n], axis=0)
             delta[start : start + n] = v_in[start : start + n] - v[asym_atom]
             start += n
@@ -712,12 +850,6 @@ def symmetrize_coordinates(self, xyz_sym, fractionals=True):
                 tmp = self.configuration.cell.to_cartesians(uvw, as_array=True)
                 return tmp.round(8).tolist(), delta.round(8).tolist()
 
-    def to_hall(self, name):
-        """Hall number given full spacegroup name or number."""
-        if isinstance(name, int):
-            return self.spacegroup_numbers_to_hall[name]
-        return self.spacegroup_names_to_hall[name]
-
     def update_group(self, value):
         if self.configuration.periodicity == 0:
             if value != "C1":
@@ -739,10 +871,22 @@ def vector_x_symop(self, vector, symop, translation=True):
         return xformed[0:3]
 
     def _expand(self):
-        """Setup the information for going from asymmetric cell to full cell."""
+        """Setup the information for going from asymmetric cell to full cell.
+
+        There is an issue with some CIF files, notably those from the Cambridge
+        structural database, where there are atoms in the asymmetric unit that
+        are duplicates. It appears that the CIF file is written with entire molecules
+        in some cases. So if the molecule is on a symmetry element, the atoms that are
+        related by symmetry are explicitly in the file.
+
+        To handle this we need to check for duplicates and remove them. Not a bad check
+        to have, anyway.
+        """
         operators = self.symmetry_matrices
 
+        symbols = self.configuration.atoms.asymmetric_symbols
         n_atoms = self.configuration.n_asymmetric_atoms
+        logger.info(f"Expanding {n_atoms} asymmetric atoms to full cell.")
         if n_atoms == 0:
             self._atom_generators = None
             self._symop_to_atom = None
@@ -750,7 +894,7 @@ def _expand(self):
         self._symop_to_atom = []  # The symmetry atom resulting from symmetry ops
         if self.configuration.periodicity == 3:
             if len(operators) == 1:
-                # P1 is a special case. Nothing to do!
+                # P1 is a special case. Nothing to do.
                 self._atom_generators = [[0] for i in range(n_atoms)]
                 self._symop_to_atom = [[i] for i in range(n_atoms)]
                 self._atom_to_asymmetric_atom = [*range(n_atoms)]
@@ -763,34 +907,38 @@ def _expand(self):
                         f"Mismatch of number of atoms in symmetry: {uvw0.shape[0]} != "
                         f"{n_atoms}"
                     )
-                logger.debug("Original coordinates")
-                logger.debug(uvw0)
-                logger.debug(f"{uvw0.shape=}")
+                logger.debug(
+                    f"""
+{n_atoms=}
+{uvw0.shape=}"
+Original coordinates
+{uvw0}
+                    """
+                )
 
-                logger.debug(f"{n_atoms=}")
                 uvw = np.ndarray((n_atoms, 4))
-                logger.debug(f"{uvw.shape=}")
                 uvw[:, 0:3] = uvw0[:, :]
-                logger.debug("Expanded coordinates")
-                logger.debug(uvw)
+                logger.debug(f"\nExpanded coordinates\n{uvw}")
                 uvw[:, 3] = 1
                 # logger.debug(self.symops)
-                logger.debug(f"{operators.shape=}")
+                logger.debug(f"\n{operators.shape=}")
                 # logger.debug(operators)
-                logger.debug(f"{uvw.shape=}")
-                logger.debug("Expanded coordinates")
-                logger.debug(uvw)
+                # logger.debug(f"{uvw.shape=}")
+                # logger.debug("Expanded coordinates")
+                # logger.debug(uvw)
                 xformed = np.einsum("ijk,lk", operators, uvw)
-                logger.debug(f"{xformed.shape=}")
-                # logger.debug(xformed)
+                logger.debug(f"\n{xformed.shape=}\n{xformed}")
 
                 # For comparison, bring all atoms into the cell [0..1)
                 tmp = xformed - np.floor(xformed)
 
-                logger.debug("tmp")
-                logger.debug(tmp)
+                logger.debug(f"\ntmp = Transformed and floored coordinates\n{tmp}")
+
+                # Keep track of atoms and coordinates found so can check for duplicates
+                found = {}
 
                 n_sym_atoms = 0
+                n_asymmetric_atoms = 0
                 for i in range(n_atoms):
                     values, I1, I2 = np.unique(
                         np.round(tmp[:, :3, i], 4),
@@ -798,15 +946,36 @@ def _expand(self):
                         return_index=True,
                         return_inverse=True,
                     )
-                    logger.debug(f"{i=}")
                     # pprint.pprint(np.round(tmp[:, :3, i], 4).tolist())
-                    logger.debug("values")
-                    # pprint.pprint(values.tolist())
-                    logger.debug(values)
-                    logger.debug("I1")
-                    logger.debug(I1)
-                    logger.debug("I2")
-                    logger.debug(I2)
+                    logger.debug(
+                        f"""
+{i=}: {symbols[i]}
+{np.round(tmp[:, :3, i], 4)}
+nvalues
+{values}
+I1
+{I1}
+I2
+{I2}
+                        """
+                    )
+                    # Check for duplicates
+                    if tuple(values[0]) in found:
+                        if found[tuple(values[0])] != symbols[i]:
+                            raise RuntimeError(
+                                "Duplicate atoms with different symbols: "
+                                f"{symbols[i]} != {found[tuple(values[0])]} at "
+                                f"{values[0]}"
+                            )
+                        else:
+                            # Same element, so ignore
+                            logger.debug("\nDuplicate atom, ignoring")
+                            continue
+
+                    for value in values:
+                        found[tuple(value)] = symbols[i]
+
+                    n_asymmetric_atoms += 1
                     I2 += n_sym_atoms
                     self._atom_generators.append(I1.tolist())
                     self._symop_to_atom.append(I2.tolist())
@@ -815,12 +984,7 @@ def _expand(self):
                 for i, generators in enumerate(self._atom_generators):
                     tmp.extend([i] * len(generators))
                 self._atom_to_asymmetric_atom = tmp
-        else:
-            if len(operators) == 1:
-                # P1 is a special case. Nothing to do!
-                self._atom_generators = [[0] for i in range(n_atoms)]
-                self._symop_to_atom = [[i] for i in range(n_atoms)]
-                self._atom_to_asymmetric_atom = [*range(n_atoms)]
+                logger.info(f"{n_asymmetric_atoms=} {n_sym_atoms=}")
 
     def _expand_bonds(self):
         """Expand the list of asymmetric bonds to the full list.
diff --git a/molsystem/system_db.py b/molsystem/system_db.py
index d50a1f1..3a64e2c 100644
--- a/molsystem/system_db.py
+++ b/molsystem/system_db.py
@@ -970,6 +970,21 @@ def _initialize(self):
                 "    ON velocities(atom, configuration)"
             )
 
+            # Keep the gradients separate from coordinates to keep the size down
+            table = self["gradients"]
+            table.add_attribute(
+                "configuration", coltype="int", references="configuration"
+            )
+            table.add_attribute("atom", coltype="int", references="atom")
+            table.add_attribute("gx", coltype="float")
+            table.add_attribute("gy", coltype="float")
+            table.add_attribute("gz", coltype="float")
+            self.db.execute("CREATE INDEX 'idx_gradients_atom' ON gradients (\"atom\")")
+            self.db.execute(
+                "CREATE INDEX idx_gradients_atom_configuration "
+                "    ON gradients(atom, configuration)"
+            )
+
             # The definition of the subsets -- templates
             table = self["template"]
             table.add_attribute("id", coltype="int", pk=True)
diff --git a/tests/conftest.py b/tests/conftest.py
index a18709c..2c88639 100644
--- a/tests/conftest.py
+++ b/tests/conftest.py
@@ -322,6 +322,9 @@ def vanadium(configuration):
         x=[0.0, 0.5],
         y=[0.0, 0.5],
         z=[0.0, 0.5],
+        gx=[0.0, 0.02],
+        gy=[0.01, -0.02],
+        gz=[-0.01, 0.01],
         vx=[0.0, 0.02],
         vy=[0.01, -0.02],
         vz=[-0.01, 0.01],
diff --git a/tests/test_atoms.py b/tests/test_atoms.py
index 03470f1..4c06ec5 100644
--- a/tests/test_atoms.py
+++ b/tests/test_atoms.py
@@ -23,7 +23,20 @@ def test_construction(configuration):
 
 def test_keys(atoms):
     """Test the default keys in an Atoms object"""
-    result = ["atno", "configuration", "id", "vx", "vy", "vz", "x", "y", "z"]
+    result = [
+        "atno",
+        "configuration",
+        "gx",
+        "gy",
+        "gz",
+        "id",
+        "vx",
+        "vy",
+        "vz",
+        "x",
+        "y",
+        "z",
+    ]
 
     if sorted([*atoms.keys()]) != result:
         print(sorted([*atoms.keys()]))
@@ -80,7 +93,21 @@ def test_append_error(atoms):
 
 def test_add_attribute(atoms):
     """Test adding an attribute"""
-    result = ["atno", "configuration", "id", "name", "vx", "vy", "vz", "x", "y", "z"]
+    result = [
+        "atno",
+        "configuration",
+        "gx",
+        "gy",
+        "gz",
+        "id",
+        "name",
+        "vx",
+        "vy",
+        "vz",
+        "x",
+        "y",
+        "z",
+    ]
     with atoms as tmp:
         tmp.add_attribute("name")
     assert sorted([*atoms.keys()]) == result
@@ -88,7 +115,21 @@ def test_add_attribute(atoms):
 
 def test_add_duplicate_attribute(atoms):
     """Test duplicate adding an attribute"""
-    result = ["atno", "configuration", "id", "name", "vx", "vy", "vz", "x", "y", "z"]
+    result = [
+        "atno",
+        "configuration",
+        "gx",
+        "gy",
+        "gz",
+        "id",
+        "name",
+        "vx",
+        "vy",
+        "vz",
+        "x",
+        "y",
+        "z",
+    ]
     with atoms as tmp:
         tmp.add_attribute("name")
     with pytest.raises(RuntimeError) as e:
@@ -103,13 +144,34 @@ def test_add_duplicate_attribute(atoms):
 
 def test_add_coordinates_attribute(atoms):
     """Test adding an attribute"""
-    result = ["atno", "configuration", "id", "spin", "vx", "vy", "vz", "x", "y", "z"]
+    correct = [
+        "atno",
+        "configuration",
+        "gx",
+        "gy",
+        "gz",
+        "id",
+        "spin",
+        "vx",
+        "vy",
+        "vz",
+        "x",
+        "y",
+        "z",
+    ]
     with atoms as tmp:
         tmp.add_attribute("spin", configuration_dependent=True)
-    assert sorted([*atoms.keys()]) == result
+    result = sorted([*atoms.keys()])
+    if result != correct:
+        print(result)
+    assert result == correct
+
     del atoms["spin"]
-    result.remove("spin")
-    assert sorted([*atoms.keys()]) == result
+    correct.remove("spin")
+    result = sorted([*atoms.keys()])
+    if result != correct:
+        print(result)
+    assert result == correct
 
 
 def test_add_attribute_with_values(atoms):
@@ -186,6 +248,9 @@ def test_deleting_column(atoms):
         del tmp["atno"]
     assert sorted([*atoms.keys()]) == [
         "configuration",
+        "gx",
+        "gy",
+        "gz",
         "id",
         "vx",
         "vy",
@@ -519,3 +584,52 @@ def test_periodic_velocities_cartesians(vanadium):
 
     vxyz = configuration.atoms.get_velocities(fractionals=False)
     assert np.allclose(vxyz, vxyz0)
+
+
+def test_have_gradients(AceticAcid):
+    """Test whether there are gradients."""
+    assert not AceticAcid.atoms.have_gradients
+
+
+def test_set_gradients(AceticAcid):
+    """Test setting gradients."""
+    configuration = AceticAcid
+
+    xs = [1.08, 0.58, 0.72, 0.71, 0.57, -0.13, 0.98, 2.17]
+    ys = [0.02, 3.14, -0.67, -0.31, 1.39, 1.71, 2.30, 0.02]
+    zs = [-0.02, 0.28, -0.79, 0.95, -0.32, -1.26, 0.59, -0.03]
+    gxyz0 = [[x, y, z] for x, y, z in zip(xs, ys, zs)]
+
+    configuration.atoms.set_gradients(gxyz0)
+
+    gxyz = configuration.atoms.gradients
+    if not np.allclose(gxyz, gxyz0):
+        print(gxyz)
+        print(" != ", gxyz0)
+    assert np.allclose(gxyz, gxyz0)
+
+
+def test_periodic_gradients(vanadium):
+    """Test getting gradients."""
+    configuration = vanadium
+
+    gxs = [0.0, 0.02]
+    gys = [0.01, -0.02]
+    gzs = [-0.01, 0.01]
+    gxyz0 = [[gx, gy, gz] for gx, gy, gz in zip(gxs, gys, gzs)]
+
+    gxyz = configuration.atoms.gradients
+    assert np.allclose(gxyz, gxyz0)
+
+
+def test_periodic_gradients_cartesians(vanadium):
+    """Test getting gradients in Cartesian gradients."""
+    configuration = vanadium
+
+    gxs = [0.0, 0.02]
+    gys = [0.01, -0.02]
+    gzs = [-0.01, 0.01]
+    gxyz0 = [[gx * 3.03, gy * 3.03, gz * 3.03] for gx, gy, gz in zip(gxs, gys, gzs)]
+
+    gxyz = configuration.atoms.get_gradients(fractionals=False)
+    assert np.allclose(gxyz, gxyz0)
diff --git a/tests/test_symmetry.py b/tests/test_symmetry.py
index 90784f6..73ec1b2 100644
--- a/tests/test_symmetry.py
+++ b/tests/test_symmetry.py
@@ -24,7 +24,7 @@ def test_spacegroup_names(symmetry):
     path = Path(__file__).parent / "data" / "spacegroup_names.json"
     with path.open() as fd:
         correct = json.load(fd)
-    answer = symmetry.spacegroup_names_to_hall
+    answer = symmetry.spacegroup_names_to_hall()
     if answer != correct:
         print(json.dumps(answer, indent=3))
     assert answer == correct
diff --git a/tests/test_templates.py b/tests/test_templates.py
index 00c8e51..916dd88 100644
--- a/tests/test_templates.py
+++ b/tests/test_templates.py
@@ -177,7 +177,7 @@ def test_atoms_data(gly):
         ],
     }
     for column in gly.atoms:
-        if column in ("vx", "vy", "vz"):
+        if column in ("gx", "gy", "gz", "vx", "vy", "vz"):
             continue
         result = gly.atoms.get_column_data(column)
         if column not in answer: