diff --git a/src/openqdc/__init__.py b/src/openqdc/__init__.py
index 1432923..87fc457 100644
--- a/src/openqdc/__init__.py
+++ b/src/openqdc/__init__.py
@@ -7,9 +7,32 @@
 
 # Dictionary of objects to lazily import; maps the object's name to its module path
 
-_lazy_imports_obj = {}
+_lazy_imports_obj = {
+    "ANI1": "openqdc.datasets.ani",
+    "ANI1CCX": "openqdc.datasets.ani",
+    "ANI1X": "openqdc.datasets.ani",
+    "Spice": "openqdc.datasets.spice",
+    "GEOM": "openqdc.datasets.geom",
+    "QMugs": "openqdc.datasets.qmugs",
+    "ISO17": "openqdc.datasets.iso_17",
+    "COMP6": "openqdc.datasets.comp6",
+    "GDML": "openqdc.datasets.gdml",
+    "Molecule3D": "openqdc.datasets.molecule3d",
+    "OrbnetDenali": "openqdc.datasets.orbnet_denali",
+    "SN2RXN": "openqdc.datasets.sn2_rxn",
+    "QM7X": "openqdc.datasets.qm7x",
+    "DESS": "openqdc.datasets.dess",
+    "NablaDFT": "openqdc.datasets.nabladft",
+    "SolvatedPeptides": "openqdc.datasets.solvated_peptides",
+    "WaterClusters": "openqdc.datasets.waterclusters3_30",
+    "TMQM": "openqdc.datasets.tmqm",
+    "Dummy": "openqdc.datasets.dummy",
+    "PCQM_B3LYP": "openqdc.datasets.pcqm",
+    "PCQM_PM6": "openqdc.datasets.pcqm",
+    "Transition1X": "openqdc.datasets.transition1x",
+}
 
-_lazy_imports_mod = {"datasets": "openqdc.datamodule", "utils": "openqdc.utils"}
+_lazy_imports_mod = {"datasets": "openqdc.datasets", "utils": "openqdc.utils"}
 
 
 def __getattr__(name):
diff --git a/src/openqdc/datasets/ani.py b/src/openqdc/datasets/ani.py
index 913fb8a..3f1b92b 100644
--- a/src/openqdc/datasets/ani.py
+++ b/src/openqdc/datasets/ani.py
@@ -39,6 +39,12 @@ class ANI1(BaseDataset):
     def root(self):
         return p_join(get_local_cache(), "ani")
 
+    def __smiles_converter__(self, x):
+        """util function to convert string to smiles: useful if the smiles is
+        encoded in a different format than its display format
+        """
+        return "-".join(x.decode("ascii").split("-")[:-1])
+
     @property
     def preprocess_path(self):
         path = p_join(self.root, "preprocessed", self.__name__)
@@ -89,6 +95,12 @@ class ANI1CCX(ANI1):
     __force_methods__ = []
     force_target_names = []
 
