diff --git a/doc/src/Tools.txt b/doc/src/Tools.txt index 6c41524d20..90f598890f 100644 --- a/doc/src/Tools.txt +++ b/doc/src/Tools.txt @@ -76,10 +76,10 @@ Post-processing tools :h3 "pymol_asphere"_#pymol, "python"_#pythontools, "reax"_#reax_tool, +"replica"_#replica, "smd"_#smd, "spin"_#spin, -"xmgrace"_#xmgrace, -"replica"_#replica :tb(c=6,ea=c,a=l) +"xmgrace"_#xmgrace :tb(c=6,ea=c,a=l) Miscellaneous tools :h3 @@ -486,6 +486,21 @@ README for more info on Pizza.py and how to use these scripts. :line +replica tool :h4,link(replica) + +The tools/replica directory contains the reorder_remd_traj python script which +can be used to reorder the replica trajectories (resulting from the use of the +temper command) according to temperature. This will produce discontinuous +trajectories with all frames at the same temperature in each trajectory. +Additional options can be used to calculate the canonical configurational +log-weight for each frame at each temperature using the pymbar package. See +the README.md file for further details. Try out the peptide example provided. + +This tool was written by (and is maintained by) Tanmoy Sanyal, +while at the Shell lab at UC Santa Barbara. (tanmoy dot 7989 at gmail.com) + +:line + reax tool :h4,link(reax_tool) The reax sub-directory contains stand-alone codes that can @@ -551,17 +566,3 @@ See the README file for details. These files were provided by Vikas Varshney (vv0210 at gmail.com) -:line - -replica tool :h4,link(replica) - -The tools/replica directory contains the reorder_remd_traj python script which -can be used to reorder the replica trajectories (resulting from the use of the -temper command) according to temperature. This will produce discontinuous -trajectories with all frames at the same temperature in each trajectory. -Additional options can be used to calculate the canonical configurational -log-weight for each frame at each temperature using the pymbar package. See -the README.md file for further details. Try out the peptide example provided. - -This tool was written by Tanmoy Sanyal, while at the Shell lab -at UC Santa Barbara. (tanmoy dot 7989 at gmail.com) diff --git a/doc/src/temper.txt b/doc/src/temper.txt index edd578fbc9..6a61dfa6dd 100644 --- a/doc/src/temper.txt +++ b/doc/src/temper.txt @@ -110,7 +110,13 @@ the information from the log.lammps file. E.g. you could produce one dump file with snapshots at 300K (from all replicas), another with snapshots at 310K, etc. Note that these new dump files will not contain "continuous trajectories" for individual atoms, because two -successive snapshots (in time) may be from different replicas. +successive snapshots (in time) may be from different replicas. The +reorder_remd_traj python script can do the reordering for you +(and additionally also calculated configurational log-weights of +trajectory snapshots in the canonical ensemble). The script can be found +in the tools/replica directory while instructions on how to use it is +available in doc/Tools (in brief) and as a README file in tools/replica +(in detail). The last argument {index} in the temper command is optional and is used when restarting a tempering run from a set of restart files (one diff --git a/tools/README b/tools/README index 54f8d86898..b20e82c53e 100644 --- a/tools/README +++ b/tools/README @@ -38,6 +38,7 @@ polybond Python tool for programmable polymer bonding 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 +replica tool to reorder LAMMPS replica trajectories according to temperature smd convert Smooth Mach Dynamics triangles to VTK spin perform a cubic polynomial interpolation of a GNEB MEP vim add-ons to VIM editor for editing LAMMPS input scripts diff --git a/tools/replica/example/in.peptide b/tools/replica/example/in.peptide index ada6941af1..5d321e34c3 100644 --- a/tools/replica/example/in.peptide +++ b/tools/replica/example/in.peptide @@ -14,7 +14,7 @@ dihedral_style charmm improper_style harmonic kspace_style pppm 0.0001 -read_data ../../../examples/peptide/data.peptide +read_data data.peptide neighbor 2.0 bin neigh_modify delay 5 diff --git a/tools/replica/reorder_remd_traj.py b/tools/replica/reorder_remd_traj.py index 5934a7fdd6..66e9489913 100644 --- a/tools/replica/reorder_remd_traj.py +++ b/tools/replica/reorder_remd_traj.py @@ -38,11 +38,11 @@ StringIO (or io if in Python 3.x) -import os, sys, numpy as np, argparse, time, pickle +import os, numpy as np, argparse, time, pickle from scipy.special import logsumexp from mpi4py import MPI -from tqdm import tqdm, trange +from tqdm import tqdm import gzip, bz2 try: # python-2 @@ -78,12 +78,10 @@ def _get_nearest_temp(temps, query_temp): """ if isinstance(temps, list): temps = np.array(temps) - idx = np.argmin(abs(temps - query_temp)) - out_temp = temps[idx] - return out_temp + return temps[np.argmin(np.abs(temps-query_temp))] -def readwrite(trajfn, mode = "rb"): +def readwrite(trajfn, mode): """ Helper function for input/output LAMMPS traj files. Trajectories may be plain text, .gz or .bz2 compressed. @@ -96,11 +94,14 @@ def readwrite(trajfn, mode = "rb"): """ if trajfn.endswith(".gz"): - return gzip.GzipFile(trajfn, mode) + of = gzip.open(trajfn, mode) + #return gzip.GzipFile(trajfn, mode) elif trajfn.endswith(".bz2"): - return bz2.BZ2File(trajfn, mode) + of = bz2.open(trajfn, mode) + #return bz2.BZ2File(trajfn, mode) else: - return file(trajfn, mode) + of = open(trajfn, mode) + return of def get_replica_frames(logfn, temps, nswap, writefreq): @@ -163,7 +164,7 @@ def get_byte_index(rep_inds, byteindfns, intrajfns): if os.path.isfile(byteindfns[n]): continue # extract bytes - fobj = readwrite(intrajfns[n]) + fobj = readwrite(intrajfns[n], "rb") byteinds = [ [0,0] ] # place file pointer at first line @@ -243,7 +244,7 @@ def write_reordered_traj(temp_inds, byte_inds, outtemps, temps, for n in temp_inds: # open string-buffer and file buf = IOBuffer() - of = readwrite(outtrajfns[n], mode = "wb") + of = readwrite(outtrajfns[n], "wb") # get frames abs_temp_ind = np.argmin( abs(temps - outtemps[n]) ) @@ -281,7 +282,7 @@ def write_reordered_traj(temp_inds, byte_inds, outtemps, temps, def get_canonical_logw(enefn, frametuple_dict, temps, nprod, writefreq, - kB = 0.001987): + kB): """ Gets configurational log-weights (logw) for each frame and at each temp. from the REMD simulation. ONLY WRITTEN FOR THE CANONICAL (NVT) ensemble. @@ -349,25 +350,25 @@ def get_canonical_logw(enefn, frametuple_dict, temps, nprod, writefreq, #3) get reduced energies (*ONLY FOR THE CANONICAL ENSEMBLE*) u_kln = np.zeros([ntemps, ntemps, nframes], float) for k in range(ntemps): - for l in range(ntemps): - u_kln[ k, l, 0:nframes_k[k] ] = beta_k[l] * u_kn[k, 0:nframes_k[k]] - + u_kln[k] = np.outer(beta_k, u_kn[k]) + # run pymbar and extract the free energies print("\nRunning pymbar...") mbar = pymbar.mbar.MBAR(u_kln, nframes_k, verbose = True) - f_k = mbar.f_k + f_k = mbar.f_k # (1 x k array) # calculate the log-weights print("\nExtracting log-weights...") log_nframes = np.log(nframes) logw = dict( (k, np.zeros([ntemps, nframes], float)) for k in range(ntemps) ) - for l in range(ntemps): - # get log-weights to reweight to this temp. - for k in range(ntemps): - for n in range(nframes): - num = -beta_k[k] * u_kn[k,n] - denom = f_k - beta_k[k] * u_kn[k,n] + # get log-weights to reweight to this temp. + for k in range(ntemps): + for n in range(nframes): + num = -beta_k[k] * u_kn[k,n] + denom = f_k - beta_k[k] * u_kn[k,n] + for l in range(ntemps): logw[l][k,n] = num - logsumexp(denom) - log_nframes + return logw @@ -516,8 +517,8 @@ if __name__ == "__main__": comm.barrier() # open all replica files for reading - infobjs = [readwrite(i) for i in intrajfns] - + infobjs = [readwrite(i, "rb") for i in intrajfns] + # open all byteindex files byte_inds = dict( (i, np.loadtxt(fn)) for i, fn in enumerate(byteindfns) )