diff --git a/.gitignore b/.gitignore
index 02fbeb97..5988c155 100644
--- a/.gitignore
+++ b/.gitignore
@@ -27,6 +27,9 @@ tests/echo*
 tests/fingerprints/
 tmp_*
 *_kliff_trained/
+tests/uq/*.pkl
+tests/uq/*.json
+tests/uq/kliff_saved_model
 
 # dataset
 Si_training_set_4_configs
diff --git a/docs/source/auto_examples/auto_examples_jupyter.zip b/docs/source/auto_examples/auto_examples_jupyter.zip
index 09072aba..da180fb1 100644
Binary files a/docs/source/auto_examples/auto_examples_jupyter.zip and b/docs/source/auto_examples/auto_examples_jupyter.zip differ
diff --git a/docs/source/auto_examples/auto_examples_python.zip b/docs/source/auto_examples/auto_examples_python.zip
index e749cbd7..8be3ecb3 100644
Binary files a/docs/source/auto_examples/auto_examples_python.zip and b/docs/source/auto_examples/auto_examples_python.zip differ
diff --git a/docs/source/auto_examples/example_uq_bootstrap.ipynb b/docs/source/auto_examples/example_uq_bootstrap.ipynb
new file mode 100644
index 00000000..4f1eb182
--- /dev/null
+++ b/docs/source/auto_examples/example_uq_bootstrap.ipynb
@@ -0,0 +1,158 @@
+{
+  "cells": [
+    {
+      "cell_type": "code",
+      "execution_count": null,
+      "metadata": {
+        "collapsed": false
+      },
+      "outputs": [],
+      "source": [
+        "%matplotlib inline"
+      ]
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "\n\n# Bootstrapping\n\nIn this example, we demonstrate how to perform uncertainty quantification (UQ) using\nbootstrap method. We use a Stillinger-Weber (SW) potential for silicon that is archived\nin OpenKIM_.\n\nFor simplicity, we only set the energy-scaling parameters, i.e., ``A`` and ``lambda`` as\nthe tunable parameters. These parameters will be calibrated to energies and forces of a\nsmall dataset, consisting of 4 compressed and stretched configurations of diamond silicon\nstructure.\n"
+      ]
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "To start, let's first install the SW model::\n\n   $ kim-api-collections-management install user SW_StillingerWeber_1985_Si__MO_405512056662_006\n\n.. seealso::\n   This installs the model and its driver into the ``User Collection``. See\n   `install_model` for more information about installing KIM models.\n\n"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": null,
+      "metadata": {
+        "collapsed": false
+      },
+      "outputs": [],
+      "source": [
+        "import matplotlib.pyplot as plt\nimport numpy as np\n\nfrom kliff.calculators import Calculator\nfrom kliff.dataset import Dataset\nfrom kliff.loss import Loss\nfrom kliff.models import KIMModel\nfrom kliff.uq.bootstrap import BootstrapEmpiricalModel\nfrom kliff.utils import download_dataset"
+      ]
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "Before running bootstrap, we need to define a loss function and train the model. More\ndetail information about this step can be found in `tut_kim_sw` and\n`tut_params_transform`.\n\n"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": null,
+      "metadata": {
+        "collapsed": false
+      },
+      "outputs": [],
+      "source": [
+        "# Create the model\nmodel = KIMModel(model_name=\"SW_StillingerWeber_1985_Si__MO_405512056662_006\")\n\n# Set the tunable parameters and the initial guess\nopt_params = {\"A\": [[\"default\"]], \"lambda\": [[\"default\"]]}\n\nmodel.set_opt_params(**opt_params)\nmodel.echo_opt_params()\n\n# Get the dataset\ndataset_path = download_dataset(dataset_name=\"Si_training_set_4_configs\")\n# Read the dataset\ntset = Dataset(dataset_path)\nconfigs = tset.get_configs()\n\n# Create calculator\ncalc = Calculator(model)\n# Only use the forces data\nca = calc.create(configs, use_energy=False, use_forces=True)\n\n# Instantiate the loss function\nresidual_data = {\"normalize_by_natoms\": False}\nloss = Loss(calc, residual_data=residual_data)"
+      ]
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "To perform UQ by bootstrapping, the general workflow starts by instantiating\n:class:`~kliff.uq.bootstrap.BootstrapEmpiricalModel`, or\n:class:`~kliff.uq.bootstrap.BootstrapNeuralNetworkModel` if using a neural network\npotential.\n\n"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": null,
+      "metadata": {
+        "collapsed": false
+      },
+      "outputs": [],
+      "source": [
+        "# Instantiate bootstrap class object\nBS = BootstrapEmpiricalModel(loss, seed=1717)"
+      ]
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "Then, we generate some bootstrap compute arguments. This is equivalent to generating\nbootstrap data. Typically, we just need to specify how many bootstrap data samples to\ngenerate. Additionally, if we call ``generate_bootstrap_compute_arguments`` multiple\ntimes, the new generated data samples will be appended to the previously generated data\nsamples. This is also the behavior if we read the data samples from the previously\nexported file.\n\n"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": null,
+      "metadata": {
+        "collapsed": false
+      },
+      "outputs": [],
+      "source": [
+        "# Generate bootstrap compute arguments\nBS.generate_bootstrap_compute_arguments(100)"
+      ]
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "Finally, we will iterate over these bootstrap data samples and train the potential\nusing each data sample. The resulting optimal parameters from each data sample give a\nsingle sample of parameters. By iterating over all data samples, then we will get an\nensemble of parameters.\n\nNote that the mapping from the bootstrap dataset to the parameters involve optimization.\nWe suggest to use the same mapping, i.e., the same optimizer setting, in each iteration.\nThis includes using the same set of initial parameter guess. In the case when the loss\nfunction has multiple local minima, we don't want the parameter ensemble to be biased\non the results of the other optimizations. For neural network model, we need to reset\nthe initial parameter value, which is done internally.\n\n"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": null,
+      "metadata": {
+        "collapsed": false
+      },
+      "outputs": [],
+      "source": [
+        "# Run bootstrap\nmin_kwargs = dict(method=\"lm\")  # Optimizer setting\ninitial_guess = calc.get_opt_params()  # Initial guess in the optimization\nBS.run(min_kwargs=min_kwargs, initial_guess=initial_guess)"
+      ]
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "The resulting parameter ensemble can be accessed in `BS.samples` as a `np.ndarray`.\nThen, we can plot the distribution of the parameters, as an example, or propagate the\nerror to the target quantities we want to study.\n\n"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": null,
+      "metadata": {
+        "collapsed": false
+      },
+      "outputs": [],
+      "source": [
+        "# Plot the distribution of the parameters\nplt.figure()\nplt.plot(*(BS.samples.T), \".\", alpha=0.5)\nparam_names = list(opt_params.keys())\nplt.xlabel(param_names[0])\nplt.ylabel(param_names[1])\nplt.show()"
+      ]
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "\n"
+      ]
+    }
+  ],
+  "metadata": {
+    "kernelspec": {
+      "display_name": "Python 3",
+      "language": "python",
+      "name": "python3"
+    },
+    "language_info": {
+      "codemirror_mode": {
+        "name": "ipython",
+        "version": 3
+      },
+      "file_extension": ".py",
+      "mimetype": "text/x-python",
+      "name": "python",
+      "nbconvert_exporter": "python",
+      "pygments_lexer": "ipython3",
+      "version": "3.8.10"
+    }
+  },
+  "nbformat": 4,
+  "nbformat_minor": 0
+}
\ No newline at end of file
diff --git a/docs/source/auto_examples/example_uq_bootstrap.py b/docs/source/auto_examples/example_uq_bootstrap.py
new file mode 100644
index 00000000..f253fc84
--- /dev/null
+++ b/docs/source/auto_examples/example_uq_bootstrap.py
@@ -0,0 +1,123 @@
+"""
+.. _tut_bootstrap:
+
+Bootstrapping
+=============
+
+In this example, we demonstrate how to perform uncertainty quantification (UQ) using
+bootstrap method. We use a Stillinger-Weber (SW) potential for silicon that is archived
+in OpenKIM_.
+
+For simplicity, we only set the energy-scaling parameters, i.e., ``A`` and ``lambda`` as
+the tunable parameters. These parameters will be calibrated to energies and forces of a
+small dataset, consisting of 4 compressed and stretched configurations of diamond silicon
+structure.
+"""
+
+
+##########################################################################################
+# To start, let's first install the SW model::
+#
+#    $ kim-api-collections-management install user SW_StillingerWeber_1985_Si__MO_405512056662_006
+#
+# .. seealso::
+#    This installs the model and its driver into the ``User Collection``. See
+#    :ref:`install_model` for more information about installing KIM models.
+
+
+import matplotlib.pyplot as plt
+import numpy as np
+
+from kliff.calculators import Calculator
+from kliff.dataset import Dataset
+from kliff.loss import Loss
+from kliff.models import KIMModel
+from kliff.uq.bootstrap import BootstrapEmpiricalModel
+from kliff.utils import download_dataset
+
+##########################################################################################
+# Before running bootstrap, we need to define a loss function and train the model. More
+# detail information about this step can be found in :ref:`tut_kim_sw` and
+# :ref:`tut_params_transform`.
+
+# Create the model
+model = KIMModel(model_name="SW_StillingerWeber_1985_Si__MO_405512056662_006")
+
+# Set the tunable parameters and the initial guess
+opt_params = {"A": [["default"]], "lambda": [["default"]]}
+
+model.set_opt_params(**opt_params)
+model.echo_opt_params()
+
+# Get the dataset
+dataset_path = download_dataset(dataset_name="Si_training_set_4_configs")
+# Read the dataset
+tset = Dataset(dataset_path)
+configs = tset.get_configs()
+
+# Create calculator
+calc = Calculator(model)
+# Only use the forces data
+ca = calc.create(configs, use_energy=False, use_forces=True)
+
+# Instantiate the loss function
+residual_data = {"normalize_by_natoms": False}
+loss = Loss(calc, residual_data=residual_data)
+
+##########################################################################################
+# To perform UQ by bootstrapping, the general workflow starts by instantiating
+# :class:`~kliff.uq.bootstrap.BootstrapEmpiricalModel`, or
+# :class:`~kliff.uq.bootstrap.BootstrapNeuralNetworkModel` if using a neural network
+# potential.
+
+
+# Instantiate bootstrap class object
+BS = BootstrapEmpiricalModel(loss, seed=1717)
+
+##########################################################################################
+# Then, we generate some bootstrap compute arguments. This is equivalent to generating
+# bootstrap data. Typically, we just need to specify how many bootstrap data samples to
+# generate. Additionally, if we call ``generate_bootstrap_compute_arguments`` multiple
+# times, the new generated data samples will be appended to the previously generated data
+# samples. This is also the behavior if we read the data samples from the previously
+# exported file.
+
+
+# Generate bootstrap compute arguments
+BS.generate_bootstrap_compute_arguments(100)
+
+##########################################################################################
+# Finally, we will iterate over these bootstrap data samples and train the potential
+# using each data sample. The resulting optimal parameters from each data sample give a
+# single sample of parameters. By iterating over all data samples, then we will get an
+# ensemble of parameters.
+#
+# Note that the mapping from the bootstrap dataset to the parameters involve optimization.
+# We suggest to use the same mapping, i.e., the same optimizer setting, in each iteration.
+# This includes using the same set of initial parameter guess. In the case when the loss
+# function has multiple local minima, we don't want the parameter ensemble to be biased
+# on the results of the other optimizations. For neural network model, we need to reset
+# the initial parameter value, which is done internally.
+
+
+# Run bootstrap
+min_kwargs = dict(method="lm")  # Optimizer setting
+initial_guess = calc.get_opt_params()  # Initial guess in the optimization
+BS.run(min_kwargs=min_kwargs, initial_guess=initial_guess)
+
+##########################################################################################
+# The resulting parameter ensemble can be accessed in `BS.samples` as a `np.ndarray`.
+# Then, we can plot the distribution of the parameters, as an example, or propagate the
+# error to the target quantities we want to study.
+
+
+# Plot the distribution of the parameters
+plt.figure()
+plt.plot(*(BS.samples.T), ".", alpha=0.5)
+param_names = list(opt_params.keys())
+plt.xlabel(param_names[0])
+plt.ylabel(param_names[1])
+plt.show()
+
+##########################################################################################
+# .. _OpenKIM: https://openkim.org
diff --git a/docs/source/auto_examples/example_uq_bootstrap.py.md5 b/docs/source/auto_examples/example_uq_bootstrap.py.md5
new file mode 100644
index 00000000..13f3a770
--- /dev/null
+++ b/docs/source/auto_examples/example_uq_bootstrap.py.md5
@@ -0,0 +1 @@
+d16579e397f3c9e5d2537a623bd65313
\ No newline at end of file
diff --git a/docs/source/auto_examples/example_uq_bootstrap.rst b/docs/source/auto_examples/example_uq_bootstrap.rst
new file mode 100644
index 00000000..20eb5467
--- /dev/null
+++ b/docs/source/auto_examples/example_uq_bootstrap.rst
@@ -0,0 +1,679 @@
+
+.. DO NOT EDIT.
+.. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY.
+.. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE:
+.. "auto_examples/example_uq_bootstrap.py"
+.. LINE NUMBERS ARE GIVEN BELOW.
+
+.. only:: html
+
+    .. note::
+        :class: sphx-glr-download-link-note
+
+        :ref:`Go to the end <sphx_glr_download_auto_examples_example_uq_bootstrap.py>`
+        to download the full example code
+
+.. rst-class:: sphx-glr-example-title
+
+.. _sphx_glr_auto_examples_example_uq_bootstrap.py:
+
+
+.. _tut_bootstrap:
+
+Bootstrapping
+=============
+
+In this example, we demonstrate how to perform uncertainty quantification (UQ) using
+bootstrap method. We use a Stillinger-Weber (SW) potential for silicon that is archived
+in OpenKIM_.
+
+For simplicity, we only set the energy-scaling parameters, i.e., ``A`` and ``lambda`` as
+the tunable parameters. These parameters will be calibrated to energies and forces of a
+small dataset, consisting of 4 compressed and stretched configurations of diamond silicon
+structure.
+
+.. GENERATED FROM PYTHON SOURCE LINES 19-26
+
+To start, let's first install the SW model::
+
+   $ kim-api-collections-management install user SW_StillingerWeber_1985_Si__MO_405512056662_006
+
+.. seealso::
+   This installs the model and its driver into the ``User Collection``. See
+   :ref:`install_model` for more information about installing KIM models.
+
+.. GENERATED FROM PYTHON SOURCE LINES 26-38
+
+.. code-block:: default
+
+
+
+    import matplotlib.pyplot as plt
+    import numpy as np
+
+    from kliff.calculators import Calculator
+    from kliff.dataset import Dataset
+    from kliff.loss import Loss
+    from kliff.models import KIMModel
+    from kliff.uq.bootstrap import BootstrapEmpiricalModel
+    from kliff.utils import download_dataset
+
+
+
+
+
+
+
+
+.. GENERATED FROM PYTHON SOURCE LINES 39-42
+
+Before running bootstrap, we need to define a loss function and train the model. More
+detail information about this step can be found in :ref:`tut_kim_sw` and
+:ref:`tut_params_transform`.
+
+.. GENERATED FROM PYTHON SOURCE LINES 42-67
+
+.. code-block:: default
+
+
+    # Create the model
+    model = KIMModel(model_name="SW_StillingerWeber_1985_Si__MO_405512056662_006")
+
+    # Set the tunable parameters and the initial guess
+    opt_params = {"A": [["default"]], "lambda": [["default"]]}
+
+    model.set_opt_params(**opt_params)
+    model.echo_opt_params()
+
+    # Get the dataset
+    dataset_path = download_dataset(dataset_name="Si_training_set_4_configs")
+    # Read the dataset
+    tset = Dataset(dataset_path)
+    configs = tset.get_configs()
+
+    # Create calculator
+    calc = Calculator(model)
+    # Only use the forces data
+    ca = calc.create(configs, use_energy=False, use_forces=True)
+
+    # Instantiate the loss function
+    residual_data = {"normalize_by_natoms": False}
+    loss = Loss(calc, residual_data=residual_data)
+
+
+
+
+
+.. rst-class:: sphx-glr-script-out
+
+ .. code-block:: none
+
+    #================================================================================
+    # Model parameters that are optimized.
+    # Note that the parameters are in the transformed space if 
+    # `params_transform` is provided when instantiating the model.
+    #================================================================================
+
+    A 1
+      1.5284847919791400e+01 
+
+    lambda 1
+      4.5532200000000003e+01 
+
+
+    2023-04-03 19:16:23.548 | INFO     | kliff.dataset.dataset:_read:398 - 4 configurations read from /home/yonatank/modules/kliff/examples/Si_training_set_4_configs
+    2023-04-03 19:16:23.550 | INFO     | kliff.calculators.calculator:create:107 - Create calculator for 4 configurations.
+
+
+
+
+.. GENERATED FROM PYTHON SOURCE LINES 68-72
+
+To perform UQ by bootstrapping, the general workflow starts by instantiating
+:class:`~kliff.uq.bootstrap.BootstrapEmpiricalModel`, or
+:class:`~kliff.uq.bootstrap.BootstrapNeuralNetworkModel` if using a neural network
+potential.
+
+.. GENERATED FROM PYTHON SOURCE LINES 72-77
+
+.. code-block:: default
+
+
+
+    # Instantiate bootstrap class object
+    BS = BootstrapEmpiricalModel(loss, seed=1717)
+
+
+
+
+
+
+
+
+.. GENERATED FROM PYTHON SOURCE LINES 78-84
+
+Then, we generate some bootstrap compute arguments. This is equivalent to generating
+bootstrap data. Typically, we just need to specify how many bootstrap data samples to
+generate. Additionally, if we call ``generate_bootstrap_compute_arguments`` multiple
+times, the new generated data samples will be appended to the previously generated data
+samples. This is also the behavior if we read the data samples from the previously
+exported file.
+
+.. GENERATED FROM PYTHON SOURCE LINES 84-89
+
+.. code-block:: default
+
+
+
+    # Generate bootstrap compute arguments
+    BS.generate_bootstrap_compute_arguments(100)
+
+
+
+
+
+
+
+
+.. GENERATED FROM PYTHON SOURCE LINES 90-101
+
+Finally, we will iterate over these bootstrap data samples and train the potential
+using each data sample. The resulting optimal parameters from each data sample give a
+single sample of parameters. By iterating over all data samples, then we will get an
+ensemble of parameters.
+
+Note that the mapping from the bootstrap dataset to the parameters involve optimization.
+We suggest to use the same mapping, i.e., the same optimizer setting, in each iteration.
+This includes using the same set of initial parameter guess. In the case when the loss
+function has multiple local minima, we don't want the parameter ensemble to be biased
+on the results of the other optimizations. For neural network model, we need to reset
+the initial parameter value, which is done internally.
+
+.. GENERATED FROM PYTHON SOURCE LINES 101-108
+
+.. code-block:: default
+
+
+
+    # Run bootstrap
+    min_kwargs = dict(method="lm")  # Optimizer setting
+    initial_guess = calc.get_opt_params()  # Initial guess in the optimization
+    BS.run(min_kwargs=min_kwargs, initial_guess=initial_guess)
+
+
+
+
+
+.. rst-class:: sphx-glr-script-out
+
+ .. code-block:: none
+
+    2023-04-03 19:16:23.553 | INFO     | kliff.loss:minimize:312 - Start minimization using method: lm.
+    2023-04-03 19:16:23.554 | INFO     | kliff.loss:_scipy_optimize:431 - Running in serial mode.
+    2023-04-03 19:16:23.596 | INFO     | kliff.loss:minimize:314 - Finish minimization using method: lm.
+    2023-04-03 19:16:23.597 | INFO     | kliff.loss:minimize:312 - Start minimization using method: lm.
+    2023-04-03 19:16:23.597 | INFO     | kliff.loss:_scipy_optimize:431 - Running in serial mode.
+    2023-04-03 19:16:23.638 | INFO     | kliff.loss:minimize:314 - Finish minimization using method: lm.
+    2023-04-03 19:16:23.638 | INFO     | kliff.loss:minimize:312 - Start minimization using method: lm.
+    2023-04-03 19:16:23.639 | INFO     | kliff.loss:_scipy_optimize:431 - Running in serial mode.
+    2023-04-03 19:16:23.679 | INFO     | kliff.loss:minimize:314 - Finish minimization using method: lm.
+    2023-04-03 19:16:23.679 | INFO     | kliff.loss:minimize:312 - Start minimization using method: lm.
+    2023-04-03 19:16:23.680 | INFO     | kliff.loss:_scipy_optimize:431 - Running in serial mode.
+    2023-04-03 19:16:23.721 | INFO     | kliff.loss:minimize:314 - Finish minimization using method: lm.
+    2023-04-03 19:16:23.722 | INFO     | kliff.loss:minimize:312 - Start minimization using method: lm.
+    2023-04-03 19:16:23.722 | INFO     | kliff.loss:_scipy_optimize:431 - Running in serial mode.
+    2023-04-03 19:16:23.762 | INFO     | kliff.loss:minimize:314 - Finish minimization using method: lm.
+    2023-04-03 19:16:23.762 | INFO     | kliff.loss:minimize:312 - Start minimization using method: lm.
+    2023-04-03 19:16:23.763 | INFO     | kliff.loss:_scipy_optimize:431 - Running in serial mode.
+    2023-04-03 19:16:23.805 | INFO     | kliff.loss:minimize:314 - Finish minimization using method: lm.
+    2023-04-03 19:16:23.806 | INFO     | kliff.loss:minimize:312 - Start minimization using method: lm.
+    2023-04-03 19:16:23.806 | INFO     | kliff.loss:_scipy_optimize:431 - Running in serial mode.
+    2023-04-03 19:16:23.848 | INFO     | kliff.loss:minimize:314 - Finish minimization using method: lm.
+    2023-04-03 19:16:23.849 | INFO     | kliff.loss:minimize:312 - Start minimization using method: lm.
+    2023-04-03 19:16:23.849 | INFO     | kliff.loss:_scipy_optimize:431 - Running in serial mode.
+    2023-04-03 19:16:23.892 | INFO     | kliff.loss:minimize:314 - Finish minimization using method: lm.
+    2023-04-03 19:16:23.892 | INFO     | kliff.loss:minimize:312 - Start minimization using method: lm.
+    2023-04-03 19:16:23.893 | INFO     | kliff.loss:_scipy_optimize:431 - Running in serial mode.
+    2023-04-03 19:16:23.937 | INFO     | kliff.loss:minimize:314 - Finish minimization using method: lm.
+    2023-04-03 19:16:23.938 | INFO     | kliff.loss:minimize:312 - Start minimization using method: lm.
+    2023-04-03 19:16:23.938 | INFO     | kliff.loss:_scipy_optimize:431 - Running in serial mode.
+    2023-04-03 19:16:23.978 | INFO     | kliff.loss:minimize:314 - Finish minimization using method: lm.
+    2023-04-03 19:16:23.979 | INFO     | kliff.loss:minimize:312 - Start minimization using method: lm.
+    2023-04-03 19:16:23.980 | INFO     | kliff.loss:_scipy_optimize:431 - Running in serial mode.
+    2023-04-03 19:16:24.019 | INFO     | kliff.loss:minimize:314 - Finish minimization using method: lm.
+    2023-04-03 19:16:24.020 | INFO     | kliff.loss:minimize:312 - Start minimization using method: lm.
+    2023-04-03 19:16:24.020 | INFO     | kliff.loss:_scipy_optimize:431 - Running in serial mode.
+    2023-04-03 19:16:24.061 | INFO     | kliff.loss:minimize:314 - Finish minimization using method: lm.
+    2023-04-03 19:16:24.062 | INFO     | kliff.loss:minimize:312 - Start minimization using method: lm.
+    2023-04-03 19:16:24.062 | INFO     | kliff.loss:_scipy_optimize:431 - Running in serial mode.
+    2023-04-03 19:16:24.105 | INFO     | kliff.loss:minimize:314 - Finish minimization using method: lm.
+    2023-04-03 19:16:24.105 | INFO     | kliff.loss:minimize:312 - Start minimization using method: lm.
+    2023-04-03 19:16:24.106 | INFO     | kliff.loss:_scipy_optimize:431 - Running in serial mode.
+    2023-04-03 19:16:24.147 | INFO     | kliff.loss:minimize:314 - Finish minimization using method: lm.
+    2023-04-03 19:16:24.148 | INFO     | kliff.loss:minimize:312 - Start minimization using method: lm.
+    2023-04-03 19:16:24.148 | INFO     | kliff.loss:_scipy_optimize:431 - Running in serial mode.
+    2023-04-03 19:16:24.188 | INFO     | kliff.loss:minimize:314 - Finish minimization using method: lm.
+    2023-04-03 19:16:24.189 | INFO     | kliff.loss:minimize:312 - Start minimization using method: lm.
+    2023-04-03 19:16:24.189 | INFO     | kliff.loss:_scipy_optimize:431 - Running in serial mode.
+    2023-04-03 19:16:24.229 | INFO     | kliff.loss:minimize:314 - Finish minimization using method: lm.
+    2023-04-03 19:16:24.230 | INFO     | kliff.loss:minimize:312 - Start minimization using method: lm.
+    2023-04-03 19:16:24.231 | INFO     | kliff.loss:_scipy_optimize:431 - Running in serial mode.
+    2023-04-03 19:16:24.274 | INFO     | kliff.loss:minimize:314 - Finish minimization using method: lm.
+    2023-04-03 19:16:24.274 | INFO     | kliff.loss:minimize:312 - Start minimization using method: lm.
+    2023-04-03 19:16:24.275 | INFO     | kliff.loss:_scipy_optimize:431 - Running in serial mode.
+    2023-04-03 19:16:24.316 | INFO     | kliff.loss:minimize:314 - Finish minimization using method: lm.
+    2023-04-03 19:16:24.317 | INFO     | kliff.loss:minimize:312 - Start minimization using method: lm.
+    2023-04-03 19:16:24.317 | INFO     | kliff.loss:_scipy_optimize:431 - Running in serial mode.
+    2023-04-03 19:16:24.357 | INFO     | kliff.loss:minimize:314 - Finish minimization using method: lm.
+    2023-04-03 19:16:24.358 | INFO     | kliff.loss:minimize:312 - Start minimization using method: lm.
+    2023-04-03 19:16:24.358 | INFO     | kliff.loss:_scipy_optimize:431 - Running in serial mode.
+    2023-04-03 19:16:24.401 | INFO     | kliff.loss:minimize:314 - Finish minimization using method: lm.
+    2023-04-03 19:16:24.402 | INFO     | kliff.loss:minimize:312 - Start minimization using method: lm.
+    2023-04-03 19:16:24.402 | INFO     | kliff.loss:_scipy_optimize:431 - Running in serial mode.
+    2023-04-03 19:16:24.444 | INFO     | kliff.loss:minimize:314 - Finish minimization using method: lm.
+    2023-04-03 19:16:24.445 | INFO     | kliff.loss:minimize:312 - Start minimization using method: lm.
+    2023-04-03 19:16:24.445 | INFO     | kliff.loss:_scipy_optimize:431 - Running in serial mode.
+    2023-04-03 19:16:24.490 | INFO     | kliff.loss:minimize:314 - Finish minimization using method: lm.
+    2023-04-03 19:16:24.491 | INFO     | kliff.loss:minimize:312 - Start minimization using method: lm.
+    2023-04-03 19:16:24.491 | INFO     | kliff.loss:_scipy_optimize:431 - Running in serial mode.
+    2023-04-03 19:16:24.534 | INFO     | kliff.loss:minimize:314 - Finish minimization using method: lm.
+    2023-04-03 19:16:24.535 | INFO     | kliff.loss:minimize:312 - Start minimization using method: lm.
+    2023-04-03 19:16:24.535 | INFO     | kliff.loss:_scipy_optimize:431 - Running in serial mode.
+    2023-04-03 19:16:24.581 | INFO     | kliff.loss:minimize:314 - Finish minimization using method: lm.
+    2023-04-03 19:16:24.582 | INFO     | kliff.loss:minimize:312 - Start minimization using method: lm.
+    2023-04-03 19:16:24.583 | INFO     | kliff.loss:_scipy_optimize:431 - Running in serial mode.
+    2023-04-03 19:16:24.625 | INFO     | kliff.loss:minimize:314 - Finish minimization using method: lm.
+    2023-04-03 19:16:24.626 | INFO     | kliff.loss:minimize:312 - Start minimization using method: lm.
+    2023-04-03 19:16:24.627 | INFO     | kliff.loss:_scipy_optimize:431 - Running in serial mode.
+    2023-04-03 19:16:24.668 | INFO     | kliff.loss:minimize:314 - Finish minimization using method: lm.
+    2023-04-03 19:16:24.669 | INFO     | kliff.loss:minimize:312 - Start minimization using method: lm.
+    2023-04-03 19:16:24.669 | INFO     | kliff.loss:_scipy_optimize:431 - Running in serial mode.
+    2023-04-03 19:16:24.714 | INFO     | kliff.loss:minimize:314 - Finish minimization using method: lm.
+    2023-04-03 19:16:24.715 | INFO     | kliff.loss:minimize:312 - Start minimization using method: lm.
+    2023-04-03 19:16:24.716 | INFO     | kliff.loss:_scipy_optimize:431 - Running in serial mode.
+    2023-04-03 19:16:24.753 | INFO     | kliff.loss:minimize:314 - Finish minimization using method: lm.
+    2023-04-03 19:16:24.754 | INFO     | kliff.loss:minimize:312 - Start minimization using method: lm.
+    2023-04-03 19:16:24.754 | INFO     | kliff.loss:_scipy_optimize:431 - Running in serial mode.
+    2023-04-03 19:16:24.790 | INFO     | kliff.loss:minimize:314 - Finish minimization using method: lm.
+    2023-04-03 19:16:24.790 | INFO     | kliff.loss:minimize:312 - Start minimization using method: lm.
+    2023-04-03 19:16:24.791 | INFO     | kliff.loss:_scipy_optimize:431 - Running in serial mode.
+    2023-04-03 19:16:24.833 | INFO     | kliff.loss:minimize:314 - Finish minimization using method: lm.
+    2023-04-03 19:16:24.834 | INFO     | kliff.loss:minimize:312 - Start minimization using method: lm.
+    2023-04-03 19:16:24.834 | INFO     | kliff.loss:_scipy_optimize:431 - Running in serial mode.
+    2023-04-03 19:16:24.874 | INFO     | kliff.loss:minimize:314 - Finish minimization using method: lm.
+    2023-04-03 19:16:24.874 | INFO     | kliff.loss:minimize:312 - Start minimization using method: lm.
+    2023-04-03 19:16:24.875 | INFO     | kliff.loss:_scipy_optimize:431 - Running in serial mode.
+    2023-04-03 19:16:24.916 | INFO     | kliff.loss:minimize:314 - Finish minimization using method: lm.
+    2023-04-03 19:16:24.917 | INFO     | kliff.loss:minimize:312 - Start minimization using method: lm.
+    2023-04-03 19:16:24.918 | INFO     | kliff.loss:_scipy_optimize:431 - Running in serial mode.
+    2023-04-03 19:16:24.958 | INFO     | kliff.loss:minimize:314 - Finish minimization using method: lm.
+    2023-04-03 19:16:24.959 | INFO     | kliff.loss:minimize:312 - Start minimization using method: lm.
+    2023-04-03 19:16:24.959 | INFO     | kliff.loss:_scipy_optimize:431 - Running in serial mode.
+    2023-04-03 19:16:25.006 | INFO     | kliff.loss:minimize:314 - Finish minimization using method: lm.
+    2023-04-03 19:16:25.006 | INFO     | kliff.loss:minimize:312 - Start minimization using method: lm.
+    2023-04-03 19:16:25.008 | INFO     | kliff.loss:_scipy_optimize:431 - Running in serial mode.
+    2023-04-03 19:16:25.050 | INFO     | kliff.loss:minimize:314 - Finish minimization using method: lm.
+    2023-04-03 19:16:25.050 | INFO     | kliff.loss:minimize:312 - Start minimization using method: lm.
+    2023-04-03 19:16:25.051 | INFO     | kliff.loss:_scipy_optimize:431 - Running in serial mode.
+    2023-04-03 19:16:25.092 | INFO     | kliff.loss:minimize:314 - Finish minimization using method: lm.
+    2023-04-03 19:16:25.092 | INFO     | kliff.loss:minimize:312 - Start minimization using method: lm.
+    2023-04-03 19:16:25.093 | INFO     | kliff.loss:_scipy_optimize:431 - Running in serial mode.
+    2023-04-03 19:16:25.138 | INFO     | kliff.loss:minimize:314 - Finish minimization using method: lm.
+    2023-04-03 19:16:25.139 | INFO     | kliff.loss:minimize:312 - Start minimization using method: lm.
+    2023-04-03 19:16:25.139 | INFO     | kliff.loss:_scipy_optimize:431 - Running in serial mode.
+    2023-04-03 19:16:25.181 | INFO     | kliff.loss:minimize:314 - Finish minimization using method: lm.
+    2023-04-03 19:16:25.182 | INFO     | kliff.loss:minimize:312 - Start minimization using method: lm.
+    2023-04-03 19:16:25.182 | INFO     | kliff.loss:_scipy_optimize:431 - Running in serial mode.
+    2023-04-03 19:16:25.221 | INFO     | kliff.loss:minimize:314 - Finish minimization using method: lm.
+    2023-04-03 19:16:25.222 | INFO     | kliff.loss:minimize:312 - Start minimization using method: lm.
+    2023-04-03 19:16:25.222 | INFO     | kliff.loss:_scipy_optimize:431 - Running in serial mode.
+    2023-04-03 19:16:25.266 | INFO     | kliff.loss:minimize:314 - Finish minimization using method: lm.
+    2023-04-03 19:16:25.267 | INFO     | kliff.loss:minimize:312 - Start minimization using method: lm.
+    2023-04-03 19:16:25.267 | INFO     | kliff.loss:_scipy_optimize:431 - Running in serial mode.
+    2023-04-03 19:16:25.308 | INFO     | kliff.loss:minimize:314 - Finish minimization using method: lm.
+    2023-04-03 19:16:25.308 | INFO     | kliff.loss:minimize:312 - Start minimization using method: lm.
+    2023-04-03 19:16:25.309 | INFO     | kliff.loss:_scipy_optimize:431 - Running in serial mode.
+    2023-04-03 19:16:25.351 | INFO     | kliff.loss:minimize:314 - Finish minimization using method: lm.
+    2023-04-03 19:16:25.352 | INFO     | kliff.loss:minimize:312 - Start minimization using method: lm.
+    2023-04-03 19:16:25.353 | INFO     | kliff.loss:_scipy_optimize:431 - Running in serial mode.
+    2023-04-03 19:16:25.395 | INFO     | kliff.loss:minimize:314 - Finish minimization using method: lm.
+    2023-04-03 19:16:25.395 | INFO     | kliff.loss:minimize:312 - Start minimization using method: lm.
+    2023-04-03 19:16:25.396 | INFO     | kliff.loss:_scipy_optimize:431 - Running in serial mode.
+    2023-04-03 19:16:25.437 | INFO     | kliff.loss:minimize:314 - Finish minimization using method: lm.
+    2023-04-03 19:16:25.437 | INFO     | kliff.loss:minimize:312 - Start minimization using method: lm.
+    2023-04-03 19:16:25.438 | INFO     | kliff.loss:_scipy_optimize:431 - Running in serial mode.
+    2023-04-03 19:16:25.479 | INFO     | kliff.loss:minimize:314 - Finish minimization using method: lm.
+    2023-04-03 19:16:25.480 | INFO     | kliff.loss:minimize:312 - Start minimization using method: lm.
+    2023-04-03 19:16:25.480 | INFO     | kliff.loss:_scipy_optimize:431 - Running in serial mode.
+    2023-04-03 19:16:25.520 | INFO     | kliff.loss:minimize:314 - Finish minimization using method: lm.
+    2023-04-03 19:16:25.521 | INFO     | kliff.loss:minimize:312 - Start minimization using method: lm.
+    2023-04-03 19:16:25.521 | INFO     | kliff.loss:_scipy_optimize:431 - Running in serial mode.
+    2023-04-03 19:16:25.561 | INFO     | kliff.loss:minimize:314 - Finish minimization using method: lm.
+    2023-04-03 19:16:25.562 | INFO     | kliff.loss:minimize:312 - Start minimization using method: lm.
+    2023-04-03 19:16:25.564 | INFO     | kliff.loss:_scipy_optimize:431 - Running in serial mode.
+    2023-04-03 19:16:25.610 | INFO     | kliff.loss:minimize:314 - Finish minimization using method: lm.
+    2023-04-03 19:16:25.611 | INFO     | kliff.loss:minimize:312 - Start minimization using method: lm.
+    2023-04-03 19:16:25.612 | INFO     | kliff.loss:_scipy_optimize:431 - Running in serial mode.
+    2023-04-03 19:16:25.655 | INFO     | kliff.loss:minimize:314 - Finish minimization using method: lm.
+    2023-04-03 19:16:25.655 | INFO     | kliff.loss:minimize:312 - Start minimization using method: lm.
+    2023-04-03 19:16:25.656 | INFO     | kliff.loss:_scipy_optimize:431 - Running in serial mode.
+    2023-04-03 19:16:25.697 | INFO     | kliff.loss:minimize:314 - Finish minimization using method: lm.
+    2023-04-03 19:16:25.697 | INFO     | kliff.loss:minimize:312 - Start minimization using method: lm.
+    2023-04-03 19:16:25.698 | INFO     | kliff.loss:_scipy_optimize:431 - Running in serial mode.
+    2023-04-03 19:16:25.738 | INFO     | kliff.loss:minimize:314 - Finish minimization using method: lm.
+    2023-04-03 19:16:25.739 | INFO     | kliff.loss:minimize:312 - Start minimization using method: lm.
+    2023-04-03 19:16:25.739 | INFO     | kliff.loss:_scipy_optimize:431 - Running in serial mode.
+    2023-04-03 19:16:25.782 | INFO     | kliff.loss:minimize:314 - Finish minimization using method: lm.
+    2023-04-03 19:16:25.782 | INFO     | kliff.loss:minimize:312 - Start minimization using method: lm.
+    2023-04-03 19:16:25.783 | INFO     | kliff.loss:_scipy_optimize:431 - Running in serial mode.
+    2023-04-03 19:16:25.825 | INFO     | kliff.loss:minimize:314 - Finish minimization using method: lm.
+    2023-04-03 19:16:25.826 | INFO     | kliff.loss:minimize:312 - Start minimization using method: lm.
+    2023-04-03 19:16:25.826 | INFO     | kliff.loss:_scipy_optimize:431 - Running in serial mode.
+    2023-04-03 19:16:25.867 | INFO     | kliff.loss:minimize:314 - Finish minimization using method: lm.
+    2023-04-03 19:16:25.867 | INFO     | kliff.loss:minimize:312 - Start minimization using method: lm.
+    2023-04-03 19:16:25.868 | INFO     | kliff.loss:_scipy_optimize:431 - Running in serial mode.
+    2023-04-03 19:16:25.909 | INFO     | kliff.loss:minimize:314 - Finish minimization using method: lm.
+    2023-04-03 19:16:25.909 | INFO     | kliff.loss:minimize:312 - Start minimization using method: lm.
+    2023-04-03 19:16:25.910 | INFO     | kliff.loss:_scipy_optimize:431 - Running in serial mode.
+    2023-04-03 19:16:25.950 | INFO     | kliff.loss:minimize:314 - Finish minimization using method: lm.
+    2023-04-03 19:16:25.951 | INFO     | kliff.loss:minimize:312 - Start minimization using method: lm.
+    2023-04-03 19:16:25.951 | INFO     | kliff.loss:_scipy_optimize:431 - Running in serial mode.
+    2023-04-03 19:16:25.997 | INFO     | kliff.loss:minimize:314 - Finish minimization using method: lm.
+    2023-04-03 19:16:25.998 | INFO     | kliff.loss:minimize:312 - Start minimization using method: lm.
+    2023-04-03 19:16:25.999 | INFO     | kliff.loss:_scipy_optimize:431 - Running in serial mode.
+    2023-04-03 19:16:26.042 | INFO     | kliff.loss:minimize:314 - Finish minimization using method: lm.
+    2023-04-03 19:16:26.043 | INFO     | kliff.loss:minimize:312 - Start minimization using method: lm.
+    2023-04-03 19:16:26.043 | INFO     | kliff.loss:_scipy_optimize:431 - Running in serial mode.
+    2023-04-03 19:16:26.084 | INFO     | kliff.loss:minimize:314 - Finish minimization using method: lm.
+    2023-04-03 19:16:26.085 | INFO     | kliff.loss:minimize:312 - Start minimization using method: lm.
+    2023-04-03 19:16:26.086 | INFO     | kliff.loss:_scipy_optimize:431 - Running in serial mode.
+    2023-04-03 19:16:26.128 | INFO     | kliff.loss:minimize:314 - Finish minimization using method: lm.
+    2023-04-03 19:16:26.128 | INFO     | kliff.loss:minimize:312 - Start minimization using method: lm.
+    2023-04-03 19:16:26.129 | INFO     | kliff.loss:_scipy_optimize:431 - Running in serial mode.
+    2023-04-03 19:16:26.174 | INFO     | kliff.loss:minimize:314 - Finish minimization using method: lm.
+    2023-04-03 19:16:26.175 | INFO     | kliff.loss:minimize:312 - Start minimization using method: lm.
+    2023-04-03 19:16:26.176 | INFO     | kliff.loss:_scipy_optimize:431 - Running in serial mode.
+    2023-04-03 19:16:26.212 | INFO     | kliff.loss:minimize:314 - Finish minimization using method: lm.
+    2023-04-03 19:16:26.213 | INFO     | kliff.loss:minimize:312 - Start minimization using method: lm.
+    2023-04-03 19:16:26.213 | INFO     | kliff.loss:_scipy_optimize:431 - Running in serial mode.
+    2023-04-03 19:16:26.254 | INFO     | kliff.loss:minimize:314 - Finish minimization using method: lm.
+    2023-04-03 19:16:26.255 | INFO     | kliff.loss:minimize:312 - Start minimization using method: lm.
+    2023-04-03 19:16:26.255 | INFO     | kliff.loss:_scipy_optimize:431 - Running in serial mode.
+    2023-04-03 19:16:26.296 | INFO     | kliff.loss:minimize:314 - Finish minimization using method: lm.
+    2023-04-03 19:16:26.297 | INFO     | kliff.loss:minimize:312 - Start minimization using method: lm.
+    2023-04-03 19:16:26.297 | INFO     | kliff.loss:_scipy_optimize:431 - Running in serial mode.
+    2023-04-03 19:16:26.338 | INFO     | kliff.loss:minimize:314 - Finish minimization using method: lm.
+    2023-04-03 19:16:26.338 | INFO     | kliff.loss:minimize:312 - Start minimization using method: lm.
+    2023-04-03 19:16:26.339 | INFO     | kliff.loss:_scipy_optimize:431 - Running in serial mode.
+    2023-04-03 19:16:26.380 | INFO     | kliff.loss:minimize:314 - Finish minimization using method: lm.
+    2023-04-03 19:16:26.381 | INFO     | kliff.loss:minimize:312 - Start minimization using method: lm.
+    2023-04-03 19:16:26.382 | INFO     | kliff.loss:_scipy_optimize:431 - Running in serial mode.
+    2023-04-03 19:16:26.424 | INFO     | kliff.loss:minimize:314 - Finish minimization using method: lm.
+    2023-04-03 19:16:26.426 | INFO     | kliff.loss:minimize:312 - Start minimization using method: lm.
+    2023-04-03 19:16:26.427 | INFO     | kliff.loss:_scipy_optimize:431 - Running in serial mode.
+    2023-04-03 19:16:26.466 | INFO     | kliff.loss:minimize:314 - Finish minimization using method: lm.
+    2023-04-03 19:16:26.466 | INFO     | kliff.loss:minimize:312 - Start minimization using method: lm.
+    2023-04-03 19:16:26.467 | INFO     | kliff.loss:_scipy_optimize:431 - Running in serial mode.
+    2023-04-03 19:16:26.506 | INFO     | kliff.loss:minimize:314 - Finish minimization using method: lm.
+    2023-04-03 19:16:26.506 | INFO     | kliff.loss:minimize:312 - Start minimization using method: lm.
+    2023-04-03 19:16:26.507 | INFO     | kliff.loss:_scipy_optimize:431 - Running in serial mode.
+    2023-04-03 19:16:26.545 | INFO     | kliff.loss:minimize:314 - Finish minimization using method: lm.
+    2023-04-03 19:16:26.546 | INFO     | kliff.loss:minimize:312 - Start minimization using method: lm.
+    2023-04-03 19:16:26.546 | INFO     | kliff.loss:_scipy_optimize:431 - Running in serial mode.
+    2023-04-03 19:16:26.589 | INFO     | kliff.loss:minimize:314 - Finish minimization using method: lm.
+    2023-04-03 19:16:26.590 | INFO     | kliff.loss:minimize:312 - Start minimization using method: lm.
+    2023-04-03 19:16:26.590 | INFO     | kliff.loss:_scipy_optimize:431 - Running in serial mode.
+    2023-04-03 19:16:26.631 | INFO     | kliff.loss:minimize:314 - Finish minimization using method: lm.
+    2023-04-03 19:16:26.631 | INFO     | kliff.loss:minimize:312 - Start minimization using method: lm.
+    2023-04-03 19:16:26.632 | INFO     | kliff.loss:_scipy_optimize:431 - Running in serial mode.
+    2023-04-03 19:16:26.673 | INFO     | kliff.loss:minimize:314 - Finish minimization using method: lm.
+    2023-04-03 19:16:26.674 | INFO     | kliff.loss:minimize:312 - Start minimization using method: lm.
+    2023-04-03 19:16:26.674 | INFO     | kliff.loss:_scipy_optimize:431 - Running in serial mode.
+    2023-04-03 19:16:26.717 | INFO     | kliff.loss:minimize:314 - Finish minimization using method: lm.
+    2023-04-03 19:16:26.717 | INFO     | kliff.loss:minimize:312 - Start minimization using method: lm.
+    2023-04-03 19:16:26.718 | INFO     | kliff.loss:_scipy_optimize:431 - Running in serial mode.
+    2023-04-03 19:16:26.757 | INFO     | kliff.loss:minimize:314 - Finish minimization using method: lm.
+    2023-04-03 19:16:26.758 | INFO     | kliff.loss:minimize:312 - Start minimization using method: lm.
+    2023-04-03 19:16:26.758 | INFO     | kliff.loss:_scipy_optimize:431 - Running in serial mode.
+    2023-04-03 19:16:26.802 | INFO     | kliff.loss:minimize:314 - Finish minimization using method: lm.
+    2023-04-03 19:16:26.803 | INFO     | kliff.loss:minimize:312 - Start minimization using method: lm.
+    2023-04-03 19:16:26.803 | INFO     | kliff.loss:_scipy_optimize:431 - Running in serial mode.
+    2023-04-03 19:16:26.843 | INFO     | kliff.loss:minimize:314 - Finish minimization using method: lm.
+    2023-04-03 19:16:26.844 | INFO     | kliff.loss:minimize:312 - Start minimization using method: lm.
+    2023-04-03 19:16:26.844 | INFO     | kliff.loss:_scipy_optimize:431 - Running in serial mode.
+    2023-04-03 19:16:26.886 | INFO     | kliff.loss:minimize:314 - Finish minimization using method: lm.
+    2023-04-03 19:16:26.887 | INFO     | kliff.loss:minimize:312 - Start minimization using method: lm.
+    2023-04-03 19:16:26.887 | INFO     | kliff.loss:_scipy_optimize:431 - Running in serial mode.
+    2023-04-03 19:16:26.929 | INFO     | kliff.loss:minimize:314 - Finish minimization using method: lm.
+    2023-04-03 19:16:26.930 | INFO     | kliff.loss:minimize:312 - Start minimization using method: lm.
+    2023-04-03 19:16:26.930 | INFO     | kliff.loss:_scipy_optimize:431 - Running in serial mode.
+    2023-04-03 19:16:26.972 | INFO     | kliff.loss:minimize:314 - Finish minimization using method: lm.
+    2023-04-03 19:16:26.972 | INFO     | kliff.loss:minimize:312 - Start minimization using method: lm.
+    2023-04-03 19:16:26.972 | INFO     | kliff.loss:_scipy_optimize:431 - Running in serial mode.
+    2023-04-03 19:16:27.014 | INFO     | kliff.loss:minimize:314 - Finish minimization using method: lm.
+    2023-04-03 19:16:27.015 | INFO     | kliff.loss:minimize:312 - Start minimization using method: lm.
+    2023-04-03 19:16:27.016 | INFO     | kliff.loss:_scipy_optimize:431 - Running in serial mode.
+    2023-04-03 19:16:27.060 | INFO     | kliff.loss:minimize:314 - Finish minimization using method: lm.
+    2023-04-03 19:16:27.061 | INFO     | kliff.loss:minimize:312 - Start minimization using method: lm.
+    2023-04-03 19:16:27.061 | INFO     | kliff.loss:_scipy_optimize:431 - Running in serial mode.
+    2023-04-03 19:16:27.105 | INFO     | kliff.loss:minimize:314 - Finish minimization using method: lm.
+    2023-04-03 19:16:27.105 | INFO     | kliff.loss:minimize:312 - Start minimization using method: lm.
+    2023-04-03 19:16:27.106 | INFO     | kliff.loss:_scipy_optimize:431 - Running in serial mode.
+    2023-04-03 19:16:27.149 | INFO     | kliff.loss:minimize:314 - Finish minimization using method: lm.
+    2023-04-03 19:16:27.150 | INFO     | kliff.loss:minimize:312 - Start minimization using method: lm.
+    2023-04-03 19:16:27.150 | INFO     | kliff.loss:_scipy_optimize:431 - Running in serial mode.
+    2023-04-03 19:16:27.193 | INFO     | kliff.loss:minimize:314 - Finish minimization using method: lm.
+    2023-04-03 19:16:27.194 | INFO     | kliff.loss:minimize:312 - Start minimization using method: lm.
+    2023-04-03 19:16:27.195 | INFO     | kliff.loss:_scipy_optimize:431 - Running in serial mode.
+    2023-04-03 19:16:27.236 | INFO     | kliff.loss:minimize:314 - Finish minimization using method: lm.
+    2023-04-03 19:16:27.236 | INFO     | kliff.loss:minimize:312 - Start minimization using method: lm.
+    2023-04-03 19:16:27.237 | INFO     | kliff.loss:_scipy_optimize:431 - Running in serial mode.
+    2023-04-03 19:16:27.279 | INFO     | kliff.loss:minimize:314 - Finish minimization using method: lm.
+    2023-04-03 19:16:27.279 | INFO     | kliff.loss:minimize:312 - Start minimization using method: lm.
+    2023-04-03 19:16:27.279 | INFO     | kliff.loss:_scipy_optimize:431 - Running in serial mode.
+    2023-04-03 19:16:27.321 | INFO     | kliff.loss:minimize:314 - Finish minimization using method: lm.
+    2023-04-03 19:16:27.321 | INFO     | kliff.loss:minimize:312 - Start minimization using method: lm.
+    2023-04-03 19:16:27.322 | INFO     | kliff.loss:_scipy_optimize:431 - Running in serial mode.
+    2023-04-03 19:16:27.362 | INFO     | kliff.loss:minimize:314 - Finish minimization using method: lm.
+    2023-04-03 19:16:27.363 | INFO     | kliff.loss:minimize:312 - Start minimization using method: lm.
+    2023-04-03 19:16:27.365 | INFO     | kliff.loss:_scipy_optimize:431 - Running in serial mode.
+    2023-04-03 19:16:27.405 | INFO     | kliff.loss:minimize:314 - Finish minimization using method: lm.
+    2023-04-03 19:16:27.406 | INFO     | kliff.loss:minimize:312 - Start minimization using method: lm.
+    2023-04-03 19:16:27.406 | INFO     | kliff.loss:_scipy_optimize:431 - Running in serial mode.
+    2023-04-03 19:16:27.446 | INFO     | kliff.loss:minimize:314 - Finish minimization using method: lm.
+    2023-04-03 19:16:27.447 | INFO     | kliff.loss:minimize:312 - Start minimization using method: lm.
+    2023-04-03 19:16:27.448 | INFO     | kliff.loss:_scipy_optimize:431 - Running in serial mode.
+    2023-04-03 19:16:27.488 | INFO     | kliff.loss:minimize:314 - Finish minimization using method: lm.
+    2023-04-03 19:16:27.488 | INFO     | kliff.loss:minimize:312 - Start minimization using method: lm.
+    2023-04-03 19:16:27.489 | INFO     | kliff.loss:_scipy_optimize:431 - Running in serial mode.
+    2023-04-03 19:16:27.532 | INFO     | kliff.loss:minimize:314 - Finish minimization using method: lm.
+    2023-04-03 19:16:27.533 | INFO     | kliff.loss:minimize:312 - Start minimization using method: lm.
+    2023-04-03 19:16:27.533 | INFO     | kliff.loss:_scipy_optimize:431 - Running in serial mode.
+    2023-04-03 19:16:27.577 | INFO     | kliff.loss:minimize:314 - Finish minimization using method: lm.
+    2023-04-03 19:16:27.577 | INFO     | kliff.loss:minimize:312 - Start minimization using method: lm.
+    2023-04-03 19:16:27.578 | INFO     | kliff.loss:_scipy_optimize:431 - Running in serial mode.
+    2023-04-03 19:16:27.619 | INFO     | kliff.loss:minimize:314 - Finish minimization using method: lm.
+    2023-04-03 19:16:27.620 | INFO     | kliff.loss:minimize:312 - Start minimization using method: lm.
+    2023-04-03 19:16:27.620 | INFO     | kliff.loss:_scipy_optimize:431 - Running in serial mode.
+    2023-04-03 19:16:27.662 | INFO     | kliff.loss:minimize:314 - Finish minimization using method: lm.
+    2023-04-03 19:16:27.662 | INFO     | kliff.loss:minimize:312 - Start minimization using method: lm.
+    2023-04-03 19:16:27.663 | INFO     | kliff.loss:_scipy_optimize:431 - Running in serial mode.
+    2023-04-03 19:16:27.704 | INFO     | kliff.loss:minimize:314 - Finish minimization using method: lm.
+    2023-04-03 19:16:27.705 | INFO     | kliff.loss:minimize:312 - Start minimization using method: lm.
+    2023-04-03 19:16:27.705 | INFO     | kliff.loss:_scipy_optimize:431 - Running in serial mode.
+    2023-04-03 19:16:27.750 | INFO     | kliff.loss:minimize:314 - Finish minimization using method: lm.
+    2023-04-03 19:16:27.750 | INFO     | kliff.loss:minimize:312 - Start minimization using method: lm.
+    2023-04-03 19:16:27.751 | INFO     | kliff.loss:_scipy_optimize:431 - Running in serial mode.
+    2023-04-03 19:16:27.793 | INFO     | kliff.loss:minimize:314 - Finish minimization using method: lm.
+    2023-04-03 19:16:27.794 | INFO     | kliff.loss:minimize:312 - Start minimization using method: lm.
+    2023-04-03 19:16:27.794 | INFO     | kliff.loss:_scipy_optimize:431 - Running in serial mode.
+    2023-04-03 19:16:27.836 | INFO     | kliff.loss:minimize:314 - Finish minimization using method: lm.
+
+    array([[15.28616743, 45.53413524],
+           [15.28539627, 45.53724361],
+           [15.28585555, 45.53650355],
+           [15.28603291, 45.53612439],
+           [15.28571892, 45.53489973],
+           [15.28539627, 45.53724361],
+           [15.28570725, 45.5312858 ],
+           [15.28584281, 45.532998  ],
+           [15.28616881, 45.53755689],
+           [15.28559134, 45.53694894],
+           [15.28546948, 45.53511245],
+           [15.28515489, 45.53371441],
+           [15.28584281, 45.532998  ],
+           [15.28559028, 45.53312488],
+           [15.28578432, 45.53278061],
+           [15.28630964, 45.5355221 ],
+           [15.28559028, 45.53312488],
+           [15.28559028, 45.53312488],
+           [15.28571892, 45.53489973],
+           [15.28559028, 45.53312488],
+           [15.28526585, 45.53538115],
+           [15.28590466, 45.53453664],
+           [15.28571892, 45.53489973],
+           [15.28571892, 45.53489973],
+           [15.28551882, 45.5352354 ],
+           [15.28584281, 45.532998  ],
+           [15.28603291, 45.53612439],
+           [15.28538128, 45.53343707],
+           [15.2864703 , 45.53508198],
+           [15.28570725, 45.5312858 ],
+           [15.28559028, 45.53312488],
+           [15.28572139, 45.53860548],
+           [15.28515489, 45.53371441],
+           [15.28616881, 45.53755689],
+           [15.28551882, 45.5352354 ],
+           [15.28526585, 45.53538115],
+           [15.28603291, 45.53612439],
+           [15.28570725, 45.5312858 ],
+           [15.28633634, 45.53370352],
+           [15.28571892, 45.53489973],
+           [15.28559134, 45.53694894],
+           [15.28572139, 45.53860548],
+           [15.28590466, 45.53453664],
+           [15.28584281, 45.532998  ],
+           [15.28526585, 45.53538115],
+           [15.28603291, 45.53612439],
+           [15.28609088, 45.5304963 ],
+           [15.28493205, 45.53767105],
+           [15.28546948, 45.53511245],
+           [15.28616743, 45.53413524],
+           [15.28633634, 45.53370352],
+           [15.28526585, 45.53538115],
+           [15.28571892, 45.53489973],
+           [15.28515489, 45.53371441],
+           [15.28578432, 45.53278061],
+           [15.28549314, 45.53163801],
+           [15.28572139, 45.53860548],
+           [15.28559134, 45.53694894],
+           [15.28526585, 45.53538115],
+           [15.28603291, 45.53612439],
+           [15.28590466, 45.53453664],
+           [15.28570725, 45.5312858 ],
+           [15.28571892, 45.53489973],
+           [15.28603275, 45.53259999],
+           [15.28598621, 45.53454737],
+           [15.28551882, 45.5352354 ],
+           [15.28590466, 45.53453664],
+           [15.28598621, 45.53454737],
+           [15.28549314, 45.53163801],
+           [15.28598621, 45.53454737],
+           [15.28539627, 45.53724361],
+           [15.28571892, 45.53489973],
+           [15.28585555, 45.53650355],
+           [15.2864703 , 45.53508198],
+           [15.28559028, 45.53312488],
+           [15.28603275, 45.53259999],
+           [15.28526585, 45.53538115],
+           [15.28515489, 45.53371441],
+           [15.28551882, 45.5352354 ],
+           [15.28590466, 45.53453664],
+           [15.28590466, 45.53453664],
+           [15.28603275, 45.53259999],
+           [15.28571892, 45.53489973],
+           [15.28504553, 45.53561156],
+           [15.28559028, 45.53312488],
+           [15.28559134, 45.53694894],
+           [15.28616743, 45.53413524],
+           [15.28481349, 45.53991774],
+           [15.28538128, 45.53343707],
+           [15.28620979, 45.53217886],
+           [15.28538128, 45.53343707],
+           [15.28570725, 45.5312858 ],
+           [15.28603291, 45.53612439],
+           [15.28546948, 45.53511245],
+           [15.28559134, 45.53694894],
+           [15.2859058 , 45.53090459],
+           [15.28571892, 45.53489973],
+           [15.28590466, 45.53453664],
+           [15.28546948, 45.53511245],
+           [15.28493205, 45.53767105]])
+
+
+
+.. GENERATED FROM PYTHON SOURCE LINES 109-112
+
+The resulting parameter ensemble can be accessed in `BS.samples` as a `np.ndarray`.
+Then, we can plot the distribution of the parameters, as an example, or propagate the
+error to the target quantities we want to study.
+
+.. GENERATED FROM PYTHON SOURCE LINES 112-122
+
+.. code-block:: default
+
+
+
+    # Plot the distribution of the parameters
+    plt.figure()
+    plt.plot(*(BS.samples.T), ".", alpha=0.5)
+    param_names = list(opt_params.keys())
+    plt.xlabel(param_names[0])
+    plt.ylabel(param_names[1])
+    plt.show()
+
+
+
+
+.. image-sg:: /auto_examples/images/sphx_glr_example_uq_bootstrap_001.png
+   :alt: example uq bootstrap
+   :srcset: /auto_examples/images/sphx_glr_example_uq_bootstrap_001.png
+   :class: sphx-glr-single-img
+
+
+
+
+
+.. GENERATED FROM PYTHON SOURCE LINES 123-124
+
+.. _OpenKIM: https://openkim.org
+
+
+.. rst-class:: sphx-glr-timing
+
+   **Total running time of the script:** ( 0 minutes  5.912 seconds)
+
+
+.. _sphx_glr_download_auto_examples_example_uq_bootstrap.py:
+
+.. only:: html
+
+  .. container:: sphx-glr-footer sphx-glr-footer-example
+
+
+
+
+    .. container:: sphx-glr-download sphx-glr-download-python
+
+      :download:`Download Python source code: example_uq_bootstrap.py <example_uq_bootstrap.py>`
+
+    .. container:: sphx-glr-download sphx-glr-download-jupyter
+
+      :download:`Download Jupyter notebook: example_uq_bootstrap.ipynb <example_uq_bootstrap.ipynb>`
+
+
+.. only:: html
+
+ .. rst-class:: sphx-glr-signature
+
+    `Gallery generated by Sphinx-Gallery <https://sphinx-gallery.github.io>`_
diff --git a/docs/source/auto_examples/example_uq_bootstrap_codeobj.pickle b/docs/source/auto_examples/example_uq_bootstrap_codeobj.pickle
new file mode 100644
index 00000000..672a4f67
Binary files /dev/null and b/docs/source/auto_examples/example_uq_bootstrap_codeobj.pickle differ
diff --git a/docs/source/auto_examples/images/sphx_glr_example_uq_bootstrap_001.png b/docs/source/auto_examples/images/sphx_glr_example_uq_bootstrap_001.png
new file mode 100644
index 00000000..6ec7608b
Binary files /dev/null and b/docs/source/auto_examples/images/sphx_glr_example_uq_bootstrap_001.png differ
diff --git a/docs/source/auto_examples/images/thumb/sphx_glr_example_uq_bootstrap_thumb.png b/docs/source/auto_examples/images/thumb/sphx_glr_example_uq_bootstrap_thumb.png
new file mode 100644
index 00000000..8a9b9499
Binary files /dev/null and b/docs/source/auto_examples/images/thumb/sphx_glr_example_uq_bootstrap_thumb.png differ
diff --git a/docs/source/auto_examples/index.rst b/docs/source/auto_examples/index.rst
index baebcafb..6c9cfad8 100644
--- a/docs/source/auto_examples/index.rst
+++ b/docs/source/auto_examples/index.rst
@@ -37,7 +37,7 @@ More examples can be found at `<https://github.com/openkim/kliff/tree/master/exa
 .. only:: html
 
   .. image:: /auto_examples/images/thumb/sphx_glr_example_linear_regression_thumb.png
-    :alt: Train a linear regression potential
+    :alt:
 
   :ref:`sphx_glr_auto_examples_example_linear_regression.py`
 
@@ -54,7 +54,7 @@ More examples can be found at `<https://github.com/openkim/kliff/tree/master/exa
 .. only:: html
 
   .. image:: /auto_examples/images/thumb/sphx_glr_example_lennard_jones_thumb.png
-    :alt: Train a Lennard-Jones potential
+    :alt:
 
   :ref:`sphx_glr_auto_examples_example_lennard_jones.py`
 
@@ -71,7 +71,7 @@ More examples can be found at `<https://github.com/openkim/kliff/tree/master/exa
 .. only:: html
 
   .. image:: /auto_examples/images/thumb/sphx_glr_example_kim_SW_Si_thumb.png
-    :alt: Train a Stillinger-Weber potential
+    :alt:
 
   :ref:`sphx_glr_auto_examples_example_kim_SW_Si.py`
 
@@ -88,7 +88,7 @@ More examples can be found at `<https://github.com/openkim/kliff/tree/master/exa
 .. only:: html
 
   .. image:: /auto_examples/images/thumb/sphx_glr_example_parameter_transform_thumb.png
-    :alt: Parameter transformation for the Stillinger-Weber potential
+    :alt:
 
   :ref:`sphx_glr_auto_examples_example_parameter_transform.py`
 
@@ -105,7 +105,7 @@ More examples can be found at `<https://github.com/openkim/kliff/tree/master/exa
 .. only:: html
 
   .. image:: /auto_examples/images/thumb/sphx_glr_example_nn_Si_thumb.png
-    :alt: Train a neural network potential
+    :alt:
 
   :ref:`sphx_glr_auto_examples_example_nn_Si.py`
 
@@ -115,6 +115,23 @@ More examples can be found at `<https://github.com/openkim/kliff/tree/master/exa
     </div>
 
 
+.. raw:: html
+
+    <div class="sphx-glr-thumbcontainer" tooltip="In this example, we demonstrate how to perform uncertainty quantification (UQ) using bootstrap ...">
+
+.. only:: html
+
+  .. image:: /auto_examples/images/thumb/sphx_glr_example_uq_bootstrap_thumb.png
+    :alt:
+
+  :ref:`sphx_glr_auto_examples_example_uq_bootstrap.py`
+
+.. raw:: html
+
+      <div class="sphx-glr-thumbnail-title">Bootstrapping</div>
+    </div>
+
+
 .. raw:: html
 
     <div class="sphx-glr-thumbcontainer" tooltip="In this tutorial, we train a neural network (NN) potential for a system containing two species:...">
@@ -122,7 +139,7 @@ More examples can be found at `<https://github.com/openkim/kliff/tree/master/exa
 .. only:: html
 
   .. image:: /auto_examples/images/thumb/sphx_glr_example_nn_SiC_thumb.png
-    :alt: Train a neural network potential for SiC
+    :alt:
 
   :ref:`sphx_glr_auto_examples_example_nn_SiC.py`
 
@@ -139,7 +156,7 @@ More examples can be found at `<https://github.com/openkim/kliff/tree/master/exa
 .. only:: html
 
   .. image:: /auto_examples/images/thumb/sphx_glr_example_uq_mcmc_thumb.png
-    :alt: MCMC sampling
+    :alt:
 
   :ref:`sphx_glr_auto_examples_example_uq_mcmc.py`
 
@@ -162,6 +179,7 @@ More examples can be found at `<https://github.com/openkim/kliff/tree/master/exa
    /auto_examples/example_kim_SW_Si
    /auto_examples/example_parameter_transform
    /auto_examples/example_nn_Si
+   /auto_examples/example_uq_bootstrap
    /auto_examples/example_nn_SiC
    /auto_examples/example_uq_mcmc
 
diff --git a/docs/source/auto_examples/sg_execution_times.rst b/docs/source/auto_examples/sg_execution_times.rst
index f4c78eeb..5cce2dac 100644
--- a/docs/source/auto_examples/sg_execution_times.rst
+++ b/docs/source/auto_examples/sg_execution_times.rst
@@ -3,12 +3,15 @@
 
 .. _sphx_glr_auto_examples_sg_execution_times:
 
+
 Computation times
 =================
-**00:36.524** total execution time for **auto_examples** files:
+**00:05.912** total execution time for **auto_examples** files:
 
 +---------------------------------------------------------------------------------------------------+-----------+--------+
-| :ref:`sphx_glr_auto_examples_example_kim_SW_Si.py` (``example_kim_SW_Si.py``)                     | 00:36.524 | 0.0 MB |
+| :ref:`sphx_glr_auto_examples_example_uq_bootstrap.py` (``example_uq_bootstrap.py``)               | 00:05.912 | 0.0 MB |
++---------------------------------------------------------------------------------------------------+-----------+--------+
+| :ref:`sphx_glr_auto_examples_example_kim_SW_Si.py` (``example_kim_SW_Si.py``)                     | 00:00.000 | 0.0 MB |
 +---------------------------------------------------------------------------------------------------+-----------+--------+
 | :ref:`sphx_glr_auto_examples_example_lennard_jones.py` (``example_lennard_jones.py``)             | 00:00.000 | 0.0 MB |
 +---------------------------------------------------------------------------------------------------+-----------+--------+
diff --git a/docs/source/contributing_guide.rst b/docs/source/contributing_guide.rst
index f2a559fc..d7d5d955 100644
--- a/docs/source/contributing_guide.rst
+++ b/docs/source/contributing_guide.rst
@@ -9,7 +9,7 @@ Code style
 
   .. code-block:: bash
 
-    pre-commit run pre-commit run --all-files --show-diff-on-failure
+    pre-commit run --all-files --show-diff-on-failure
 
 - The docstring of **KLIFF** follows the `Google` style, which can be found at googledoc_.
 
diff --git a/docs/source/img/bootstrap.png b/docs/source/img/bootstrap.png
new file mode 100755
index 00000000..b9259dff
Binary files /dev/null and b/docs/source/img/bootstrap.png differ
diff --git a/docs/source/img/uq_cartoon.png b/docs/source/img/uq_cartoon.png
new file mode 100755
index 00000000..d9ea35da
Binary files /dev/null and b/docs/source/img/uq_cartoon.png differ
diff --git a/docs/source/modules/uq.rst b/docs/source/modules/uq.rst
index bced7e74..472435aa 100644
--- a/docs/source/modules/uq.rst
+++ b/docs/source/modules/uq.rst
@@ -5,18 +5,34 @@ Uncertainty Quantification (UQ)
 ===============================
 
 Uncertainty quantification (UQ) is an emerging field in applied mathematics that aims to
-quantify uncertainties in mathematical models as a result from error propagation in the
+quantify uncertainties in mathematical models as a result of error propagation in the
 modeling process. This is especially important since we use the model, i.e., the
 potential, to predict material properties that are not used in the
 training process. Thus, UQ process is especially important to assess the reliability of
 these out-of-sample predictions.
 
-In UQ process, we first quantify the uncertainty of the model parameters. Having found
-the parametric uncertainty, then we can propagate the uncertainty of the parameters and
-get the uncertainty of the material properties of interest, e.g., by evaluating the
-ensemble that is obtained from sampling the distribution of the parameters. As the first
-uncertainty propagation is more involved, KLIFF implements tools to quantify the
-uncertainty of the parameters.
+.. figure:: ../img/uq_cartoon.png
+   :alt: Illustration of general uncertainty quantification process.
+   :align: left
+
+   Illustration of general uncertainty quantification process. The error from the
+   training data (represented by a dashed ellipse on the bottom left plot) is
+   propagated to the uncertainty of the model's parameters (the ellipse in the middle
+   plot). Then, the uncertainty of the parameters is propagated forward further to
+   the uncertainty of the material properties of interest (the ellipse on the bottom
+   right plot).
+
+In UQ process, we first quantify the uncertainty of the model parameters (represented by
+the dashed ellipse on the middle plot of the figure above). Having found the parametric
+uncertainty, then we can propagate the uncertainty of the parameters and get the
+uncertainty of the material properties of interest, e.g., by evaluating the ensemble that
+is obtained from sampling the distribution of the parameters. As the first uncertainty
+propagation is more involved, KLIFF implements tools to quantify the uncertainty of the
+parameters.
+
+In KLIFF, the UQ tools are implemented in :mod:`~kliff.uq`. In the current version, there
+are 2 methods implemented: Bayesian MCMC sampling and bootstrapping, with the integration
+of other UQ methods will be added in the future.
 
 
 MCMC
@@ -43,7 +59,7 @@ The likelihood function is given by
 
 The inclusion of the sampling temperature :math:`T` is to account for model inadequacy,
 or bias, in the potential [Kurniawan2022]_. Frederiksen et al. (2004) [Frederiksen2004]_
-suggest to estimate the bias by setting the temperature to
+suggest estimating the bias by setting the temperature to
 
 .. math::
 
@@ -71,16 +87,14 @@ parameters and :math:`N` is the number of tunable parameters.
 .. [KurniawanKLIFFUQ]
    Kurniawan, Y., Petrie, C.L., Transtrum, M.K., Tadmor, E.B., Elliott, R.S., Karls,
    D.S., Wen, M., 2022. Extending OpenKIM with an Uncertainty Quantification Toolkit for
-   Molecular Modeling. arXiv:2206.00578 [physics.comp-ph]
+   Molecular Modeling, in: 2022 IEEE 18th International Conference on E-Science
+   (e-Science). Presented at the 2022 IEEE 18th International Conference on e-Science
+   (e-Science), pp. 367–377. https://doi.org/10.1109/eScience55777.2022.00050
 
 
 
-MCMC implementation
-===================
-
-In KLIFF, the UQ tools are implemented in :mod:`~kliff.uq`. In the currect version, only
-MCMC sampling is implemented, with the integration of other UQ methods will be added in
-the future.
+Implementation
+--------------
 
 For the MCMC sampling, KLIFF adopts parallel-tempered MCMC (PTMCMC) methods, via the
 ptemcee_ Python package, as a way to perform MCMC sampling with several different
@@ -109,9 +123,9 @@ instantiate the sampler are:
 * ``nwalkers`` specifies the number of parallel walkers to run for each sampling
   temperature. As a default, this quantity is set to twice the number of parameters in
   the model.
-* ``logprior_fn`` argument allows user to specify the prior distribution
+* ``logprior_fn`` argument allows the user to specify the prior distribution
   :math:`\pi(\theta)` to use. The function should accept an array of parameter values as
-  an input and compute the logarithm of the prior distribution. Note that the prior
+  input and compute the logarithm of the prior distribution. Note that the prior
   distribution doesn't need to be normalized. The default prior is a uniform distribution
   over a finite range. See the next argument on how to set the boundaries of the uniform
   prior.
@@ -145,7 +159,7 @@ The required arguments are:
 * ``p0``, which is a :math:`K \times L \times N` array containing the position of each
   walker for each temperature in parameter space, where :math:`K`, :math:`L`, and
   :math:`N` are the number of temperatures, walkers, and parameters, respectively.
-* ``iterations`` specifies the number of MCMC step to take. Since the position is step
+* ``iterations`` specifies the number of MCMC steps to take. Since the position of step
   :math:`i` in Markov chain only depends on step :math:`(i-1)`, it is possible to break
   up the MCMC run into smaller batches, with the note that the initial positions of the
   current run needs to be set to the last positions of the previous run.
@@ -153,13 +167,13 @@ The required arguments are:
 .. seealso::
    For other possible arguments, see also ptemcee.Sampler.run_mcmc_.
 
-The resulting chain can be retrieved from via ``sampler.chain`` as a
+The resulting chain can be retrieved via ``sampler.chain`` as a
 :math:`K \times L \times M \times N` array, where :math:`M` is the total number of
-iteratiions.
+iterations.
 
 
 Parallelization
----------------
+^^^^^^^^^^^^^^^
 
 In principle, parallelization for the MCMC run can be done in 2 places: in the likelihood
 (or loss function) evaluation for each parameter set (see :ref:`run_in_parallel`) and in
@@ -167,14 +181,14 @@ the likelihood evaluation across different walkers. In the current implementatio
 supports OpenMP-style parallelization in the loss evaluation and both OpenMP and MPI for
 the sampling for different walkers when running MCMC sampling.
 
-In general, parallelization in the sampling process is done by declaring a pool and set
-it to ``sampler.pool`` prior to running MCMC, for example:
+In general, parallelization in the sampling process is done by declaring a pool and
+setting it to ``sampler.pool`` prior to running MCMC, for example:
 
 .. code-block:: python
 
    from multiprocessing import Pool
    
-   sampler.pool = Pool(nprocs)  # nprocs is the number of parallel process to use
+   sampler.pool = Pool(nprocs)  # nprocs is the number of parallel processes to use
    sampler.run_mcmc(p0, iterations, *args, **kwargs)
 
 To do parallelization with MPI, we can utilize ``MPIPool`` from schwimmbad_:
@@ -189,7 +203,7 @@ To do parallelization with MPI, we can utilize ``MPIPool`` from schwimmbad_:
 and run the Python script with ``mpiexec`` bash command.
 
 If enough compute resources are available, we can also employ a hybrid parallelization,
-for example, using ``multiprocessing`` in the loss evaluation (by specifying argument
+for example, using ``multiprocessing`` in the loss evaluation (by specifying the argument
 ``nprocs`` > 1) and MPI in the likelihood evaluation across different walkers. Then, we
 can run the Python script as follows.
 
@@ -206,9 +220,9 @@ can run the Python script as follows.
 
 
 MCMC analysis
-=============
+-------------
 
-The chains from the MCMC simulation needs to be processed. In a nutshell, the steps to
+The chains from the MCMC simulation need to be processed. In a nutshell, the steps to
 take are
 
 * Estimate the burn-in time and discard it from the beginning of the chain,
@@ -218,11 +232,11 @@ take are
 
 
 Burn-in time
-------------
+^^^^^^^^^^^^
 
-First we need to discard the first few iterations in the beginning of each chain as a
+First, we need to discard the first few iterations at the beginning of each chain as a
 burn-in time. This is similar to the equilibration time in a molecular dynamics
-simulation before the measurement. This action also ensure that the result is independent
+simulation before the measurement. This action also ensures that the result is independent
 of the initial positions of the walkers.
 
 KLIFF provides a function to estimate the burn-in time, based on the Marginal Standard
@@ -231,13 +245,13 @@ Error Rule (MSER). This can calculation can be done using the function
 performed for each temperature, walker, and parameter dimension separately.
 
 Autocorrelation length
-----------------------
+^^^^^^^^^^^^^^^^^^^^^^
 
-In Markov chain, the position at step :math:`i` is not independent from the previous step.
-However, after several iterations (denote this number by :math:`\tau`, which is the
+In the Markov chain, the position at step :math:`i` is not independent of the previous
+step. However, after several iterations (denote this number by :math:`\tau`, which is the
 autocorrelation length), the walker will "forget" where it started, i.e., the position at
 step :math:`i` is independent from step :math:`(i + \tau)`. Thus, we need to only keep
-every :math:`\tau \text{-th}` step to obtain the independent, unceorrelated samples.
+every :math:`\tau \text{-th}` step to obtain the independent, uncorrelated samples.
 
 The estimation of the autocorrelation length in KLIFF is done via the function
 :func:`~kliff.uq.mcmc_utils.autocorr`, which wraps over
@@ -248,7 +262,7 @@ respectively, and :math:`\tilde{M}` is the remaining number of iterations after
 discarding the burn-in time.
 
 Convergence
------------
+^^^^^^^^^^^
 
 Finally, after a sufficient number of iterations, the distribution of the MCMC samples
 will converge to the posterior. For multi-chain MCMC simulation, the convergence can be
@@ -261,14 +275,139 @@ In KLIFF, the function :func:`~kliff.uq.rhat` computes :math:`\hat{R}^p` for one
 temperature. The required input argument is a :math:`L \times \tilde{M}^* \times N` array
 of independent samples (:math:`\tilde{M}^*` is the number of independent samples in each
 walker). When the resulting :math:`\hat{R}^p` values are larger than the threshold
-(e.g., 1.1), then the MCMC simulation should be continued until this criteria is
+(e.g., 1.1), then the MCMC simulation should be continued until this criterion is
 satisfied.
 
 .. note::
-   Some sampling temperatures might converge at slower rates compared to the others.
+   Some sampling temperatures might converge at slower rates compared to others.
    So, user can terminate the MCMC simulation as long as the samples at the target
    temperatures, e.g., :math:`T_0`, have converged.
 
 
 .. seealso::
    See the tutorial for running MCMC in :ref:`tut_mcmc`.
+
+
+
+Bootstrap
+=========
+
+In general, the training dataset contains some random noise. When the data collection
+process is repeated, we will not get the exact same values, but instead, we will get
+(slightly) different values, where the deviation comes from the random noise. If we
+train the model to fit different realizations of the training dataset, we will get a
+distribution of the parameters. The uncertainty of the parameters from this distribution
+gives how the error in the training data is propagated to the uncertainty of the
+parameters. However, oftentimes we don't have the luxury to repeat the data collection.
+A suggestion, in this case, is to generate artificial datasets and train the model to fit
+these artificial datasets.
+
+.. figure:: ../img/bootstrap.png
+   :alt: Illustration on how bootstrapping works.
+
+Bootstrapping is a way to generate artificial datasets. We assume that the original
+dataset contains :math:`N` *independent and identically distributed (iid)* data points.
+An artificial, bootstrap dataset is generated by sample :math:`M` points from the original
+dataset with replacement. Note that this means that there are some data points that are
+repeated, while some other data points are not sampled, thus the bootstrap dataset is not
+the same as the original dataset. The difference between the datasets gives a sense of
+probability in data.
+
+.. note::
+   Although usually :math:`M` is set to be the same as :math:`N`, in principle it
+   doesn't need to be.
+
+
+
+Implementation
+--------------
+
+Bootstrapping is implemented in :class:`~kliff.uq.Bootstrap`. A general workflow for this
+calculation is
+
+1. Instantiate :class:`~kliff.uq.Bootstrap` class instance.
+
+   This process is straightforward. The only required argument is the :class:`~kliff.loss.Loss`
+   instance.
+
+   .. code-block:: python
+
+      from kliff.uq import Bootstrap
+
+      loss = ...  # define the loss function
+      # Train the potential
+      min_kwargs = ...  # Optimizer setting
+      loss.minimize(**min_kwargs)
+
+      bs = Bootstrap(loss, *args, **kwargs)
+
+   When instantiating the parent class :class:`~kliff.uq.Bootstrap`, it will return either
+   an instance of :class:`~kliff.uq.BootstrapEmpiricalModel` or
+   :class:`~kliff.uq.BootstrapNeuralNetworkModel`, depending on whether we have a
+   physics-based (empirical) model or a neural network model, respectively. When a neural
+   network model is used, the user can specify an additional argument `orig_state_filename`,
+   which specified the name and path of the file to use to export the initial state of the
+   model prior to running bootstrap. This is to reset the state of the model at the end
+   of performing bootstrap UQ.
+
+2. Generate bootstrap datasets.
+
+   In this implementation, we assume that the training dataset consists of many atomic
+   configurations and the corresponding quantities. Note that the quantities corresponding
+   to a single atomic configuration are **not** independent of each other. Thus, the
+   resampling process to generate a bootstrap dataset should not be done at the data point
+   level. Instead, we should generate a bootstrap dataset by resampling the atomic
+   configurations.
+
+   The built-in bootstrap dataset generator function was set up to perform this type of
+   resampling. Note that atomic configurations here are referred to as compute arguments,
+   which also contain the type of data and weights to use.
+
+   .. code-block:: python
+
+      nsamples = ...  # Number of samples to generate
+      bs.generate_bootstrap_compute_arguments(nsamples)
+
+   When an empirical model with multiple calculators is used, the resampling is done to
+   the combined list of the compute arguments across all calculators. Then, an internal
+   function will automatically assign back the bootstrap compute arguments to their
+   respective calculators. This means that the number of compute arguments in each
+   calculator when we do bootstrapping is more likely to be different than the original
+   number of compute arguments per calculator, although the total number of compute
+   arguments is still the same.
+
+   Also, note that the built-in bootstrap compute arguments generator assumes that the
+   configurations are independent of each other. In the case where this is not satisfied,
+   then a more sophisticated resampling method should be used. This can be done by
+   defining a custom bootstrap compute arguments generator function. The only required
+   argument for this function is the requested number of samples.
+
+3. Run the optimization for each bootstrap dataset.
+
+   After a set of bootstrap compute arguments is generated, then we need to iterate over
+   each of them, and train the potential to fit each bootstrap dataset.
+
+   .. code-block:: python
+
+      bs.run(min_kwargs=min_kwargs)
+
+   There are 2 arguments that are the same to run the optimization stage of bootstrapping,
+   regardless if we use an empirical or neural network model. These arguments are:
+
+   * ``min_kwargs``, which is a dictionary containing the keyword arguments that will be
+     passed into the optimizer. This argument can be thought of as the optimizer setting.
+
+     .. note::
+	Since the mapping from the bootstrap dataset to the inferred parameters contains
+	optimization, then it is recommended to use the same optimizer setting when we
+	iterate over each bootstrap compute argument and train the potential.
+	Additionally, the optimizer setting should also be the same as the setting used
+	in the initial training, when we use the original set of compute arguments to
+	train the potential.
+
+   * ``callback``, which is an option to specify a function that will be called in each
+     iteration. This can be used as a debugging tool or to monitor convergence.
+
+   For other additional arguments, please refer to the respective function documentation,
+   i.e., :meth:`~kliff.uq.BootstrapEmpiricalModel.run` for an empirical model or
+   :meth:`~kliff.uq.BootstrapNeuralNetworkModel.run` for neural network model.
diff --git a/examples/example_uq_bootstrap.py b/examples/example_uq_bootstrap.py
new file mode 100644
index 00000000..f253fc84
--- /dev/null
+++ b/examples/example_uq_bootstrap.py
@@ -0,0 +1,123 @@
+"""
+.. _tut_bootstrap:
+
+Bootstrapping
+=============
+
+In this example, we demonstrate how to perform uncertainty quantification (UQ) using
+bootstrap method. We use a Stillinger-Weber (SW) potential for silicon that is archived
+in OpenKIM_.
+
+For simplicity, we only set the energy-scaling parameters, i.e., ``A`` and ``lambda`` as
+the tunable parameters. These parameters will be calibrated to energies and forces of a
+small dataset, consisting of 4 compressed and stretched configurations of diamond silicon
+structure.
+"""
+
+
+##########################################################################################
+# To start, let's first install the SW model::
+#
+#    $ kim-api-collections-management install user SW_StillingerWeber_1985_Si__MO_405512056662_006
+#
+# .. seealso::
+#    This installs the model and its driver into the ``User Collection``. See
+#    :ref:`install_model` for more information about installing KIM models.
+
+
+import matplotlib.pyplot as plt
+import numpy as np
+
+from kliff.calculators import Calculator
+from kliff.dataset import Dataset
+from kliff.loss import Loss
+from kliff.models import KIMModel
+from kliff.uq.bootstrap import BootstrapEmpiricalModel
+from kliff.utils import download_dataset
+
+##########################################################################################
+# Before running bootstrap, we need to define a loss function and train the model. More
+# detail information about this step can be found in :ref:`tut_kim_sw` and
+# :ref:`tut_params_transform`.
+
+# Create the model
+model = KIMModel(model_name="SW_StillingerWeber_1985_Si__MO_405512056662_006")
+
+# Set the tunable parameters and the initial guess
+opt_params = {"A": [["default"]], "lambda": [["default"]]}
+
+model.set_opt_params(**opt_params)
+model.echo_opt_params()
+
+# Get the dataset
+dataset_path = download_dataset(dataset_name="Si_training_set_4_configs")
+# Read the dataset
+tset = Dataset(dataset_path)
+configs = tset.get_configs()
+
+# Create calculator
+calc = Calculator(model)
+# Only use the forces data
+ca = calc.create(configs, use_energy=False, use_forces=True)
+
+# Instantiate the loss function
+residual_data = {"normalize_by_natoms": False}
+loss = Loss(calc, residual_data=residual_data)
+
+##########################################################################################
+# To perform UQ by bootstrapping, the general workflow starts by instantiating
+# :class:`~kliff.uq.bootstrap.BootstrapEmpiricalModel`, or
+# :class:`~kliff.uq.bootstrap.BootstrapNeuralNetworkModel` if using a neural network
+# potential.
+
+
+# Instantiate bootstrap class object
+BS = BootstrapEmpiricalModel(loss, seed=1717)
+
+##########################################################################################
+# Then, we generate some bootstrap compute arguments. This is equivalent to generating
+# bootstrap data. Typically, we just need to specify how many bootstrap data samples to
+# generate. Additionally, if we call ``generate_bootstrap_compute_arguments`` multiple
+# times, the new generated data samples will be appended to the previously generated data
+# samples. This is also the behavior if we read the data samples from the previously
+# exported file.
+
+
+# Generate bootstrap compute arguments
+BS.generate_bootstrap_compute_arguments(100)
+
+##########################################################################################
+# Finally, we will iterate over these bootstrap data samples and train the potential
+# using each data sample. The resulting optimal parameters from each data sample give a
+# single sample of parameters. By iterating over all data samples, then we will get an
+# ensemble of parameters.
+#
+# Note that the mapping from the bootstrap dataset to the parameters involve optimization.
+# We suggest to use the same mapping, i.e., the same optimizer setting, in each iteration.
+# This includes using the same set of initial parameter guess. In the case when the loss
+# function has multiple local minima, we don't want the parameter ensemble to be biased
+# on the results of the other optimizations. For neural network model, we need to reset
+# the initial parameter value, which is done internally.
+
+
+# Run bootstrap
+min_kwargs = dict(method="lm")  # Optimizer setting
+initial_guess = calc.get_opt_params()  # Initial guess in the optimization
+BS.run(min_kwargs=min_kwargs, initial_guess=initial_guess)
+
+##########################################################################################
+# The resulting parameter ensemble can be accessed in `BS.samples` as a `np.ndarray`.
+# Then, we can plot the distribution of the parameters, as an example, or propagate the
+# error to the target quantities we want to study.
+
+
+# Plot the distribution of the parameters
+plt.figure()
+plt.plot(*(BS.samples.T), ".", alpha=0.5)
+param_names = list(opt_params.keys())
+plt.xlabel(param_names[0])
+plt.ylabel(param_names[1])
+plt.show()
+
+##########################################################################################
+# .. _OpenKIM: https://openkim.org
diff --git a/kliff/calculators/calculator.py b/kliff/calculators/calculator.py
index c0d25923..e8e11958 100644
--- a/kliff/calculators/calculator.py
+++ b/kliff/calculators/calculator.py
@@ -295,13 +295,7 @@ class _WrapperCalculator(object):
 
     # TODO methods like get_prediction needs to be implemented. The prediction of
     #  energy (force) should be the sum of each calculator.
-    def __init__(self, calculators=List[Calculator]):
-        """
-        Parameters
-        ----------
-
-        calculators: instance of Calculator
-        """
+    def __init__(self, calculators: List[Calculator]):
         self._calculators = calculators
         self._start_end = self._set_start_end()
 
@@ -324,11 +318,26 @@ def _set_start_end(self):
             i += n
         return start_end
 
-    def get_compute_arguments(self):
+    def get_compute_arguments(self, flat: bool = True) -> np.array:
+        """
+        Retrieve list of compute arguments. If ``flat=True``, then a 1D list will be
+        returned. Otherwise, nested lists (2D) will be returned, where each nested list
+        contain the compute arguments for each calculator.
+
+        Args:
+            flat: Whether return a flat, 1D array or nested lists. If nested lists is
+                returned, each nested list contains compute arguments for each calculator.
+
+        Returns:
+            all_cas: Compute arguments of all calculators.
+        """
         all_cas = []
         for calc in self._calculators:
             cas = calc.get_compute_arguments()
-            all_cas.extend(cas)
+            if flat:
+                all_cas.extend(cas)
+            else:
+                all_cas.append(cas)
         return all_cas
 
     def get_num_opt_params(self):
@@ -364,6 +373,20 @@ def get_calculator_list(self):
             calc_list.extend([calc] * N)
         return calc_list
 
+    def has_opt_params_bounds(self) -> bool:
+        """
+        Return a boolean to indicate whether there are parameters whose bounds are
+        provided.
+
+        Returns:
+            True if all calculators have bounds on the parameters, otherwise it is false.
+        """
+        calc_list = self._calculators
+        has_opt_params_bounds_per_calc = [
+            calc.model.has_opt_params_bounds() for calc in calc_list
+        ]
+        return all(has_opt_params_bounds_per_calc)
+
 
 class CalculatorError(Exception):
     def __init__(self, msg):
diff --git a/kliff/calculators/calculator_torch.py b/kliff/calculators/calculator_torch.py
index 9f255c28..c535b871 100644
--- a/kliff/calculators/calculator_torch.py
+++ b/kliff/calculators/calculator_torch.py
@@ -1,6 +1,7 @@
 from pathlib import Path
-from typing import Dict, List, Optional, Union
+from typing import Any, Dict, List, Optional, Tuple, Union
 
+import numpy as np
 import torch
 from loguru import logger
 from torch.utils.data import DataLoader
@@ -121,18 +122,38 @@ def create(
                 nprocs,
             )
 
+        # Finally, assign fingerprints dataset property as a FingerprintsDataset instance
+        self.fingerprints_dataset = FingerprintsDataset(self.fingerprints_path)
+
+    def get_fingerprints(self) -> List[dict]:
+        """
+        Return a list of fingerprints of the configurations.
+        """
+        return self.fingerprints_dataset.fp
+
     def get_compute_arguments(self, batch_size: int = 1):
         """
         Return the dataloader with batch size set to `batch_size`.
         """
-        fname = self.fingerprints_path
-        fp = FingerprintsDataset(fname)
         loader = DataLoader(
-            dataset=fp, batch_size=batch_size, collate_fn=fingerprints_collate_fn
+            dataset=self.fingerprints_dataset,
+            batch_size=batch_size,
+            collate_fn=fingerprints_collate_fn,
         )
 
         return loader
 
+    def set_fingerprints(self, fingerprints: List[dict]):
+        """
+        Update the fingerprints of the calculator. The fingerprints input argument should
+        be in the same format as the output of `meth:~kliff.descriptors.descriptor.load_fingerprints`,
+        which is a list of dictionaries.
+
+        Args:
+            fingerprints: A list of fingerprints.
+        """
+        self.fingerprints_dataset.fp = fingerprints
+
     def fit(self):
         path = self.fingerprints_path
         self.model.fit(path)
@@ -253,6 +274,85 @@ def _compute_stress(denergy_dzeta, dzetadr, volume):
         forces = torch.tensordot(denergy_dzeta, dzetadr, dims=([0, 1], [0, 1])) / volume
         return forces
 
+    def get_size_opt_params(self) -> Tuple[List[int], List[int], int]:
+        """
+        Return the size of the parameters.
+
+        Returns:
+            sizes: Each element in the list gives the shape of each type of parameter
+                tensors, containing, e.g., weights and biases, for each layer.
+            nelements: Number of elements of each parameter tensor.
+            nparams: Total number of parameters
+        """
+        sizes = []  # Size of each parameter tensor
+        nelements = []  # The number of elements for each tensor
+        for param in self.model.parameters():
+            sizes.append(param.size())
+            nelements.append(np.prod(param.size()))
+        nparams = sum(nelements)
+        return sizes, nelements, nparams
+
+    def get_num_opt_params(self) -> int:
+        """
+        Return the total number of parameters.
+        """
+        return self.get_size_opt_params()[-1]
+
+    def get_opt_params(self, flat: bool = True) -> Union[List, np.array]:
+        """
+        Retrieve the parameters, i.e., weights and biases.
+
+        Args:
+            flat: A flag to return a flat, 1D array.
+
+        Returns:
+            Parameters, i.e., weights and biases. If ``flat=True``, a 1D np.ndarray will
+            be returned. Otherwise, nested lists will be returned, where each list contain
+            the weights and biases for each layer.
+        """
+        parameters = []
+        for param in self.model.parameters():
+            if flat:
+                # Make sure that the parameters are stored in host memory
+                param_host = param.data.cpu()
+                parameters = np.append(parameters, param_host.numpy().flatten())
+            else:
+                parameters.append(param)
+        return parameters
+
+    def update_model_params(self, parameters: np.array):
+        """
+        Update the model parameters from a 1D array.
+
+        Args:
+            parameters: New parameter values to set. It needs to be a 1D array.
+        """
+        # Convert to the right format
+        parameters = self._convert_parameters_from_1d_array(parameters)
+        # Update the weights and biases
+        for ii, param in enumerate(self.model.parameters()):
+            param.data = parameters[ii]
+
+    def _convert_parameters_from_1d_array(self, flat_params: np.array) -> List:
+        """
+        Convert the parameters from a 1D array format to nested lists format.
+
+        Args:
+            flat_params: A 1D array containing weights and biases of the model.
+
+        Returns:
+            parameters: Parameters (weiths and biases) in nested lists format.
+        """
+        sizes, nelems, _ = self.get_size_opt_params()
+        # Indices to index the flat array to get the appropriate portion of each parameter
+        # tensor
+        idx = np.append(0.0, np.cumsum(nelems)).astype(int)
+        parameters = []
+        for ii, size in enumerate(sizes):
+            params = flat_params[idx[ii] : idx[ii + 1]]
+            parameters.append(torch.Tensor(params.reshape(size)))
+        return parameters
+
 
 class CalculatorTorchSeparateSpecies(CalculatorTorch):
     """
diff --git a/kliff/dataset/weight.py b/kliff/dataset/weight.py
index 66b7baf0..dfd2727d 100644
--- a/kliff/dataset/weight.py
+++ b/kliff/dataset/weight.py
@@ -144,8 +144,8 @@ def compute_weight(self, config):
             # Use the absolute value of the energy
             energy_norm = np.abs(energy)
             self._energy_weight = self._compute_weight_one_property(
-                energy_norm, self._weight_params["energy_weight_params"], "energy"
-            )
+                [energy_norm], self._weight_params["energy_weight_params"], "energy"
+            )[0]
         # Forces
         if forces is not None:
             # Use the magnitude of the force vector
@@ -163,8 +163,8 @@ def compute_weight(self, config):
             shear_stress_norm = np.linalg.norm(stress[3:])
             stress_norm = np.sqrt(normal_stress_norm**2 + 2 * shear_stress_norm**2)
             self._stress_weight = self._compute_weight_one_property(
-                stress_norm, self._weight_params["stress_weight_params"], "stress"
-            )
+                [stress_norm], self._weight_params["stress_weight_params"], "stress"
+            )[0]
 
         self._check_compute_flag(config)
 
@@ -174,7 +174,7 @@ def _compute_weight_one_property(data_norm, property_weight_params, property_typ
         Compute the weight based for one property.
         """
         c1, c2 = property_weight_params
-        sigma = np.linalg.norm([c1, (c2 * data_norm)])
+        sigma = np.array([np.linalg.norm([c1, c2 * dn]) for dn in data_norm])
         weight = 1 / sigma
         if np.any(sigma < 1e-12):
             logger.warning(
diff --git a/kliff/uq/__init__.py b/kliff/uq/__init__.py
index d44e3fb7..9a5009f5 100644
--- a/kliff/uq/__init__.py
+++ b/kliff/uq/__init__.py
@@ -1,3 +1,4 @@
+from .bootstrap import Bootstrap, BootstrapEmpiricalModel, BootstrapNeuralNetworkModel
 from .mcmc import MCMC, get_T0
 from .mcmc_utils import autocorr, mser, rhat
 
@@ -7,4 +8,7 @@
     "mser",
     "autocorr",
     "rhat",
+    "Bootstrap",
+    "BootstrapEmpiricalModel",
+    "BootstrapNeuralNetworkModel",
 ]
diff --git a/kliff/uq/bootstrap.py b/kliff/uq/bootstrap.py
new file mode 100644
index 00000000..eb32c55c
--- /dev/null
+++ b/kliff/uq/bootstrap.py
@@ -0,0 +1,773 @@
+import copy
+import json
+import os
+from pathlib import Path
+from typing import Any, Callable, List, Optional, Union
+
+import numpy as np
+
+from kliff.calculators.calculator import Calculator, _WrapperCalculator
+from kliff.calculators.calculator_torch import CalculatorTorchSeparateSpecies
+from kliff.loss import (
+    Loss,
+    LossNeuralNetworkModel,
+    LossPhysicsMotivatedModel,
+    energy_forces_residual,
+    energy_residual,
+    forces_residual,
+)
+
+
+class _BaseBootstrap:
+    """
+    Base class for bootstrap sampler.
+
+    Args:
+        loss: Loss function class instance from :class:`~kliff.loss.Loss`.
+        seed: Random number generator seed.
+    """
+
+    def __init__(self, loss: Loss, seed: Optional[int] = 1717):
+        self.loss = loss
+        self.calculator = loss.calculator
+        # Cache the original parameter values
+        self.orig_params = copy.copy(self.calculator.get_opt_params())
+        # Initiate the bootstrap configurations property
+        self.bootstrap_compute_arguments = {}
+        self.samples = np.empty((0, self.calculator.get_num_opt_params()))
+
+        # Set the random state
+        self.random_state = np.random.RandomState(seed)
+
+    @property
+    def _nsamples_done(self) -> int:
+        """
+        Returns the number of bootstrap compute arguments have been used in to train the
+        model.
+
+        For each bootstrap compute arguments sample, we need to train the potential using
+        these compute arguments. This function returns how many bootstrap compute
+        arguments samples have been used in the training. This is to help so that we can
+        continue the calculation midway.
+        """
+        return len(self.samples)
+
+    @property
+    def _nsamples_prepared(self) -> int:
+        """
+        Get how many bootstrap compute arguments are prepared.
+
+        This also include those compute arguments set that have been evaluated.
+        """
+        return len(self.bootstrap_compute_arguments)
+
+    def _generate_bootstrap_compute_arguments(
+        self, nsamples: int, bootstrap_cas_generator_fn: Callable, **kwargs
+    ):
+        """
+        Generate bootstrap compute arguments samples.
+
+        If this function is called multiple, say, K times, then it will in total
+        generate: math: `K \times nsamples` bootstrap compute arguments samples. That is,
+        consecutive call of this function will append the generated compute arguments
+        samples.
+
+        Args:
+            nsamples: Number of bootstrap samples to generate.
+            bootstrap_cas_generator_fn: A function to generate bootstrap compute argument
+                samples. The default function combine the compute arguments across all
+                calculators and do sampling with replacement from the combined list.
+                Another possible convention is to do sampling with replacement on the
+                compute arguments list of each calculator separately, in which case a
+                custom function needs to be defined and used.
+            kwargs: Additional keyword arguments to ``bootstrap_cas_generator_fn``.
+        """
+        # Generate a new bootstrap configurations
+        new_bootstrap_compute_arguments = bootstrap_cas_generator_fn(nsamples, **kwargs)
+        # Update the old list with this new list
+        self.bootstrap_compute_arguments = self._update_bootstrap_compute_arguments(
+            new_bootstrap_compute_arguments
+        )
+
+    def _update_bootstrap_compute_arguments(
+        self, new_bootstrap_compute_arguments: dict
+    ) -> dict:
+        """
+        Append the new generated bootstrap compute arguments samples to the old list.
+
+        Args:
+            new_bootstrap_compute_arguments: Newly generated bootstrap compute arguments
+                to append to the old dictionary of bootstrap compute arguments.
+
+        Returns:
+            Updated dictionary containing all generated compute arguments, including the
+            newly generated ones.
+        """
+        bootstrap_compute_arguments = copy.copy(self.bootstrap_compute_arguments)
+        for ii, vals in new_bootstrap_compute_arguments.items():
+            iteration = ii + self._nsamples_prepared
+            bootstrap_compute_arguments.update({iteration: vals})
+
+        return bootstrap_compute_arguments
+
+    def _save_bootstrap_compute_arguments(
+        self, filename: Union[Path, str], identifiers_converter_fn: Callable
+    ):
+        """
+        Export the generated bootstrap compute arguments as a json file.
+
+        The json file will contain the identifier of the compute arguments for each
+        sample.
+
+        Args:
+            filename: Where to export the bootstrap compute arguments samples
+            identifiers_conferter_fn: A function to convert the compute arguments to their
+                identifiers.
+        """
+        # We cannot directly export the bootstrap compute arguments. Instead, we will
+        # first convert it and list the identifiers and export the identifiers.
+        # Convert to identifiers
+        bootstrap_compute_arguments_identifiers = {}
+        for ii in self.bootstrap_compute_arguments:
+            bootstrap_compute_arguments_identifiers.update(
+                {ii: identifiers_converter_fn(self.bootstrap_compute_arguments[ii])}
+            )
+
+        with open(filename, "w") as f:
+            json.dump(bootstrap_compute_arguments_identifiers, f, indent=4)
+
+    def reset(self):
+        """
+        Reset the bootstrap sampler.
+        """
+        self.restore_loss()
+        self.bootstrap_compute_arguments = {}
+        self.samples = np.empty((0, self.calculator.get_num_opt_params()))
+
+
+class Bootstrap:
+    """
+    Bootstrap sampler class for uncertainty quantification.
+
+    This is a wrapper over :class:`BootstrapEmpiricalModel` and
+    :class:`BootstrapNeuralNetworkModel` to provide a united interface. You can use the
+    two classes directly.
+
+    Args:
+        loss: Loss function class instance from :class:`~kliff.loss.Loss`.
+        seed: Random number generator seed.
+        args, kwargs: Additional positional and keyword arguments for instantiating
+            :class:`BootstrapEmpiricalModel` or :class:`BootstrapNeuralNetworkModel`.
+    """
+
+    def __new__(self, loss: Loss, seed: Optional[int] = 1717, *args, **kwargs):
+        if isinstance(loss, LossPhysicsMotivatedModel):
+            return BootstrapEmpiricalModel(loss, seed, *args, **kwargs)
+        elif isinstance(loss, LossNeuralNetworkModel):
+            return BootstrapNeuralNetworkModel(loss, seed, *args, **kwargs)
+
+
+def bootstrap_cas_generator_empirical(
+    nsamples: int,
+    orig_cas: List,
+    ncas: Optional[int] = None,
+    rstate: Optional[np.random.RandomState] = None,
+) -> dict:
+    """
+    Default class to generate bootstrap compute arguments for empirical, physics-based
+    model.
+
+    The compute arguments from all calculators will be combined, then the bootstrap
+    sample configurations will be generated from the combined list. Afterwards, the
+    configurations will be split into their respective calculators.
+
+    Args:
+        nsamples: Number of the bootstrap compute arguments requested.
+        orig_cas: The original list of compute arguments. The bootstrap compute arguments
+            will be generated from this list. The format of this input is given below::
+
+                orig_cas = [
+                    [calc0_ca0, calc0_ca1, ...],
+                    [calc1_ca0, calc1_ca1, ...],
+                    ...
+                ]
+
+        ncas: Number of compute arguments to have in each sample. If not specified, the
+            function will generate the same number of compute arguments sample as the
+            number of the original compute argument list.
+        rstate: The state of random number generator.
+
+
+    Returns:
+        A set of bootstrap compute arguments, written in a dictionary format, where the
+        keys index the bootstrap samples compute arguments::
+
+            bootstrap_configs = {
+                0: [[calc0_cas], [calc1_cas]],
+                1: [[calc0_cas], [calc1_cas]]
+            }
+
+    """
+    if rstate is None:
+        state = np.random.get_state()  # Get the state of global random number generator
+        rstate = np.random.RandomState()  # Instantiate a local random state
+        rstate.set_state(state)  # Set the state
+
+    ncalc = len(orig_cas)  # Number of calculators
+    # Number of compute args per calc
+    ncas_per_calc = [len(cas) for cas in orig_cas]
+    if ncas is None:
+        ncas = sum(ncas_per_calc)  # Total number of compute arguments
+    # Combine the compute arguments
+    comb_orig_cas = np.concatenate((orig_cas))
+    # Index of which calculator each ca correspond to
+    calc_idx = np.concatenate([[ii] * nc for ii, nc in enumerate(ncas_per_calc)])
+    bootstrap_cas = {}
+    for ii in range(nsamples):
+        # Generate a bootstrap sample configuration
+        # Generate the bootstrap indices
+        bootstrap_idx = rstate.choice(range(ncas), size=ncas, replace=True)
+        # From the indices, get bootstrap compute arguments
+        comb_bootstrap_cas = [comb_orig_cas[ii] for ii in bootstrap_idx]
+        # We also need to deal with the calculator index
+        comb_bootstrap_calc_idx = calc_idx[bootstrap_idx]
+
+        # Split the bootstrap cas into separate calculators
+        bootstrap_cas_single_sample = [[] for _ in range(ncalc)]
+        for idx, ca in zip(comb_bootstrap_calc_idx, comb_bootstrap_cas):
+            bootstrap_cas_single_sample[idx].append(ca)
+
+        # Update the bootstrap compute arguments dictionary
+        bootstrap_cas.update({ii: bootstrap_cas_single_sample})
+    return bootstrap_cas
+
+
+def get_identifiers_from_compute_arguments(compute_arguments: List) -> List[List[str]]:
+    """
+    Retrieve the identifiers from a list of compute arguments.
+
+    Args:
+        compute_arguments: A list of :class:`~kliff.models.model.ComputeArguments`.
+
+    Returns:
+        A list of compute arguments' identifiers, which shows the paths to the xyz files.
+    """
+    identifiers = []
+    for cas in compute_arguments:
+        # Iterate over compute arguments corresponding to each calculator
+        identifiers.append([ca.conf.identifier for ca in cas])
+    return identifiers
+
+
+def default_callback(*args):
+    """
+    Default callback function that does nothing.
+    """
+    return False
+
+
+class BootstrapEmpiricalModel(_BaseBootstrap):
+    """
+    Bootstrap sampler class for empirical, physics-based potentials.
+
+    Args:
+        loss: Loss function class instance from :class:`~kliff.loss.Loss`.
+        seed: Random number generator seed.
+    """
+
+    def __init__(self, loss: Loss, seed: Optional[int] = 1717):
+        super().__init__(loss, seed)
+        # Cache the original compute arguments
+        if isinstance(self.calculator, Calculator):
+            self.orig_compute_arguments = [
+                copy.copy(self.calculator.get_compute_arguments())
+            ]
+            self.use_multi_calc = False
+        elif isinstance(self.calculator, _WrapperCalculator):
+            self.orig_compute_arguments = copy.copy(
+                self.calculator.get_compute_arguments(flat=False)
+            )
+            self.use_multi_calc = True
+        self._orig_compute_arguments_identifiers = (
+            get_identifiers_from_compute_arguments(self.orig_compute_arguments)
+        )
+
+    def generate_bootstrap_compute_arguments(
+        self,
+        nsamples: int,
+        bootstrap_cas_generator_fn: Optional[Callable] = None,
+        **kwargs,
+    ):
+        """
+        Generate bootstrap compute arguments samples.
+
+        If this function is called multiple, say, K times, then it will in total
+        generate: math: `K \times nsamples` bootstrap compute arguments samples. That is,
+        consecutive call of this function will append the generated compute arguments
+        samples.
+
+        Args:
+            nsamples: Number of bootstrap samples to generate.
+            bootstrap_cas_generator_fn: A function to generate bootstrap compute argument
+                samples. The default function combine the compute arguments across all
+                calculators and do sampling with replacement from the combined list.
+                Another possible convention is to do sampling with replacement on the
+                compute arguments list of each calculator separately, in which case a
+                custom function needs to be defined and used. The required argument for
+                the custom generator functions is the requested number of samples.
+            kwargs: Additional keyword arguments to ``bootstrap_cas_generator_fn``.
+        """
+        # Function to generate bootstrap configurations
+        if bootstrap_cas_generator_fn is None:
+            bootstrap_cas_generator_fn = bootstrap_cas_generator_empirical
+            kwargs.update(
+                {"orig_cas": self.orig_compute_arguments, "rstate": self.random_state}
+            )
+        self._generate_bootstrap_compute_arguments(
+            nsamples, bootstrap_cas_generator_fn, **kwargs
+        )
+
+    def save_bootstrap_compute_arguments(self, filename: Union[Path, str]):
+        """
+        Export the generated bootstrap compute arguments as a json file.
+
+        The json file will contain the identifier of the compute arguments for each
+        sample.
+
+        Args:
+            filename: Where to export the bootstrap compute arguments samples
+        """
+        self._save_bootstrap_compute_arguments(
+            filename, get_identifiers_from_compute_arguments
+        )
+
+    def load_bootstrap_compute_arguments(self, filename: Union[Path, str]) -> dict:
+        """
+        Load the bootstrap compute arguments from a json file.
+
+        If a list of bootstrap compute arguments samples exists prior to this function
+        call, then the samples read from this file will be appended to the old list.
+
+        Args:
+            filename: Name or path of json file to read.
+
+        Returns:
+            Dictionary read from the json file.
+        """
+        # Load the json file
+        with open(filename, "r") as f:
+            new_bootstrap_compute_arguments_identifiers = json.load(f)
+
+        # The information stored in the json file are the identifiers. We need to
+        # convert it back to compute arguments.
+        keys = [int(key) for key in new_bootstrap_compute_arguments_identifiers.keys()]
+        new_bootstrap_compute_arguments = {}
+        # Iterate over sample
+        for ii in keys:
+            # List of identifier for step ii
+            identifiers_ii = new_bootstrap_compute_arguments_identifiers[str(ii)]
+            # Iterate over the calculator
+            cas_ii = []
+            for jj, identifiers_calc in enumerate(identifiers_ii):
+                reference = self._orig_compute_arguments_identifiers[jj]
+                cas_calc = [
+                    self.orig_compute_arguments[jj][reference.index(ss)]
+                    for ss in identifiers_calc
+                ]
+                cas_ii.append(cas_calc)
+            # Update the new bootstrap compute arguments dictionary
+            new_bootstrap_compute_arguments.update({ii: cas_ii})
+
+        # Update the old list with this new list
+        self.bootstrap_compute_arguments = self._update_bootstrap_compute_arguments(
+            new_bootstrap_compute_arguments
+        )
+        return new_bootstrap_compute_arguments_identifiers
+
+    def run(
+        self,
+        min_kwargs: Optional[dict] = None,
+        initial_guess: Optional[np.ndarray] = None,
+        residual_fn_list: Optional[List] = None,
+        callback: Optional[Callable] = None,
+    ) -> np.ndarray:
+        """
+        Iterate over the generated bootstrap compute arguments samples and train the
+        potential using each compute arguments sample.
+
+        Args:
+            min_kwargs: Keyword arguments for :meth:`~kliff.loss.Loss.minimize`.
+            initial_guess: (ndim,) Initial guess of parameters to use for the
+                minimization. It is recommended to use the same values as used in the
+                training process if such step is done prior to running bootstrap.
+            residual_fn_list: List of residual function to use in each calculator.
+                Currently, this only affect the case when multiple calculators are used.
+                If there is only a single calculator, don't worry about this argument.
+            callback: Called after each iteration. The arguments for this function are
+                the bootstrap instance and and output of
+                :meth:`~kliff.loss.Loss.minimize`. This function can also be used to
+                break the run, by returning boolean `True`.
+
+        Returns:
+            (nsamples, ndim,) Parameter samples from bootstrapping.
+
+        Raises:
+            BootstrapError: If there is no bootstrap compute areguments generated prior to
+                calling this method.
+            ValueError: If the calculators use neither the energy nor forces.
+        """
+        if self._nsamples_prepared == 0:
+            # Bootstrap compute arguments have not been generated
+            raise BootstrapError("Please generate a bootstrap compute arguments first")
+
+        # Optimizer setting
+        if min_kwargs is None:
+            min_kwargs = {}
+
+        # Callback function
+        if callback is None:
+            callback = default_callback
+
+        # Train the model using each bootstrap compute arguments
+        for ii in range(self._nsamples_done, self._nsamples_prepared):
+            # Update the compute arguments
+            if self.use_multi_calc:
+                # There are multiple calculators used
+                for jj, calc in enumerate(self.calculator.calculators):
+                    calc.compute_arguments = self.bootstrap_compute_arguments[ii][jj]
+            else:
+                self.calculator.compute_arguments = self.bootstrap_compute_arguments[
+                    ii
+                ][0]
+
+            # Set the initial parameter guess
+            if initial_guess is None:
+                initial_guess = self.calculator.get_opt_params().copy()
+            self.calculator.update_model_params(initial_guess)
+            # TODO This assumes that we use the built-in residual functions
+            if self.use_multi_calc:
+                if residual_fn_list is None:
+                    # If multiple calculators are used, we need to update the residual
+                    # function used for each configuration. This is to ensure that we use
+                    # the correct residual function for each configuration.
+                    calc_list = self.calculator.get_calculator_list()
+                    residual_fn_list = []
+                    for calculator in calc_list:
+                        if calculator.use_energy and calculator.use_forces:
+                            residual_fn = energy_forces_residual
+                        elif calculator.use_energy:
+                            residual_fn = energy_residual
+                        elif calculator.use_forces:
+                            residual_fn = forces_residual
+                        else:
+                            raise ValueError(
+                                "Calculator does not use energy or forces."
+                            )
+                        residual_fn_list.append(residual_fn)
+                self.loss.residual_fn = residual_fn_list
+            # Minimization
+            opt_res = self.loss.minimize(**min_kwargs)
+
+            # Append the parameters to the samples
+            self.samples = np.row_stack(
+                (self.samples, self.loss.calculator.get_opt_params())
+            )
+
+            # Callback
+            if callback(self, opt_res):
+                break
+
+        # Finishing up
+        self.restore_loss()  # Restore the loss function
+        return self.samples
+
+    def restore_loss(self):
+        """
+        Restore the loss function: revert back the compute arguments and the parameters
+        to the original state.
+        """
+        # Restore the parameters and configurations back
+        self.calculator.compute_arguments = self.orig_compute_arguments
+        self.calculator.update_model_params(self.orig_params)
+
+
+def bootstrap_cas_generator_neuralnetwork(
+    nsamples: int,
+    orig_fingerprints: List,
+    nfingerprints: Optional[int] = None,
+    rstate: Optional[np.random.RandomState] = None,
+) -> dict:
+    """
+    Default class to generate bootstrap compute arguments (fingerprints) for neural
+    network model.
+
+    When models for separate species are used, we assume that the compute arguments used
+    for those models are the same. This assumption is valid since usually the atomic
+    configurations contain multiple atom species, which are the same as the species used
+    for those models.
+
+    Args:
+        nsamples: Number of the bootstrap compute arguments requested.
+        orig_fingerprints: The original list of compute arguments(fingerprints). The
+            bootstrap compute arguments will be generated from this list. The format of
+            this input is given below: :
+
+                orig_fingerprints = [ca0, ca1, ...]
+
+        nfingerprints: Number of compute arguments to have in each sample. If not specified, the
+            function will generate the same number of compute arguments sample as the
+            number of the original compute argument list.
+        rstate: The state of random number generator
+
+    Returns:
+       A set of bootstrap compute arguments(fingerprints), written in a dictionary
+       format, where the keys index the bootstrap samples compute arguments::
+
+            bootstrap_configs = {
+                0: [ca0, ca1, ...],
+                1: [ca0, ca1, ...],
+            }
+
+    """
+    if rstate is None:
+        state = np.random.get_state()  # Get the state of global random number generator
+        rstate = np.random.RandomState()  # Instantiate a local random state
+        rstate.set_state(state)  # Set the state
+
+    bootstrap_fingerprints = {}
+    if nfingerprints is None:
+        nfingerprints = len(orig_fingerprints)
+    for ii in range(nsamples):
+        # Get 1 sample of bootstrap fingerprints
+        bootstrap_fingerprints_single_sample = rstate.choice(
+            orig_fingerprints, size=nfingerprints, replace=True
+        )
+        bootstrap_fingerprints.update({ii: bootstrap_fingerprints_single_sample})
+    return bootstrap_fingerprints
+
+
+def get_identifiers_from_fingerprints(fingerprints: List) -> List[str]:
+    """
+    Retrieve the identifiers of a list of fingerprints.
+
+    Args:
+        fingerprints: A list of fingerprints.
+
+    Returns:
+        A list of fingerprints' identifiers, which shows the paths to the xyz files.
+    """
+    identifiers = [fp["configuration"].identifier for fp in fingerprints]
+    return identifiers
+
+
+class BootstrapNeuralNetworkModel(_BaseBootstrap):
+    """
+    Bootstrap sampler class for neural network potentials.
+
+    Args:
+        loss: Loss function class instance from :class:`~kliff.loss.Loss`.
+        seed: Random number generator seed.
+        orig_state_filename: Name of the file in which the initial state of the model
+            prior to bootstrapping will be stored. This is to use at the end of the
+            bootstrap run to reset the model to the initial state.
+    """
+
+    def __init__(
+        self,
+        loss: Loss,
+        seed: Optional[int] = 1717,
+        orig_state_filename: Optional[Union[Path, str]] = "orig_model.pkl",
+    ):
+        super().__init__(loss, seed)
+        # Check if the calculator uses separate species
+        if isinstance(self.calculator, CalculatorTorchSeparateSpecies):
+            self._calc_separate_species = True
+            self.model = [model[1] for model in self.calculator.models.items()]
+            self._species = self.calculator.models
+        else:
+            self._calc_separate_species = False
+            self.model = [self.calculator.model]
+        # Cache the original fingerprints
+        self.orig_compute_arguments = self.calculator.get_fingerprints()
+        self._orig_compute_arguments_identifiers = get_identifiers_from_fingerprints(
+            self.orig_compute_arguments
+        )
+
+        # Save the original state of the model before running bootstrap
+        if self._calc_separate_species:
+            self.orig_state_filename = []
+            for sp, model in zip(self._species, self.model):
+                splitted_path = os.path.splitext(orig_state_filename)
+                path_with_species = splitted_path[0] + f"_{sp}" + splitted_path[1]
+                self.orig_state_filename.append(path_with_species)
+                model.save(path_with_species)
+        else:
+            self.orig_state_filename = [orig_state_filename]
+            self.model[0].save(orig_state_filename)
+
+    def generate_bootstrap_compute_arguments(
+        self,
+        nsamples: int,
+        bootstrap_cas_generator_fn: Optional[Callable] = None,
+        **kwargs,
+    ):
+        """
+        Generate bootstrap compute arguments samples.
+
+        If this function is called multiple, say, K times, then it will in total
+        generate: math: `K \times nsamples` bootstrap compute arguments samples. That is,
+        consecutive call of this function will append the generated compute arguments
+        samples.
+
+        Args:
+            nsamples: Number of bootstrap samples to generate.
+            bootstrap_cas_generator_fn: A function to generate bootstrap compute argument
+                samples. The default function combine the compute arguments across all
+                calculators and do sampling with replacement from the combined list.
+                Another possible convention is to do sampling with replacement on the
+                compute arguments list of each calculator separately, in which case a
+                custom function needs to be defined and used.
+            kwargs: Additional keyword arguments to ``bootstrap_cas_generator_fn``.
+        """
+        # Function to generate bootstrap configurations
+        if bootstrap_cas_generator_fn is None:
+            bootstrap_cas_generator_fn = bootstrap_cas_generator_neuralnetwork
+            kwargs.update(
+                {
+                    "orig_fingerprints": self.orig_compute_arguments,
+                    "rstate": self.random_state,
+                }
+            )
+        self._generate_bootstrap_compute_arguments(
+            nsamples, bootstrap_cas_generator_fn, **kwargs
+        )
+
+    def save_bootstrap_compute_arguments(self, filename: Union[Path, str]):
+        """
+        Export the generated bootstrap compute arguments as a json file.
+
+        The json file will contain the identifier of the compute arguments for each
+        sample.
+
+        Args:
+            filename: Where to export the bootstrap compute arguments samples
+        """
+        self._save_bootstrap_compute_arguments(
+            filename, get_identifiers_from_fingerprints
+        )
+
+    def load_bootstrap_compute_arguments(self, filename: Union[Path, str]) -> dict:
+        """
+        Load the bootstrap compute arguments from a json file.
+
+        If a list of bootstrap compute arguments samples exists prior to this function
+        call, then the samples read from this file will be appended to the old list.
+
+        Args:
+            filename: Name or path of json file to read.
+
+        Returns:
+            Dictionary read from the json file.
+        """
+        # Load the json file
+        with open(filename, "r") as f:
+            new_bootstrap_compute_arguments_identifiers = json.load(f)
+
+        # The information stored in the json file are the identifiers. We need to
+        # convert it back to fingerprints.
+        keys = [int(key) for key in new_bootstrap_compute_arguments_identifiers.keys()]
+        new_bootstrap_compute_arguments = {}
+        # Iterate over sample
+        for ii in keys:
+            # List of identifier for step ii
+            identifiers_ii = new_bootstrap_compute_arguments_identifiers[str(ii)]
+            reference = self._orig_compute_arguments_identifiers
+            fp_ii = [
+                self.orig_compute_arguments[reference.index(ss)]
+                for ss in identifiers_ii
+            ]
+            # Update the new bootstrap fingerprints dictionary
+            new_bootstrap_compute_arguments.update({ii: fp_ii})
+
+        # Update the old list with this new list
+        self.bootstrap_compute_arguments = self._update_bootstrap_compute_arguments(
+            new_bootstrap_compute_arguments
+        )
+        return new_bootstrap_compute_arguments_identifiers
+
+    def run(
+        self, min_kwargs: Optional[dict] = None, callback: Optional[Callable] = None
+    ) -> np.ndarray:
+        """
+        Iterate over the generated bootstrap compute arguments samples and train the
+        potential using each compute arguments sample.
+
+        Args:
+            min_kwargs: Keyword arguments for :meth:`~kliff.loss.Loss.minimize`.
+            callback: Called after each iteration. The arguments for this function are
+                the bootstrap instance and and output of
+                :meth:`~kliff.loss.Loss.minimize`. This function can also be used to
+                break the run, by returning boolean `True`.
+
+        Returns:
+            (nsamples, ndim,) Parameter samples from bootstrapping.
+
+        Raises:
+            BootstrapError: If there is no bootstrap compute areguments generated prior to
+                calling this method.
+        """
+        if self._nsamples_prepared == 0:
+            # Bootstrap fingerprints have not been generated
+            raise BootstrapError("Please generate a bootstrap compute_arguments first")
+
+        # Optimizer setting
+        if min_kwargs is None:
+            min_kwargs = {}
+
+        # Callback function
+        if callback is None:
+            callback = default_callback
+
+        # Train the model using each bootstrap fingerprints
+        for ii in range(self._nsamples_done, self._nsamples_prepared):
+            # Update the fingerprints
+            self.calculator.set_fingerprints(self.bootstrap_compute_arguments[ii])
+
+            for model in self.model:
+                # Reset the initial parameters
+                for layer in model.layers:
+                    try:
+                        layer.reset_parameters()
+                    except AttributeError:
+                        pass
+            # Minimization
+            self.loss.minimize(**min_kwargs)
+
+            # Append the parameters to the samples
+            self.samples = np.row_stack(
+                (self.samples, self.loss.calculator.get_opt_params())
+            )
+
+            # Callback
+            if callback(self):
+                break
+
+        # Finishing up, restore the state
+        self.restore_loss()
+        return self.samples
+
+    def restore_loss(self):
+        """
+        Restore the loss function: revert back the compute arguments and the parameters
+        to the original state.
+        """
+        # Restore the parameters and configurations back by loading the original state
+        # back.
+        for model, fname in zip(self.model, self.orig_state_filename):
+            model.load(fname)
+
+
+class BootstrapError(Exception):
+    def __init__(self, msg: str):
+        super(BootstrapError, self).__init__(msg)
+        self.msg = msg
diff --git a/kliff/uq/mcmc.py b/kliff/uq/mcmc.py
index 75d8a98b..cad4fba3 100644
--- a/kliff/uq/mcmc.py
+++ b/kliff/uq/mcmc.py
@@ -22,20 +22,17 @@
 
 
 def logprior_uniform(x: np.ndarray, bounds: np.ndarray) -> float:
-    """Logarithm of the non-normalized joint uniform prior. This is the default prior
-    distribution.
-
-    Parameters
-    ----------
-    x : np.ndarray
-        Parameter values to evaluate.
-    bounds : np.ndarray (nparams, 2,)
-        An array containing the boundaries of the uniform prior. The first column of the
-        array contains the lower bounds and the second column contains the upper bounds.
-
-    Returns
-    -------
-    float:
+    """Logarithm of the non-normalized joint uniform prior.
+
+    This is the default prior distribution.
+
+    Args:
+        x: (ndim,) Parameter values to evaluate.
+        bounds: (ndim, 2,) An array containing the boundaries of the uniform prior. The
+            first column of the array contains the lower bounds and the second column
+            contains the upper bounds.
+
+    Returns:
         Logarithm of the non-normalized joint uniform prior evaluated at parameter ``x``.
     """
     l_bounds, u_bounds = bounds.T
@@ -47,18 +44,15 @@ def logprior_uniform(x: np.ndarray, bounds: np.ndarray) -> float:
     return ret
 
 
-def get_T0(loss: Loss):
-    """Compute the natural temperature. The minimum loss is the loss value at the optimal
-    parameters.
+def get_T0(loss: Loss) -> float:
+    """Compute the natural temperature.
+
+    The minimum loss is the loss value at the optimal parameters.
 
-    Parameters
-    ----------
-    loss : Loss
-        Loss function class from :class:`~kliff.loss.Loss`.
+    Args:
+        loss: Loss function class from :class:`~kliff.loss.Loss`.
 
-    Returns
-    -------
-    float
+    Returns:
         Value of the natural temperature.
     """
     xopt = loss.calculator.get_opt_params()
@@ -73,29 +67,21 @@ class MCMC:
     This is a wrapper over :class:`PtemceeSampler` and
     :class:`EmceeSampler`. Currently, only these  2 samplers implemented.
 
-    Parameters
-    ----------
-    loss : Loss
-        Loss function class from :class:`~kliff.loss.Loss`.
-    nwalkers : Optional[int]
-        Number of walkers to simulate. The minimum number of walkers
-        is twice the number of parameters. It defaults to this minimum
-        value.
-    logprior_fn : Optional[Callable]
-        A function that evaluate logarithm of the prior
-        distribution. The prior doesn't need to be normalized. It
-        defaults to a uniform prior over a finite range.
-    logprior_args : Optional[tuple]
-        Additional positional arguments of the ``logprior_fn``. If the
-        default ``logprior_fn`` is used, then the boundaries of the
-        uniform prior can be specified here.
-    sampler : Optional[str] or sampler instance
-        An argument that specifies the MCMC sampler to use. The value can be one of the
-        strings ``"ptemcee"`` (the default value) or ``"emcee"``, or a sampler class
-        instance. If ``"ptemcee"`` or ``"emcee"`` is given, a respective internal sampler
-        class will be uses.
-    **kwargs : Optional[dict]
-        Additional keyword arguments for ``ptemcee.Sampler`` or ``emcee.EnsembleSampler``.
+    Args:
+        loss: Loss function class from :class:`~kliff.loss.Loss`.
+        nwalkers: Number of walkers to simulate. The minimum number of walkers is twice
+        the number of parameters. It defaults to this minimum value.
+        logprior_fn: A function that evaluate logarithm of the prior distribution. The
+            prior doesn't need to be normalized. It defaults to a uniform prior over a
+            finite range.
+        logprior_args: Additional positional arguments of the ``logprior_fn``. If the
+            default ``logprior_fn`` is used, then the boundaries of the uniform prior can
+            be specified here.
+        sampler: An argument that specifies the MCMC sampler to use. The value can be one
+            of the strings ``"ptemcee"`` (the default value) or ``"emcee"``, or a sampler
+            class instance. If ``"ptemcee"`` or ``"emcee"`` is given, a respective
+            internal sampler class will be uses.
+        **kwargs: Additional keyword arguments for ``ptemcee.Sampler`` or ``emcee.EnsembleSampler``.
     """
 
     builtin_samplers = ["ptemcee", "emcee"]
@@ -166,45 +152,30 @@ def __new__(
     class PtemceeSampler(ptemcee.Sampler):
         """Sampler class for PTMCMC via ``ptemcee`` Python package.
 
-        Parameters
-        ----------
-        loss : Loss
-            Loss function instance from :class:`~kliff.loss.Loss`.
-        nwalkers : Optional[int]
-            Number of walkers to simulate. The minimum number of walkers
-            is twice the number of parameters. It defaults to this minimum
-            value.
-        ntemps: Optional[int]
-            Number of temperatures to simulate. It defaults to 10.
-        Tmax_ratio: Optional[float]
-            The ratio between the highest temperature to use and the natural temperature.
-            Higher value means that the maximum temperature is higher than :math:`T_0`
-            [Frederiksen2004]_. It defaults to 1.0.
-        Tladder: Optional[List]
-            A list containing the temperature values to use. The values nedd to be
-            monotonically increasing or decreasing. It defaults to ``None``, which will
-            be generated from ``ntemps`` and ``Tmax_ratio``.
-        logprior_fn : Optional[Callable]
-            A function that evaluate logarithm of the prior
-            distribution. The prior doesn't need to be normalized. It
-            defaults to a uniform prior over a finite range.
-        logprior_args : Optional[tuple]
-            Additional positional arguments of the ``logprior_fn``. If the
-            default ``logprior_fn`` is used, then the boundaries of the
-            uniform prior can be specified here.
-        **kwargs : Optional[dict]
-            Additional keyword arguments for ``ptemcee.Sampler``.
-
-        Attributes
-        ----------
-        loss: Loss
-            Loss function instance from :class:`~kliff.loss.Loss`
-        T0: float
-            Values of the natural temperature, :math:`T_0` [Frederiksen2004]_.
-        Tladder: np.ndarray
-            An array containing the values of the sampling temperatures.
-        sampler: ptemcee.Sampler
-            Sampler instance from ``ptemcee.Sampler``
+        Args:
+            loss: Loss function instance from :class:`~kliff.loss.Loss`.
+            nwalkers: Number of walkers to simulate. The minimum number of walkers is
+                twice the number of parameters. It defaults to this minimum value.
+            ntemps: Number of temperatures to simulate. It defaults to 10.
+            Tmax_ratio: The ratio between the highest temperature to use and the natural
+                temperature. Higher value means that the maximum temperature is higher
+                than :math:`T_0` [Frederiksen2004]_. It defaults to 1.0.
+            Tladder: A list containing the temperature values to use. The values nedd to
+                be monotonically increasing or decreasing. It defaults to ``None``, which
+                will be generated from ``ntemps`` and ``Tmax_ratio``.
+            logprior_fn: A function that evaluate logarithm of the prior distribution.
+                The prior doesn't need to be normalized. It defaults to a uniform prior
+                over a finite range.
+            logprior_args: Additional positional arguments of the ``logprior_fn``. If the
+                default ``logprior_fn`` is used, then the boundaries of the uniform prior
+                can be specified here.
+            **kwargs: Additional keyword arguments for ``ptemcee.Sampler``.
+
+        Attributes:
+            loss: Loss function instance from :class:`~kliff.loss.Loss`
+            T0: Values of the natural temperature, :math:`T_0` [Frederiksen2004]_.
+            Tladder: An array containing the values of the sampling temperatures.
+            sampler: Sampler instance from ``ptemcee.Sampler``
         """
 
         def __init__(
@@ -277,51 +248,37 @@ def _generate_temp_ladder(self, ntemps, Tmax_ratio, Tladder):
     class EmceeSampler(emcee.EnsembleSampler):
         """Sampler class for affine invariant MCMC via ``emcee`` Python package.
 
-        Parameters
-        ----------
-        loss : Loss
-            Loss function instance from :class:`~kliff.loss.Loss`.
-        nwalkers : Optional[int]
-            Number of walkers to simulate. The minimum number of walkers
-            is twice the number of parameters. It defaults to this minimum
-            value.
-        T: Optional[float]
-            Sampling temperatures, used to inflate the likelihood function in the MCMC
-            sampling. It defaults to the natural temperature :math:`T_0` [Frederiksen2004]_.
-        logprior_fn : Optional[Callable]
-            A function that evaluate logarithm of the prior
-            distribution. The prior doesn't need to be normalized. It
-            defaults to a uniform prior over a finite range.
-        logprior_args : Optional[tuple]
-            Additional positional arguments of the ``logprior_fn``. If the
-            default ``logprior_fn`` is used, then the boundaries of the
-            uniform prior can be specified here.
-        **kwargs : Optional[dict]
-            Additional keyword arguments for ``emcee.EnsembleSampler``.
-
-        Attributes
-        ----------
-        loss: Loss
-            Loss function instance from :class:`~kliff.loss.Loss`
-        T: float
-            Values of the sampling temperature.
-        sampler: emcee.EnsembleSampler
-            Sampler instance from ``emcee.EnsembleSampler``
-
-        Notes
-        -----
-        As a convention, KLIFF inflates the likelihood by some sampling temperature, i.e.,
-        :math:`L(\\theta) \propto \exp(-C(\\theta) / T)`. As a default, the sampling
-        temperature is set to the natural temperature. To use the untempered likelihood
-        (:math:`T=1`), user should specify the argument ``T=1``.
-
-
-        References
-        ----------
-        .. [Frederiksen2004] S. L. Frederiksen, K. W. Jacobsen, K. S. Brown, and J. P.
-           Sethna, “Bayesian Ensemble Approach to Error Estimation of Interatomic
-           Potentials,” Phys. Rev. Lett., vol. 93, no. 16, p. 165501, Oct. 2004,
-           doi: 10.1103/PhysRevLett.93.165501.
+        Args:
+            loss: Loss function instance from :class:`~kliff.loss.Loss`.
+            nwalkers: Number of walkers to simulate. The minimum number of walkers is
+                twice the number of parameters. It defaults to this minimum value.
+            T: Sampling temperatures, used to inflate the likelihood function in the MCMC
+                sampling. It defaults to the natural temperature :math:`T_0` [Frederiksen2004]_.
+            logprior_fn: A function that evaluate logarithm of the prior distribution.
+                The prior doesn't need to be normalized. It defaults to a uniform prior
+                over a finite range.
+            logprior_args: Additional positional arguments of the ``logprior_fn``. If the
+                default ``logprior_fn`` is used, then the boundaries of the uniform prior
+                can be specified here.
+            **kwargs: Additional keyword arguments for ``emcee.EnsembleSampler``.
+
+        Attributes:
+            loss: Loss function instance from :class:`~kliff.loss.Loss`
+            T: Values of the sampling temperature.
+            sampler: Sampler instance from ``emcee.EnsembleSampler``
+
+        Notes:
+            As a convention, KLIFF inflates the likelihood by some sampling temperature,
+            i.e., :math:`L(\\theta) \propto \exp(-C(\\theta) / T)`. As a default, the
+            sampling temperature is set to the natural temperature. To use the untempered
+            likelihood (:math:`T=1`), user should specify the argument ``T=1``.
+
+
+        References:
+            .. [Frederiksen2004] S. L. Frederiksen, K. W. Jacobsen, K. S. Brown, and J. P.
+               Sethna, “Bayesian Ensemble Approach to Error Estimation of Interatomic
+               Potentials,” Phys. Rev. Lett., vol. 93, no. 16, p. 165501, Oct. 2004,
+               doi: 10.1103/PhysRevLett.93.165501.
         """
 
         def __init__(
diff --git a/kliff/uq/mcmc_utils.py b/kliff/uq/mcmc_utils.py
index 0f05c531..b2ea4ad8 100644
--- a/kliff/uq/mcmc_utils.py
+++ b/kliff/uq/mcmc_utils.py
@@ -1,4 +1,4 @@
-from typing import Optional
+from typing import Optional, Tuple, Union
 
 import numpy as np
 
@@ -19,30 +19,23 @@ def mser(
     dstep: Optional[int] = 10,
     dmax: Optional[int] = -1,
     full_output: Optional[bool] = False,
-) -> int:
-    """Estimate the equilibration time using marginal standard error rule (MSER). This is
-    done by calculating the standard error (square) of ``chain_d``, where ``chain_d``
-    contains the last :math:`n-d` element of the chain (n is the total number of
-    iterations for each chain), for progresively larger d values, starting from ``dmin``
-    upto ``dmax``, incremented by ``dstep``. The SE values are stored in a list. Then we
-    search the minimum element in the list and return the index of that element.
-
-    Parameters
-    ----------
-    chain: 1D np.ndarray
-        Array containing the time series.
-    dmin: int
-        Index where to start the search in the time series.
-    dstep: int
-        How much to increment the search is done.
-    dmax: int
-        Index where to stop the search in the time series.
-    full_output: bool
-        A flag to return the list of squared standard error.
-
-    Returns
-    -------
-    dstar: int or dict
+) -> Union[int, dict]:
+    """Estimate the equilibration time using marginal standard error rule (MSER).
+
+    This is done by calculating the standard error (square) of ``chain_d``, where
+    ``chain_d`` contains the last :math:`n-d` element of the chain (n is the total number
+    of iterations for each chain), for progresively larger d values, starting from
+    ``dmin`` upto ``dmax``, incremented by ``dstep``. The SE values are stored in a list.
+    Then we search the minimum element in the list and return the index of that element.
+
+    Args:
+        chain: (nsteps,) Array containing the time series.
+        dmin: Index where to start the search in the time series.
+        dstep: How much to increment the search is done.
+        dmax: Index where to stop the search in the time series.
+        full_output: A flag to return the list of squared standard error.
+
+    Returns:
         Estimate of the equilibration time using MSER. If ``full_output=True``, then a
         dictionary containing the estimated equilibration time and the list of squared
         standard errors will be returned.
@@ -65,21 +58,15 @@ def mser(
 
 
 # Estimate autocorrelation length
-def autocorr(chain: np.ndarray, *args, **kwargs):
+def autocorr(chain: np.ndarray, *args, **kwargs) -> np.ndarray:
     """Use ``emcee`` package to estimate the autocorrelation length.
 
-    Parameters
-    ----------
-    chain: np.ndarray (nwalkers, nsteps, ndim,)
-        Chains from the MCMC simulation. The shape of the chains needs to be
-        (nsteps, nwalkers, ndim). Note that the burn-in time needs to be discarded prior
-        to this calculation
-    args, kwargs
-        Additional positional and keyword arguments of ``emcee.autocorr.integrated_time``.
-
-    Returns
-    -------
-    float or array:
+    Args:
+        chain: (nwalkers, nsteps, ndim,) Chains from the MCMC simulation. Note that the
+            burn-in time needs to be discarded prior to this calculation
+        args, kwargs: Additional positional and keyword arguments of ``emcee.autocorr.integrated_time``.
+
+    Returns:
         Estimate of the autocorrelation length for each parameter.
     """
     if emcee_avail:
@@ -90,38 +77,32 @@ def autocorr(chain: np.ndarray, *args, **kwargs):
 
 
 # Assess convergence
-def rhat(chain, time_axis: int = 1, return_WB: bool = False):
-    """Compute the value of :math:`\\hat{r}` proposed by Brooks and Gelman
-    [BrooksGelman1998]_. If the samples come from PTMCMC simulation, then the chain needs
-    to be from one of the temperature only.
-
-    Parameters
-    ----------
-    chain: ndarray
-        The MCMC chain as a ndarray, preferrably with the shape
-        (nwalkers, nsteps, ndims). However, the shape can also be
-        (nsteps, nwalkers, ndims), but the argument time_axis needs to be set
-        to 0.
-    time_axis: int (optional)
-        Axis in which the time series is stored (0 or 1). For emcee results,
-        the time series is stored in axis 0, but for ptemcee for a given
-        temperature, the time axis is 1.
-    return_WB: bool (optional)
-        A flag to return covariance matrices within and between chains.
-
-    Returns
-    -------
-    r: float
-        The value of rhat.
-    W, B: 2d ndarray
-        Matrices of covariance within and between the chains.
-
-    References
-    ----------
-    .. [BrooksGelman1998]
-       Brooks, S.P., Gelman, A., 1998. General Methods for Monitoring Convergence of
-       Iterative Simulations. Journal of Computational and Graphical Statistics 7,
-       434455. https://doi.org/10.1080/10618600.1998.10474787
+def rhat(
+    chain: np.ndarray, time_axis: Optional[int] = 1, return_WB: Optional[bool] = False
+) -> Union[float, Tuple[float, np.ndarray, np.ndarray]]:
+    """Compute the value of :math:`\\hat{r}` proposed by Brooks and Gelman [BrooksGelman1998]_.
+
+    If the samples come from PTMCMC simulation, then the chain needs to be from one of
+    the temperature only.
+
+    Args:
+        chain: The MCMC chain as a ndarray, preferrably with the shape
+            (nwalkers, nsteps, ndim,). However, the shape can also be
+            (nsteps, nwalkers, ndim,), but the argument time_axis needs to be set to 0.
+        time_axis: Axis in which the time series is stored (0 or 1). For emcee results,
+            the time series is stored in axis 0, but for ptemcee for a given temperature,
+            the time axis is 1.
+        return_WB: A flag to return covariance matrices within and between chains.
+
+    Returns:
+        The value of rhat. if ``return_WB=True``, also returns matrices of
+        covariance within and between the chains.
+
+    References:
+        .. [BrooksGelman1998]
+           Brooks, S.P., Gelman, A., 1998. General Methods for Monitoring Convergence of
+           Iterative Simulations. Journal of Computational and Graphical Statistics 7,
+           434455. https://doi.org/10.1080/10618600.1998.10474787
     """
     if time_axis == 1:
         # Reshape the chain so that the time axis is in axis 1
diff --git a/tests/calculators/test_calculator_torch.py b/tests/calculators/test_calculator_torch.py
new file mode 100644
index 00000000..419f72dc
--- /dev/null
+++ b/tests/calculators/test_calculator_torch.py
@@ -0,0 +1,132 @@
+from pathlib import Path
+
+import numpy as np
+import torch
+from torch import Tensor
+
+from kliff import nn
+from kliff.calculators import CalculatorTorch
+from kliff.dataset import Dataset
+from kliff.descriptors import SymmetryFunction
+from kliff.models import NeuralNetwork
+
+# model
+descriptor = SymmetryFunction(
+    cut_name="cos", cut_dists={"Si-Si": 5.0}, hyperparams="set30", normalize=True
+)
+
+N1 = np.random.randint(5, 10)
+N2 = np.random.randint(5, 10)
+model = NeuralNetwork(descriptor)
+model.add_layers(
+    # first hidden layer
+    nn.Linear(descriptor.get_size(), N1),
+    nn.Tanh(),
+    # second hidden layer
+    nn.Linear(N1, N2),
+    nn.Tanh(),
+    # output layer
+    nn.Linear(N2, 1),
+)
+
+# training set
+path = Path(__file__).absolute().parents[1].joinpath("configs_extxyz/Si_4")
+data = Dataset(path)
+configs = data.get_configs()
+
+# calculator
+calc = CalculatorTorch(model, gpu=False)
+_ = calc.create(configs, reuse=False)
+loader = calc.get_compute_arguments(batch_size=100)
+
+# data on parameter sizes
+exp_sizes = [
+    torch.Size([N1, 30]),
+    torch.Size([N1]),
+    torch.Size([N2, N1]),
+    torch.Size([N2]),
+    torch.Size([1, N2]),
+    torch.Size([1]),
+]
+exp_nparams_per_layer = [N1 * 30, N1, N2 * N1, N2, N2, 1]
+exp_nparams_total = np.sum(exp_nparams_per_layer)
+
+# parameters to try
+p0 = np.zeros(exp_nparams_total)
+p1 = np.ones(exp_nparams_total)
+
+
+# Test if the functions to update parameters work
+
+
+def test_get_parameters_sizes():
+    """
+    Test if the function to get parameters sizes works.
+
+    Given the descriptor and the number of nodes in each layer, we can in principle find
+    the expected size of the parameters in each layer.
+    """
+    sizes, nparams_per_layer, nparams_total = calc.get_size_opt_params()
+    assert sizes == exp_sizes, "The sizes retrieved are incorrect"
+    assert (
+        nparams_per_layer == exp_nparams_per_layer
+    ), "The numbers of parameters per layer are incorrect"
+    assert nparams_total == exp_nparams_total, "Total number of parameters is incorrect"
+
+
+def test_parameter_values():
+    """
+    Test if the parameter values are updated.
+
+    This is done by calling `calc.update_model_params` and set the parameter values to
+    zeros. Then, when we retrieve the parameters, we should also get zeros.
+    """
+    # Update the parameters
+    calc.update_model_params(p0)
+    # Check if the parameters retrieved are all zeros
+    assert np.all(
+        calc.get_opt_params() == p0
+    ), "Either `update_model_params` or `get_opt_params` not working"
+
+
+def test_predictions_change():
+    """
+    Test if changing parameters affect the predictions.
+
+    There are two steps of this test. The first one, if we set all the parameters to be
+    zero, then the (forces) predictions should also be zero.
+
+    Then, if we change the parameters to some other values, the predictions should change
+    and they should not be zero, unless there is something special with the
+    configurations.
+    """
+    # Test if predictions are zeros when all parameters are zeros
+    calc.update_model_params(p0)
+    for batch in loader:
+        calc.compute(batch)
+    # We will only look at the forces
+    forces0 = calc.results["forces"]
+    all_zeros = []
+    for f0 in forces0:
+        all_zeros.append(Tensor.all(f0 == 0.0))
+    assert np.all(all_zeros), (
+        "Problem in predicitons calculation: "
+        + "there are non-zero forces when all parameters are zero"
+    )
+
+    # Test if predictions change when we change parameters
+    calc.update_model_params(p1)
+    for batch in loader:
+        calc.compute(batch)
+    forces1 = calc.results["forces"]
+    change = []
+    for f0, f1 in zip(forces0, forces1):
+        change.append(not Tensor.all(f0 - f1 == 0.0))
+    # Use any since there might be special configurations
+    assert np.any(change), "Changing parameters doesn't change predictions"
+
+
+if __name__ == "__main__":
+    test_get_parameters_sizes()
+    test_parameter_values()
+    test_predictions_change()
diff --git a/tests/configs_extxyz/SiC_4/T300_step5000.xyz b/tests/configs_extxyz/SiC_4/T300_step5000.xyz
new file mode 100644
index 00000000..78fde58b
--- /dev/null
+++ b/tests/configs_extxyz/SiC_4/T300_step5000.xyz
@@ -0,0 +1,66 @@
+64
+Lattice="10.862 0 0 0 10.862 0 0 0 10.862 " PBC="1 1 1" Properties=species:S:1:pos:R:3:force:R:3 Energy="-272.98168"
+Si    1.08706e+01  -2.29312e-02   1.09753e+01  1.51360e-01   1.29077e+00  -1.52733e+00
+Si    2.73391e+00   2.70065e+00  -9.06342e-03  6.12039e-01  -5.38713e-01   1.00476e+00
+Si    1.08199e+01   2.78258e+00   2.68842e+00  4.50514e-02  -2.23850e-01   4.51355e-01
+Si    2.73507e+00   4.06389e-02   2.68650e+00  1.21082e-01  -1.83448e-01   9.32376e-01
+C    1.35173e+00   1.41080e+00   1.44068e+00  5.85101e-02  -1.56081e-01  -1.02932e+00
+Si    4.15516e+00   4.11626e+00   1.27052e+00 -7.17150e-01  -1.96152e-01   6.29056e-03
+Si    1.29010e+00   4.08558e+00   4.08204e+00  9.74443e-01   3.78931e-01  -4.65320e-01
+Si    4.03287e+00   1.51702e+00   4.07380e+00  2.41665e-01  -1.78825e+00   7.27809e-03
+Si    5.50415e+00   1.08921e+01   1.09033e+01 -4.59559e-01   2.16510e-01  -5.99524e-01
+Si    8.08230e+00   2.73974e+00   1.09527e+01 -2.35809e-01  -6.83882e-01  -1.05173e+00
+Si    5.40989e+00   2.80707e+00   2.73622e+00  4.84655e-01  -2.29879e-01  -3.77130e-01
+Si    8.16743e+00   1.08698e+01   2.75252e+00  2.97603e-02  -4.54039e-01   4.38043e-01
+Si    6.80158e+00   1.35230e+00   1.46696e+00 -5.10308e-01   5.66508e-01  -2.55501e-01
+Si    9.37403e+00   4.05788e+00   1.40696e+00  7.13992e-01   1.55634e+00   2.00805e-01
+Si    6.77232e+00   4.15275e+00   4.05958e+00 -1.01318e+00   5.19067e-01   1.06870e+00
+C    9.53541e+00   1.41765e+00   4.04694e+00 -1.84262e-01  -6.32007e-01   2.67453e-01
+Si    1.08894e+01   5.36357e+00   1.09343e+01 -9.53305e-01   1.48779e-01  -6.62152e-01
+Si    2.71363e+00   8.08649e+00  -5.10257e-02  2.45951e-01   5.33613e-01   5.70378e-01
+Si    1.28814e-02   8.10487e+00   2.79206e+00 -1.03290e-01   6.97331e-01  -1.83173e+00
+Si    2.71837e+00   5.43762e+00   2.60184e+00  4.58488e-01  -1.00881e+00   1.60455e+00
+Si    1.36878e+00   6.74831e+00   1.37776e+00 -5.62436e-01   9.50964e-01  -5.70933e-01
+Si    4.09700e+00   9.55798e+00   1.31785e+00 -8.10593e-01  -4.46633e-01  -4.65980e-01
+Si    1.40897e+00   9.42544e+00   4.05412e+00  6.32419e-01   5.64328e-01  -2.83858e-01
+Si    4.10437e+00   6.78860e+00   3.98012e+00 -8.11835e-01   3.40268e-01   4.72911e-01
+Si    5.51453e+00   5.43427e+00  -9.65005e-02 -9.93830e-01   1.12890e+00   5.18954e-01
+Si    8.17637e+00   8.19218e+00  -5.02392e-02 -7.05163e-02  -3.16919e-01   1.93546e-01
+Si    5.39754e+00   8.20448e+00   2.65464e+00  6.31803e-01   2.49553e-01   9.23714e-01
+Si    8.03543e+00   5.54799e+00   2.74274e+00  1.13916e+00  -4.08040e-01  -4.64795e-01
+Si    6.75763e+00   6.87600e+00   1.34293e+00  2.15844e-01  -1.83011e-01  -9.80749e-01
+C    9.51024e+00   9.53541e+00   1.33235e+00 -1.59666e-01  -1.17020e+00   1.59728e+00
+Si    6.78385e+00   9.49177e+00   4.15971e+00 -2.23425e-01  -8.74425e-02  -1.08958e+00
+Si    9.46887e+00   6.86346e+00   4.07202e+00 -7.17944e-02  -1.61371e+00   5.93209e-01
+Si    5.84751e-02  -6.96244e-02   5.31462e+00 -4.35583e-01   7.03869e-01   6.18022e-01
+Si    2.66905e+00   2.72471e+00   5.47343e+00 -5.99184e-01   1.89938e-01  -2.35962e-01
+Si    1.08286e+01   2.64356e+00   8.15069e+00  6.21847e-01  -5.83702e-01   1.93889e-01
+Si    2.69535e+00   8.89280e-02   8.17185e+00  8.62088e-02  -1.84017e+00   7.31498e-01
+Si    1.36007e+00   1.22777e+00   6.78619e+00 -8.70852e-01   1.76346e+00  -6.16356e-01
+Si    4.08996e+00   4.07930e+00   6.72275e+00  9.52266e-02  -8.19671e-01   1.24282e+00
+Si    1.35234e+00   3.94019e+00   9.57792e+00 -7.59233e-01   8.33868e-01  -8.85000e-01
+Si    4.09350e+00   1.35894e+00   9.49499e+00  7.23679e-01   3.51873e-01   3.55377e-01
+Si    5.35071e+00  -1.90501e-02   5.42690e+00  3.81758e-02   7.68951e-01  -3.13179e-01
+Si    8.09918e+00   2.74572e+00   5.36865e+00  5.74354e-01  -4.71762e-01   3.67605e-01
+C    5.50996e+00   2.70071e+00   8.12823e+00 -7.11421e-01  -4.34436e-01  -8.98137e-01
+Si    8.16609e+00  -1.57494e-02   8.13037e+00  6.11616e-01  -3.91077e-01  -1.77962e-01
+Si    6.84987e+00   1.25605e+00   6.72079e+00 -8.53129e-01   1.10434e+00   2.26867e-01
+Si    9.53345e+00   3.98061e+00   6.77859e+00 -8.01219e-01   4.83960e-01  -1.17999e-01
+Si    6.78385e+00   4.06483e+00   9.48696e+00  1.20934e+00  -2.38081e-01  -1.84825e-01
+Si    9.42467e+00   1.36164e+00   9.54232e+00  7.38793e-01  -3.21302e-01   7.41018e-01
+Si    1.07950e+01   5.40225e+00   5.44638e+00  5.99360e-01   8.67526e-01  -2.12099e-02
+Si    2.74569e+00   8.10777e+00   5.45158e+00 -8.36755e-03  -3.87478e-01  -3.98498e-01
+C    1.09358e+01   8.07708e+00   8.06821e+00  8.65545e-02   8.63913e-01   9.71933e-01
+Si    2.70664e+00   5.26522e+00   8.17694e+00 -3.64082e-02   3.60222e-01  -1.25221e-01
+Si    1.44321e+00   6.69986e+00   6.83633e+00 -5.91952e-01  -4.50647e-01  -1.06598e+00
+Si    3.98057e+00   9.55196e+00   6.80313e+00  6.37694e-01  -3.77256e-01   2.35405e-01
+Si    1.33122e+00   9.49749e+00   9.53191e+00  2.55721e-01   1.25510e-01  -4.81926e-01
+Si    4.03213e+00   6.73500e+00   9.39925e+00  5.93527e-01  -7.41217e-02   8.77090e-01
+Si    5.36552e+00   5.55381e+00   5.51075e+00  1.45805e+00  -5.74016e-01  -9.49913e-01
+Si    8.11682e+00   8.07091e+00   5.47149e+00  5.96012e-01   1.59499e+00   6.45565e-01
+Si    5.39132e+00   8.17028e+00   8.13755e+00 -2.02935e-01   6.83964e-01   3.80599e-01
+Si    8.20385e+00   5.39076e+00   8.12896e+00  7.73413e-02   6.10681e-02   6.28329e-01
+Si    6.81014e+00   6.82820e+00   6.92450e+00 -3.39664e-01  -7.26523e-01  -1.05084e+00
+Si    9.61890e+00   9.50268e+00   6.77065e+00 -1.25284e+00  -5.28147e-01  -3.82532e-01
+Si    6.79896e+00   9.53872e+00   9.48647e+00 -2.82715e-01  -3.10328e-01  -2.85606e-02
+Si    9.57379e+00   6.81824e+00   9.47305e+00 -1.33267e-01  -1.04954e+00   5.22115e-01
diff --git a/tests/configs_extxyz/SiC_4/T300_step5100.xyz b/tests/configs_extxyz/SiC_4/T300_step5100.xyz
new file mode 100644
index 00000000..b6c8799d
--- /dev/null
+++ b/tests/configs_extxyz/SiC_4/T300_step5100.xyz
@@ -0,0 +1,66 @@
+64
+Lattice="10.862 0 0 0 10.862 0 0 0 10.862 " PBC="1 1 1" Properties=species:S:1:pos:R:3:force:R:3 Energy="-272.89139"
+Si    1.08845e+01   2.52820e-02   1.07067e+01 -1.55709e-01  -9.91907e-01   2.00533e+00
+Si    2.71776e+00   2.76860e+00  -1.32744e-02 -8.31924e-01   1.64169e-01  -7.88431e-02
+Si    1.09379e+01   2.62039e+00   2.73174e+00 -2.54929e-01   7.15361e-01  -4.41133e-01
+Si    2.69115e+00  -2.42743e-02   2.75971e+00  5.95271e-01   1.75718e-01  -4.00328e-01
+Si    1.42098e+00   1.30004e+00   1.29411e+00 -5.80240e-01   2.82915e-01   6.53391e-02
+Si    3.99713e+00   4.01849e+00   1.45235e+00  8.05363e-01   3.41239e-01  -7.29688e-01
+Si    1.43968e+00   4.06293e+00   4.07204e+00 -8.23487e-01   1.56047e-02   2.32339e-01
+C    4.09772e+00   1.26348e+00   4.15078e+00 -4.59120e-01   1.17055e+00  -9.85557e-01
+Si    5.30203e+00   1.08121e+01   1.08010e+01  2.03758e-01  -3.34273e-02   1.25131e+00
+Si    8.15766e+00   2.68522e+00   1.07367e+01  2.15767e-01   6.26227e-01   1.60105e+00
+Si    5.43191e+00   2.60115e+00   2.63080e+00 -3.29769e-02   5.80762e-01   1.55207e+00
+Si    8.15723e+00   1.09088e+01   2.67657e+00 -2.50201e-01  -1.03760e+00  -7.25650e-01
+Si    6.77706e+00   1.33575e+00   1.23187e+00  5.79189e-02  -5.78415e-01  -1.33541e-01
+Si    9.55998e+00   4.03799e+00   1.32652e+00  6.30999e-01  -1.10591e+00  -7.11305e-01
+Si    6.80799e+00   4.04956e+00   4.04657e+00  1.36754e-01  -1.24603e+00  -9.28454e-01
+Si    9.51456e+00   1.29292e+00   4.06791e+00 -1.94886e-01   1.36174e+00   2.65320e-01
+Si    1.09028e+01   5.47673e+00   1.08430e+01 -5.01832e-01  -1.58907e-02   1.29413e-01
+Si    2.76501e+00   8.14566e+00   9.20678e-03 -5.86744e-01   1.37888e-01  -2.30226e-01
+Si   -3.84732e-02   8.18388e+00   2.69552e+00  7.67934e-02  -4.44829e-01   1.11181e+00
+Si    2.74392e+00   5.45465e+00   2.74603e+00 -7.22323e-01   9.41142e-01  -1.81837e-01
+Si    1.31714e+00   6.89960e+00   1.34397e+00  9.43403e-01  -1.33540e+00   3.84611e-01
+Si    4.02086e+00   9.40438e+00   1.45884e+00  9.46069e-01   9.20571e-01  -4.77190e-01
+Si    1.32523e+00   9.56153e+00   4.06581e+00 -7.35982e-01  -1.80193e-01   9.36450e-01
+Si    4.10570e+00   6.80635e+00   4.13803e+00 -4.13159e-02  -7.53051e-01  -6.48603e-02
+Si    5.33588e+00   5.43998e+00   1.02977e-01  4.64495e-01  -1.39599e+00  -6.89798e-01
+C    8.05450e+00   8.01440e+00   3.96543e-02  2.97398e-01   1.11844e+00  -4.82488e-01
+Si    5.48115e+00   8.06435e+00   2.74360e+00 -6.14091e-01   8.90909e-02  -3.64582e-01
+Si    8.24974e+00   5.32983e+00   2.65956e+00 -1.15574e+00   6.51296e-01   8.01292e-01
+Si    6.75686e+00   6.67442e+00   1.38256e+00  5.79101e-01   1.51782e-02   7.33616e-01
+Si    9.47038e+00   9.41757e+00   1.33478e+00 -5.13355e-02   1.60530e+00  -1.10661e+00
+Si    6.81425e+00   9.46199e+00   3.97398e+00  3.14263e-01   5.47434e-01   1.21095e+00
+Si    9.53131e+00   6.70081e+00   4.08706e+00  2.28776e-01   8.35848e-01  -5.26611e-01
+Si   -1.54559e-02   1.19025e-01   5.51970e+00  3.99788e-01  -1.45334e+00  -6.74839e-01
+Si    2.72432e+00   2.71878e+00   5.46968e+00  1.09726e+00   2.62405e-01  -3.29329e-01
+Si    1.09166e+01   2.78752e+00   8.11984e+00 -3.58638e-01   8.38963e-01   1.15483e-01
+Si    2.72101e+00  -4.53199e-02   8.09781e+00  6.87526e-02   1.12171e+00   1.50566e-01
+Si    1.40906e+00   1.43990e+00   6.80092e+00 -4.44491e-01  -8.07655e-01   4.62781e-01
+Si    4.02536e+00   4.09780e+00   6.90435e+00 -3.62960e-01   5.03655e-01  -1.37867e+00
+Si    1.43125e+00   4.23303e+00   9.44403e+00 -5.77487e-01  -1.20483e+00   6.04055e-01
+Si    3.95862e+00   1.36408e+00   9.51456e+00  6.78546e-01  -2.10807e-01  -6.48376e-01
+Si    5.45688e+00  -3.94292e-02   5.42946e+00  8.51698e-01  -4.65555e-01   6.30219e-01
+C    8.14387e+00   2.70161e+00   5.44906e+00  1.39489e-01   1.79348e-02  -6.06797e-02
+Si    5.32666e+00   2.73475e+00   8.16122e+00  8.94350e-01  -2.92816e-01   1.15479e+00
+Si    8.13051e+00   2.91742e-03   8.22023e+00 -3.02150e-01   4.71536e-02  -4.21941e-01
+Si    6.74908e+00   1.39077e+00   6.86706e+00  1.54134e-01  -4.65330e-01  -6.25526e-01
+Si    9.46059e+00   4.17368e+00   6.78752e+00  6.28939e-01  -3.60888e-01   1.88558e-01
+Si    6.76179e+00   4.07666e+00   9.54715e+00 -1.03592e+00   5.35820e-01  -4.00457e-01
+Si    9.55474e+00   1.41418e+00   9.46134e+00 -6.15028e-02  -2.78742e-01  -9.18477e-01
+Si    1.09755e+01   5.49048e+00   5.46197e+00 -9.81819e-01  -1.21302e+00  -5.15269e-01
+Si    2.71459e+00   8.20234e+00   5.42475e+00  2.60725e-01   2.48513e-01   1.43278e-01
+Si    1.08378e+01   8.24858e+00   8.19737e+00 -3.64675e-01  -7.55907e-01  -8.13539e-01
+Si    2.70117e+00   5.61978e+00   8.14061e+00  8.38475e-01  -4.75096e-01  -5.94859e-02
+Si    1.33665e+00   6.89002e+00   6.76893e+00  1.16703e-01   8.35196e-01   7.16284e-01
+Si    4.12901e+00   9.48700e+00   6.85365e+00 -5.33692e-01   2.08416e-01  -1.45737e+00
+Si    1.44729e+00   9.51860e+00   9.47408e+00 -5.14330e-01  -3.19187e-01   3.49375e-01
+C    4.14774e+00   6.82813e+00   9.54229e+00 -4.99405e-01   1.85776e-01  -4.21043e-01
+Si    5.44590e+00   5.29075e+00   5.41548e+00 -6.19045e-01   1.28420e+00   5.71548e-01
+Si    8.14498e+00   8.09044e+00   5.40948e+00  2.09383e-01   5.15762e-01  -7.82864e-01
+Si    5.46611e+00   8.15702e+00   8.11275e+00  2.27699e-01  -9.35350e-01   3.12738e-01
+Si    8.08406e+00   5.51159e+00   8.20811e+00 -1.37995e-01  -1.00705e+00  -9.70478e-01
+Si    6.78248e+00   6.77078e+00   6.70459e+00 -4.84048e-01  -2.95199e-01   7.91892e-01
+Si    9.43142e+00   9.53599e+00   6.75313e+00  4.91186e-01   8.46936e-02   7.20599e-01
+Si    6.71557e+00   9.47026e+00   9.53376e+00  9.40735e-01  -4.68722e-02   8.42826e-02
+Si    9.46240e+00   6.78073e+00   9.49901e+00  7.76993e-01   7.19624e-01   4.59691e-01
diff --git a/tests/configs_extxyz/SiC_4/T300_step5200.xyz b/tests/configs_extxyz/SiC_4/T300_step5200.xyz
new file mode 100644
index 00000000..5d5c1ed5
--- /dev/null
+++ b/tests/configs_extxyz/SiC_4/T300_step5200.xyz
@@ -0,0 +1,66 @@
+64
+Lattice="10.862 0 0 0 10.862 0 0 0 10.862 " PBC="1 1 1" Properties=species:S:1:pos:R:3:force:R:3 Energy="-272.66106"
+Si    1.08500e+01   1.50218e-02   1.09351e+01  1.37951e-01  -3.63970e-01  -6.01600e-01
+Si    2.78990e+00   2.64781e+00   1.52983e-03 -2.24560e-01   5.57739e-01  -1.77788e-01
+Si    1.07368e+01   2.75442e+00   2.75108e+00  8.07987e-01  -3.04998e-01  -3.44536e-01
+Si    2.73246e+00  -3.13199e-02   2.64303e+00 -7.42985e-01   3.96681e-01   1.91314e-01
+Si    1.30717e+00   1.41582e+00   1.35084e+00  4.23929e-01  -2.09336e-01   8.52803e-01
+Si    4.15721e+00   4.13718e+00   1.28701e+00 -5.93303e-01  -6.13105e-01   8.16028e-01
+Si    1.34147e+00   4.05016e+00   4.05116e+00 -3.60850e-01   1.80264e-01  -1.00550e-01
+Si    4.01788e+00   1.38969e+00   3.99910e+00  5.63324e-01  -4.98191e-01   7.55043e-01
+C    5.47260e+00   1.09204e+01   1.08825e+01  1.04538e+00   1.47966e-01  -1.00618e+00
+Si    8.20178e+00   2.75892e+00   1.09000e+01 -5.19996e-01  -2.44467e-01   4.90860e-01
+Si    5.44597e+00   2.80811e+00   2.80388e+00 -3.52853e-01  -3.02536e-01  -1.28776e+00
+Si    8.08822e+00   1.07893e+01   2.81612e+00  6.17405e-01   1.48291e+00  -1.74563e-01
+Si    6.82551e+00   1.41686e+00   1.47388e+00 -1.38338e-01  -1.15831e-01  -5.30033e-01
+Si    9.46729e+00   4.21398e+00   1.39564e+00 -7.44595e-01  -6.16684e-01   5.39416e-01
+Si    6.79011e+00   4.16358e+00   4.06550e+00  9.59780e-01   3.10570e-01   6.71096e-01
+Si    9.51643e+00   1.31122e+00   4.17078e+00 -6.76982e-02  -4.29290e-01  -1.09132e+00
+Si    1.07955e+01   5.34993e+00   1.08007e+01  1.36769e+00   4.00385e-01   3.63401e-01
+C    2.70360e+00   8.20123e+00  -2.70794e-03  1.64992e-03  -1.18661e+00  -3.77687e-01
+Si    3.15399e-02   8.08872e+00   2.73366e+00  1.42304e-01   3.58890e-01  -3.85927e-01
+Si    2.71030e+00   5.41285e+00   2.71716e+00  7.16255e-01  -3.17519e-01  -3.75396e-01
+Si    1.40762e+00   6.67857e+00   1.28321e+00 -1.01684e+00   9.83024e-01   2.49248e-01
+Si    4.10083e+00   9.54890e+00   1.22343e+00 -1.70180e-01  -1.16276e-01   1.32773e+00
+Si    1.36632e+00   9.51490e+00   4.12845e+00  7.71684e-01  -9.21802e-01  -1.79143e+00
+Si    4.07926e+00   6.80732e+00   4.04160e+00  1.31572e-01   2.47988e-01   1.37575e-01
+Si    5.49996e+00   5.47535e+00  -6.21377e-02 -3.68148e-01   5.51398e-01   3.08779e-01
+Si    8.23831e+00   8.21751e+00   7.22925e-04 -2.67836e-01  -4.18642e-01   2.17844e-01
+Si    5.43616e+00   8.23737e+00   2.65392e+00 -9.39410e-01  -6.92007e-01   2.56881e-02
+Si    8.08502e+00   5.66369e+00   2.71861e+00  1.29180e-01  -1.95997e+00  -8.82086e-02
+Si    6.79073e+00   6.90676e+00   1.30578e+00 -1.41962e-01   2.32075e-01  -4.04148e-01
+Si    9.49680e+00   9.56468e+00   1.42794e+00  5.26830e-01  -9.39669e-01   2.96755e-01
+C    6.67433e+00   9.52979e+00   4.08986e+00  2.67495e-01   8.59416e-02   1.92806e-01
+Si    9.46617e+00   6.81709e+00   4.11299e+00  1.52958e-01   9.66763e-01  -1.19885e-01
+Si    8.82243e-02  -5.88158e-02   5.42540e+00 -1.07820e+00   8.29193e-01   8.07927e-01
+Si    2.68050e+00   2.69960e+00   5.39226e+00 -3.24540e-01  -4.69384e-01   1.50724e-01
+Si    1.08246e+01   2.56520e+00   8.18553e+00 -2.44310e-01   4.31206e-01  -5.91504e-01
+Si    2.71053e+00  -2.53842e-02   8.16700e+00 -6.06134e-02  -4.57677e-02  -2.73448e-01
+Si    1.35037e+00   1.31237e+00   6.80342e+00  5.30404e-01   7.52295e-02  -1.93649e-01
+Si    4.06665e+00   4.00750e+00   6.70389e+00  6.36182e-01   2.37870e-01   5.26170e-01
+Si    1.27435e+00   3.88460e+00   9.51839e+00  1.28990e+00   1.11149e+00  -1.02869e-01
+Si    4.12977e+00   1.34987e+00   9.41267e+00 -2.20313e-01  -2.54574e-01   1.04648e+00
+Si    5.37004e+00   9.40110e-03   5.45678e+00 -4.14944e-01   5.17623e-01  -4.20950e-01
+Si    8.21121e+00   2.73062e+00   5.40566e+00 -1.20319e+00   2.71764e-01   3.76665e-01
+C    5.45879e+00   2.72717e+00   8.13075e+00  4.55415e-01   1.93965e-01  -4.16698e-01
+Si    8.21042e+00   1.51257e-03   8.06978e+00 -3.97037e-01  -4.25845e-01   7.07289e-01
+Si    6.80913e+00   1.34596e+00   6.76022e+00  5.88401e-02   4.55501e-01  -1.27984e-01
+Si    9.56343e+00   4.02563e+00   6.79316e+00 -5.48750e-01  -2.94608e-02   3.64358e-01
+Si    6.83620e+00   4.08271e+00   9.53239e+00 -4.73002e-02  -3.03613e-01  -7.62904e-01
+Si    9.46968e+00   1.29213e+00   9.57023e+00 -2.71012e-03   2.12112e-01  -2.09882e-01
+Si    1.07894e+01   5.38890e+00   5.36466e+00  8.13404e-01   1.46918e-02   6.99181e-01
+Si    2.67426e+00   8.11814e+00   5.44518e+00  5.18297e-02  -6.60089e-01   2.28637e-01
+Si    1.09115e+01   8.07910e+00   8.20900e+00 -3.20277e-01   3.09072e-01  -7.75741e-01
+Si    2.75067e+00   5.27998e+00   8.13101e+00 -1.17503e+00   2.94381e-01   5.74669e-01
+Si    1.31150e+00   6.66295e+00   6.81311e+00  6.14627e-01   2.85463e-01  -3.44692e-01
+Si    4.02667e+00   9.47664e+00   6.78571e+00  2.98798e-01  -1.23702e-01   4.29538e-01
+Si    1.33967e+00   9.49917e+00   9.52668e+00 -1.40720e-01   4.99384e-01   2.41528e-01
+Si    4.02096e+00   6.78067e+00   9.48982e+00  3.08738e-02  -4.21122e-01   1.50020e-01
+Si    5.44704e+00   5.47260e+00   5.47833e+00 -8.30015e-01  -3.40489e-01  -9.09213e-01
+Si    8.11628e+00   8.27857e+00   5.38872e+00 -2.66023e-01  -1.56537e+00   1.20915e+00
+Si    5.39589e+00   8.11574e+00   8.18960e+00 -2.95816e-02   9.29825e-01  -9.04100e-02
+Si    8.14247e+00   5.42072e+00   8.10463e+00  4.57692e-01   3.05447e-01   4.29430e-01
+Si    6.71072e+00   6.81334e+00   6.83035e+00  1.64719e+00   5.79584e-01  -6.17128e-01
+Si    9.58860e+00   9.47856e+00   6.78129e+00 -2.64272e-02  -1.88738e-01   1.71331e-01
+Si    6.87705e+00   9.47829e+00   9.48873e+00 -9.10584e-01   9.29916e-01  -1.59752e-01
+Si    9.52926e+00   6.78779e+00   9.53584e+00 -7.58404e-01  -3.07252e-01  -6.95646e-01
diff --git a/tests/configs_extxyz/SiC_4/T300_step5300.xyz b/tests/configs_extxyz/SiC_4/T300_step5300.xyz
new file mode 100644
index 00000000..e2914575
--- /dev/null
+++ b/tests/configs_extxyz/SiC_4/T300_step5300.xyz
@@ -0,0 +1,66 @@
+64
+Lattice="10.862 0 0 0 10.862 0 0 0 10.862 " PBC="1 1 1" Properties=species:S:1:pos:R:3:force:R:3 Energy="-272.42889"
+Si    1.08699e+01  -4.20199e-02   1.07996e+01  2.20512e-02   1.51149e+00  -3.88977e-02
+Si    2.58926e+00   2.78052e+00  -1.39593e-03  8.10141e-01  -1.24719e+00  -2.08758e-01
+Si    1.09340e+01   2.79775e+00   2.68485e+00 -2.80682e-01  -1.35182e+00   2.85280e-01
+Si    2.78443e+00  -2.48951e-03   2.79580e+00 -7.82129e-01   7.88056e-02  -2.03905e-01
+Si    1.35503e+00   1.28984e+00   1.39591e+00  1.88528e-01   4.17038e-01  -4.20778e-01
+Si    4.01058e+00   4.05717e+00   1.31825e+00  1.28453e-01   4.12616e-01   6.79169e-01
+Si    1.31134e+00   4.05657e+00   4.18647e+00  1.35227e+00   8.42998e-02  -9.26991e-01
+Si    4.08316e+00   1.36895e+00   4.14670e+00  5.56300e-01  -1.26922e-01  -2.25283e-01
+Si    5.41222e+00   1.07896e+01   1.09244e+01 -1.09324e+00   1.04236e-01  -4.52515e-01
+Si    8.16418e+00   2.65424e+00   1.09087e+01 -4.82632e-01   5.83148e-03  -2.07760e+00
+Si    5.36683e+00   2.72239e+00   2.70035e+00  4.79537e-01  -1.06699e+00  -4.45521e-01
+C    8.19408e+00   1.08616e+01   2.64395e+00 -2.67026e-01  -3.92955e-01   5.16582e-03
+Si    6.78254e+00   1.27326e+00   1.31466e+00 -4.46977e-01   1.05058e+00   1.01459e+00
+Si    9.51490e+00   3.97871e+00   1.33263e+00  3.25397e-01   1.25051e+00  -9.27610e-02
+Si    6.79497e+00   3.95420e+00   4.06542e+00 -1.42953e+00   3.59371e-01   6.19683e-01
+Si    9.50184e+00   1.39267e+00   3.93446e+00  1.06430e-01   5.03365e-01   1.23346e+00
+Si    1.08311e+01   5.53410e+00   1.09143e+01 -2.19641e-01  -9.44005e-01  -2.20972e-01
+Si    2.67895e+00   8.15965e+00   1.34848e-02  3.66671e-01   7.90716e-01   8.25052e-01
+Si   -4.09161e-02   8.22395e+00   2.69041e+00 -2.39025e-01  -4.46686e-01  -2.48683e-01
+Si    2.70615e+00   5.41372e+00   2.77430e+00 -1.47685e-01   1.87959e-01  -7.38639e-02
+Si    1.35516e+00   6.82529e+00   1.39230e+00  1.04146e-02  -4.35720e-03   5.06950e-01
+Si    4.11047e+00   9.41791e+00   1.53019e+00 -4.72843e-01   3.60300e-01  -1.80520e+00
+Si    1.26876e+00   9.59244e+00   4.06695e+00  7.74796e-01  -3.69443e-01   8.50232e-01
+C    4.13904e+00   6.81374e+00   4.06377e+00 -9.37815e-01  -1.52482e-01  -4.80170e-02
+Si    5.46350e+00   5.40901e+00   5.61699e-02 -5.30576e-01   8.40196e-02  -4.48148e-01
+Si    8.05160e+00   8.04856e+00   5.00546e-02  5.65470e-01   1.00978e+00  -3.45266e-01
+Si    5.52644e+00   8.05991e+00   2.71424e+00  4.71998e-01   6.14675e-01   1.09958e+00
+Si    8.09539e+00   5.18038e+00   2.73493e+00  1.50899e+00   1.65899e+00  -3.32183e-01
+Si    6.85115e+00   6.65866e+00   1.46806e+00 -5.01094e-01  -1.99957e-01  -7.67642e-01
+Si    9.49164e+00   9.54394e+00   1.22918e+00 -4.28349e-01  -1.14523e+00   7.75558e-01
+Si    6.93129e+00   9.47006e+00   4.07182e+00 -5.38255e-01  -3.12750e-01  -7.06112e-01
+Si    9.45474e+00   6.74002e+00   3.97755e+00  6.39942e-01  -7.15496e-01   3.04976e-01
+Si   -9.26225e-02   9.82621e-02   5.46668e+00  1.13188e-01  -5.65146e-01  -1.00007e+00
+Si    2.72276e+00   2.70205e+00   5.50377e+00 -1.91892e-01   2.37605e-01  -4.15003e-01
+Si    1.09006e+01   2.81115e+00   8.09420e+00  4.77521e-01  -3.04592e-01   8.76114e-01
+Si    2.63815e+00   9.25758e-02   8.08371e+00  8.37312e-01  -1.09962e+00   3.12564e-01
+Si    1.29031e+00   1.39708e+00   6.77271e+00  2.56466e-01   4.08914e-01  -8.03143e-03
+Si    4.10800e+00   4.07545e+00   6.78062e+00 -3.21328e-01   6.71510e-01   4.81479e-01
+Si    1.33655e+00   4.19463e+00   9.51522e+00 -3.74252e-01  -2.88609e-03  -2.01544e-01
+Si    4.02017e+00   1.33863e+00   9.50812e+00  1.22600e-02   4.46657e-01   7.11490e-01
+Si    5.58925e+00  -9.90293e-03   5.38052e+00 -1.21168e+00   2.98443e-01   8.25907e-02
+Si    8.12405e+00   2.69049e+00   5.42817e+00  1.03062e+00  -1.90342e-01   1.15034e-02
+Si    5.41311e+00   2.71611e+00   8.17005e+00 -9.40033e-01  -4.69003e-01  -4.56167e-01
+C    8.13746e+00   1.56631e-02   8.18234e+00 -5.08745e-01   8.08146e-02  -6.65697e-01
+Si    6.73839e+00   1.41857e+00   6.78444e+00  1.09425e+00  -9.01860e-01   5.89977e-01
+Si    9.52958e+00   4.04451e+00   6.81241e+00 -5.53528e-01   3.03582e-01  -1.17501e+00
+Si    6.71623e+00   4.04331e+00   9.44917e+00  1.69517e+00   2.23959e-01   1.56232e+00
+Si    9.58375e+00   1.40600e+00   9.38541e+00 -5.82259e-01  -2.65509e-01   1.49451e+00
+Si    1.08565e+01   5.41033e+00   5.44115e+00 -5.67697e-01   1.23435e+00   2.33092e-01
+Si    2.78228e+00   8.21961e+00   5.41274e+00 -8.40222e-01   1.28399e-01  -7.28743e-01
+Si    1.07908e+01   8.19900e+00   8.08132e+00  8.84873e-01  -5.34081e-01   1.27021e+00
+Si    2.76769e+00   5.50318e+00   8.19283e+00 -3.13686e-01   9.02925e-02  -8.08844e-01
+Si    1.37193e+00   6.88303e+00   6.75033e+00  1.20661e-01  -7.80477e-01   3.57705e-01
+Si    4.10913e+00   9.60385e+00   6.72173e+00  7.26362e-02  -1.19276e-01   8.08844e-01
+C    1.35137e+00   9.53981e+00   9.48303e+00  1.38700e-01  -7.01077e-01  -8.33977e-01
+Si    4.07326e+00   6.84492e+00   9.60865e+00  8.06560e-01  -3.96432e-01  -7.79410e-01
+Si    5.43447e+00   5.43639e+00   5.36088e+00  1.05125e+00  -5.88982e-01   8.42410e-01
+Si    8.19853e+00   8.05968e+00   5.42050e+00  4.75263e-01   8.69562e-01  -2.07440e-01
+Si    5.42052e+00   8.16010e+00   8.11935e+00  2.43288e-01  -2.87820e-01  -1.90204e-01
+Si    8.23261e+00   5.36515e+00   8.17716e+00 -1.18266e+00   3.80891e-01  -1.14765e-01
+Si    6.90282e+00   6.74996e+00   6.75362e+00 -1.47083e+00  -4.26810e-01   8.54180e-01
+Si    9.47885e+00   9.51360e+00   6.78234e+00 -7.58140e-01   4.70274e-01  -4.16790e-01
+Si    6.74641e+00   9.49460e+00   9.54064e+00  4.19989e-01  -3.71642e-01  -2.19099e-01
+Si    9.45233e+00   6.81490e+00   9.52829e+00  5.77054e-01   1.52024e-01  -3.78803e-01
diff --git a/tests/uq/test_bootstrap_empirical.py b/tests/uq/test_bootstrap_empirical.py
new file mode 100644
index 00000000..ce70579e
--- /dev/null
+++ b/tests/uq/test_bootstrap_empirical.py
@@ -0,0 +1,200 @@
+from pathlib import Path
+
+import numpy as np
+import pytest
+from scipy.optimize import OptimizeResult
+
+from kliff.calculators.calculator import Calculator, _WrapperCalculator
+from kliff.dataset import Dataset
+from kliff.loss import Loss
+from kliff.models import KIMModel
+from kliff.uq.bootstrap import (
+    Bootstrap,
+    BootstrapEmpiricalModel,
+    BootstrapError,
+    bootstrap_cas_generator_empirical,
+)
+
+seed = 1717
+np.random.seed(seed)
+
+# model
+modelname = "SW_StillingerWeber_1985_Si__MO_405512056662_006"
+model = KIMModel(modelname)
+model.set_opt_params(A=[["default"]])
+
+# training set
+FILE_DIR = Path(__file__).absolute().parent  # Directory of test file
+path = FILE_DIR.parent.joinpath("configs_extxyz/Si_4")
+data = Dataset(path)
+configs = data.get_configs()
+
+# calculators
+# forces
+calc_forces = Calculator(model)
+cas_forces = calc_forces.create(configs, use_energy=False, use_forces=True)
+ncas_forces = len(calc_forces.get_compute_arguments())
+calc_energy = Calculator(model)
+cas_energy = calc_energy.create(configs, use_energy=True, use_forces=False)
+ncas_energy = len(calc_energy.get_compute_arguments())
+calc_comb = _WrapperCalculator(calculators=[calc_energy, calc_forces])
+
+
+# Single calculator
+# loss
+min_kwargs = dict(method="lm")
+loss_forces = Loss(calc_forces)
+loss_forces.minimize(**min_kwargs)
+# initial state
+orig_params = calc_forces.get_opt_params()
+
+# Bootstrap class
+BS_1calc = Bootstrap(loss_forces)
+nsamples = np.random.randint(1, 5)
+
+
+# Multiple calculators
+# loss
+loss_comb = Loss(calc_comb)
+loss_comb.minimize(**min_kwargs)
+
+# Bootstrap class
+BS_2calc = Bootstrap(loss_comb)
+
+
+def test_wrapper():
+    """Test if the Bootstrap class wrapper instantiate the correct class."""
+    assert isinstance(
+        BS_1calc, BootstrapEmpiricalModel
+    ), "Wrapper should instantiate BootstrapEmpiricalModel"
+
+
+def test_error():
+    """Test if BootstrapError is raised when we try to call run before generating
+    bootstrap compute arguments.
+    """
+    with pytest.raises(BootstrapError):
+        BS_1calc.run(min_kwargs=min_kwargs)
+
+
+def test_bootstrap_cas_generator():
+    """Test the shape of generated bootstrap compute arguments."""
+    BS_1calc.generate_bootstrap_compute_arguments(nsamples)
+    bootstrap_cas = BS_1calc.bootstrap_compute_arguments
+    # Test the shape of bootstrap cas samples with default arguments
+    assert (
+        len(bootstrap_cas) == nsamples
+    ), "The number of generated cas is not the same as requested, check the generator"
+    assert np.all(
+        [
+            [len(bs_cas) == ncas_forces for bs_cas in bootstrap_cas[ii]]
+            for ii in range(nsamples)
+        ]
+    ), "For each sample, generator should generate the same number of cas as the original"
+    assert (
+        BS_1calc._nsamples_prepared == nsamples
+    ), "`_nsamples_prepared` property doesn't work"
+
+    # Test the shape of bootstrap cas samples if we specify the number of cas to generate
+    ncas = ncas_forces - 1
+    bootstrap_cas_2 = bootstrap_cas_generator_empirical(
+        nsamples, [cas_forces], ncas=ncas
+    )
+    assert np.all(
+        [
+            [len(bs_cas) == ncas for bs_cas in bootstrap_cas_2[ii]]
+            for ii in range(nsamples)
+        ]
+    ), "Generator doesn't generate the same number of cas as requested for each sample"
+
+
+def test_callback():
+    """Test if callback function works and can break the loop in run method."""
+
+    def callback(_, opt):
+        assert isinstance(opt, OptimizeResult), "Callback cannot capture the run"
+        return True
+
+    BS_1calc.run(min_kwargs=min_kwargs, callback=callback)
+    assert (
+        len(BS_1calc.samples) == 1
+    ), "Callback function cannot break the loop in run method."
+
+
+def test_run():
+    """Test the method to run the calculation."""
+    BS_1calc.run(min_kwargs=min_kwargs)
+    # Test the shape of ensembles
+    samples = BS_1calc.samples
+    shape = samples.shape
+    exp_shape = (nsamples, len(orig_params))
+    assert (
+        shape == exp_shape
+    ), f"Samples doesn't have the right shape; expected {exp_shape}, got {shape}"
+    assert BS_1calc._nsamples_done == shape[0]
+    # Assert restore loss
+    assert np.all(
+        BS_1calc.calculator.get_opt_params() == orig_params
+    ), "Parameters are not restored to initial state"
+
+
+def test_appending_cas():
+    """If we call the generate method again, test if it appends the newly generated
+    bootstrap compute arguments to the previously generated list.
+    """
+    BS_1calc.generate_bootstrap_compute_arguments(nsamples)
+    assert (
+        BS_1calc._nsamples_prepared == 2 * nsamples
+    ), "The newly generated cas is not appended to the old list"
+
+
+def test_save_load_cas():
+    """Test the function to load bootstrap compute arguments."""
+    filename = FILE_DIR / "bootstrap_cas.json"
+    BS_1calc.save_bootstrap_compute_arguments(filename)
+    BS_1calc.load_bootstrap_compute_arguments(filename)
+    filename.unlink()
+    assert (
+        BS_1calc._nsamples_prepared == 4 * nsamples
+    ), "Problem with function to load bootstrap cas"
+
+
+def test_reset():
+    """Test the reset method."""
+    BS_1calc.reset()
+    # Check reset parameters
+    assert np.all(
+        BS_1calc.calculator.get_opt_params() == orig_params
+    ), "Parameters are not restored to initial state"
+    # Check reset bootstrap samples
+    assert BS_1calc._nsamples_prepared == 0, "Reset bootstrap cas failed"
+    assert BS_1calc._nsamples_done == 0, "Reset ensembles failed"
+
+
+def test_multi_calc_cas_generator():
+    """Test the shape of generated bootstrap compute arguments when we use multiple
+    calculators.
+    """
+    BS_2calc.generate_bootstrap_compute_arguments(nsamples)
+    bootstrap_cas = BS_2calc.bootstrap_compute_arguments
+    assert (
+        len(bootstrap_cas) == nsamples
+    ), "The number of generated cas is not the same as requested, check the generator"
+    assert (
+        len(bootstrap_cas[0]) == 2
+    ), "For each sample, the generator should generate cas for each calculator"
+    assert (
+        len(bootstrap_cas[0][0]) + len(bootstrap_cas[0][1]) == ncas_energy + ncas_forces
+    ), "For each sample, generator should generate the same number of cas in total as the original"
+
+
+if __name__ == "__main__":
+    test_wrapper()
+    test_error()
+    test_bootstrap_cas_generator()
+    test_callback()
+    test_run()
+    test_appending_cas()
+    test_save_load_cas()
+    test_reset()
+    test_multi_calc_cas_generator()
diff --git a/tests/uq/test_bootstrap_nn.py b/tests/uq/test_bootstrap_nn.py
new file mode 100644
index 00000000..45d05722
--- /dev/null
+++ b/tests/uq/test_bootstrap_nn.py
@@ -0,0 +1,159 @@
+"""
+TODO:
+- [] Test if the implementation works with model with multiple elements.
+"""
+
+from pathlib import Path
+
+import numpy as np
+import pytest
+
+from kliff import nn
+from kliff.calculators import CalculatorTorch
+from kliff.dataset import Dataset
+from kliff.descriptors import SymmetryFunction
+from kliff.loss import Loss
+from kliff.models import NeuralNetwork
+from kliff.uq.bootstrap import (
+    Bootstrap,
+    BootstrapError,
+    BootstrapNeuralNetworkModel,
+    bootstrap_cas_generator_neuralnetwork,
+)
+
+seed = 1717
+np.random.seed(seed)
+
+# descriptor
+descriptor = SymmetryFunction(
+    cut_name="cos", cut_dists={"Si-Si": 5.0}, hyperparams="set30", normalize=True
+)
+
+# model
+N1 = np.random.randint(5, 10)
+model = NeuralNetwork(descriptor)
+model.add_layers(
+    nn.Linear(descriptor.get_size(), N1),
+    nn.Tanh(),
+    nn.Linear(N1, 1),
+)
+
+# training set
+FILE_DIR = Path(__file__).absolute().parent  # Directory of test file
+path = FILE_DIR.parent.joinpath("configs_extxyz/Si_4")
+data = Dataset(path)
+configs = data.get_configs()
+
+# calculators
+calc = CalculatorTorch(model)
+_ = calc.create(configs, use_energy=False, use_forces=True)
+fingerprints = calc.get_fingerprints()
+nfingerprints = len(fingerprints)
+
+loss = Loss(calc)
+min_kwargs = dict(method="Adam", num_epochs=10, batch_size=100, lr=0.001)
+result = loss.minimize(**min_kwargs)
+
+orig_state_filename = FILE_DIR / "orig_model.pkl"
+BS = Bootstrap(loss, orig_state_filename=orig_state_filename)
+nsamples = np.random.randint(1, 5)
+
+
+def test_wrapper():
+    """Test if the Bootstrap class wrapper instantiate the correct class."""
+    assert isinstance(
+        BS, BootstrapNeuralNetworkModel
+    ), "Wrapper should instantiate BootstrapNeuralNetworkModel"
+
+
+def test_error():
+    """Test if BootstrapError is raised when we try to call run before generating
+    bootstrap compute arguments.
+    """
+    with pytest.raises(BootstrapError):
+        BS.run(min_kwargs=min_kwargs)
+
+
+def test_original_state():
+    """Test we export the original state when we instantiate the class."""
+    assert orig_state_filename.exists(), "Original state is not exported"
+
+
+def test_bootstrap_cas_generator():
+    """Test the generator function generate the same number of bootstrap compute arguments
+    samples as requested.
+    """
+    # Test the shape of bootstrap cas samples with default arguments
+    BS.generate_bootstrap_compute_arguments(nsamples)
+    bootstrap_cas = BS.bootstrap_compute_arguments
+    assert (
+        len(bootstrap_cas) == nsamples
+    ), "The number of generated cas is not the same as requested, check the generator"
+    assert np.all(
+        [len(bs_cas) == nfingerprints for _, bs_cas in bootstrap_cas.items()]
+    ), "For each sample, generator should generate the same number of cas as the original"
+    assert (
+        BS._nsamples_prepared == nsamples
+    ), "`_nsamples_prepared` property doesn't work"
+
+    # Test the shape of bootstrap cas samples if we specify the number of cas to generate
+    nfp = nfingerprints - 1
+    bootstrap_cas_2 = bootstrap_cas_generator_neuralnetwork(
+        nsamples, fingerprints, nfingerprints=nfp
+    )
+    assert np.all(
+        [len(bs_cas) == nfp for _, bs_cas in bootstrap_cas_2.items()]
+    ), "Generator doesn't generate the same number of cas as requested for each sample"
+
+
+def test_run():
+    """Test the method to run the calculation."""
+    BS.run(min_kwargs=min_kwargs)
+    # Test the shape of ensembles
+    samples = BS.samples
+    shape = samples.shape
+    exp_shape = (nsamples, calc.get_size_opt_params()[-1])
+    assert (
+        shape == exp_shape
+    ), f"Samples doesn't have the right shape; expected {exp_shape}, got {shape}"
+    assert BS._nsamples_done == nsamples
+
+
+def test_appending_cas():
+    """If we call the generate method again, test if it appends the newly generated
+    bootstrap compute arguments to the previously generated list.
+    """
+    BS.generate_bootstrap_compute_arguments(nsamples)
+    assert (
+        BS._nsamples_prepared == 2 * nsamples
+    ), "The newly generated cas is not appended to the old list"
+
+
+def test_save_load_cas():
+    """Test the function to load bootstrap compute arguments."""
+    filename = FILE_DIR / "bootstrap_cas.json"
+    BS.save_bootstrap_compute_arguments(filename)
+    BS.load_bootstrap_compute_arguments(filename)
+    filename.unlink()
+    assert (
+        BS._nsamples_prepared == 4 * nsamples
+    ), "Problem with function to load bootstrap cas"
+
+
+def test_reset():
+    """Test the reset method."""
+    BS.reset()
+    # Check reset bootstrap samples
+    assert BS._nsamples_prepared == 0, "Reset bootstrap cas failed"
+    assert BS._nsamples_done == 0, "Reset ensembles failed"
+
+
+if __name__ == "__main__":
+    test_wrapper()
+    test_error()
+    test_original_state()
+    test_bootstrap_cas_generator()
+    test_run()
+    test_appending_cas()
+    test_save_load_cas()
+    test_reset()
diff --git a/tests/uq/test_bootstrap_nn_separate_species.py b/tests/uq/test_bootstrap_nn_separate_species.py
new file mode 100644
index 00000000..e839df71
--- /dev/null
+++ b/tests/uq/test_bootstrap_nn_separate_species.py
@@ -0,0 +1,109 @@
+import os
+from pathlib import Path
+
+import numpy as np
+
+from kliff import nn
+from kliff.calculators.calculator_torch import CalculatorTorchSeparateSpecies
+from kliff.dataset import Dataset
+from kliff.descriptors import SymmetryFunction
+from kliff.loss import Loss
+from kliff.models import NeuralNetwork
+from kliff.uq.bootstrap import BootstrapNeuralNetworkModel
+
+seed = 1717
+np.random.seed(seed)
+
+# descriptor
+descriptor = SymmetryFunction(
+    cut_name="cos",
+    cut_dists={"Si-Si": 5.0, "C-C": 5.0, "Si-C": 5.0},
+    hyperparams="set30",
+    normalize=True,
+)
+
+# models
+# Si
+N1 = np.random.randint(5, 10)
+model_si = NeuralNetwork(descriptor)
+model_si.add_layers(
+    nn.Linear(descriptor.get_size(), N1),
+    nn.Tanh(),
+    nn.Linear(N1, 1),
+)
+
+# C
+model_c = NeuralNetwork(descriptor)
+model_c.add_layers(
+    nn.Linear(descriptor.get_size(), N1),
+    nn.Tanh(),
+    nn.Linear(N1, 1),
+)
+
+# training set
+FILE_DIR = Path(__file__).absolute().parent  # Directory of test file
+path = FILE_DIR.parent.joinpath("configs_extxyz/SiC_4")
+data = Dataset(path)
+configs = data.get_configs()
+
+# calculators
+calc = CalculatorTorchSeparateSpecies({"Si": model_si, "C": model_c})
+_ = calc.create(configs, use_energy=False, use_forces=True)
+
+loss = Loss(calc)
+min_kwargs = dict(method="Adam", num_epochs=10, batch_size=100, lr=0.001)
+result = loss.minimize(**min_kwargs)
+
+orig_state_filename = FILE_DIR / "orig_model.pkl"
+BS = BootstrapNeuralNetworkModel(loss, orig_state_filename=orig_state_filename)
+nsamples = np.random.randint(1, 5)
+
+
+def test_model():
+    """Test the model property."""
+    # models
+    assert len(BS.model) == 2, "There should be 2 elements in BS.model"
+    # species
+    assert np.all(
+        [sp in ["Si", "C"] for sp in BS._species]
+    ), "There are species not recorded"
+
+
+def test_original_state():
+    """Test we export the original states for all models when we instantiate the class."""
+    # check shape
+    assert (
+        len(BS.orig_state_filename) == 2
+    ), "There should be 2 elements in orig_state_filename"
+    # check elements
+    splitted_path = os.path.splitext(orig_state_filename)
+    fnames = [
+        Path(splitted_path[0] + f"_{el}" + splitted_path[1]) for el in ["Si", "C"]
+    ]
+    assert np.all(
+        [str(fname) in BS.orig_state_filename for fname in fnames]
+    ), "Not all original state filename are listed"
+    # check if initial states are exported
+    assert np.all(
+        [fname.exists() for fname in fnames]
+    ), "Not all original state is not exported"
+
+
+def test_run():
+    """Test the method to run the calculation."""
+    BS.generate_bootstrap_compute_arguments(nsamples)
+    BS.run(min_kwargs=min_kwargs)
+    # Test the shape of ensembles
+    samples = BS.samples
+    shape = samples.shape
+    exp_shape = (nsamples, calc.get_size_opt_params()[-1])
+    assert (
+        shape == exp_shape
+    ), f"Samples doesn't have the right shape; expected {exp_shape}, got {shape}"
+    assert BS._nsamples_done == nsamples
+
+
+if __name__ == "__main__":
+    test_model()
+    test_original_state()
+    test_run()
diff --git a/tests/uq/test_mcmc.py b/tests/uq/test_mcmc.py
index fbec073c..60038d3f 100644
--- a/tests/uq/test_mcmc.py
+++ b/tests/uq/test_mcmc.py
@@ -116,9 +116,9 @@ def test_dimensionality():
     if emcee_avail:
         p0 = np.random.uniform(0, 10, (nwalkers, ndim))
         sampler.run_mcmc(initial_state=p0, nsteps=nsteps)
-        assert sampler.chain.shape == (
-            nwalkers,
+        assert sampler.get_chain().shape == (
             nsteps,
+            nwalkers,
             ndim,
         ), "Dimensionality from the emcee wrapper is not right"
     else: