diff --git a/tools/tabulate/pair_bi_tabulate.py b/tools/tabulate/pair_bi_tabulate.py
new file mode 100755
index 0000000000..910bedf64e
--- /dev/null
+++ b/tools/tabulate/pair_bi_tabulate.py
@@ -0,0 +1,117 @@
+#!/usr/bin/env python3
+
+from tabulate import PairTabulate
+import sys
+import argparse
+import numpy as np
+from scipy.signal import savgol_filter
+from scipy.optimize import curve_fit
+
+"""
+   This script gives an example on how to make tabulated forces from radial
+   distribution function using tabulate.py.
+   Required: python3, numpy, scipy.
+   BI stands for Boltzmann Inversion.
+"""
+###############################################################################
+
+
+class BI(PairTabulate):
+    def __init__(self, units=None, comment=None, T=1):
+        super(PairTabulate, self).__init__("pair", self.energy, units, comment)
+        self.parser.add_argument(
+            "--eshift",
+            "-e",
+            dest="eshift",
+            default=False,
+            action="store_true",
+            help="Shift potential energy to be zero at outer cutoff",
+        )
+        self.parser.add_argument(
+            "--rdffile", default="rdf.dat", help="Rdf file to be read."
+        )
+        try:
+            self.args = self.parser.parse_args()
+        except argparse.ArgumentError:
+            sys.exit()
+
+        kb = 1
+        # Add more kb units if you need
+        if units == "si":
+            kb = 1.380649e-23  # J/K
+        elif units == "metal":
+            kb = 8.617333e-5  # eV/K
+        elif units == "real":
+            kb = 1.987204e-3  # kcal/mol/K
+        else:
+            sys.stdout.write("WARNING: Unknown or lj units, using kb=1\n")
+        self.kbT = kb * T
+        self.r, self.e, self.f = self.read_rdf(self.args.rdffile)
+
+    # This function assumes LAMMPS format for rdf with a single entry
+    def read_rdf(self, rdffile):
+        data = np.loadtxt(rdffile, skiprows=4)
+        r = data[:, 1]
+        g = data[:, 2]
+
+        # savgol_filter is an example of smoothing.
+        # Other filters/functions can be used.
+        g = savgol_filter(g, 10, 5)
+        return self.inversion(r, g)
+
+    def inversion(self, r, g):
+        r = r
+        e = -self.kbT * np.log(g)
+        e = self.complete_exponential(r, e)
+        f = -np.gradient(e, r)
+        return r, e, f
+
+    def complete_exponential(self, r, e):
+        r_temp = r[e != np.inf]
+        e_temp = e[e != np.inf]
+
+        # Optimising the parameter for a function for derivation
+        # to be continuous.
+        # Here a gaussian function, can be anything relevant defined in func.
+        popt, pcov = curve_fit(self.func, r_temp[:2], e_temp[:2])
+        for i, _ in enumerate(e):
+            if e[i] == np.inf:
+                e[i] = self.func(r[i], *popt)
+        return e
+
+    def func(self, x, K, s):
+        return K * np.exp(-0.5 * (x / s) ** 2) / (s * np.sqrt(2 * np.pi))
+
+    def energy(self, x):
+        e = self.e
+        r = self.r
+        # Force estimation at minimum distance.
+        # Should not be that useful
+        f0 = (e[1] - e[0]) / (r[1] - r[0])
+
+        minr = min(r)
+        maxr = max(r)
+        # Note that you might want OOB to return an error.
+        if x >= maxr:
+            return 0
+        if x < minr:
+            dx = minr - x
+            return -f0 * dx
+        else:
+            # Linear interpolation between points.
+            for i, ri in enumerate(r):
+                if r[i] < x:
+                    r1, e1 = r[i], e[i]
+                    r2, e2 = r[i + 1], e[i + 1]
+            dr12 = r2 - r1
+            dr = x - r1
+            de = (e2 - e1) / (r2 - r1)
+            return e1 + (de * dr / dr12)
+
+
+###############################################################################
+
+
+if __name__ == "__main__":
+    ptable = BI()
+    ptable.run("BI")