diff --git a/tools/README b/tools/README
index b611b6c97f..853ae50c6b 100644
--- a/tools/README
+++ b/tools/README
@@ -23,6 +23,7 @@ eam_database	       one tool to generate EAM alloy potential files
 eam_generate	       2nd tool to generate EAM alloy potential files
 eff		       scripts for working with the eFF (electron force field)
 emacs		       add-ons to EMACS editor for editing LAMMPS input scripts
+fep		       scripts for free-energy perturbation with USER-FEP pkg
 ipp		       input pre-processor Perl tool for creating input scripts
 lmp2arc		       convert LAMMPS output to Accelrys Insight format
 lmp2cfg		       convert LAMMPS output to CFG files for AtomEye viz
diff --git a/tools/fep/README b/tools/fep/README
new file mode 100644
index 0000000000..0663963752
--- /dev/null
+++ b/tools/fep/README
@@ -0,0 +1,7 @@
+These are utility scripts provided as part of the USER-FEP package for
+free energy perturbation simulations with soft-core pair potentials in
+LAMMPS.
+
+The person who created these tools is Agilio Padua at Université
+Blaise Pascal Clermont-Ferrand (agilio.padua at univ-bpclermont.fr)
+Contact him directly if you have questions.
diff --git a/tools/fep/bar.py b/tools/fep/bar.py
new file mode 100755
index 0000000000..a2c51a7164
--- /dev/null
+++ b/tools/fep/bar.py
@@ -0,0 +1,88 @@
+#!/usr/bin/env python
+# bar.py - Bennet's acceptance ratio method for free energy calculation
+
+import sys, math
+
+if len(sys.argv) < 4:
+    print "Bennet's acceptance ratio method"
+    print "usage: bar.py temperature datafile01 datafile10 [delf_lo delf_hi]"
+    print "  datafile01 contains (U_1 - U_0)_0 in 2nd column"
+    print "  datafile10 contains (U_0 - U_1)_1 in 2nd column"
+    print "  (first column is index, time step, etc. and is ignored)"
+    print "  delf_lo and delf_hi are optional guesses bracketing the solution"
+    sys.exit()
+
+if len(sys.argv) == 6:
+    delf_lo = float(sys.argv[4])
+    delf_hi = float(sys.argv[5])
+else:
+    delf_lo = -50.0
+    delf_hi =  50.0
+    
+def fermi(x):
+    if x > 100:
+        return 0.0
+    else:
+        return 1.0 / (1.0 + math.exp(x))
+
+def avefermi(eng, delf):
+    ave = 0.0
+    n = 0
+    for du in eng:
+        ave += fermi((du + delf) / rt)
+        n += 1
+    return ave / n
+
+def bareq(delf):
+    ave0 = avefermi(eng01, -delf)
+    ave1 = avefermi(eng10,  delf)
+    return ave1 - ave0
+
+def bisect(func, xlo, xhi, xtol = 1.0e-4, maxit = 20):
+    if xlo > xhi:
+        aux = xhi
+        xhi = xlo
+        xlo = aux
+    if func(xlo) * func(xhi) > 0.0:
+        print("error: root not bracketed by interval")
+        sys.exit(2)
+    for i in range(maxit):
+        sys.stdout.write('.')
+        sys.stdout.flush()
+        xmid = (xlo + xhi) / 2.0
+        if func(xlo) * func(xmid) < 0.0:
+            xhi = xmid
+        else:
+            xlo = xmid
+        if xhi - xlo < xtol:
+            return xmid
+    return xmid
+
+print "Bennet's acceptance ratio method"
+print sys.argv[1], " K"
+rt = 0.008314 / 4.184 * float(sys.argv[1])
+
+eng01 = []                                # read datafiles
+with open(sys.argv[2], 'r') as f:
+    for line in f:
+        if line.startswith('#'):
+            continue
+        eng01.append(float(line.split()[1]))
+
+eng10 = []
+with open(sys.argv[3], 'r') as f:
+    for line in f:
+        if line.startswith('#'):
+            continue
+        eng10.append(float(line.split()[1]))
+
+sys.stdout.write("solving")
+delf = bisect(bareq, delf_lo, delf_hi)
+print "."
+
+ave0 = avefermi(eng01, -delf)
+ave1 = avefermi(eng10,  delf)
+
+print "<...>0 = ", ave0
+print "<...>1 = ", ave1
+print "deltaA = ", delf
diff --git a/tools/fep/fdti.py b/tools/fep/fdti.py
new file mode 100755
index 0000000000..c17eb5a037
--- /dev/null
+++ b/tools/fep/fdti.py
@@ -0,0 +1,36 @@
+#!/usr/bin/env python
+# fdti.py - integrate compute fep results using the trapezoidal rule
+
+import sys, math
+
+if len(sys.argv) < 3:
+    print "Finite Difference Thermodynamic Integration (Mezei 1987)"
+    print "Trapezoidal integration of compute_fep results at equally-spaced points"
+    print "usage: fdti.py temperature hderiv < fep.lmp"
+    sys.exit()
+
+rt = 0.008314 / 4.184 * float(sys.argv[1])
+hderiv = float(sys.argv[2])
+
+line = sys.stdin.readline()
+while line.startswith("#"):
+    line = sys.stdin.readline()
+
+v = 1.0
+tok = line.split()
+if len(tok) == 4:
+    v = float(tok[3])
+lo = -rt * math.log(float(tok[2]) / v)
+    
+i = 0
+sum = 0.0
+for line in sys.stdin:
+    tok = line.split()
+    if len(tok) == 4:
+        v = float(tok[3])
+    hi = - rt * math.log(float(tok[2]) / v)
+    sum += (hi + lo) / (2 * hderiv)
+    lo = hi
+    i += 1
+
+print sum / i      # int_0^1: divide by i == multiply by delta
diff --git a/tools/fep/fep.py b/tools/fep/fep.py
new file mode 100755
index 0000000000..3ed9290ea4
--- /dev/null
+++ b/tools/fep/fep.py
@@ -0,0 +1,23 @@
+#!/usr/bin/env python
+# fep.py - calculate free energy from compute fep results
+
+import sys, math
+
+if len(sys.argv) < 2:
+    print "Free Energy Perturbation"
+    print "usage: fep.py temperature < fep.lmp"
+    sys.exit()
+
+rt = 0.008314 / 4.184 * float(sys.argv[1])
+
+v = 1.0
+sum = 0.0
+for line in sys.stdin:
+    if line.startswith("#"):
+        continue
+    tok = line.split()
+    if len(tok) == 4:
+        v = float(tok[3])
+    sum += math.log(float(tok[2]) / v)
+
+print -rt * sum
diff --git a/tools/fep/nti.py b/tools/fep/nti.py
new file mode 100755
index 0000000000..716dea7cb3
--- /dev/null
+++ b/tools/fep/nti.py
@@ -0,0 +1,30 @@
+#!/usr/bin/env python
+# nti.py - integrate compute fep results using the trapezoidal rule
+
+import sys, math
+
+if len(sys.argv) < 3:
+    print "Thermodynamic Integration with Numerical Derivative"
+    print "Trapezoidal integration of compute_fep results at equally-spaced points"
+    print "usage: nti.py temperature hderiv < fep.lmp"
+    sys.exit()
+
+hderiv = float(sys.argv[2])
+
+line = sys.stdin.readline()
+while line.startswith("#"):
+    line = sys.stdin.readline()
+
+tok = line.split()
+lo = float(tok[1])
+
+i = 0
+sum = 0.0
+for line in sys.stdin:
+    tok = line.split()
+    hi = float(tok[1])
+    sum += (hi + lo) / (2 * hderiv)
+    lo = hi
+    i += 1
+
+print sum / i      # int_0^1: divide by i == multiply by delta