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
129 changed files with 4717 additions and 484 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/commands.txt
+++ b/doc/src/commands.txt
@ -37,6 +37,7 @@ Commands :h1
   dump_image
   dump_modify
   dump_molfile
+   dump_nc
   echo
   fix
   fix_modify
--- 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
@ -0,0 +1,66 @@
+"LAMMPS WWW Site"_lws - "LAMMPS Documentation"_ld - "LAMMPS Commands"_lc :c
+
+:link(lws,http://lammps.sandia.gov)
+:link(ld,Manual.html)
+:link(lc,Section_commands.html#comm)
+
+:line
+
+dump nc command :h3
+dump nc/mpiio command :h3
+
+[Syntax:]
+
+dump ID group-ID nc N file.nc args
+dump ID group-ID nc/mpiio N file.nc args :pre
+
+ID = user-assigned name for the dump :ulb,l
+group-ID = ID of the group of atoms to be imaged :l
+{nc} or {nc/mpiio}  = style of dump command (other styles {atom} or {cfg} or {dcd} or {xtc} or {xyz} or {local} or {custom} are discussed on the "dump"_dump.html doc page) :l
+N = dump every this many timesteps :l
+file.nc = name of file to write to :l
+args = list of per atom data elements to dump, same as for the 'custom' dump style. :l,ule
+
+[Examples:]
+
+dump 1 all nc 100 traj.nc type x y z vx vy vz
+dump_modify 1 append yes at -1 global c_thermo_pe c_thermo_temp c_thermo_press :pre
+
+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.
+
+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.
+
+:line
+
+[Related commands:]
+
+"dump"_dump.html, "dump_modify"_dump_modify.html, "undump"_undump.html
+
--- 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/fixes.txt
+++ b/doc/src/fixes.txt
@ -48,6 +48,7 @@ Fixes :h1
   fix_gld
   fix_gle
   fix_gravity
+   fix_halt
   fix_heat
   fix_imd
   fix_indent
--- a/doc/src/lammps.book
+++ b/doc/src/lammps.book
@ -170,6 +170,7 @@ fix_gcmc.html
 fix_gld.html
 fix_gle.html
 fix_gravity.html
+fix_halt.html
 fix_heat.html
 fix_imd.html
 fix_indent.html
--- 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_gauss.txt
+++ b/doc/src/pair_gauss.txt
@ -106,8 +106,31 @@ more instructions on how to use the accelerated styles effectively.

 [Mixing, shift, table, tail correction, restart, rRESPA info]:

-These pair styles do not support mixing. Thus, coefficients for all
-I,J pairs must be specified explicitly.
+For atom type pairs I,J and I != J, the A, B, H, sigma_h, r_mh
+parameters, and the cutoff distance for these pair styles can be mixed:
+A (energy units)
+sqrt(1/B) (distance units, see below)
+H (energy units)
+sigma_h (distance units)
+r_mh (distance units)
+cutoff (distance units):ul
+
+The default mix value is {geometric}.
+Only {arithmetic} and {geometric} mix values are supported.
+See the "pair_modify" command for details.
+
+The A and H parameters are mixed using the same rules normally
+used to mix the "epsilon" parameter in a Lennard Jones interaction.
+The sigma_h, r_mh, and the cutoff distance are mixed using the same
+rules used to mix the "sigma" parameter in a Lennard Jones interaction.
+The B parameter is converted to a distance (sigma), before mixing
+(using sigma=B^-0.5), and converted back to a coefficient
+afterwards (using B=sigma^2).
+Negative A values are converted to positive A values (using abs(A))
+before mixing, and converted back after mixing
+(by multiplying by sign(Ai)*sign(Aj)).
+This way, if either particle is repulsive (if Ai<0 or Aj<0),
+then the default interaction between both particles will be repulsive.

 The {gauss} style does not support the "pair_modify"_pair_modify.html
 shift option. There is no effect due to the Gaussian well beyond the
--- 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/lib/netcdf/Makefile.lammps
+++ b/lib/netcdf/Makefile.lammps
@ -0,0 +1,19 @@
+# Settings that the LAMMPS build will import when a package using the
+# netCDF library is installed. This tries to automate configuration
+# via the nc-config tool, which is part of the netCDF installation.
+
+netcdf_SYSINC =
+netcdf_SYSLIB =
+netcdf_SYSPATH =
+
+ifneq ($(shell which nc-config 2>> /dev/null),)
+netcdf_SYSINC += -DLMP_HAS_NETCDF $(shell nc-config --cflags)
+netcdf_SYSLIB += $(shell nc-config --libs)
+endif
+
+ifneq ($(shell which ncmpidump 2>> /dev/null),)
+PNETCDF_PATH = $(shell which ncmpidump | sed -e 's,bin/ncmpidump,,')
+netcdf_SYSINC += -DLMP_HAS_PNETCDF -I$(PNETCDF_PATH)/include
+netcdf_SYSLIB += -lpnetcdf
+netcdf_SYSPATH += -L$(PNETCDF_PATH)/lib
+endif
--- a/lib/netcdf/README
+++ b/lib/netcdf/README
@ -0,0 +1,43 @@
+The Makefile.lammps file in this directory is used when building
+LAMMPS with packages that make use of the NetCDF library or its
+parallel version.  The file has several settings needed to compile
+and link LAMMPS with the NetCDF and parallel NetCDF support.
+For any regular NetCDF installation, all required flags should be
+autodetected. Please note that parallel NetCDF support is
+beneficial only when you run on a machine with very many processors
+like an IBM BlueGene or Cray. For most people regular NetCDF
+support should be sufficient and not cause any performance
+penalties.
+
+If you have problems compiling or linking, you may have to set
+the flags manually. There are three makefile variables
+
+1) netcdf_SYSINC
+This is for setting preprocessor options and include file paths.
+Set -DLMP_HAS_NETCDF, if you have NetCDF installed.
+Set -DLMP_HAS_PNETCDF, if you have parallel NetCDF installed.
+You can have either or both defines set. If none of these are
+set, LAMMPS will compile, but the NetCDF enabled functionality
+will not be available.
+In addition you may have to point to the folder with the include
+with -I/path/to/netcdf/include
+
+Example for a Fedora 24 machine with serial NetCDF installed as
+netcdf-devel-4.4.0-3.fc24.x86_64 RPM package:
+
+netcdf_SYSINC = -DLMP_HAS_NETCDF -I/usr/include -I/usr/include/hdf
+
+2) netcdf_SYSLIB
+This is the setting for all required libraries that need to be linked to.
+
+Example for a Fedora 24 machine with serial NetCDF installed as
+netcdf-devel-4.4.0-3.fc24.x86_64 RPM package:
+
+netcdf_SYSLIB = -lnetcdf
+
+3) netcdf_SYSPATH
+This is the setting for the path of directories with the NetCDF libraries.
+Typically, this will be of the form -L/path/to/netcdf/lib
+
+In the example from above, it can be left empty, because the Linux
+distribution provided libraries are installed in a system library location.
--- 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/.gitignore
+++ b/src/.gitignore
@ -271,6 +271,10 @@
 /dump_custom_vtk.h
 /dump_h5md.cpp
 /dump_h5md.h
