diff --git a/pyxtal/__init__.py b/pyxtal/__init__.py
index 17248535..4b5ce2dc 100644
--- a/pyxtal/__init__.py
+++ b/pyxtal/__init__.py
@@ -1386,11 +1386,22 @@ def optimize_lattice(self, iterations=5, force=False, standard=False):
                 msg = "Cannot find the standard setting"
                 raise RuntimeError(msg)
 
+    def update_wyckoffs(self):
+        """
+        rescale the coordinates in the wyckoff site
+        """
+
+        sites = self.mol_sites if self.molecular else self.atom_sites
+
+        for site in sites:
+            site.update()
+
     def get_std_representation(self, trans):
         """
         Perform cell transformation so that the symmetry operations
         follow standard space group notation
         """
+        pass
 
     def get_1D_representation(self, standard=False):
         """
diff --git a/pyxtal/interface/charmm.py b/pyxtal/interface/charmm.py
index d761f08d..d416f0f2 100644
--- a/pyxtal/interface/charmm.py
+++ b/pyxtal/interface/charmm.py
@@ -179,6 +179,12 @@ def write(self):
                             j + 1 + count, i + 1, res_name, label, *coord
                         )
                     )
+                    # quickly check if 
+                    if abs(coord).max() > 500.0:
+                        print("Unexpectedly large input coordinates, stop and debug")
+                        print(self.structure)
+                        self.structure.to_file('bug.cif')
+                        import sys; sys.exit()
 
             f.write(f"write psf card name {self.psf:s}\n")
             f.write(f"write coor crd card name {self.crd:s}\n")
@@ -282,6 +288,7 @@ def read(self):
                     count += len(site.molecule.mol)
                 # print("after relaxation  : ", self.structure.lattice, "iter: ", self.structure.iter)
                 self.structure.optimize_lattice()
+                self.structure.update_wyckoffs()
                 # print("after latticeopt  : ", self.structure.lattice, self.structure.check_distance()); import sys; sys.exit()
             except:
                 # molecular connectivity or lattice optimization
diff --git a/pyxtal/wyckoff_site.py b/pyxtal/wyckoff_site.py
index 05d277c3..fd55ca2b 100644
--- a/pyxtal/wyckoff_site.py
+++ b/pyxtal/wyckoff_site.py
@@ -750,13 +750,18 @@ def show_molecule_in_box(self, id=0):
         cell, vertices, center = self.molecule.get_box_coordinates(mol.cart_coords)
         return display_molecule(mol, center, cell)
 
-    def update(self, coords, lattice=None, absolute=False, update_mol=True):
+    def update(self, coords=None, lattice=None, absolute=False, update_mol=True):
         """
         After the geometry relaxation, the returned atomic coordinates
         maybe rescaled to [0, 1] bound. In this case, we need to refind
         the molecular coordinates according to the original neighbor list.
         If the list does not change, we return the new coordinates
         otherwise, terminate the calculation.
+
+        Args:
+            coords (float): input xyz coordinates for the given molecule
+            absolute (bool): whether the input is absolute coordinates
+            update_mol (bool): whether update the molecule (used for substitution)
         """
         from pyxtal.molecule import Orientation, compare_mol_connectivity
 
@@ -768,13 +773,20 @@ def update(self, coords, lattice=None, absolute=False, update_mol=True):
 
         if lattice is not None:
             self.lattice = lattice
+
+        if coords is None:
+            coords, _ = self._get_coords_and_species(absolute=False, first=True, unitcell=True)
+            absolute = False
+
         if not absolute:
             coords = coords.dot(self.lattice.matrix)
+
         # mol = Molecule(self.symbols, coords-np.mean(coords, axis=0))
         center = self.molecule.get_center(coords)
         mol = Molecule(self.symbols, coords - center)
 
         # match, _ = compare_mol_connectivity(mol, self.mol, True)
+        # Update the orientation matrix
         match, _ = compare_mol_connectivity(mol, self.molecule.mol)
         if match:
             # position = np.mean(coords, axis=0).dot(self.lattice.inv_matrix)