diff --git a/scripts/match.py b/scripts/match.py
new file mode 100644
index 00000000..cf1f6c5c
--- /dev/null
+++ b/scripts/match.py
@@ -0,0 +1,127 @@
+
+import pymatgen as mg
+import numpy as np
+from pyxtal.util import parse_cif
+from optparse import OptionParser
+import pymatgen.analysis.structure_matcher as sm
+from pyxtal import pyxtal
+from pyxtal.XRD import Similarity
+from pyxtal.optimize.base import GlobalOptimize
+import warnings
+warnings.filterwarnings("ignore")
+
+parser = OptionParser()
+parser.add_option("-f", dest="cif", default="WFS-gaff.cif", help="input cif file")
+parser.add_option("-r", dest="ref", help="reference")
+parser.add_option("-o", dest="out", default="Matched.cif", help="output")
+parser.add_option("--emin", dest="emin", type=float, default=0,
+                  help="minimum energy, default 0")
+parser.add_option("--emax", dest="emax", type=float, default=100,
+                  help="maximum energy, default 100")
+parser.add_option("--early_stop", dest="early", action="store_true", default=False,
+                  help="stop when the first match is found")
+parser.add_option("--XRD", dest="xrd",
+                  action="store_true", default=False,
+                  help="Compare XRD")
+parser.add_option("--smin", dest="smin", type=float, default=0.80,
+                  help="min similarity for XRD")
+
+(options, args) = parser.parse_args()
+matcher = sm.StructureMatcher(ltol=0.3, stol=0.4, angle_tol=5.0)
+
+with open(options.cif, 'r') as f:
+    lines = f.readlines()
+    smiles = []
+    for l in lines:
+        if 'smile' in l:
+            smile_str = l.split(':')[1].strip()
+            smiles = [s + '.smi' for s in smile_str.split('.')]
+            break
+print(smiles)
+with open(options.out, 'w') as f: f.write(f'smiles: {smile_str}\n')
+
+if options.ref is None:
+    raise ValueError("Reference structure is required.")
+else:
+    pmg_ref = mg.core.Structure.from_file(options.ref)
+    xtal = pyxtal(molecular=True)
+    xtal.from_seed(pmg_ref, molecules = smiles)
+    print(f"Reference Structure loaded from {options.ref} {pmg_ref.density:.3f}")
+    pmg_ref.remove_species("H")
+    print(xtal)
+
+    if options.xrd:
+        thetas = [0, 35.0]
+        xrd = xtal.get_XRD(thetas=thetas)
+        p_ref = xrd.get_profile(res=0.15, user_kwargs={"FWHM": 0.25})
+
+cifs, engs = parse_cif(options.cif, eng=True)
+print("Total Number of Structures:", len(cifs))
+engs = np.array(engs)
+ids = np.argsort(engs)
+cifs = [cifs[id] for id in ids]
+engs = engs[ids]
+print(f"Min energy in eV: {engs.min()+options.emin/96.485:.4f} {engs.min()+options.emax/96.485:.4f}")
+engs_norm = engs - engs.min()  # Normalize energies to the lowest one
+engs_norm *= 96.485
+
+# Find the id of energy that is between [options.emin, options.emax]
+n1 = np.searchsorted(engs_norm, options.emin, side='left')
+n2 = np.searchsorted(engs_norm, options.emax, side='right')
+engs = engs[n1:n2]
+engs_norm = engs_norm[n1:n2]
+cifs = [cifs[id] for id in range(n1, n2)]
+ids = ids[n1:n2]
+count = 0
+xtal = pyxtal(molecular=True)
+for id, cif in enumerate(cifs):
+    pmg = mg.core.Structure.from_str(cif, fmt='cif')
+    try:
+        xtal.from_seed(pmg, molecules = smiles)
+        #xtal.energy = engs_norm[id]
+        spg = xtal.group.number
+        den = xtal.get_density()
+        raw_eng = engs[id]
+        norm_eng = engs_norm[id]
+        match = False
+        strs = f"Struc {ids[id]:6d}: {spg:3d} {raw_eng:.3f} kJ/mol, {den:.3f} g/cm^3, {norm_eng:.3f}"
+        if options.xrd:
+            p1 = xtal.get_XRD(thetas=thetas).get_profile(res=0.15, user_kwargs={"FWHM": 0.25})
+            sim = Similarity(p1, p_ref, x_range=thetas).value
+            if sim > options.smin: match = True
+            strs += f'{sim:12.3f} in PXRD similarity'
+        else:
+            pmg.remove_species("H")
+            if abs(pmg.density-pmg_ref.density) <= 0.35:
+                strs += '****'
+                if matcher.fit(pmg, pmg_ref):
+                    match = True
+                    
+                    pmg_gen = xtal.to_pymatgen()
+                    pmg_gen.remove_species("H")
+
+                    ref_copy = pmg_ref.copy()
+                    ref_copy.remove_species("H")
+                    print(f"Generated structure sites: {len(pmg_gen.sites)}, Reference sites: {len(ref_copy.sites)}")
+                    print(f"Species (gen): {[str(sp) for sp in pmg_gen.species]}")
+                    print(f"Species (ref): {[str(sp) for sp in ref_copy.species]}")
+                    rmsd = matcher.get_rms_dist( ref_copy, pmg_gen)
+                    if rmsd is not None:
+                        rms_lat, rms_cart = rmsd
+                        print(f"RMSD – Lattice: {rms_lat:.3f} Å, Cartesian: {rms_cart:.3f} Å")
+                    else:
+                        print("RMSD calculation failed: structures not comparable")
+                    
+    except:
+        continue
+
+    if match:
+        count += 1
+        strs += '+++++++++++'
+        label = f"{count}-d{den:.3f}-spg{spg}-e{norm_eng:.3f}"
+        if options.xrd: label += f"-s{sim:.3f}"
+        with open(options.out, 'a+') as f: f.writelines(xtal.to_file(header=label))
+        if options.early: break
+    print(strs)
+print(f"Found {count} matches")
+
diff --git a/scripts/rank.py b/scripts/rank.py
new file mode 100644
index 00000000..55a21050
--- /dev/null
+++ b/scripts/rank.py
@@ -0,0 +1,148 @@
+import pymatgen as mg
+import numpy as np
+from pyxtal import pyxtal
+from pyxtal.util import parse_cif
+from optparse import OptionParser
+import pymatgen.analysis.structure_matcher as sm
+
+import warnings
+warnings.filterwarnings("ignore")
+
+def new_struc(xtal, xtals, max_num=100):
+    """
+    check if this is a new structure
+
+    Args:
+        xtal: input structure
+        xtals: list of reference structures
+
+    Return:
+        `None` or the id of matched structure
+    """
+    if len(xtals) > max_num:
+        start = len(xtals) - max_num
+    else:
+        start = 0
+    sg1 = xtal.group.number
+    pmg_s1 = xtal.to_pymatgen()
+    pmg_s1.remove_species("H")
+    vol1 = pmg_s1.lattice.volume
+
+    for xtal2 in xtals[start:]:
+        sg2 = xtal2.group.number
+        if sg1 == sg2:
+            pmg_s2 = xtal2.to_pymatgen()
+            vol2 = pmg_s2.lattice.volume
+            if abs(vol1-vol2)/vol1<5e-2:
+                pmg_s2.remove_species("H")
+                if sm.StructureMatcher().fit(pmg_s1, pmg_s2):
+                    return False
+    return True
+
+
+parser = OptionParser()
+parser.add_option("-f", "--cif", dest="cif", default="WFS-gaff.cif",
+                  help="cif file name, optional")
+parser.add_option("-r", "--rank", dest="rank", default='energy',
+                  help="ranking criteria: default is energy")
+parser.add_option("-s", "--n1", dest="n1", type=int, default=0,
+                  help="starting id, optional")
+parser.add_option("-e", "--n2", dest="n2", type=int, default=-1,
+                  help="ennding id, optional")
+parser.add_option("-c", "--cut", dest="cut", type=int,
+                  help="cutoff number, optional")
+parser.add_option("--dmax", dest="dmax", type=float, default=10.0,
+                  help="maximum density in g/cm^3, optional")
+
+(options, args) = parser.parse_args()
+rank = options.rank
+n1 = options.n1
+n2 = options.n2
+output1 = 'Ranked.cif'
+
+"""
+Read the smile from the following contents
+-------Global Crystal Structure Prediction------
+smile     : CC(=O)OC1=CC=CC=C1C(=O)O
+"""
+
+with open(options.cif, 'r') as f:
+    lines = f.readlines()
+    smiles = []
+    for l in lines:
+        if 'smile' in l:
+            smile_str = l.split(':')[1].strip()
+            smiles = [smile_str + '.smi']
+            break
+print(smiles)
+
+cifs, engs = parse_cif(options.cif, eng=True)
+print("Total Number of Structures:", len(cifs))
+if options.cut is None:
+    cut = len(cifs)
+else:
+    cut = options.cut
+
+if options.rank == 'energy':
+    engs = np.array(engs)
+    ids = np.argsort(engs)
+else:
+    cifs, sims = parse_cif(options.cif, sim=True)
+    sims = np.array(sims)
+    ids = np.argsort(-1*sims)
+    sims = [sims[id] for id in ids]
+
+cifs = [cifs[id] for id in ids]
+engs = engs[ids]
+eng0 = engs.min()
+
+if n2 == -1:
+    cifs = cifs[n1:]
+    engs = engs[n1:]
+    ids = ids[n1:]
+else:
+    n2 = min(n2, len(cifs))
+    cifs = cifs[n1:n2]
+    engs = engs[n1:n2]
+    ids = ids[n1:n2]
+print("Index", n1, n2, cut, len(cifs))
+with open(output1, 'w') as f: f.write(l)
+
+xtals = []
+count = 0
+with open(output1, 'a+') as f:
+    for id, cif in enumerate(cifs):
+        pmg = mg.core.Structure.from_str(cif, fmt='cif')
+        try:
+            xtal = pyxtal(molecular=True)
+            xtal.from_seed(pmg, molecules = smiles)
+            xtal.energy = engs[id]
+            if new_struc(xtal, xtals, 100):
+                xtals.append(xtal)
+                spg = xtal.group.number
+                den = xtal.get_density()
+                eng = engs[id]
+                label = f"{count}-d{den:.3f}-spg{spg}-e{eng:.3f}"
+                if den < options.dmax:
+                    cif_lines = xtal.to_file(header=label).splitlines()
+                    # Extract energy from label
+                    energy_from_label = float(label.split('-e')[1])
+                    energy_line = f"#Energy: {energy_from_label:.3f} eV/cell\n"
+                    output_lines = []
+                    energy_replaced = False
+                    for line in cif_lines:
+                        stripped = line.strip()
+                        # Replace any existing energy line
+                        if stripped.startswith("#Energy:"):
+                            output_lines.append(energy_line)
+                            energy_replaced = True
+                            continue
+                        output_lines.append(line + '\n')
+                    f.writelines(output_lines)
+                    print(f"{ids[id]:6d} {label} {(eng - eng0)*96.485:6.2f}")
+                count += 1
+                if count == cut:
+                    print("Stop", count)
+                    break
+        except:
+            print("Problem in reading")