+/dump_nc.cpp
+/dump_nc.h
+/dump_nc_mpiio.cpp
+/dump_nc_mpiio.h
 /dump_xtc.cpp
 /dump_xtc.h
 /dump_xyz_mpiio.cpp
--- 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/Makefile
+++ b/src/Makefile
@ -50,14 +50,14 @@ PACKAGE = asphere body class2 colloid compress coreshell dipole gpu \

 PACKUSER = user-atc user-awpmd user-cg-cmm user-colvars \
 	   user-diffraction user-dpd user-drude user-eff user-fep user-h5md \
-	   user-intel user-lb user-manifold user-mgpt \
-	   user-misc user-molfile user-omp user-phonon user-qmmm user-qtb \
+	   user-intel user-lb user-manifold user-mgpt user-misc user-molfile \
+	   user-nc-dump user-omp user-phonon user-qmmm user-qtb \
 	   user-quip user-reaxc user-smd user-smtbq user-sph user-tally \
 	   user-vtk

 PACKLIB = compress gpu kim kokkos meam mpiio poems python reax voronoi \
 	  user-atc user-awpmd user-colvars user-h5md user-lb user-molfile \
-	  user-qmmm user-quip user-vtk
+	  user-nc-dump user-qmmm user-quip user-smd user-vtk

 PACKALL = $(PACKAGE) $(PACKUSER)