+    def __smiles_converter__(self, x):
+        """util function to convert string to smiles: useful if the smiles is
+        encoded in a different format than its display format
+        """
+        return x
+
 
 class ANI1X(ANI1):
     """
@@ -145,3 +157,9 @@ class ANI1X(ANI1):
 
     def convert_forces(self, x):
         return super().convert_forces(x) * 0.529177249  # correct the Dataset error
+
+    def __smiles_converter__(self, x):
+        """util function to convert string to smiles: useful if the smiles is
+        encoded in a different format than its display format
+        """
+        return x
diff --git a/src/openqdc/datasets/base.py b/src/openqdc/datasets/base.py
index 1de6ff1..d7c8234 100644
--- a/src/openqdc/datasets/base.py
+++ b/src/openqdc/datasets/base.py
@@ -1,4 +1,6 @@
 import os
+import pickle as pkl
+from copy import deepcopy
 from os.path import join as p_join
 from typing import Dict, List, Optional, Union
 
@@ -14,17 +16,27 @@
     IsolatedAtomEnergyFactory,
     chemical_symbols,
 )
-from openqdc.utils.constants import NB_ATOMIC_FEATURES
+from openqdc.utils.constants import (
+    NB_ATOMIC_FEATURES,
+    NOT_DEFINED,
+    POSSIBLE_NORMALIZATION,
+)
+from openqdc.utils.exceptions import (
+    DatasetNotAvailableError,
+    NormalizationNotAvailableError,
+    StatisticsNotAvailableError,
+)
 from openqdc.utils.io import (
     copy_exists,
     dict_to_atoms,
     get_local_cache,
     load_hdf5_file,
+    load_pkl,
     pull_locally,
     push_remote,
     set_cache_dir,
 )
-from openqdc.utils.molecule import atom_table
+from openqdc.utils.molecule import atom_table, z_to_formula
 from openqdc.utils.package_utils import requires_package
 from openqdc.utils.units import get_conversion
 
@@ -43,7 +55,7 @@ def extract_entry(
 
     res = dict(
         name=np.array([df["name"][i]]),
-        subset=np.array([subset]),
+        subset=np.array([subset if subset is not None else z_to_formula(x)]),
         energies=energies.reshape((1, -1)).astype(np.float32),
         atomic_inputs=np.concatenate((xs, positions), axis=-1, dtype=np.float32),
         n_atoms=np.array([x.shape[0]], dtype=np.int32),
@@ -64,8 +76,8 @@ def read_qc_archive_h5(
 ) -> List[Dict[str, np.ndarray]]:
     data = load_hdf5_file(raw_path)
     data_t = {k2: data[k1][k2][:] for k1 in data.keys() for k2 in data[k1].keys()}
-    n = len(data_t["molecule_id"])
 
+    n = len(data_t["molecule_id"])
     samples = [extract_entry(data_t, i, subset, energy_target_names, force_target_names) for i in tqdm(range(n))]
     return samples
 
@@ -83,6 +95,8 @@ class BaseDataset(torch.utils.data.Dataset):
     __fn_energy__ = lambda x: x
     __fn_distance__ = lambda x: x
     __fn_forces__ = lambda x: x
+    __average_nb_atoms__ = None
+    __stats__ = {}
 
     def __init__(
         self,
@@ -93,26 +107,104 @@ def __init__(
     ) -> None:
         set_cache_dir(cache_dir)
         self.data = None
-        self._set_units(energy_unit, distance_unit)
         if not self.is_preprocessed():
-            logger.info("This dataset not available. Please open an issue on Github for the team to look into it.")
-            # entries = self.read_raw_entries()
-            # res = self.collate_list(entries)
-            # self.save_preprocess(res)
+            raise DatasetNotAvailableError(self.__name__)
         else:
             self.read_preprocess(overwrite_local_cache=overwrite_local_cache)
-            self._set_isolated_atom_energies()
+        self._post_init(overwrite_local_cache, energy_unit, distance_unit)
+
+    def _post_init(
+        self,
+        overwrite_local_cache: bool = False,
+        energy_unit: Optional[str] = None,
+        distance_unit: Optional[str] = None,
+    ) -> None:
+        self._set_units(None, None)
+        self._set_isolated_atom_energies()
+        self._precompute_statistics(overwrite_local_cache=overwrite_local_cache)
+        self._set_units(energy_unit, distance_unit)
+        self._convert_data()
+        self._set_isolated_atom_energies()
+
+    def _convert_data(self):
+        logger.info(
+            f"Converting {self.__name__} data to the following units:\n\
+                     Energy: {self.energy_unit},\n\
+                     Distance: {self.distance_unit},\n\
+                     Forces: {self.force_unit if self.__force_methods__ else 'None'}"
+        )
+        for key in self.data_keys:
+            self.data[key] = self._convert_on_loading(self.data[key], key)
+
+    def _precompute_statistics(self, overwrite_local_cache: bool = False):
+        local_path = p_join(self.preprocess_path, "stats.pkl")
+        if self.is_preprocessed_statistics() and not overwrite_local_cache:
+            stats = load_pkl(local_path)
+            logger.info("Loaded precomputed statistics")
+        else:
+            logger.info("Precomputing relevant statistics")
+            (formation_E_mean, formation_E_std, total_E_mean, total_E_std) = self._precompute_E()
+            forces_dict = self._precompute_F()
+            stats = {
+                "formation": {"energy": {"mean": formation_E_mean, "std": formation_E_std}, "forces": forces_dict},
+                "total": {"energy": {"mean": total_E_mean, "std": total_E_std}, "forces": forces_dict},
+            }
+            with open(local_path, "wb") as f:
+                pkl.dump(stats, f)
+        self._compute_average_nb_atoms()
+        self.__stats__ = stats
+
+    def _compute_average_nb_atoms(self):
+        self.__average_nb_atoms__ = np.mean(self.data["n_atoms"])
+
+    def _precompute_E(self):
+        splits_idx = self.data["position_idx_range"][:, 1]
+        s = np.array(self.data["atomic_inputs"][:, :2], dtype=int)
+        s[:, 1] += IsolatedAtomEnergyFactory.max_charge
+        matrixs = [matrix[s[:, 0], s[:, 1]] for matrix in self.__isolated_atom_energies__]
+        converted_energy_data = self.convert_energy(self.data["energies"])
+        # calculation per molecule formation energy statistics
+        E = []
+        for i, matrix in enumerate(matrixs):
+            c = np.cumsum(np.append([0], matrix))[splits_idx]
+            c[1:] = c[1:] - c[:-1]
+            E.append(converted_energy_data[:, i] - c)
+        E = np.array(E).T
+        formation_E_mean = np.nanmean(E, axis=0)
+        formation_E_std = np.nanstd(E, axis=0)
+        total_E_mean = np.nanmean(converted_energy_data, axis=0)
+        total_E_std = np.nanstd(converted_energy_data, axis=0)
+
+        return (
+            np.atleast_2d(formation_E_mean),
+            np.atleast_2d(formation_E_std),
+            np.atleast_2d(total_E_mean),
+            np.atleast_2d(total_E_std),
+        )
+
+    def _precompute_F(self):
+        if len(self.__force_methods__) == 0:
+            return NOT_DEFINED
+        converted_force_data = self.convert_forces(self.data["forces"])
+        force_mean = np.nanmean(converted_force_data, axis=0)
+        force_std = np.nanstd(converted_force_data, axis=0)
+        force_rms = np.sqrt(np.nanmean(converted_force_data**2, axis=0))
+        return {
+            "mean": np.atleast_2d(force_mean.mean(axis=0)),
+            "std": np.atleast_2d(force_std.mean(axis=0)),
+            "components": {"rms": force_rms, "std": force_std, "mean": force_mean},
+        }
 
     @property
     def numbers(self):
         if hasattr(self, "_numbers"):
             return self._numbers
-        self._numbers = np.array(list(set(self.data["atomic_inputs"][..., 0])), dtype=np.int32)
+        self._numbers = pd.unique(self.data["atomic_inputs"][..., 0]).astype(np.int32)
         return self._numbers
 
     @property
     def chemical_species(self):
-        return [chemical_symbols[z] for z in self.numbers]
+        return np.array(chemical_symbols)[self.numbers]
 
     @property
     def energy_unit(self):
@@ -211,10 +303,11 @@ def collate_list(self, list_entries):
         # concatenate entries
         res = {key: np.concatenate([r[key] for r in list_entries if r is not None], axis=0) for key in list_entries[0]}
 
-        csum = np.cumsum(res.pop("n_atoms"))
+        csum = np.cumsum(res.get("n_atoms"))
         x = np.zeros((csum.shape[0], 2), dtype=np.int32)
         x[1:, 0], x[:, 1] = csum[:-1], csum
         res["position_idx_range"] = x
+
         return res
 
     def save_preprocess(self, data_dict):
@@ -228,12 +321,25 @@ def save_preprocess(self, data_dict):
             push_remote(local_path, overwrite=True)
 
         # save smiles and subset
+        local_path = p_join(self.preprocess_path, "props.pkl")
         for key in ["name", "subset"]:
-            local_path = p_join(self.preprocess_path, f"{key}.npz")
-            uniques, inv_indices = np.unique(data_dict[key], return_inverse=True)
-            with open(local_path, "wb") as f:
-                np.savez_compressed(f, uniques=uniques, inv_indices=inv_indices)
-            push_remote(local_path)
+            data_dict[key] = np.unique(data_dict[key], return_inverse=True)
+
+        with open(local_path, "wb") as f:
+            pkl.dump(data_dict, f)
+        push_remote(local_path, overwrite=True)
+
+    def _convert_on_loading(self, x, key):
+        if key == "energies":
+            return self.convert_energy(x)
+        elif key == "forces":
+            return self.convert_forces(x)
+        elif key == "atomic_inputs":
+            x = np.array(x, dtype=np.float32)
+            x[:, -3:] = self.convert_distance(x[:, -3:])
+            return x
+        else:
+            return x
 
     def read_preprocess(self, overwrite_local_cache=False):
         logger.info("Reading preprocessed data")
@@ -247,44 +353,45 @@ def read_preprocess(self, overwrite_local_cache=False):
         for key in self.data_keys:
             filename = p_join(self.preprocess_path, f"{key}.mmap")
             pull_locally(filename, overwrite=overwrite_local_cache)
-            self.data[key] = np.memmap(
-                filename,
-                mode="r",
-                dtype=self.data_types[key],
-            ).reshape(self.data_shapes[key])
+            self.data[key] = np.memmap(filename, mode="r", dtype=self.data_types[key]).reshape(self.data_shapes[key])
+
+        filename = p_join(self.preprocess_path, "props.pkl")
+        pull_locally(filename, overwrite=overwrite_local_cache)
+        with open(filename, "rb") as f:
+            tmp = pkl.load(f)
+            for key in ["name", "subset", "n_atoms"]:
+                x = tmp.pop(key)
+                if len(x) == 2:
+                    self.data[key] = x[0][x[1]]
+                else:
+                    self.data[key] = x
 
         for key in self.data:
-            print(f"Loaded {key} with shape {self.data[key].shape}, dtype {self.data[key].dtype}")
-
-        for key in ["name", "subset"]:
-            filename = p_join(self.preprocess_path, f"{key}.npz")
-            pull_locally(filename)
-            self.data[key] = dict()
-            with open(filename, "rb") as f:
-                tmp = np.load(f)
-                for k in tmp:
-                    self.data[key][k] = tmp[k]
-                    print(f"Loaded {key}_{k} with shape {self.data[key][k].shape}, dtype {self.data[key][k].dtype}")
+            logger.info(f"Loaded {key} with shape {self.data[key].shape}, dtype {self.data[key].dtype}")
 
     def is_preprocessed(self):
         predicats = [copy_exists(p_join(self.preprocess_path, f"{key}.mmap")) for key in self.data_keys]
-        predicats += [copy_exists(p_join(self.preprocess_path, f"{x}.npz")) for x in ["name", "subset"]]
+        predicats += [copy_exists(p_join(self.preprocess_path, "props.pkl"))]
         return all(predicats)
 
-    def preprocess(self):
-        if not self.is_preprocessed():
+    def is_preprocessed_statistics(self):
+        return bool(copy_exists(p_join(self.preprocess_path, "stats.pkl")))
+
+    def preprocess(self, overwrite=False):
+        if overwrite or not self.is_preprocessed():
             entries = self.read_raw_entries()
             res = self.collate_list(entries)
             self.save_preprocess(res)
 
-    def save_xyz(self, idx: int, path: Optional[str] = None):
+    def save_xyz(self, idx: int, path: Optional[str] = None, name=None):
         """
         Save the entry at index idx as an extxyz file.
         """
         if path is None:
             path = os.getcwd()
         at = self.get_ase_atoms(idx, ext=True)
-        name = at.info["name"]
+        if name is not None:
+            name = at.info["name"]
         write_extxyz(p_join(path, f"{name}.xyz"), at)
 
     def get_ase_atoms(self, idx: int, ext=True):
@@ -305,7 +412,7 @@ def get_ase_atoms(self, idx: int, ext=True):
 
     @requires_package("dscribe")
     @requires_package("datamol")
-    def chemical_space(
+    def soap_descriptors(
         self,
         n_samples: Optional[Union[List[int], int]] = None,
         return_idxs: bool = True,
@@ -350,7 +457,7 @@ def chemical_space(
             idxs = list(range(len(self)))
         elif isinstance(n_samples, int):
             idxs = np.random.choice(len(self), size=n_samples, replace=False)
-        elif isinstance(n_samples, list):
+        else:  # list, set, np.ndarray
             idxs = n_samples
         datum = {}
         r_cut = soap_kwargs.pop("r_cut", 5.0)
@@ -383,7 +490,7 @@ def wrapper(idx):
             entry = self.get_ase_atoms(idx, ext=False)
             return soap.create(entry, centers=entry.positions)
 
-        descr = dm.parallelized(wrapper, idxs, progress=progress, scheduler="threads")
+        descr = dm.parallelized(wrapper, idxs, progress=progress, scheduler="threads", n_jobs=-1)
         datum["soap"] = np.vstack(descr)
         if return_idxs:
             datum["idxs"] = idxs
@@ -392,6 +499,12 @@ def wrapper(idx):
     def __len__(self):
         return self.data["energies"].shape[0]
 
+    def __smiles_converter__(self, x):
+        """util function to convert string to smiles: useful if the smiles is
+        encoded in a different format than its display format
+        """
+        return x
+
     def __getitem__(self, idx: int):
         shift = IsolatedAtomEnergyFactory.max_charge
         p_start, p_end = self.data["position_idx_range"][idx]
@@ -399,14 +512,14 @@ def __getitem__(self, idx: int):
         z, c, positions, energies = (
             np.array(input[:, 0], dtype=np.int32),
             np.array(input[:, 1], dtype=np.int32),
-            self.convert_distance(np.array(input[:, -3:], dtype=np.float32)),
-            self.convert_energy(np.array(self.data["energies"][idx], dtype=np.float32)),
+            np.array(input[:, -3:], dtype=np.float32),
+            np.array(self.data["energies"][idx], dtype=np.float32),
         )
-        name = self.data["name"]["uniques"][self.data["name"]["inv_indices"][idx]]
-        subset = self.data["subset"]["uniques"][self.data["subset"]["inv_indices"][idx]]
+        name = self.__smiles_converter__(self.data["name"][idx])
+        subset = self.data["subset"][idx]
 
         if "forces" in self.data:
-            forces = self.convert_forces(np.array(self.data["forces"][p_start:p_end], dtype=np.float32))
+            forces = np.array(self.data["forces"][p_start:p_end], dtype=np.float32)
         else:
             forces = None
         return Bunch(
@@ -425,3 +538,61 @@ def __str__(self):
 
     def __repr__(self):
         return f"{self.__name__}"
+
+    @property
+    def _stats(self):
+        return self.__stats__
+
+    @property
+    def average_n_atoms(self):
+        """
+        Average number of atoms in a molecule in the dataset.
+        """
+        if self.__average_nb_atoms__ is None:
+            raise StatisticsNotAvailableError(self.__name__)
+        return self.__average_nb_atoms__
+
+    def get_statistics(self, normalization: str = "formation", return_none: bool = True):
+        """
+        Get the statistics of the dataset.
+        normalization : str, optional
+            Type of energy, by default "formation", must be one of ["formation", "total"]
+        return_none : bool, optional
+            Whether to return None if the statistics for the forces are not available, by default True
+            Otherwise, the statistics for the forces are set to 0.0
+        """
+        stats = deepcopy(self._stats)
+        if len(stats) == 0:
+            raise StatisticsNotAvailableError(self.__name__)
+        if normalization not in POSSIBLE_NORMALIZATION:
+            raise NormalizationNotAvailableError(normalization)
+        selected_stats = stats[normalization]
+        if len(self.__force_methods__) == 0 and not return_none:
+            selected_stats.update(
+                {
+                    "forces": {
+                        "mean": np.array([0.0]),
+                        "std": np.array([0.0]),
+                        "components": {
+                            "mean": np.array([[0.0], [0.0], [0.0]]),
+                            "std": np.array([[0.0], [0.0], [0.0]]),
+                            "rms": np.array([[0.0], [0.0], [0.0]]),
+                        },
+                    }
+                }
+            )
+        # cycle trough dict to convert units
+        for key in selected_stats:
+            if key == "forces":
+                for key2 in selected_stats[key]:
+                    if key2 != "components":
+                        selected_stats[key][key2] = self.convert_forces(selected_stats[key][key2])
+                    else:
+                        for key2 in selected_stats[key]["components"]:
+                            selected_stats[key]["components"][key2] = self.convert_forces(
+                                selected_stats[key]["components"][key2]
+                            )
+            else:
+                for key2 in selected_stats[key]:
+                    selected_stats[key][key2] = self.convert_energy(selected_stats[key][key2])
+        return selected_stats
diff --git a/src/openqdc/datasets/comp6.py b/src/openqdc/datasets/comp6.py
index c95ec17..7b6890b 100644
--- a/src/openqdc/datasets/comp6.py
+++ b/src/openqdc/datasets/comp6.py
@@ -35,8 +35,8 @@ class COMP6(BaseDataset):
         "pbe-d3bj/def2-tzvp",
         "pbe/def2-tzvp",
         "svwn/def2-tzvp",
-        "wb97m-d3bj/def2-tzvp",
-        "wb97m/def2-tzvp",
+        # "wb97m-d3bj/def2-tzvp",
+        # "wb97m/def2-tzvp",
     ]
 
     energy_target_names = [
@@ -47,8 +47,8 @@ class COMP6(BaseDataset):
         "PBE-D3M(BJ):def2-tzvp",
         "PBE:def2-tzvp",
         "SVWN:def2-tzvp",
-        "WB97M-D3(BJ):def2-tzvp",
-        "WB97M:def2-tzvp",
+        # "WB97M-D3(BJ):def2-tzvp",
+        # "WB97M:def2-tzvp",
     ]
 
     __force_methods__ = [
@@ -59,6 +59,12 @@ class COMP6(BaseDataset):
         "Gradient",
     ]
 
+    def __smiles_converter__(self, x):
+        """util function to convert string to smiles: useful if the smiles is
+        encoded in a different format than its display format
+        """
+        return "-".join(x.decode("ascii").split("_")[:-1])
+
     def read_raw_entries(self):
         samples = []
         for subset in ["ani_md", "drugbank", "gdb7_9", "gdb10_13", "s66x8", "tripeptides"]:
diff --git a/src/openqdc/datasets/dummy.py b/src/openqdc/datasets/dummy.py
index 4e1ff17..99b5106 100644
--- a/src/openqdc/datasets/dummy.py
+++ b/src/openqdc/datasets/dummy.py
@@ -1,7 +1,10 @@
 import numpy as np  # noqa
+from numpy import array
 from sklearn.utils import Bunch
 
 from openqdc.datasets.base import BaseDataset
+from openqdc.utils.atomization_energies import IsolatedAtomEnergyFactory
+from openqdc.utils.constants import NOT_DEFINED
 
 
 class Dummy(BaseDataset):
@@ -10,38 +13,66 @@ class Dummy(BaseDataset):
     """
 
     __name__ = "dummy"
