From db78f042efe3108e3147d04857f2de5c9da9492b Mon Sep 17 00:00:00 2001
From: sjplimp <sjplimp@f3b2605a-c512-4ea7-a41b-209d697bcdaa>
Date: Fri, 9 May 2014 15:34:57 +0000
Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@11955
f3b2605a-c512-4ea7-a41b-209d697bcdaa
---
tools/README | 1 +
tools/fep/README | 7 ++++
tools/fep/bar.py | 88 +++++++++++++++++++++++++++++++++++++++++++++++
tools/fep/fdti.py | 36 +++++++++++++++++++
tools/fep/fep.py | 23 +++++++++++++
tools/fep/nti.py | 30 ++++++++++++++++
6 files changed, 185 insertions(+)
create mode 100644 tools/fep/README
create mode 100755 tools/fep/bar.py
create mode 100755 tools/fep/fdti.py
create mode 100755 tools/fep/fep.py
create mode 100755 tools/fep/nti.py
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
--
GitLab