new library functions

Merge pull request #239 from akohlmey/static-analysis-fixes
Static analysis fixes
2016-10-27 09:34:04 -06:00 · 2016-10-25 10:32:55 -06:00 · 2016-10-25 10:32:41 -06:00 · 2016-10-25 10:30:34 -06:00 · 2016-10-25 12:12:48 -04:00 · 2016-10-25 11:58:42 -04:00
102 changed files with 1850 additions and 464 deletions
--- a/doc/.gitignore
+++ b/doc/.gitignore
@ -1 +1,5 @@
 /html
+/LAMMPS.epub
+/LAMMPS.mobi
+/Manual.pdf
+/Developer.pdf
--- 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 html pdf old venv
+.PHONY: help clean-all clean epub html pdf old venv

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

@ -32,6 +32,7 @@ help:
 	@echo "  pdf        create Manual.pdf and Developer.pdf in this dir"
 	@echo "  old        create old-style HTML doc pages in old dir"
 	@echo "  fetch      fetch HTML and PDF files from LAMMPS web site"
+	@echo "  epub       create ePUB format manual for e-book readers"
 	@echo "  clean      remove all intermediate RST files"
 	@echo "  clean-all  reset the entire build environment"
 	@echo "  txt2html   build txt2html tool"
@ -63,6 +64,20 @@ html: $(OBJECTS)
 	@rm -rf html/USER/*/*.[sg]*
 	@echo "Build finished. The HTML pages are in doc/html."

+epub: $(OBJECTS)
+	@mkdir -p epub
+	@rm -f LAMMPS.epub
+	@cp src/JPG/lammps-logo.png epub/
+	@(\
+		. $(VENV)/bin/activate ;\
+		cp -r src/* $(RSTDIR)/ ;\
+		sphinx-build -j 8 -b epub -c utils/sphinx-config -d $(BUILDDIR)/doctrees $(RSTDIR) epub ;\
+		deactivate ;\
+	)
+	@mv  epub/LAMMPS.epub .
+	@rm -rf epub
+	@echo "Build finished. The ePUB manual file is created."
+
 pdf: utils/txt2html/txt2html.exe
 	@(\
 		cd src; \
--- a/doc/README
+++ b/doc/README
@ -1,13 +1,14 @@
 LAMMPS Documentation

 Depending on how you obtained LAMMPS, this directory has 2 or 3
-sub-directories and optionally 2 PDF files:
+sub-directories and optionally 2 PDF files and an ePUB file:

 src             content files for LAMMPS documentation
 html            HTML version of the LAMMPS manual (see html/Manual.html)
 tools           tools and settings for building the documentation
 Manual.pdf      large PDF version of entire manual
 Developer.pdf   small PDF with info about how LAMMPS is structured
+LAMMPS.epub     Manual in ePUB format

 If you downloaded LAMMPS as a tarball from the web site, all these
 directories and files should be included.
@ -49,6 +50,7 @@ make pdf          # generate 2 PDF files (Manual.pdf,Developer.pdf)
 make old          # generate old-style HTML pages in old dir via txt2html
 make fetch        # fetch HTML doc pages and 2 PDF files from web site
                  #   as a tarball and unpack into html dir and 2 PDFs
+make epub         # generate LAMMPS.epub in ePUB format using Sphinx 
 make clean        # remove intermediate RST files created by HTML build
 make clean-all    # remove entire build folder and any cached data

@ -91,3 +93,23 @@ This will install virtualenv from the Python Package Index.
 ----------------

 Installing prerequisites for PDF build
+
+[TBA]
+
+----------------
+
+Installing prerequisites for epub build
+
+## ePUB
+
+Same as for HTML. This uses the same tools and configuration
+files as the HTML tree.
+
+For converting the generated ePUB file to a mobi format file
+(for e-book readers like Kindle, that cannot read ePUB), you
+also need to have the 'ebook-convert' tool from the "calibre"
+software installed. http://calibre-ebook.com/
+You first create the ePUB file with 'make epub' and then do:
+
+ebook-convert LAMMPS.epub LAMMPS.mobi
+
--- a/doc/src/JPG/lammps-logo.png
+++ b/doc/src/JPG/lammps-logo.png
--- 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="18 Oct 2016 version">
+<META NAME="docnumber" CONTENT="27 Oct 2016 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
-18 Oct 2016 version :c,h4
+27 Oct 2016 version :c,h4

 Version info: :h4

--- a/doc/src/PDF/colvars-refman-lammps.pdf
+++ b/doc/src/PDF/colvars-refman-lammps.pdf
--- a/doc/src/Section_commands.txt
+++ b/doc/src/Section_commands.txt
@ -106,7 +106,7 @@ the $.  Thus $\{myTemp\} and $x refer to variable names "myTemp" and
 "x".

 How the variable is converted to a text string depends on what style
-of variable it is; see the "variable"_variable doc page for details.
+of variable it is; see the "variable"_variable.html doc page for details.
 It can be a variable that stores multiple text strings, and return one
 of them.  The returned text string can be multiple "words" (space
 separated) which will then be interpreted as multiple arguments in the
@ -528,6 +528,8 @@ These are additional commands in USER packages, which can be used if
 package"_Section_start.html#start_3.

 "dump custom/vtk"_dump_custom_vtk.html,
+"dump nc"_dump_nc.html,
+"dump nc/mpiio"_dump_nc.html,
 "group2ndx"_group2ndx.html,
 "ndx2group"_group2ndx.html :tb(c=3,ea=c)

--- a/doc/src/Section_errors.txt
+++ b/doc/src/Section_errors.txt
@ -8116,11 +8116,11 @@ boundary of a processor's sub-domain has moved more than 1/2 the
 rebuilt and atoms being migrated to new processors.  This also means
 you may be missing pairwise interactions that need to be computed.
 The solution is to change the re-neighboring criteria via the
-"neigh_modify"_neigh_modify command.  The safest settings are "delay 0
-every 1 check yes".  Second, it may mean that an atom has moved far
-outside a processor's sub-domain or even the entire simulation box.
-This indicates bad physics, e.g. due to highly overlapping atoms, too
-large a timestep, etc. :dd
+"neigh_modify"_neigh_modify.html command.  The safest settings are
+"delay 0 every 1 check yes".  Second, it may mean that an atom has
+moved far outside a processor's sub-domain or even the entire
+simulation box. This indicates bad physics, e.g. due to highly
+overlapping atoms, too large a timestep, etc. :dd

 {Out of range atoms - cannot compute PPPM} :dt

@ -8132,11 +8132,11 @@ boundary of a processor's sub-domain has moved more than 1/2 the
 rebuilt and atoms being migrated to new processors.  This also means
 you may be missing pairwise interactions that need to be computed.
 The solution is to change the re-neighboring criteria via the
-"neigh_modify"_neigh_modify command.  The safest settings are "delay 0
-every 1 check yes".  Second, it may mean that an atom has moved far
-outside a processor's sub-domain or even the entire simulation box.
-This indicates bad physics, e.g. due to highly overlapping atoms, too
-large a timestep, etc. :dd
+"neigh_modify"_neigh_modify.html command.  The safest settings are
+"delay 0 every 1 check yes".  Second, it may mean that an atom has
+moved far outside a processor's sub-domain or even the entire
+simulation box. This indicates bad physics, e.g. due to highly
+overlapping atoms, too large a timestep, etc. :dd

 {Out of range atoms - cannot compute PPPMDisp} :dt

@ -8148,11 +8148,11 @@ boundary of a processor's sub-domain has moved more than 1/2 the
 rebuilt and atoms being migrated to new processors.  This also means
 you may be missing pairwise interactions that need to be computed.
 The solution is to change the re-neighboring criteria via the
-"neigh_modify"_neigh_modify command.  The safest settings are "delay 0
-every 1 check yes".  Second, it may mean that an atom has moved far
-outside a processor's sub-domain or even the entire simulation box.
-This indicates bad physics, e.g. due to highly overlapping atoms, too
-large a timestep, etc. :dd
+"neigh_modify"_neigh_modify.html command.  The safest settings are
+"delay 0 every 1 check yes".  Second, it may mean that an atom has
+moved far outside a processor's sub-domain or even the entire
+simulation box. This indicates bad physics, e.g. due to highly
+overlapping atoms, too large a timestep, etc. :dd

 {Overflow of allocated fix vector storage} :dt

--- a/doc/src/Section_howto.txt
+++ b/doc/src/Section_howto.txt
@ -1854,13 +1854,19 @@ internal LAMMPS operations.  Note that LAMMPS classes are defined
 within a LAMMPS namespace (LAMMPS_NS) if you use them from another C++
 application.

-Library.cpp contains these 5 basic functions:
+Library.cpp contains these functions for creating and destroying an
+instance of LAMMPS and sending it commands to execute.  See the
+documentation in the src/library.cpp file for details:

 void lammps_open(int, char **, MPI_Comm, void **)
+void lammps_open_no_mpi(int, char **, void **)
 void lammps_close(void *)
 int lammps_version(void *)
 void lammps_file(void *, char *)
-char *lammps_command(void *, char *) :pre
+char *lammps_command(void *, char *)
+void lammps_commands_list(void *, int, char **)
+void lammps_commands_string(void *, char *)
+void lammps_free(void *) :pre

 The lammps_open() function is used to initialize LAMMPS, passing in a
 list of strings as if they were "command-line
@ -1880,6 +1886,10 @@ half to the other code and run both codes simultaneously before
 syncing them up periodically.  Or it might instantiate multiple
 instances of LAMMPS to perform different calculations.

+The lammps_open_no_mpi() function is similar except that no MPI
+communicator is passed from the caller.  Instead, MPI_COMM_WORLD is
+used to instantiate LAMMPS, and MPI is initialzed if necessary.
+
 The lammps_close() function is used to shut down an instance of LAMMPS
 and free all its memory.

@ -1891,44 +1901,93 @@ changes to the LAMMPS command syntax between versions. The returned
 LAMMPS version code is an integer (e.g. 2 Sep 2015 results in
 20150902) that grows with every new LAMMPS version.

-The lammps_file() and lammps_command() functions are used to pass a
-file or string to LAMMPS as if it were an input script or single
-command in an input script.  Thus the calling code can read or
-generate a series of LAMMPS commands one line at a time and pass it
-thru the library interface to setup a problem and then run it,
-interleaving the lammps_command() calls with other calls to extract
-information from LAMMPS, perform its own operations, or call another
-code's library.
+The lammps_file(), lammps_command(), lammps_commands_list(), and
+lammps_commands_string() functions are used to pass one or more
+commands to LAMMPS to execute, the same as if they were coming from an
+input script.

-Other useful functions are also included in library.cpp.  For example:
+Via these functions, the calling code can read or generate a series of
+LAMMPS commands one or multiple at a time and pass it thru the library
+interface to setup a problem and then run it in stages.  The caller
+can interleave the command function calls with operations it performs,
+calls to extract information from or set information within LAMMPS, or
+calls to another code's library.
+
+The lammps_file() function passes the filename of an input script.
+The lammps_command() function passes a single command as a string.
+The lammps_commands_list() function passes multiple commands in a
+char** list.  In both lammps_command() and lammps_commands_list(),
+individual commands may or may not have a trailing newline.  The
+lammps_commands_string() function passes multiple commands
+concatenated into one long string, separated by newline characters.
+In both lammps_commands_list() and lammps_commands_string(), a single
+command can be spread across multiple lines, if the last printable
+character of all but the last line is "&", the same as if the lines
+appeared in an input script.
+
+The lammps_free() function is a clean-up function to free memory that
+the library allocated previously via other function calls.  See
+comments in src/library.cpp file for which other functions need this
+clean-up.
+
+Library.cpp also contains these functions for extracting information
+from LAMMPS and setting value within LAMMPS.  Again, see the
+documentation in the src/library.cpp file for details, including
+which quantities can be queried by name:

 void *lammps_extract_global(void *, char *)
 void *lammps_extract_atom(void *, char *)
 void *lammps_extract_compute(void *, char *, int, int)
 void *lammps_extract_fix(void *, char *, int, int, int, int)
-void *lammps_extract_variable(void *, char *, char *)
+void *lammps_extract_variable(void *, char *, char *) :pre
+
 int lammps_set_variable(void *, char *, char *)
+double lammps_get_thermo(void *, char *) :pre
+
 int lammps_get_natoms(void *)
-void lammps_get_coords(void *, double *)
-void lammps_put_coords(void *, double *) :pre
+void lammps_gather_atoms(void *, double *)
+void lammps_scatter_atoms(void *, double *) :pre
+void lammps_create_atoms(void *, int, tagint *, int *, double *, double *) :pre

-These can extract various global or per-atom quantities from LAMMPS as
-well as values calculated by a compute, fix, or variable.  The
-"set_variable" function can set an existing string-style variable to a
-new value, so that subsequent LAMMPS commands can access the variable.
-The "get" and "put" operations can retrieve and reset atom
-coordinates.  See the library.cpp file and its associated header file
-library.h for details.
+The extract functions return a pointer to various global or per-atom
+quantities stored in LAMMPS or to values calculated by a compute, fix,
+or variable.  The pointer returned by the extract_global() function
+can be used as a permanent reference to a value which may change.  For
+the other extract functions, the underlying storage may be reallocated
+as LAMMPS runs, so you need to re-call the function to assure a
+current pointer or returned value(s).

-The key idea of the library interface is that you can write any
-functions you wish to define how your code talks to LAMMPS and add
-them to src/library.cpp and src/library.h, as well as to the "Python
-interface"_Section_python.html.  The routines you add can access or
-change any LAMMPS data you wish.  The examples/COUPLE and python
-directories have example C++ and C and Python codes which show how a
-driver code can link to LAMMPS as a library, run LAMMPS on a subset of
-processors, grab data from LAMMPS, change it, and put it back into
-LAMMPS.
+The lammps_set_variable() function can set an existing string-style
+variable to a new string value, so that subsequent LAMMPS commands can
+access the variable.  The lammps_get_thermo() function returns the
+current value of a thermo keyword as a double.
+
+The lammps_get_natoms() function returns the total number of atoms in
+the system and can be used by the caller to allocate space for the
+lammps_gather_atoms() and lammps_scatter_atoms() functions.  The
+gather function collects atom info of the requested type (atom coords,
+types, forces, etc) from all procsesors, orders them by atom ID, and
+returns a full list to each calling processor.  The scatter function
+does the inverse.  It distributes the same kinds of values, 
+passed by the caller, to each atom owned by individual processors.
+
+The lammps_create_atoms() function takes a list of N atoms as input
+with atom types and coords (required), an optionally atom IDs and
+velocities.  It uses the coords of each atom to assign it as a new
+atom to the processor that owns it.  Additional properties for the new
+atoms can be assigned via the lammps_scatter_atoms() or
+lammps_extract_atom() functions.
+
+The examples/COUPLE and python directories have example C++ and C and
+Python codes which show how a driver code can link to LAMMPS as a
+library, run LAMMPS on a subset of processors, grab data from LAMMPS,
+change it, and put it back into LAMMPS.
+
+NOTE: You can write code for additional functions as needed to define
+how your code talks to LAMMPS and add them to src/library.cpp and
+src/library.h, as well as to the "Python
+interface"_Section_python.html.  The added functions can access or
+change any LAMMPS data you wish.

 :line

--- a/doc/src/Section_packages.txt
+++ b/doc/src/Section_packages.txt
@ -1153,6 +1153,7 @@ Package, Description, Author(s), Doc page, Example, Pic/movie, Library
 "USER-MISC"_#USER-MISC, single-file contributions, USER-MISC/README, USER-MISC/README, -, -, -
 "USER-MANIFOLD"_#USER-MANIFOLD, motion on 2d surface, Stefan Paquay (Eindhoven U of Technology), "fix manifoldforce"_fix_manifoldforce.html, USER/manifold, "manifold"_manifold, -
 "USER-MOLFILE"_#USER-MOLFILE, "VMD"_VMD molfile plug-ins, Axel Kohlmeyer (Temple U), "dump molfile"_dump_molfile.html, -, -, VMD-MOLFILE
+"USER-NC-DUMP"_#USER-NC-DUMP, dump output via NetCDF, Lars Pastewka (Karlsruhe Institute of Technology, KIT), "dump nc, dump nc/mpiio"_dump_nc.html, -, -, lib/netcdf
 "USER-OMP"_#USER-OMP, OpenMP threaded styles, Axel Kohlmeyer (Temple U), "Section 5.3.4"_accelerate_omp.html, -, -, -
 "USER-PHONON"_#USER-PHONON, phonon dynamical matrix, Ling-Ti Kong (Shanghai Jiao Tong U), "fix phonon"_fix_phonon.html, USER/phonon, -, -
 "USER-QMMM"_#USER-QMMM, QM/MM coupling, Axel Kohlmeyer (Temple U), "fix qmmm"_fix_qmmm.html, USER/qmmm, -, lib/qmmm
@ -1598,6 +1599,29 @@ The person who created this package is Axel Kohlmeyer at Temple U

 :line

+USER-NC-DUMP package :link(USER-NC-DUMP),h5
+
+Contents: Dump styles for writing NetCDF format files.  NetCDF is a binary,
+portable, self-describing file format on top of HDF5. The file format
+contents follow the AMBER NetCDF trajectory conventions
+(http://ambermd.org/netcdf/nctraj.xhtml), but include extensions to this
+convention. This package implements a "dump nc"_dump_nc.html command
+and a "dump nc/mpiio"_dump_nc.html command to output LAMMPS snapshots
+in this format.  See src/USER-NC-DUMP/README for more details.
+
+NetCDF files can be directly visualized with the following tools:
+Ovito (http://www.ovito.org/). Ovito supports the AMBER convention
+  and all of the above extensions. :ulb,l
+VMD (http://www.ks.uiuc.edu/Research/vmd/) :l
+AtomEye (http://www.libatoms.org/). The libAtoms version of AtomEye contains
+  a NetCDF reader that is not present in the standard distribution of AtomEye :l,ule
+
+The person who created these files is Lars Pastewka at
+Karlsruhe Institute of Technology (lars.pastewka at kit.edu).
+Contact him directly if you have questions.
+
+:line
+
 USER-OMP package :link(USER-OMP),h5

 Supporting info:
--- a/doc/src/Section_python.txt
+++ b/doc/src/Section_python.txt
@ -534,10 +534,11 @@ from lammps import lammps :pre
 These are the methods defined by the lammps module.  If you look at
 the files src/library.cpp and src/library.h you will see that they
 correspond one-to-one with calls you can make to the LAMMPS library
-from a C++ or C or Fortran program.
+from a C++ or C or Fortran program, and which are described in
+"Section 6.19"_Section_howto.html#howto_19 of the manual.

 lmp = lammps()           # create a LAMMPS object using the default liblammps.so library
-                         4 optional args are allowed: name, cmdargs, ptr, comm
+                         # 4 optional args are allowed: name, cmdargs, ptr, comm
 lmp = lammps(ptr=lmpptr) # use lmpptr as previously created LAMMPS object
 lmp = lammps(comm=split) # create a LAMMPS object with a custom communicator, requires mpi4py 2.0.0 or later
 lmp = lammps(name="g++")   # create a LAMMPS object using the liblammps_g++.so library
@ -549,6 +550,8 @@ version = lmp.version()  # return the numerical version id, e.g. LAMMPS 2 Sep 20

 lmp.file(file)           # run an entire input script, file = "in.lj"
 lmp.command(cmd)         # invoke a single LAMMPS command, cmd = "run 100" :pre
+lmp.commands_list(cmdlist)     # invoke commands in cmdlist = ["run 10", "run 20"]
+lmp.commands_string(multicmd)  # invoke commands in multicmd = "run 10\nrun 20"

 xlo = lmp.extract_global(name,type)  # extract a global quantity
                                     # name = "boxxlo", "nlocal", etc
@ -580,6 +583,8 @@ var = lmp.extract_variable(name,group,flag)  # extract value(s) from a variable
                                             #        1 = atom-style variable :pre

 flag = lmp.set_variable(name,value)       # set existing named string-style variable to value, flag = 0 if successful
+value = lmp.get_thermo(name)              # return current value of a thermo keyword
+
 natoms = lmp.get_natoms()                 # total # of atoms as int
 data = lmp.gather_atoms(name,type,count)  # return atom attribute of all atoms gathered into data, ordered by atom ID
                                          # name = "x", "charge", "type", etc
@ -599,9 +604,10 @@ create an instance of LAMMPS, wrapped in a Python class by the lammps
 Python module, and return an instance of the Python class as lmp.  It
 is used to make all subequent calls to the LAMMPS library.

-Additional arguments can be used to tell Python the name of the shared
-library to load or to pass arguments to the LAMMPS instance, the same
-as if LAMMPS were launched from a command-line prompt.
+Additional arguments to lammps() can be used to tell Python the name
+of the shared library to load or to pass arguments to the LAMMPS
+instance, the same as if LAMMPS were launched from a command-line
+prompt.

 If the ptr argument is set like this:

@ -626,8 +632,9 @@ lmp2 = lammps()
 lmp1.file("in.file1")
 lmp2.file("in.file2") :pre

-The file() and command() methods allow an input script or single
-commands to be invoked.
+The file(), command(), commands_list(), commands_string() methods
+allow an input script, a single command, or multiple commands to be
+invoked.

 The extract_global(), extract_atom(), extract_compute(),
 extract_fix(), and extract_variable() methods return values or
--- a/doc/src/Section_start.txt
+++ b/doc/src/Section_start.txt
@ -1601,9 +1601,9 @@ implementations, either by environment variables that specify how to
 order physical processors, or by config files that specify what
 physical processors to assign to each MPI rank.  The -reorder switch
 simply gives you a portable way to do this without relying on MPI
-itself.  See the "processors out"_processors command for how to output
-info on the final assignment of physical processors to the LAMMPS
-simulation domain.
+itself.  See the "processors out"_processors.html command for how
+to output info on the final assignment of physical processors to
+the LAMMPS simulation domain.

 -screen file :pre

--- a/doc/src/accelerate_intel.txt
+++ b/doc/src/accelerate_intel.txt
@ -151,7 +151,7 @@ can start running so that the CPU pipeline is still being used
 efficiently. Although benefits can be seen by launching a MPI task
 for every hardware thread, for multinode simulations, we recommend
 that OpenMP threads are used for SMT instead, either with the
-USER-INTEL package, "USER-OMP package"_accelerate_omp.html", or
+USER-INTEL package, "USER-OMP package"_accelerate_omp.html, or
 "KOKKOS package"_accelerate_kokkos.html. In the example above, up
 to 36X speedups can be observed by using all 36 physical cores with
 LAMMPS. By using all 72 hardware threads, an additional 10-30%
@ -343,7 +343,7 @@ when using offload.

 Not all styles are supported in the USER-INTEL package. You can mix
 the USER-INTEL package with styles from the "OPT"_accelerate_opt.html
-package or the "USER-OMP package"_accelerate_omp.html". Of course,
+package or the "USER-OMP package"_accelerate_omp.html. Of course,
 this requires that these packages were installed at build time. This
 can performed automatically by using "-sf hybrid intel opt" or
 "-sf hybrid intel omp" command-line options. Alternatively, the "opt"
--- a/doc/src/atom_style.txt
+++ b/doc/src/atom_style.txt
@ -166,7 +166,7 @@ stores a per-particle mass and size and orientation (i.e. the corner
 points of the triangle).

 The {template} style allows molecular topolgy (bonds,angles,etc) to be
-defined via a molecule template using the "molecule"_molecule.txt
+defined via a molecule template using the "molecule"_molecule.html
 command.  The template stores one or more molecules with a single copy
 of the topology info (bonds,angles,etc) of each.  Individual atoms
 only store a template index and template atom to identify which
--- a/doc/src/bond_fene.txt
+++ b/doc/src/bond_fene.txt
@ -73,7 +73,7 @@ This bond style can only be used if LAMMPS was built with the
 MOLECULE package (which it is by default).  See the "Making
 LAMMPS"_Section_start.html#start_3 section for more info on packages.

-You typically should specify "special_bonds fene"_special_bonds.html"
+You typically should specify "special_bonds fene"_special_bonds.html
 or "special_bonds lj/coul 0 1 1"_special_bonds.html to use this bond
 style.  LAMMPS will issue a warning it that's not the case.

--- a/doc/src/bond_fene_expand.txt
+++ b/doc/src/bond_fene_expand.txt
@ -76,7 +76,7 @@ This bond style can only be used if LAMMPS was built with the
 MOLECULE package (which it is by default).  See the "Making
 LAMMPS"_Section_start.html#start_3 section for more info on packages.

-You typically should specify "special_bonds fene"_special_bonds.html"
+You typically should specify "special_bonds fene"_special_bonds.html
 or "special_bonds lj/coul 0 1 1"_special_bonds.html to use this bond
 style.  LAMMPS will issue a warning it that's not the case.

--- a/doc/src/compute_fep.txt
+++ b/doc/src/compute_fep.txt
@ -236,7 +236,7 @@ LAMMPS"_Section_start.html#start_3 section for more info.
 [Related commands:]

 "fix adapt/fep"_fix_adapt_fep.html, "fix ave/time"_fix_ave_time.html,
-"pair_lj_soft_coul_soft"_pair_lj_soft_coul_soft.txt
+"pair_style lj/soft/coul/soft"_pair_lj_soft.html

 [Default:]

--- a/doc/src/compute_orientorder_atom.txt
+++ b/doc/src/compute_orientorder_atom.txt
@ -15,7 +15,7 @@ compute ID group-ID orientorder/atom keyword values ... :pre
 ID, group-ID are documented in "compute"_compute.html command :ulb,l
 orientorder/atom = style name of this compute command :l
 one or more keyword/value pairs may be appended :l
-keyword = {cutoff} or {nnn} or {ql}
+keyword = {cutoff} or {nnn} or {degrees}
  {cutoff} value = distance cutoff
  {nnn} value = number of nearest neighbors
  {degrees} values = nlvalues, l1, l2,...  :pre
@ -111,7 +111,7 @@ options.

 [Default:]

-The option defaults are {cutoff} = pair style cutoff, {nnn} = 12, {degrees} = 5 4 6 8 9 10 12 i.e. {Q}4, {Q}6, {Q}8, {Q}10, and {Q}12.
+The option defaults are {cutoff} = pair style cutoff, {nnn} = 12, {degrees} = 5 4 6 8 10 12 i.e. {Q}4, {Q}6, {Q}8, {Q}10, and {Q}12.

 :line

--- a/doc/src/compute_property_local.txt
+++ b/doc/src/compute_property_local.txt
@ -78,7 +78,7 @@ defined by the "pair_style"_pair_style.html command for the types of
 the two atoms is used.  For the {radius} setting, the sum of the radii
 of the two particles is used as a cutoff.  For example, this is
 appropriate for granular particles which only interact when they are
-overlapping, as computed by "granular pair styles"_pair_gran.txt.
+overlapping, as computed by "granular pair styles"_pair_gran.html.

 If the inputs are bond, angle, etc attributes, the local data is
 generated by looping over all the atoms owned on a processor and
--- a/doc/src/dump.txt
+++ b/doc/src/dump.txt
@ -12,6 +12,7 @@ dump command :h3
 "dump image"_dump_image.html command :h3
 "dump movie"_dump_image.html command :h3
 "dump molfile"_dump_molfile.html command :h3
+"dump nc"_dump_nc.html command :h3

 [Syntax:]

@ -43,7 +44,9 @@ args = list of arguments for a particular style :l

  {movie} args = discussed on "dump image"_dump_image.html doc page :pre

-  {molfile} args = discussed on "dump molfile"_dump_molfile.html doc page :pre
+  {molfile} args = discussed on "dump molfile"_dump_molfile.html doc page
+
+  {nc} args = discussed on "dump nc"_dump_nc.html doc page :pre

  {local} args = list of local attributes
    possible attributes = index, c_ID, c_ID\[I\], f_ID, f_ID\[I\]
--- a/doc/src/dump_nc.txt
+++ b/doc/src/dump_nc.txt
@ -31,30 +31,32 @@ dump 1 all nc/mpiio 1000 traj.nc id type x y z :pre
 [Description:]

 Dump a snapshot of atom coordinates every N timesteps in Amber-style
-NetCDF file format. NetCDF files are binary, portable and self-describing.
-This dump style will write only one file on the root node.  The dump
-style {nc} uses the "standard NetCDF library"_netcdf-home all data is
-collected on one processor and then written to the dump file. Dump style
-{nc/mpiio} used the "parallel NetCDF library"_pnetcdf-home and MPI-IO;
-it has better performance on a larger number of processors. Note that
-'nc' outputs all atoms sorted by atom tag while 'nc/mpiio' outputs in
-order of the MPI rank.
+NetCDF file format. NetCDF files are binary, portable and
+self-describing.  This dump style will write only one file on the root
+node.  The dump style {nc} uses the "standard NetCDF
+library"_netcdf-home all data is collected on one processor and then
+written to the dump file. Dump style {nc/mpiio} used the "parallel
+NetCDF library"_pnetcdf-home and MPI-IO; it has better performance on
+a larger number of processors. Note that 'nc' outputs all atoms sorted
+by atom tag while 'nc/mpiio' outputs in order of the MPI rank.

-In addition to per-atom data, also global (i.e. not per atom, but per frame)
-quantities can be included in the dump file. This can be variables, output
-from computes or fixes data prefixed with v_, c_ and f_, respectively.
-These properties are included via "dump_modify"_dump_modify.html {global}.
+In addition to per-atom data, also global (i.e. not per atom, but per
+frame) quantities can be included in the dump file. This can be
+variables, output from computes or fixes data prefixed with v_, c_ and
+f_, respectively.  These properties are included via
+"dump_modify"_dump_modify.html {global}.
+
+:link(netcdf-home,http://www.unidata.ucar.edu/software/netcdf/)
+:link(pnetcdf-home,http://trac.mcs.anl.gov/projects/parallel-netcdf/)

 :line

 [Restrictions:]

-The {nc} and {nc/mpiio} dump styles are part of the USER-NC-DUMP 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.
-
-:link(netcdf-home,http://www.unidata.ucar.edu/software/netcdf/)
-:link(pnetcdf-home,http://trac.mcs.anl.gov/projects/parallel-netcdf/)
+The {nc} and {nc/mpiio} dump styles are part of the USER-NC-DUMP
+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.

 :line

--- a/doc/src/fix.txt
+++ b/doc/src/fix.txt
@ -190,6 +190,7 @@ of "this page"_Section_commands.html#cmd_5.
 "gcmc"_fix_gcmc.html - grand canonical insertions/deletions
 "gld"_fix_gcmc.html - generalized Langevin dynamics integrator
 "gravity"_fix_gravity.html - add gravity to atoms in a granular simulation
+"halt"_fix_halt.html - terminate a dynamics run or minimization
 "heat"_fix_heat.html - add/subtract momentum-conserving heat
 "indent"_fix_indent.html - impose force due to an indenter
 "langevin"_fix_langevin.html - Langevin temperature control
--- a/doc/src/fix_phonon.txt
+++ b/doc/src/fix_phonon.txt
@ -184,7 +184,7 @@ This fix requires LAMMPS be built with an FFT library.  See the
 [Default:]

 The option defaults are sysdim = the same dimemsion as specified by
-the "dimension"_dimension command, and nasr = 20.
+the "dimension"_dimension.html command, and nasr = 20.

 :line

--- a/doc/src/fix_property_atom.txt
+++ b/doc/src/fix_property_atom.txt
@ -177,7 +177,7 @@ their values.  This means that the values can be output via the "dump
 custom"_dump.html command, accessed by fixes like "fix
 ave/atom"_fix_ave_atom.html, accessed by other computes like "compute
 reduce"_compute_reduce.html, or used in "atom-style
-variables"_variables.
+variables"_variable.html.

 For example, these commands will output two new properties to a custom
 dump file:
--- a/doc/src/fix_rigid.txt
+++ b/doc/src/fix_rigid.txt
@ -620,7 +620,7 @@ rigid styles for the rigid bodies. :l

 Use "fix press/berendsen"_fix_press_berendsen.html to compute the
 pressure and change the box dimensions.  Use one of the 4 NVE or 2 NVT
-rigid styles for the rigid bodies.  Use "fix nvt"_fix_nh.thml (or any
+rigid styles for the rigid bodies.  Use "fix nvt"_fix_nh.html (or any
 other thermostat) for the non-rigid particles. :l
 :ule

--- a/doc/src/fix_ti_spring.txt
+++ b/doc/src/fix_ti_spring.txt
@ -104,7 +104,7 @@ the Nose-Hoover thermostat ("fix nvt"_fix_nh.html) is {NOT}
 recommended due to its well documented issues with the canonical
 sampling of harmonic degrees of freedom (notice that the {chain}
 option will {NOT} solve this problem). The Langevin thermostat ("fix
-langevin"_fix_langevin.html") correctly thermostats the system and we
+langevin"_fix_langevin.html) correctly thermostats the system and we
 advise its usage with ti/spring command.

 [Restart, fix_modify, output, run start/stop, minimize info:]
--- a/doc/src/neb.txt
+++ b/doc/src/neb.txt
@ -83,9 +83,9 @@ replica.  Conceptually, the non-NEB atoms provide a background force
 field for the NEB atoms.  They can be allowed to move during the NEB
 minimiation procedure (which will typically induce different
 coordinates for non-NEB atoms in different replicas), or held fixed
-using other LAMMPS commands such as "fix setforce"_fix_setforce.  Note
-that the "partition"_partition.html command can be used to invoke a
-command on a subset of the replicas, e.g. if you wish to hold NEB or
+using other LAMMPS commands such as "fix setforce"_fix_setforce.html.
+Note that the "partition"_partition.html command can be used to invoke
+a command on a subset of the replicas, e.g. if you wish to hold NEB or
 non-NEB atoms fixed in only the end-point replicas.

 The initial atomic configuration for each of the replicas can be
--- a/doc/src/pair_hbond_dreiding.txt
+++ b/doc/src/pair_hbond_dreiding.txt
@ -138,8 +138,8 @@ angle cutoff (degrees) :ul
 A single hydrogen atom type K can be specified, or a wild-card
 asterisk can be used in place of or in conjunction with the K
 arguments to select multiple types as hydrogens.  This takes the form
-"*" or "*n" or "n*" or "m*n".  See the "pair_coeff"_pair_coeff command
-doc page for details.
+"*" or "*n" or "n*" or "m*n".  See the "pair_coeff"_pair_coeff.html
+command doc page for details.

 If the donor flag is {i}, then the atom of type I in the pair_coeff
 command is treated as the donor, and J is the acceptor.  If the donor
--- a/doc/src/pair_line_lj.txt
+++ b/doc/src/pair_line_lj.txt
@ -60,8 +60,8 @@ pair_style command or overridden with an optional argument in the
 pair_coeff command for a type pair as discussed below.  The distance
 between the centers of 2 line segments, or the center of a line
 segment and a point particle, must be less than this distance (plus
-the neighbor skin; see the "neighbor"_neighbor command), for the pair
-of particles to be included in the neighbor list.
+the neighbor skin; see the "neighbor"_neighbor.html command), for
+the pair of particles to be included in the neighbor list.

 NOTE: This means that a too-short value for the {cutoff} setting can
 exclude a pair of particles from the neighbor list even if pairs of
--- a/doc/src/pair_lubricateU.txt
+++ b/doc/src/pair_lubricateU.txt
@ -119,7 +119,7 @@ of walls (whether moving or stationary) will affect the volume
 fraction available to colloidal particles. This is currently accounted
 for with the following types of walls: "wall/lj93"_fix_wall.html,
 "wall/lj126"_fix_wall.html, "wall/colloid"_fix_wall.html, and
-"wall/harmonic_fix_wall.html".  For these wall styles, the correct
+"wall/harmonic"_fix_wall.html.  For these wall styles, the correct
 volume fraction will be used when walls do not coincide with the box
 boundary, as well as when walls move and thereby cause a change in the
 volume fraction. To use these wall styles with pair_style {lubricateU}
--- a/doc/src/pair_thole.txt
+++ b/doc/src/pair_thole.txt
@ -180,7 +180,7 @@ package.
 langevin/drude"_fix_langevin_drude.html, "fix
 drude/transform"_fix_drude_transform.html, "compute
 temp/drude"_compute_temp_drude.html
-"pair_style lj/cut/coul/long"_pair_lj_cut_coul_long
+"pair_style lj/cut/coul/long"_pair_lj.html

 [Default:] none

--- a/doc/src/read_data.txt
+++ b/doc/src/read_data.txt
@ -431,8 +431,8 @@ Atoms # sphere
 Pair Coeffs # lj/cut :pre

 will check if the currently-defined "atom_style"_atom_style.html is
-{sphere}, and the current "pair_style"_pair_style is {lj/cut}.  If
-not, LAMMPS will issue a warning to indicate that the data file
+{sphere}, and the current "pair_style"_pair_style.html is {lj/cut}.
+If not, LAMMPS will issue a warning to indicate that the data file
 section likely does not contain the correct number or type of
 parameters expected for the currently-defined style.

--- a/doc/src/read_dump.txt
+++ b/doc/src/read_dump.txt
@ -322,3 +322,6 @@ They are only enabled if LAMMPS was built with that packages.  See the

 The option defaults are box = yes, replace = yes, purge = no, trim =
 no, add = no, scaled = no, wrapped = yes, and format = native.
+
+:link(vmd,http://www.ks.uiuc.edu/Research/vmd)
+
--- a/doc/src/read_restart.txt
+++ b/doc/src/read_restart.txt
@ -141,11 +141,11 @@ these settings after the restart file is read.
 "units"_units.html
 "newton bond"_newton.html (see discussion of newton command below)
 "atom style"_atom_style.html and "atom_modify"_atom_modify.html settings id, map, sort
-"comm style"_comm_style.html and "comm_modify"_comm_modify settings mode, cutoff, vel
+"comm style"_comm_style.html and "comm_modify"_comm_modify.html settings mode, cutoff, vel
 "timestep"_timestep.html
 simulation box size and shape and "boundary"_boundary.html settings
 atom "group"_group.html definitions
-per-type atom settings such as "mass"_mass.thml
+per-type atom settings such as "mass"_mass.html
 per-atom attributes including their group assignments and molecular topology attributes (bonds, angles, etc)
 force field styles ("pair"_pair_style.html, "bond"_bond_style.html, "angle"_angle_style.html, etc)
 force field coefficients ("pair"_pair_coeff.html, "bond"_bond_coeff.html, "angle"_angle_coeff.html, etc) in some cases (see below)
--- a/doc/src/thermo_style.txt
+++ b/doc/src/thermo_style.txt
@ -108,7 +108,7 @@ Style {custom} is the most general setting and allows you to specify
 which of the keywords listed above you want printed on each
 thermodynamic timestep.  Note that the keywords c_ID, f_ID, v_name are
 references to "computes"_compute.html, "fixes"_fix.html, and
-equal-style "variables"_variable.html" that have been defined
+equal-style "variables"_variable.html that have been defined
 elsewhere in the input script or can even be new styles which users
 have added to LAMMPS (see the "Section 10"_Section_modify.html
 section of the documentation).  Thus the {custom} style provides a
--- a/doc/src/write_dump.txt
+++ b/doc/src/write_dump.txt
@ -26,7 +26,7 @@ write_dump all atom dump.atom
 write_dump subgroup atom dump.run.bin
 write_dump all custom dump.myforce.* id type x y vx fx
 write_dump flow custom dump.%.myforce id type c_myF\[3\] v_ke modify sort id
-write_dump all xyz system.xyz modify sort id elements O H
+write_dump all xyz system.xyz modify sort id element O H
 write_dump all image snap*.jpg type type size 960 960 modify backcolor white
 write_dump all image snap*.jpg element element &
   bond atom 0.3 shiny 0.1 ssao yes 6345 0.2 size 1600 1600  &
--- a/doc/utils/sphinx-config/conf.py
+++ b/doc/utils/sphinx-config/conf.py
@ -276,4 +276,27 @@ texinfo_documents = [
 # If true, do not generate a @detailmenu in the "Top" node's menu.
 #texinfo_no_detailmenu = False

+# -- Options for ePUB output ----------------------------------------------
+
+epub_title = 'LAMMPS Documentation - ' + get_lammps_version()
+
+epub_cover = ('lammps-logo.png', '')
+
+epub_description = """
+This is the Manual for the LAMMPS software package.
+
+LAMMPS stands for Large-scale Atomic/Molecular Massively Parallel
+Simulator and is a classical molecular dynamics simulation code
+designed to run efficiently on parallel computers.  It was developed
+at Sandia National Laboratories, a US Department of Energy facility,
+with funding from the DOE.  It is an open-source code, distributed
+freely under the terms of the GNU Public License (GPL).
+
+The primary author of the code is Steve Plimpton, who can be emailed
+at sjplimp@sandia.gov.  The LAMMPS WWW Site at lammps.sandia.gov has
+more information about the code and its uses.
+"""
+
+epub_author = 'The LAMMPS Developers'
+

--- a/examples/COUPLE/simple/in.lj
+++ b/examples/COUPLE/simple/in.lj
@ -20,4 +20,6 @@ neigh_modify	delay 0 every 20 check no

 fix		1 all nve

+variable        fx atom fx
+
 run		10
--- a/examples/COUPLE/simple/simple.c
+++ b/examples/COUPLE/simple/simple.c
@ -71,8 +71,8 @@ int main(int narg, char **arg)
     (could just send it to proc 0 of comm_lammps and let it Bcast)
     all LAMMPS procs call lammps_command() on the line */