-    __energy_methods__ = ["I_solved_the_schrodinger_equation_by_hand"]
-    __force_methods__ = ["I_made_up_random_forces"]
+    __energy_methods__ = ["I_solved_the_schrodinger_equation_by_hand", "PM6"]
+    __force_methods__ = ["I_made_up_random_forces", "writing_1_to_every_coordinate"]
     __energy_unit__ = "kcal/mol"
     __distance_unit__ = "ang"
     __forces_unit__ = "kcal/mol/ang"
 
-    energy_target_names = ["energy"]
+    energy_target_names = [f"energy{i}" for i in range(len(__energy_methods__))]
 
-    force_target_names = ["forces"]
+    force_target_names = [f"forces{i}" for i in range(len(__force_methods__))]
+    __isolated_atom_energies__ = []
+    __average_n_atoms__ = 20
+
+    @property
+    def _stats(self):
+        return {
+            "formation": {
+                "energy": {
+                    "mean": array([-12.94348027, -9.83037297]),
+                    "std": array([4.39971409, 3.3574188]),
+                },
+                "forces": NOT_DEFINED,
+            },
+            "total": {
+                "energy": {
+                    "mean": array([-89.44242, -1740.5336]),
+                    "std": array([29.599571, 791.48663]),
+                },
+                "forces": NOT_DEFINED,
+            },
+        }
 
     def __init__(self, energy_unit=None, distance_unit=None, cache_dir=None) -> None:
         try:
             super().__init__(energy_unit=energy_unit, distance_unit=distance_unit, cache_dir=cache_dir)
