diff --git a/tools/README b/tools/README index b611b6c97..853ae50c6 100644 --- a/tools/README +++ b/tools/README @@ -1,47 +1,48 @@ LAMMPS pre- and post-processing tools This directory contains several stand-alone tools for creating LAMMPS input files and massaging LAMMPS output data. Instructions on how to use the tools are discussed in the "Additional Tools" section of the LAMMPS documentation. Tools that are single source files in this directory have additional comments that may be useful at the top of the source file. Tools that reside in their own sub-directories have README files you should look at. These are the included tools: amber2lmp python scripts for using AMBER to setup LAMMPS input binary2txt convert a LAMMPS dump file from binary to ASCII text ch2lmp convert CHARMM files to LAMMPS input chain create a data file of bead-spring chains colvars post-process output of the fix colvars command createatoms generate lattices of atoms within a geometry data2xmovie convert a data file to a snapshot that xmovie can viz 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 lmp2traj convert LAMMPS output to contour, density profiles lmp2vmd tools for visualizing and analyzing LAMMPS data with VMD matlab MatLab scripts for post-processing LAMMPS output micelle2d create a data file of small lipid chains in solvent msi2lmp use Accelrys Insight code to setup LAMMPS input phonon post-process output of the fix phonon command pymol_asphere convert LAMMPS output of ellipsoids to PyMol format python Python scripts for post-processing LAMMPS output reax Tools for analyzing output of ReaxFF simulations vim add-ons to VIM editor for editing LAMMPS input scripts xmovie a quick/simple viz package (2d projections of 3d) xmgrace a collection of scripts to generate xmgrace plots For tools that are single C, C++, or Fortran files, a Makefile for building them is included in this directory. You may need to edit it for the compilers and paths on your system. For tools in their own sub-directories, see their README file for info on how to build and use it. diff --git a/tools/fep/README b/tools/fep/README new file mode 100644 index 000000000..066396375 --- /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 000000000..a2c51a716 --- /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 000000000..c17eb5a03 --- /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 000000000..3ed9290ea --- /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 000000000..716dea7cb --- /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