From bbb0f5740efd2c4ea6734a29e1fe54271f3d2f92 Mon Sep 17 00:00:00 2001 From: "tanmoy.7989" Date: Fri, 6 Sep 2019 11:18:33 -0700 Subject: [PATCH 1/4] link to data.peptide was deleted by me by mistake. Now that it's re-added, I revoked (un-necessary) changes I made since to the in.peptide input script --- tools/replica/example/in.peptide | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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 From f41a1f83030f343aa1df647f269405e917adbc6c Mon Sep 17 00:00:00 2001 From: "tanmoy.7989" Date: Sun, 8 Sep 2019 10:43:22 -0700 Subject: [PATCH 2/4] vectorized in parts and made changes as suggested by evoyiatzis --- tools/replica/reorder_remd_traj.py | 47 +++++++++++++++--------------- 1 file changed, 24 insertions(+), 23 deletions(-) diff --git a/tools/replica/reorder_remd_traj.py b/tools/replica/reorder_remd_traj.py index c8c467ffe0..cef417cd0b 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. @@ -348,25 +349,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 @@ -515,7 +516,7 @@ 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) ) From 91a19719773e02f05de9b44c418a6d93cc2a143c Mon Sep 17 00:00:00 2001 From: "tanmoy.7989" Date: Mon, 9 Sep 2019 02:04:59 -0700 Subject: [PATCH 3/4] added a line to tools/README and fixed the alphabetical ordering in docs --- doc/src/Tools.txt | 29 +++++++++++++++-------------- tools/README | 1 + 2 files changed, 16 insertions(+), 14 deletions(-) diff --git a/doc/src/Tools.txt b/doc/src/Tools.txt index 6c41524d20..e1042f3a7e 100644 --- a/doc/src/Tools.txt +++ b/doc/src/Tools.txt @@ -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 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/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 From cde16580c086e98917042ea5648d3675b64e8d02 Mon Sep 17 00:00:00 2001 From: "tanmoy.7989" Date: Mon, 9 Sep 2019 14:15:05 -0700 Subject: [PATCH 4/4] fixed alphabetical ordering in Tools.txt and added a line highlighting the tool in temper.txt --- doc/src/Tools.txt | 8 ++++---- doc/src/temper.txt | 8 +++++++- 2 files changed, 11 insertions(+), 5 deletions(-) diff --git a/doc/src/Tools.txt b/doc/src/Tools.txt index e1042f3a7e..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 @@ -496,8 +496,8 @@ 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) +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 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