diff --git a/pysisyphus/Geometry.py b/pysisyphus/Geometry.py
index 79eff2daa..b2b8ef09a 100644
--- a/pysisyphus/Geometry.py
+++ b/pysisyphus/Geometry.py
@@ -177,6 +177,8 @@ def __init__(
         coord_kwargs=None,
         isotopes=None,
         freeze_atoms=None,
+        remove_com=False,
+        remove_centroid=False,
         comment="",
         name="",
     ):
@@ -208,6 +210,10 @@ def __init__(
             Given a float, this float will be directly used as mass.
         freeze_atoms : iterable of integers
             Specifies which atoms should remain fixed at their initial positions.
+        remove_com : bool, optional
+            Move center of mass to the origin.
+        remove_centroid : bool, optional
+            Move centroid to the origin.
         comment : str, optional
             Comment string.
         name : str, optional
@@ -236,10 +242,17 @@ def __init__(
         elif type(freeze_atoms) is str:
             freeze_atoms = full_expand(freeze_atoms)
         self.freeze_atoms = np.array(freeze_atoms, dtype=int)
+        self.remove_com = bool(remove_com)
+        self.remove_centroid = bool(remove_centroid)
         self.comment = comment
         self.name = name
 
         self._masses = None
+        if self.remove_com:
+            self.coords3d = self.coords3d - self.center_of_mass[None, :]
+        if self.remove_centroid:
+            self.coords3d = self.coords3d - self.centroid[None, :]
+
         self._energy = None
         self._forces = None
         self._hessian = None
@@ -551,7 +564,7 @@ def set_coords(self, coords, cartesian=False, update_constraints=False):
         coords = np.array(coords).flatten()
 
         # Do Internal->Cartesian backtransformation if internal coordinates are used.
-        if self.internal:
+        if hasattr(self, "internal") and self.internal:
             # When internal coordinates are employed it may happen, that the underlying
             # Cartesian coordinates are updated, e.g. from the IPIServer calculator, which
             # may yield different internal coordinates.
diff --git a/scripts/modes3d.py b/scripts/modes3d.py
index 2802ddaf1..4cf171259 100755
--- a/scripts/modes3d.py
+++ b/scripts/modes3d.py
@@ -42,6 +42,9 @@
     "o": 1.25,
     "s": 1.25,
 }
+GEOM_KWARGS = {
+    "remove_com": True,
+}
 
 
 def get_bond_data(coords3d, bonds):
@@ -134,7 +137,7 @@ def palindrome_cycler(iterable):
 
 
 def from_h5_hessian(fn):
-    geom = geom_from_hessian(fn, coord_type="redund")
+    geom = geom_from_hessian(fn, coord_type="redund", **GEOM_KWARGS)
     geoms = (geom,)
 
     # Cartesian displacements
@@ -158,7 +161,7 @@ def from_h5_hessian(fn):
 
 
 def from_geom(fn):
-    geoms = geom_loader(fn, coord_type="redund", iterable=True)
+    geoms = geom_loader(fn, coord_type="redund", iterable=True, **GEOM_KWARGS)
     trj_coords = [geom.coords3d for geom in geoms]
     cycler = palindrome_cycler(trj_coords)
     return geoms, cycler, trj_coords