+
         except:  # noqa
             pass
+        self._set_isolated_atom_energies()
+
+    def is_preprocessed(self):
+        return True
 
     def read_raw_entries(self):
         pass
 
     def __len__(self):
-        return 999999999
+        return 9999
 
     def __getitem__(self, idx: int):
-        size = np.random.randint(1, 250)
+        shift = IsolatedAtomEnergyFactory.max_charge
+        size = np.random.randint(1, 100)
         z = np.random.randint(1, 100, size)
+        c = np.random.randint(-1, 2, size)
         return Bunch(
             positions=np.random.rand(size, 3) * 10,
             atomic_numbers=z,
-            charges=np.random.randint(-1, 2, size),
-            e0=np.zeros(size),
-            energies=np.random.rand(1) * 100,
+            charges=c,
+            e0=self.__isolated_atom_energies__[..., z, c + shift].T,
+            energies=np.random.randn(len(self.__energy_methods__)),
             name="dummy_{}".format(idx),
             subset="dummy",
-            forces=np.random.rand(size, 3) * 100,
+            forces=(np.random.randn(size, 3, len(self.__force_methods__)) * 100),
         )
diff --git a/src/openqdc/datasets/gdml.py b/src/openqdc/datasets/gdml.py
index 789f84a..e40b3fa 100644
--- a/src/openqdc/datasets/gdml.py
+++ b/src/openqdc/datasets/gdml.py
@@ -32,7 +32,7 @@ class GDML(BaseDataset):
     __energy_methods__ = [
         "ccsd/cc-pvdz",
         "ccsd(t)/cc-pvdz",
-        "pbe/mbd",  # MD22
+        # "pbe/mbd",  # MD22
         # "pbe+mbd/tight", #MD22
         "pbe/vdw-ts",  # MD17
     ]
@@ -46,7 +46,7 @@ class GDML(BaseDataset):
     __force_methods__ = [
         "ccsd/cc-pvdz",
         "ccsd(t)/cc-pvdz",
-        "pbe/mbd",  # MD22
+        # "pbe/mbd",  # MD22
         # "pbe+mbd/tight", #MD22
         "pbe/vdw-ts",  # MD17
     ]
diff --git a/src/openqdc/datasets/iso_17.py b/src/openqdc/datasets/iso_17.py
index 735ae67..4553ec1 100644
--- a/src/openqdc/datasets/iso_17.py
+++ b/src/openqdc/datasets/iso_17.py
@@ -43,6 +43,12 @@ class ISO17(BaseDataset):
     __distance_unit__ = "bohr"  # bohr
     __forces_unit__ = "ev/bohr"
 
+    def __smiles_converter__(self, x):
+        """util function to convert string to smiles: useful if the smiles is
+        encoded in a different format than its display format
+        """
+        return "-".join(x.decode("ascii").split("_")[:-1])
+
     def read_raw_entries(self):
         raw_path = p_join(self.root, "iso_17.h5")
         samples = read_qc_archive_h5(raw_path, "iso_17", self.energy_target_names, self.force_target_names)
diff --git a/src/openqdc/datasets/nabladft.py b/src/openqdc/datasets/nabladft.py
index e7d9eb8..0555cdc 100644
--- a/src/openqdc/datasets/nabladft.py
+++ b/src/openqdc/datasets/nabladft.py
@@ -4,30 +4,32 @@
 
 import datamol as dm
 import numpy as np