-  void *ptr;
-  if (lammps == 1) lammps_open(0,NULL,comm_lammps,&ptr);
+  void *lmp = NULL;
+  if (lammps == 1) lammps_open(0,NULL,comm_lammps,&lmp);

  int n;
  char line[1024];
@ -85,7 +85,7 @@ int main(int narg, char **arg)
    MPI_Bcast(&n,1,MPI_INT,0,MPI_COMM_WORLD);
    if (n == 0) break;
    MPI_Bcast(line,n,MPI_CHAR,0,MPI_COMM_WORLD);
-    if (lammps == 1) lammps_command(ptr,line);
+    if (lammps == 1) lammps_command(lmp,line);
  }

  /* run 10 more steps
@ -94,23 +94,72 @@ int main(int narg, char **arg)
     put coords back into LAMMPS
     run a single step with changed coords */

-  if (lammps == 1) {
-    lammps_command(ptr,"run 10");
+  double *x = NULL;
+  double *v = NULL;

-    int natoms = lammps_get_natoms(ptr);
-    double *x = (double *) malloc(3*natoms*sizeof(double));
-    lammps_gather_atoms(ptr,"x",1,3,x);
+  if (lammps == 1) {
+    lammps_command(lmp,"run 10");
+
+    int natoms = lammps_get_natoms(lmp);
+    x = (double *) malloc(3*natoms*sizeof(double));
+    lammps_gather_atoms(lmp,"x",1,3,x);
+    v = (double *) malloc(3*natoms*sizeof(double));
+    lammps_gather_atoms(lmp,"v",1,3,v);
    double epsilon = 0.1;
    x[0] += epsilon;
-    lammps_scatter_atoms(ptr,"x",1,3,x);
-    free(x);
+    lammps_scatter_atoms(lmp,"x",1,3,x);

-    lammps_command(ptr,"run 1");
+    lammps_command(lmp,"run 1");
  }

-  if (lammps == 1) lammps_close(ptr);
+  // extract force on single atom two different ways
+
+  if (lammps == 1) {
+    double **f = (double **) lammps_extract_atom(lmp,"f");
+    printf("Force on 1 atom via extract_atom: %g\n",f[0][0]);
+
+    double *fx = (double *) lammps_extract_variable(lmp,"fx","all");
+    printf("Force on 1 atom via extract_variable: %g\n",fx[0]);
+  }
+
+  /* use commands_string() and commands_list() to invoke more commands */
+
+  char *strtwo = "run 10\nrun 20";
+  if (lammps == 1) lammps_commands_string(lmp,strtwo);
+
+  char *cmds[2];
+  cmds[0] = "run 10";
+  cmds[1] = "run 20";
+  if (lammps == 1) lammps_commands_list(lmp,2,cmds);
+
+  /* delete all atoms
+     create_atoms() to create new ones with old coords, vels
+     initial thermo should be same as step 20 */
+
+  int *type = NULL;
+
+  if (lammps == 1) {
+    int natoms = lammps_get_natoms(lmp);
+    type = (int *) malloc(natoms*sizeof(double));
+    int i;
+    for (i = 0; i < natoms; i++) type[i] = 1;
+
+    lammps_command(lmp,"delete_atoms group all");
+    lammps_create_atoms(lmp,natoms,NULL,type,x,v);
+    lammps_command(lmp,"run 10");
+  }
+
+  if (x) free(x);
+  if (v) free(v);
+  if (type) free(type);
+
+  // close down LAMMPS
+
+  lammps_close(lmp);

  /* close down MPI */

+  if (lammps == 1) MPI_Comm_free(&comm_lammps);
+  MPI_Barrier(MPI_COMM_WORLD);
  MPI_Finalize();
 }
