Start update of fep examples and doc

examples/USER/fep/CH4hyd/fdti01/fdti01.fep

@@ -0,0 +1,23 @@
+# Time-averaged data for fix FEP
+# TimeStep c_FEP[1] c_FEP[2]
+100000 0.000675923 0.998872
+200000 0.00676971 0.988725
+300000 0.012229 0.979748
+400000 0.0195902 0.967815
+500000 0.0264713 0.956934
+600000 0.0311625 0.949904
+700000 0.0149588 0.975649
+800000 0.00933401 0.984811
+900000 0.00651534 0.989335
+1000000 0.00224647 0.996372
+1100000 0.00227072 0.996327
+1200000 -0.000223228 1.00045
+1300000 -0.000553474 1.00099
+1400000 -0.00183338 1.00313
+1500000 -0.00231784 1.00394
+1600000 -0.00338448 1.00572
+1700000 -0.00364256 1.00615
+1800000 -0.00432816 1.0073
+1900000 -0.00528456 1.00891
+2000000 -0.00581585 1.00981
+2100000 -0.00657668 1.0111

examples/USER/fep/CH4hyd/fdti01/fdti01.lmp

@@ -1,22 +0,0 @@
-# Time-averaged data for fix FEP
-# TimeStep c_FEP[1] c_FEP[2]
-100000 0.000289416 0.999521
-200000 0.00590502 0.990163
-300000 0.0115179 0.980934
-400000 0.0216052 0.96457
-500000 0.0222451 0.963797
-600000 0.0200038 0.967603
-700000 0.0152292 0.975176
-800000 0.00896315 0.985384
-900000 0.00585213 0.990434
-1000000 0.00327599 0.99467
-1100000 0.00159845 0.997437
-1200000 0.000171108 0.999804
-1300000 -0.000313183 1.00059
-1400000 -0.00148 1.00254
-1500000 -0.00297976 1.00504
-1600000 -0.00287926 1.00487
-1700000 -0.00430671 1.00727
-1800000 -0.00453729 1.00765
-1900000 -0.00507356 1.00856
-2000000 -0.00586954 1.0099

examples/USER/fep/CH4hyd/fdti01/in.fdti01.lmp

@@ -72,12 +72,24 @@ compute FEP all fep ${TK} &
   atom charge 1 v_dq1 &
   atom charge 2 v_dq2
 
-fix FEP all ave/time 20 4000 100000 c_FEP[1] c_FEP[2] file fdti01.lmp
+fix FEP all ave/time 20 4000 100000 c_FEP[1] c_FEP[2] file fdti01.fep
 
 dump TRAJ all custom 20000 dump.lammpstrj id mol type element x y z ix iy iz
 dump_modify TRAJ element C H O H
 
 run 2000000
 
+unfix ADAPT
+
+pair_coeff    1    3  lj/cut/tip4p/long/soft   0.1036   3.3279  1.0  # C4H Ow
+pair_coeff    1    4  lj/cut/tip4p/long/soft   0.0000   1.0000  1.0  # C4H Hw
+pair_coeff    2    3  lj/cut/tip4p/long/soft   0.0699   2.8126  1.0  # H Ow
+pair_coeff    2    4  lj/cut/tip4p/long/soft   0.0000   1.0000  1.0  # H Hw
+
+set type 1 charge -0.24
+set type 2 charge 0.06
+
+run 100000
+
 # write_restart restart.*.lmp
 write_data data.*.lmp

tools/fep/README

@@ -2,6 +2,6 @@ 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)
+The person who created these tools is Agilio Padua at ENS de Lyon
+(agilio.padua at ens-lyon.fr)
 Contact him directly if you have questions.

tools/fep/bar.py

@@ -1,15 +1,16 @@
 #!/usr/bin/env python
 # bar.py - Bennet's acceptance ratio method for free energy calculation
 
-import sys, math
+import sys
+import 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"
+    print("Bennet 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:
@@ -58,8 +59,8 @@ def bisect(func, xlo, xhi, xtol = 1.0e-4, maxit = 20):
             return xmid
     return xmid
 
-print "Bennet's acceptance ratio method"
-print sys.argv[1], " K"
+print("Bennet acceptance ratio method")
+print(sys.argv[1], " K")
 rt = 0.008314 / 4.184 * float(sys.argv[1])
 
 eng01 = []                                # read datafiles
@@ -78,11 +79,11 @@ with open(sys.argv[3], 'r') as f:
 
 sys.stdout.write("solving")
 delf = bisect(bareq, delf_lo, delf_hi)
-print "."
+print(".")
 
 ave0 = avefermi(eng01, -delf)
 ave1 = avefermi(eng10,  delf)
 
-print "<...>0 = ", ave0
-print "<...>1 = ", ave1
-print "deltaA = ", delf
+print("<...>0 = ", ave0)
+print("<...>1 = ", ave1)
+print("deltaA = ", delf)

tools/fep/fdti.py

@@ -1,12 +1,13 @@
 #!/usr/bin/env python
 # fdti.py - integrate compute fep results using the trapezoidal rule
 
-import sys, math
+import sys
+import 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"
+    print("Finite Difference Thermodynamic Integration (Mezei 1987)")
+    print("Trapezoidal integration of compute_fep results at equally-spaced points")
+    print("usage: fdti.py temperature hderiv < out.fep")
     sys.exit()
 
 rt = 0.008314 / 4.184 * float(sys.argv[1])
@@ -33,4 +34,4 @@ for line in sys.stdin:
     lo = hi
     i += 1
 
-print sum / i      # int_0^1: divide by i == multiply by delta
+print(sum/(i - 1))    # int_0^1: divide by i - 1 == multiply by delta

tools/fep/fep.py

@@ -1,11 +1,12 @@
 #!/usr/bin/env python
 # fep.py - calculate free energy from compute fep results
 
-import sys, math
+import sys
+import math
 
 if len(sys.argv) < 2:
-    print "Free Energy Perturbation"
-    print "usage: fep.py temperature < fep.lmp"
+    print("Free Energy Perturbation")
+    print("usage: fep.py temperature < out.fep")
     sys.exit()
 
 rt = 0.008314 / 4.184 * float(sys.argv[1])
@@ -20,4 +21,4 @@ for line in sys.stdin:
         v = float(tok[3])
     sum += math.log(float(tok[2]) / v)
 
-print -rt * sum
+print(-rt * sum)

tools/fep/nti.py

@@ -1,12 +1,13 @@
 #!/usr/bin/env python
 # nti.py - integrate compute fep results using the trapezoidal rule
 
-import sys, math
+import sys
+import 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"
+    print("Thermodynamic Integration with Numerical Derivative")
+    print("Trapezoidal integration of compute_fep results at equally-spaced points")
+    print("usage: nti.py temperature hderiv < out.fep")
     sys.exit()
 
 hderiv = float(sys.argv[2])
@@ -27,4 +28,4 @@ for line in sys.stdin:
     lo = hi
     i += 1
 
-print sum / i      # int_0^1: divide by i == multiply by delta
+print(sum/(i - 1))    # int_0^1: divide by i - 1 == multiply by delta