-from tqdm import tqdm
+import pandas as pd
 
 from openqdc.datasets.base import BaseDataset
+from openqdc.utils.molecule import z_to_formula
 from openqdc.utils.package_utils import requires_package
 
 
-def to_mol(entry) -> Dict[str, np.ndarray]:
+def to_mol(entry, metadata) -> Dict[str, np.ndarray]:
     Z, R, E, F = entry[:4]
     C = np.zeros_like(Z)
+    E[0] = metadata["DFT TOTAL ENERGY"]
 
     res = dict(
         atomic_inputs=np.concatenate((Z[:, None], C[:, None], R), axis=-1).astype(np.float32),
-        name=np.array([""]),
+        name=np.array([metadata["SMILES"]]),
         energies=E[:, None].astype(np.float32),
         forces=F[:, :, None].astype(np.float32),
         n_atoms=np.array([Z.shape[0]], dtype=np.int32),
-        subset=np.array(["nabla"]),
+        subset=np.array([z_to_formula(Z)]),
     )
 
     return res
 
 
 @requires_package("nablaDFT")
-def read_chunk_from_db(raw_path, start_idx, stop_idx, step_size=1000):
+def read_chunk_from_db(raw_path, start_idx, stop_idx, labels, step_size=1000):
     from nablaDFT.dataset import HamiltonianDatabase
 
     print(f"Loading from {start_idx} to {stop_idx}")
@@ -35,7 +37,13 @@ def read_chunk_from_db(raw_path, start_idx, stop_idx, step_size=1000):
     idxs = list(np.arange(start_idx, stop_idx))
     n, s = len(idxs), step_size
 
-    samples = [to_mol(entry) for i in tqdm(range(0, n, s)) for entry in db[idxs[i : i + s]]]
+    cursor = db._get_connection().cursor()
+    data_idxs = cursor.execute("""SELECT * FROM dataset_ids WHERE id IN (""" + str(idxs)[1:-1] + ")").fetchall()
+    c_idxs = [tuple(x[1:]) for x in data_idxs]
+
+    samples = [
+        to_mol(entry, labels[c_idxs[i + j]]) for i in range(0, n, s) for j, entry in enumerate(db[idxs[i : i + s]])
+    ]
     return samples
 
 
@@ -68,12 +76,16 @@ class NablaDFT(BaseDataset):
     def read_raw_entries(self):
         from nablaDFT.dataset import HamiltonianDatabase
 
+        label_path = p_join(self.root, "summary.csv")
+        df = pd.read_csv(label_path, usecols=["MOSES id", "CONFORMER id", "SMILES", "DFT TOTAL ENERGY"])
+        labels = df.set_index(keys=["MOSES id", "CONFORMER id"]).to_dict("index")
+
         raw_path = p_join(self.root, "dataset_full.db")
         train = HamiltonianDatabase(raw_path)
         n, c = len(train), 20
         step_size = int(np.ceil(n / os.cpu_count()))
 
-        fn = lambda i: read_chunk_from_db(raw_path, i * step_size, min((i + 1) * step_size, n))
+        fn = lambda i: read_chunk_from_db(raw_path, i * step_size, min((i + 1) * step_size, n), labels=labels)
         samples = dm.parallelized(
             fn, list(range(c)), n_jobs=c, progress=False, scheduler="threads"
         )  # don't use more than 1 job
diff --git a/src/openqdc/datasets/pcqm.py b/src/openqdc/datasets/pcqm.py
index 505eef1..d1a344c 100644
--- a/src/openqdc/datasets/pcqm.py
+++ b/src/openqdc/datasets/pcqm.py
@@ -11,7 +11,7 @@
 from loguru import logger
 
 from openqdc.datasets.base import BaseDataset
-from openqdc.utils.io import get_local_cache
+from openqdc.utils.io import get_local_cache, push_remote
 
 
 def flatten_dict(d, sep: str = "."):
@@ -80,27 +80,83 @@ def __init__(self, energy_unit=None, distance_unit=None) -> None:
     def root(self):
         return p_join(get_local_cache(), "pubchemqc")
 
-    def collate_list(self, list_entries, partial=False):
-        # default partial=False is necessary for compatibility with the base class
-        if partial:
-            predicat = list_entries is not None and len(list_entries) > 0
-            list_entries = [x for x in list_entries if x is not None]
-            return super().collate_list(list_entries) if predicat else None
+    @property
+    def preprocess_path(self):
+        path = p_join(self.root, "preprocessed", self.__name__)
+        os.makedirs(path, exist_ok=True)
+        return path
+
+    def collate_list(self, list_entries):
+        predicat = list_entries is not None and len(list_entries) > 0
+        list_entries = [x for x in list_entries if x is not None]
+        if predicat:
+            res = super().collate_list(list_entries)
         else:
-            n = 0
-            for i in range(len(list_entries)):
-                list_entries[i]["position_idx_range"] += n
-                n += list_entries[i]["position_idx_range"].max()
-            res = {key: np.concatenate([r[key] for r in list_entries], axis=0) for key in list_entries[0]}
-            return res
+            res = None
+        return res
 
     def read_raw_entries(self):
         arxiv_paths = glob(p_join(self.root, f"{self.__energy_methods__[0]}", "*.pkl"))
-        f = lambda x: self.collate_list(read_preprocessed_archive(x), partial=True)
+        f = lambda x: self.collate_list(read_preprocessed_archive(x))
         samples = dm.parallelized(f, arxiv_paths, n_jobs=1, progress=True)
         samples = [x for x in samples if x is not None]
         return samples
 
+    def preprocess(self, overwrite=False):
+        if overwrite or not self.is_preprocessed():
+            logger.info("Preprocessing data and saving it to cache.")
+            logger.info(
+                f"Dataset {self.__name__} data with the following units:\n"
+                f"Energy: {self.energy_unit}, Distance: {self.distance_unit}, "
+                f"Forces: {self.force_unit if self.__force_methods__ else 'None'}"
+            )
+            entries = self.read_raw_entries()
+            self.collate_and_save_list(entries)
+
+    def collate_and_save_list(self, list_entries):
+        n_molecules, n_atoms = 0, 0
+        for i in range(len(list_entries)):
+            list_entries[i]["position_idx_range"] += n_atoms
+            n_atoms += list_entries[i]["position_idx_range"].max()
+            n_molecules += list_entries[i]["position_idx_range"].shape[0]
+
+        for key in self.data_keys:
+            first = list_entries[0][key]
+            shape = (n_molecules, *first.shape[1:])
+            local_path = p_join(self.preprocess_path, f"{key}.mmap")
+            out = np.memmap(local_path, mode="w+", dtype=first.dtype, shape=shape)
+
+            start = 0
+            for i in range(len(list_entries)):
+                x = list_entries[i].pop(key)
+                n = x.shape[0]
+                out[start : start + n] = x
+                out.flush()
+            push_remote(local_path, overwrite=True)
+
+        # save smiles and subset
+        tmp, n = dict(name=[]), len(list_entries)
+        local_path = p_join(self.preprocess_path, "props.pkl")
+        names = [list_entries[i].pop("name") for i in range(n)]
+        f = lambda xs: [dm.to_inchikey(x) for x in xs]
+        res = dm.parallelized(f, names, n_jobs=-1, progress=False)
+        for x in res:
+            tmp["name"] += x
+        for key in ["subset", "n_atoms"]:
+            tmp[key] = []
+            for i in range(n):
+                tmp[key] += list(list_entries[i].pop(key))
+        with open(local_path, "wb") as f:
+            pkl.dump(tmp, f)
+        push_remote(local_path, overwrite=True)
+        # for key in ["name", "subset"]:
+        #     local_path = p_join(self.preprocess_path, f"{key}.npz")
+        #     x = [el for i in range(len(list_entries)) for el in list_entries[i].pop(key)]
+        #     uniques, inv_indices = np.unique(x, return_inverse=True)
+        #     with open(local_path, "wb") as f:
+        #         np.savez_compressed(f, uniques=uniques, inv_indices=inv_indices)
+        #     push_remote(local_path, overwrite=True)
+
 
 class PCQM_B3LYP(PCQM_PM6):
     __name__ = "pubchemqc_b3lyp"