--- a/src/REPLICA/prd.cpp
+++ b/src/REPLICA/prd.cpp
@ -795,7 +795,7 @@ void PRD::replicate(int ireplica)

  if (cmode == SINGLE_PROC_DIRECT) {
    MPI_Bcast(atom->x[0],3*atom->nlocal,MPI_DOUBLE,ireplica,comm_replica);
-    MPI_Bcast(atom->image,atom->nlocal,MPI_INT,ireplica,comm_replica);
+    MPI_Bcast(atom->image,atom->nlocal,MPI_LMP_IMAGEINT,ireplica,comm_replica);
    return;
  }

@ -815,9 +815,9 @@ void PRD::replicate(int ireplica)
      memcpy(imageall,image,nlocal*sizeof(imageint));
    }

-    MPI_Bcast(tagall,natoms,MPI_INT,ireplica,comm_replica);
+    MPI_Bcast(tagall,natoms,MPI_LMP_TAGINT,ireplica,comm_replica);
    MPI_Bcast(xall[0],3*natoms,MPI_DOUBLE,ireplica,comm_replica);
-    MPI_Bcast(imageall,natoms,MPI_INT,ireplica,comm_replica);
+    MPI_Bcast(imageall,natoms,MPI_LMP_IMAGEINT,ireplica,comm_replica);

    for (i = 0; i < nlocal; i++) {
      m = atom->map(tagall[i]);
@ -844,8 +844,8 @@ void PRD::replicate(int ireplica)
      displacements[i+1] = displacements[i] + counts[i];
    MPI_Gatherv(atom->tag,atom->nlocal,MPI_LMP_TAGINT,
                tagall,counts,displacements,MPI_LMP_TAGINT,0,world);
-    MPI_Gatherv(atom->image,atom->nlocal,MPI_LMP_TAGINT,
-                imageall,counts,displacements,MPI_LMP_TAGINT,0,world);
+    MPI_Gatherv(atom->image,atom->nlocal,MPI_LMP_IMAGEINT,
+                imageall,counts,displacements,MPI_LMP_IMAGEINT,0,world);
    for (i = 0; i < nprocs; i++) counts[i] *= 3;
    for (i = 0; i < nprocs-1; i++)
      displacements[i+1] = displacements[i] + counts[i];
@ -854,13 +854,13 @@ void PRD::replicate(int ireplica)
  }
  
  if (me == 0) {
-    MPI_Bcast(tagall,natoms,MPI_INT,ireplica,comm_replica);
-    MPI_Bcast(imageall,natoms,MPI_INT,ireplica,comm_replica);
+    MPI_Bcast(tagall,natoms,MPI_LMP_TAGINT,ireplica,comm_replica);
+    MPI_Bcast(imageall,natoms,MPI_LMP_IMAGEINT,ireplica,comm_replica);
    MPI_Bcast(xall[0],3*natoms,MPI_DOUBLE,ireplica,comm_replica);
  }
  
-  MPI_Bcast(tagall,natoms,MPI_INT,0,world);
-  MPI_Bcast(imageall,natoms,MPI_INT,0,world);
+  MPI_Bcast(tagall,natoms,MPI_LMP_TAGINT,0,world);
+  MPI_Bcast(imageall,natoms,MPI_LMP_IMAGEINT,0,world);
  MPI_Bcast(xall[0],3*natoms,MPI_DOUBLE,0,world);
  
  double **x = atom->x;
--- 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-MISC/pair_gauss_cut.cpp
+++ b/src/USER-MISC/pair_gauss_cut.cpp
@ -39,6 +39,7 @@ using namespace MathConst;
 PairGaussCut::PairGaussCut(LAMMPS *lmp) : Pair(lmp)
 {
  respa_enable = 0;
+  writedata = 1;
 }

 /* ---------------------------------------------------------------------- */