--- a/examples/COUPLE/simple/simple.cpp
+++ b/examples/COUPLE/simple/simple.cpp
@ -77,7 +77,7 @@ int main(int narg, char **arg)
  // (could just send it to proc 0 of comm_lammps and let it Bcast)
  // all LAMMPS procs call input->one() on the line
  
-  LAMMPS *lmp;
+  LAMMPS *lmp = NULL;
  if (lammps == 1) lmp = new LAMMPS(0,NULL,comm_lammps);

  int n;
@ -91,7 +91,7 @@ int main(int narg, char **arg)
    MPI_Bcast(&n,1,MPI_INT,0,MPI_COMM_WORLD);
    if (n == 0) break;
    MPI_Bcast(line,n,MPI_CHAR,0,MPI_COMM_WORLD);
-    if (lammps == 1) lmp->input->one(line);
+    if (lammps == 1) lammps_command(lmp,line);
  }

  // run 10 more steps
@ -100,23 +100,74 @@ int main(int narg, char **arg)
  // put coords back into LAMMPS
  // run a single step with changed coords

+  double *x = NULL;
+  double *v = NULL;
+
  if (lammps == 1) {
    lmp->input->one("run 10");

    int natoms = static_cast<int> (lmp->atom->natoms);
-    double *x = new double[3*natoms];
+    x = new double[3*natoms];
+    v = new double[3*natoms];
    lammps_gather_atoms(lmp,"x",1,3,x);
+    lammps_gather_atoms(lmp,"v",1,3,v);
    double epsilon = 0.1;
    x[0] += epsilon;
    lammps_scatter_atoms(lmp,"x",1,3,x);
-    delete [] x;

+    // these 2 lines are the same
+
+    // lammps_command(lmp,"run 1");
    lmp->input->one("run 1");
  }

-  if (lammps == 1) delete lmp;
+  // extract force on single atom two different ways
+
+  if (lammps == 1) {
+    double **f = (double **) lammps_extract_atom(lmp,"f");
+    printf("Force on 1 atom via extract_atom: %g\n",f[0][0]);
+
+    double *fx = (double *) lammps_extract_variable(lmp,"fx","all");
+    printf("Force on 1 atom via extract_variable: %g\n",fx[0]);
+  }
+
+  // use commands_string() and commands_list() to invoke more commands
+
+  char *strtwo = "run 10\nrun 20";
+  if (lammps == 1) lammps_commands_string(lmp,strtwo);
+
+  char *cmds[2];
+  cmds[0] = "run 10";
+  cmds[1] = "run 20";
+  if (lammps == 1) lammps_commands_list(lmp,2,cmds);
+
+  // delete all atoms
+  // create_atoms() to create new ones with old coords, vels
+  // initial thermo should be same as step 20
+
+  int *type = NULL;
+
+  if (lammps == 1) {
+    int natoms = static_cast<int> (lmp->atom->natoms);
+    type = new int[natoms];
+    for (int i = 0; i < natoms; i++) type[i] = 1;
+
+    lmp->input->one("delete_atoms group all");
+    lammps_create_atoms(lmp,natoms,NULL,type,x,v);
+    lmp->input->one("run 10");
+  }
+
+  delete [] x;
+  delete [] v;
+  delete [] type;
+
+  // close down LAMMPS
+
+  delete lmp;

  // close down MPI

+  if (lammps == 1) MPI_Comm_free(&comm_lammps);
+  MPI_Barrier(MPI_COMM_WORLD);
  MPI_Finalize();
 }
--- a/lib/colvars/colvar.cpp
+++ b/lib/colvars/colvar.cpp
@ -1144,9 +1144,9 @@ cvm::real colvar::update_forces_energy()
    // For a periodic colvar, both walls may be applicable at the same time
    // in which case we pick the closer one
    if ( (!is_enabled(f_cv_upper_wall)) ||
-         (this->dist2(x, lower_wall) < this->dist2(x, upper_wall)) ) {
+         (this->dist2(x_reported, lower_wall) < this->dist2(x_reported, upper_wall)) ) {

-      cvm::real const grad = this->dist2_lgrad(x, lower_wall);
+      cvm::real const grad = this->dist2_lgrad(x_reported, lower_wall);
      if (grad < 0.0) {
        fw = -0.5 * lower_wall_k * grad;
        f += fw;
@ -1157,7 +1157,7 @@ cvm::real colvar::update_forces_energy()

    } else {

-      cvm::real const grad = this->dist2_lgrad(x, upper_wall);
+      cvm::real const grad = this->dist2_lgrad(x_reported, upper_wall);
      if (grad > 0.0) {
        fw = -0.5 * upper_wall_k * grad;
        f += fw;
@ -1177,15 +1177,21 @@ cvm::real colvar::update_forces_energy()
    // atoms only feel the harmonic force
    // fr: bias force on extended variable (without harmonic spring), for output in trajectory
    // f_ext: total force on extended variable (including harmonic spring)
-    // f: - initially, external biasing force
+    // f: - initially, external biasing force (including wall forces)
    //    - after this code block, colvar force to be applied to atomic coordinates, ie. spring force
    fr    = f;
    f_ext = f + (-0.5 * ext_force_k) * this->dist2_lgrad(xr, x);
    f     =     (-0.5 * ext_force_k) * this->dist2_rgrad(xr, x);

-    // The total force acting on the extended variable is f_ext
-    // This will be used in the next timestep
-    ft_reported = f_ext;
+    if (is_enabled(f_cv_subtract_applied_force)) {
+      // Report a "system" force without the biases on this colvar
+      // that is, just the spring force
+      ft_reported = (-0.5 * ext_force_k) * this->dist2_lgrad(xr, x);
+    } else {
+      // The total force acting on the extended variable is f_ext
+      // This will be used in the next timestep
+      ft_reported = f_ext;
+    }

    // leapfrog: starting from x_i, f_i, v_(i-1/2)
    vr  += (0.5 * dt) * f_ext / ext_mass;
--- a/lib/colvars/colvaratoms.cpp
+++ b/lib/colvars/colvaratoms.cpp
@ -367,6 +367,7 @@ int cvm::atom_group::parse(std::string const &conf)
        cvm::error("Error: atomsColValue, if provided, must be non-zero.\n", INPUT_ERROR);
      }

+      // NOTE: calls to add_atom() and/or add_atom_id() are in the proxy-implemented function
      cvm::load_atoms(atoms_file_name.c_str(), *this, atoms_col, atoms_col_value);
    }
  }
@ -403,11 +404,21 @@ int cvm::atom_group::parse(std::string const &conf)
    }
  }

+  // We need to know the fitting options to decide whether the group is scalable
+  parse_error |= parse_fitting_options(group_conf);
+
+  if (is_available(f_ag_scalable_com) && !b_rotate) {
+    enable(f_ag_scalable_com);
+    enable(f_ag_scalable);
+  }
+
  if (is_enabled(f_ag_scalable) && !b_dummy) {
+    cvm::log("Enabling scalable calculation for group \""+this->key+"\".\n");
    index = (cvm::proxy)->init_atom_group(atoms_ids);
  }

-  parse_error |= parse_fitting_options(group_conf);
+  bool b_print_atom_ids = false;
+  get_keyval(group_conf, "printAtomIDs", b_print_atom_ids, false, colvarparse::parse_silent);

  // TODO move this to colvarparse object
  check_keywords(group_conf, key.c_str());
@ -427,6 +438,10 @@ int cvm::atom_group::parse(std::string const &conf)
 	    cvm::to_str(total_mass)+", total charge = "+
            cvm::to_str(total_charge)+".\n");

+  if (b_print_atom_ids) {
+    cvm::log("Internal definition of the atom group:\n");
+  }
+
  cvm::decrease_depth();

  return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
@ -583,6 +598,21 @@ int cvm::atom_group::add_atom_name_residue_range(std::string const &psf_segid,
 }


+std::string const cvm::atom_group::print_atom_ids() const
+{
+  size_t line_count = 0;
+  std::ostringstream os("");
+  for (size_t i = 0; i < atoms_ids.size(); i++) {
+    os << " " << std::setw(9) << atoms_ids[i];
+    if (++line_count == 7) {
+      os << "\n";
+      line_count = 0;
+    }
+  }
+  return os.str();
+}
+
+
 int cvm::atom_group::parse_fitting_options(std::string const &group_conf)
 {
  bool b_defined_center = get_keyval(group_conf, "centerReference", b_center, false);
@ -1118,8 +1148,7 @@ void cvm::atom_group::apply_colvar_force(cvm::real const &force)
    log("Communicating a colvar force from atom group to the MD engine.\n");
  }