diff --git a/src/openqdc/datasets/sn2_rxn.py b/src/openqdc/datasets/sn2_rxn.py
index 3e75e91..abcbd62 100644
--- a/src/openqdc/datasets/sn2_rxn.py
+++ b/src/openqdc/datasets/sn2_rxn.py
@@ -25,8 +25,38 @@ class SN2RXN(BaseDataset):
         "DSD-BLYP-D3(BJ):def2-TZVP Gradient",
     ]
 
+    def __smiles_converter__(self, x):
+        """util function to convert string to smiles: useful if the smiles is
+        encoded in a different format than its display format
+        """
+        return "-".join(x.decode("ascii").split("_")[:-1])
+
     def read_raw_entries(self):
         raw_path = p_join(self.root, "sn2_rxn.h5")
+
+        # raw_path = p_join(self.root, "sn2_reactions.npz")
+        # data = np.load(raw_path)
+
+        # # as example for accessing individual entries, print the data for entry idx=0
+        # idx = 0
+        # print("Data for entry " + str(idx)+":")
+        # print("Number of atoms")
+        # print(data["N"][idx])
+        # print("Energy [eV]")
+        # print(data["E"][idx])
+        # print("Total charge")
+        # print(data["Q"][idx])
+        # print("Dipole moment vector (with respect to [0.0 0.0 0.0]) [eA]")
+        # print(data["D"][idx,:])
+        # print("Nuclear charges")
+        # print(data["Z"][idx,:data["N"][idx]])
+        # print("Cartesian coordinates [A]")
+        # print(data["R"][idx,:data["N"][idx],:])
+        # print("Forces [eV/A]")
+        # print(data["F"][idx,:data["N"][idx],:])
+
+        # exit()
+
         samples = read_qc_archive_h5(raw_path, "sn2_rxn", self.energy_target_names, self.force_target_names)
 
         return samples
diff --git a/src/openqdc/datasets/solvated_peptides.py b/src/openqdc/datasets/solvated_peptides.py
index 9846bdf..216ecdd 100644
--- a/src/openqdc/datasets/solvated_peptides.py
+++ b/src/openqdc/datasets/solvated_peptides.py
@@ -27,6 +27,12 @@ class SolvatedPeptides(BaseDataset):
     __distance_unit__ = "bohr"
     __forces_unit__ = "hartree/bohr"
 
+    def __smiles_converter__(self, x):
+        """util function to convert string to smiles: useful if the smiles is
+        encoded in a different format than its display format
+        """
+        return "_".join(x.decode("ascii").split("_")[:-1])
+
     def read_raw_entries(self):
         raw_path = p_join(self.root, "solvated_peptides.h5")
         samples = read_qc_archive_h5(raw_path, "solvated_peptides", self.energy_target_names, self.force_target_names)
diff --git a/src/openqdc/raws/config_factory.py b/src/openqdc/raws/config_factory.py
index 38bec86..c8dddba 100644
--- a/src/openqdc/raws/config_factory.py
+++ b/src/openqdc/raws/config_factory.py
@@ -37,7 +37,7 @@ class DataConfigFactory:
 
     sn2_rxn = dict(
         dataset_name="sn2_rxn",
-        links={"sn2_rxn.hdf5.gz": "https://zenodo.org/record/3585800/files/212.hdf5.gz"},
+        links={"sn2_rxn.hdf5.gz": "https://zenodo.org/records/2605341/files/sn2_reactions.npz"},
     )
 
     # FROM: https://sites.uw.edu/wdbase/database-of-water-clusters/
diff --git a/src/openqdc/utils/atomization_energies.py b/src/openqdc/utils/atomization_energies.py
index 40d0d13..6a1a638 100644
--- a/src/openqdc/utils/atomization_energies.py
+++ b/src/openqdc/utils/atomization_energies.py
@@ -2,124 +2,126 @@
 
 import numpy as np
 from loguru import logger
+from rdkit import Chem
 
 from openqdc.utils.constants import MAX_ATOMIC_NUMBER
 
+atom_table = Chem.GetPeriodicTable()
+
 __all__ = ["chemical_symbols", "atomic_numbers", "IsolatedAtomEnergyFactory"]
 
 EF_KEY: TypeAlias = Tuple[str, int]
 
-ATOM_SPECIES = "H", "Li", "B", "C", "N", "O", "F", "Na", "Mg", "Si", "P", "S", "Cl", "K", "Ca", "Br", "I"
-# Energy in atomic unit/ Hartree / Ang
-
 # didn t calculate for Pd, Pt, Mo, Ni, Fe, Cu, see DESS
 atomic_numbers = {}