@ -96,7 +97,6 @@ void PairGaussCut::compute(int eflag, int vflag)

    for (jj = 0; jj < jnum; jj++) {
      j = jlist[jj];
-
      factor_lj = special_lj[sbmask(j)];
      j &= NEIGHMASK;

@ -221,7 +221,21 @@ void PairGaussCut::coeff(int narg, char **arg)

 double PairGaussCut::init_one(int i, int j)
 {
-  if (setflag[i][j] == 0) error->all(FLERR,"All pair coeffs are not set");
+  if (setflag[i][j] == 0) {
+    hgauss[i][j] = mix_energy(fabs(hgauss[i][i]), fabs(hgauss[j][j]),
+			      fabs(sigmah[i][i]), fabs(sigmah[j][j]));
+
+    // If either of the particles is repulsive (ie, if hgauss > 0),
+    // then the interaction between both is repulsive.
+    double sign_hi = (hgauss[i][i] >= 0.0) ? 1.0 : -1.0;
+    double sign_hj = (hgauss[j][j] >= 0.0) ? 1.0 : -1.0;
+    hgauss[i][j] *= MAX(sign_hi, sign_hj);
+
+    sigmah[i][j] = mix_distance(sigmah[i][i], sigmah[j][j]);
+    rmh[i][j] = mix_distance(rmh[i][i], rmh[j][j]);
+    cut[i][j] = mix_distance(cut[i][i], cut[j][j]);
+  }
+
  pgauss[i][j] = hgauss[i][j] / sqrt(MY_2PI) / sigmah[i][j];

  if (offset_flag) {
@ -334,6 +348,27 @@ void PairGaussCut::read_restart_settings(FILE *fp)
  MPI_Bcast(&mix_flag,1,MPI_INT,0,world);
 }

+/* ----------------------------------------------------------------------
+   proc 0 writes to data file
+------------------------------------------------------------------------- */
+
+void PairGaussCut::write_data(FILE *fp)
+{
+  for (int i = 1; i <= atom->ntypes; i++)
+    fprintf(fp,"%d %g %g %g\n",i,hgauss[i][i],rmh[i][i],sigmah[i][i]);
+}
+
+/* ----------------------------------------------------------------------
+   proc 0 writes all pairs to data file
+------------------------------------------------------------------------- */
+
+void PairGaussCut::write_data_all(FILE *fp)
+{
+  for (int i = 1; i <= atom->ntypes; i++)
+    for (int j = i; j <= atom->ntypes; j++)
+      fprintf(fp,"%d %d %g %g %g %g\n",i,j,hgauss[i][j],rmh[i][j],sigmah[i][j],cut[i][j]);
+}
+
 /* ---------------------------------------------------------------------- */

 double PairGaussCut::single(int i, int j, int itype, int jtype, double rsq,
--- a/src/USER-MISC/pair_gauss_cut.h
+++ b/src/USER-MISC/pair_gauss_cut.h
@ -42,6 +42,8 @@ class PairGaussCut : public Pair {
  virtual void read_restart(FILE *);
  virtual void write_restart_settings(FILE *);
  virtual void read_restart_settings(FILE *);
+  virtual void write_data(FILE *fp);
+  virtual void write_data_all(FILE *fp);

  virtual double memory_usage();

--- a/src/USER-NC-DUMP/Install.sh
+++ b/src/USER-NC-DUMP/Install.sh
@ -0,0 +1,63 @@
+# Install/unInstall package files in LAMMPS
+# mode = 0/1/2 for uninstall/install/update
+
+mode=$1
+
+# enforce using portable C locale
+LC_ALL=C
+export LC_ALL
+
+# arg1 = file, arg2 = file it depends on
+
+action () {
+  if (test $mode = 0) then
+    rm -f ../$1
+  elif (! cmp -s $1 ../$1) then
+    if (test -z "$2" || test -e ../$2) then
+      cp $1 ..
+      if (test $mode = 2) then
+        echo "  updating src/$1"
+      fi
+    fi
+  elif (test -n "$2") then
+    if (test ! -e ../$2) then
+      rm -f ../$1
+    fi
+  fi
+}
+
+for file in *.cpp *.h; do
+  action $file
+done
+
+# edit 2 Makefile.package files to include/exclude package info
+
+if (test $1 = 1) then
+
+  if (test -e ../Makefile.package) then
+    sed -i -e 's/[^ \t]*netcdf[^ \t]* //g' ../Makefile.package
+    sed -i -e 's|^PKG_SYSINC =[ \t]*|&$(netcdf_SYSINC) |' ../Makefile.package
+    sed -i -e 's|^PKG_SYSLIB =[ \t]*|&$(netcdf_SYSLIB) |' ../Makefile.package
+    sed -i -e 's|^PKG_SYSPATH =[ \t]*|&$(netcdf_SYSPATH) |' ../Makefile.package
+  fi
+
+  if (test -e ../Makefile.package.settings) then
+    sed -i -e '/^include.*netcdf.*$/d' ../Makefile.package.settings
+    # multiline form needed for BSD sed on Macs
+    sed -i -e '4 i \
+include ..\/..\/lib\/netcdf\/Makefile.lammps
+' ../Makefile.package.settings
+
+  fi
+
+elif (test $1 = 0) then
+
+  if (test -e ../Makefile.package) then
+    sed -i -e 's/[^ \t]*netcdf[^ \t]* //g' ../Makefile.package
+  fi
+
+  if (test -e ../Makefile.package.settings) then
+    sed -i -e '/^include.*netcdf.*$/d' ../Makefile.package.settings
+  fi
+
+fi
--- a/src/USER-NC-DUMP/README
+++ b/src/USER-NC-DUMP/README
@ -0,0 +1,39 @@
+USER-NC-DUMP
+============
+
+This package provides the nc and (optionally) the nc/mpiio dump styles.
+See the doc page for dump nc or dump nc/mpiio command for how to use them.
+Compiling these dump styles requires having the netCDF library installed
+on your system. See lib/netcdf/README for additional details.
+
+PACKAGE DESCRIPTION
+-------------------
+
+This is a LAMMPS (http://lammps.sandia.gov/) dump style for output into a NetCDF
+database. The database format follows the AMBER NetCDF trajectory convention
+(http://ambermd.org/netcdf/nctraj.xhtml), but includes extensions to this
+convention. These extension are:
+* A variable "cell_origin" (of dimension "frame", "cell_spatial") that contains
+  the bottom left corner of the simulation cell.
+* Any number of additional variables corresponding to per atom scalar, vector
+  or tensor quantities available within LAMMPS. Tensor quantities are written in
+  Voigt notation. An additional dimension "Voigt" of length 6 is created for
+  this purpose.
+* Possibility to output to an HDF5 database.
+
+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.
+* VMD (http://www.ks.uiuc.edu/Research/vmd/).
+* AtomEye (http://www.libatoms.org/). The libAtoms version of AtomEye contains
+  a NetCDF reader that is not present in the standard distribution of AtomEye.
+
+The person who created these files is Lars Pastewka at
+Karlsruhe Institute of Technology (lars.pastewka@kit.edu).
+Contact him directly if you have questions.
+
+Lars Pastewka
+Institute for Applied Materials (IAM)
+Karlsruhe Institute of Technology (KIT)
+Kaiserstrasse 12, 76131 Karlsruhe
+e-mail: lars.pastewka@kit.edu
--- a/src/USER-NC-DUMP/dump_nc.cpp
+++ b/src/USER-NC-DUMP/dump_nc.cpp
--- a/src/USER-NC-DUMP/dump_nc.h
+++ b/src/USER-NC-DUMP/dump_nc.h
@ -0,0 +1,143 @@
+/* ======================================================================
+   LAMMPS NetCDF dump style
+   https://github.com/pastewka/lammps-netcdf
+   Lars Pastewka, lars.pastewka@kit.edu
+
+   Copyright (2011-2013) Fraunhofer IWM
+   Copyright (2014) Karlsruhe Institute of Technology
+
+   This program is free software: you can redistribute it and/or modify
+   it under the terms of the GNU General Public License as published by
+   the Free Software Foundation, either version 2 of the License, or
+   (at your option) any later version.
+
+   This program is distributed in the hope that it will be useful,
+   but WITHOUT ANY WARRANTY; without even the implied warranty of
+   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+   GNU General Public License for more details.
+
+   You should have received a copy of the GNU General Public License
+   along with this program.  If not, see <http://www.gnu.org/licenses/>.
+   ====================================================================== */
+
+/* ----------------------------------------------------------------------
+   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.
+------------------------------------------------------------------------- */
+#if defined(LMP_HAS_NETCDF)
+
+#ifdef DUMP_CLASS
+
+DumpStyle(nc,DumpNC)
+
+#else
+
+#ifndef LMP_DUMP_NC_H
+#define LMP_DUMP_NC_H
+
+#include "dump_custom.h"
+
+namespace LAMMPS_NS {
+
+const int NC_FIELD_NAME_MAX = 100;
+const int DUMP_NC_MAX_DIMS  = 100;
+
+class DumpNC : public DumpCustom {
+ public:
+  DumpNC(class LAMMPS *, int, char **);
+  virtual ~DumpNC();
+  virtual void write();
+
+ private:
+  // per-atoms quantities (positions, velocities, etc.)
+  struct nc_perat_t {
+    int dims;                     // number of dimensions
+    int field[DUMP_NC_MAX_DIMS];  // field indices corresponding to the dim.
+    char name[NC_FIELD_NAME_MAX]; // field name
+    int var;                      // NetCDF variable
+
+    bool constant;                // is this property per file (not per frame)
+    int ndumped;                  // number of enties written for this prop.
+  };
+
+  typedef void (DumpNC::*funcptr_t)(void *);
+
+  // per-frame quantities (variables, fixes or computes)
+  struct nc_perframe_t {
+    char name[NC_FIELD_NAME_MAX]; // field name
+    int var;                      // NetCDF variable
+    int type;                     // variable, fix, compute or callback
+    int index;                    // index in fix/compute list
+    funcptr_t compute;            // compute function
+    int dim;                      // dimension
+    char id[NC_FIELD_NAME_MAX];   // variable id
+
+    bigint bigint_data;           // actual data
+    double double_data;           // actual data
+  };
+
+  int framei;                  // current frame index
+  int blocki;                  // current block index
+  int ndata;                   // number of data blocks to expect
+
+  bigint ntotalgr;             // # of atoms
+
+  int n_perat;                 // # of netcdf per-atom properties
+  nc_perat_t *perat;           // per-atom properties
+
+  int n_perframe;              // # of global netcdf (not per-atom) fix props
+  nc_perframe_t *perframe;     // global properties
+
+  bool double_precision;       // write everything as double precision
+
+  bigint n_buffer;             // size of buffer
+  int *int_buffer;             // buffer for passing data to netcdf
+  double *double_buffer;       // buffer for passing data to netcdf
+
+  int ncid;
+
+  int frame_dim;
+  int spatial_dim;
+  int Voigt_dim;
+  int atom_dim;
+  int cell_spatial_dim;
+  int cell_angular_dim;
+  int label_dim;
+
+  int spatial_var;
+  int cell_spatial_var;
+  int cell_angular_var;
+
+  int time_var;
+  int cell_origin_var;
+  int cell_lengths_var;
+  int cell_angles_var;
+
+  virtual void openfile();
+  void closefile();
+  virtual void write_header(bigint);
+  virtual void write_data(int, double *);
+  void write_prmtop();
+
+  virtual int modify_param(int, char **);
+
+  void ncerr(int, const char *, int);
+
+  void compute_step(void *);
+  void compute_elapsed(void *);
+  void compute_elapsed_long(void *);
+};
+
+}
+
+#endif
+#endif
+#endif /* defined(LMP_HAS_NETCDF) */
--- a/src/USER-NC-DUMP/dump_nc_mpiio.cpp
+++ b/src/USER-NC-DUMP/dump_nc_mpiio.cpp
--- a/src/USER-NC-DUMP/dump_nc_mpiio.h
+++ b/src/USER-NC-DUMP/dump_nc_mpiio.h
@ -0,0 +1,140 @@
+/* ======================================================================
+   LAMMPS NetCDF dump style
+   https://github.com/pastewka/lammps-netcdf
+   Lars Pastewka, lars.pastewka@kit.edu
+
+   Copyright (2011-2013) Fraunhofer IWM
+   Copyright (2014) Karlsruhe Institute of Technology
+
+   This program is free software: you can redistribute it and/or modify
+   it under the terms of the GNU General Public License as published by
+   the Free Software Foundation, either version 2 of the License, or
+   (at your option) any later version.
+
+   This program is distributed in the hope that it will be useful,
+   but WITHOUT ANY WARRANTY; without even the implied warranty of
+   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+   GNU General Public License for more details.
+
+   You should have received a copy of the GNU General Public License
+   along with this program.  If not, see <http://www.gnu.org/licenses/>.
+   ====================================================================== */
+
+/* ----------------------------------------------------------------------
+   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.
+------------------------------------------------------------------------- */
+#if defined(LMP_HAS_PNETCDF)
+
+#ifdef DUMP_CLASS
+
+DumpStyle(nc/mpiio,DumpNCMPIIO)
+
+#else
+
+#ifndef LMP_DUMP_NC_MPIIO_H
+#define LMP_DUMP_NC_MPIIO_H
+
+#include "dump_custom.h"
+
+namespace LAMMPS_NS {
+
+const int NC_MPIIO_FIELD_NAME_MAX = 100;
+const int DUMP_NC_MPIIO_MAX_DIMS  = 100;
+
+class DumpNCMPIIO : public DumpCustom {
+ public:
+  DumpNCMPIIO(class LAMMPS *, int, char **);
+  virtual ~DumpNCMPIIO();
+  virtual void write();
+
+ private:
+  // per-atoms quantities (positions, velocities, etc.)
+  struct nc_perat_t {
+    int dims;                           // number of dimensions
+    int field[DUMP_NC_MPIIO_MAX_DIMS];  // field indices corresponding to the dim.
+    char name[NC_MPIIO_FIELD_NAME_MAX]; // field name
+    int var;                            // NetCDF variable
+  };
+
+  typedef void (DumpNCMPIIO::*funcptr_t)(void *);
+
+  // per-frame quantities (variables, fixes or computes)
+  struct nc_perframe_t {
+    char name[NC_MPIIO_FIELD_NAME_MAX]; // field name
+    int var;                            // NetCDF variable
+    int type;                           // variable, fix, compute or callback
+    int index;                          // index in fix/compute list
+    funcptr_t compute;                  // compute function
+    int dim;                            // dimension
+    char id[NC_MPIIO_FIELD_NAME_MAX];   // variable id
+
+    bigint bigint_data;                 // actual data
+    double double_data;                 // actual data
+  };
+
+  int framei;                           // current frame index
+  int blocki;                           // current block index
+  int ndata;                            // number of data blocks to expect
+
+  bigint ntotalgr;                      // # of atoms
+
+  int n_perat;                          // # of netcdf per-atom properties
+  nc_perat_t *perat;                    // per-atom properties
+
+  int n_perframe;                       // # of global netcdf (not per-atom) fix props
+  nc_perframe_t *perframe;              // global properties
+
+  bool double_precision;                // write everything as double precision
+
+  bigint n_buffer;                      // size of buffer
+  int *int_buffer;                      // buffer for passing data to netcdf
+  double *double_buffer;                // buffer for passing data to netcdf
+
+  int ncid;
+
+  int frame_dim;
+  int spatial_dim;
+  int Voigt_dim;
+  int atom_dim;
+  int cell_spatial_dim;
+  int cell_angular_dim;
+  int label_dim;
+
+  int spatial_var;
+  int cell_spatial_var;
+  int cell_angular_var;
+
+  int time_var;
+  int cell_origin_var;
+  int cell_lengths_var;
+  int cell_angles_var;
+
+  virtual void openfile();
+  void closefile();
+  void write_time_and_cell();
+  virtual void write_data(int, double *);
+  void write_prmtop();
+
+  virtual int modify_param(int, char **);
+
+  void ncerr(int, const char *, int);
+
+  void compute_step(void *);
+  void compute_elapsed(void *);
+  void compute_elapsed_long(void *);
+};
+
+}
+
+#endif
+#endif
+#endif /* defined(LMP_HAS_PNETCDF) */
--- 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/compute_bond_local.cpp
+++ b/src/compute_bond_local.cpp
@ -148,7 +148,7 @@ int ComputeBondLocal::compute_bonds(int flag)
  double vvib,vrotsq;
  double inertia,omegasq;
  double mvv2e;
-  double engpot,engtrans,engvib,engrot,engtot,fbond;
+  double engpot,engtrans,engvib,engrot,fbond;
  double *ptr;

  double **x = atom->x;
--- 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/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
Steve Plimpton	3d254780de	Merge branch 'master' of github.com:lammps/lammps	2016-10-19 08:58:22 -06:00
sjplimp	b0cf1ded38	Merge pull request #229 from akohlmey/small-corrections collected various small corrections and updates	2016-10-19 08:57:54 -06:00
sjplimp	0891ed83a5	Merge pull request #228 from akohlmey/correct-prd-for-bigbig make PRD code compatible with -DLAMMPS_BIGBIG	2016-10-19 08:56:19 -06:00
Axel Kohlmeyer	ffb8eb36da	Revert "recognize any word starting with FP as FPRIME style table" This reverts commit `86f5b51133`.	2016-10-19 10:44:38 -04:00
Axel Kohlmeyer	1a10857244	make PRD code compatible with -DLAMMPS_BIGBIG	2016-10-19 10:43:25 -04:00
Axel Kohlmeyer	2634468676	dead code removal	2016-10-19 10:42:02 -04:00
Steve Plimpton	58141e0559	text edit	2016-10-19 08:41:44 -06:00
sjplimp	393337e7cf	Merge pull request #226 from akohlmey/pair-gauss-mixing add mixing support for pair styles gauss and gauss/cut	2016-10-19 08:37:32 -06:00
sjplimp	e0ed218cc4	Merge pull request #224 from akohlmey/get-vars provide an accessor for names of LAMMPS variables	2016-10-19 08:36:52 -06:00
sjplimp	efc7b23bd8	Merge pull request #216 from akohlmey/user-nc-dump USER-NC-DUMP package	2016-10-19 08:33:59 -06:00
Axel Kohlmeyer	566effc9c9	account for new fix halt properly in PDF and online docs	2016-10-19 09:08:20 -04:00
Axel Kohlmeyer	86f5b51133	recognize any word starting with FP as FPRIME style table	2016-10-19 09:05:25 -04:00
Axel Kohlmeyer	ccca642b3a	fix typo	2016-10-19 08:08:52 -04:00
Axel Kohlmeyer	019bc0ba0b	fix typos in Makefile.lammps for NetCDF support	2016-10-19 07:39:03 -04:00
Axel Kohlmeyer	7b7cb297e6	add mixing for pair styles gauss and gauss/cut contributed by andrew jewett. also add support for write_data.	2016-10-18 22:26:47 -04:00
Axel Kohlmeyer	5a6dea7a22	Merge branch 'master' into user-nc-dump	2016-10-18 15:36:18 -04:00
Axel Kohlmeyer	27796b431e	avoid distracting error output	2016-10-18 15:32:20 -04:00
Axel Kohlmeyer	6e54443d8c	make plain nc dump style a compile time option as well according to the package author, lars pastewka, some machines may have only parallel netcdf and not the serial netcdf installed. hence both netcdf dump styles neet to be optional depending on the compile time setting in lib/netcdf/Makefile.lammps The makefile has been set up to auto-detect the necessary flags for most regular and correct installations.	2016-10-18 15:11:22 -04:00
Axel Kohlmeyer	2b47fa6653	update Makefile.lammps and add README according to suggestions by lars pastewka	2016-10-18 15:05:30 -04:00
Axel Kohlmeyer	ae77664bdd	forgot to include USER-NC-DUMP in PACKUSER variable	2016-10-18 14:45:25 -04:00
Axel Kohlmeyer	4502d3276b	updated Makefile.lammps for NetCDF support based on suggestions from lars pastewka	2016-10-18 14:45:06 -04:00
Axel Kohlmeyer	c5513c4f75	provide an accessor for names of LAMMPS variables	2016-10-17 19:56:18 -04:00
Axel Kohlmeyer	2e8edbd2b9	USER-SMD was missing from PACKLIB variable	2016-10-12 06:53:37 -04:00
Axel Kohlmeyer	da66c1e649	update #include directives to match current LAMMPS coding conventions	2016-10-12 00:07:52 -04:00
Axel Kohlmeyer	e10ea91c5f	add doc file to toctree	2016-10-11 23:51:42 -04:00
Axel Kohlmeyer	69cb831705	import dump nc and dump nc/mpiio from lammps-netcdf as USER-NC-DUMP package by lars pastewka	2016-10-11 23:44:34 -04:00
Axel Kohlmeyer	ecb03dd2df	import lammps-netcdf code as USER-NC-DUMP package	2016-10-11 23:03:18 -04:00