patch 17Mar17

Merge pull request #424 from timattox/dpd_restart_fix
Remove unneeded restart_peratom flags in USER-DPD package
2017-03-17 11:29:23 -06:00 · 2017-03-17 11:19:04 -06:00 · 2017-03-17 11:18:52 -06:00 · 2017-03-17 11:18:08 -06:00 · 2017-03-17 11:17:28 -06:00 · 2017-03-17 11:16:58 -06:00
56 changed files with 17241 additions and 152 deletions
--- a/doc/Makefile
+++ b/doc/Makefile
@ -22,7 +22,7 @@ endif
 SOURCES=$(wildcard src/*.txt)
 OBJECTS=$(SOURCES:src/%.txt=$(RSTDIR)/%.rst)

-.PHONY: help clean-all clean epub html pdf old venv spelling
+.PHONY: help clean-all clean epub html pdf old venv spelling anchor_check

 # ------------------------------------------

@ -36,6 +36,7 @@ help:
 	@echo "  clean      remove all intermediate RST files"
 	@echo "  clean-all  reset the entire build environment"
 	@echo "  txt2html   build txt2html tool"
+	@echo "  anchor_check  scan for duplicate anchor labels"

 # ------------------------------------------

@ -54,6 +55,9 @@ html: $(OBJECTS)
 		. $(VENV)/bin/activate ;\
 		cp -r src/* $(RSTDIR)/ ;\
 		sphinx-build -j 8 -b html -c utils/sphinx-config -d $(BUILDDIR)/doctrees $(RSTDIR) html ;\
+		echo "############################################" ;\
+		doc_anchor_check src/*.txt ;\
+		echo "############################################" ;\
 		deactivate ;\
 	)
 	-rm html/searchindex.js
@ -127,6 +131,13 @@ fetch:

 txt2html: utils/txt2html/txt2html.exe

+anchor_check : $(TXT2RST)
+	@(\
+		. $(VENV)/bin/activate ;\
+		doc_anchor_check src/*.txt ;\
+		deactivate ;\
+	)
+
 # ------------------------------------------

 utils/txt2html/txt2html.exe: utils/txt2html/txt2html.cpp
--- a/doc/src/Manual.txt
+++ b/doc/src/Manual.txt
@ -1,7 +1,7 @@
 <!-- HTML_ONLY -->
 <HEAD>
 <TITLE>LAMMPS Users Manual</TITLE>
-<META NAME="docnumber" CONTENT="10 Mar 2017 version">
+<META NAME="docnumber" CONTENT="17 Mar 2017 version">
 <META NAME="author" CONTENT="http://lammps.sandia.gov - Sandia National Laboratories">
 <META NAME="copyright" CONTENT="Copyright (2003) Sandia Corporation.  This software and manual is distributed under the GNU General Public License.">
 </HEAD>
@ -21,7 +21,7 @@
 <H1></H1>

 LAMMPS Documentation :c,h3
-10 Mar 2017 version :c,h4
+17 Mar 2017 version :c,h4

 Version info: :h4

--- a/doc/src/PDF/USER-CGDNA-overview.pdf
+++ b/doc/src/PDF/USER-CGDNA-overview.pdf
--- a/doc/src/Section_commands.txt
+++ b/doc/src/Section_commands.txt
@ -687,6 +687,7 @@ package"_Section_start.html#start_3.
 "eos/cv"_fix_eos_cv.html,
 "eos/table"_fix_eos_table.html,
 "eos/table/rx"_fix_eos_table_rx.html,
+"filter/corotate"_fix_filter_corotate.html,
 "flow/gauss"_fix_flow_gauss.html,
 "gle"_fix_gle.html,
 "grem"_fix_grem.html,
--- a/doc/src/Section_howto.txt
+++ b/doc/src/Section_howto.txt
@ -1032,6 +1032,10 @@ profile consistent with the applied shear strain rate.
 An alternative method for calculating viscosities is provided via the
 "fix viscosity"_fix_viscosity.html command.

+NEMD simulations can also be used to measure transport properties of a fluid
+through a pore or channel. Simulations of steady-state flow can be performed
+using the "fix flow/gauss"_fix_flow_gauss.html command.
+
 :line

 6.14 Finite-size spherical and aspherical particles :link(howto_14),h4
--- a/doc/src/fix_filter_corotate.txt
+++ b/doc/src/fix_filter_corotate.txt
@ -0,0 +1,87 @@
+"LAMMPS WWW Site"_lws - "LAMMPS Documentation"_ld - "LAMMPS Commands"_lc :c
+
+:link(lws,http://lammps.sandia.gov)
+:link(ld,Manual.html)
+:link(lc,Section_commands.html#comm)
+
+:line
+
+fix filter/corotate command :h3
+
+[Syntax:]
+
+fix ID group-ID filter/corotate keyword value ... :pre
+
+ID, group-ID are documented in "fix"_fix.html command :ulb,l
+one or more constraint/value pairs are appended :l
+constraint = {b} or {a} or {t} or {m} :l
+  {b} values = one or more bond types
+  {a} values = one or more angle types
+  {t} values = one or more atom types
+  {m} value = one or more mass values :pre
+:ule
+
+[Examples:]
+
+timestep 8
+run_style respa 3 2 8 bond 1 pair 2 kspace 3
+fix cor all filter/corotate m 1.0 :pre
+
+fix cor all filter/corotate b 4 19 a 3 5 2 :pre
+
+[Description:]
+
+This fix implements a corotational filter for a mollified impulse
+method. In biomolecular simulations, it allows the usage of larger
+timesteps for long-range electrostatic interactions.  For details, see
+"(Fath)"_#Fath2017.
+
+When using "run_style respa"_run_style.html for a biomolecular
+simulation with high-frequency covalent bonds, the outer time-step is
+restricted to below ~ 4fs due to resonance problems. This fix filters
+the outer stage of the respa and thus a larger (outer) time-step can
+be used. Since in large biomolecular simulations the computation of
+the long-range electrostatic contributions poses a major bottleneck,
+this can significantly accelerate the simulation.
+
+The filter computes a cluster decomposition of the molecular structure
+following the criteria indicated by the options a, b, t and m. This
+process is similar to the approach in "fix shake"_fix_shake.html,
+however, the clusters are not kept contrained. Instead, the position
+is slightly modified only for the computation of long-range forces. A
+good cluster decomposition constitutes in building clusters which
+contain the fastest covalent bonds inside clusters.
+
+If the clusters are chosen suitably, the "run_style
+respa"_run_style.html is stable for outer time-steps of at least 8fs.
+
+:line
+
+[Restart, fix_modify, output, run start/stop, minimize info:]
+
+No information about these fixes is written to "binary restart
+files"_restart.html.  None of the "fix_modify"_fix_modify.html options
+are relevant to these fixes.  No global or per-atom quantities are
+stored by these fixes for access by various "output
+commands"_Section_howto.html#howto_15.  No parameter of these fixes
+can be used with the {start/stop} keywords of the "run"_run.html
+command.  These fixes are not invoked during "energy
+minimization"_minimize.html.
+
+[Restrictions:]
+
+This fix is part of the USER-MISC package. It is only enabled if 
+LAMMPS was built with that package. See the "Making 
+LAMMPS"_Section_start.html#start_3 section for more info.
+
+Currently, it does not support "molecule templates"_molecule.html.
+
+[Related commands:]
+
+
+[Default:] none
+
+:line
+
+:link(Fath2017)
+[(Fath)] Fath, Hochbruck, Singh, J Comp Phys, 333, 180-198 (2017).
--- a/doc/src/fix_halt.txt
+++ b/doc/src/fix_halt.txt
@ -121,7 +121,7 @@ halt ID, so that the same condition is not immediately triggered in a
 subsequent run.

 The optional {message} keyword determines whether a message is printed
-to the screen and logfile when the half condition is triggered.  If
+to the screen and logfile when the halt condition is triggered.  If
 {message} is set to yes, a one line message with the values that
 triggered the halt is printed.  If {message} is set to no, no message
 is printed; the run simply exits.  The latter may be desirable for
--- a/doc/src/fix_qeq_reax.txt
+++ b/doc/src/fix_qeq_reax.txt
@ -36,7 +36,7 @@ uses charges on each atom.  The "fix qeq/comb"_fix_qeq_comb.html
 command should be used to perform charge equilibration with the "COMB
 potential"_pair_comb.html.  For more technical details about the
 charge equilibration performed by fix qeq/reax, see the
-"(Aktulga)"_#Aktulga paper.
+"(Aktulga)"_#qeq-Aktulga paper.

 The QEq method minimizes the electrostatic energy of the system by
 adjusting the partial charge on individual atoms based on interactions
@ -119,6 +119,6 @@ be used for periodic cell dimensions less than 10 angstroms.
 :link(Nakano)
 [(Nakano)] Nakano, Computer Physics Communications, 104, 59-69 (1997).

-:link(Aktulga)
-(Aktulga) Aktulga, Fogarty, Pandit, Grama, Parallel Computing, 38,
+:link(qeq-Aktulga)
+[(Aktulga)] Aktulga, Fogarty, Pandit, Grama, Parallel Computing, 38,
 245-259 (2012).
--- a/doc/src/fixes.txt
+++ b/doc/src/fixes.txt
@ -42,6 +42,7 @@ Fixes :h1
   fix_eos_table_rx
   fix_evaporate
   fix_external
+   fix_filter_corotate
   fix_flow_gauss
   fix_freeze
   fix_gcmc
--- a/doc/src/info.txt
+++ b/doc/src/info.txt
@ -12,7 +12,7 @@ info command :h3

 info args :pre

-args = one or more of the following keywords: {out}, {all}, {system}, {communication}, {computes}, {dumps}, {fixes}, {groups}, {regions}, {variables}, {styles}, {time}, or {configuration}
+args = one or more of the following keywords: {out}, {all}, {system}, {memory}, {communication}, {computes}, {dumps}, {fixes}, {groups}, {regions}, {variables}, {styles}, {time}, or {configuration}
     {out} values = {screen}, {log}, {append} filename, {overwrite} filename
     {styles} values = {all}, {angle}, {atom}, {bond}, {compute}, {command}, {dump}, {dihedral}, {fix}, {improper}, {integrate}, {kspace}, {minimize}, {pair}, {region} :ul

@ -40,6 +40,17 @@ to that file, which is either appended to or overwritten, respectively.

 The {all} flag activates printing all categories listed below.

+The {configuration} category prints some information about the
+LAMMPS version as well as architecture and OS it is run on.
+
+The {memory} category prints some information about the current
+memory allocation of MPI rank 0 (this the amount of dynamically
+allocated memory reported by LAMMPS classes). Where supported,
+also some OS specific information about the size of the reserved
+memory pool size (this is where malloc() and the new operator
+request memory from) and the maximum resident set size is reported
+(this is the maximum amount of physical memory occupied so far).
+
 The {system} category prints a general system overview listing.  This
 includes the unit style, atom style, number of atoms, bonds, angles,
 dihedrals, and impropers and the number of the respective types, box
@ -93,11 +104,6 @@ region :ul
 The {time} category prints the accumulated CPU and wall time for the
 process that writes output (usually MPI rank 0).

-The {configuration} command prints some information about the LAMMPS
-version and architecture and OS it is run on. Where supported, also
-information about the memory consumption provided by the OS is
-reported.
-
 [Restrictions:] none

 [Related commands:]
--- a/doc/src/lammps.book
+++ b/doc/src/lammps.book
@ -168,6 +168,7 @@ fix_eos_table.html
 fix_eos_table_rx.html
 fix_evaporate.html
 fix_external.html
+fix_filter_corotate.html
 fix_flow_gauss.html
 fix_freeze.html
 fix_gcmc.html
--- a/doc/src/package.txt
+++ b/doc/src/package.txt
@ -62,12 +62,13 @@ args = arguments specific to the style :l
      {no_affinity} values = none
  {kokkos} args = keyword value ...
    zero or more keyword/value pairs may be appended
-    keywords = {neigh} or {newton} or {binsize} or {comm} or {comm/exchange} or {comm/forward}
-      {neigh} value = {full} or {half} or {n2} or {full/cluster}
+    keywords = {neigh} or {neigh/qeq} or {newton} or {binsize} or {comm} or {comm/exchange} or {comm/forward}
+      {neigh} value = {full} or {half}
+        full = full neighbor list
+        half = half neighbor list built in thread-safe manner
+      {neigh/qeq} value = {full} or {half}
        full = full neighbor list
        half = half neighbor list built in thread-safe manner
-        n2 = non-binning neighbor list build, O(N^2) algorithm
-        full/cluster = full neighbor list with clustered groups of atoms
      {newton} = {off} or {on}
        off = set Newton pairwise and bonded flags off (default)
        on = set Newton pairwise and bonded flags on
@ -392,10 +393,7 @@ default value as listed below.

 The {neigh} keyword determines how neighbor lists are built.  A value
 of {half} uses a thread-safe variant of half-neighbor lists,
-the same as used by most pair styles in LAMMPS.  A value of
-{n2} uses an O(N^2) algorithm to build the neighbor list without
-binning, where N = # of atoms on a processor.  It is typically slower
-than the other methods, which use binning.
+the same as used by most pair styles in LAMMPS.

 A value of {full} uses a full neighbor lists and is the default.  This
 performs twice as much computation as the {half} option, however that
@ -403,15 +401,9 @@ is often a win because it is thread-safe and doesn't require atomic
 operations in the calculation of pair forces.  For that reason, {full}
 is the default setting.  However, when running in MPI-only mode with 1
 thread per MPI task, {half} neighbor lists will typically be faster,
-just as it is for non-accelerated pair styles.
-
-A value of {full/cluster} is an experimental neighbor style, where
-particles interact with all particles within a small cluster, if at
-least one of the clusters particles is within the neighbor cutoff
-range.  This potentially allows for better vectorization on
-architectures such as the Intel Phi.  If also reduces the size of the
-neighbor list by roughly a factor of the cluster size, thus reducing
-the total memory footprint considerably.
+just as it is for non-accelerated pair styles. Similarly, the {neigh/qeq}
+keyword determines how neighbor lists are built for "fix qeq/reax/kk"_fix_qeq_reax.html.
+If not explicitly set, the value of {neigh/qeq} will match {neigh}.

 The {newton} keyword sets the Newton flags for pairwise and bonded
 interactions to {off} or {on}, the same as the "newton"_newton.html
@ -582,9 +574,9 @@ is used.  If it is not used, you must invoke the package intel
 command in your input script or or via the "-pk intel" "command-line
 switch"_Section_start.html#start_7.

-For the KOKKOS package, the option defaults neigh = full, newton =
-off, binsize = 0.0, and comm = device.  These settings are made
-automatically by the required "-k on" "command-line
+For the KOKKOS package, the option defaults neigh = full, neigh/qeq
+ = full, newton = off, binsize = 0.0, and comm = device.  These settings
+ are made automatically by the required "-k on" "command-line
 switch"_Section_start.html#start_7.  You can change them bu using the
 package kokkos command in your input script or via the "-pk kokkos"
 "command-line switch"_Section_start.html#start_7.
--- a/doc/src/pair_polymorphic.txt
+++ b/doc/src/pair_polymorphic.txt
@ -52,7 +52,7 @@ to Stillinger-Weber potential ("SW"_#SW) if we set
 :c,image(Eqs/polymorphic4.jpg)

 The potential reduces to Tersoff types of potential
-("Tersoff"_#Tersoff or "Albe"_#Albe) if we set
+("Tersoff"_#Tersoff or "Albe"_#poly-Albe) if we set

 :c,image(Eqs/polymorphic5.jpg)
 :c,image(Eqs/polymorphic6.jpg)
@ -63,7 +63,7 @@ The potential reduces to Rockett-Tersoff ("Wang"_#Wang) type if we set
 :c,image(Eqs/polymorphic6.jpg)
 :c,image(Eqs/polymorphic8.jpg)

-The potential becomes embedded atom method ("Daw"_#Daw) if we set
+The potential becomes embedded atom method ("Daw"_#poly-Daw) if we set

 :c,image(Eqs/polymorphic9.jpg)

@ -218,12 +218,12 @@ F. P. Doty, J. Mater. Sci. Res., 4, 15 (2015).
 :link(Tersoff)
 [(Tersoff)] J. Tersoff, Phys. Rev. B, 39, 5566 (1989).

-:link(Albe)
+:link(poly-Albe)
 [(Albe)] K. Albe, K. Nordlund, J. Nord, and A. Kuronen, Phys. Rev. B,
 66, 035205 (2002).

 :link(Wang)
 [(Wang)] J. Wang, and A. Rockett, Phys. Rev. B, 43, 12571 (1991).

-:link(Daw)
+:link(poly-Daw)
 [(Daw)] M. S. Daw, and M. I. Baskes, Phys. Rev. B, 29, 6443 (1984).
--- a/doc/src/pair_tersoff_zbl.txt
+++ b/doc/src/pair_tersoff_zbl.txt
@ -23,9 +23,9 @@ pair_coeff * * SiC.tersoff.zbl Si C Si :pre
 [Description:]

 The {tersoff/zbl} style computes a 3-body Tersoff potential
-"(Tersoff_1)"_#Tersoff_1 with a close-separation pairwise modification
+"(Tersoff_1)"_#zbl-Tersoff_1 with a close-separation pairwise modification
 based on a Coulomb potential and the Ziegler-Biersack-Littmark
-universal screening function "(ZBL)"_#ZBL, giving the energy E of a
+universal screening function "(ZBL)"_#zbl-ZBL, giving the energy E of a
 system of atoms as

 :c,image(Eqs/pair_tersoff_zbl.jpg)
@ -146,16 +146,16 @@ be set to 0.0 if desired.
 Note that the twobody parameters in entries such as SiCC and CSiSi
 are often the same, due to the common use of symmetric mixing rules,
 but this is not always the case. For example, the beta and n parameters in
-Tersoff_2 "(Tersoff_2)"_#Tersoff_2 are not symmetric.
+Tersoff_2 "(Tersoff_2)"_#zbl-Tersoff_2 are not symmetric.

 We chose the above form so as to enable users to define all commonly
 used variants of the Tersoff portion of the potential.  In particular,
 our form reduces to the original Tersoff form when m = 3 and gamma =
-1, while it reduces to the form of "Albe et al."_#Albe when beta = 1
+1, while it reduces to the form of "Albe et al."_#zbl-Albe when beta = 1
 and m = 1.  Note that in the current Tersoff implementation in LAMMPS,
 m must be specified as either 3 or 1.  Tersoff used a slightly
 different but equivalent form for alloys, which we will refer to as
-Tersoff_2 potential "(Tersoff_2)"_#Tersoff_2.
+Tersoff_2 potential "(Tersoff_2)"_#zbl-Tersoff_2.

 LAMMPS parameter values for Tersoff_2 can be obtained as follows:
 gamma = omega_ijk, lambda3 = 0 and the value of
@ -253,16 +253,16 @@ units.

 :line

-:link(Tersoff_1)
+:link(zbl-Tersoff_1)
 [(Tersoff_1)] J. Tersoff, Phys Rev B, 37, 6991 (1988).

-:link(ZBL)
+:link(zbl-ZBL)
 [(ZBL)] J.F. Ziegler, J.P. Biersack, U. Littmark, 'Stopping and Ranges
 of Ions in Matter' Vol 1, 1985, Pergamon Press.

-:link(Albe)
+:link(zbl-Albe)
 [(Albe)] J. Nord, K. Albe, P. Erhart and K. Nordlund, J. Phys.:
 Condens. Matter, 15, 5649(2003).

-:link(Tersoff_2)
+:link(zbl-Tersoff_2)
 [(Tersoff_2)] J. Tersoff, Phys Rev B, 39, 5566 (1989); errata (PRB 41, 3248)
--- a/doc/utils/converters/lammpsdoc/doc_anchor_check.py
+++ b/doc/utils/converters/lammpsdoc/doc_anchor_check.py
@ -0,0 +1,60 @@
+#! /usr/bin/env python3
+# LAMMPS Documentation Utilities
+#
+# Scan for duplicate anchor labels in documentation files
+#
+# Copyright (C) 2017 Richard Berger
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program.  If not, see <http://www.gnu.org/licenses/>.
+import re
+import sys
+import argparse
+
+def main():
+    parser = argparse.ArgumentParser(description='scan for duplicate anchor labels in documentation files')
+    parser.add_argument('files',  metavar='file', nargs='+', help='one or more files to scan')
+    parsed_args = parser.parse_args()
+
+    anchor_pattern = re.compile(r'^:link\(([^,\)]*)\)')
+    anchors = {}
+
+    for filename in parsed_args.files:
+        with open(filename, 'rt') as f:
+            for line_number, line in enumerate(f):
+                m = anchor_pattern.match(line)
+                if m:
+                    label = m.group(1)
+                    if label in anchors:
+                        anchors[label].append((filename, line_number+1))
+                    else:
+                        anchors[label] = [(filename, line_number+1)]
+
+    count = 0
+
+    for label in sorted(anchors.keys()):
+        if len(anchors[label]) > 1:
+            print(label)
+            count += 1
+            for filename, line_number in anchors[label]:
+                print(" - %s:%d" % (filename, line_number))
+
+
+    if count > 0:
+        print("Found %d anchor label errors." % count)
+        sys.exit(1)
+    else:
+        print("No anchor label errors.")
+
+if __name__ == "__main__":
+    main()
--- a/doc/utils/converters/setup.py
+++ b/doc/utils/converters/setup.py
@ -12,6 +12,7 @@ setup(name='LAMMPS Documentation Utilities',
      tests_require=['nose'],
      entry_points = {
          "console_scripts": ['txt2html = lammpsdoc.txt2html:main',
-                              'txt2rst  = lammpsdoc.txt2rst:main']
+                              'txt2rst  = lammpsdoc.txt2rst:main',
+                              'doc_anchor_check = lammpsdoc.doc_anchor_check:main ']
      },
 )
--- a/examples/USER/cgdna/util/generate.py
+++ b/examples/USER/cgdna/util/generate.py
@ -0,0 +1,630 @@
+#!/usr/bin/env python
+"""
+/* ----------------------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   http://lammps.sandia.gov, Sandia National Laboratories
+   Steve Plimpton, sjplimp@sandia.gov
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+/* ----------------------------------------------------------------------
+   Contributing author: Oliver Henrich (EPCC, University of Edinburgh)
+------------------------------------------------------------------------- */
+"""
+
+
+"""
+Import basic modules
+"""
+import sys, os, timeit
+
+from timeit import default_timer as timer
+start_time = timer()
+"""
+Try to import numpy; if failed, import a local version mynumpy 
+which needs to be provided
+"""
+try:
+    import numpy as np
+except:
+    print >> sys.stderr, "numpy not found. Exiting."
+    sys.exit(1)
+
+"""
+Check that the required arguments (box offset and size in simulation units 
+and the sequence file were provided
+"""
+try:
+    box_offset = float(sys.argv[1])
+    box_length = float(sys.argv[2])
+    infile = sys.argv[3]
+except:
+    print >> sys.stderr, "Usage: %s <%s> <%s> <%s>" % (sys.argv[0], \
+	"box offset", "box length", "file with sequences")
+    sys.exit(1)
+box = np.array ([box_length, box_length, box_length])
+
+"""
+Try to open the file and fail gracefully if file cannot be opened
+"""
+try:
+    inp = open (infile, 'r')
+    inp.close()
+except:
+    print >> sys.stderr, "Could not open file '%s' for reading. \
+					      Aborting." % infile
+    sys.exit(2)
+
+# return parts of a string
+def partition(s, d):
+    if d in s:
+        sp = s.split(d, 1)
+        return sp[0], d, sp[1]
+    else:
+        return s, "", ""
+
+"""
+Define the model constants
+"""
+# set model constants
+PI = np.pi
+POS_BASE =  0.4
+POS_BACK = -0.4
+EXCL_RC1 =  0.711879214356
+EXCL_RC2 =  0.335388426126
+EXCL_RC3 =  0.52329943261
+
+"""
+Define auxillary variables for the construction of a helix
+"""
+# center of the double strand
+CM_CENTER_DS = POS_BASE + 0.2
+
+# ideal distance between base sites of two nucleotides 
+# which are to be base paired in a duplex
+BASE_BASE = 0.3897628551303122
+
+# cutoff distance for overlap check
+RC2 = 16
+
+# squares of the excluded volume distances for overlap check
+RC2_BACK = EXCL_RC1**2
+RC2_BASE = EXCL_RC2**2
+RC2_BACK_BASE = EXCL_RC3**2
+
+# enumeration to translate from letters to numbers and vice versa
+number_to_base = {1 : 'A', 2 : 'C', 3 : 'G', 4 : 'T'}
+base_to_number = {'A' : 1, 'a' : 1, 'C' : 2, 'c' : 2,
+                  'G' : 3, 'g' : 3, 'T' : 4, 't' : 4}
+
+# auxillary arrays
+positions = []
+a1s = []
+a3s = []
+quaternions = []
+
+newpositions = []
+newa1s = []
+newa3s = []
+
+basetype  = []
+strandnum = []
+
+bonds = []
+
+""" 
+Convert local body frame to quaternion DOF
+"""
+def exyz_to_quat (mya1, mya3):
+
+    mya2 = np.cross(mya3, mya1)
+    myquat = [1,0,0,0]
+
+    q0sq = 0.25 * (mya1[0] + mya2[1] + mya3[2] + 1.0)
+    q1sq = q0sq - 0.5 * (mya2[1] + mya3[2])
+    q2sq = q0sq - 0.5 * (mya1[0] + mya3[2])
+    q3sq = q0sq - 0.5 * (mya1[0] + mya2[1])
+
+    # some component must be greater than 1/4 since they sum to 1
+    # compute other components from it
+
+    if q0sq >= 0.25:
+	myquat[0] = np.sqrt(q0sq)
+	myquat[1] = (mya2[2] - mya3[1]) / (4.0*myquat[0])
+	myquat[2] = (mya3[0] - mya1[2]) / (4.0*myquat[0])
+	myquat[3] = (mya1[1] - mya2[0]) / (4.0*myquat[0])
+    elif q1sq >= 0.25:
+	myquat[1] = np.sqrt(q1sq)
+	myquat[0] = (mya2[2] - mya3[1]) / (4.0*myquat[1])
+	myquat[2] = (mya2[0] + mya1[1]) / (4.0*myquat[1])
+	myquat[3] = (mya1[2] + mya3[0]) / (4.0*myquat[1])
+    elif q2sq >= 0.25:
+	myquat[2] = np.sqrt(q2sq)
+	myquat[0] = (mya3[0] - mya1[2]) / (4.0*myquat[2])
+	myquat[1] = (mya2[0] + mya1[1]) / (4.0*myquat[2])
+	myquat[3] = (mya3[1] + mya2[2]) / (4.0*myquat[2])
+    elif q3sq >= 0.25:
+	myquat[3] = np.sqrt(q3sq)
+	myquat[0] = (mya1[1] - mya2[0]) / (4.0*myquat[3])
+	myquat[1] = (mya3[0] + mya1[2]) / (4.0*myquat[3])
+	myquat[2] = (mya3[1] + mya2[2]) / (4.0*myquat[3])
+
+    norm = 1.0/np.sqrt(myquat[0]*myquat[0] + myquat[1]*myquat[1] + \
+			  myquat[2]*myquat[2] + myquat[3]*myquat[3])
+    myquat[0] *= norm
+    myquat[1] *= norm
+    myquat[2] *= norm
+    myquat[3] *= norm
+
+    return np.array([myquat[0],myquat[1],myquat[2],myquat[3]])
+
+"""
+Adds a strand to the system by appending it to the array of previous strands
+"""
+def add_strands (mynewpositions, mynewa1s, mynewa3s):
+    overlap = False
+	
+    # This is a simple check for each of the particles where for previously 
+    # placed particles i we check whether it overlaps with any of the 
+    # newly created particles j
+
+    print >> sys.stdout, "## Checking for overlaps"
+
+    for i in xrange(len(positions)):
+
+	p = positions[i]
+	pa1 = a1s[i]
+
+	for j in xrange (len(mynewpositions)):
+
+	    q = mynewpositions[j]
+	    qa1 = mynewa1s[j]
+
+	    # skip particles that are anyway too far away
+	    dr = p - q
+	    dr -= box * np.rint (dr / box)
+	    if np.dot(dr, dr) > RC2:
+		continue
+
+	    # base site and backbone site of the two particles
+            p_pos_back = p + pa1 * POS_BACK
+            p_pos_base = p + pa1 * POS_BASE
+            q_pos_back = q + qa1 * POS_BACK
+            q_pos_base = q + qa1 * POS_BASE
+
+	    # check for no overlap between the two backbone sites
+            dr = p_pos_back - q_pos_back
+            dr -= box * np.rint (dr / box)
+            if np.dot(dr, dr) < RC2_BACK:
+                overlap = True
+
+	    # check for no overlap between the two base sites
+            dr = p_pos_base -  q_pos_base
+            dr -= box * np.rint (dr / box)
+            if np.dot(dr, dr) < RC2_BASE:
+                overlap = True
+
+	    # check for no overlap between backbone site of particle p 
+	    # with base site of particle q
+            dr = p_pos_back - q_pos_base
+            dr -= box * np.rint (dr / box)
+            if np.dot(dr, dr) < RC2_BACK_BASE:
+                overlap = True
+
+	    # check for no overlap between base site of particle p and 
+	    # backbone site of particle q
+            dr = p_pos_base - q_pos_back
+            dr -= box * np.rint (dr / box)
+            if np.dot(dr, dr) < RC2_BACK_BASE:
+                overlap = True
+
+	    # exit if there is an overlap
+            if overlap:
+                return False
+
+    # append to the existing list if no overlap is found
+    if not overlap:
+
+        for p in mynewpositions:
+            positions.append(p)
+        for p in mynewa1s:
+            a1s.append (p)
+        for p in mynewa3s:
+            a3s.append (p)
+	# calculate quaternion from local body frame and append
+	for ia in xrange(len(mynewpositions)):
+	    mynewquaternions = exyz_to_quat(mynewa1s[ia],mynewa3s[ia])
+	    quaternions.append(mynewquaternions)
+
+    return True
+
+
+"""
+Returns the rotation matrix defined by an axis and angle
+"""
+def get_rotation_matrix(axis, anglest):
+    # The argument anglest can be either an angle in radiants
+    # (accepted types are float, int or np.float64 or np.float64)
+    # or a tuple [angle, units] where angle is a number and
+    # units is a string. It tells the routine whether to use degrees,
+    # radiants (the default) or base pairs turns.
+    if not isinstance (anglest, (np.float64, np.float32, float, int)):
+        if len(anglest) > 1:
+            if anglest[1] in ["degrees", "deg", "o"]:
+                #angle = np.deg2rad (anglest[0])
+                angle = (np.pi / 180.) * (anglest[0])
+            elif anglest[1] in ["bp"]:
+                angle = int(anglest[0]) * (np.pi / 180.) * (35.9)
+            else:
+                angle = float(anglest[0])
+        else:
+            angle = float(anglest[0])
+    else:
+        angle = float(anglest) # in degrees (?)
+
+    axis = np.array(axis)
+    axis /= np.sqrt(np.dot(axis, axis))
+
+    ct = np.cos(angle)
+    st = np.sin(angle)
+    olc = 1. - ct
+    x, y, z = axis
+
+    return np.array([[olc*x*x+ct, olc*x*y-st*z, olc*x*z+st*y],
+                    [olc*x*y+st*z, olc*y*y+ct, olc*y*z-st*x],
+                    [olc*x*z-st*y, olc*y*z+st*x, olc*z*z+ct]])
+
+"""
+Generates the position and orientation vectors of a 
+(single or double) strand from a sequence string
+"""
+def generate_strand(bp, sequence=None, start_pos=np.array([0, 0, 0]), \
+	  dir=np.array([0, 0, 1]), perp=False, double=True, rot=0.):
+    # generate empty arrays
+    mynewpositions, mynewa1s, mynewa3s = [], [], []
+
+    # cast the provided start_pos array into a numpy array
+    start_pos = np.array(start_pos, dtype=float)
+
+    # overall direction of the helix
+    dir = np.array(dir, dtype=float)
+    if sequence == None:
+	sequence = np.random.randint(1, 5, bp)
+
+    # the elseif here is most likely redundant 
+    elif len(sequence) != bp:
+	n = bp - len(sequence)
+	sequence += np.random.randint(1, 5, n)
+	print >> sys.stderr, "sequence is too short, adding %d random bases" % n
+
+    # normalize direction
+    dir_norm = np.sqrt(np.dot(dir,dir))
+    if dir_norm < 1e-10:
+	print >> sys.stderr, "direction must be a valid vector, \
+			      defaulting to (0, 0, 1)"
+	dir = np.array([0, 0, 1])
+    else: dir /= dir_norm
+
+    # find a vector orthogonal to dir to act as helix direction,
+    # if not provided switch off random orientation
+    if perp is None or perp is False:
+	v1 = np.random.random_sample(3)
+	v1 -= dir * (np.dot(dir, v1))
+	v1 /= np.sqrt(sum(v1*v1))
+    else:
+	v1 = perp;
+
+    # generate rotational matrix representing the overall rotation of the helix
+    R0 = get_rotation_matrix(dir, rot)
+	    
+    # rotation matrix corresponding to one step along the helix
+    R = get_rotation_matrix(dir, [1, "bp"])
+
+    # set the vector a1 (backbone to base) to v1 
+    a1 = v1
+    
+    # apply the global rotation to a1 
+    a1 = np.dot(R0, a1)
+    
+    # set the position of the fist backbone site to start_pos
+    rb = np.array(start_pos)
+	    
+    # set a3 to the direction of the helix
+    a3 = dir
+    for i in range(bp):
+    # work out the position of the centre of mass of the nucleotide
+	rcdm = rb - CM_CENTER_DS * a1
+	
+	# append to newpositions
+	mynewpositions.append(rcdm)
+	mynewa1s.append(a1)
+	mynewa3s.append(a3)
+	
+	# if we are not at the end of the helix, we work out a1 and rb for the 
+	# next nucleotide along the helix
+	if i != bp - 1:
+	    a1 = np.dot(R, a1)
+	    rb += a3 * BASE_BASE
+
+    # if we are working on a double strand, we do a cycle similar 
+    # to the previous one but backwards
+    if double == True:
+	a1 = -a1
+	a3 = -dir
+	R = R.transpose()
+	for i in range(bp):
+	    rcdm = rb - CM_CENTER_DS * a1
+	    mynewpositions.append (rcdm)
+	    mynewa1s.append (a1)
+	    mynewa3s.append (a3)
+	    a1 = np.dot(R, a1)
+	    rb += a3 * BASE_BASE
+
+    assert (len (mynewpositions) > 0)
+
+    return [mynewpositions, mynewa1s, mynewa3s]
+
+
+"""
+Main function for this script.
+Reads a text file with the following format:
+- Each line contains the sequence for a single strand (A,C,G,T)
+- Lines beginning with the keyword 'DOUBLE' produce double-stranded DNA
+
+Ex: Two ssDNA (single stranded DNA)
+ATATATA
+GCGCGCG
+
+Ex: Two strands, one double stranded, the other single stranded.
+DOUBLE AGGGCT
+CCTGTA
+
+"""
+
+def read_strands(filename):
+    try:
+        infile = open (filename)
+    except:
+        print >> sys.stderr, "Could not open file '%s'. Aborting." % filename
+        sys.exit(2)
+
+    # This block works out the number of nucleotides and strands by reading 
+    # the number of non-empty lines in the input file and the number of letters,
+    # taking the possible DOUBLE keyword into account.
+    nstrands, nnucl, nbonds = 0, 0, 0
+    lines = infile.readlines()
+    for line in lines:
+        line = line.upper().strip()
+        if len(line) == 0:
+            continue
+        if line[:6] == 'DOUBLE':
+            line = line.split()[1]
+            length = len(line)
+            print >> sys.stdout, "## Found duplex of %i base pairs" % length
+            nnucl += 2*length
+            nstrands += 2
+	    nbonds += (2*length-2)
+        else:
+            line = line.split()[0]
+            length = len(line)
+            print >> sys.stdout, \
+		    "## Found single strand of %i bases" % length
+            nnucl += length
+            nstrands += 1
+	    nbonds += length-1
+    # rewind the sequence input file
+    infile.seek(0)
+
+    print >> sys.stdout, "## nstrands, nnucl = ", nstrands, nnucl
+
+    # generate the data file in LAMMPS format
+    try:
+        out = open ("data.oxdna", "w")
+    except:
+        print >> sys.stderr, "Could not open data file for writing. Aborting."
+        sys.exit(2)
+	
+    lines = infile.readlines()
+    nlines = len(lines)
+    i = 1
+    myns = 0
+    noffset = 1
+
+    for line in lines:
+        line = line.upper().strip()
+
+        # skip empty lines
+        if len(line) == 0: 
+	    i += 1
+	    continue
+
+	# block for duplexes: last argument of the generate function 
+	# is set to 'True'
+        if line[:6] == 'DOUBLE':
+            line = line.split()[1]
+            length = len(line)
+            seq = [(base_to_number[x]) for x in line]
+
+	    myns += 1
+	    for b in xrange(length):
+		basetype.append(seq[b])
+		strandnum.append(myns)
+
+	    for b in xrange(length-1):
+		bondpair = [noffset + b, noffset + b + 1]
+		bonds.append(bondpair)
+	    noffset += length
+
+	    # create the sequence of the second strand as made of 
+	    # complementary bases
+	    seq2 = [5-s for s in seq]
+	    seq2.reverse()
+
+	    myns += 1
+	    for b in xrange(length):
+		basetype.append(seq2[b])
+		strandnum.append(myns)
+
+	    for b in xrange(length-1):
+		bondpair = [noffset + b, noffset + b + 1]
+		bonds.append(bondpair)
+	    noffset += length
+ 
+            print >> sys.stdout, "## Created duplex of %i bases" % (2*length)
+
+	    # generate random position of the first nucleotide
+            cdm = box_offset + np.random.random_sample(3) * box
+
+            # generate the random direction of the helix 
+            axis = np.random.random_sample(3)
+            axis /= np.sqrt(np.dot(axis, axis))
+
+            # use the generate function defined above to create 
+	    # the position and orientation vector of the strand 
+            newpositions, newa1s, newa3s = generate_strand(len(line), \
+		    sequence=seq, dir=axis, start_pos=cdm, double=True)
+
+            # generate a new position for the strand until it does not overlap
+	    # with anything already present
+	    start = timer()
+            while not add_strands(newpositions, newa1s, newa3s):
+                cdm = box_offset + np.random.random_sample(3) * box
+                axis = np.random.random_sample(3)
+                axis /= np.sqrt(np.dot(axis, axis))
+                newpositions, newa1s, newa3s = generate_strand(len(line), \
+		      sequence=seq, dir=axis, start_pos=cdm, double=True)
+                print >> sys.stdout, "## Trying %i" % i
+	    end = timer()
+            print >> sys.stdout, "## Added duplex of %i bases (line %i/%i) in %.2fs, now at %i/%i" % \
+				      (2*length, i, nlines, end-start, len(positions), nnucl)
+
+	# block for single strands: last argument of the generate function 
+	# is set to 'False'
+        else:
+            length = len(line)
+            seq = [(base_to_number[x]) for x in line]
+
+	    myns += 1
+	    for b in xrange(length):
+		basetype.append(seq[b])
+		strandnum.append(myns)
+
+	    for b in xrange(length-1):
+		bondpair = [noffset + b, noffset + b + 1]
+		bonds.append(bondpair)
+	    noffset += length
+
+	    # generate random position of the first nucleotide
+            cdm = box_offset + np.random.random_sample(3) * box
+
+            # generate the random direction of the helix 
+            axis = np.random.random_sample(3)
+            axis /= np.sqrt(np.dot(axis, axis))
+
+            print >> sys.stdout, \
+		      "## Created single strand of %i bases" % length
+
+            newpositions, newa1s, newa3s = generate_strand(length, \
+		      sequence=seq, dir=axis, start_pos=cdm, double=False)
+	    start = timer()
+            while not add_strands(newpositions, newa1s, newa3s):
+                cdm = box_offset + np.random.random_sample(3) * box
+                axis = np.random.random_sample(3)
+		axis /= np.sqrt(np.dot(axis, axis))
+                newpositions, newa1s, newa3s = generate_strand(length, \
+			  sequence=seq, dir=axis, start_pos=cdm, double=False)
+                print >> sys.stdout, "## Trying  %i" % (i)
+	    end = timer()
+            print >> sys.stdout, "## Added single strand of %i bases (line %i/%i) in %.2fs, now at %i/%i" % \
+				      (length, i, nlines, end-start,len(positions), nnucl)
+
+        i += 1
+
+    # sanity check
+    if not len(positions) == nnucl:
+        print len(positions), nnucl
+        raise AssertionError
+
+    out.write('# LAMMPS data file\n')
+    out.write('%d atoms\n' % nnucl)
+    out.write('%d ellipsoids\n' % nnucl)
+    out.write('%d bonds\n' % nbonds)
+    out.write('\n')
+    out.write('4 atom types\n')
+    out.write('1 bond types\n')
+    out.write('\n')
+    out.write('# System size\n')
+    out.write('%f %f xlo xhi\n' % (box_offset,box_offset+box_length))
+    out.write('%f %f ylo yhi\n' % (box_offset,box_offset+box_length))
+    out.write('%f %f zlo zhi\n' % (box_offset,box_offset+box_length))
+
+    out.write('\n')
+    out.write('Masses\n')
+    out.write('\n')
+    out.write('1 3.1575\n')
+    out.write('2 3.1575\n')
+    out.write('3 3.1575\n')
+    out.write('4 3.1575\n')
+
+    # for each nucleotide print a line under the headers
+    # Atoms, Velocities, Ellipsoids and Bonds
+    out.write('\n')
+    out.write(\
+      '# Atom-ID, type, position, molecule-ID, ellipsoid flag, density\n')
+    out.write('Atoms\n')
+    out.write('\n')
+
+    for i in xrange(nnucl):
+	out.write('%d %d %22.15le %22.15le %22.15le %d 1 1\n' \
+		  % (i+1, basetype[i], \
+		     positions[i][0], positions[i][1], positions[i][2], \
+		     strandnum[i]))
+
+    out.write('\n')
+    out.write('# Atom-ID, translational, rotational velocity\n')
+    out.write('Velocities\n')
+    out.write('\n')
+
+    for i in xrange(nnucl):
+	out.write("%d %22.15le %22.15le %22.15le %22.15le %22.15le %22.15le\n" \
+		  % (i+1,0.0,0.0,0.0,0.0,0.0,0.0))
+
+    out.write('\n')
+    out.write('# Atom-ID, shape, quaternion\n')
+    out.write('Ellipsoids\n')
+    out.write('\n')
+
+    for i in xrange(nnucl):
+	out.write(\
+    "%d %22.15le %22.15le %22.15le %22.15le %22.15le %22.15le %22.15le\n"  \
+      % (i+1,1.1739845031423408,1.1739845031423408,1.1739845031423408, \
+	quaternions[i][0],quaternions[i][1], quaternions[i][2],quaternions[i][3]))
+ 
+    out.write('\n')
+    out.write('# Bond topology\n')
+    out.write('Bonds\n')
+    out.write('\n')
+
+    for i in xrange(nbonds):
+	out.write("%d  %d  %d  %d\n" % (i+1,1,bonds[i][0],bonds[i][1]))
+
+    out.close()
+
+    print >> sys.stdout, "## Wrote data to 'data.oxdna'"
+    print >> sys.stdout, "## DONE"
+
+# call the above main() function, which executes the program
+read_strands (infile)
+
+end_time=timer()
+runtime = end_time-start_time
+hours = runtime/3600
+minutes = (runtime-np.rint(hours)*3600)/60
+seconds = (runtime-np.rint(hours)*3600-np.rint(minutes)*60)%60
+print >> sys.stdout, "## Total runtime %ih:%im:%.2fs" % (hours,minutes,seconds)
--- a/examples/USER/cgdna/util/generate_simple.py
+++ b/examples/USER/cgdna/util/generate_simple.py
@ -29,7 +29,7 @@ def single():

  strandstart=len(nucleotide)+1

-  for letter in strand[2]:
+  for letter in strand[1]:
    temp=[]

    temp.append(nt2num[letter])
@ -58,7 +58,7 @@ def single_helix():
  strand = inp[1].split(':')

  com_start=strand[0].split(',')
-  twist=float(strand[1])
+  twist=0.6

  posx = float(com_start[0])
  posy = float(com_start[1])
@ -79,7 +79,7 @@ def single_helix():
  qrot2=0
  qrot3=math.sin(0.5*twist)

-  for letter in strand[2]:
+  for letter in strand[1]:
    temp=[]

    temp.append(nt2num[letter])
@ -114,7 +114,7 @@ def duplex():
  strand = inp[1].split(':')

  com_start=strand[0].split(',')
-  twist=float(strand[1])
+  twist=0.6

  compstrand=[]
  comptopo=[]
@ -145,7 +145,7 @@ def duplex():
  qrot2=0
  qrot3=math.sin(0.5*twist)

-  for letter in strand[2]:
+  for letter in strand[1]:
    temp1=[]
    temp2=[]

@ -189,7 +189,7 @@ def duplex():

    if (len(nucleotide)+1 > strandstart):
      topology.append([1,len(nucleotide),len(nucleotide)+1])
-      comptopo.append([1,len(nucleotide)+len(strand[2]),len(nucleotide)+len(strand[2])+1])
+      comptopo.append([1,len(nucleotide)+len(strand[1]),len(nucleotide)+len(strand[1])+1])

    nucleotide.append(temp1)
    compstrand.append(temp2)
@ -208,7 +208,7 @@ def duplex_array():
  strand = inp[1].split(':')
  number=strand[0].split(',')
  posz1_0 = float(strand[1])
-  twist=float(strand[2])
+  twist=0.6

  nx = int(number[0])
  ny = int(number[1])
@ -249,7 +249,7 @@ def duplex_array():
      qrot2=0
      qrot3=math.sin(0.5*twist)

-      for letter in strand[3]:
+      for letter in strand[2]:
 	temp1=[]
 	temp2=[]

@ -293,7 +293,7 @@ def duplex_array():

 	if (len(nucleotide)+1 > strandstart):
 	  topology.append([1,len(nucleotide),len(nucleotide)+1])
-	  comptopo.append([1,len(nucleotide)+len(strand[3]),len(nucleotide)+len(strand[3])+1])
+	  comptopo.append([1,len(nucleotide)+len(strand[2]),len(nucleotide)+len(strand[2])+1])

 	nucleotide.append(temp1)
 	compstrand.append(temp2)
--- a/examples/USER/cgdna/util/sequence.txt
+++ b/examples/USER/cgdna/util/sequence.txt
@ -1,4 +1,3 @@
-single 0,0,0:0.6:AAAAA
-single_helix 0,0,0:0.6:AAAAA
-duplex 0,0,0:0.6:AAAAA
-duplex_array 10,10:-112.0:0.6:AAAAA
+DOUBLE ACGTA
+
+ACGTA
--- a/examples/USER/cgdna/util/sequence_simple.txt
+++ b/examples/USER/cgdna/util/sequence_simple.txt
@ -0,0 +1,4 @@
+single 0,0,0:AAAAA
+single_helix 0,0,0:AAAAA
+duplex 0,0,0:AAAAA
+duplex_array 10,10:-112.0:AAAAA
--- a/examples/USER/misc/filter-corotate/data.bpti
+++ b/examples/USER/misc/filter-corotate/data.bpti
--- a/examples/USER/misc/filter-corotate/data.peptide
+++ b/examples/USER/misc/filter-corotate/data.peptide
--- a/examples/USER/misc/filter-corotate/in.bpti
+++ b/examples/USER/misc/filter-corotate/in.bpti
@ -0,0 +1,43 @@
+
+units           real
+
+atom_style      full
+bond_style      harmonic
+angle_style     charmm
+dihedral_style  charmm
+improper_style  harmonic
+
+pair_style      lj/charmm/coul/long 8 10
+pair_modify     mix arithmetic
+kspace_style    pppm 1e-4
+
+read_data       data.bpti
+
+special_bonds   charmm
+neigh_modify    delay 2 every 1
+
+
+# ------------- MINIMIZE ----------
+
+minimize 	1e-4 1e-6 1000 10000
+reset_timestep  0
+
+# ------------- RUN ---------------
+
+thermo          100
+thermo_style    multi
+timestep	8
+
+run_style respa 3 2 8 bond 1 pair 2 kspace 3
+
+velocity        all create 200.0 12345678 dist uniform
+#dump            dump1 all atom 100 4pti.dump
+
+fix             1 all nvt temp 200 300 25
+fix             cor all filter/corotate m 1.0
+
+run             1000
+
+unfix           cor
+unfix           1
+
--- a/examples/USER/misc/filter-corotate/in.peptide
+++ b/examples/USER/misc/filter-corotate/in.peptide
@ -0,0 +1,32 @@
+# Solvated 5-mer peptide, run for 8ps in NVT
+
+units           real
+atom_style      full
+
+pair_style      lj/charmm/coul/long 8.0 10.0 10.0
+bond_style      harmonic
+angle_style     charmm
+dihedral_style  charmm
+improper_style  harmonic
+kspace_style    pppm 0.0001
+
+read_data       data.peptide
+
+neighbor        2.0 bin
+neigh_modify    delay 5
+
+thermo          50
+#dump            dump1 all atom 100 peptide.dump
+
+timestep        8
+
+run_style respa 3 2 8 bond 1 pair 2 kspace 3
+
+fix             1 all nvt temp 250.0 250.0 100.0 tchain 1
+fix             cor all filter/corotate m 1.0
+run             1000
+unfix           cor
+unfix           1
+
+
+
--- a/examples/USER/misc/filter-corotate/log.10Mar2017.bpti.g++.1
+++ b/examples/USER/misc/filter-corotate/log.10Mar2017.bpti.g++.1
@ -0,0 +1,240 @@
+LAMMPS (10 Mar 2017)
+  using 1 OpenMP thread(s) per MPI task
+
+units           real
+
+atom_style      full
+bond_style      harmonic
+angle_style     charmm
+dihedral_style  charmm
+improper_style  harmonic
+
+pair_style      lj/charmm/coul/long 8 10
+pair_modify     mix arithmetic
+kspace_style    pppm 1e-4
+
+read_data       data.bpti
+  orthogonal box = (-10 -10 -30) to (50 50 30)
+  1 by 1 by 1 MPI processor grid
+  reading atoms ...
+  892 atoms
+  scanning bonds ...
+  4 = max bonds/atom
+  scanning angles ...
+  6 = max angles/atom
+  scanning dihedrals ...
+  18 = max dihedrals/atom
+  scanning impropers ...
+  2 = max impropers/atom
+  reading bonds ...
+  906 bonds
+  reading angles ...
+  1626 angles
+  reading dihedrals ...
+  2501 dihedrals
+  reading impropers ...
+  137 impropers
+  4 = max # of 1-2 neighbors
+  9 = max # of 1-3 neighbors
+  19 = max # of 1-4 neighbors
+  21 = max # of special neighbors
+
+special_bonds   charmm
+neigh_modify    delay 2 every 1
+
+
+# ------------- MINIMIZE ----------
+
+minimize 	1e-4 1e-6 1000 10000
+WARNING: Resetting reneighboring criteria during minimization (../min.cpp:168)
+PPPM initialization ...
+WARNING: System is not charge neutral, net charge = 6 (../kspace.cpp:302)
+WARNING: Using 12-bit tables for long-range coulomb (../kspace.cpp:321)
+  G vector (1/distance) = 0.203272
+  grid = 16 16 16
+  stencil order = 5
+  estimated absolute RMS force accuracy = 0.0316399
+  estimated relative force accuracy = 9.52826e-05
+  using double precision FFTs
+  3d grid and FFT values/proc = 9261 4096
+Neighbor list info ...
+  update every 1 steps, delay 0 steps, check yes
+  max neighbors/atom: 2000, page size: 100000
+  master list distance cutoff = 12
+  ghost atom cutoff = 12
+  binsize = 6, bins = 10 10 10
+  1 neighbor lists, perpetual/occasional/extra = 1 0 0
+  (1) pair lj/charmm/coul/long, perpetual
+      attributes: half, newton on
+      pair build: half/bin/newton
+      stencil: half/bin/3d/newton
+      bin: standard
+Per MPI rank memory usage (min/avg/max) = 17.8596/1/0 Mbytes
+Step Temp E_pair E_mol TotEng Press 
+       0            0   -3075.6498    943.91164   -2131.7381   -380.67776 
+     241            0    -4503.313    749.58662   -3753.7264   -29.045104 
+Loop time of 3.35722 on 1 procs for 241 steps with 892 atoms
+
+99.7% CPU use with 1 MPI tasks x 1 OpenMP threads
+
+Minimization stats:
+  Stopping criterion = energy tolerance
+  Energy initial, next-to-last, final = 
+        -2131.73812515     -3753.43984087     -3753.72636847
+  Force two-norm initial, final = 1086.21 26.3688
+  Force max component initial, final = 310.811 3.92748
+  Final line search alpha, max atom move = 0.00596649 0.0234333
+  Iterations, force evaluations = 241 463
+
+MPI task timing breakdown:
+Section |  min time  |  avg time  |  max time  |%varavg| %total
+---------------------------------------------------------------
+Pair    | 2.5003     | 2.5003     | 2.5003     |   0.0 | 74.48
+Bond    | 0.24287    | 0.24287    | 0.24287    |   0.0 |  7.23
+Kspace  | 0.53428    | 0.53428    | 0.53428    |   0.0 | 15.91
+Neigh   | 0.069765   | 0.069765   | 0.069765   |   0.0 |  2.08
+Comm    | 0.00065374 | 0.00065374 | 0.00065374 |   0.0 |  0.02
+Output  | 0          | 0          | 0          |   0.0 |  0.00
+Modify  | 0          | 0          | 0          |   0.0 |  0.00
+Other   |            | 0.009358   |            |       |  0.28
+
+Nlocal:    892 ave 892 max 892 min
+Histogram: 1 0 0 0 0 0 0 0 0 0
+Nghost:    31 ave 31 max 31 min
+Histogram: 1 0 0 0 0 0 0 0 0 0
+Neighs:    148891 ave 148891 max 148891 min
+Histogram: 1 0 0 0 0 0 0 0 0 0
+
+Total # of neighbors = 148891
+Ave neighs/atom = 166.918
+Ave special neighs/atom = 10.9395
+Neighbor list builds = 15
+Dangerous builds = 0
+reset_timestep  0
+
+# ------------- RUN ---------------
+
+thermo          100
+thermo_style    multi
+timestep	8
+
+run_style respa 3 2 8 bond 1 pair 2 kspace 3
+Respa levels:
+  1 = bond angle dihedral improper
+  2 = pair
+  3 = kspace
+
+velocity        all create 200.0 12345678 dist uniform
+#dump            dump1 all atom 100 4pti.dump
+
+fix             1 all nvt temp 200 300 25
+fix             cor all filter/corotate m 1.0
+  163 = # of size 2 clusters
+  0 = # of size 3 clusters
+  25 = # of size 4 clusters
+  0 = # of size 5 clusters
+  100 = # of frozen angles
+
+run             1000
+PPPM initialization ...
+WARNING: Using 12-bit tables for long-range coulomb (../kspace.cpp:321)
+  G vector (1/distance) = 0.203272
+  grid = 16 16 16
+  stencil order = 5
+  estimated absolute RMS force accuracy = 0.0316399
+  estimated relative force accuracy = 9.52826e-05
+  using double precision FFTs
+  3d grid and FFT values/proc = 9261 4096
+Per MPI rank memory usage (min/avg/max) = 19.5425/1/0 Mbytes
+---------------- Step        0 ----- CPU =      0.0000 (sec) ----------------
+TotEng   =     -3220.3378 KinEng   =       531.1804 Temp     =       200.0000 
+PotEng   =     -3751.5181 E_bond   =        42.2810 E_angle  =       345.2592 
+E_dihed  =       337.8361 E_impro  =        24.2103 E_vdwl   =      -288.5339 
+E_coul   =      -886.3622 E_long   =     -3326.2088 Press    =        83.2283 
+---------------- Step      100 ----- CPU =      3.9414 (sec) ----------------
+TotEng   =     -2718.8970 KinEng   =       538.6206 Temp     =       202.8014 
+PotEng   =     -3257.5176 E_bond   =       203.3367 E_angle  =       566.5317 
+E_dihed  =       397.6202 E_impro  =        34.6623 E_vdwl   =      -248.7451 
+E_coul   =      -874.5122 E_long   =     -3336.4111 Press    =       135.8662 
+---------------- Step      200 ----- CPU =      7.9028 (sec) ----------------
+TotEng   =     -2660.1406 KinEng   =       626.3319 Temp     =       235.8265 
+PotEng   =     -3286.4725 E_bond   =       209.5147 E_angle  =       591.7773 
+E_dihed  =       388.9591 E_impro  =        29.4992 E_vdwl   =      -243.5808 
+E_coul   =      -923.5115 E_long   =     -3339.1306 Press    =        88.9000 
+---------------- Step      300 ----- CPU =     11.8246 (sec) ----------------
+TotEng   =     -2673.8090 KinEng   =       616.7924 Temp     =       232.2346 
+PotEng   =     -3290.6014 E_bond   =       202.8254 E_angle  =       568.6860 
+E_dihed  =       378.4182 E_impro  =        38.2399 E_vdwl   =      -221.3236 
+E_coul   =      -915.3004 E_long   =     -3342.1468 Press    =        78.8527 
+---------------- Step      400 ----- CPU =     15.7990 (sec) ----------------
+TotEng   =     -2614.9416 KinEng   =       649.3474 Temp     =       244.4922 
+PotEng   =     -3264.2890 E_bond   =       211.6116 E_angle  =       617.2026 
+E_dihed  =       399.8744 E_impro  =        40.2678 E_vdwl   =      -211.7790 
+E_coul   =      -978.1624 E_long   =     -3343.3041 Press    =        -4.1958 
+---------------- Step      500 ----- CPU =     19.8146 (sec) ----------------
+TotEng   =     -2588.6772 KinEng   =       660.1424 Temp     =       248.5568 
+PotEng   =     -3248.8196 E_bond   =       218.4786 E_angle  =       620.8605 
+E_dihed  =       390.3220 E_impro  =        41.6794 E_vdwl   =      -226.3657 
+E_coul   =      -953.1676 E_long   =     -3340.6269 Press    =        99.3200 
+---------------- Step      600 ----- CPU =     23.8587 (sec) ----------------
+TotEng   =     -2550.4618 KinEng   =       693.3384 Temp     =       261.0557 
+PotEng   =     -3243.8002 E_bond   =       232.3563 E_angle  =       606.2922 
+E_dihed  =       396.2469 E_impro  =        37.1980 E_vdwl   =      -235.8425 
+E_coul   =      -937.1208 E_long   =     -3342.9303 Press    =       -21.7737 
+---------------- Step      700 ----- CPU =     27.8381 (sec) ----------------
+TotEng   =     -2554.4355 KinEng   =       692.8951 Temp     =       260.8888 
+PotEng   =     -3247.3306 E_bond   =       216.3395 E_angle  =       637.7785 
+E_dihed  =       391.5940 E_impro  =        43.1426 E_vdwl   =      -187.6159 
+E_coul   =     -1008.1694 E_long   =     -3340.3998 Press    =        75.1484 
+---------------- Step      800 ----- CPU =     31.8039 (sec) ----------------
+TotEng   =     -2508.3551 KinEng   =       699.0766 Temp     =       263.2163 
+PotEng   =     -3207.4317 E_bond   =       241.9936 E_angle  =       641.3631 
+E_dihed  =       386.2198 E_impro  =        43.7793 E_vdwl   =      -217.7523 
+E_coul   =      -964.6070 E_long   =     -3338.4282 Press    =      -127.7337 
+---------------- Step      900 ----- CPU =     35.7700 (sec) ----------------
+TotEng   =     -2452.7644 KinEng   =       762.1842 Temp     =       286.9776 
+PotEng   =     -3214.9485 E_bond   =       243.9191 E_angle  =       649.8664 
+E_dihed  =       382.4351 E_impro  =        39.0029 E_vdwl   =      -221.3389 
+E_coul   =      -970.8965 E_long   =     -3337.9366 Press    =       122.7720 
+---------------- Step     1000 ----- CPU =     39.7695 (sec) ----------------
+TotEng   =     -2386.6805 KinEng   =       799.0253 Temp     =       300.8490 
+PotEng   =     -3185.7058 E_bond   =       265.3649 E_angle  =       661.7543 
+E_dihed  =       374.6843 E_impro  =        38.6877 E_vdwl   =      -229.2030 
+E_coul   =      -960.7041 E_long   =     -3336.2899 Press    =       -17.9910 
+Loop time of 39.7695 on 1 procs for 1000 steps with 892 atoms
+
+Performance: 17.380 ns/day, 1.381 hours/ns, 25.145 timesteps/s
+99.6% CPU use with 1 MPI tasks x 1 OpenMP threads
+
+MPI task timing breakdown:
+Section |  min time  |  avg time  |  max time  |%varavg| %total
+---------------------------------------------------------------
+Pair    | 29.169     | 29.169     | 29.169     |   0.0 | 73.34
+Bond    | 7.6249     | 7.6249     | 7.6249     |   0.0 | 19.17
+Kspace  | 1.1525     | 1.1525     | 1.1525     |   0.0 |  2.90
+Neigh   | 0.87606    | 0.87606    | 0.87606    |   0.0 |  2.20
+Comm    | 0.01563    | 0.01563    | 0.01563    |   0.0 |  0.04
+Output  | 0.00048423 | 0.00048423 | 0.00048423 |   0.0 |  0.00
+Modify  | 0.80446    | 0.80446    | 0.80446    |   0.0 |  2.02
+Other   |            | 0.1266     |            |       |  0.32
+
+Nlocal:    892 ave 892 max 892 min
+Histogram: 1 0 0 0 0 0 0 0 0 0
+Nghost:    27 ave 27 max 27 min
+Histogram: 1 0 0 0 0 0 0 0 0 0
+Neighs:    146206 ave 146206 max 146206 min
+Histogram: 1 0 0 0 0 0 0 0 0 0
+
+Total # of neighbors = 146206
+Ave neighs/atom = 163.908
+Ave special neighs/atom = 10.9395
+Neighbor list builds = 186
+Dangerous builds = 0
+
+unfix           cor
+unfix           1
+
+
+Please see the log.cite file for references relevant to this simulation
+
+Total wall time: 0:00:43
--- a/examples/USER/misc/filter-corotate/log.10Mar2017.bpti.g++.4
+++ b/examples/USER/misc/filter-corotate/log.10Mar2017.bpti.g++.4
@ -0,0 +1,240 @@
+LAMMPS (10 Mar 2017)
+  using 1 OpenMP thread(s) per MPI task
+
+units           real
+
+atom_style      full
+bond_style      harmonic
+angle_style     charmm
+dihedral_style  charmm
+improper_style  harmonic
+
+pair_style      lj/charmm/coul/long 8 10
+pair_modify     mix arithmetic
+kspace_style    pppm 1e-4
+
+read_data       data.bpti
+  orthogonal box = (-10 -10 -30) to (50 50 30)
+  1 by 2 by 2 MPI processor grid
+  reading atoms ...
+  892 atoms
+  scanning bonds ...
+  4 = max bonds/atom
+  scanning angles ...
+  6 = max angles/atom
+  scanning dihedrals ...
+  18 = max dihedrals/atom
+  scanning impropers ...
+  2 = max impropers/atom
+  reading bonds ...
+  906 bonds
+  reading angles ...
+  1626 angles
+  reading dihedrals ...
+  2501 dihedrals
+  reading impropers ...
+  137 impropers
+  4 = max # of 1-2 neighbors
+  9 = max # of 1-3 neighbors
+  19 = max # of 1-4 neighbors
+  21 = max # of special neighbors
+
+special_bonds   charmm
+neigh_modify    delay 2 every 1
+
+
+# ------------- MINIMIZE ----------
+
+minimize 	1e-4 1e-6 1000 10000
+WARNING: Resetting reneighboring criteria during minimization (../min.cpp:168)
+PPPM initialization ...
+WARNING: System is not charge neutral, net charge = 6 (../kspace.cpp:302)
+WARNING: Using 12-bit tables for long-range coulomb (../kspace.cpp:321)
+  G vector (1/distance) = 0.203272
+  grid = 16 16 16
+  stencil order = 5
+  estimated absolute RMS force accuracy = 0.0316399
+  estimated relative force accuracy = 9.52826e-05
+  using double precision FFTs
+  3d grid and FFT values/proc = 3549 1024
+Neighbor list info ...
+  update every 1 steps, delay 0 steps, check yes
+  max neighbors/atom: 2000, page size: 100000
+  master list distance cutoff = 12
+  ghost atom cutoff = 12
+  binsize = 6, bins = 10 10 10
+  1 neighbor lists, perpetual/occasional/extra = 1 0 0
+  (1) pair lj/charmm/coul/long, perpetual
+      attributes: half, newton on
+      pair build: half/bin/newton
+      stencil: half/bin/3d/newton
+      bin: standard
+Per MPI rank memory usage (min/avg/max) = 16.9693/0.981879/0 Mbytes
+Step Temp E_pair E_mol TotEng Press 
+       0            0   -3075.6498    943.91164   -2131.7381   -380.67776 
+     241            0   -4503.3131    749.58666   -3753.7264   -29.045153 
+Loop time of 1.26594 on 4 procs for 241 steps with 892 atoms
+
+99.0% CPU use with 4 MPI tasks x 1 OpenMP threads
+
+Minimization stats:
+  Stopping criterion = energy tolerance
+  Energy initial, next-to-last, final = 
+        -2131.73812515     -3753.43983927     -3753.72640137
+  Force two-norm initial, final = 1086.21 26.3688
+  Force max component initial, final = 310.811 3.92751
+  Final line search alpha, max atom move = 0.00596649 0.0234334
+  Iterations, force evaluations = 241 463
+
+MPI task timing breakdown:
+Section |  min time  |  avg time  |  max time  |%varavg| %total
+---------------------------------------------------------------
+Pair    | 0.34267    | 0.63792    | 0.90268    |  25.2 | 50.39
+Bond    | 0.025776   | 0.063318   | 0.095631   |  10.8 |  5.00
+Kspace  | 0.21904    | 0.51601    | 0.84895    |  31.3 | 40.76
+Neigh   | 0.023185   | 0.023363   | 0.023538   |   0.1 |  1.85
+Comm    | 0.012025   | 0.014189   | 0.016335   |   1.4 |  1.12
+Output  | 0          | 0          | 0          |   0.0 |  0.00
+Modify  | 0          | 0          | 0          |   0.0 |  0.00
+Other   |            | 0.01114    |            |       |  0.88
+
+Nlocal:    223 ave 323 max 89 min
+Histogram: 1 0 0 0 1 0 0 0 1 1
+Nghost:    613 ave 675 max 557 min
+Histogram: 1 0 0 1 0 1 0 0 0 1
+Neighs:    37222.8 ave 50005 max 20830 min
+Histogram: 1 0 0 0 1 0 0 1 0 1
+
+Total # of neighbors = 148891
+Ave neighs/atom = 166.918
+Ave special neighs/atom = 10.9395
+Neighbor list builds = 15
+Dangerous builds = 0
+reset_timestep  0
+
+# ------------- RUN ---------------
+
+thermo          100
+thermo_style    multi
+timestep	8
+
+run_style respa 3 2 8 bond 1 pair 2 kspace 3
+Respa levels:
+  1 = bond angle dihedral improper
+  2 = pair
+  3 = kspace
+
+velocity        all create 200.0 12345678 dist uniform
+#dump            dump1 all atom 100 4pti.dump
+
+fix             1 all nvt temp 200 300 25
+fix             cor all filter/corotate m 1.0
+  163 = # of size 2 clusters
+  0 = # of size 3 clusters
+  25 = # of size 4 clusters
+  0 = # of size 5 clusters
+  100 = # of frozen angles
+
+run             1000
+PPPM initialization ...
+WARNING: Using 12-bit tables for long-range coulomb (../kspace.cpp:321)
+  G vector (1/distance) = 0.203272
+  grid = 16 16 16
+  stencil order = 5
+  estimated absolute RMS force accuracy = 0.0316399
+  estimated relative force accuracy = 9.52826e-05
+  using double precision FFTs
+  3d grid and FFT values/proc = 3549 1024
+Per MPI rank memory usage (min/avg/max) = 17.142/0.97212/0 Mbytes
+---------------- Step        0 ----- CPU =      0.0000 (sec) ----------------
+TotEng   =     -3220.3378 KinEng   =       531.1804 Temp     =       200.0000 
+PotEng   =     -3751.5182 E_bond   =        42.2810 E_angle  =       345.2592 
+E_dihed  =       337.8361 E_impro  =        24.2103 E_vdwl   =      -288.5339 
+E_coul   =      -886.3622 E_long   =     -3326.2088 Press    =        83.2282 
+---------------- Step      100 ----- CPU =      1.5457 (sec) ----------------
+TotEng   =     -2718.9184 KinEng   =       538.6205 Temp     =       202.8014 
+PotEng   =     -3257.5389 E_bond   =       203.3365 E_angle  =       566.5311 
+E_dihed  =       397.6202 E_impro  =        34.6621 E_vdwl   =      -248.7451 
+E_coul   =      -874.5326 E_long   =     -3336.4111 Press    =       135.8435 
+---------------- Step      200 ----- CPU =      3.0720 (sec) ----------------
+TotEng   =     -2660.1146 KinEng   =       626.3474 Temp     =       235.8323 
+PotEng   =     -3286.4620 E_bond   =       209.5168 E_angle  =       591.7735 
+E_dihed  =       388.9615 E_impro  =        29.5000 E_vdwl   =      -243.5840 
+E_coul   =      -923.4998 E_long   =     -3339.1299 Press    =        88.8857 
+---------------- Step      300 ----- CPU =      4.5597 (sec) ----------------
+TotEng   =     -2669.7442 KinEng   =       619.3625 Temp     =       233.2023 
+PotEng   =     -3289.1067 E_bond   =       203.4405 E_angle  =       569.5281 
+E_dihed  =       378.3314 E_impro  =        38.2880 E_vdwl   =      -221.1904 
+E_coul   =      -915.3396 E_long   =     -3342.1646 Press    =        79.3780 
+---------------- Step      400 ----- CPU =      5.9808 (sec) ----------------
+TotEng   =     -2618.9975 KinEng   =       644.6145 Temp     =       242.7102 
+PotEng   =     -3263.6119 E_bond   =       209.5864 E_angle  =       618.8954 
+E_dihed  =       401.3798 E_impro  =        39.9064 E_vdwl   =      -212.1271 
+E_coul   =      -977.1589 E_long   =     -3344.0940 Press    =        -7.8938 
+---------------- Step      500 ----- CPU =      7.4159 (sec) ----------------
+TotEng   =     -2579.7486 KinEng   =       666.4643 Temp     =       250.9371 
+PotEng   =     -3246.2129 E_bond   =       219.2549 E_angle  =       620.3474 
+E_dihed  =       388.4395 E_impro  =        41.4499 E_vdwl   =      -225.9686 
+E_coul   =      -949.3689 E_long   =     -3340.3672 Press    =       113.2543 
+---------------- Step      600 ----- CPU =      8.9252 (sec) ----------------
+TotEng   =     -2535.8235 KinEng   =       708.5919 Temp     =       266.7990 
+PotEng   =     -3244.4154 E_bond   =       243.9451 E_angle  =       606.0866 
+E_dihed  =       400.0562 E_impro  =        33.9708 E_vdwl   =      -223.1319 
+E_coul   =      -964.9940 E_long   =     -3340.3482 Press    =      -102.4475 
+---------------- Step      700 ----- CPU =     10.4022 (sec) ----------------
+TotEng   =     -2552.6681 KinEng   =       702.3080 Temp     =       264.4330 
+PotEng   =     -3254.9761 E_bond   =       250.8834 E_angle  =       639.0977 
+E_dihed  =       386.4014 E_impro  =        42.3004 E_vdwl   =      -224.4816 
+E_coul   =     -1011.8551 E_long   =     -3337.3222 Press    =        10.6424 
+---------------- Step      800 ----- CPU =     11.8699 (sec) ----------------
+TotEng   =     -2423.5415 KinEng   =       772.1254 Temp     =       290.7206 
+PotEng   =     -3195.6670 E_bond   =       238.5831 E_angle  =       640.9180 
+E_dihed  =       377.7994 E_impro  =        40.3135 E_vdwl   =      -216.5705 
+E_coul   =      -935.1087 E_long   =     -3341.6019 Press    =       -38.2479 
+---------------- Step      900 ----- CPU =     13.3548 (sec) ----------------
+TotEng   =     -2394.4779 KinEng   =       766.6895 Temp     =       288.6739 
+PotEng   =     -3161.1673 E_bond   =       284.8428 E_angle  =       671.0959 
+E_dihed  =       380.3406 E_impro  =        51.2975 E_vdwl   =      -219.5211 
+E_coul   =      -990.6305 E_long   =     -3338.5925 Press    =       -15.2279 
+---------------- Step     1000 ----- CPU =     14.7908 (sec) ----------------
+TotEng   =     -2340.1471 KinEng   =       799.0198 Temp     =       300.8469 
+PotEng   =     -3139.1669 E_bond   =       271.0389 E_angle  =       683.8278 
+E_dihed  =       407.0795 E_impro  =        39.6209 E_vdwl   =      -230.5355 
+E_coul   =      -974.2981 E_long   =     -3335.9003 Press    =       -94.3420 
+Loop time of 14.7909 on 4 procs for 1000 steps with 892 atoms
+
+Performance: 46.732 ns/day, 0.514 hours/ns, 67.609 timesteps/s
+99.1% CPU use with 4 MPI tasks x 1 OpenMP threads
+
+MPI task timing breakdown:
+Section |  min time  |  avg time  |  max time  |%varavg| %total
+---------------------------------------------------------------
+Pair    | 4.4184     | 7.5543     | 10.133     |  74.2 | 51.07
+Bond    | 0.94027    | 1.9781     | 2.7492     |  54.4 | 13.37
+Kspace  | 0.45487    | 0.45887    | 0.46343    |   0.4 |  3.10
+Neigh   | 0.28145    | 0.28339    | 0.28539    |   0.3 |  1.92
+Comm    | 0.7515     | 4.1484     | 8.3861     | 135.5 | 28.05
+Output  | 0.00049973 | 0.00055474 | 0.00066924 |   0.0 |  0.00
+Modify  | 0.26165    | 0.31142    | 0.35023    |   6.7 |  2.11
+Other   |            | 0.05572    |            |       |  0.38
+
+Nlocal:    223 ave 313 max 122 min
+Histogram: 1 0 0 1 0 0 0 1 0 1
+Nghost:    584.5 ave 605 max 553 min
+Histogram: 1 0 0 0 0 1 0 0 0 2
+Neighs:    35448 ave 42093 max 25175 min
+Histogram: 1 0 0 0 0 0 1 1 0 1
+
+Total # of neighbors = 141792
+Ave neighs/atom = 158.96
+Ave special neighs/atom = 10.9395
+Neighbor list builds = 186
+Dangerous builds = 0
+
+unfix           cor
+unfix           1
+
+
+Please see the log.cite file for references relevant to this simulation
+
+Total wall time: 0:00:16
--- a/examples/USER/misc/filter-corotate/log.10Mar2017.peptide.g++.1
+++ b/examples/USER/misc/filter-corotate/log.10Mar2017.peptide.g++.1
@ -0,0 +1,146 @@
+LAMMPS (10 Mar 2017)
+  using 1 OpenMP thread(s) per MPI task
+# Solvated 5-mer peptide, run for 8ps in NVT
+
+units           real
+atom_style      full
+
+pair_style      lj/charmm/coul/long 8.0 10.0 10.0
+bond_style      harmonic
+angle_style     charmm
+dihedral_style  charmm
+improper_style  harmonic
+kspace_style    pppm 0.0001
+
+read_data       data.peptide
+  orthogonal box = (36.8402 41.0137 29.7681) to (64.2116 68.3851 57.1395)
+  1 by 1 by 1 MPI processor grid
+  reading atoms ...
+  2004 atoms
+  reading velocities ...
+  2004 velocities
+  scanning bonds ...
+  3 = max bonds/atom
+  scanning angles ...
+  6 = max angles/atom
+  scanning dihedrals ...
+  14 = max dihedrals/atom
+  scanning impropers ...
+  1 = max impropers/atom
+  reading bonds ...
+  1365 bonds
+  reading angles ...
+  786 angles
+  reading dihedrals ...
+  207 dihedrals
+  reading impropers ...
+  12 impropers
+  4 = max # of 1-2 neighbors
+  7 = max # of 1-3 neighbors
+  14 = max # of 1-4 neighbors
+  18 = max # of special neighbors
+
+neighbor        2.0 bin
+neigh_modify    delay 5
+
+thermo          50
+#dump            dump1 all atom 100 peptide.dump
+
+timestep        8
+
+run_style respa 3 2 8 bond 1 pair 2 kspace 3
+Respa levels:
+  1 = bond angle dihedral improper
+  2 = pair
+  3 = kspace
+
+fix             1 all nvt temp 250.0 250.0 100.0 tchain 1
+fix             cor all filter/corotate m 1.0
+  19 = # of size 2 clusters
+  0 = # of size 3 clusters
+  3 = # of size 4 clusters
+  0 = # of size 5 clusters
+  646 = # of frozen angles
+run             1000
+PPPM initialization ...
+WARNING: Using 12-bit tables for long-range coulomb (../kspace.cpp:321)
+  G vector (1/distance) = 0.268725
+  grid = 15 15 15
+  stencil order = 5
+  estimated absolute RMS force accuracy = 0.0228209
+  estimated relative force accuracy = 6.87243e-05
+  using double precision FFTs
+  3d grid and FFT values/proc = 10648 3375
+Neighbor list info ...
+  update every 1 steps, delay 5 steps, check yes
+  max neighbors/atom: 2000, page size: 100000
+  master list distance cutoff = 12
+  ghost atom cutoff = 12
+  binsize = 6, bins = 5 5 5
+  1 neighbor lists, perpetual/occasional/extra = 1 0 0
+  (1) pair lj/charmm/coul/long, perpetual
+      attributes: half, newton on
+      pair build: half/bin/newton
+      stencil: half/bin/3d/newton
+      bin: standard
+Per MPI rank memory usage (min/avg/max) = 22.6706/1/0 Mbytes
+Step Temp E_pair E_mol TotEng Press 
+       0     190.0857   -6785.6785    70.391457   -5580.3684    19434.821 
+      50    239.46028   -7546.5667    1092.8874   -5023.9668   -24643.891 
+     100    242.81799   -7125.5527     416.0788   -5259.7139    15525.465 
+     150    235.97108   -7531.9334    932.35464   -5190.6987   -14838.489 
+     200    252.06415   -7195.6011    568.02993   -5122.6064     8841.332 
+     250    249.99431   -7586.5092    881.83491   -5212.0676    -9330.345 
+     300     240.3382   -7333.0933    633.29951   -5264.8395    5137.9757 
+     350    255.34529   -7568.2413    856.46371   -5187.2226    -6206.063 
+     400    242.99276   -7419.9031    713.23943   -5255.8602    2447.0091 
+     450    251.10653    -7622.061    844.20584   -5278.6079   -4906.6559 
+     500    255.59314    -7439.253    710.84907   -5202.3691    1571.0032 
+     550     253.2025   -7660.5101    823.05373    -5325.695    -4551.399 
+     600    249.05313   -7509.6729    741.48104   -5281.2046       992.87 
+     650    251.75984   -7593.6589    847.08244   -5243.4286   -3510.1176 
+     700    249.25027   -7601.9112     794.0912   -5319.6557    305.76021 
+     750      255.415   -7602.2674    822.98524   -5254.3109    -2333.421 
+     800    241.99621   -7643.8878    796.53352   -5402.5008   -298.66565 
+     850     253.6428   -7598.3764    816.45457   -5267.5316   -1905.3478 
+     900    247.20231   -7690.2806    789.75999   -5424.5838   -1331.7228 
+     950    255.92583   -7634.7505    831.18272   -5275.5466   -2186.5117 
+    1000     253.2126   -7647.9526    823.93602    -5312.195   -1189.9659 
+Loop time of 150.664 on 1 procs for 1000 steps with 2004 atoms
+
+Performance: 4.588 ns/day, 5.231 hours/ns, 6.637 timesteps/s
+99.7% CPU use with 1 MPI tasks x 1 OpenMP threads
+
+MPI task timing breakdown:
+Section |  min time  |  avg time  |  max time  |%varavg| %total
+---------------------------------------------------------------
+Pair    | 135.81     | 135.81     | 135.81     |   0.0 | 90.14
+Bond    | 2.5889     | 2.5889     | 2.5889     |   0.0 |  1.72
+Kspace  | 2.0379     | 2.0379     | 2.0379     |   0.0 |  1.35
+Neigh   | 5.893      | 5.893      | 5.893      |   0.0 |  3.91
+Comm    | 1.6998     | 1.6998     | 1.6998     |   0.0 |  1.13
+Output  | 0.00077915 | 0.00077915 | 0.00077915 |   0.0 |  0.00
+Modify  | 2          | 2          | 2          |   0.0 |  1.33
+Other   |            | 0.6352     |            |       |  0.42
+
+Nlocal:    2004 ave 2004 max 2004 min
+Histogram: 1 0 0 0 0 0 0 0 0 0
+Nghost:    11197 ave 11197 max 11197 min
+Histogram: 1 0 0 0 0 0 0 0 0 0
+Neighs:    707779 ave 707779 max 707779 min
+Histogram: 1 0 0 0 0 0 0 0 0 0
+
+Total # of neighbors = 707779
+Ave neighs/atom = 353.183
+Ave special neighs/atom = 2.34032
+Neighbor list builds = 200
+Dangerous builds = 200
+unfix           cor
+unfix           1
+
+
+
+
+Please see the log.cite file for references relevant to this simulation
+
+Total wall time: 0:02:30
--- a/examples/USER/misc/filter-corotate/log.10Mar2017.peptide.g++.4
+++ b/examples/USER/misc/filter-corotate/log.10Mar2017.peptide.g++.4
@ -0,0 +1,146 @@
+LAMMPS (10 Mar 2017)
+  using 1 OpenMP thread(s) per MPI task
+# Solvated 5-mer peptide, run for 8ps in NVT
+
+units           real
+atom_style      full
+
+pair_style      lj/charmm/coul/long 8.0 10.0 10.0
+bond_style      harmonic
+angle_style     charmm
+dihedral_style  charmm
+improper_style  harmonic
+kspace_style    pppm 0.0001
+
+read_data       data.peptide
+  orthogonal box = (36.8402 41.0137 29.7681) to (64.2116 68.3851 57.1395)
+  1 by 2 by 2 MPI processor grid
+  reading atoms ...
+  2004 atoms
+  reading velocities ...
+  2004 velocities
+  scanning bonds ...
+  3 = max bonds/atom
+  scanning angles ...
+  6 = max angles/atom
+  scanning dihedrals ...
+  14 = max dihedrals/atom
+  scanning impropers ...
+  1 = max impropers/atom
+  reading bonds ...
+  1365 bonds
+  reading angles ...
+  786 angles
+  reading dihedrals ...
+  207 dihedrals
+  reading impropers ...
+  12 impropers
+  4 = max # of 1-2 neighbors
+  7 = max # of 1-3 neighbors
+  14 = max # of 1-4 neighbors
+  18 = max # of special neighbors
+
+neighbor        2.0 bin
+neigh_modify    delay 5
+
+thermo          50
+#dump            dump1 all atom 100 peptide.dump
+
+timestep        8
+
+run_style respa 3 2 8 bond 1 pair 2 kspace 3
+Respa levels:
+  1 = bond angle dihedral improper
+  2 = pair
+  3 = kspace
+
+fix             1 all nvt temp 250.0 250.0 100.0 tchain 1
+fix             cor all filter/corotate m 1.0
+  19 = # of size 2 clusters
+  0 = # of size 3 clusters
+  3 = # of size 4 clusters
+  0 = # of size 5 clusters
+  646 = # of frozen angles
+run             1000
+PPPM initialization ...
+WARNING: Using 12-bit tables for long-range coulomb (../kspace.cpp:321)
+  G vector (1/distance) = 0.268725
+  grid = 15 15 15
+  stencil order = 5
+  estimated absolute RMS force accuracy = 0.0228209
+  estimated relative force accuracy = 6.87243e-05
+  using double precision FFTs
+  3d grid and FFT values/proc = 4312 960
+Neighbor list info ...
+  update every 1 steps, delay 5 steps, check yes
+  max neighbors/atom: 2000, page size: 100000
+  master list distance cutoff = 12
+  ghost atom cutoff = 12
+  binsize = 6, bins = 5 5 5
+  1 neighbor lists, perpetual/occasional/extra = 1 0 0
+  (1) pair lj/charmm/coul/long, perpetual
+      attributes: half, newton on
+      pair build: half/bin/newton
+      stencil: half/bin/3d/newton
+      bin: standard
+Per MPI rank memory usage (min/avg/max) = 16.8394/0.98826/0 Mbytes
+Step Temp E_pair E_mol TotEng Press 
+       0     190.0857   -6785.6785    70.391457   -5580.3684    19434.821 
+      50    239.46028   -7546.5668    1092.8874   -5023.9668   -24643.891 
+     100    242.81819   -7125.5629    416.08082   -5259.7209    15525.244 
+     150    235.94928   -7531.9186    932.50658   -5190.6621   -14842.431 
+     200    255.85551   -7254.4065     568.8803   -5157.9249    8936.8651 
+     250     247.8705   -7607.4583    858.06087   -5269.4711   -9926.0442 
+     300    257.64176    -7267.424     618.5573   -5110.6004    5173.3307 
+     350    251.65439   -7572.3806    821.15745   -5248.7049    -7092.327 
+     400    256.87927   -7414.2145    655.33178    -5225.169    4119.4095 
+     450    257.12393   -7576.5541    853.39773   -5187.9819   -5224.8823 
+     500    242.42371    -7524.705    705.75357   -5371.5455    2111.3878 
+     550    248.97188    -7541.076    792.86994   -5261.7038   -2278.4185 
+     600    249.81862   -7592.0499    767.17722   -5333.3149   -1149.4759 
+     650    253.31349   -7578.2665    813.75975   -5252.0827   -2915.5706 
+     700    256.61152   -7588.1475    761.03356   -5294.9988   -747.88089 
+     750     248.3606    -7660.457    837.71615   -5339.8883   -3072.8311 
+     800    253.81464   -7638.6089     782.4229   -5340.7698    -1025.909 
+     850    245.69185   -7660.9036    795.66792   -5398.3172   -2717.5851 
+     900    249.13156   -7589.4769    806.43464   -5295.5867   -761.63361 
+     950    251.11482   -7691.4981    869.34937    -5322.852   -3282.3031 
+    1000     241.9195   -7630.9899    828.59107   -5358.0033   -95.962685 
+Loop time of 45.5507 on 4 procs for 1000 steps with 2004 atoms
+
+Performance: 15.174 ns/day, 1.582 hours/ns, 21.954 timesteps/s
+99.4% CPU use with 4 MPI tasks x 1 OpenMP threads
+
+MPI task timing breakdown:
+Section |  min time  |  avg time  |  max time  |%varavg| %total
+---------------------------------------------------------------
+Pair    | 35.545     | 36.674     | 38.004     |  15.8 | 80.51
+Bond    | 0.51302    | 0.67796    | 0.86345    |  18.6 |  1.49
+Kspace  | 0.66031    | 0.68459    | 0.70506    |   2.1 |  1.50
+Neigh   | 1.5605     | 1.5627     | 1.5649     |   0.1 |  3.43
+Comm    | 3.4611     | 4.9841     | 6.294      |  47.2 | 10.94
+Output  | 0.00079799 | 0.00086641 | 0.0010369  |   0.0 |  0.00
+Modify  | 0.67341    | 0.69059    | 0.71186    |   1.7 |  1.52
+Other   |            | 0.2762     |            |       |  0.61
+
+Nlocal:    501 ave 523 max 473 min
+Histogram: 1 0 0 0 0 0 2 0 0 1
+Nghost:    6643.25 ave 6708 max 6566 min
+Histogram: 1 1 0 0 0 0 0 0 0 2
+Neighs:    176977 ave 185765 max 164931 min
+Histogram: 1 0 0 0 1 0 0 0 1 1
+
+Total # of neighbors = 707908
+Ave neighs/atom = 353.248
+Ave special neighs/atom = 2.34032
+Neighbor list builds = 200
+Dangerous builds = 200
+unfix           cor
+unfix           1
+
+
+
+
+Please see the log.cite file for references relevant to this simulation
+
+Total wall time: 0:00:45
--- a/python/examples/matplotlib_plot.py
+++ b/python/examples/matplotlib_plot.py
@ -0,0 +1,93 @@
+#!/usr/bin/env python -i
+# preceding line should have path for Python on your machine
+
+# matplotlib_plot.py
+# Purpose: plot Temp of running LAMMPS simulation via matplotlib
+# Syntax:  plot.py in.lammps Nfreq Nsteps compute-ID
+#          in.lammps = LAMMPS input script
+#          Nfreq = plot data point every this many steps
+#          Nsteps = run for this many steps
+#          compute-ID = ID of compute that calculates temperature
+#                       (or any other scalar quantity)
+
+from __future__ import print_function
+import sys
+sys.path.append("./pizza")
+import matplotlib
+matplotlib.use('tkagg')
+import matplotlib.pyplot as plt
+
+# parse command line
+
+argv = sys.argv
+if len(argv) != 5:
+  print("Syntax: plot.py in.lammps Nfreq Nsteps compute-ID")
+  sys.exit()
+
+infile = sys.argv[1]
+nfreq = int(sys.argv[2])
+nsteps = int(sys.argv[3])
+compute = sys.argv[4]
+
+me = 0
+# uncomment if running in parallel via Pypar
+#import pypar
+#me = pypar.rank()
+#nprocs = pypar.size()
+
+from lammps import lammps
+lmp = lammps()
+
+# run infile all at once
+# assumed to have no run command in it
+
+lmp.file(infile)
+lmp.command("thermo %d" % nfreq)
+
+# initial 0-step run to generate initial 1-point plot
+
+lmp.command("run 0 pre yes post no")
+value = lmp.extract_compute(compute,0,0)
+ntimestep = 0
+xaxis = [ntimestep]
+yaxis = [value]
+
+# create matplotlib plot
+# just proc 0 handles plotting
+
+if me == 0:
+  fig = plt.figure()
+  line, = plt.plot(xaxis, yaxis)
+  plt.xlim([0, nsteps])
+  plt.title(compute)
+  plt.xlabel("Timestep")
+  plt.ylabel("Temperature")
+  plt.show(block=False)
+
+# run nfreq steps at a time w/out pre/post, query compute, refresh plot
+import time
+
+while ntimestep < nsteps:
+  lmp.command("run %d pre no post no" % nfreq)
+  ntimestep += nfreq
+  value = lmp.extract_compute(compute,0,0)
+  xaxis.append(ntimestep)
+  yaxis.append(value)
+  if me == 0:
+    line.set_xdata(xaxis)
+    line.set_ydata(yaxis)
+    ax = plt.gca()
+    ax.relim()
+    ax.autoscale_view(True, True, True)
+    fig.canvas.draw()
+
+lmp.command("run 0 pre no post yes")
+
+# uncomment if running in parallel via Pypar
+#print("Proc %d out of %d procs has" % (me,nprocs), lmp)
+#pypar.finalize()
+
+if sys.version_info[0] == 3:
+    input("Press Enter to exit...")
+else:
+    raw_input("Press Enter to exit...")
--- a/python/examples/pizza/dump.py
+++ b/python/examples/pizza/dump.py
@ -146,6 +146,8 @@ d.extra(data)				   extract bond/tri/line list from data
 # History
 #   8/05, Steve Plimpton (SNL): original version
 #   12/09, David Hart (SNL): allow use of NumPy or Numeric
+#   03/17, Richard Berger (Temple U): improve Python 3 compatibility,
+#                                     simplify read_snapshot by using reshape

 # ToDo list
 #   try to optimize this line in read_snap: words += f.readline().split()
@ -224,7 +226,7 @@ class dump:
    self.flist = []
    for word in words: self.flist += glob.glob(word)
    if len(self.flist) == 0 and len(list) == 1:
-      raise StandardError("no dump file specified")
+      raise Exception("no dump file specified")
    
    if len(list) == 1:
      self.increment = 0
@ -299,7 +301,7 @@ class dump:

  def next(self):

-    if not self.increment: raise StandardError("cannot read incrementally")
+    if not self.increment: raise Exception("cannot read incrementally")

    # read next snapshot in current file using eof as pointer
    # if fail, try next file
@ -344,13 +346,13 @@ class dump:
    try:
      snap = Snap()
      item = f.readline()
-      snap.time = int(f.readline().split()[0])    # just grab 1st field
+      snap.time = int(f.readline().decode().split()[0])    # just grab 1st field
      item = f.readline()
-      snap.natoms = int(f.readline())
+      snap.natoms = int(f.readline().decode())

      snap.aselect = np.zeros(snap.natoms)

-      item = f.readline()
+      item = f.readline().decode()
      words = f.readline().split()
      snap.xlo,snap.xhi = float(words[0]),float(words[1])
      words = f.readline().split()
@ -358,7 +360,7 @@ class dump:
      words = f.readline().split()
      snap.zlo,snap.zhi = float(words[0]),float(words[1])

-      item = f.readline()
+      item = f.readline().decode()
      if len(self.names) == 0:
        words = item.split()[2:]
        if len(words):
@ -372,24 +374,22 @@ class dump:
            else: self.names[words[i]] = i

      if snap.natoms:
-        words = f.readline().split()
+        words = f.readline().decode().split()
        ncol = len(words)
        for i in range(1,snap.natoms):
-          words += f.readline().split()
+          words += f.readline().decode().split()
        floats = map(float,words)
-        if oldnumeric: atoms = np.zeros((snap.natoms,ncol),np.Float)
-        else: atoms = np.zeros((snap.natoms,ncol),np.float)
-        start = 0
-        stop = ncol
-        for i in range(snap.natoms):
-          atoms[i] = floats[start:stop]
-          start = stop
-          stop += ncol
-      else: atoms = None
-      snap.atoms = atoms
+        if oldnumeric:
+          atom_data = np.array(list(floats),np.Float)
+        else:
+          atom_data = np.array(list(floats),np.float)
+
+        snap.atoms = atom_data.reshape((snap.natoms, ncol))
+      else:
+        snap.atoms = None
      return snap
    except:
-      return 0
+      return None

  # --------------------------------------------------------------------
  # decide if snapshot i is scaled/unscaled from coords of first and last atom
@ -417,7 +417,7 @@ class dump:
  
  def map(self,*pairs):
    if len(pairs) % 2 != 0:
-      raise StandardError("dump map() requires pairs of mappings")
+      raise Exception("dump map() requires pairs of mappings")
    for i in range(0,len(pairs),2):
      j = i + 1
      self.names[pairs[j]] = pairs[i]-1
@ -734,7 +734,7 @@ class dump:
    for snap in self.snaps:
      if not snap.tselect: continue
      if snap.nselect != len(vec):
-        raise StandardError("vec length does not match # of selected atoms")
+        raise Exception("vec length does not match # of selected atoms")
      atoms = snap.atoms
      m = 0
      for i in range(snap.natoms):
@ -800,7 +800,7 @@ class dump:

  def atom(self,n,*list):
    if len(list) == 0:
-      raise StandardError("no columns specified")
+      raise Exception("no columns specified")
    columns = []
    values = []
    for name in list:
@ -816,7 +816,7 @@ class dump:
      for i in range(snap.natoms):
        if atoms[i][id] == n: break
      if atoms[i][id] != n:
-        raise StandardError("could not find atom ID in snapshot")
+        raise Exception("could not find atom ID in snapshot")
      for j in range(ncol):
        values[j][m] = atoms[i][columns[j]]
      m += 1
@ -831,7 +831,7 @@ class dump:
    snap = self.snaps[self.findtime(n)]
    
    if len(list) == 0:
-      raise StandardError("no columns specified")
+      raise Exception("no columns specified")
    columns = []
    values = []
    for name in list:
@ -957,9 +957,9 @@ class dump:
  # --------------------------------------------------------------------

  def findtime(self,n):
-    for i in range(self.nsnaps):
-      if self.snaps[i].time == n: return i
-    raise StandardError("no step %d exists" % n)
+    for i, snap in enumerate(self.snaps):
+      if snap.time == n: return i
+    raise Exception("no step %d exists" % n)

  # --------------------------------------------------------------------
  # return maximum box size across all selected snapshots
@ -1008,7 +1008,7 @@ class dump:
        nbonds = int(f.readline())
        item = f.readline()
        if not re.search("BONDS",item):
-          raise StandardError("could not read bonds from dump file")
+          raise Exception("could not read bonds from dump file")

        words = f.readline().split()
        ncol = len(words)
@ -1031,7 +1031,7 @@ class dump:
          self.bondflag = 1
          self.bondlist = bondlist
      except:
-        raise StandardError("could not read from bond dump file")
+        raise Exception("could not read from bond dump file")
      
    # request bonds from data object
    
@ -1047,7 +1047,7 @@ class dump:
          self.bondflag = 1
          self.bondlist = bondlist
      except:
-        raise StandardError("could not extract bonds from data object")
+        raise Exception("could not extract bonds from data object")

    # request tris/lines from cdata object
    
@ -1061,7 +1061,7 @@ class dump:
          self.lineflag = 1
          self.linelist = lines
      except:
-        raise StandardError("could not extract tris/lines from cdata object")
+        raise Exception("could not extract tris/lines from cdata object")

    # request tris from mdump object
    
@ -1070,10 +1070,10 @@ class dump:
        self.triflag = 2
        self.triobj = arg
      except:
-        raise StandardError("could not extract tris from mdump object")
+        raise Exception("could not extract tris from mdump object")

    else:
-      raise StandardError("unrecognized argument to dump.extra()")
+      raise Exception("unrecognized argument to dump.extra()")
      
  # --------------------------------------------------------------------

--- a/python/examples/pizza/gl.py
+++ b/python/examples/pizza/gl.py
@ -638,7 +638,7 @@ class gl:
                          fraction*(self.scale_stop - self.scale_start)
          self.viewupright()

-	if n == nstart or self.panflag: self.center = compute_center(box)
+        if n == nstart or self.panflag: self.center = compute_center(box)

        if bonds: self.bonds_augment(bonds)

--- a/python/examples/pizza/pdbfile.py
+++ b/python/examples/pizza/pdbfile.py
@ -54,6 +54,7 @@ index,time,flag = p.iterator(1)

 # History
 #   8/05, Steve Plimpton (SNL): original version
+#   3/17, Richard Berger (Temple U): improve Python 3 compatibility

 # ToDo list
 #   for generic PDB file (no template) from a LJ unit system,
@ -68,6 +69,12 @@ index,time,flag = p.iterator(1)
 # Imports and external programs

 import sys, types, glob, urllib
+PY3 = sys.version_info[0] == 3
+
+if PY3:
+    string_types = str,
+else:
+    string_types = basestring

 # Class definition

@ -77,7 +84,7 @@ class pdbfile:

  def __init__(self,*args):
    if len(args) == 1:
-      if type(args[0]) is types.StringType:
+      if type(args[0]) is string_types:
        filestr = args[0]
        self.data = None
      else:
--- a/python/examples/viz_atomeye.py
+++ b/python/examples/viz_atomeye.py
@ -42,7 +42,7 @@ lmp = lammps()

 lmp.file(infile)
 lmp.command("thermo %d" % nfreq)
-lmp.command("dump python all cfg %d tmp.cfg.* id type xs ys zs" % nfreq)
+lmp.command("dump python all cfg %d tmp.cfg.* mass type xs ys zs id" % nfreq)

 # initial 0-step run to generate dump file and image

--- a/python/examples/vizplotgui_atomeye.py
+++ b/python/examples/vizplotgui_atomeye.py
@ -73,7 +73,7 @@ lmp = lammps()

 lmp.file(infile)
 lmp.command("thermo %d" % nfreq)
-lmp.command("dump python all cfg %d tmp.cfg.* id type xs ys zs" % nfreq)
+lmp.command("dump python all cfg %d tmp.cfg.* mass type xs ys zs id" % nfreq)

 # initial 0-step run to generate initial 1-point plot, dump file, and image

--- a/python/lammps.py
+++ b/python/lammps.py
@ -190,12 +190,15 @@ class lammps(object):
  # send a list of commands

  def commands_list(self,cmdlist):
-    args = (c_char_p * len(cmdlist))(*cmdlist)
+    cmds = [x.encode() for x in cmdlist if type(x) is str]
+    args = (c_char_p * len(cmdlist))(*cmds)
    self.lib.lammps_commands_list(self.lmp,len(cmdlist),args)
    
  # send a string of commands

  def commands_string(self,multicmd):
+    if type(multicmd) is str:
+        multicmd = multicmd.encode()
    self.lib.lammps_commands_string(self.lmp,c_char_p(multicmd))
    
  # extract global info
@ -359,19 +362,27 @@ class lammps(object):
  #   e.g. for Python list or NumPy, etc
  #   ditto for gather_atoms() above

-  def create_atoms(self,n,id,type,x,v):
+  def create_atoms(self,n,id,type,x,v,image=None,shrinkexceed=False):
    if id:
      id_lmp = (c_int * n)()
      id_lmp[:] = id
-    else: id_lmp = id
+    else:
+      id_lmp = id
+
+    if image:
+      image_lmp = (c_int * n)()
+      image_lmp[:] = image
+    else:
+      image_lmp = image
+
    type_lmp = (c_int * n)()
    type_lmp[:] = type
-    self.lib.lammps_create_atoms(self.lmp,n,id_lmp,type_lmp,x,v)
+    self.lib.lammps_create_atoms(self.lmp,n,id_lmp,type_lmp,x,v,image_lmp,shrinkexceed)
+

-  # document this?
-    
  @property
  def uses_exceptions(self):
+    """ Return whether the LAMMPS shared library was compiled with C++ exceptions handling enabled """
    try:
      if self.lib.lammps_has_error:
        return True
--- a/src/.gitignore
+++ b/src/.gitignore
@ -328,6 +328,8 @@
 /fix_eos_table.h
 /fix_evaporate.cpp
 /fix_evaporate.h
+/fix_filter_corotate.cpp
+/fix_filter_corotate.h
 /fix_viscosity.cpp
 /fix_viscosity.h
 /fix_ehex.cpp
@ -653,6 +655,8 @@
 /pair_hbond_dreiding_lj.h
 /pair_hbond_dreiding_morse.cpp
 /pair_hbond_dreiding_morse.h
+/pair_kolmogorov_crespi_z.cpp
+/pair_kolmogorov_crespi_z.h
 /pair_lcbop.cpp
 /pair_lcbop.h
 /pair_line_lj.cpp
--- a/src/KOKKOS/fix_qeq_reax_kokkos.cpp
+++ b/src/KOKKOS/fix_qeq_reax_kokkos.cpp
@ -82,7 +82,7 @@ void FixQEqReaxKokkos<DeviceType>::init()

  FixQEqReax::init();

-  neighflag = lmp->kokkos->neighflag;
+  neighflag = lmp->kokkos->neighflag_qeq;
  int irequest = neighbor->nrequest - 1;
  
  neighbor->requests[irequest]->
--- a/src/KOKKOS/kokkos.cpp
+++ b/src/KOKKOS/kokkos.cpp
@ -119,6 +119,8 @@ KokkosLMP::KokkosLMP(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp)
  // default settings for package kokkos command

  neighflag = FULL;
+  neighflag_qeq = FULL;
+  neighflag_qeq_set = 0;
  exchange_comm_classic = 0;
  forward_comm_classic = 0;
  exchange_comm_on_host = 0;
@ -152,6 +154,8 @@ void KokkosLMP::accelerator(int narg, char **arg)
  // defaults

  neighflag = FULL;
+  neighflag_qeq = FULL;
+  neighflag_qeq_set = 0;
  int newtonflag = 0;
  double binsize = 0.0;
  exchange_comm_classic = forward_comm_classic = 0;
@ -169,6 +173,19 @@ void KokkosLMP::accelerator(int narg, char **arg)
          neighflag = HALF;
      } else if (strcmp(arg[iarg+1],"n2") == 0) neighflag = N2;
      else error->all(FLERR,"Illegal package kokkos command");
+      if (!neighflag_qeq_set) neighflag_qeq = neighflag;
+      iarg += 2;
+    } else if (strcmp(arg[iarg],"neigh/qeq") == 0) {
+      if (iarg+2 > narg) error->all(FLERR,"Illegal package kokkos command");
+      if (strcmp(arg[iarg+1],"full") == 0) neighflag_qeq = FULL;
+      else if (strcmp(arg[iarg+1],"half") == 0) {
+        if (num_threads > 1 || ngpu > 0)
+          neighflag_qeq = HALFTHREAD;
+        else 
+          neighflag_qeq = HALF;
+      } else if (strcmp(arg[iarg+1],"n2") == 0) neighflag_qeq = N2;
+      else error->all(FLERR,"Illegal package kokkos command");
+      neighflag_qeq_set = 1;
      iarg += 2;
    } else if (strcmp(arg[iarg],"binsize") == 0) {
      if (iarg+2 > narg) error->all(FLERR,"Illegal package kokkos command");
--- a/src/KOKKOS/kokkos.h
+++ b/src/KOKKOS/kokkos.h
@ -23,6 +23,8 @@ class KokkosLMP : protected Pointers {
 public:
  int kokkos_exists;
  int neighflag;
+  int neighflag_qeq;
+  int neighflag_qeq_set;
  int exchange_comm_classic;
  int forward_comm_classic;
  int exchange_comm_on_host;
--- a/src/RIGID/fix_shake.cpp
+++ b/src/RIGID/fix_shake.cpp
@ -66,6 +66,7 @@ FixShake::FixShake(LAMMPS *lmp, int narg, char **arg) :
  virial_flag = 1;
  create_attribute = 1;
  dof_flag = 1;
+
  // error check

  molecular = atom->molecular;
--- a/src/STUBS/mpi.c
+++ b/src/STUBS/mpi.c
@ -190,6 +190,13 @@ int MPI_Type_size(MPI_Datatype datatype, int *size)

 /* ---------------------------------------------------------------------- */

+int MPI_Request_free(MPI_Request *request)
+{
+  return 0;
+}
+
+/* ---------------------------------------------------------------------- */
+
 int MPI_Send(const void *buf, int count, MPI_Datatype datatype,
             int dest, int tag, MPI_Comm comm)
 {
--- a/src/STUBS/mpi.h
+++ b/src/STUBS/mpi.h
@ -89,6 +89,7 @@ int MPI_Finalize();
 double MPI_Wtime();

 int MPI_Type_size(int, int *);
+int MPI_Request_free(MPI_Request *request);

 int MPI_Send(const void *buf, int count, MPI_Datatype datatype,
             int dest, int tag, MPI_Comm comm);
--- a/src/USER-CGDNA/mf_oxdna.h
+++ b/src/USER-CGDNA/mf_oxdna.h
@ -253,8 +253,8 @@ inline double MFOxdna::DF5(double x, double a, double x_ast,
 }

 /* ----------------------------------------------------------------------
-   test for directionality by projecting base normal n onto delr,
-   returns 1 if nucleotide a to nucleotide b is 3' to 5', otherwise -1
+   test for directionality by projecting base normal n onto delr = a - b,
+   returns 1 if nucleotide b to nucleotide a is 3' to 5', otherwise -1
   ------------------------------------------------------------------------- */
 inline double MFOxdna::is_3pto5p(const double * delr, const double * n)
 {
--- a/src/USER-DPD/fix_eos_cv.cpp
+++ b/src/USER-DPD/fix_eos_cv.cpp
@ -34,7 +34,6 @@ FixEOScv::FixEOScv(LAMMPS *lmp, int narg, char **arg) :
  cvEOS = force->numeric(FLERR,arg[3]);
  if(cvEOS <= 0.0) error->all(FLERR,"EOS cv must be > 0.0");

-  restart_peratom = 1;
  nevery = 1;

  if (atom->dpd_flag != 1)
--- a/src/USER-DPD/fix_eos_table.cpp
+++ b/src/USER-DPD/fix_eos_table.cpp
@ -34,7 +34,6 @@ FixEOStable::FixEOStable(LAMMPS *lmp, int narg, char **arg) :
  Fix(lmp, narg, arg), ntables(0), tables(NULL)
 {
  if (narg != 7) error->all(FLERR,"Illegal fix eos/table command");
-  restart_peratom = 1;
  nevery = 1;

  if (strcmp(arg[3],"linear") == 0) tabstyle = LINEAR;
--- a/src/USER-DPD/fix_eos_table_rx.cpp
+++ b/src/USER-DPD/fix_eos_table_rx.cpp
@ -44,7 +44,6 @@ FixEOStableRX::FixEOStableRX(LAMMPS *lmp, int narg, char **arg) :
  tables2(NULL), dHf(NULL), eosSpecies(NULL)
 {
  if (narg != 8 && narg != 10) error->all(FLERR,"Illegal fix eos/table/rx command");
-  restart_peratom = 1;
  nevery = 1;

  rx_flag = false;
--- a/src/USER-MISC/README
+++ b/src/USER-MISC/README
@ -38,6 +38,7 @@ dihedral_style spherical, Andrew Jewett, jewett.aij@gmail.com, 15 Jul 16
 dihedral_style table, Andrew Jewett, jewett.aij@gmail.com, 10 Jan 12
 fix addtorque, Laurent Joly (U Lyon), ljoly.ulyon at gmail.com, 8 Aug 11
 fix ave/correlate/long, Jorge Ramirez (UPM Madrid), jorge.ramirez at upm.es, 21 Oct 2015
+fix filter/corotate, Lukas Fath (KIT), lukas.fath at kit.edu, 15 Mar 2017
 fix flow/gauss, Joel Eaves (CU Boulder), Joel.Eaves@Colorado.edu, 23 Aug 2016
 fix gle, Michele Ceriotti (EPFL Lausanne), michele.ceriotti at gmail.com, 24 Nov 2014
 fix grem, David Stelter, dstelter@bu.edu, 22 Nov 16
--- a/src/USER-MISC/fix_filter_corotate.cpp
+++ b/src/USER-MISC/fix_filter_corotate.cpp
--- a/src/USER-MISC/fix_filter_corotate.h
+++ b/src/USER-MISC/fix_filter_corotate.h
@ -0,0 +1,139 @@
+/* ----------------------------------------------------------------------
+ LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+ http://lammps.sandia.gov, Sandia National Laboratories
+ Steve Plimpton, sjplimp@sandia.gov
+
+ Copyright (2003) Sandia Corporation.  Under the terms of Contract
+ DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+ certain rights in this software.  This software is distributed under
+ the GNU General Public License.
+
+ See the README file in the top-level LAMMPS directory.
+ ------------------------------------------------------------------------- */
+
+/* ----------------------------------------------------------------------
+    Contributing author: Lukas Fath (KIT)
+    some subroutines are from fix_shake.cpp
+  ------------------------------------------------------------------------- */
+
+#ifdef FIX_CLASS
+
+FixStyle(filter/corotate,FixFilterCorotate)
+
+#else
+
+#ifndef LMP_FIX_FILTER_COROTATE_H
+#define LMP_FIX_FILTER_COROTATE_H
+
+#include "fix.h"
+
+namespace LAMMPS_NS
+{
+
+  class FixFilterCorotate : public Fix
+  {
+  public:
+
+    FixFilterCorotate(class LAMMPS *, int, char **);
+    ~FixFilterCorotate();
+		void setup(int);
+    void setup_pre_neighbor();
+    void pre_neighbor();
+    void setup_pre_force_respa(int,int);
+//    void setup_post_force_respa(int,int);
+    void pre_force_respa(int, int, int);
+    void post_force_respa(int, int, int);
+
+    void init();
+    int setmask();
+
+    double compute_array(int,int);
+
+    int pack_forward_comm(int, int *, double *, int, int *);
+    void unpack_forward_comm(int, int, double *);
+    void grow_arrays(int );
+    double memory_usage();
+
+    void copy_arrays(int, int, int);
+    void set_arrays(int);
+    void update_arrays(int, int);
+
+    int pack_exchange(int, double *);
+    int unpack_exchange(int, double *);
+
+  protected:
+
+    int me,nprocs;
+
+    int flevel;          //filtered respa level
+
+    double **help2;      //temp derivative
+    double **x_store;    //temp for atom->x
+    double *g;             //temp for derivative
+
+    double*n1,*n2,*n3, *del1, *del2,*del3;
+
+    double**dn1dx,**dn2dx,**dn3dx;
+
+    int *bond_flag,*angle_flag;            // bond/angle types to constrain
+    int *type_flag;                        // constrain bonds to these types
+    double *mass_list;                     // constrain bonds to these masses
+    int nmass;                             // # of masses in mass_list
+
+    int molecular;                         // copy of atom->molecular
+    double *bond_distance,*angle_distance; // constraint distances
+
+    int nlevels_respa;                     // copies of needed rRESPA variables
+
+    double **x,**v,**f;                    // local ptrs to atom class quantities
+    double *mass,*rmass;
+    int *type;
+    int nlocal;
+    // atom-based arrays
+    int *shake_flag;                       // 0 if atom not in SHAKE cluster
+    // 1 = size 3 angle cluster
+    // 2,3,4 = size of bond-only cluster
+    tagint **shake_atom;                   // global IDs of atoms in cluster
+    // central atom is 1st
+    // lowest global ID is 1st for size 2
+    int **shake_type;                      // bondtype of each bond in cluster
+    // for angle cluster, 3rd value
+    //   is angletype
+    int *nshake;                           // count
+
+    int *list;                            // list of clusters to SHAKE
+    int nlist,maxlist;                    // size and max-size of list
+    double ***clist_derv;      //stores derivative
+    double **clist_q0;        //stores reference config
+    int *clist_nselect1, *clist_nselect2; //stores length of each selec. list
+    int **clist_select1, **clist_select2; //stores selection lists
+
+    void find_clusters();
+    int masscheck(double);
+
+    void filter_inner();
+    void filter_outer();
+
+    void general_cluster(int,int);
+
+    void stats();
+    int bondtype_findset(int, tagint, tagint, int);
+    int angletype_findset(int, tagint, tagint, int);
+
+    // static variable for ring communication callback to access class data
+    // callback functions for ring communication
+
+    static FixFilterCorotate *fsptr;
+    static void ring_bonds(int, char *);
+    static void ring_nshake(int, char *);
+    static void ring_shake(int, char *);
+
+    int sgn(double val) {
+      return (0 < val) - (val < 0);
+    };
+  };
+
+}
+
+#endif
+#endif
--- a/src/fix_store_state.cpp
+++ b/src/fix_store_state.cpp
@ -38,7 +38,8 @@ enum{KEYWORD,COMPUTE,FIX,VARIABLE,DNAME,INAME};

 FixStoreState::FixStoreState(LAMMPS *lmp, int narg, char **arg) :
  Fix(lmp, narg, arg),
-  nvalues(0), which(NULL), argindex(NULL), value2index(NULL), ids(NULL), values(NULL),
+  nvalues(0), which(NULL), argindex(NULL), value2index(NULL), 
+  ids(NULL), values(NULL),
  vbuf(NULL), pack_choice(NULL)
 {
  if (narg < 5) error->all(FLERR,"Illegal fix store/state command");
@ -1046,7 +1047,8 @@ void FixStoreState::pack_ysu(int n)

  for (int i = 0; i < nlocal; i++) {
    if (mask[i] & groupbit) 
-      vbuf[n] = (x[i][1]-boxylo)*invyprd + (image[i] >> IMGBITS & IMGMASK) - IMGMAX;
+      vbuf[n] = (x[i][1]-boxylo)*invyprd + 
+        (image[i] >> IMGBITS & IMGMASK) - IMGMAX;
    else vbuf[n] = 0.0;
    n += nvalues;
  }
@ -1158,7 +1160,8 @@ void FixStoreState::pack_iy(int n)
  int nlocal = atom->nlocal;

  for (int i = 0; i < nlocal; i++) {
-    if (mask[i] & groupbit) vbuf[n] = (image[i] >> IMGBITS & IMGMASK) - IMGMAX;
+    if (mask[i] & groupbit) 
+      vbuf[n] = (image[i] >> IMGBITS & IMGMASK) - IMGMAX;
    else vbuf[n] = 0.0;
    n += nvalues;
  }
--- a/src/info.cpp
+++ b/src/info.cpp
@ -72,32 +72,33 @@ enum {COMPUTES=1<<0,
      REGIONS=1<<4,
      CONFIG=1<<5,
      TIME=1<<6,
-      VARIABLES=1<<7,
-      SYSTEM=1<<8,
-      COMM=1<<9,
-      ATOM_STYLES=1<<10,
-      INTEGRATE_STYLES=1<<11,
-      MINIMIZE_STYLES=1<<12,
-      PAIR_STYLES=1<<13,
-      BOND_STYLES=1<<14,
-      ANGLE_STYLES=1<<15,
-      DIHEDRAL_STYLES=1<<16,
-      IMPROPER_STYLES=1<<17,
-      KSPACE_STYLES=1<<18,
-      FIX_STYLES=1<<19,
-      COMPUTE_STYLES=1<<20,
-      REGION_STYLES=1<<21,
-      DUMP_STYLES=1<<22,
-      COMMAND_STYLES=1<<23,
+      MEMORY=1<<7,
+      VARIABLES=1<<8,
+      SYSTEM=1<<9,
+      COMM=1<<10,
+      ATOM_STYLES=1<<11,
+      INTEGRATE_STYLES=1<<12,
+      MINIMIZE_STYLES=1<<13,
+      PAIR_STYLES=1<<14,
+      BOND_STYLES=1<<15,
+      ANGLE_STYLES=1<<16,
+      DIHEDRAL_STYLES=1<<17,
+      IMPROPER_STYLES=1<<18,
+      KSPACE_STYLES=1<<19,
+      FIX_STYLES=1<<20,
+      COMPUTE_STYLES=1<<21,
+      REGION_STYLES=1<<22,
+      DUMP_STYLES=1<<23,
+      COMMAND_STYLES=1<<24,
      ALL=~0};

-static const int STYLES = ATOM_STYLES | INTEGRATE_STYLES | MINIMIZE_STYLES | PAIR_STYLES | BOND_STYLES | \
-                         ANGLE_STYLES | DIHEDRAL_STYLES | IMPROPER_STYLES | KSPACE_STYLES | FIX_STYLES | \
-                         COMPUTE_STYLES | REGION_STYLES | DUMP_STYLES | COMMAND_STYLES;
-
+static const int STYLES = ATOM_STYLES | INTEGRATE_STYLES | MINIMIZE_STYLES
+                        | PAIR_STYLES | BOND_STYLES | ANGLE_STYLES
+                        | DIHEDRAL_STYLES | IMPROPER_STYLES | KSPACE_STYLES
+                        | FIX_STYLES | COMPUTE_STYLES | REGION_STYLES
+                        | DUMP_STYLES | COMMAND_STYLES;
 }

-
 static const char *varstyles[] = {
  "index", "loop", "world", "universe", "uloop", "string", "getenv",
  "file", "atomfile", "format", "equal", "atom", "python", "(unknown)"};
@ -174,6 +175,9 @@ void Info::command(int narg, char **arg)
    } else if (strncmp(arg[idx],"time",3) == 0) {
      flags |= TIME;
      ++idx;
+    } else if (strncmp(arg[idx],"memory",3) == 0) {
+      flags |= MEMORY;
+      ++idx;
    } else if (strncmp(arg[idx],"variables",3) == 0) {
      flags |= VARIABLES;
      ++idx;
@ -293,27 +297,42 @@ void Info::command(int narg, char **arg)
    fprintf(out,"\nOS information: %s %s on %s\n",
            ut.sysname, ut.release, ut.machine);
 #endif
+  }
+  
+  if (flags & MEMORY) {

-    fprintf(out,"\nMemory allocation information (MPI rank 0)\n");
+    fprintf(out,"\nMemory allocation information (MPI rank 0):\n\n");
+
+    bigint bytes = 0;
+    bytes += atom->memory_usage();
+    bytes += neighbor->memory_usage();
+    bytes += comm->memory_usage();
+    bytes += update->memory_usage();
+    bytes += force->memory_usage();
+    bytes += modify->memory_usage();
+    for (int i = 0; i < output->ndump; i++)
+      bytes += output->dump[i]->memory_usage();
+    double mbytes = bytes/1024.0/1024.0;
+    fprintf(out,"Total dynamically allocated memory: %.4g Mbyte\n",mbytes);

 #if defined(_WIN32)
    HANDLE phandle = GetCurrentProcess();
    PROCESS_MEMORY_COUNTERS_EX pmc;
    GetProcessMemoryInfo(phandle,(PROCESS_MEMORY_COUNTERS *)&pmc,sizeof(pmc));
-    fprintf(out,"Non-shared memory use: %.3g Mbyte\n",
+    fprintf(out,"Non-shared memory use: %.4g Mbyte\n",
            (double)pmc.PrivateUsage/1048576.0);
-    fprintf(out,"Maximum working set size: %.3g Mbyte\n",
+    fprintf(out,"Maximum working set size: %.4g Mbyte\n",
            (double)pmc.PeakWorkingSetSize/1048576.0);
 #else
 #if defined(__linux)
    struct mallinfo mi;
    mi = mallinfo();
-    fprintf(out,"Total dynamically allocated memory: %.3g Mbyte\n",
-            (double)mi.uordblks/1048576.0);
+    fprintf(out,"Current reserved memory pool size: %.4g Mbyte\n",
+            (double)mi.uordblks/1048576.0+(double)mi.hblkhd/1048576.0);
 #endif
    struct rusage ru;
    if (getrusage(RUSAGE_SELF, &ru) == 0) {
-      fprintf(out,"Maximum resident set size: %.3g Mbyte\n",
+      fprintf(out,"Maximum resident set size: %.4g Mbyte\n",
              (double)ru.ru_maxrss/1024.0);
    }
 #endif
--- a/src/irregular.cpp
+++ b/src/irregular.cpp
@ -395,7 +395,9 @@ int Irregular::create_atom(int n, int *sizes, int *proclist, int sortflag)

  sendmax_proc = 0;
  for (i = 0; i < nsend_proc; i++) {
-    MPI_Send(&length_send[i],1,MPI_INT,proc_send[i],0,world);
+    MPI_Request tmpReq; // Use non-blocking send to avoid possible deadlock
+    MPI_Isend(&length_send[i],1,MPI_INT,proc_send[i],0,world,&tmpReq);
+    MPI_Request_free(&tmpReq); // the MPI_Barrier below marks completion
    sendmax_proc = MAX(sendmax_proc,length_send[i]);
  }

@ -641,7 +643,9 @@ int Irregular::create_data(int n, int *proclist, int sortflag)

  sendmax_proc = 0;
  for (i = 0; i < nsend_proc; i++) {
-    MPI_Send(&num_send[i],1,MPI_INT,proc_send[i],0,world);
+    MPI_Request tmpReq; // Use non-blocking send to avoid possible deadlock
+    MPI_Isend(&num_send[i],1,MPI_INT,proc_send[i],0,world,&tmpReq);
+    MPI_Request_free(&tmpReq); // the MPI_Barrier below marks completion
    sendmax_proc = MAX(sendmax_proc,num_send[i]);
  }

--- a/src/library.cpp
+++ b/src/library.cpp
@ -971,11 +971,9 @@ void lammps_create_atoms(void *ptr, int n, tagint *id, int *type,
      xdata[0] = x[3*i];
      xdata[1] = x[3*i+1];
      xdata[2] = x[3*i+2];
-      if (image) {
-        if (!domain->ownatom(id[i],xdata,&image[i],shrinkexceed)) continue;
-      } else {
-        if (!domain->ownatom(id[i],xdata,NULL,shrinkexceed)) continue;
-      }
+      imageint * img = image ? &image[i] : NULL;
+      tagint     tag = id    ? id[i]     : -1;
+      if (!domain->ownatom(tag, xdata, img, shrinkexceed)) continue;
  
      atom->avec->create_atom(type[i],xdata);
      if (id) atom->tag[nlocal] = id[i];
--- a/src/output.cpp
+++ b/src/output.cpp
@ -812,9 +812,9 @@ void Output::create_restart(int narg, char **arg)
   sum and print memory usage
   result is only memory on proc 0, not averaged across procs
 ------------------------------------------------------------------------- */
-
 void Output::memory_usage()
 {
+
  bigint bytes = 0;
  bytes += atom->memory_usage();
  bytes += neighbor->memory_usage();
@ -825,11 +825,18 @@ void Output::memory_usage()
  for (int i = 0; i < ndump; i++) bytes += dump[i]->memory_usage();

  double mbytes = bytes/1024.0/1024.0;
+  double mbavg,mbmin,mbmax;
+  MPI_Reduce(&mbytes,&mbavg,1,MPI_DOUBLE,MPI_SUM,0,world);
+  MPI_Reduce(&mbytes,&mbmin,1,MPI_DOUBLE,MPI_MIN,0,world);
+  MPI_Reduce(&mbytes,&mbmax,1,MPI_DOUBLE,MPI_MAX,0,world);
+  mbavg /= comm->nprocs;

  if (comm->me == 0) {
    if (screen)
-      fprintf(screen,"Memory usage per processor = %g Mbytes\n",mbytes);
+      fprintf(screen,"Per MPI rank memory allocation (min/avg/max) = "
+              "%.4g | %.4g | %.4g Mbytes\n",mbmin,mbavg,mbmax);
    if (logfile)
-      fprintf(logfile,"Memory usage per processor = %g Mbytes\n",mbytes);
+      fprintf(logfile,"Per MPI rank memory allocation (min/avg/max) = "
+              "%.4g | %.4g | %.4g Mbytes\n",mbmin,mbavg,mbmax);
  }
 }
--- a/src/rcb.cpp
+++ b/src/rcb.cpp
@ -1176,7 +1176,6 @@ void median_merge(void *in, void *inout, int *len, MPI_Datatype *dptr)

 void RCB::invert(int sortflag)
 {
-
  // only create Irregular if not previously created
  // allows Irregular to persist for multiple RCB calls by fix balance

--- a/src/version.h
+++ b/src/version.h
@ -1 +1 @@
-#define LAMMPS_VERSION "10 Mar 2017"
+#define LAMMPS_VERSION "17 Mar 2017"
Author	SHA1	Message	Date
Steve Plimpton	4a90bca7a3	patch 17Mar17	2017-03-17 11:29:23 -06:00
sjplimp	9f35b764f8	Merge pull request #424 from timattox/dpd_restart_fix Remove unneeded restart_peratom flags in USER-DPD package	2017-03-17 11:19:04 -06:00
sjplimp	7ca5dce2f5	Merge pull request #423 from timattox/master_bugfix_irregular bugfix for irregular::create_atom and create_data	2017-03-17 11:18:52 -06:00
sjplimp	fcc3b3bd36	Merge pull request #421 from akohlmey/memory-output improve memory usage reporting in info and output class	2017-03-17 11:18:08 -06:00
sjplimp	53a3877c3d	Merge pull request #420 from rbberger/doc_anchor_check Add utility to detect duplicate anchors in docs	2017-03-17 11:17:28 -06:00
sjplimp	a936b7b2ab	Merge pull request #419 from rbberger/python_fixes Library interface fixes, Python example fixes and better Python 3 compatibility	2017-03-17 11:16:58 -06:00
sjplimp	a91b851f3d	Merge pull request #416 from lukin17/filter_corotate Added fix filter/corotate.	2017-03-17 11:14:59 -06:00
Stan Moore	d31c591b60	Remove unneeded restart_peratom flags	2017-03-17 09:40:39 -06:00
Axel Kohlmeyer	ae5ebf6001	add support for MPI_Request_free() to MPI STUBS library	2017-03-17 11:40:09 -04:00
Axel Kohlmeyer	7fb741d53d	Revert "In irregular.cpp use simpler and slightly faster MPI_Reduce_scatter_block()" This reverts commit `8e75616c14`.	2017-03-17 11:35:59 -04:00
Tim Mattox	8e75616c14	In irregular.cpp use simpler and slightly faster MPI_Reduce_scatter_block()	2017-03-17 03:02:42 -04:00
Tim Mattox	411c069ba6	BUGFIX: Prevent possible deadlock in Irregular::create_atom and create_data	2017-03-17 03:02:19 -04:00
Axel Kohlmeyer	ac82d041cc	ignore package pair style	2017-03-16 23:12:49 -04:00
Richard Berger	621d7d5ce0	Correct off-by-one error in line number reported	2017-03-16 23:05:03 -04:00
Richard Berger	1bb9c7da42	Remove some duplicate anchors in documentation	2017-03-16 22:36:13 -04:00
Richard Berger	f893104b18	Add anchor_check to doc Makefile	2017-03-16 22:21:55 -04:00
Richard Berger	efb2a942e0	Add utility to detect duplicate anchors in documentation files	2017-03-16 22:21:12 -04:00
Axel Kohlmeyer	070ce33a13	improve memory usage reporting in info and output class	2017-03-16 18:35:04 -04:00
Axel Kohlmeyer	f604f86cfc	add fix filter/corotate to some administrative files	2017-03-16 15:08:17 -04:00
Axel Kohlmeyer	bed288339e	simplify and shorten examples for fix filter/corotate and combine into a single folder	2017-03-16 15:02:50 -04:00
Axel Kohlmeyer	1995f434f3	fix some more code formatting issues, add newline at EOF	2017-03-16 14:24:28 -04:00
Axel Kohlmeyer	db0281b4df	Merge branch 'filter_corotate' of https://github.com/lukin17/lammps into pull-416	2017-03-16 12:14:09 -04:00
lukin17	2f5e711acd	Merge remote-tracking branch 'upstream/master' into filter_corotate	2017-03-16 10:00:12 +01:00
lukin17	bdb7669e27	Fixed coding style.	2017-03-16 09:44:07 +01:00
Richard Berger	cda8213892	Added Python matplotlib plot example	2017-03-16 01:41:28 -04:00
Richard Berger	ef940d226c	Improve Python 3 compatibility of pizza tools and simplify read_snapshot code	2017-03-16 01:38:05 -04:00
Richard Berger	36da9223ec	Fix dump cfg in vizplotgui_atomeye.py example	2017-03-15 22:55:51 -04:00
Richard Berger	eb29ef32b1	Fix space/tab error in pizza/gl.py	2017-03-15 22:43:00 -04:00
Richard Berger	29550d472d	Fix dump cfg in viz_atomeye.py example	2017-03-15 22:31:14 -04:00
Richard Berger	79cae51156	Document property 'uses_exceptions' of Python interface	2017-03-15 22:20:30 -04:00
Richard Berger	a210867025	Fixes lammps_create_atoms library function and its Python interface variant The interface of that function has changed and includes two additional parameters, which haven't been added to the Python interface either. This showed up by trying to run the simple.py example.	2017-03-15 22:13:06 -04:00
Richard Berger	0262a54ecf	Fix Python 3 compatibility by encoding strings passed as c_char_p	2017-03-15 22:00:43 -04:00
Axel Kohlmeyer	0d8f74f0c5	Merge branch 'filter_corotate' of https://github.com/lukin17/lammps into pull-416	2017-03-15 18:54:41 -04:00
sjplimp	3a2da51a82	Merge pull request #413 from ohenrich/user-cgdna User cgdna	2017-03-15 13:12:43 -06:00
sjplimp	b1c59126f7	Merge pull request #415 from stanmoore1/kk_qeq Add neigh/qeq option to Kokkos package	2017-03-15 13:12:08 -06:00
sjplimp	4c77838514	Merge pull request #414 from sstrong99/flow-gauss-doc-addition flow/gauss documentation update	2017-03-15 13:11:26 -06:00
sjplimp	f9468f46f5	Merge pull request #412 from timattox/master_typofix Correct a typo in the fix_halt.txt documentation.	2017-03-15 13:10:58 -06:00
lukin17	c3ce3747e0	Added fix filter/corotate.	2017-03-15 11:34:01 +01:00
Stan Moore	fdc390ad05	Tweaking docs for Kokkos package	2017-03-14 14:08:14 -06:00
Stan Moore	580f6b567b	Add neigh/qeq option to Kokkos	2017-03-14 10:44:31 -06:00
Steven Strong	27b1c33a16	updated the NEMD discussion in the how-to documentation about flow/gauss	2017-03-14 10:39:06 -06:00
Oliver Henrich	7a75cd111c	Minor updates in documentation and setup tool, merge before upgrade to oxDNA2	2017-03-14 11:39:09 +00:00
Oliver Henrich	23b8287933	Updated documentation and simple setup tool	2017-03-14 11:36:44 +00:00
Tim Mattox	4cfe623bc1	Correct a typo in the fix_halt.txt documentation.	2017-03-10 21:30:03 -05:00
Oliver Henrich	3626496c7c	Corrected comment in 3' to 5' directionality check	2017-02-22 20:06:49 +00:00
Oliver Henrich	458b6749e7	Corrected comment in 3' to 5' directionality check.	2017-02-22 20:03:41 +00:00
Oliver Henrich	0efd209480	Merge branch 'master' into user-cgdna	2017-02-09 11:50:03 +00:00
Oliver Henrich	ed0a347fbf	Merge branch 'master' into user-cgdna	2017-01-30 10:31:50 +00:00
Oliver Henrich	149f37e764	Corrected reference to Fig.1	2017-01-26 19:08:59 +00:00