-chemical_symbols = [
-    "X",
-    "H",
-    "He",
-    "Li",
-    "Be",
-    "B",
-    "C",
-    "N",
-    "O",
-    "F",
-    "Ne",
-    "Na",
-    "Mg",
-    "Al",
-    "Si",
-    "P",
-    "S",
-    "Cl",
-    "Ar",
-    "K",
-    "Ca",
-    "Sc",
-    "Ti",
-    "V",
-    "Cr",
-    "Mn",
-    "Fe",
-    "Co",
-    "Ni",
-    "Cu",
-    "Zn",
-    "Ga",
-    "Ge",
-    "As",
-    "Se",
-    "Br",
-    "Kr",
-    "Rb",
-    "Sr",
-    "Y",
-    "Zr",
-    "Nb",
-    "Mo",
-    "Tc",
-    "Ru",
-    "Rh",
-    "Pd",
-    "Ag",
-    "Cd",
-    "In",
-    "Sn",
-    "Sb",
-    "Te",
-    "I",
-    "Xe",
-    "Cs",
-    "Ba",
-    "La",
-    "Ce",
-    "Pr",
-    "Nd",
-    "Pm",
-    "Sm",
-    "Eu",
-    "Gd",
-    "Tb",
-    "Dy",
-    "Ho",
-    "Er",
-    "Tm",
-    "Yb",
-    "Lu",
-    "Hf",
-    "Ta",
-    "W",
-    "Re",
-    "Os",
-    "Ir",
-    "Pt",
-    "Au",
-    "Hg",
-    "Tl",
-    "Pb",
-    "Bi",
-    "Po",
-    "At",
-    "Rn",
-    "Fr",
-    "Ra",
-    "Ac",
-    "Th",
-    "Pa",
-    "U",
-    "Np",
-    "Pu",
-    "Am",
-    "Cm",
-    "Bk",
-    "Cf",
-    "Es",
-    "Fm",
-    "Md",
-    "No",
-    "Lr",
-]
+chemical_symbols = np.array(
+    [
+        "X",
+        "H",
+        "He",
+        "Li",
+        "Be",
+        "B",
+        "C",
+        "N",
+        "O",
+        "F",
+        "Ne",
+        "Na",
+        "Mg",
+        "Al",
+        "Si",
+        "P",
+        "S",
+        "Cl",
+        "Ar",
+        "K",
+        "Ca",
+        "Sc",
+        "Ti",
+        "V",
+        "Cr",
+        "Mn",
+        "Fe",
+        "Co",
+        "Ni",
+        "Cu",
+        "Zn",
+        "Ga",
+        "Ge",
+        "As",
+        "Se",
+        "Br",
+        "Kr",
+        "Rb",
+        "Sr",
+        "Y",
+        "Zr",
+        "Nb",
+        "Mo",
+        "Tc",
+        "Ru",
+        "Rh",
+        "Pd",
+        "Ag",
+        "Cd",
+        "In",
+        "Sn",
+        "Sb",
+        "Te",
+        "I",
+        "Xe",
+        "Cs",
+        "Ba",
+        "La",
+        "Ce",
+        "Pr",
+        "Nd",
+        "Pm",
+        "Sm",
+        "Eu",
+        "Gd",
+        "Tb",
+        "Dy",
+        "Ho",
+        "Er",
+        "Tm",
+        "Yb",
+        "Lu",
+        "Hf",
+        "Ta",
+        "W",
+        "Re",
+        "Os",
+        "Ir",
+        "Pt",
+        "Au",
+        "Hg",
+        "Tl",
+        "Pb",
+        "Bi",
+        "Po",
+        "At",
+        "Rn",
+        "Fr",
+        "Ra",
+        "Ac",
+        "Th",
+        "Pa",
+        "U",
+        "Np",
+        "Pu",
+        "Am",
+        "Cm",
+        "Bk",
+        "Cf",
+        "Es",
+        "Fm",
+        "Md",
+        "No",
+        "Lr",
+    ]
+)
 
 
 for Z, symbol in enumerate(chemical_symbols):