-  if (b_dummy)
-    return;
+  if (b_dummy) return;

  if (noforce) {
    cvm::error("Error: sending a force to a group that has "
@ -1161,17 +1190,21 @@ void cvm::atom_group::apply_colvar_force(cvm::real const &force)

 void cvm::atom_group::apply_force(cvm::rvector const &force)
 {
-  if (b_dummy)
-    return;
+  if (cvm::debug()) {
+    log("Communicating a colvar force from atom group to the MD engine.\n");
+  }
+
+  if (b_dummy) return;

  if (noforce) {
    cvm::error("Error: sending a force to a group that has "
-               "\"disableForces\" defined.\n");
+               "\"enableForces\" set to off.\n");
    return;
  }

  if (is_enabled(f_ag_scalable)) {
    (cvm::proxy)->apply_atom_group_force(index, force);
+    return;
  }

  if (b_rotate) {
--- a/lib/colvars/colvaratoms.h
+++ b/lib/colvars/colvaratoms.h
@ -253,6 +253,8 @@ public:
    return atoms.size();
  }

+  std::string const print_atom_ids() const;
+
  /// \brief If this option is on, this group merely acts as a wrapper
  /// for a fixed position; any calls to atoms within or to
  /// functions that return disaggregated data will fail
--- a/lib/colvars/colvarbias_abf.cpp
+++ b/lib/colvars/colvarbias_abf.cpp
@ -80,6 +80,7 @@ int colvarbias_abf::init(std::string const &conf)

  if (update_bias) {
  // Request calculation of total force (which also checks for availability)
+  // TODO - change this to a dependency - needs ABF-specific features
    if(enable(f_cvb_get_total_force)) return cvm::get_error();
  }

@ -133,6 +134,10 @@ int colvarbias_abf::init(std::string const &conf)

  // Data for eABF z-based estimator
  if (b_extended) {
+    // CZAR output files for stratified eABF
+    get_keyval(conf, "writeCZARwindowFile", b_czar_window_file, false,
+               colvarparse::parse_silent);
+
    z_bin.assign(colvars.size(), 0);
    z_samples   = new colvar_grid_count(colvars);
    z_samples->request_actual_value();
@ -241,7 +246,7 @@ int colvarbias_abf::update()

      for (size_t i = 0; i < colvars.size(); i++) {
        // get total forces (lagging by 1 timestep) from colvars
-        // and subtract previous ABF force
+        // and subtract previous ABF force if necessary
        update_system_force(i);
      }
      gradients->acc_force(force_bin, system_force);
@ -457,28 +462,30 @@ void colvarbias_abf::write_gradients_samples(const std::string &prefix, bool app

  if (z_gradients) {
    // Write eABF-related quantities
+
    std::string  z_samples_out_name = prefix + ".zcount";
-    std::string  z_gradients_out_name = prefix + ".zgrad";
-    std::string  czar_gradients_out_name = prefix + ".czar";
    cvm::ofstream z_samples_os;
-    cvm::ofstream z_gradients_os;
-    cvm::ofstream czar_gradients_os;

    if (!append) cvm::backup_file(z_samples_out_name.c_str());
    z_samples_os.open(z_samples_out_name.c_str(), mode);
    if (!z_samples_os.is_open()) {
-      cvm::error("Error opening ABF z sample file " + z_samples_out_name + " for writing");
+      cvm::error("Error opening eABF z-histogram file " + z_samples_out_name + " for writing");
    }
    z_samples->write_multicol(z_samples_os);
    z_samples_os.close();

-    if (!append) cvm::backup_file(z_gradients_out_name.c_str());
-    z_gradients_os.open(z_gradients_out_name.c_str(), mode);
-    if (!z_gradients_os.is_open()) {
-      cvm::error("Error opening ABF z gradient file " + z_gradients_out_name + " for writing");
+    if (b_czar_window_file) {
+      std::string  z_gradients_out_name = prefix + ".zgrad";
+      cvm::ofstream z_gradients_os;
+
+      if (!append) cvm::backup_file(z_gradients_out_name.c_str());
+      z_gradients_os.open(z_gradients_out_name.c_str(), mode);
+      if (!z_gradients_os.is_open()) {
+        cvm::error("Error opening eABF z-gradient file " + z_gradients_out_name + " for writing");
+      }
+      z_gradients->write_multicol(z_gradients_os);
+      z_gradients_os.close();
    }
-    z_gradients->write_multicol(z_gradients_os);
-    z_gradients_os.close();

    // Calculate CZAR estimator of gradients
    for (std::vector<int> ix = czar_gradients->new_index();
@ -490,6 +497,9 @@ void colvarbias_abf::write_gradients_samples(const std::string &prefix, bool app
      }
    }

+    std::string  czar_gradients_out_name = prefix + ".czar.grad";
+    cvm::ofstream czar_gradients_os;
+
    if (!append) cvm::backup_file(czar_gradients_out_name.c_str());
    czar_gradients_os.open(czar_gradients_out_name.c_str(), mode);
    if (!czar_gradients_os.is_open()) {
@ -499,17 +509,7 @@ void colvarbias_abf::write_gradients_samples(const std::string &prefix, bool app
    czar_gradients_os.close();

    if (colvars.size() == 1) {
-      std::string  z_pmf_out_name = prefix + ".zpmf";
-      if (!append) cvm::backup_file(z_pmf_out_name.c_str());
-      cvm::ofstream z_pmf_os;
-      // Do numerical integration and output a PMF
-      z_pmf_os.open(z_pmf_out_name.c_str(), mode);
-      if (!z_pmf_os.is_open())  cvm::error("Error opening z pmf file " + z_pmf_out_name + " for writing");
-      z_gradients->write_1D_integral(z_pmf_os);
-      z_pmf_os << std::endl;
-      z_pmf_os.close();
-
-      std::string  czar_pmf_out_name = prefix + ".czarpmf";
+      std::string  czar_pmf_out_name = prefix + ".czar.pmf";
      if (!append) cvm::backup_file(czar_pmf_out_name.c_str());
      cvm::ofstream czar_pmf_os;
      // Do numerical integration and output a PMF
@ -520,8 +520,6 @@ void colvarbias_abf::write_gradients_samples(const std::string &prefix, bool app
      czar_pmf_os.close();
    }
  }
-
-
  return;
 }

@ -559,7 +557,7 @@ void colvarbias_abf::read_gradients_samples()

    std::ifstream is;

-    cvm::log("Reading sample count from " + samples_in_name + " and gradients from " + gradients_in_name);
+    cvm::log("Reading sample count from " + samples_in_name + " and gradient from " + gradients_in_name);
    is.open(samples_in_name.c_str());
    if (!is.is_open()) cvm::error("Error opening ABF samples file " + samples_in_name + " for reading");
    samples->read_multicol(is, true);
@ -572,17 +570,18 @@ void colvarbias_abf::read_gradients_samples()
    is.close();

    if (z_gradients) {
-      cvm::log("Reading z sample count from " + z_samples_in_name + " and z gradients from " + z_gradients_in_name);
+      // Read eABF z-averaged data for CZAR
+      cvm::log("Reading z-histogram from " + z_samples_in_name + " and z-gradient from " + z_gradients_in_name);

      is.clear();
      is.open(z_samples_in_name.c_str());
-      if (!is.is_open())  cvm::error("Error opening ABF z samples file " + z_samples_in_name + " for reading");
+      if (!is.is_open())  cvm::error("Error opening eABF z-histogram file " + z_samples_in_name + " for reading");
      z_samples->read_multicol(is, true);
      is.close();
      is.clear();

      is.open(z_gradients_in_name.c_str());
-      if (!is.is_open())  cvm::error("Error opening ABF z gradient file " + z_gradients_in_name + " for reading");
+      if (!is.is_open())  cvm::error("Error opening eABF z-gradient file " + z_gradients_in_name + " for reading");
      z_gradients->read_multicol(is, true);
      is.close();
    }
--- a/lib/colvars/colvarbias_abf.h
+++ b/lib/colvars/colvarbias_abf.h
@ -40,6 +40,8 @@ private:
  int		output_freq;
  /// Write combined files with a history of all output data?
  bool      b_history_files;
+  /// Write CZAR output file for stratified eABF (.zgrad)
+  bool      b_czar_window_file;
  size_t    history_freq;

  /// Cap applied biasing force?
--- a/lib/colvars/colvarbias_meta.cpp
+++ b/lib/colvars/colvarbias_meta.cpp
@ -159,6 +159,23 @@ int colvarbias_meta::init(std::string const &conf)
    cvm::log("The bias temperature is "+cvm::to_str(bias_temperature)+".\n");
  }

+
+  // for ebmeta
+  target_dist = NULL;
+  get_keyval(conf, "ebMeta", ebmeta, false);
+  if(ebmeta){
+    target_dist = new colvar_grid_scalar();
+    target_dist->init_from_colvars(colvars);
+    get_keyval(conf, "targetdistfile", target_dist_file);
+    std::ifstream targetdiststream(target_dist_file.c_str());
+    target_dist->read_multicol(targetdiststream);
+    // normalize target distribution and multiply by effective volume = exp(differential entropy)
+    target_dist->multiply_constant(1.0/target_dist->integral());
+    cvm::real volume = std::exp(target_dist->entropy());
+    target_dist->multiply_constant(volume);
+    get_keyval(conf, "ebMetaEquilSteps", ebmeta_equil_steps, 0);
+  }
+
  if (cvm::debug())
    cvm::log("Done initializing the metadynamics bias \""+this->name+"\""+
             ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+".\n");
@ -186,6 +203,11 @@ colvarbias_meta::~colvarbias_meta()
  if (hills_traj_os.is_open())
    hills_traj_os.close();

+  if(target_dist) {
+    delete target_dist;
+    target_dist = NULL;
+  }
+
  if (cvm::n_meta_biases > 0)
    cvm::n_meta_biases -= 1;
 }
@ -221,9 +243,7 @@ colvarbias_meta::create_hill(colvarbias_meta::hill const &h)
  // output to trajectory (if specified)
  if (hills_traj_os.is_open()) {
    hills_traj_os << (hills.back()).output_traj();
-    if (cvm::debug()) {
-      hills_traj_os.flush();
-    }
+    hills_traj_os.flush();
  }

  has_data = true;
@ -258,8 +278,7 @@ colvarbias_meta::delete_hill(hill_iter &h)
    hills_traj_os << "# DELETED this hill: "
                  << (hills.back()).output_traj()
                  << "\n";
-    if (cvm::debug())
-      hills_traj_os.flush();
+    hills_traj_os.flush();
  }

  return hills.erase(h);
@ -381,38 +400,37 @@ int colvarbias_meta::update()
               ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+
               ": adding a new hill at step "+cvm::to_str(cvm::step_absolute())+".\n");

+    cvm::real hills_scale=1.0;
+
+    if (ebmeta) {
+       hills_scale *= 1.0/target_dist->value(target_dist->get_colvars_index());
+       if(cvm::step_absolute() <= ebmeta_equil_steps) {
+         cvm::real const hills_lambda=(cvm::real(ebmeta_equil_steps - cvm::step_absolute()))/(cvm::real(ebmeta_equil_steps));
+         hills_scale = hills_lambda + (1-hills_lambda)*hills_scale;
+       }
+    }
+
+    if (well_tempered) {
+      cvm::real hills_energy_sum_here = 0.0;
+      if (use_grids) {
+        std::vector<int> curr_bin = hills_energy->get_colvars_index();
+        hills_energy_sum_here = hills_energy->value(curr_bin);
+      } else {
+        calc_hills(new_hills_begin, hills.end(), hills_energy_sum_here);
+      }
+      hills_scale *= std::exp(-1.0*hills_energy_sum_here/(bias_temperature*cvm::boltzmann()));
+    }
+
    switch (comm) {

    case single_replica:
-      if (well_tempered) {
-        cvm::real hills_energy_sum_here = 0.0;
-        if (use_grids) {
-          std::vector<int> curr_bin = hills_energy->get_colvars_index();
-          hills_energy_sum_here = hills_energy->value(curr_bin);
-        } else {
-          calc_hills(new_hills_begin, hills.end(), hills_energy_sum_here);
-        }
-        cvm::real const exp_weight = std::exp(-1.0*hills_energy_sum_here/(bias_temperature*cvm::boltzmann()));
-        create_hill(hill((hill_weight*exp_weight), colvars, hill_width));
-      } else {
-        create_hill(hill(hill_weight, colvars, hill_width));
-      }
+
+      create_hill(hill(hill_weight*hills_scale, colvars, hill_width));
+
      break;

    case multiple_replicas:
-      if (well_tempered) {
-        cvm::real hills_energy_sum_here = 0.0;
-        if (use_grids) {
-          std::vector<int> curr_bin = hills_energy->get_colvars_index();
-          hills_energy_sum_here = hills_energy->value(curr_bin);
-        } else {
-          calc_hills(new_hills_begin, hills.end(), hills_energy_sum_here);
-        }
-        cvm::real const exp_weight = std::exp(-1.0*hills_energy_sum_here/(bias_temperature*cvm::boltzmann()));
-        create_hill(hill((hill_weight*exp_weight), colvars, hill_width, replica_id));
-      } else {
-        create_hill(hill(hill_weight, colvars, hill_width, replica_id));
-      }
+      create_hill(hill(hill_weight*hills_scale, colvars, hill_width, replica_id));
      if (replica_hills_os.is_open()) {
        replica_hills_os << hills.back();
      } else {
@ -1507,7 +1525,9 @@ int colvarbias_meta::setup_output()
                                       ("."+replica_id) :
                                       ("") )+
                                     ".hills.traj");
-    hills_traj_os.open(traj_file_name.c_str());
+    if (!hills_traj_os.is_open()) {
+      hills_traj_os.open(traj_file_name.c_str());
+    }
    if (!hills_traj_os.is_open())
      cvm::error("Error: in opening hills output file \"" +
                 traj_file_name+"\".\n", FILE_ERROR);
--- a/lib/colvars/colvarbias_meta.h
+++ b/lib/colvars/colvarbias_meta.h
@ -147,6 +147,14 @@ protected:
  /// \brief Biasing temperature in well-tempered metadynamics
  cvm::real  bias_temperature;

+  // EBmeta parameters
+  bool       ebmeta;
+  colvar_grid_scalar* target_dist;
+  std::string target_dist_file;
+  cvm::real target_dist_volume;
+  size_t ebmeta_equil_steps;
+
+
  /// \brief Try to read the restart information by allocating new
  /// grids before replacing the current ones (used e.g. in
  /// multiple_replicas)
--- a/lib/colvars/colvarcomp.cpp
+++ b/lib/colvars/colvarcomp.cpp
@ -84,15 +84,12 @@ cvm::atom_group *colvar::cvc::parse_group(std::string const &conf,
      if (is_available(f_cvc_scalable_com) && is_available(f_cvc_com_based)) {
        enable(f_cvc_scalable_com);
        enable(f_cvc_scalable);
-        group->enable(f_ag_scalable_com);
-        group->enable(f_ag_scalable);
+        // The CVC makes the feature available;
+        // the atom group will enable it unless it needs to compute a rotational fit
+        group->provide(f_ag_scalable_com);
      }

      // TODO check for other types of parallelism here
-
-      if (is_enabled(f_cvc_scalable)) {
-        cvm::log("Will enable scalable calculation for group \""+group->key+"\".\n");
-      }
    }

    if (group->parse(conf) == COLVARS_OK) {
--- a/lib/colvars/colvardeps.cpp
+++ b/lib/colvars/colvardeps.cpp
@ -449,8 +449,10 @@ void colvardeps::init_ag_requires() {

  // Features that are implemented (or not) by all atom groups
  feature_states[f_ag_active]->available = true;
-  feature_states[f_ag_scalable_com]->available = (cvm::proxy->scalable_group_coms() == COLVARS_OK);
-  feature_states[f_ag_scalable]->available = feature_states[f_ag_scalable_com]->available;
+  // f_ag_scalable_com is provided by the CVC iff it is COM-based
+  feature_states[f_ag_scalable_com]->available = false;
+  // TODO make f_ag_scalable depend on f_ag_scalable_com (or something else)
+  feature_states[f_ag_scalable]->available = true;
 }


--- a/lib/colvars/colvarmodule.cpp
+++ b/lib/colvars/colvarmodule.cpp
@ -156,6 +156,12 @@ int colvarmodule::parse_global_params(std::string const &conf)
    read_index_file(index_file_name.c_str());
  }

+  if (parse->get_keyval(conf, "smp", proxy->b_smp_active, proxy->b_smp_active)) {
+    if (proxy->b_smp_active == false) {
+      cvm::log("SMP parallelism has been disabled.\n");
+    }
+  }
+
  parse->get_keyval(conf, "analysis", b_analysis, b_analysis);

  parse->get_keyval(conf, "debugGradientsStepSize", debug_gradients_step_size,
@ -810,6 +816,7 @@ int colvarmodule::analyze()

 int colvarmodule::setup()
 {
+  if (this->size() == 0) return cvm::get_error();
  // loop over all components of all colvars to reset masses of all groups
  for (std::vector<colvar *>::iterator cvi = colvars.begin();
       cvi != colvars.end();  cvi++) {
@ -867,25 +874,35 @@ int colvarmodule::reset()

 int colvarmodule::setup_input()
 {
-  // name of input state file
-  restart_in_name = proxy->input_prefix().size() ?
-    std::string(proxy->input_prefix()+".colvars.state") :
-    std::string("") ;
+  if (this->size() == 0) return cvm::get_error();
+
+  std::string restart_in_name("");

  // read the restart configuration, if available
-  if (restart_in_name.size()) {
+  if (proxy->input_prefix().size()) {
    // read the restart file
+    restart_in_name = proxy->input_prefix();
    std::ifstream input_is(restart_in_name.c_str());
    if (!input_is.good()) {
-      cvm::error("Error: in opening restart file \""+
+      // try by adding the suffix
+      input_is.clear();
+      restart_in_name = restart_in_name+std::string(".colvars.state");
+      input_is.open(restart_in_name.c_str());
+    }
+
+    if (!input_is.good()) {
+      cvm::error("Error: in opening input file \""+
                 std::string(restart_in_name)+"\".\n",
                 FILE_ERROR);
      return COLVARS_ERROR;
    } else {
+      cvm::log(cvm::line_marker);
      cvm::log("Restarting from file \""+restart_in_name+"\".\n");
      read_restart(input_is);
      if (cvm::get_error() != COLVARS_OK) {
        return COLVARS_ERROR;
+      } else {
+        proxy->input_prefix().clear();
      }
      cvm::log(cvm::line_marker);
    }
@ -897,7 +914,9 @@ int colvarmodule::setup_input()

 int colvarmodule::setup_output()
 {
-  int error_code = 0;
+  if (this->size() == 0) return cvm::get_error();
+
+  int error_code = COLVARS_OK;

  // output state file (restart)
  restart_out_name = proxy->restart_output_prefix().size() ?
@ -1545,11 +1564,7 @@ std::list<std::string> colvarmodule::index_group_names;
 std::list<std::vector<int> > colvarmodule::index_groups;
 bool     colvarmodule::use_scripted_forces = false;
 bool     colvarmodule::scripting_after_biases = true;
-
-// file name prefixes
 std::string colvarmodule::output_prefix = "";
-std::string colvarmodule::restart_in_name = "";
-

 // i/o constants
 size_t const colvarmodule::it_width = 12;
--- a/lib/colvars/colvarmodule.h
+++ b/lib/colvars/colvarmodule.h
@ -4,7 +4,7 @@
 #define COLVARMODULE_H

 #ifndef COLVARS_VERSION
-#define COLVARS_VERSION "2016-10-05"
+#define COLVARS_VERSION "2016-10-21"
 #endif

 #ifndef COLVARS_DEBUG
@ -12,6 +12,9 @@
 #endif

 /*! \mainpage Main page
+This is the Developer's documentation for the Collective Variables Module.
+
+You can browse the class hierarchy or the list of source files.
 */

 /// \file colvarmodule.h
@ -154,9 +157,6 @@ public:
  /// Prefix for all output files for this run
  static std::string output_prefix;

-  /// input restart file name (determined from input_prefix)
-  static std::string restart_in_name;
-

  /// Array of collective variables
  static std::vector<colvar *>     colvars;
@ -197,6 +197,11 @@ public:
    return COLVARS_DEBUG;
  }

+  /// \brief How many objects are configured yet?
+  inline size_t const size() const
+  {
+    return colvars.size() + biases.size();
+  }

  /// \brief Constructor \param config_name Configuration file name
  /// \param restart_name (optional) Restart file name
--- a/lib/colvars/colvarparse.cpp
+++ b/lib/colvars/colvarparse.cpp
@ -644,9 +644,9 @@ bool colvarparse::key_lookup(std::string const &conf,

      // find the matching closing brace

-      if (cvm::debug()) {
-        cvm::log("Multi-line value, config is now \""+line+"\".\n");
-      }
+//       if (cvm::debug()) {
+//         cvm::log("Multi-line value, config is now \""+line+"\".\n");
+//       }

      int brace_count = 1;

@ -689,9 +689,9 @@ bool colvarparse::key_lookup(std::string const &conf,
            line_end = nl;
          line.append(conf, line_begin, (line_end-line_begin));

-          if (cvm::debug()) {
-            cvm::log("Added a new line, config is now \""+line+"\".\n");
-          }
+//           if (cvm::debug()) {
+//             cvm::log("Added a new line, config is now \""+line+"\".\n");
+//           }
        }

        if (brace_count < 0) {
--- a/lib/colvars/colvarproxy.h
+++ b/lib/colvars/colvarproxy.h
@ -25,7 +25,7 @@ public:
  colvarmodule *colvars;

  /// Default constructor
-  inline colvarproxy() : script(NULL) {}
+  inline colvarproxy() : script(NULL), b_smp_active(true) {}

  /// Default destructor
  virtual ~colvarproxy() {}
@ -80,7 +80,7 @@ public:
  /// configuration)
  std::string input_prefix_str, output_prefix_str, restart_output_prefix_str;

-  inline std::string input_prefix()
+  inline std::string & input_prefix()
  {
    return input_prefix_str;
  }
@ -116,12 +116,15 @@ public:

  // ***************** SHARED-MEMORY PARALLELIZATION *****************

-  /// Whether or not threaded parallelization is available
+  /// Whether threaded parallelization is available (TODO: make this a cvm::deps feature)
  virtual int smp_enabled()
  {
    return COLVARS_NOT_IMPLEMENTED;
  }

+  /// Whether threaded parallelization should be used (TODO: make this a cvm::deps feature)
+  bool b_smp_active;
+
  /// Distribute calculation of colvars (and their components) across threads
  virtual int smp_colvars_loop()
  {
--- a/lib/colvars/colvarscript.cpp
+++ b/lib/colvars/colvarscript.cpp
@ -14,6 +14,36 @@ colvarscript::colvarscript(colvarproxy *p)
 {
 }

+
+extern "C" {
+
+  // Generic hooks; NAMD and VMD have Tcl-specific versions in the respective proxies
+
+  int run_colvarscript_command(int argc, const char **argv)
+  {
+    colvarproxy *cvp = cvm::proxy;
+    if (!cvp) {
+      return -1;
+    }
+    if (!cvp->script) {
+      cvm::error("Called run_colvarscript_command without a script object initialized.\n");
+      return -1;
+    }
+    return cvp->script->run(argc, argv);
+  }
+
+  const char * get_colvarscript_result()
+  {
+    colvarproxy *cvp = cvm::proxy;
+    if (!cvp->script) {
+      cvm::error("Called run_colvarscript_command without a script object initialized.\n");
+      return "";
+    }
+    return cvp->script->result.c_str();
+  }
+}
+
+
 /// Run method based on given arguments
 int colvarscript::run(int argc, char const *argv[]) {

@ -125,7 +155,7 @@ int colvarscript::run(int argc, char const *argv[]) {
      result = "Missing arguments\n" + help_string();
      return COLVARSCRIPT_ERROR;
    }
-    proxy->input_prefix_str = argv[2];
+    proxy->input_prefix() = argv[2];
    if (colvars->setup_input() == COLVARS_OK) {
      return COLVARS_OK;
    } else {
--- a/python/examples/simple.py
+++ b/python/examples/simple.py
@ -13,6 +13,8 @@

 from __future__ import print_function
 import sys
+import numpy as np
+import ctypes

 # parse command line

@ -51,17 +53,39 @@ for line in lines: lmp.command(line)

 lmp.command("run 10")
 x = lmp.gather_atoms("x",1,3)
+v = lmp.gather_atoms("v",1,3)
 epsilon = 0.1
 x[0] += epsilon
 lmp.scatter_atoms("x",1,3,x)
 lmp.command("run 1");

+# extract force on single atom two different ways
+
 f = lmp.extract_atom("f",3)
 print("Force on 1 atom via extract_atom: ",f[0][0])

 fx = lmp.extract_variable("fx","all",1)
 print("Force on 1 atom via extract_variable:",fx[0])

+# use commands_string() and commands_list() to invoke more commands
+
+strtwo = "run 10\nrun 20"
+lmp.commands_string(strtwo)
+
+cmds = ["run 10","run 20"]
+lmp.commands_list(cmds)
+
+# delete all atoms
+# create_atoms() to create new ones with old coords, vels
+# initial thermo should be same as step 20
+
+natoms = lmp.get_natoms()
+type = natoms*[1]
+
+lmp.command("delete_atoms group all");
+lmp.create_atoms(natoms,None,type,x,v);
+lmp.command("run 10");
+
 # uncomment if running in parallel via Pypar
 #print("Proc %d out of %d procs has" % (me,nprocs), lmp)
 #pypar.finalize()
--- a/python/lammps.py
+++ b/python/lammps.py
@ -172,6 +172,8 @@ class lammps(object):
    if file: file = file.encode()
    self.lib.lammps_file(self.lmp,file)

+  # send a single command
+    
  def command(self,cmd):
    if cmd: cmd = cmd.encode()
    self.lib.lammps_command(self.lmp,cmd)
@ -185,6 +187,19 @@ class lammps(object):
        raise MPIAbortException(error_msg)
      raise Exception(error_msg)

+  # send a list of commands
+
+  def commands_list(self,cmdlist):
+    args = (c_char_p * len(cmdlist))(*cmdlist)
+    self.lib.lammps_commands_list(self.lmp,len(cmdlist),args)
+    
+  # send a string of commands
+
+  def commands_string(self,multicmd):
+    self.lib.lammps_commands_string(self.lmp,c_char_p(multicmd))
+    
+  # extract global info
+    
  def extract_global(self,name,type):
    if name: name = name.encode()
    if type == 0:
@ -195,6 +210,8 @@ class lammps(object):
    ptr = self.lib.lammps_extract_global(self.lmp,name)
    return ptr[0]

+  # extract per-atom info
+  
  def extract_atom(self,name,type):
    if name: name = name.encode()
    if type == 0:
@ -209,6 +226,8 @@ class lammps(object):
    ptr = self.lib.lammps_extract_atom(self.lmp,name)
    return ptr

+  # extract compute info
+  
  def extract_compute(self,id,style,type):
    if id: id = id.encode()
    if type == 0:
@ -226,6 +245,7 @@ class lammps(object):
      return ptr
    return None

+  # extract fix info
  # in case of global datum, free memory for 1 double via lammps_free()
  # double was allocated by library interface function

@ -249,6 +269,7 @@ class lammps(object):
    else:
      return None

+  # extract variable info
  # free memory for 1 double or 1 vector of doubles via lammps_free()
  # for vector, must copy nlocal returned values to local c_double vector
  # memory was allocated by library interface function
@ -296,7 +317,14 @@ class lammps(object):
    return self.lib.lammps_get_natoms(self.lmp)

  # return vector of atom properties gathered across procs, ordered by atom ID
-
+  # name = atom property recognized by LAMMPS in atom->extract()
+  # type = 0 for integer values, 1 for double values
+  # count = number of per-atom valus, 1 for type or charge, 3 for x or f
+  # returned data is a 1d vector - doc how it is ordered?
+  # NOTE: how could we insure are converting to correct Python type
+  #   e.g. for Python list or NumPy, etc
+  #   ditto for extact_atom() above
+  
  def gather_atoms(self,name,type,count):
    if name: name = name.encode()
    natoms = self.lib.lammps_get_natoms(self.lmp)
@ -310,12 +338,38 @@ class lammps(object):
    return data

  # scatter vector of atom properties across procs, ordered by atom ID
-  # assume vector is of correct type and length, as created by gather_atoms()
-
+  # name = atom property recognized by LAMMPS in atom->extract()
+  # type = 0 for integer values, 1 for double values
+  # count = number of per-atom valus, 1 for type or charge, 3 for x or f
+  # assume data is of correct type and length, as created by gather_atoms()
+  # NOTE: how could we insure are passing correct type to LAMMPS
+  # e.g. for Python list or NumPy, etc
+    
  def scatter_atoms(self,name,type,count,data):
    if name: name = name.encode()
    self.lib.lammps_scatter_atoms(self.lmp,name,type,count,data)

+  # create N atoms on all procs
+  # N = global number of atoms
+  # id = ID of each atom (optional, can be None)
+  # type = type of each atom (1 to Ntypes) (required)
+  # x = coords of each atom as (N,3) array (required)
+  # v = velocity of each atom as (N,3) array (optional, can be None)
+  # NOTE: how could we insure are passing correct type to LAMMPS
+  #   e.g. for Python list or NumPy, etc
+  #   ditto for gather_atoms() above
+
+  def create_atoms(self,n,id,type,x,v):
+    if id:
+      id_lmp = (c_int * n)()
+      id_lmp[:] = id
+    else: id_lmp = id
+    type_lmp = (c_int * n)()
+    type_lmp[:] = type
+    self.lib.lammps_create_atoms(self.lmp,n,id_lmp,type_lmp,x,v)
+
+  # document this?
+    
  @property
  def uses_exceptions(self):
    try:
--- a/src/KOKKOS/Install.sh
+++ b/src/KOKKOS/Install.sh
@ -83,6 +83,10 @@ action fix_nvt_kokkos.cpp
 action fix_nvt_kokkos.h
 action fix_qeq_reax_kokkos.cpp fix_qeq_reax.cpp
 action fix_qeq_reax_kokkos.h fix_qeq_reax.h
+action fix_reaxc_bonds_kokkos.cpp fix_reaxc_bonds.cpp
+action fix_reaxc_bonds_kokkos.h fix_reaxc_bonds.h
+action fix_reaxc_species_kokkos.cpp fix_reaxc_species.cpp
+action fix_reaxc_species_kokkos.h fix_reaxc_species.h
 action fix_setforce_kokkos.cpp
 action fix_setforce_kokkos.h
 action fix_wall_reflect_kokkos.cpp
--- a/src/KOKKOS/fix_qeq_reax_kokkos.cpp
+++ b/src/KOKKOS/fix_qeq_reax_kokkos.cpp
@ -431,8 +431,8 @@ void FixQEqReaxKokkos<DeviceType>::compute_h_item(int ii, int &m_fill, const boo
      const F_FLOAT rsq = delx*delx + dely*dely + delz*delz;
      if (rsq > cutsq) continue;

-      const F_FLOAT r = sqrt(rsq);
      if (final) {
+        const F_FLOAT r = sqrt(rsq);
        d_jlist(m_fill) = j;
        const F_FLOAT shldij = d_shield(itype,jtype);
        d_val(m_fill) = calculate_H_k(r,shldij);
--- a/src/KOKKOS/fix_reaxc_bonds_kokkos.cpp
+++ b/src/KOKKOS/fix_reaxc_bonds_kokkos.cpp
@ -0,0 +1,124 @@
+/* ----------------------------------------------------------------------
+   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: Stan Moore (Sandia)
+------------------------------------------------------------------------- */
+
+#include <stdlib.h>
+#include <string.h>
+#include "fix_ave_atom.h"
+#include "fix_reaxc_bonds_kokkos.h"
+#include "atom.h"
+#include "update.h"
+#include "pair_reax_c_kokkos.h"
+#include "modify.h"
+#include "neighbor.h"
+#include "neigh_list.h"
+#include "neigh_request.h"
+#include "comm.h"
+#include "force.h"
+#include "compute.h"
+#include "input.h"
+#include "variable.h"
+#include "memory.h"
+#include "error.h"
+#include "reaxc_list.h"
+#include "reaxc_types.h"
+#include "reaxc_defs.h"
+#include "atom_masks.h"
+
+using namespace LAMMPS_NS;
+using namespace FixConst;
+
+/* ---------------------------------------------------------------------- */
+
+FixReaxCBondsKokkos::FixReaxCBondsKokkos(LAMMPS *lmp, int narg, char **arg) :
+  FixReaxCBonds(lmp, narg, arg)
+{
+  kokkosable = 1;
+  atomKK = (AtomKokkos *) atom;
+
+  datamask_read = EMPTY_MASK;
+  datamask_modify = EMPTY_MASK;
+}
+
+/* ---------------------------------------------------------------------- */
+
+FixReaxCBondsKokkos::~FixReaxCBondsKokkos()
+{
+
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixReaxCBondsKokkos::init()
+{
+  Pair *pair_kk = force->pair_match("reax/c/kk",1);
+  if (pair_kk == NULL) error->all(FLERR,"Cannot use fix reax/c/bonds without "
+                  "pair_style reax/c/kk");
+
+  FixReaxCBonds::init();
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixReaxCBondsKokkos::Output_ReaxC_Bonds(bigint ntimestep, FILE *fp)
+
+{
+  int i, j;
+  int nbuf_local;
+  int nlocal_max, numbonds, numbonds_max;
+  double *buf;
+  DAT::tdual_ffloat_1d k_buf;
+
+  int nlocal = atom->nlocal;
+  int nlocal_tot = static_cast<int> (atom->natoms);
+
+  numbonds = 0;
+  if (reaxc->execution_space == Device)
+    ((PairReaxCKokkos<LMPDeviceType>*) reaxc)->FindBond(numbonds);
+  else
+    ((PairReaxCKokkos<LMPHostType>*) reaxc)->FindBond(numbonds);
+
+  // allocate a temporary buffer for the snapshot info
+  MPI_Allreduce(&numbonds,&numbonds_max,1,MPI_INT,MPI_MAX,world);
+  MPI_Allreduce(&nlocal,&nlocal_max,1,MPI_INT,MPI_MAX,world);
+
+  nbuf = 1+(numbonds_max*2+10)*nlocal_max;
+  memory->create_kokkos(k_buf,buf,nbuf,"reax/c/bonds:buf");
+
+  // Pass information to buffer
+  if (reaxc->execution_space == Device)
+    ((PairReaxCKokkos<LMPDeviceType>*) reaxc)->PackBondBuffer(k_buf,nbuf_local);
+  else
+    ((PairReaxCKokkos<LMPHostType>*) reaxc)->PackBondBuffer(k_buf,nbuf_local);
+
+  // Receive information from buffer for output
+  RecvBuffer(buf, nbuf, nbuf_local, nlocal_tot, numbonds_max);
+
+  memory->destroy_kokkos(k_buf,buf);
+}
+
+double FixReaxCBondsKokkos::memory_usage()
+{
+  double bytes;
+
+  bytes = nbuf*sizeof(double);
+  // These are accounted for in PairReaxCKokkos:
+  //bytes += nmax*sizeof(int);
+  //bytes += 1.0*nmax*MAXREAXBOND*sizeof(double);
+  //bytes += 1.0*nmax*MAXREAXBOND*sizeof(int);
+
+  return bytes;
+}
--- a/src/KOKKOS/fix_reaxc_bonds_kokkos.h
+++ b/src/KOKKOS/fix_reaxc_bonds_kokkos.h
@ -0,0 +1,42 @@
+/* -*- c++ -*- ----------------------------------------------------------
+   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.
+------------------------------------------------------------------------- */
+
+#ifdef FIX_CLASS
+
+FixStyle(reax/c/bonds/kk,FixReaxCBondsKokkos)
+
+#else
+
+#ifndef LMP_FIX_REAXC_BONDS_KOKKOS_H
+#define LMP_FIX_REAXC_BONDS_KOKKOS_H
+
+#include "fix_reaxc_bonds.h"
+#include "kokkos_type.h"
+
+namespace LAMMPS_NS {
+
+class FixReaxCBondsKokkos : public FixReaxCBonds {
+ public:
+  FixReaxCBondsKokkos(class LAMMPS *, int, char **);
+  virtual ~FixReaxCBondsKokkos();
+  void init();
+
+ private:
+  double nbuf;
+  void Output_ReaxC_Bonds(bigint, FILE *);
+  double memory_usage();
+};
+}
+
+#endif
+#endif
--- a/src/KOKKOS/fix_reaxc_species_kokkos.cpp
+++ b/src/KOKKOS/fix_reaxc_species_kokkos.cpp
@ -0,0 +1,157 @@
+/* ----------------------------------------------------------------------
+   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 authors: Stan Moore (Sandia)
+------------------------------------------------------------------------- */
+
+#include <stdlib.h>
+#include <math.h>
+#include "atom.h"
+#include <string.h>
+#include "fix_ave_atom.h"
+#include "fix_reaxc_species_kokkos.h"
+#include "domain.h"
+#include "update.h"
+#include "pair_reax_c_kokkos.h"
+#include "modify.h"
+#include "neighbor.h"
+#include "neigh_list.h"
+#include "neigh_request.h"
+#include "comm.h"
+#include "force.h"
+#include "compute.h"
+#include "input.h"
+#include "variable.h"
+#include "memory.h"
+#include "error.h"
+#include "reaxc_list.h"
+#include "atom_masks.h"
+
+using namespace LAMMPS_NS;
+using namespace FixConst;
+
+/* ---------------------------------------------------------------------- */
+
+FixReaxCSpeciesKokkos::FixReaxCSpeciesKokkos(LAMMPS *lmp, int narg, char **arg) :
+  FixReaxCSpecies(lmp, narg, arg)
+{
+  kokkosable = 1;
+  atomKK = (AtomKokkos *) atom;
+  
+  // NOTE: Could improve performance if a Kokkos version of ComputeSpecAtom is added
+
+  datamask_read = X_MASK | V_MASK | Q_MASK | MASK_MASK;
+  datamask_modify = EMPTY_MASK;
+}
+
+/* ---------------------------------------------------------------------- */
+
+FixReaxCSpeciesKokkos::~FixReaxCSpeciesKokkos()
+{
+
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixReaxCSpeciesKokkos::init()
+{
+  Pair* pair_kk = force->pair_match("reax/c/kk",1);
+  if (pair_kk == NULL) error->all(FLERR,"Cannot use fix reax/c/species/kk without "
+                  "pair_style reax/c/kk");
+
+  FixReaxCSpecies::init();
+
+  int irequest = neighbor->request(this,instance_me);
+  neighbor->requests[irequest]->pair = 0;
+  neighbor->requests[irequest]->fix = 1;
+  neighbor->requests[irequest]->newton = 2;
+  neighbor->requests[irequest]->ghost = 1;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixReaxCSpeciesKokkos::FindMolecule()
+{
+  int i,j,ii,jj,inum,itype,jtype,loop,looptot;
+  int change,done,anychange;
+  int *mask = atom->mask;
+  int *ilist;
+  double bo_tmp,bo_cut;
+  double **spec_atom = f_SPECBOND->array_atom;
+
+  inum = list->inum;
+  ilist = list->ilist;
+
+  for (ii = 0; ii < inum; ii++) {
+    i = ilist[ii];
+    if (mask[i] & groupbit) {
+      clusterID[i] = atom->tag[i];
+      x0[i].x = spec_atom[i][1];
+      x0[i].y = spec_atom[i][2];
+      x0[i].z = spec_atom[i][3];
+    }
+    else clusterID[i] = 0.0;
+  }
+
+  loop = 0;
+  while (1) {
+    comm->forward_comm_fix(this);
+    loop ++;
+
+    change = 0;
+    while (1) {
+      done = 1;
+
+      for (ii = 0; ii < inum; ii++) {
+      	i = ilist[ii];
+      	if (!(mask[i] & groupbit)) continue;
+
+      	itype = atom->type[i];
+
+        for (jj = 0; jj < MAXSPECBOND; jj++) {
+      	  j = reaxc->tmpid[i][jj];
+
+      	  if (j < i) continue;
+      	  if (!(mask[j] & groupbit)) continue;
+
+      	  if (clusterID[i] == clusterID[j] && PBCconnected[i] == PBCconnected[j]
+	    && x0[i].x == x0[j].x && x0[i].y == x0[j].y && x0[i].z == x0[j].z) continue;
+
+          jtype = atom->type[j];
+      	  bo_cut = BOCut[itype][jtype];
+      	  bo_tmp = spec_atom[i][jj+7];
+
+      	  if (bo_tmp > bo_cut) {
+            clusterID[i] = clusterID[j] = MIN(clusterID[i], clusterID[j]);
+            PBCconnected[i] = PBCconnected[j] = MAX(PBCconnected[i], PBCconnected[j]);
+            x0[i] = x0[j] = chAnchor(x0[i], x0[j]);
+            if ((fabs(spec_atom[i][1] - spec_atom[j][1]) > reaxc->control->bond_cut)
+             || (fabs(spec_atom[i][2] - spec_atom[j][2]) > reaxc->control->bond_cut)
+             || (fabs(spec_atom[i][3] - spec_atom[j][3]) > reaxc->control->bond_cut))
+              PBCconnected[i] = PBCconnected[j] = 1;
+      	    done = 0;
+      	  }
+      	}
+      }
+      if (!done) change = 1;
+      if (done) break;
+    }
+    MPI_Allreduce(&change,&anychange,1,MPI_INT,MPI_MAX,world);
+    if (!anychange) break;
+
+    MPI_Allreduce(&loop,&looptot,1,MPI_INT,MPI_SUM,world);
+    if (looptot >= 400*nprocs) break;
+
+  }
+}
--- a/src/KOKKOS/fix_reaxc_species_kokkos.h
+++ b/src/KOKKOS/fix_reaxc_species_kokkos.h
@ -0,0 +1,41 @@
+/* -*- c++ -*- ----------------------------------------------------------
+   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.
+------------------------------------------------------------------------- */
+
+#ifdef FIX_CLASS
+
+FixStyle(reax/c/species/kk,FixReaxCSpeciesKokkos)
+
+#else
+
+#ifndef LMP_FIX_REAXC_SPECIES_KOKKOS_H
+#define LMP_FIX_REAXC_SPECIES_KOKKOS_H
+
+#include "fix_reaxc_species.h"
+
+#define BUFLEN 1000
+
+namespace LAMMPS_NS {
+
+class FixReaxCSpeciesKokkos : public FixReaxCSpecies {
+ public:
+  FixReaxCSpeciesKokkos(class LAMMPS *, int, char **);
+  virtual ~FixReaxCSpeciesKokkos();
+  void init();
+
+ private:
+  void FindMolecule();
+};
+}
+
+#endif
+#endif
--- a/src/KOKKOS/neighbor_kokkos.cpp
+++ b/src/KOKKOS/neighbor_kokkos.cpp
@ -181,6 +181,9 @@ int NeighborKokkos::init_lists_kokkos()

 void NeighborKokkos::init_list_flags1_kokkos(int i)
 {
+  if (style != BIN)
+    error->all(FLERR,"KOKKOS package only supports 'bin' neighbor lists");
+
  if (lists_host[i]) {
    lists_host[i]->buildflag = 1;
    if (pair_build_host[i] == NULL) lists_host[i]->buildflag = 0;
--- a/src/KOKKOS/neighbor_kokkos.h
+++ b/src/KOKKOS/neighbor_kokkos.h
@ -434,6 +434,10 @@ class NeighborKokkos : public Neighbor {

 /* ERROR/WARNING messages:

+E: KOKKOS package only supports 'bin' neighbor lists
+
+Self-explanatory.
+
 E: Too many local+ghost atoms for neighbor list

 The number of nlocal + nghost atoms on a processor
--- a/src/KOKKOS/pair_kokkos.h
+++ b/src/KOKKOS/pair_kokkos.h
@ -611,7 +611,7 @@ struct PairComputeFunctor<PairStyle,N2,STACKPARAMS,Specialisation>  {
 // pair_compute_neighlist and pair_compute_fullcluster will match - either the dummy version
 // or the real one further below.
 template<class PairStyle, unsigned NEIGHFLAG, class Specialisation>
-EV_FLOAT pair_compute_neighlist (PairStyle* fpair, typename Kokkos::Impl::enable_if<!(NEIGHFLAG&PairStyle::EnabledNeighFlags), NeighListKokkos<typename PairStyle::device_type>*>::type list) {
+EV_FLOAT pair_compute_neighlist (PairStyle* fpair, typename Kokkos::Impl::enable_if<!((NEIGHFLAG&PairStyle::EnabledNeighFlags) != 0), NeighListKokkos<typename PairStyle::device_type>*>::type list) {
  EV_FLOAT ev;
  (void) fpair;
  (void) list;
@ -620,7 +620,7 @@ EV_FLOAT pair_compute_neighlist (PairStyle* fpair, typename Kokkos::Impl::enable
 }

 template<class PairStyle, class Specialisation>
-EV_FLOAT pair_compute_fullcluster (PairStyle* fpair, typename Kokkos::Impl::enable_if<!(FULLCLUSTER&PairStyle::EnabledNeighFlags), NeighListKokkos<typename PairStyle::device_type>*>::type list) {
+EV_FLOAT pair_compute_fullcluster (PairStyle* fpair, typename Kokkos::Impl::enable_if<!((FULLCLUSTER&PairStyle::EnabledNeighFlags) != 0), NeighListKokkos<typename PairStyle::device_type>*>::type list) {
  EV_FLOAT ev;
  (void) fpair;
  (void) list;
@ -630,7 +630,7 @@ EV_FLOAT pair_compute_fullcluster (PairStyle* fpair, typename Kokkos::Impl::enab

 // Submit ParallelFor for NEIGHFLAG=HALF,HALFTHREAD,FULL,N2
 template<class PairStyle, unsigned NEIGHFLAG, class Specialisation>
-EV_FLOAT pair_compute_neighlist (PairStyle* fpair, typename Kokkos::Impl::enable_if<NEIGHFLAG&PairStyle::EnabledNeighFlags, NeighListKokkos<typename PairStyle::device_type>*>::type list) {
+EV_FLOAT pair_compute_neighlist (PairStyle* fpair, typename Kokkos::Impl::enable_if<(NEIGHFLAG&PairStyle::EnabledNeighFlags) != 0, NeighListKokkos<typename PairStyle::device_type>*>::type list) {
  EV_FLOAT ev;
  if(fpair->atom->ntypes > MAX_TYPES_STACKPARAMS) {
    PairComputeFunctor<PairStyle,NEIGHFLAG,false,Specialisation > ff(fpair,list);
@ -646,7 +646,7 @@ EV_FLOAT pair_compute_neighlist (PairStyle* fpair, typename Kokkos::Impl::enable

 // Submit ParallelFor for NEIGHFLAG=FULLCLUSTER
 template<class PairStyle, class Specialisation>
-EV_FLOAT pair_compute_fullcluster (PairStyle* fpair, typename Kokkos::Impl::enable_if<FULLCLUSTER&PairStyle::EnabledNeighFlags, NeighListKokkos<typename PairStyle::device_type>*>::type list) {
+EV_FLOAT pair_compute_fullcluster (PairStyle* fpair, typename Kokkos::Impl::enable_if<(FULLCLUSTER&PairStyle::EnabledNeighFlags) != 0, NeighListKokkos<typename PairStyle::device_type>*>::type list) {
  EV_FLOAT ev;
  if(fpair->atom->ntypes > MAX_TYPES_STACKPARAMS) {
    typedef PairComputeFunctor<PairStyle,FULLCLUSTER,false,Specialisation >
--- a/src/KOKKOS/pair_reax_c_kokkos.cpp
+++ b/src/KOKKOS/pair_reax_c_kokkos.cpp
@ -44,7 +44,7 @@

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

-namespace LAMMPS_NS{
+namespace LAMMPS_NS {
 using namespace MathConst;
 using namespace MathSpecial;

@ -69,6 +69,9 @@ PairReaxCKokkos<DeviceType>::PairReaxCKokkos(LAMMPS *lmp) : PairReaxC(lmp)
  nmax = 0;
  maxbo = 1;
  maxhb = 1;
+
+  k_error_flag = DAT::tdual_int_scalar("pair:error_flag");
+  k_nbuf_local = DAT::tdual_int_scalar("pair:nbuf_local");
 }

 /* ---------------------------------------------------------------------- */
@ -76,10 +79,15 @@ PairReaxCKokkos<DeviceType>::PairReaxCKokkos(LAMMPS *lmp) : PairReaxC(lmp)
 template<class DeviceType>
 PairReaxCKokkos<DeviceType>::~PairReaxCKokkos()
 {
-  if (!copymode) {
-    memory->destroy_kokkos(k_eatom,eatom);
-    memory->destroy_kokkos(k_vatom,vatom);
-  }
+  if (copymode) return;
+
+  memory->destroy_kokkos(k_eatom,eatom);
+  memory->destroy_kokkos(k_vatom,vatom);
+
+  memory->destroy_kokkos(k_tmpid,tmpid);
+  tmpid = NULL;
+  memory->destroy_kokkos(k_tmpbo,tmpbo);
+  tmpbo = NULL;
 }

 /* ---------------------------------------------------------------------- */
@ -306,6 +314,11 @@ void PairReaxCKokkos<DeviceType>::setup()
  bo_cut = 0.01 * gp[29];
  thb_cut = control->thb_cut;
  thb_cutsq = 0.000010; //thb_cut*thb_cut;
+
+  if (atom->nmax > nmax) {
+    nmax = atom->nmax;
+    allocate_array();
+  }
 }

 /* ---------------------------------------------------------------------- */
@ -980,6 +993,9 @@ void PairReaxCKokkos<DeviceType>::compute(int eflag_in, int vflag_in)
    k_vatom.template sync<LMPHostType>();
  }

+  if (fixspecies_flag)
+    FindBondSpecies();
+
  copymode = 0;
 }

@ -1346,6 +1362,19 @@ void PairReaxCKokkos<DeviceType>::allocate_array()
  d_Delta_lp_temp = typename AT::t_ffloat_1d("reax/c/kk:Delta_lp_temp",nmax);
  d_CdDelta = typename AT::t_ffloat_1d("reax/c/kk:CdDelta",nmax);
  d_sum_ovun = typename AT::t_ffloat_2d_dl("reax/c/kk:sum_ovun",nmax,3);
+
+  // FixReaxCSpecies
+  if (fixspecies_flag) {
+    memory->destroy_kokkos(k_tmpid,tmpid);
+    memory->destroy_kokkos(k_tmpbo,tmpbo);
+    memory->create_kokkos(k_tmpid,tmpid,nmax,MAXSPECBOND,"pair:tmpid");
+    memory->create_kokkos(k_tmpbo,tmpbo,nmax,MAXSPECBOND,"pair:tmpbo");
+  }
+
+  // FixReaxCBonds
+  d_abo = typename AT::t_ffloat_2d("reax/c/kk:abo",nmax,maxbo);
+  d_neighid = typename AT::t_tagint_2d("reax/c/kk:neighid",nmax,maxbo);
+  d_numneigh_bonds = typename AT::t_int_1d("reax/c/kk:numneigh_bonds",nmax);
 }

 /* ---------------------------------------------------------------------- */
@ -3954,11 +3983,214 @@ double PairReaxCKokkos<DeviceType>::memory_usage()
  bytes += nmax*17*sizeof(F_FLOAT);
  bytes += maxbo*nmax*34*sizeof(F_FLOAT);

+  // FixReaxCSpecies
+  if (fixspecies_flag) {
+    bytes += MAXSPECBOND*nmax*sizeof(tagint);
+    bytes += MAXSPECBOND*nmax*sizeof(F_FLOAT);
+  }
+
+  // FixReaxCBonds
+  bytes += maxbo*nmax*sizeof(tagint);
+  bytes += maxbo*nmax*sizeof(F_FLOAT);
+  bytes += nmax*sizeof(int);
+
  return bytes;
 }

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

+template<class DeviceType>
+void PairReaxCKokkos<DeviceType>::FindBond(int &numbonds)
+{
+  copymode = 1;
+  Kokkos::parallel_for(Kokkos::RangePolicy<DeviceType, PairReaxFindBondZero>(0,nmax),*this);
+  DeviceType::fence();
+
+  bo_cut_bond = control->bg_cut;
+
+  atomKK->sync(execution_space,TAG_MASK);
+  tag = atomKK->k_tag.view<DeviceType>();
+
+  const int inum = list->inum;
+  NeighListKokkos<DeviceType>* k_list = static_cast<NeighListKokkos<DeviceType>*>(list);
+  d_ilist = k_list->d_ilist;
+  k_list->clean_copy();
+
+  numbonds = 0;
+  PairReaxCKokkosFindBondFunctor<DeviceType> find_bond_functor(this);
+  Kokkos::parallel_reduce(inum,find_bond_functor,numbonds);
+  DeviceType::fence();
+  copymode = 0;
+}
+
+template<class DeviceType>
+KOKKOS_INLINE_FUNCTION
+void PairReaxCKokkos<DeviceType>::operator()(PairReaxFindBondZero, const int &i) const {
+  d_numneigh_bonds[i] = 0;
+  for (int j = 0; j < maxbo; j++) {
+    d_neighid(i,j) = 0;
+    d_abo(i,j) = 0.0;
+  }
+}
+
+template<class DeviceType>
+KOKKOS_INLINE_FUNCTION
+void PairReaxCKokkos<DeviceType>::calculate_find_bond_item(int ii, int &numbonds) const
+{
+  const int i = d_ilist[ii];
+  int nj = 0;
+
+  const int j_start = d_bo_first[i];
+  const int j_end = j_start + d_bo_num[i];
+  for (int jj = j_start; jj < j_end; jj++) {
+    int j = d_bo_list[jj];
+    j &= NEIGHMASK;
+    const tagint jtag = tag[j];
+    const int j_index = jj - j_start;
+    double bo_tmp = d_BO(i,j_index);
+
+    if (bo_tmp > bo_cut_bond) {
+      d_neighid(i,nj) = jtag;
+      d_abo(i,nj) = bo_tmp;
+      nj++;
+    }
+  }
+  d_numneigh_bonds[i] = nj;
+  if (nj > numbonds) numbonds = nj;
+}
+
+/* ---------------------------------------------------------------------- */
+
+template<class DeviceType>
+void PairReaxCKokkos<DeviceType>::PackBondBuffer(DAT::tdual_ffloat_1d k_buf, int &nbuf_local)
+{
+  d_buf = k_buf.view<DeviceType>();
+  k_params_sing.template sync<DeviceType>();
+  atomKK->sync(execution_space,TAG_MASK|TYPE_MASK|Q_MASK|MOLECULE_MASK);
+
+  tag = atomKK->k_tag.view<DeviceType>();
+  type = atomKK->k_type.view<DeviceType>();
+  q = atomKK->k_q.view<DeviceType>();
+  if (atom->molecule)
+    molecule = atomKK->k_molecule.view<DeviceType>();
+
+  copymode = 1;
+  nlocal = atomKK->nlocal;
+  PairReaxCKokkosPackBondBufferFunctor<DeviceType> pack_bond_buffer_functor(this);
+  Kokkos::parallel_scan(nlocal,pack_bond_buffer_functor);
+  DeviceType::fence();
+  copymode = 0;
+
+  k_buf.modify<DeviceType>();
+  k_nbuf_local.modify<DeviceType>();
+
+  k_buf.sync<LMPHostType>();
+  k_nbuf_local.sync<LMPHostType>();
+  nbuf_local = k_nbuf_local.h_view();
+}
+
+template<class DeviceType>
+KOKKOS_INLINE_FUNCTION
+void PairReaxCKokkos<DeviceType>::pack_bond_buffer_item(int i, int &j, const bool &final) const
+{
+  if (i == 0)
+    j += 2;
+
+  if (final) {
+    d_buf[j-1] = tag[i];
+    d_buf[j+0] = type[i];
+    d_buf[j+1] = d_total_bo[i];
+    d_buf[j+2] = paramssing(type[i]).nlp_opt - d_Delta_lp[i];
+    d_buf[j+3] = q[i];
+    d_buf[j+4] = d_numneigh_bonds[i];
+  }
+  const int numbonds = d_numneigh_bonds[i];
+
+  if (final) {
+    for (int k = 5; k < 5+numbonds; k++) {
+      d_buf[j+k] = d_neighid(i,k-5);
+    }
+  }
+  j += (5+numbonds);
+
+  if (final) {
+    if (!molecule.data()) d_buf[j] = 0.0;
+    else d_buf[j] = molecule[i];
+  }
+  j++;
+
+  if (final) {
+    for (int k = 0; k < numbonds; k++) {
+      d_buf[j+k] = d_abo(i,k);
+    }
+  }
+  j += (1+numbonds);
+
+  if (final && i == nlocal-1)
+    k_nbuf_local.view<DeviceType>()() = j - 1;
+}
+
+/* ---------------------------------------------------------------------- */
+
+template<class DeviceType>
+void PairReaxCKokkos<DeviceType>::FindBondSpecies()
+{
+  copymode = 1;
+  Kokkos::parallel_for(Kokkos::RangePolicy<DeviceType, PairReaxFindBondSpeciesZero>(0,nmax),*this);
+  DeviceType::fence();
+
+  nlocal = atomKK->nlocal;
+  Kokkos::parallel_for(Kokkos::RangePolicy<DeviceType, PairReaxFindBondSpecies>(0,nlocal),*this);
+  DeviceType::fence();
+  copymode = 0;
+
+  // NOTE: Could improve performance if a Kokkos version of ComputeSpecAtom is added
+
+  k_tmpbo.modify<DeviceType>();
+  k_tmpid.modify<DeviceType>();
+  k_error_flag.modify<DeviceType>();
+
+  k_tmpbo.sync<LMPHostType>();
+  k_tmpid.sync<LMPHostType>();
+  k_error_flag.sync<LMPHostType>();
+
+  if (k_error_flag.h_view())
+    error->all(FLERR,"Increase MAXSPECBOND in reaxc_defs.h");
+}
+
+template<class DeviceType>
+KOKKOS_INLINE_FUNCTION
+void PairReaxCKokkos<DeviceType>::operator()(PairReaxFindBondSpeciesZero, const int &i) const {
+  for (int j = 0; j < MAXSPECBOND; j++) {
+    k_tmpbo.view<DeviceType>()(i,j) = 0.0;
+    k_tmpid.view<DeviceType>()(i,j) = 0;
+  }
+}
+
+template<class DeviceType>
+KOKKOS_INLINE_FUNCTION
+void PairReaxCKokkos<DeviceType>::operator()(PairReaxFindBondSpecies, const int &i) const {
+  int nj = 0;
+
+  const int j_start = d_bo_first[i];
+  const int j_end = j_start + d_bo_num[i];
+  for (int jj = j_start; jj < j_end; jj++) {
+    int j = d_bo_list[jj];
+    j &= NEIGHMASK;
+    if (j < i) continue;
+    const int j_index = jj - j_start;
+  
+    double bo_tmp = d_BO(i,j_index);
+  
+    if (bo_tmp >= 0.10 ) { // Why is this a hardcoded value?
+      k_tmpid.view<DeviceType>()(i,nj) = j;
+      k_tmpbo.view<DeviceType>()(i,nj) = bo_tmp;
+      nj++;
+      if (nj > MAXSPECBOND) k_error_flag.view<DeviceType>()() = 1;
+    }
+  }
+}
+
 template class PairReaxCKokkos<LMPDeviceType>;
 #ifdef KOKKOS_HAVE_CUDA
 template class PairReaxCKokkos<LMPHostType>;
--- a/src/KOKKOS/pair_reax_c_kokkos.h
+++ b/src/KOKKOS/pair_reax_c_kokkos.h
@ -115,6 +115,13 @@ struct PairReaxComputeTorsion{};
 template<int NEIGHFLAG, int EVFLAG>
 struct PairReaxComputeHydrogen{};

+struct PairReaxFindBondZero{};
+
+struct PairReaxFindBondSpeciesZero{};
+
+struct PairReaxFindBondSpecies{};
+
+
 template<class DeviceType>
 class PairReaxCKokkos : public PairReaxC {
 public:
@ -132,6 +139,9 @@ class PairReaxCKokkos : public PairReaxC {
  void *extract(const char *, int &);
  void init_style();
  double memory_usage();
+  void FindBond(int &);
+  void PackBondBuffer(DAT::tdual_ffloat_1d, int &);
+  void FindBondSpecies();

  template<int NEIGHFLAG, int EVFLAG>
  KOKKOS_INLINE_FUNCTION
@ -245,6 +255,21 @@ class PairReaxCKokkos : public PairReaxC {
  KOKKOS_INLINE_FUNCTION
  void operator()(PairReaxComputeHydrogen<NEIGHFLAG,EVFLAG>, const int&) const;

+  KOKKOS_INLINE_FUNCTION
+  void operator()(PairReaxFindBondZero, const int&) const;
+
+  KOKKOS_INLINE_FUNCTION
+  void calculate_find_bond_item(int, int&) const;
+
+  KOKKOS_INLINE_FUNCTION
+  void pack_bond_buffer_item(int, int&, const bool&) const;
+
+  KOKKOS_INLINE_FUNCTION
+  void operator()(PairReaxFindBondSpeciesZero, const int&) const;
+
+  KOKKOS_INLINE_FUNCTION
+  void operator()(PairReaxFindBondSpecies, const int&) const;
+
  struct params_sing{
    KOKKOS_INLINE_FUNCTION
    params_sing(){mass=0;chi=0;eta=0;r_s=0;r_pi=0;r_pi2=0;valency=0;valency_val=0;valency_e=0;valency_boc=0;nlp_opt=0;
@ -359,8 +384,9 @@ class PairReaxCKokkos : public PairReaxC {
  typename AT::t_x_array_randomread x;
  typename AT::t_f_array f;
  typename AT::t_int_1d_randomread type;
-  typename AT::t_tagint_1d tag;
+  typename AT::t_tagint_1d_randomread tag;
  typename AT::t_float_1d_randomread q;
+  typename AT::t_tagint_1d_randomread molecule;

  DAT::tdual_efloat_1d k_eatom;
  typename AT::t_efloat_1d v_eatom;
@ -405,6 +431,7 @@ class PairReaxCKokkos : public PairReaxC {
  int neighflag,newton_pair, maxnumneigh, maxhb, maxbo;
  int nlocal,nall,eflag,vflag;
  F_FLOAT cut_nbsq, cut_hbsq, cut_bosq, bo_cut, thb_cut, thb_cutsq;
+  F_FLOAT bo_cut_bond;

  int vdwflag, lgflag;
  F_FLOAT gp[39], p_boc1, p_boc2;
@ -418,6 +445,49 @@ class PairReaxCKokkos : public PairReaxC {

  tdual_LR_lookup_table_kk_2d k_LR;
  t_LR_lookup_table_kk_2d d_LR;
+
+  DAT::tdual_int_2d k_tmpid;
+  DAT::tdual_ffloat_2d k_tmpbo;
+  DAT::tdual_int_scalar k_error_flag;
+
+  typename AT::t_int_1d d_numneigh_bonds;
+  typename AT::t_tagint_2d d_neighid;
+  typename AT::t_ffloat_2d d_abo;
+
+  typename AT::t_ffloat_1d d_buf;
+  DAT::tdual_int_scalar k_nbuf_local;
+};
+
+template <class DeviceType>
+struct PairReaxCKokkosFindBondFunctor  {
+  typedef DeviceType device_type;
+  typedef int value_type;
+  PairReaxCKokkos<DeviceType> c;
+  PairReaxCKokkosFindBondFunctor(PairReaxCKokkos<DeviceType>* c_ptr):c(*c_ptr) {};
+
+  KOKKOS_INLINE_FUNCTION
+  void join(volatile int &dst,
+             const volatile int &src) const {
+    dst = MAX(dst,src);
+  }
+
+  KOKKOS_INLINE_FUNCTION
+  void operator()(const int ii, int &numbonds) const {
+    c.calculate_find_bond_item(ii,numbonds);
+  }
+};
+
+template <class DeviceType>
+struct PairReaxCKokkosPackBondBufferFunctor  {
+  typedef DeviceType device_type;
+  typedef int value_type;
+  PairReaxCKokkos<DeviceType> c;
+  PairReaxCKokkosPackBondBufferFunctor(PairReaxCKokkos<DeviceType>* c_ptr):c(*c_ptr) {};
+
+  KOKKOS_INLINE_FUNCTION
+  void operator()(const int ii, int &j, const bool &final) const {
+    c.pack_bond_buffer_item(ii,j,final);
+  }
 };

 }
--- a/src/MOLECULE/bond_table.cpp
+++ b/src/MOLECULE/bond_table.cpp
@ -585,17 +585,26 @@ double BondTable::splint(double *xa, double *ya, double *y2a, int n, double x)

 /* ----------------------------------------------------------------------
   calculate potential u and force f at distance x
-   insure x is between bond min/max
+   insure x is between bond min/max, exit with error if not
 ------------------------------------------------------------------------- */

 void BondTable::uf_lookup(int type, double x, double &u, double &f)
 {
  int itable;
  double fraction,a,b;
+  char estr[128];

  Table *tb = &tables[tabindex[type]];
-  x = MAX(x,tb->lo);
-  x = MIN(x,tb->hi);
+  if (x < tb->lo) {
+    sprintf(estr,"Bond length < table inner cutoff: "
+            "type %d length %g",type,x);
+    error->one(FLERR,estr);
+  }
+  if (x > tb->hi) {
+    sprintf(estr,"Bond length > table outer cutoff: "
+            "type %d length %g",type,x);
+    error->one(FLERR,estr);
+  }

  if (tabstyle == LINEAR) {
    itable = static_cast<int> ((x - tb->lo) * tb->invdelta);
--- a/src/USER-COLVARS/colvarproxy_lammps.cpp
+++ b/src/USER-COLVARS/colvarproxy_lammps.cpp
@ -333,7 +333,10 @@ int colvarproxy_lammps::backup_file(char const *filename)

 int colvarproxy_lammps::smp_enabled()
 {
-  return COLVARS_OK;
+  if (b_smp_active) {
+    return COLVARS_OK;
+  }
+  return COLVARS_ERROR;
 }


--- a/src/USER-REAXC/compute_spec_atom.cpp
+++ b/src/USER-REAXC/compute_spec_atom.cpp
@ -44,6 +44,8 @@ ComputeSpecAtom::ComputeSpecAtom(LAMMPS *lmp, int narg, char **arg) :

  // Initiate reaxc
  reaxc = (PairReaxC *) force->pair_match("reax/c",1);
+  if (reaxc == NULL)
+    reaxc = (PairReaxC *) force->pair_match("reax/c/kk",1);

  pack_choice = new FnPtrPack[nvalues];

--- a/src/USER-REAXC/fix_reaxc_bonds.cpp
+++ b/src/USER-REAXC/fix_reaxc_bonds.cpp
@ -107,11 +107,10 @@ void FixReaxCBonds::setup(int vflag)

 void FixReaxCBonds::init()
 {
-  Pair *pair_kk = force->pair_match("reax/c/kk",1);
-  if (pair_kk != NULL) error->all(FLERR,"Cannot (yet) use fix reax/c/bonds with "
-                  "pair_style reax/c/kk");
-
  reaxc = (PairReaxC *) force->pair_match("reax/c",1);
+  if (reaxc == NULL)
+    reaxc = (PairReaxC *) force->pair_match("reax/c/kk",1);
+
  if (reaxc == NULL) error->all(FLERR,"Cannot use fix reax/c/bonds without "
                  "pair_style reax/c");

--- a/src/USER-REAXC/fix_reaxc_bonds.h
+++ b/src/USER-REAXC/fix_reaxc_bonds.h
@ -29,13 +29,13 @@ namespace LAMMPS_NS {
 class FixReaxCBonds : public Fix {
 public:
  FixReaxCBonds(class LAMMPS *, int, char **);
-  ~FixReaxCBonds();
+  virtual ~FixReaxCBonds();
  int setmask();
-  void init();
+  virtual void init();
  void setup(int);
  void end_of_step();

- private:
+ protected:
  int me, nprocs, nmax, ntypes, maxsize;
  int *numneigh;
  tagint **neighid;
@ -44,12 +44,12 @@ class FixReaxCBonds : public Fix {

  void allocate();
  void destroy();
-  void Output_ReaxC_Bonds(bigint, FILE *);
+  virtual void Output_ReaxC_Bonds(bigint, FILE *);
  void FindBond(struct _reax_list*, int &);
  void PassBuffer(double *, int &);
  void RecvBuffer(double *, int, int, int, int);
  int nint(const double &);
-  double memory_usage();
+  virtual double memory_usage();

  bigint nvalid, nextvalid();
  struct _reax_list *lists;
--- a/src/USER-REAXC/fix_reaxc_species.cpp
+++ b/src/USER-REAXC/fix_reaxc_species.cpp
@ -284,11 +284,10 @@ void FixReaxCSpecies::init()
  if (atom->tag_enable == 0)
    error->all(FLERR,"Cannot use fix reax/c/species unless atoms have IDs");

-  Pair *pair_kk = force->pair_match("reax/c/kk",1);
-  if (pair_kk != NULL) error->all(FLERR,"Cannot (yet) use fix reax/c/species with "
-                  "pair_style reax/c/kk");
-
  reaxc = (PairReaxC *) force->pair_match("reax/c",1);
+  if (reaxc == NULL)
+    reaxc = (PairReaxC *) force->pair_match("reax/c/kk",1);
+
  if (reaxc == NULL) error->all(FLERR,"Cannot use fix reax/c/species without "
 		  "pair_style reax/c");

@ -485,7 +484,7 @@ void FixReaxCSpecies::Output_ReaxC_Bonds(bigint ntimestep, FILE *fp)

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

-AtomCoord chAnchor(AtomCoord in1, AtomCoord in2)
+AtomCoord FixReaxCSpecies::chAnchor(AtomCoord in1, AtomCoord in2)
 {
  if (in1.x < in2.x)
    return in1;
--- a/src/USER-REAXC/fix_reaxc_species.h
+++ b/src/USER-REAXC/fix_reaxc_species.h
@ -38,15 +38,15 @@ typedef struct {
 class FixReaxCSpecies : public Fix {
 public:
  FixReaxCSpecies(class LAMMPS *, int, char **);
-  ~FixReaxCSpecies();
+  virtual ~FixReaxCSpecies();
  int setmask();
-  void init();
+  virtual void init();
  void init_list(int, class NeighList *);
  void setup(int);
  void post_integrate();
  double compute_vector(int);

- private:
+ protected:
  int me, nprocs, nmax, nlocal, ntypes, ntotal;
  int nrepeat, nfreq, posfreq;
  int Nmoltype, vector_nmole, vector_nspec;
@ -67,7 +67,8 @@ class FixReaxCSpecies : public Fix {
  void Output_ReaxC_Bonds(bigint, FILE *);
  void create_compute();
  void create_fix();
-  void FindMolecule();
+  AtomCoord chAnchor(AtomCoord, AtomCoord);
+  virtual void FindMolecule();
  void SortMolecule(int &);
  void FindSpecies(int, int &);
  void WriteFormulas(int, int);
--- a/src/atom.cpp
+++ b/src/atom.cpp
@ -26,11 +26,14 @@
 #include "force.h"
 #include "modify.h"
 #include "fix.h"
+#include "compute.h"
 #include "output.h"
 #include "thermo.h"
 #include "update.h"
 #include "domain.h"
 #include "group.h"
+#include "input.h"
+#include "variable.h"
 #include "molecule.h"
 #include "atom_masks.h"
 #include "math_const.h"
@ -1388,6 +1391,33 @@ void Atom::data_bodies(int n, char *buf, AtomVecBody *avec_body,
  delete [] dvalues;
 }

+/* ----------------------------------------------------------------------
+   init per-atom fix/compute/variable values for newly created atoms
+   called from create_atoms, read_data, read_dump,
+     lib::lammps_create_atoms()
+   fixes, computes, variables may or may not exist when called
+------------------------------------------------------------------------- */
+
+void Atom::data_fix_compute_variable(int nprev, int nnew)
+{
+  for (int m = 0; m < modify->nfix; m++) {
+    Fix *fix = modify->fix[m];
+    if (fix->create_attribute)
+      for (int i = nprev; i < nnew; i++)
+        fix->set_arrays(i);
+  }
+
+  for (int m = 0; m < modify->ncompute; m++) {
+    Compute *compute = modify->compute[m];
+    if (compute->create_attribute)
+      for (int i = nprev; i < nnew; i++)
+        compute->set_arrays(i);
+  }
+
+  for (int i = nprev; i < nnew; i++)
+    input->variable->set_arrays(i);
+}
+
 /* ----------------------------------------------------------------------
   allocate arrays of length ntypes
   only done after ntypes is set
--- a/src/atom.h
+++ b/src/atom.h
@ -225,15 +225,14 @@ class Atom : protected Pointers {

  void data_atoms(int, char *, tagint, int, int, double *);
  void data_vels(int, char *, tagint);
-
  void data_bonds(int, char *, int *, tagint, int);
  void data_angles(int, char *, int *, tagint, int);
  void data_dihedrals(int, char *, int *, tagint, int);
  void data_impropers(int, char *, int *, tagint, int);
-
  void data_bonus(int, char *, class AtomVec *, tagint);
  void data_bodies(int, char *, class AtomVecBody *, tagint);
-
+  void data_fix_compute_variable(int, int);
+  
  virtual void allocate_type_arrays();
  void set_mass(const char *, int);
  void set_mass(int, double);
--- a/src/create_atoms.cpp
+++ b/src/create_atoms.cpp
@ -359,30 +359,9 @@ void CreateAtoms::command(int narg, char **arg)
  else if (style == RANDOM) add_random();
  else add_lattice();

-  // invoke set_arrays() for fixes/computes/variables
-  //   that need initialization of attributes of new atoms
-  // don't use modify->create_attributes() since would be inefficient
-  //   for large number of atoms
-  // note that for typical early use of create_atoms,
-  //   no fixes/computes/variables exist yet
-
-  int nlocal = atom->nlocal;
-  for (int m = 0; m < modify->nfix; m++) {
-    Fix *fix = modify->fix[m];
-    if (fix->create_attribute)
-      for (int i = nlocal_previous; i < nlocal; i++)
-        fix->set_arrays(i);
-  }
-
-  for (int m = 0; m < modify->ncompute; m++) {
-    Compute *compute = modify->compute[m];
-    if (compute->create_attribute)
-      for (int i = nlocal_previous; i < nlocal; i++)
-        compute->set_arrays(i);
-  }
-
-  for (int i = nlocal_previous; i < nlocal; i++)
-    input->variable->set_arrays(i);
+  // init per-atom fix/compute/variable values for created atoms
+  
+  atom->data_fix_compute_variable(nlocal_previous,atom->nlocal);

  // set new total # of atoms and error check

--- a/src/create_box.cpp
+++ b/src/create_box.cpp
@ -159,6 +159,7 @@ void CreateBox::command(int narg, char **arg)
    } else if (strcmp(arg[iarg],"extra/special/per/atom") == 0) {
      if (iarg+2 > narg) error->all(FLERR,"Illegal create_box command");
      force->special_extra = force->inumeric(FLERR,arg[iarg+1]);
+      atom->maxspecial += force->special_extra;
      iarg += 2;
    } else error->all(FLERR,"Illegal create_box command");
  }
--- a/src/delete_atoms.cpp
+++ b/src/delete_atoms.cpp
@ -59,7 +59,9 @@ void DeleteAtoms::command(int narg, char **arg)
  bigint ndihedrals_previous = atom->ndihedrals;
  bigint nimpropers_previous = atom->nimpropers;

-  // delete the atoms
+  // flag atoms for deletion
+
+  allflag = 0;

  if (strcmp(arg[0],"group") == 0) delete_group(narg,arg);
  else if (strcmp(arg[0],"region") == 0) delete_region(narg,arg);
@ -67,36 +69,44 @@ void DeleteAtoms::command(int narg, char **arg)
  else if (strcmp(arg[0],"porosity") == 0) delete_porosity(narg,arg);
  else error->all(FLERR,"Illegal delete_atoms command");

-  // optionally delete additional bonds or atoms in molecules
+  // if allflag = 1, just reset atom->nlocal
+  // else delete atoms one by one

-  if (bond_flag) delete_bond();
-  if (mol_flag) delete_molecule();
+  if (allflag) atom->nlocal = 0;
+  else {

-  // delete local atoms flagged in dlist
-  // reset nlocal
+    // optionally delete additional bonds or atoms in molecules

-  AtomVec *avec = atom->avec;
-  int nlocal = atom->nlocal;
+    if (bond_flag) delete_bond();
+    if (mol_flag) delete_molecule();

-  int i = 0;
-  while (i < nlocal) {
-    if (dlist[i]) {
-      avec->copy(nlocal-1,i,1);
-      dlist[i] = dlist[nlocal-1];
-      nlocal--;
-    } else i++;
+    // delete local atoms flagged in dlist
+    // reset nlocal
+
+    AtomVec *avec = atom->avec;
+    int nlocal = atom->nlocal;
+
+    int i = 0;
+    while (i < nlocal) {
+      if (dlist[i]) {
+	avec->copy(nlocal-1,i,1);
+	dlist[i] = dlist[nlocal-1];
+	nlocal--;
+      } else i++;
+    }
+    
+    atom->nlocal = nlocal;
+    memory->destroy(dlist);
  }
-
-  atom->nlocal = nlocal;
-  memory->destroy(dlist);
-
+  
  // if non-molecular system and compress flag set,
  // reset atom tags to be contiguous
  // set all atom IDs to 0, call tag_extend()

  if (atom->molecular == 0 && compress_flag) {
    tagint *tag = atom->tag;
-    for (i = 0; i < nlocal; i++) tag[i] = 0;
+    int nlocal = atom->nlocal;
+    for (int i = 0; i < nlocal; i++) tag[i] = 0;
    atom->tag_extend();
  }

@ -185,6 +195,13 @@ void DeleteAtoms::delete_group(int narg, char **arg)
  if (igroup == -1) error->all(FLERR,"Could not find delete_atoms group ID");
  options(narg-2,&arg[2]);

+  // check for special case of group = all
+
+  if (strcmp(arg[1],"all") == 0) {
+    allflag = 1;
+    return;
+  }
+  
  // allocate and initialize deletion list

  int nlocal = atom->nlocal;
--- a/src/delete_atoms.h
+++ b/src/delete_atoms.h
@ -32,7 +32,7 @@ class DeleteAtoms : protected Pointers {

 private:
  int *dlist;
-  int compress_flag,bond_flag,mol_flag;
+  int allflag,compress_flag,bond_flag,mol_flag;
  std::map<tagint,int> *hash;

  void delete_group(int, char **);
--- a/src/domain.cpp
+++ b/src/domain.cpp
@ -1496,6 +1496,28 @@ void Domain::image_flip(int m, int n, int p)
  }
 }

+/* ----------------------------------------------------------------------
+   return 1 if this proc owns atom with coords x, else return 0
+   x is returned remapped into periodic box
+------------------------------------------------------------------------- */
+
+int Domain::ownatom(double *x)
+{
+  double lamda[3];
+  double *coord;
+  
+  remap(x);
+  if (triclinic) {
+    x2lamda(x,lamda);
+    coord = lamda;
+  } else coord = x;
+  
+  if (coord[0] >= sublo[0] && coord[0] < subhi[0] &&
+      coord[1] >= sublo[1] && coord[1] < subhi[1] &&
+      coord[2] >= sublo[2] && coord[2] < subhi[2]) return 1;
+  return 0;
+}
+
 /* ----------------------------------------------------------------------
   create a lattice
 ------------------------------------------------------------------------- */
@ -1929,5 +1951,3 @@ void Domain::lamda_box_corners(double *lo, double *hi)
  corners[7][0] = hi[0]; corners[7][1] = hi[1]; corners[7][2] = subhi_lamda[2];
  lamda2x(corners[7],corners[7]);
 }
-
-
--- a/src/domain.h
+++ b/src/domain.h
@ -121,6 +121,8 @@ class Domain : protected Pointers {
  void unmap(double *, imageint);
  void unmap(double *, imageint, double *);
  void image_flip(int, int, int);
+  int ownatom(double *);
+  
  void set_lattice(int, char **);
  void add_region(int, char **);
  void delete_region(int, char **);
--- a/src/fix_halt.h
+++ b/src/fix_halt.h
@ -36,7 +36,7 @@ class FixHalt : public Fix {

 private:
  int attribute,operation,eflag,ivar;
-  double attvalue,value;
+  double value;
  char *idvar;

  double bondmax();
--- a/src/fix_nh.cpp
+++ b/src/fix_nh.cpp
@ -51,11 +51,10 @@ enum{ISO,ANISO,TRICLINIC};
 ---------------------------------------------------------------------- */

 FixNH::FixNH(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg),
-tstat_flag(0), pstat_flag(0),
 rfix(NULL), id_dilate(NULL), irregular(NULL), id_temp(NULL), id_press(NULL),
-tcomputeflag(0), pcomputeflag(0), eta(NULL), eta_dot(NULL), eta_dotdot(NULL),
+eta(NULL), eta_dot(NULL), eta_dotdot(NULL),
 eta_mass(NULL), etap(NULL), etap_dot(NULL), etap_dotdot(NULL), 
-etap_mass(NULL), mpchain(0)
+etap_mass(NULL)
 {
  if (narg < 4) error->all(FLERR,"Illegal fix nvt/npt/nph command");

@ -478,8 +477,10 @@ etap_mass(NULL), mpchain(0)

    // pre_exchange only required if flips can occur due to shape changes

-    if (flipflag && (p_flag[3] || p_flag[4] || p_flag[5])) pre_exchange_flag = 1;
-    if (flipflag && (domain->yz != 0.0 || domain->xz != 0.0 || domain->xy != 0.0))
+    if (flipflag && (p_flag[3] || p_flag[4] || p_flag[5])) 
+      pre_exchange_flag = 1;
+    if (flipflag && (domain->yz != 0.0 || domain->xz != 0.0 || 
+                     domain->xy != 0.0))
      pre_exchange_flag = 1;
  }

--- a/src/input.cpp
+++ b/src/input.cpp
@ -278,7 +278,8 @@ void Input::file(const char *filename)
 }

 /* ----------------------------------------------------------------------
-   copy command in single to line, parse and execute it
+   invoke one command in single
+   first copy to line, then parse, then execute it
   return command name to caller
 ------------------------------------------------------------------------- */

--- a/src/library.cpp
+++ b/src/library.cpp
@ -18,9 +18,11 @@
 #include <string.h>
 #include <stdlib.h>
 #include "library.h"
+#include "lmptype.h"
 #include "lammps.h"
 #include "universe.h"
 #include "input.h"
+#include "atom_vec.h"
 #include "atom.h"
 #include "domain.h"
 #include "update.h"
@ -38,8 +40,12 @@

 using namespace LAMMPS_NS;

+// ----------------------------------------------------------------------
+// utility macros
+// ----------------------------------------------------------------------
+
 /* ----------------------------------------------------------------------
-   Utility macros for optional code path which captures all exceptions
+   macros for optional code path which captures all exceptions
   and stores the last error message. These assume there is a variable lmp
   which is a pointer to the current LAMMPS instance.

@ -76,6 +82,37 @@ using namespace LAMMPS_NS;
 #define END_CAPTURE
 #endif

+// ----------------------------------------------------------------------
+// helper functions, not in library API
+// ----------------------------------------------------------------------
+
+/* ----------------------------------------------------------------------
+   concatenate one or more LAMMPS input lines starting at ptr
+   removes NULL terminator when last printable char of line = '&'
+     by replacing both NULL and '&' with space character
+   repeat as many times as needed
+   on return, ptr now points to longer line
+------------------------------------------------------------------------- */
+
+void concatenate_lines(char *ptr)
+{
+  int nend = strlen(ptr);
+  int n = nend-1;
+  while (n && isspace(ptr[n])) n--;
+  while (ptr[n] == '&') {
+    ptr[nend] = ' ';
+    ptr[n] = ' ';
+    strtok(ptr,"\n");
+    nend = strlen(ptr);
+    n = nend-1;
+    while (n && isspace(ptr[n])) n--;
+  }
+}
+
+// ----------------------------------------------------------------------
+// library API functions to create/destroy an instance of LAMMPS
+//   and communicate commands to it
+// ----------------------------------------------------------------------

 /* ----------------------------------------------------------------------
   create an instance of LAMMPS and return pointer to it
@ -172,12 +209,14 @@ void lammps_file(void *ptr, char *str)

 /* ----------------------------------------------------------------------
   process a single input command in str
+   does not matter if str ends in newline
+   return command name to caller
 ------------------------------------------------------------------------- */

 char *lammps_command(void *ptr, char *str)
 {
  LAMMPS *lmp = (LAMMPS *) ptr;
-  char * result = NULL;
+  char *result = NULL;

  BEGIN_CAPTURE
  {
@ -188,6 +227,69 @@ char *lammps_command(void *ptr, char *str)
  return result;
 }

+/* ----------------------------------------------------------------------
+   process multiple input commands in cmds = list of strings
+   does not matter if each string ends in newline
+   create long contatentated string for processing by commands_string()
+   insert newlines in concatenated string as needed
+------------------------------------------------------------------------- */
+
+void lammps_commands_list(void *ptr, int ncmd, char **cmds)
+{
+  LAMMPS *lmp = (LAMMPS *) ptr;
+
+  int n = ncmd+1;
+  for (int i = 0; i < ncmd; i++) n += strlen(cmds[i]);
+
+  char *str = (char *) lmp->memory->smalloc(n,"lib/commands/list:str");
+  str[0] = '\0';
+  n = 0;
+
+  for (int i = 0; i < ncmd; i++) {
+    strcpy(&str[n],cmds[i]);
+    n += strlen(cmds[i]);
+    if (str[n-1] != '\n') {
+      str[n] = '\n';
+      str[n+1] = '\0';
+      n++;
+    }
+  }
+
+  lammps_commands_string(ptr,str);
+  lmp->memory->sfree(str);
+}
+
+/* ----------------------------------------------------------------------
+   process multiple input commands in single long str, separated by newlines
+   single command can span multiple lines via continuation characters 
+   multi-line commands enabled by triple quotes will not work
+------------------------------------------------------------------------- */
+
+void lammps_commands_string(void *ptr, char *str)
+{
+  LAMMPS *lmp = (LAMMPS *) ptr;
+
+  // make copy of str so can strtok() it
+
+  int n = strlen(str) + 1;
+  char *copy = new char[n];
+  strcpy(copy,str);
+
+  BEGIN_CAPTURE
+  {
+    char *ptr = strtok(copy,"\n");
+    if (ptr) concatenate_lines(ptr);
+    while (ptr) {
+      lmp->input->one(ptr);
+      ptr = strtok(NULL,"\n");
+      if (ptr) concatenate_lines(ptr);
+    }
+  }
+  END_CAPTURE
+
+  delete [] copy;
+}
+
 /* ----------------------------------------------------------------------
   clean-up function to free memory allocated by lib and returned to caller
 ------------------------------------------------------------------------- */
@ -197,6 +299,10 @@ void lammps_free(void *ptr)
  free(ptr);
 }

+// ----------------------------------------------------------------------
+// library API functions to extract info from LAMMPS or set info in LAMMPS
+// ----------------------------------------------------------------------
+
 /* ----------------------------------------------------------------------
   add LAMMPS-specific library functions
   all must receive LAMMPS pointer as argument
@ -209,6 +315,8 @@ void lammps_free(void *ptr)
   returns a void pointer to the entity
     which the caller can cast to the proper data type
   returns a NULL if name not listed below
+   this function need only be invoked once
+     the returned pointer is a permanent valid reference to the quantity
   customize by adding names
 ------------------------------------------------------------------------- */

@ -232,14 +340,16 @@ void *lammps_extract_global(void *ptr, char *name)
  if (strcmp(name,"ndihedrals") == 0) return (void *) &lmp->atom->ndihedrals;
  if (strcmp(name,"nimpropers") == 0) return (void *) &lmp->atom->nimpropers;
  if (strcmp(name,"nlocal") == 0) return (void *) &lmp->atom->nlocal;
+  if (strcmp(name,"nghost") == 0) return (void *) &lmp->atom->nghost;
+  if (strcmp(name,"nmax") == 0) return (void *) &lmp->atom->nmax;
  if (strcmp(name,"ntimestep") == 0) return (void *) &lmp->update->ntimestep;

-  // NOTE: we cannot give access to the thermo "time" data by reference,
-  // as that is a recomputed property.  Only "atime" can be provided as pointer.
-  // please use lammps_get_thermo() defined below to access all supported
-  // thermo keywords by value.
+  // update atime can be referenced as a pointer
+  // thermo "timer" data cannot be, since it is computed on request
+  // lammps_get_thermo() can access all thermo keywords by value

  if (strcmp(name,"atime") == 0) return (void *) &lmp->update->atime;
+  if (strcmp(name,"atimestep") == 0) return (void *) &lmp->update->atimestep;

  return NULL;
 }
@ -250,6 +360,8 @@ void *lammps_extract_global(void *ptr, char *name)
   returns a void pointer to the entity
     which the caller can cast to the proper data type
   returns a NULL if Atom::extract() does not recognize the name
+   the returned pointer is not a permanent valid reference to the
+     per-atom quantity, since LAMMPS may reallocate per-atom data
   customize by adding names to Atom::extract()
 ------------------------------------------------------------------------- */

@ -261,6 +373,7 @@ void *lammps_extract_atom(void *ptr, char *name)

 /* ----------------------------------------------------------------------
   extract a pointer to an internal LAMMPS compute-based entity
+   the compute is invoked if its value(s) is not current
   id = compute ID
   style = 0 for global data, 1 for per-atom data, 2 for local data
   type = 0 for scalar, 1 for vector, 2 for array
@ -275,6 +388,8 @@ void *lammps_extract_atom(void *ptr, char *name)
   returns a void pointer to the compute's internal data structure
     for the entity which the caller can cast to the proper data type
   returns a NULL if id is not recognized or style/type not supported
+   the returned pointer is not a permanent valid reference to the
+     compute data, this function should be re-invoked
   IMPORTANT: if the compute is not current it will be invoked
     LAMMPS cannot easily check here if it is valid to invoke the compute,
     so caller must insure that it is OK
@ -492,11 +607,11 @@ int lammps_set_variable(void *ptr, char *name, char *str)
 }

 /* ----------------------------------------------------------------------
-   return the current value of a thermo keyword as double.
+   return the current value of a thermo keyword as a double
   unlike lammps_extract_global() this does not give access to the
-   storage of the data in question, and thus needs to be called
-   again to retrieve an updated value. The upshot is that it allows
-   accessing information that is only computed on-the-fly.
+     storage of the data in question
+   instead it triggers the Thermo class to compute the current value
+     and returns it
 ------------------------------------------------------------------------- */

 double lammps_get_thermo(void *ptr, char *name)
@ -529,12 +644,13 @@ int lammps_get_natoms(void *ptr)

 /* ----------------------------------------------------------------------
   gather the named atom-based entity across all processors
+   atom IDs must be consecutive from 1 to N
   name = desired quantity, e.g. x or charge
   type = 0 for integer values, 1 for double values
   count = # of per-atom values, e.g. 1 for type or charge, 3 for x or f
   return atom-based values in data, ordered by count, then by atom ID
     e.g. x[0][0],x[0][1],x[0][2],x[1][0],x[1][1],x[1][2],x[2][0],...
-   data must be pre-allocated by caller to correct length
+     data must be pre-allocated by caller to correct length
 ------------------------------------------------------------------------- */

 void lammps_gather_atoms(void *ptr, char *name,
@ -547,7 +663,8 @@ void lammps_gather_atoms(void *ptr, char *name,
    // error if tags are not defined or not consecutive

    int flag = 0;
-    if (lmp->atom->tag_enable == 0 || lmp->atom->tag_consecutive() == 0) flag = 1;
+    if (lmp->atom->tag_enable == 0 || lmp->atom->tag_consecutive() == 0) 
+      flag = 1;
    if (lmp->atom->natoms > MAXSMALLINT) flag = 1;
    if (flag) {
      if (lmp->comm->me == 0)
@ -586,7 +703,7 @@ void lammps_gather_atoms(void *ptr, char *name,
          for (j = 0; j < count; j++)
            copy[offset++] = array[i][0];
        }
-
+      
      MPI_Allreduce(copy,data,count*natoms,MPI_INT,MPI_SUM,lmp->world);
      lmp->memory->destroy(copy);

@ -623,6 +740,7 @@ void lammps_gather_atoms(void *ptr, char *name,

 /* ----------------------------------------------------------------------
   scatter the named atom-based entity across all processors
+   atom IDs must be consecutive from 1 to N
   name = desired quantity, e.g. x or charge
   type = 0 for integer values, 1 for double values
   count = # of per-atom values, e.g. 1 for type or charge, 3 for x or f
@ -640,7 +758,8 @@ void lammps_scatter_atoms(void *ptr, char *name,
    // error if tags are not defined or not consecutive or no atom map

    int flag = 0;
-    if (lmp->atom->tag_enable == 0 || lmp->atom->tag_consecutive() == 0) flag = 1;
+    if (lmp->atom->tag_enable == 0 || lmp->atom->tag_consecutive() == 0) 
+      flag = 1;
    if (lmp->atom->natoms > MAXSMALLINT) flag = 1;
    if (lmp->atom->map_style == 0) flag = 1;
    if (flag) {
@ -702,9 +821,91 @@ void lammps_scatter_atoms(void *ptr, char *name,
  END_CAPTURE
 }

-#ifdef LAMMPS_EXCEPTIONS
 /* ----------------------------------------------------------------------
-   Check if a new error message
+   create N atoms and assign them to procs based on coords
+   id = atom IDs (optional, NULL if just use 1 to N)
+   type = N-length vector of atom types (required)
+   x = 3N-length vector of atom coords (required)
+   v = 3N-length vector of atom velocities (optional, NULL if just 0.0)
+   x,v = ordered by xyz, then by atom
+     e.g. x[0][0],x[0][1],x[0][2],x[1][0],x[1][1],x[1][2],x[2][0],...
+------------------------------------------------------------------------- */
+
+void lammps_create_atoms(void *ptr, int n, tagint *id, int *type,
+			 double *x, double *v)
+{
+  LAMMPS *lmp = (LAMMPS *) ptr;
+
+  BEGIN_CAPTURE
+  {
+    // error if box does not exist or tags not defined
+
+    int flag = 0;
+    if (lmp->domain->box_exist == 0) flag = 1;
+    if (lmp->atom->tag_enable == 0) flag = 1;
+    if (flag) {
+      if (lmp->comm->me == 0)
+        lmp->error->warning(FLERR,"Library error in lammps_create_atoms");
+      return;
+    }
+
+    // loop over N atoms of entire system
+    // if this proc owns it based on coords, invoke create_atom()
+    // optionally set atom tags and velocities
+
+    Atom *atom = lmp->atom;
+    Domain *domain = lmp->domain;
+    int nlocal = atom->nlocal;
+
+    int nprev = nlocal;
+    double xdata[3];
+    
+    for (int i = 0; i < n; i++) {
+      xdata[0] = x[3*i];
+      xdata[1] = x[3*i+1];
+      xdata[2] = x[3*i+2];
+      if (!domain->ownatom(xdata)) continue;
+  
+      atom->avec->create_atom(type[i],xdata);
+      if (id) atom->tag[nlocal] = id[i];
+      else atom->tag[nlocal] = i+1;
+      if (v) {
+	atom->v[nlocal][0] = v[3*i];
+	atom->v[nlocal][1] = v[3*i+1];
+	atom->v[nlocal][2] = v[3*i+2];
+      }
+      nlocal++;
+    }
+
+    // need to reset atom->natoms inside LAMMPS
+
+    bigint ncurrent = nlocal;
+    MPI_Allreduce(&ncurrent,&lmp->atom->natoms,1,MPI_LMP_BIGINT,
+                  MPI_SUM,lmp->world);
+
+    // init per-atom fix/compute/variable values for created atoms
+
+    atom->data_fix_compute_variable(nprev,nlocal);
+     
+    // if global map exists, reset it
+    // invoke map_init() b/c atom count has grown
+
+    if (lmp->atom->map_style) {
+      lmp->atom->map_init();
+      lmp->atom->map_set();
+    }
+  }
+  END_CAPTURE
+}
+
+// ----------------------------------------------------------------------
+// library API functions for error handling
+// ----------------------------------------------------------------------
+
+#ifdef LAMMPS_EXCEPTIONS
+
+/* ----------------------------------------------------------------------
+   check if a new error message
 ------------------------------------------------------------------------- */

 int lammps_has_error(void *ptr) {
@ -714,8 +915,8 @@ int lammps_has_error(void *ptr) {
 }

 /* ----------------------------------------------------------------------
-   Copy the last error message of LAMMPS into a character buffer
-   The return value encodes which type of error it is.
+   copy the last error message of LAMMPS into a character buffer
+   return value encodes which type of error:
   1 = normal error (recoverable)
   2 = abort error (non-recoverable)
 ------------------------------------------------------------------------- */
@ -732,4 +933,5 @@ int lammps_get_last_error_message(void *ptr, char * buffer, int buffer_size) {
  }
  return 0;
 }
+
 #endif
--- a/src/library.h
+++ b/src/library.h
@ -30,6 +30,8 @@ void lammps_close(void *);
 int  lammps_version(void *);
 void lammps_file(void *, char *);
 char *lammps_command(void *, char *);
+void lammps_commands_list(void *, int, char **);
+void lammps_commands_string(void *, char *);
 void lammps_free(void *);

 void *lammps_extract_global(void *, char *);
@ -40,10 +42,11 @@ void *lammps_extract_variable(void *, char *, char *);

 int lammps_set_variable(void *, char *, char *);
 double lammps_get_thermo(void *, char *);
-int lammps_get_natoms(void *);

+int lammps_get_natoms(void *);
 void lammps_gather_atoms(void *, char *, int, int, void *);
 void lammps_scatter_atoms(void *, char *, int, int, void *);
+void lammps_create_atoms(void *, int, int *, int *, double *, double *);

 #ifdef LAMMPS_EXCEPTIONS
 int lammps_has_error(void *);
--- a/src/molecule.cpp
+++ b/src/molecule.cpp
@ -595,6 +595,10 @@ void Molecule::read(int flag)
  // set maxspecial on first pass, so allocate() has a size

  if (bondflag && specialflag == 0) {
+    if (domain->box_exist == 0)
+      error->all(FLERR,"Cannot auto-generate special bonds before "
+                       "simulation box is defined");
+
    maxspecial = atom->maxspecial;
    if (flag) {
      special_generate();
--- a/src/pair_table.cpp
+++ b/src/pair_table.cpp
@ -554,7 +554,7 @@ void PairTable::spline_table(Table *tb)

 /* ----------------------------------------------------------------------
   extract attributes from parameter line in table section
-   format of line: N value R/RSQ/BITMAP lo hi FP fplo fphi
+   format of line: N value R/RSQ/BITMAP lo hi FPRIME fplo fphi
   N is required, other params are optional
 ------------------------------------------------------------------------- */

@ -578,7 +578,7 @@ void PairTable::param_extract(Table *tb, char *line)
      tb->rlo = atof(word);
      word = strtok(NULL," \t\n\r\f");
      tb->rhi = atof(word);
-    } else if (strcmp(word,"FP") == 0) {
+    } else if (strcmp(word,"FPRIME") == 0) {
      tb->fpflag = 1;
      word = strtok(NULL," \t\n\r\f");
      tb->fplo = atof(word);
--- a/src/read_data.cpp
+++ b/src/read_data.cpp
@ -703,7 +703,11 @@ void ReadData::command(int narg, char **arg)
    if (addflag == NONE) atom->deallocate_topology();
    atom->avec->grow(atom->nmax);
  }
+
+  // init per-atom fix/compute/variable values for created atoms
  
+  atom->data_fix_compute_variable(nlocal_previous,atom->nlocal);
+
  // assign atoms added by this data file to specified group

  if (groupbit) {
@ -830,6 +834,7 @@ void ReadData::command(int narg, char **arg)
  }

  // restore old styles, when reading with nocoeff flag given
+  
  if (coeffflag == 0) {
    if (force->pair) delete force->pair;
    force->pair = saved_pair;
--- a/src/read_dump.cpp
+++ b/src/read_dump.cpp
@ -908,27 +908,9 @@ void ReadDump::process_atoms(int n)
    }
  }

-  // invoke set_arrays() for fixes/computes/variables
-  //   that need initialization of attributes of new atoms
-  // same as in CreateAtoms
-  // don't use modify->create_attributes() since would be inefficient
-  //   for large number of atoms
-
-  nlocal = atom->nlocal;
-  for (int m = 0; m < modify->nfix; m++) {
-    Fix *fix = modify->fix[m];
-    if (fix->create_attribute)
-      for (i = nlocal_previous; i < nlocal; i++)
-        fix->set_arrays(i);
-  }
-  for (int m = 0; m < modify->ncompute; m++) {
-    Compute *compute = modify->compute[m];
-    if (compute->create_attribute)
-      for (i = nlocal_previous; i < nlocal; i++)
-        compute->set_arrays(i);
-  }
-  for (int i = nlocal_previous; i < nlocal; i++)
-    input->variable->set_arrays(i);
+  // init per-atom fix/compute/variable values for created atoms
+  
+  atom->data_fix_compute_variable(nlocal_previous,atom->nlocal);
 }

 /* ----------------------------------------------------------------------
--- a/src/region.cpp
+++ b/src/region.cpp
@ -486,7 +486,7 @@ void Region::set_velocity()
    else v[0] = v[1] = v[2] = 0.0;
    prev[0] = dx;
    prev[1] = dy;
-    prev[2] = dz;    
+    prev[2] = dz;
  }

  if (rotateflag) {
@ -546,14 +546,13 @@ void Region::velocity_contact(double *vwall, double *x, int ic)

 void Region::length_restart_string(int &n)
 {
-  n += sizeof(int) + strlen(id)+1 + 
+  n += sizeof(int) + strlen(id)+1 +
    sizeof(int) + strlen(style)+1 + sizeof(int) +
    size_restart*sizeof(double);
 }

 /* ----------------------------------------------------------------------
-   region writes its current style, id, number of sub-regions 
-     and position/angle
+   region writes its current style, id, number of sub-regions, position/angle
   needed by fix/wall/gran/region to compute velocity by differencing scheme
 ------------------------------------------------------------------------- */

@ -562,50 +561,36 @@ void Region::write_restart(FILE *fp)
  int sizeid = (strlen(id)+1);
  int sizestyle = (strlen(style)+1);
  fwrite(&sizeid, sizeof(int), 1, fp);
-  fwrite(id, 1, sizeid, fp);
-  fwrite(&sizestyle, sizeof(int), 1, fp);
-  fwrite(style, 1, sizestyle, fp);  
+  fwrite(id,1,sizeid,fp);
+  fwrite(&sizestyle,sizeof(int),1,fp);
+  fwrite(style,1,sizestyle,fp);
  fwrite(&nregion,sizeof(int),1,fp);
-
-  fwrite(prev, sizeof(double), size_restart, fp);  
+  fwrite(prev,sizeof(double),size_restart,fp);
 }

 /* ----------------------------------------------------------------------
   region reads style, id, number of sub-regions from restart file
-     if they match current region, also read previous position/angle
+   if they match current region, also read previous position/angle
   needed by fix/wall/gran/region to compute velocity by differencing scheme
 ------------------------------------------------------------------------- */

 int Region::restart(char *buf, int &n)
 {
-  int sizeid = buf[n];
+  int size = *((int *) (&buf[n]));
  n += sizeof(int);
-  char *restart_id = new char[sizeid];
-  for (int i = 0; i < sizeid; i++)
-    restart_id[i] = buf[n++];    
-  if (strcmp(restart_id,id) != 0) return 0;
+  if ((size <= 0) || (strcmp(&buf[n],id) != 0)) return 0;
+  n += size;

-  int sizestyle = buf[n];
+  size = *((int *) (&buf[n]));
  n += sizeof(int);
-  char *restart_style = new char[sizestyle];
-  for (int i = 0; i < sizestyle; i++)
-    restart_style[i] = buf[n++];  
-  if (strcmp(restart_style,style) != 0) return 0;    
+  if ((size <= 0) || (strcmp(&buf[n],style) != 0)) return 0;
+  n += size;

-  int restart_nreg = buf[n];
+  int restart_nreg = *((int *) (&buf[n]));
  n += sizeof(int);
  if (restart_nreg != nregion) return 0;

-  char *rlist = new char[size_restart*sizeof(double)];  
-  for (int i = 0; i < size_restart*sizeof(double); i++)
-    rlist[i] = buf[n++]; 
-  for (int i = 0; i < size_restart; i++){
-    prev[i] = ((double *)rlist)[i];
-  }
-  
-  delete [] rlist;
-  delete [] restart_id;
-  delete [] restart_style;
+  memcpy(prev,&buf[n],size_restart*sizeof(double));
  return 1;
 }

--- a/src/region.h
+++ b/src/region.h
@ -60,7 +60,7 @@ class Region : protected Pointers {
  double omega[3];	      // angular velocity
  double rprev;               // speed of time-dependent radius, if applicable
  double xcenter[3];          // translated/rotated center of cylinder/sphere (only used if varshape)
-  double prev[5];             // stores displacement (X3), angle and if 
+  double prev[5];             // stores displacement (X3), angle and if
                              //  necessary, region variable size (e.g. radius)
                              //  at previous time step
  int vel_timestep;           // store timestep at which set_velocity was called
@ -106,11 +106,12 @@ class Region : protected Pointers {
  void options(int, char **);
  void point_on_line_segment(double *, double *, double *, double *);
  void forward_transform(double &, double &, double &);
+  double point[3],runit[3];

 private:
  char *xstr,*ystr,*zstr,*tstr;
  int xvar,yvar,zvar,tvar;
-  double axis[3],point[3],runit[3];
+  double axis[3];

  void inverse_transform(double &, double &, double &);
  void rotate(double &, double &, double &, double);
--- a/src/region_block.cpp
+++ b/src/region_block.cpp
@ -344,7 +344,7 @@ int RegBlock::surface_exterior(double *x, double cutoff)
  store closest point in xc,yc,zc
 --------------------------------------------------------------------------*/

-double RegBlock::find_closest_point(int i, double *x, 
+double RegBlock::find_closest_point(int i, double *x,
                                    double &xc, double &yc, double &zc)
 {
  double dot,d2,d2min;
@ -372,7 +372,7 @@ double RegBlock::find_closest_point(int i, double *x,

  } else {
    point_on_line_segment(corners[i][0],corners[i][1],x,p);
-    d2 = (p[0]-x[0])*(p[0]-x[0]) + (p[1]-x[1])*(p[1]-x[1]) + 
+    d2 = (p[0]-x[0])*(p[0]-x[0]) + (p[1]-x[1])*(p[1]-x[1]) +
      (p[2]-x[2])*(p[2]-x[2]);
    if (d2 < d2min) {
      d2min = d2;
@ -382,7 +382,7 @@ double RegBlock::find_closest_point(int i, double *x,
    }

    point_on_line_segment(corners[i][1],corners[i][2],x,p);
-    d2 = (p[0]-x[0])*(p[0]-x[0]) + (p[1]-x[1])*(p[1]-x[1]) + 
+    d2 = (p[0]-x[0])*(p[0]-x[0]) + (p[1]-x[1])*(p[1]-x[1]) +
      (p[2]-x[2])*(p[2]-x[2]);
    if (d2 < d2min) {
      d2min = d2;
@ -392,7 +392,7 @@ double RegBlock::find_closest_point(int i, double *x,
    }

    point_on_line_segment(corners[i][2],corners[i][3],x,p);
-    d2 = (p[0]-x[0])*(p[0]-x[0]) + (p[1]-x[1])*(p[1]-x[1]) + 
+    d2 = (p[0]-x[0])*(p[0]-x[0]) + (p[1]-x[1])*(p[1]-x[1]) +
      (p[2]-x[2])*(p[2]-x[2]);
    if (d2 < d2min) {
      d2min = d2;
@ -402,7 +402,7 @@ double RegBlock::find_closest_point(int i, double *x,
    }

    point_on_line_segment(corners[i][3],corners[i][4],x,p);
-    d2 = (p[0]-x[0])*(p[0]-x[0]) + (p[1]-x[1])*(p[1]-x[1]) + 
+    d2 = (p[0]-x[0])*(p[0]-x[0]) + (p[1]-x[1])*(p[1]-x[1]) +
      (p[2]-x[2])*(p[2]-x[2]);
    if (d2 < d2min) {
      d2min = d2;
--- a/src/region_cone.cpp
+++ b/src/region_cone.cpp
@ -416,7 +416,7 @@ int RegCone::surface_exterior(double *x, double cutoff)
    // x is far enough from cone that there is no contact
    // x is interior to cone

-    if (r >= maxradius+cutoff || x[0] <= lo-cutoff || x[0] >= hi+cutoff) 
+    if (r >= maxradius+cutoff || x[0] <= lo-cutoff || x[0] >= hi+cutoff)
      return 0;
    if (r < currentradius && x[0] > lo && x[0] < hi) return 0;

--- a/src/region_cylinder.cpp
+++ b/src/region_cylinder.cpp
@ -36,7 +36,7 @@ RegCylinder::RegCylinder(LAMMPS *lmp, int narg, char **arg) :

  // check open face settings

-  if (openflag && (open_faces[2] || open_faces[3] || 
+  if (openflag && (open_faces[2] || open_faces[3] ||
                   open_faces[4] || open_faces[5]))
    error->all(FLERR,"Invalid region cylinder open setting");

@ -691,16 +691,16 @@ void RegCylinder::variable_check()


 /* ----------------------------------------------------------------------
-   Set values needed to calculate velocity due to shape changes. 
+   Set values needed to calculate velocity due to shape changes.
   These values do not depend on the contact, so this function is
   called once per timestep by fix/wall/gran/region.
-   
+
 ------------------------------------------------------------------------- */

 void RegCylinder::set_velocity_shape()
 {
  if (axis == 'x'){
-    xcenter[0] = 0; 
+    xcenter[0] = 0;
    xcenter[1] = c1;
    xcenter[2] = c2;
  }
@ -715,15 +715,15 @@ void RegCylinder::set_velocity_shape()
    xcenter[2] = 0;
  }
  forward_transform(xcenter[0], xcenter[1], xcenter[2]);
-  if (update->ntimestep > 0) rprev = prev[4];  
+  if (update->ntimestep > 0) rprev = prev[4];
  else rprev = radius;
-  prev[4] = radius;     
+  prev[4] = radius;
 }



 /* ----------------------------------------------------------------------
-   add velocity due to shape change to wall velocity 
+   add velocity due to shape change to wall velocity
 ------------------------------------------------------------------------- */

 void RegCylinder::velocity_contact_shape(double *vwall, double *xc)
@ -746,7 +746,7 @@ void RegCylinder::velocity_contact_shape(double *vwall, double *xc)
  }
  vwall[0] += delx/update->dt;
  vwall[1] += dely/update->dt;
-  vwall[2] += delz/update->dt;     
+  vwall[2] += delz/update->dt;
  //printf ("R is %g, prev %g, velocity of wall at %g %g %g is %g %g %g\n",radius,rprev,xc[0],xc[1],xc[2],vwall[0],vwall[1],vwall[2]);
 }

--- a/src/region_intersect.cpp
+++ b/src/region_intersect.cpp
@ -43,7 +43,7 @@ RegIntersect::RegIntersect(LAMMPS *lmp, int narg, char **arg) :
    idsub[nregion] = new char[m];
    strcpy(idsub[nregion],arg[iarg+3]);
    iregion = domain->find_region(idsub[nregion]);
-    if (iregion == -1) 
+    if (iregion == -1)
      error->all(FLERR,"Region intersect region ID does not exist");
    list[nregion++] = iregion;
  }
@ -124,7 +124,7 @@ void RegIntersect::init()
  int iregion;
  for (int ilist = 0; ilist < nregion; ilist++) {
    iregion = domain->find_region(idsub[ilist]);
-    if (iregion == -1) 
+    if (iregion == -1)
      error->all(FLERR,"Region union region ID does not exist");
    list[ilist] = iregion;
  }
@ -290,7 +290,7 @@ void RegIntersect::set_velocity()

 void RegIntersect::length_restart_string(int& n)
 {
-  n += sizeof(int) + strlen(id)+1 + 
+  n += sizeof(int) + strlen(id)+1 +
    sizeof(int) + strlen(style)+1 + sizeof(int);
  for (int ilist = 0; ilist < nregion; ilist++)
    domain->regions[list[ilist]]->length_restart_string(n);
@ -308,7 +308,7 @@ void RegIntersect::write_restart(FILE *fp)
  fwrite(&sizeid, sizeof(int), 1, fp);
  fwrite(id, 1, sizeid, fp);
  fwrite(&sizestyle, sizeof(int), 1, fp);
-  fwrite(style, 1, sizestyle, fp);  
+  fwrite(style, 1, sizestyle, fp);
  fwrite(&nregion,sizeof(int),1,fp);

  for (int ilist = 0; ilist < nregion; ilist++){
@ -321,33 +321,26 @@ void RegIntersect::write_restart(FILE *fp)
   needed by fix/wall/gran/region to compute velocity by differencing scheme
 ------------------------------------------------------------------------- */

-int RegIntersect::restart(char *buf, int& n)
+int RegIntersect::restart(char *buf, int &n)
 {
-  int sizeid = buf[n];
+  int size = *((int *) (&buf[n]));
  n += sizeof(int);
-  char *restart_id = new char[sizeid];
-  for (int i = 0; i < sizeid; i++)
-    restart_id[i] = buf[n++];    
-  if (strcmp(restart_id,id) != 0) return 0;
+  if ((size <= 0) || (strcmp(&buf[n],id) != 0)) return 0;
+  n += size;

-  int sizestyle = buf[n];
+  size = *((int *) (&buf[n]));
  n += sizeof(int);
+  if ((size <= 0) || (strcmp(&buf[n],style) != 0)) return 0;
+  n += size;

-  char *restart_style = new char[sizestyle];
-  for (int i = 0; i < sizestyle; i++)
-    restart_style[i] = buf[n++];  
-  if (strcmp(restart_style,style) != 0) return 0;    
-
-  int restart_nreg = buf[n];
+  int restart_nreg = *((int *) (&buf[n]));
  n += sizeof(int);
  if (restart_nreg != nregion) return 0;

-  for (int ilist = 0; ilist < nregion; ilist++){
-    if (!domain->regions[list[ilist]]->restart(buf, n)){
-      return 0;
-    }
-  }
-  return 1; 
+  for (int ilist = 0; ilist < nregion; ilist++)
+    if (!domain->regions[list[ilist]]->restart(buf,n)) return 0;
+
+  return 1;
 }

 /* ----------------------------------------------------------------------
--- a/src/region_sphere.cpp
+++ b/src/region_sphere.cpp
@ -188,10 +188,10 @@ void RegSphere::variable_check()


 /* ----------------------------------------------------------------------
-   Set values needed to calculate velocity due to shape changes. 
+   Set values needed to calculate velocity due to shape changes.
   These values do not depend on the contact, so this function is
   called once per timestep by fix/wall/gran/region.
-   
+
 ------------------------------------------------------------------------- */

 void RegSphere::set_velocity_shape()
@ -200,27 +200,27 @@ void RegSphere::set_velocity_shape()
  xcenter[1] = yc;
  xcenter[2] = zc;
  forward_transform(xcenter[0], xcenter[1], xcenter[2]);
-  if (update->ntimestep > 0) rprev = prev[4];  
+  if (update->ntimestep > 0) rprev = prev[4];
  else rprev = radius;
-  prev[4] = radius;     
+  prev[4] = radius;
 }



 /* ----------------------------------------------------------------------
-   add velocity due to shape change to wall velocity 
+   add velocity due to shape change to wall velocity
 ------------------------------------------------------------------------- */

 void RegSphere::velocity_contact_shape(double *vwall, double *xc)
 {
  double delx, dely, delz; // Displacement of contact point in x,y,z
- 
+
  delx = (xc[0] - xcenter[0])*(1 - rprev/radius);
  dely = (xc[1] - xcenter[1])*(1 - rprev/radius);
  delz = (xc[2] - xcenter[2])*(1 - rprev/radius);

  vwall[0] += delx/update->dt;
  vwall[1] += dely/update->dt;
-  vwall[2] += delz/update->dt;     
+  vwall[2] += delz/update->dt;
 }

--- a/src/region_union.cpp
+++ b/src/region_union.cpp
@ -44,7 +44,7 @@ RegUnion::RegUnion(LAMMPS *lmp, int narg, char **arg) : Region(lmp, narg, arg)
    idsub[nregion] = new char[m];
    strcpy(idsub[nregion],arg[iarg+3]);
    iregion = domain->find_region(idsub[nregion]);
-    if (iregion == -1) 
+    if (iregion == -1)
      error->all(FLERR,"Region union region ID does not exist");
    list[nregion++] = iregion;
  }
@ -118,7 +118,7 @@ void RegUnion::init()
  int iregion;
  for (int ilist = 0; ilist < nregion; ilist++) {
    iregion = domain->find_region(idsub[ilist]);
-    if (iregion == -1) 
+    if (iregion == -1)
      error->all(FLERR,"Region union region ID does not exist");
    list[ilist] = iregion;
  }
@ -171,7 +171,7 @@ int RegUnion::surface_interior(double *x, double cutoff)
      for (jlist = 0; jlist < nregion; jlist++) {
        if (jlist == ilist) continue;
        jregion = list[jlist];
-        if (regions[jregion]->match(xs,ys,zs) && 
+        if (regions[jregion]->match(xs,ys,zs) &&
            !regions[jregion]->openflag) break;
      }
      if (jlist == nregion) {
@ -284,7 +284,7 @@ void RegUnion::set_velocity()

 void RegUnion::length_restart_string(int& n)
 {
-  n += sizeof(int) + strlen(id)+1 + 
+  n += sizeof(int) + strlen(id)+1 +
    sizeof(int) + strlen(style)+1 + sizeof(int);
  for (int ilist = 0; ilist < nregion; ilist++)
    domain->regions[list[ilist]]->length_restart_string(n);
@ -302,10 +302,10 @@ void RegUnion::write_restart(FILE *fp)
  fwrite(&sizeid, sizeof(int), 1, fp);
  fwrite(id, 1, sizeid, fp);
  fwrite(&sizestyle, sizeof(int), 1, fp);
-  fwrite(style, 1, sizestyle, fp);  
+  fwrite(style, 1, sizestyle, fp);
  fwrite(&nregion,sizeof(int),1,fp);
  for (int ilist = 0; ilist < nregion; ilist++)
-    domain->regions[list[ilist]]->write_restart(fp);  
+    domain->regions[list[ilist]]->write_restart(fp);
 }

 /* ----------------------------------------------------------------------
@ -315,27 +315,23 @@ void RegUnion::write_restart(FILE *fp)

 int RegUnion::restart(char *buf, int &n)
 {
-  int sizeid = buf[n];
+  int size = *((int *) (&buf[n]));
  n += sizeof(int);
-  char *restart_id = new char[sizeid];
-  for (int i = 0; i < sizeid; i++)
-    restart_id[i] = buf[n++];    
-  if (strcmp(restart_id,id) != 0) return 0;
+  if ((size <= 0) || (strcmp(&buf[n],id) != 0)) return 0;
+  n += size;

-  int sizestyle = buf[n];
+  size = *((int *) (&buf[n]));
  n += sizeof(int);
-  char *restart_style = new char[sizestyle];
-  for (int i = 0; i < sizestyle; i++)
-    restart_style[i] = buf[n++];  
-  if (strcmp(restart_style,style) != 0) return 0;    
+  if ((size <= 0) || (strcmp(&buf[n],style) != 0)) return 0;
+  n += size;

-  int restart_nreg = buf[n];
+  int restart_nreg = *((int *) (&buf[n]));
  n += sizeof(int);
  if (restart_nreg != nregion) return 0;

-  for (int ilist = 0; ilist < nregion; ilist++){
-    if (!domain->regions[list[ilist]]->restart(buf, n)) return 0;
-  }
+  for (int ilist = 0; ilist < nregion; ilist++)
+    if (!domain->regions[list[ilist]]->restart(buf,n)) return 0;
+
  return 1;
 }

--- a/src/verlet.cpp
+++ b/src/verlet.cpp
@ -95,6 +95,9 @@ void Verlet::setup()
    timer->print_timeout(screen);
  }

+  if (lmp->kokkos)
+    error->all(FLERR,"KOKKOS package requires run_style verlet/kk");
+
  update->setupflag = 1;

  // setup domain, communication and neighboring
--- a/Show More
+++ b/Show More
Author	SHA1	Message	Date
Steve Plimpton	d9891abdf4	new library functions	2016-10-27 09:34:04 -06:00
sjplimp	efaa8feab5	Merge pull request #239 from akohlmey/static-analysis-fixes Static analysis fixes	2016-10-25 10:32:55 -06:00
sjplimp	ad5f7c4581	Merge pull request #238 from giacomofiorin/colvars-2016-10-24 Colvars fixes and small changes	2016-10-25 10:32:41 -06:00
sjplimp	6b33499135	Merge pull request #231 from akohlmey/collected-doc-fixes Collected doc fixes	2016-10-25 10:30:34 -06:00
Axel Kohlmeyer	63eada2425	fix issue with docs for orientorder/atom compute reported by @andeplane this closes #243	2016-10-25 12:12:48 -04:00
Axel Kohlmeyer	1a436bd7a9	Merge branch 'collected-doc-fixes' of github.com:akohlmey/lammps into collected-doc-fixes	2016-10-25 11:58:42 -04:00
Axel Kohlmeyer	52dd9aee5f	Merge branch 'master' into collected-doc-fixes	2016-10-25 11:55:09 -04:00
Steve Plimpton	eca96e21ef	Merge branch 'doc'	2016-10-25 09:46:07 -06:00
Steve Plimpton	9c81ad1ab6	doc page changes	2016-10-25 09:45:55 -06:00
Axel Kohlmeyer	f8367e3d0f	update documentation pdf for updated colvars lib	2016-10-24 17:49:53 -04:00
Axel Kohlmeyer	ba6d1528bb	Merge branch 'colvars-2016-10-24' of https://github.com/giacomofiorin/lammps into colvars-update	2016-10-24 17:34:28 -04:00
Giacomo Fiorin	182141b850	Make SMP parallelism for Colvars optional	2016-10-24 17:13:34 -04:00
Axel Kohlmeyer	512c413b7e	whitespace cleanup	2016-10-24 17:13:21 -04:00
Axel Kohlmeyer	7b89e47a38	apply corrections to issues reported by static code analysis	2016-10-24 17:12:28 -04:00
Giacomo Fiorin	e02505c8cc	Add ensemble-biased metadynamics (Fabrizio Marinelli, NIH)	2016-10-24 17:11:09 -04:00
Giacomo Fiorin	be2d155cef	Minor changes and fixes not relevant to LAMMPS	2016-10-24 17:10:52 -04:00
Giacomo Fiorin	c243093980	Fix wall forces and subtractAppliedForce for extended-Lagrangian ABF	2016-10-24 17:05:47 -04:00
Giacomo Fiorin	ad57a17f48	Add C-linkage wrapper for colvarscript (useful with ctypes)	2016-10-24 16:48:20 -04:00
sjplimp	477ddaf112	Merge pull request #232 from akohlmey/small-bugfixes Small bugfixes	2016-10-24 08:15:08 -06:00
sjplimp	4f69d91a99	Merge pull request #230 from akohlmey/manual-in-ebook-format generate LAMMPS manual in ebook format	2016-10-24 08:12:08 -06:00
Axel Kohlmeyer	bc44988003	correct typo in write_dump docs this closes #233	2016-10-23 15:18:25 -04:00
Axel Kohlmeyer	db36c8bcc3	stop with error, if molecule command requires special bond auto-generation before box is defined	2016-10-21 14:51:09 -04:00
Axel Kohlmeyer	991034b632	have bond style table exit when bond length is outside table range	2016-10-21 14:01:06 -04:00
Axel Kohlmeyer	607246f923	ignore mobi file as well	2016-10-21 13:25:53 -04:00
Axel Kohlmeyer	6742fb634a	remove mobi file format creation from makefile and explain it in README instead	2016-10-21 12:05:21 -04:00
Axel Kohlmeyer	ed3f02f249	ignore generated PDF and ePUB files	2016-10-21 12:04:48 -04:00
Axel Kohlmeyer	a2e34aab0a	make certain, that atom->maxspecial is incremented with extra special space	2016-10-21 11:55:36 -04:00
Axel Kohlmeyer	6cd6c106ef	Merge branch 'collected-small-changes' into collected-doc-fixes	2016-10-20 19:27:18 -04:00
Axel Kohlmeyer	a9572275ee	Revert "support generation of manual in ePUB format" This reverts commit `8c3f5cb307`.	2016-10-20 16:27:00 -04:00
Axel Kohlmeyer	2cf77ff778	Add support for ebook generation in ePUB and mobi format	2016-10-20 16:16:17 -04:00
Axel Kohlmeyer	f022f6d88a	fix various formatting and broken link issues identified by ebook-convert	2016-10-20 14:40:18 -04:00
Axel Kohlmeyer	8c3f5cb307	support generation of manual in ePUB format	2016-10-20 09:27:26 -04:00
Axel Kohlmeyer	e8359923f1	update packages section in manual with information about USER-NC-DUMP	2016-10-19 15:58:50 -04:00
Steve Plimpton	9954d5d346	forgot pair table change	2016-10-19 10:47:07 -06:00