@@ -131,7 +133,7 @@ class IsolatedAtomEnergyFactory:
     Factory method to get the isolated atom energies for a given level of theory.
     """
 
-    max_charge = 4
+    max_charge = 6
 
     def __init__(self):
         pass
@@ -207,7 +209,14 @@ def get_matrix(level_of_theory: str) -> np.ndarray:
         if tuple_hashmap is None:
             return matrix
         for key in tuple_hashmap.keys():
-            matrix[atomic_numbers[key[0]], key[1] + shift] = tuple_hashmap[key]
+            try:
+                matrix[atomic_numbers[key[0]], key[1] + shift] = tuple_hashmap[key]
+            except KeyError:
+                print(key, list(tuple_hashmap.items()))
+                print(key[0], "?", key[1], "?", shift)
+                print(matrix.shape, atomic_numbers[key[0]], key[1] + shift)
+                logger.warning(f"Isolated atom energies not found for {key} and level of theory {level_of_theory}")
+                matrix[atomic_numbers[key[0]], key[1] + shift] = 0
         return matrix
 
 
diff --git a/src/openqdc/utils/constants.py b/src/openqdc/utils/constants.py
index a8a8215..9244637 100644
--- a/src/openqdc/utils/constants.py
+++ b/src/openqdc/utils/constants.py
@@ -5,3 +5,15 @@
 HAR2EV = 27.211386246
 
 BOHR2ANG = 0.52917721092
+
+POSSIBLE_NORMALIZATION = ["formation", "total"]
+
+NOT_DEFINED = {
+    "mean": None,
+    "std": None,
+    "components": {
+        "mean": None,
+        "std": None,
+        "rms": None,
+    },
+}
diff --git a/src/openqdc/utils/exceptions.py b/src/openqdc/utils/exceptions.py
new file mode 100644
index 0000000..246d01c
--- /dev/null
+++ b/src/openqdc/utils/exceptions.py
@@ -0,0 +1,68 @@
+from typing import Final
+
+from openqdc.utils.constants import POSSIBLE_NORMALIZATION
+
+PROPERTY_NOT_AVAILABLE_ERROR: Final[
+    str
+] = """This property for this dataset not available.
+Please open an issue on Github for the team to look into it."""
+
+
+class OpenQDCException(Exception):
+    """Base exception for custom exceptions raised by the openQDC"""
+
+    def __init__(self, msg: str):
+        super().__init__(msg)
+        self.msg = msg
+
+    def __str__(self):
+        return self.msg
+
+
+class DatasetNotAvailableError(OpenQDCException):
+    """Raised when a dataset is not available"""
+
+    msg = "Dataset {dataset_name} is not available. Please open an issue on Github for the team to look into it."
+
+    def __init__(self, dataset_name):
+        super().__init__(self.msg.format(dataset_name=dataset_name))
+
+
+class StatisticsNotAvailableError(DatasetNotAvailableError):
+    """Raised when statistics are not available"""
+
+    msg = (
+        "Statistics for dataset {dataset_name} are not available."
+        + "Please open an issue on Github for the team to look into it."
+    )
+
+
+class NormalizationNotAvailableError(OpenQDCException):
+    """Raised when normalization is not available"""
+
+    def __init__(self, normalization):
+        msg = f"Normalization={normalization} is not valid. Must be one of {POSSIBLE_NORMALIZATION}"
+        super().__init__(msg)
+
+
+class ConversionNotDefinedError(OpenQDCException, ValueError):
+    """Raised when a conversion is not defined"""
+
+    _error_message = """
+    Conversion from {in_unit} to {out_unit} is not defined in the conversion registry.
+    To add a new conversion, use the following syntax or open an issue on Github for the team to look into it:
+
+    Conversion("{in_unit}", "{out_unit}", lambda x: x * conversion_factor)
+    """
+
+    def __init__(self, in_unit, out_unit):
+        super().__init__(self._error_message.format(in_unit=in_unit, out_unit=out_unit))
+
+
+class ConversionAlreadyDefined(ConversionNotDefinedError):
+    """Raised when a conversion is not defined"""
+
+    _error_message = """
+    Conversion from {in_unit} to {out_unit} is alread defined in the conversion registry.
+    To reuse the same metric, use get_conversion({in_unit}, {out_unit}).
+    """
diff --git a/src/openqdc/utils/molecule.py b/src/openqdc/utils/molecule.py
index cd2290f..82a58d2 100644
--- a/src/openqdc/utils/molecule.py
+++ b/src/openqdc/utils/molecule.py
@@ -1,9 +1,22 @@
+from typing import Any
+
 import numpy as np
+from numpy import ndarray
 from rdkit import Chem
 
+from openqdc.utils.atomization_energies import chemical_symbols
+
 atom_table = Chem.GetPeriodicTable()
 
 
+def z_to_formula(z):
+    u, c = np.unique(z, return_counts=True)
+    idxs = np.argsort(u)
+    u, c = u[idxs], c[idxs]
+
+    return "".join([f"{chemical_symbols[u[i]]}{c[i] if c[i] > 1 else ''}" for i in range(len(u))])
+
+
 def get_atomic_number(mol: Chem.Mol):
     """Returns atomic numbers for rdkit molecule"""
     return np.array([atom.GetAtomicNum() for atom in mol.GetAtoms()])
@@ -17,3 +30,76 @@ def get_atomic_charge(mol: Chem.Mol):
 def get_atomic_number_and_charge(mol: Chem.Mol):
     """Returns atoms number and charge for rdkit molecule"""
     return np.array([[atom.GetAtomicNum(), atom.GetFormalCharge()] for atom in mol.GetAtoms()])
+
+
+def rmsd(P: ndarray, Q: ndarray, **kwargs) -> float:
+    """
+    Calculate Root-mean-square deviation from two sets of vectors V and W.
+
+    Parameters
+    ----------
+    V : array
+        (N,D) matrix, where N is points and D is dimension.
+    W : array
+        (N,D) matrix, where N is points and D is dimension.
+
+    Returns
+    -------
+    rmsd : float
+        Root-mean-square deviation between the two vectors
+    """
+    diff = P - Q
+    return np.sqrt((diff * diff).sum() / P.shape[0])
+
+
+def kabsch_rmsd(
+    P: ndarray,
+    Q: ndarray,
+    translate: bool = False,
+    **kwargs: Any,
+) -> float:
+    """
+    Rotate matrix P unto Q using Kabsch algorithm and calculate the RMSD.
+
+    Parameters
+    ----------
+    P : array
+        (N,D) matrix, where N is points and D is dimension.
+    Q : array
+        (N,D) matrix, where N is points and D is dimension.
+    translate : bool
+        Use centroids to translate vector P and Q unto each other.
+
+    Returns
+    -------
+    rmsd : float
+        root-mean squared deviation
+    """
+
+    if translate:
+        Q = Q - Q.mean(axis=0)
+        P = P - P.mean(axis=0)
+
+    # Computation of the covariance matrix
+    C = np.dot(np.transpose(P), Q)
+
+    # Computation of the optimal rotation matrix
+    # This can be done using singular value decomposition (SVD)
+    # Getting the sign of the det(V)*(W) to decide
+    # whether we need to correct our rotation matrix to ensure a
+    # right-handed coordinate system.
+    # And finally calculating the optimal rotation matrix U
+    # see http://en.wikipedia.org/wiki/Kabsch_algorithm
+    V, S, W = np.linalg.svd(C)
+    d = (np.linalg.det(V) * np.linalg.det(W)) < 0.0
+
+    if d:
+        S[-1] = -S[-1]
+        V[:, -1] = -V[:, -1]
+
+    # Create Rotation matrix U
+    U = np.dot(V, W)
+
+    # Rotate P
+    P_prime = np.dot(P, U)
+    return rmsd(P_prime, Q)
diff --git a/src/openqdc/utils/preprocess.py b/src/openqdc/utils/preprocess.py
index 1142dca..b34499e 100644
--- a/src/openqdc/utils/preprocess.py
+++ b/src/openqdc/utils/preprocess.py
@@ -36,9 +36,11 @@
 def preprocess(dataset):
     if dataset not in options_map:
         dataset_id = int(dataset)
+        data_class = options[dataset_id]
+    else:
+        data_class = options_map[dataset]
 
-    data_class = options[dataset_id]
-    data_class().preprocess()
+    data_class().preprocess(overwrite=False)
     data = data_class()
     logger.info(f"Preprocessing {data.__name__}")
 
@@ -47,7 +49,7 @@ def preprocess(dataset):
         x = data[i]
         print(x.name, x.subset, end=" ")
         for k in x:
-            if x[k] is not None:
+            if isinstance(x[k], np.ndarray):
                 print(k, x[k].shape, end=" ")
         print()
 
diff --git a/src/openqdc/utils/units.py b/src/openqdc/utils/units.py
index fb895ce..f79ebce 100644
--- a/src/openqdc/utils/units.py
+++ b/src/openqdc/utils/units.py
@@ -1,5 +1,7 @@
 from typing import Callable
 
+from openqdc.utils.exceptions import ConversionAlreadyDefined, ConversionNotDefinedError
+
 CONVERSION_REGISTRY = {}
 
 
@@ -13,7 +15,7 @@ def __init__(self, in_unit: str, out_unit: str, func: Callable[[float], float]):
         name = "convert_" + in_unit.lower().strip() + "_to_" + out_unit.lower().strip()
 
         if name in CONVERSION_REGISTRY:
-            raise ValueError(f"{name} is already registered. To reuse the same metric, use Metric.get_by_name().")
+            raise ConversionAlreadyDefined(in_unit, out_unit)
         CONVERSION_REGISTRY[name] = self
 
         self.name = name
@@ -29,7 +31,7 @@ def get_conversion(in_unit: str, out_unit: str):
     if in_unit.lower().strip() == out_unit.lower().strip():
         return lambda x: x
     if name not in CONVERSION_REGISTRY:
-        raise ValueError(f"{name} is not a valid metric. Valid metrics are: {list(CONVERSION_REGISTRY.keys())}")
+        raise ConversionNotDefinedError(in_unit, out_unit)
     return CONVERSION_REGISTRY[name]
 
 
@@ -73,3 +75,10 @@ def get_conversion(in_unit: str, out_unit: str):
 Conversion("hartree/ang", "kcal/mol/ang", lambda x: get_conversion("hartree", "kcal/mol")(x))
 Conversion("hartree/ang", "hartree/bohr", lambda x: get_conversion("bohr", "ang")(x))
 Conversion("hartree/bohr", "hartree/ang", lambda x: get_conversion("ang", "bohr")(x))
+Conversion("kcal/mol/bohr", "Hartree/bohr", lambda x: get_conversion("kcal/mol", "hartree")(x))
+Conversion("ev/ang", "hartree/ang", lambda x: get_conversion("ev", "hartree")(x))
+Conversion("ev/bohr", "hartree/bohr", lambda x: get_conversion("ev", "hartree")(x))
+Conversion("ev/bohr", "ev/ang", lambda x: get_conversion("ang", "bohr")(x))
+Conversion("ev/bohr", "kcal/mol/ang", lambda x: get_conversion("ang", "bohr")(get_conversion("ev", "kcal/mol")(x)))
+Conversion("kcal/mol/bohr", "kcal/mol/ang", lambda x: get_conversion("ang", "bohr")(x))
+Conversion("ev/ang", "kcal/mol/ang", lambda x: get_conversion("ev", "kcal/mol")(x))