diff --git a/doc/Developer.pdf b/doc/Developer.pdf
deleted file mode 100644
index 55b5015e4..000000000
Binary files a/doc/Developer.pdf and /dev/null differ
diff --git a/doc/Makefile b/doc/Makefile
index 52ed2d0d7..256dcd51c 100644
--- a/doc/Makefile
+++ b/doc/Makefile
@@ -1,87 +1,126 @@
 # Makefile for LAMMPS documentation
+
+SHELL 	      = /bin/bash
 SHA1          = $(shell echo $USER-$PWD | python utils/sha1sum.py)
 BUILDDIR      = /tmp/lammps-docs-$(SHA1)
 RSTDIR        = $(BUILDDIR)/rst
 VENV          = $(BUILDDIR)/docenv
 TXT2RST       = $(VENV)/bin/txt2rst
 
 PYTHON        = $(shell which python3)
 
 ifeq ($(shell which python3 >/dev/null 2>&1; echo $$?), 1)
 $(error Python3 was not found! Please check README.md for further instructions)
 endif
 
 ifeq ($(shell which virtualenv >/dev/null 2>&1; echo $$?), 1)
 $(error virtualenv was not found! Please check README.md for further instructions)
 endif
 
 SOURCES=$(wildcard src/*.txt)
 OBJECTS=$(SOURCES:src/%.txt=$(RSTDIR)/%.rst)
 
-.PHONY: help clean-all clean html pdf venv
+.PHONY: help clean-all clean html pdf old venv
+
+# ------------------------------------------
 
 help:
 	@echo "Please use \`make <target>' where <target> is one of"
-	@echo "  html       to make HTML version of documentation using Sphinx"
-	@echo "  pdf        to make Manual.pdf"
-	@echo "  txt2html   to build txt2html tool"
-	@echo "  clean      to remove all generated RST files"
-	@echo "  clean-all  to reset the entire build environment"
+	@echo "  html       create HTML doc pages in html dir"
+	@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 "  clean      remove all intermediate RST files"
+	@echo "  clean-all  reset the entire build environment"
+	@echo "  txt2html   build txt2html tool"
+
+# ------------------------------------------
 
 clean-all:
 	rm -rf $(BUILDDIR)/* utils/txt2html/txt2html.exe
 
 clean:
 	rm -rf $(RSTDIR)
 
-txt2html: utils/txt2html/txt2html.exe
-
 html: $(OBJECTS)
 	@(\
 		. $(VENV)/bin/activate ;\
 		cp -r src/* $(RSTDIR)/ ;\
 		sphinx-build -j 8 -b html -c utils/sphinx-config -d $(BUILDDIR)/doctrees $(RSTDIR) html ;\
 		deactivate ;\
 	)
 	-rm html/searchindex.js
-	-rm -rf html/_sources
+	@rm -rf html/_sources
+	@rm -rf html/PDF
+	@rm -rf html/USER
+	@cp -r src/PDF html/PDF
+	@cp -r src/USER html/USER
+	@rm -rf html/PDF/.[sg]*
+	@rm -rf html/USER/.[sg]*
+	@rm -rf html/USER/*/.[sg]*
+	@rm -rf html/USER/*/*.[sg]*
 	@echo "Build finished. The HTML pages are in doc/html."
 
 pdf: utils/txt2html/txt2html.exe
 	@(\
 		cd src; \
 		../utils/txt2html/txt2html.exe -b *.txt; \
-		htmldoc --batch ../lammps.book;          \
+		htmldoc --batch lammps.book;          \
 		for s in `echo *.txt | sed -e 's,\.txt,\.html,g'` ; \
-			do grep -q $$s ../lammps.book || \
-			echo doc file $$s missing in lammps.book; done; \
+			do grep -q $$s lammps.book || \
+			echo doc file $$s missing in src/lammps.book; done; \
 		rm *.html; \
+		cd Developer; \
+		pdflatex developer; \
+		pdflatex developer; \
+		mv developer.pdf ../../Developer.pdf; \
 	)
 
+old: utils/txt2html/txt2html.exe
+	@rm -rf old
+	@mkdir old; mkdir old/Eqs; mkdir old/JPG; mkdir old/PDF
+	@cd src; ../utils/txt2html/txt2html.exe -b *.txt; \
+	  mv *.html ../old; \
+	  cp Eqs/*.jpg ../old/Eqs; \
+	  cp JPG/* ../old/JPG; \
+	  cp PDF/* ../old/PDF;
+
+fetch:
+	@rm -rf html_www Manual_www.pdf Developer_www.pdf
+	@curl -s -o Manual_www.pdf http://lammps.sandia.gov/doc/Manual.pdf
+	@curl -s -o Developer_www.pdf http://lammps.sandia.gov/doc/Developer.pdf
+	@curl -s -o lammps-doc.tar.gz http://lammps.sandia.gov/tars/lammps-doc.tar.gz
+	@tar xzf lammps-doc.tar.gz
+	@rm -f lammps-doc.tar.gz
+
+txt2html: utils/txt2html/txt2html.exe
+
+# ------------------------------------------
+
 utils/txt2html/txt2html.exe: utils/txt2html/txt2html.cpp
 	g++ -O -Wall -o $@ $<
 
 $(RSTDIR)/%.rst : src/%.txt $(TXT2RST)
 	@(\
 		mkdir -p $(RSTDIR) ; \
 		. $(VENV)/bin/activate ;\
 		txt2rst $< > $@ ;\
 		deactivate ;\
 	)
 
 $(VENV):
 	@( \
 		virtualenv -p $(PYTHON) $(VENV); \
 		. $(VENV)/bin/activate; \
 		pip install Sphinx; \
 		pip install sphinxcontrib-images; \
 		deactivate;\
 	)
 
 $(TXT2RST): $(VENV)
 	@( \
 		. $(VENV)/bin/activate; \
 		(cd utils/converters;\
 		python setup.py develop);\
 		deactivate;\
 	)
diff --git a/doc/Manual.pdf b/doc/Manual.pdf
deleted file mode 100644
index 0f5f3497d..000000000
Binary files a/doc/Manual.pdf and /dev/null differ
diff --git a/doc/README b/doc/README
new file mode 100644
index 000000000..686a1802a
--- /dev/null
+++ b/doc/README
@@ -0,0 +1,53 @@
+Generation of LAMMPS Documentation
+
+The generation of all documentation is managed by the Makefile in this
+dir.
+
+----------------
+
+Options:
+
+make html         # generate HTML in html dir using Sphinx
+make pdf          # generate 2 PDF files (Manual.pdf,Developer.pdf)
+                  #   in this dir via htmldoc and pdflatex
+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 clean        # remove intermediate RST files created by HTML build
+make clean-all    # remove entire build folder and any cached data
+
+----------------
+
+Installing prerequisites for HTML build
+
+To run the HTML documention build toolchain, Python 3 and virtualenv
+have to be installed.  Here are instructions for common setups:
+
+# Ubuntu
+
+sudo apt-get install python-virtualenv
+
+# Fedora (up to version 21)
+# Red Hat Enterprise Linux or CentOS (up to version 7.x)
+
+sudo yum install python3-virtualenv
+
+# Fedora (since version 22)
+
+sudo dnf install python3-virtualenv
+
+# MacOS X
+
+## Python 3
+
+Download the latest Python 3 MacOS X package from
+https://www.python.org and install it.  This will install both Python
+3 and pip3.
+
+## virtualenv
+
+Once Python 3 is installed, open a Terminal and type
+
+pip3 install virtualenv
+
+This will install virtualenv from the Python Package Index.
diff --git a/doc/README.md b/doc/README.md
deleted file mode 100644
index b843b7eab..000000000
--- a/doc/README.md
+++ /dev/null
@@ -1,48 +0,0 @@
-# Generation of LAMMPS Documentation
-
-The generation of all the documentation is managed by the Makefile inside the
-`doc/` folder.
-
-## Usage:
-
-```bash
-make html         # generate HTML using Sphinx
-make pdf          # generate PDF using htmldoc
-make clean        # remove generated RST files
-make clean-all    # remove entire build folder and any cached data
-```
-
-## Installing prerequisites
-
-To run the documention build toolchain, Python 3 and virtualenv have
-to be installed. Here are instructions for common setups:
-
-### Ubuntu
-
-```bash
-sudo apt-get install python-virtualenv
-```
-
-### Fedora (up to version 21), Red Hat Enterprise Linux or CentOS (up to version 7.x)
-
-```bash
-sudo yum install python3-virtualenv
-```
-
-### Fedora (since version 22)
-
-```bash
-sudo dnf install python3-virtualenv
-```
-
-### MacOS X
-
-## Python 3
-
-Download the latest Python 3 MacOS X package from https://www.python.org and install it.
-This will install both Python 3 and pip3.
-
-## virtualenv
-
-Once Python 3 is installed, open a Terminal and type `pip3 install virtualenv`. This will
-install virtualenv from the Python Package Index.
diff --git a/doc/src/Manual.txt b/doc/src/Manual.txt
index 3e822bd96..b276794a5 100644
--- a/doc/src/Manual.txt
+++ b/doc/src/Manual.txt
@@ -1,340 +1,340 @@
 <!-- HTML_ONLY -->
 <HEAD>
 <TITLE>LAMMPS Users Manual</TITLE>
-<META NAME="docnumber" CONTENT="20 Sep 2016 version">
+<META NAME="docnumber" CONTENT="22 Sep 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>
 
 <BODY>
 
 <!-- END_HTML_ONLY -->
 
 "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
 
 <H1></H1>
 
 LAMMPS Documentation :c,h3
-20 Sep 2016 version :c,h4
+22 Sep 2016 version :c,h4
 
 Version info: :h4
 
 The LAMMPS "version" is the date when it was released, such as 1 May
 2010. LAMMPS is updated continuously.  Whenever we fix a bug or add a
 feature, we release it immediately, and post a notice on "this page of
 the WWW site"_bug. Every 2-4 months one of the incremental releases
 is subjected to more thorough testing and labeled as a {stable} version.
 
 Each dated copy of LAMMPS contains all the
 features and bug-fixes up to and including that version date. The
 version date is printed to the screen and logfile every time you run
 LAMMPS. It is also in the file src/version.h and in the LAMMPS
 directory name created when you unpack a tarball, and at the top of
 the first page of the manual (this page).
 
 If you browse the HTML doc pages on the LAMMPS WWW site, they always
 describe the most current version of LAMMPS. :ulb,l
 
 If you browse the HTML doc pages included in your tarball, they
 describe the version you have. :l
 
 The "PDF file"_Manual.pdf on the WWW site or in the tarball is updated
 about once per month.  This is because it is large, and we don't want
 it to be part of every patch. :l
 
 There is also a "Developer.pdf"_Developer.pdf file in the doc
 directory, which describes the internal structure and algorithms of
 LAMMPS.  :l
 :ule
 
 LAMMPS stands for Large-scale Atomic/Molecular Massively Parallel
 Simulator.
 
 LAMMPS 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 current core group of LAMMPS developers is at Sandia National
 Labs and Temple University:
 
 "Steve Plimpton"_sjp, sjplimp at sandia.gov :ulb,l
 Aidan Thompson, athomps at sandia.gov :l
 Stan Moore, stamoore at sandia.gov :l
 "Axel Kohlmeyer"_ako, akohlmey at gmail.com :l
 :ule
 
 Past core developers include Paul Crozier, Ray Shan and Mark Stevens,
 all at Sandia. The [LAMMPS home page] at
 "http://lammps.sandia.gov"_http://lammps.sandia.gov has more information
 about the code and its uses. Interaction with external LAMMPS developers,
 bug reports and feature requests are mainly coordinated through the
 "LAMMPS project on GitHub."_https://github.com/lammps/lammps
 The lammps.org domain, currently hosting "public continuous integration
 testing"_https://ci.lammps.org/job/lammps/ and "precompiled Linux
 RPM and Windows installer packages"_http://rpm.lammps.org is located
 at Temple University and managed by Richard Berger,
 richard.berger at temple.edu.
 
 :link(bug,http://lammps.sandia.gov/bug.html)
 :link(sjp,http://www.sandia.gov/~sjplimp)
 :link(ako,http://goo.gl/1wk0)
 
 :line
 
 The LAMMPS documentation is organized into the following sections.  If
 you find errors or omissions in this manual or have suggestions for
 useful information to add, please send an email to the developers so
 we can improve the LAMMPS documentation.
 
 Once you are familiar with LAMMPS, you may want to bookmark "this
 page"_Section_commands.html#comm at Section_commands.html#comm since
 it gives quick access to documentation for all LAMMPS commands.
 
 "PDF file"_Manual.pdf of the entire manual, generated by
 "htmldoc"_http://freecode.com/projects/htmldoc
 
 <!-- RST
 
 .. toctree::
    :maxdepth: 2
    :numbered:
    :caption: User Documentation
    :name: userdoc
    :includehidden:
    
    Section_intro
    Section_start
    Section_commands
    Section_packages
    Section_accelerate
    Section_howto
    Section_example
    Section_perf
    Section_tools
    Section_modify
    Section_python
    Section_errors
    Section_history
 
 .. toctree::
    :caption: Index
    :name: index
    :hidden:
 
    tutorials
    commands
    fixes
    computes
    pairs
    bonds
    angles
    dihedrals
    impropers
 
 Indices and tables
 ==================
 
 * :ref:`genindex`
 * :ref:`search`
    
 END_RST -->
 
 <!-- HTML_ONLY -->
 "Introduction"_Section_intro.html :olb,l
   1.1 "What is LAMMPS"_intro_1 :ulb,b
   1.2 "LAMMPS features"_intro_2 :b
   1.3 "LAMMPS non-features"_intro_3 :b
   1.4 "Open source distribution"_intro_4 :b
   1.5 "Acknowledgments and citations"_intro_5 :ule,b
 "Getting started"_Section_start.html :l
   2.1 "What's in the LAMMPS distribution"_start_1 :ulb,b
   2.2 "Making LAMMPS"_start_2 :b
   2.3 "Making LAMMPS with optional packages"_start_3 :b
   2.4 "Building LAMMPS via the Make.py script"_start_4 :b
   2.5 "Building LAMMPS as a library"_start_5 :b
   2.6 "Running LAMMPS"_start_6 :b
   2.7 "Command-line options"_start_7 :b
   2.8 "Screen output"_start_8 :b
   2.9 "Tips for users of previous versions"_start_9 :ule,b
 "Commands"_Section_commands.html :l
   3.1 "LAMMPS input script"_cmd_1 :ulb,b
   3.2 "Parsing rules"_cmd_2 :b
   3.3 "Input script structure"_cmd_3 :b
   3.4 "Commands listed by category"_cmd_4 :b
   3.5 "Commands listed alphabetically"_cmd_5 :ule,b
 "Packages"_Section_packages.html :l
   4.1 "Standard packages"_pkg_1 :ulb,b
   4.2 "User packages"_pkg_2 :ule,b
 "Accelerating LAMMPS performance"_Section_accelerate.html :l
   5.1 "Measuring performance"_acc_1 :ulb,b
   5.2 "Algorithms and code options to boost performace"_acc_2 :b
   5.3 "Accelerator packages with optimized styles"_acc_3 :b
     5.3.1 "GPU package"_accelerate_gpu.html :ulb,b
     5.3.2 "USER-INTEL package"_accelerate_intel.html :b
     5.3.3 "KOKKOS package"_accelerate_kokkos.html :b
     5.3.4 "USER-OMP package"_accelerate_omp.html :b
     5.3.5 "OPT package"_accelerate_opt.html :ule,b
   5.4 "Comparison of various accelerator packages"_acc_4 :ule,b
 "How-to discussions"_Section_howto.html :l
   6.1 "Restarting a simulation"_howto_1 :ulb,b
   6.2 "2d simulations"_howto_2 :b
   6.3 "CHARMM and AMBER force fields"_howto_3 :b
   6.4 "Running multiple simulations from one input script"_howto_4 :b
   6.5 "Multi-replica simulations"_howto_5 :b
   6.6 "Granular models"_howto_6 :b
   6.7 "TIP3P water model"_howto_7 :b
   6.8 "TIP4P water model"_howto_8 :b
   6.9 "SPC water model"_howto_9 :b
   6.10 "Coupling LAMMPS to other codes"_howto_10 :b
   6.11 "Visualizing LAMMPS snapshots"_howto_11 :b
   6.12 "Triclinic (non-orthogonal) simulation boxes"_howto_12 :b
   6.13 "NEMD simulations"_howto_13 :b
   6.14 "Finite-size spherical and aspherical particles"_howto_14 :b
   6.15 "Output from LAMMPS (thermo, dumps, computes, fixes, variables)"_howto_15 :b
   6.16 "Thermostatting, barostatting, and compute temperature"_howto_16 :b
   6.17 "Walls"_howto_17 :b
   6.18 "Elastic constants"_howto_18 :b
   6.19 "Library interface to LAMMPS"_howto_19 :b
   6.20 "Calculating thermal conductivity"_howto_20 :b
   6.21 "Calculating viscosity"_howto_21 :b
   6.22 "Calculating a diffusion coefficient"_howto_22 :b
   6.23 "Using chunks to calculate system properties"_howto_23 :b
   6.24 "Setting parameters for pppm/disp"_howto_24 :b
   6.25 "Polarizable models"_howto_25 :b
   6.26 "Adiabatic core/shell model"_howto_26 :b
   6.27 "Drude induced dipoles"_howto_27 :ule,b
 "Example problems"_Section_example.html :l
 "Performance & scalability"_Section_perf.html :l
 "Additional tools"_Section_tools.html :l
 "Modifying & extending LAMMPS"_Section_modify.html :l
   10.1 "Atom styles"_mod_1 :ulb,b
   10.2 "Bond, angle, dihedral, improper potentials"_mod_2 :b
   10.3 "Compute styles"_mod_3 :b
   10.4 "Dump styles"_mod_4 :b
   10.5 "Dump custom output options"_mod_5 :b
   10.6 "Fix styles"_mod_6 :b
   10.7 "Input script commands"_mod_7 :b
   10.8 "Kspace computations"_mod_8 :b
   10.9 "Minimization styles"_mod_9 :b
   10.10 "Pairwise potentials"_mod_10 :b
   10.11 "Region styles"_mod_11 :b
   10.12 "Body styles"_mod_12 :b
   10.13 "Thermodynamic output options"_mod_13 :b
   10.14 "Variable options"_mod_14 :b
   10.15 "Submitting new features for inclusion in LAMMPS"_mod_15 :ule,b
 "Python interface"_Section_python.html :l
   11.1 "Overview of running LAMMPS from Python"_py_1 :ulb,b
   11.2 "Overview of using Python from a LAMMPS script"_py_2 :b
   11.3 "Building LAMMPS as a shared library"_py_3 :b
   11.4 "Installing the Python wrapper into Python"_py_4 :b
   11.5 "Extending Python with MPI to run in parallel"_py_5 :b
   11.6 "Testing the Python-LAMMPS interface"_py_6 :b
   11.7 "Using LAMMPS from Python"_py_7 :b
   11.8 "Example Python scripts that use LAMMPS"_py_8 :ule,b
 "Errors"_Section_errors.html :l
   12.1 "Common problems"_err_1 :ulb,b
   12.2 "Reporting bugs"_err_2 :b
   12.3 "Error & warning messages"_err_3 :ule,b
 "Future and history"_Section_history.html :l
   13.1 "Coming attractions"_hist_1 :ulb,b
   13.2 "Past versions"_hist_2 :ule,b
 :ole
 
 :link(intro_1,Section_intro.html#intro_1)
 :link(intro_2,Section_intro.html#intro_2)
 :link(intro_3,Section_intro.html#intro_3)
 :link(intro_4,Section_intro.html#intro_4)
 :link(intro_5,Section_intro.html#intro_5)
 
 :link(start_1,Section_start.html#start_1)
 :link(start_2,Section_start.html#start_2)
 :link(start_3,Section_start.html#start_3)
 :link(start_4,Section_start.html#start_4)
 :link(start_5,Section_start.html#start_5)
 :link(start_6,Section_start.html#start_6)
 :link(start_7,Section_start.html#start_7)
 :link(start_8,Section_start.html#start_8)
 :link(start_9,Section_start.html#start_9)
 
 :link(cmd_1,Section_commands.html#cmd_1)
 :link(cmd_2,Section_commands.html#cmd_2)
 :link(cmd_3,Section_commands.html#cmd_3)
 :link(cmd_4,Section_commands.html#cmd_4)
 :link(cmd_5,Section_commands.html#cmd_5)
 
 :link(pkg_1,Section_packages.html#pkg_1)
 :link(pkg_2,Section_packages.html#pkg_2)
 
 :link(acc_1,Section_accelerate.html#acc_1)
 :link(acc_2,Section_accelerate.html#acc_2)
 :link(acc_3,Section_accelerate.html#acc_3)
 :link(acc_4,Section_accelerate.html#acc_4)
 
 :link(howto_1,Section_howto.html#howto_1)
 :link(howto_2,Section_howto.html#howto_2)
 :link(howto_3,Section_howto.html#howto_3)
 :link(howto_4,Section_howto.html#howto_4)
 :link(howto_5,Section_howto.html#howto_5)
 :link(howto_6,Section_howto.html#howto_6)
 :link(howto_7,Section_howto.html#howto_7)
 :link(howto_8,Section_howto.html#howto_8)
 :link(howto_9,Section_howto.html#howto_9)
 :link(howto_10,Section_howto.html#howto_10)
 :link(howto_11,Section_howto.html#howto_11)
 :link(howto_12,Section_howto.html#howto_12)
 :link(howto_13,Section_howto.html#howto_13)
 :link(howto_14,Section_howto.html#howto_14)
 :link(howto_15,Section_howto.html#howto_15)
 :link(howto_16,Section_howto.html#howto_16)
 :link(howto_17,Section_howto.html#howto_17)
 :link(howto_18,Section_howto.html#howto_18)
 :link(howto_19,Section_howto.html#howto_19)
 :link(howto_20,Section_howto.html#howto_20)
 :link(howto_21,Section_howto.html#howto_21)
 :link(howto_22,Section_howto.html#howto_22)
 :link(howto_23,Section_howto.html#howto_23)
 :link(howto_24,Section_howto.html#howto_24)
 :link(howto_25,Section_howto.html#howto_25)
 :link(howto_26,Section_howto.html#howto_26)
 :link(howto_27,Section_howto.html#howto_27)
 
 :link(mod_1,Section_modify.html#mod_1)
 :link(mod_2,Section_modify.html#mod_2)
 :link(mod_3,Section_modify.html#mod_3)
 :link(mod_4,Section_modify.html#mod_4)
 :link(mod_5,Section_modify.html#mod_5)
 :link(mod_6,Section_modify.html#mod_6)
 :link(mod_7,Section_modify.html#mod_7)
 :link(mod_8,Section_modify.html#mod_8)
 :link(mod_9,Section_modify.html#mod_9)
 :link(mod_10,Section_modify.html#mod_10)
 :link(mod_11,Section_modify.html#mod_11)
 :link(mod_12,Section_modify.html#mod_12)
 :link(mod_13,Section_modify.html#mod_13)
 :link(mod_14,Section_modify.html#mod_14)
 :link(mod_15,Section_modify.html#mod_15)
 
 :link(py_1,Section_python.html#py_1)
 :link(py_2,Section_python.html#py_2)
 :link(py_3,Section_python.html#py_3)
 :link(py_4,Section_python.html#py_4)
 :link(py_5,Section_python.html#py_5)
 :link(py_6,Section_python.html#py_6)
 
 :link(err_1,Section_errors.html#err_1)
 :link(err_2,Section_errors.html#err_2)
 :link(err_3,Section_errors.html#err_3)
 
 :link(hist_1,Section_history.html#hist_1)
 :link(hist_2,Section_history.html#hist_2)
 <!-- END_HTML_ONLY -->
 
 </BODY>
diff --git a/doc/lammps.book b/doc/src/lammps.book
similarity index 99%
rename from doc/lammps.book
rename to doc/src/lammps.book
index f28bb48b5..620eee483 100644
--- a/doc/lammps.book
+++ b/doc/src/lammps.book
@@ -1,624 +1,622 @@
 #HTMLDOC 1.8.27
 -t pdf14 -f "../Manual.pdf" --book --toclevels 4 --no-numbered --toctitle "Table of Contents" --title --textcolor #000000 --linkcolor #0000ff --linkstyle plain --bodycolor #ffffff --size Universal --left 1.00in --right 0.50in --top 0.50in --bottom 0.50in --header .t. --header1 ... --footer ..1 --nup 1 --tocheader .t. --tocfooter ..i --portrait --color --no-pscommands --no-xrxcomments --compression=1 --jpeg=0 --fontsize 11.0 --fontspacing 1.2 --headingfont helvetica --bodyfont times --headfootsize 11.0 --headfootfont helvetica --charset iso-8859-1 --links --embedfonts --pagemode document --pagelayout single --firstpage c1 --pageeffect none --pageduration 10 --effectduration 1.0 --no-encryption --permissions all  --owner-password ""  --user-password "" --browserwidth 680 --no-strict --no-overflow
 
 Manual.html
 Section_intro.html
 Section_start.html
 Section_commands.html
 Section_packages.html
 Section_accelerate.html
 accelerate_gpu.html
 accelerate_intel.html
 accelerate_kokkos.html
 accelerate_omp.html
 accelerate_opt.html
 Section_howto.html
 Section_example.html
 Section_perf.html
 Section_tools.html
 Section_modify.html
 Section_python.html
 Section_errors.html
 Section_history.html
 
 tutorial_drude.html
 tutorial_github.html
 
 body.html
 manifolds.html
 
 angle_coeff.html
 angle_style.html
 atom_modify.html
 atom_style.html
 balance.html
 bond_coeff.html
 bond_style.html
 bond_write.html
 boundary.html
 box.html
 change_box.html
 clear.html
 comm_modify.html
 comm_style.html
 compute.html
 compute_modify.html
 create_atoms.html
 create_bonds.html
 create_box.html
 delete_atoms.html
 delete_bonds.html
 dielectric.html
 dihedral_coeff.html
 dihedral_style.html
 dimension.html
 displace_atoms.html
 dump.html
 dump_custom_vtk.html
 dump_h5md.html
 dump_image.html
 dump_modify.html
 dump_molfile.html
 echo.html
 fix.html
 fix_modify.html
 group.html
 group2ndx.html
 if.html
 improper_coeff.html
 improper_style.html
 include.html
 info.html
 jump.html
 kspace_modify.html
 kspace_style.html
 label.html
 lattice.html
 log.html
 mass.html
 min_modify.html
 min_style.html
 minimize.html
 molecule.html
 neb.html
 neigh_modify.html
 neighbor.html
 newton.html
 next.html
 package.html
 pair_coeff.html
 pair_modify.html
 pair_style.html
 pair_write.html
 partition.html
 prd.html
 print.html
 processors.html
 python.html
 quit.html
 read_data.html
 read_dump.html
 read_restart.html
 region.html
 replicate.html
 rerun.html
 reset_timestep.html
 restart.html
 run.html
 run_style.html
 set.html
 shell.html
 special_bonds.html
 suffix.html
 tad.html
 temper.html
 thermo.html
 thermo_modify.html
 thermo_style.html
 timer.html
 timestep.html
 uncompute.html
 undump.html
 unfix.html
 units.html
 variable.html
 velocity.html
 write_coeff.html
 write_data.html
 write_dump.html
 write_restart.html
 
 fix_adapt.html
 fix_adapt_fep.html
 fix_addforce.html
 fix_addtorque.html
 fix_append_atoms.html
 fix_atc.html
 fix_atom_swap.html
 fix_ave_atom.html
 fix_ave_chunk.html
 fix_ave_correlate.html
 fix_ave_correlate_long.html
 fix_ave_histo.html
 fix_ave_time.html
 fix_aveforce.html
 fix_balance.html
 fix_bond_break.html
 fix_bond_create.html
 fix_bond_swap.html
 fix_box_relax.html
 fix_colvars.html
 fix_controller.html
 fix_deform.html
 fix_deposit.html
 fix_drag.html
 fix_drude.html
 fix_drude_transform.html
 fix_dt_reset.html
 fix_efield.html
 fix_ehex.html
 fix_enforce2d.html
 fix_eos_cv.html
 fix_eos_table.html
 fix_eos_table_rx.html
 fix_evaporate.html
 fix_external.html
 fix_flow_gauss.html
 fix_freeze.html
 fix_gcmc.html
 fix_gld.html
 fix_gle.html
 fix_gravity.html
 fix_heat.html
 fix_imd.html
 fix_indent.html
 fix_ipi.html
 fix_langevin.html
 fix_langevin_drude.html
 fix_langevin_eff.html
 fix_lb_fluid.html
 fix_lb_momentum.html
 fix_lb_pc.html
 fix_lb_rigid_pc_sphere.html
 fix_lb_viscous.html
 fix_lineforce.html
 fix_manifoldforce.html
 fix_meso.html
 fix_meso_stationary.html
 fix_momentum.html
 fix_move.html
 fix_msst.html
 fix_neb.html
 fix_nh.html
 fix_nh_eff.html
 fix_nph_asphere.html
 fix_nph_body.html
 fix_nph_sphere.html
 fix_nphug.html
 fix_npt_asphere.html
 fix_npt_body.html
 fix_npt_sphere.html
 fix_nve.html
 fix_nve_asphere.html
 fix_nve_asphere_noforce.html
 fix_nve_body.html
 fix_nve_eff.html
 fix_nve_limit.html
 fix_nve_line.html
 fix_nve_manifold_rattle.html
 fix_nve_noforce.html
 fix_nve_sphere.html
 fix_nve_tri.html
 fix_nvt_asphere.html
 fix_nvt_body.html
 fix_nvt_manifold_rattle.html
 fix_nvt_sllod.html
 fix_nvt_sllod_eff.html
 fix_nvt_sphere.html
 fix_oneway.html
 fix_orient.html
 fix_phonon.html
 fix_pimd.html
 fix_planeforce.html
 fix_poems.html
 fix_pour.html
 fix_press_berendsen.html
 fix_print.html
 fix_property_atom.html
 fix_qbmsst.html
 fix_qeq.html
 fix_qeq_comb.html
 fix_qeq_reax.html
 fix_qmmm.html
 fix_qtb.html
 fix_reax_bonds.html
 fix_reaxc_species.html
 fix_recenter.html
 fix_restrain.html
 fix_rigid.html
 fix_rx.html
 fix_saed_vtk.html
 fix_setforce.html
 fix_shake.html
 fix_shardlow.html
 fix_smd.html
 fix_smd_adjust_dt.html
 fix_smd_integrate_tlsph.html
 fix_smd_integrate_ulsph.html
 fix_smd_move_triangulated_surface.html
 fix_smd_setvel.html
-fix_smd_tlsph_reference_configuration.html
 fix_smd_wall_surface.html
 fix_spring.html
 fix_spring_chunk.html
 fix_spring_rg.html
 fix_spring_self.html
 fix_srd.html
 fix_store_force.html
 fix_store_state.html
 fix_temp_berendsen.html
 fix_temp_csvr.html
 fix_temp_rescale.html
 fix_temp_rescale_eff.html
 fix_tfmc.html
 fix_thermal_conductivity.html
-fix_ti_rs.html
 fix_ti_spring.html
 fix_tmd.html
 fix_ttm.html
 fix_tune_kspace.html
 fix_vector.html
 fix_viscosity.html
 fix_viscous.html
 fix_wall.html
 fix_wall_gran.html
 fix_wall_piston.html
 fix_wall_reflect.html
 fix_wall_region.html
 fix_wall_srd.html
 
 compute_ackland_atom.html
 compute_angle.html
 compute_angle_local.html
 compute_angmom_chunk.html
 compute_basal_atom.html
 compute_body_local.html
 compute_bond.html
 compute_bond_local.html
 compute_centro_atom.html
 compute_chunk_atom.html
 compute_cluster_atom.html
 compute_cna_atom.html
 compute_com.html
 compute_com_chunk.html
 compute_contact_atom.html
 compute_coord_atom.html
 compute_damage_atom.html
 compute_dihedral.html
 compute_dihedral_local.html
 compute_dilatation_atom.html
 compute_dipole_chunk.html
 compute_displace_atom.html
 compute_dpd.html
 compute_dpd_atom.html
 compute_erotate_asphere.html
 compute_erotate_rigid.html
 compute_erotate_sphere.html
 compute_erotate_sphere_atom.html
 compute_event_displace.html
 compute_fep.html
 compute_group_group.html
 compute_gyration.html
 compute_gyration_chunk.html
 compute_heat_flux.html
 compute_hexorder_atom.html
 compute_improper.html
 compute_improper_local.html
 compute_inertia_chunk.html
 compute_ke.html
 compute_ke_atom.html
 compute_ke_atom_eff.html
 compute_ke_eff.html
 compute_ke_rigid.html
 compute_meso_e_atom.html
 compute_meso_rho_atom.html
 compute_meso_t_atom.html
 compute_msd.html
 compute_msd_chunk.html
 compute_msd_nongauss.html
 compute_omega_chunk.html
 compute_orientorder_atom.html
 compute_pair.html
 compute_pair_local.html
 compute_pe.html
 compute_pe_atom.html
 compute_plasticity_atom.html
 compute_pressure.html
 compute_property_atom.html
 compute_property_chunk.html
 compute_property_local.html
 compute_rdf.html
 compute_reduce.html
 compute_rigid_local.html
 compute_saed.html
 compute_slice.html
 compute_smd_contact_radius.html
 compute_smd_damage.html
 compute_smd_hourglass_error.html
 compute_smd_internal_energy.html
 compute_smd_plastic_strain.html
 compute_smd_plastic_strain_rate.html
 compute_smd_rho.html
 compute_smd_tlsph_defgrad.html
 compute_smd_tlsph_dt.html
 compute_smd_tlsph_num_neighs.html
 compute_smd_tlsph_shape.html
 compute_smd_tlsph_strain.html
 compute_smd_tlsph_strain_rate.html
 compute_smd_tlsph_stress.html
 compute_smd_triangle_mesh_vertices.html
 compute_smd_ulsph_num_neighs.html
 compute_smd_ulsph_strain.html
 compute_smd_ulsph_strain_rate.html
 compute_smd_ulsph_stress.html
 compute_smd_vol.html
 compute_sna_atom.html
 compute_stress_atom.html
 compute_tally.html
 compute_temp.html
 compute_temp_asphere.html
 compute_temp_body.html
 compute_temp_chunk.html
 compute_temp_com.html
 compute_temp_cs.html
 compute_temp_deform.html
 compute_temp_deform_eff.html
 compute_temp_drude.html
 compute_temp_eff.html
 compute_temp_partial.html
 compute_temp_profile.html
 compute_temp_ramp.html
 compute_temp_region.html
 compute_temp_region_eff.html
 compute_temp_rotate.html
 compute_temp_sphere.html
 compute_ti.html
 compute_torque_chunk.html
 compute_vacf.html
 compute_vcm_chunk.html
 compute_voronoi_atom.html
 compute_xrd.html
 
 pair_adp.html
 pair_airebo.html
 pair_awpmd.html
 pair_beck.html
 pair_body.html
 pair_bop.html
 pair_born.html
 pair_brownian.html
 pair_buck.html
 pair_buck_long.html
 pair_charmm.html
 pair_class2.html
 pair_colloid.html
 pair_comb.html
 pair_coul.html
 pair_coul_diel.html
 pair_cs.html
 pair_dipole.html
 pair_dpd.html
 pair_dpd_fdt.html
 pair_dsmc.html
 pair_eam.html
 pair_edip.html
 pair_eff.html
 pair_eim.html
 pair_exp6_rx.html
 pair_gauss.html
 pair_gayberne.html
 pair_gran.html
 pair_gromacs.html
 pair_hbond_dreiding.html
 pair_hybrid.html
 pair_kim.html
 pair_lcbop.html
 pair_line_lj.html
 pair_list.html
 pair_lj.html
 pair_lj96.html
 pair_lj_cubic.html
 pair_lj_expand.html
 pair_lj_long.html
 pair_lj_sf.html
 pair_lj_smooth.html
 pair_lj_smooth_linear.html
 pair_lj_soft.html
 pair_lubricate.html
 pair_lubricateU.html
 pair_mdf.html
 pair_meam.html
 pair_meam_spline.html
 pair_meam_sw_spline.html
 pair_mgpt.html
 pair_mie.html
 pair_morse.html
 pair_multi_lucy.html
 pair_multi_lucy_rx.html
 pair_nb3b_harmonic.html
 pair_nm.html
 pair_none.html
 pair_peri.html
 pair_polymorphic.html
 pair_quip.html
 pair_reax.html
 pair_reax_c.html
 pair_resquared.html
 pair_sdk.html
 pair_smd_hertz.html
 pair_smd_tlsph.html
 pair_smd_triangulated_surface.html
 pair_smd_ulsph.html
 pair_smtbq.html
 pair_snap.html
 pair_soft.html
 pair_sph_heatconduction.html
 pair_sph_idealgas.html
 pair_sph_lj.html
 pair_sph_rhosum.html
 pair_sph_taitwater.html
 pair_sph_taitwater_morris.html
 pair_srp.html
 pair_sw.html
 pair_table.html
 pair_table_rx.html
 pair_tersoff.html
 pair_tersoff_mod.html
 pair_tersoff_zbl.html
 pair_thole.html
 pair_tri_lj.html
 pair_vashishta.html
 pair_yukawa.html
 pair_yukawa_colloid.html
 pair_zbl.html
 pair_zero.html
 
 bond_class2.html
 bond_fene.html
 bond_fene_expand.html
 bond_harmonic.html
 bond_harmonic_shift.html
 bond_harmonic_shift_cut.html
 bond_hybrid.html
 bond_morse.html
 bond_none.html
 bond_nonlinear.html
 bond_quartic.html
 bond_table.html
 bond_zero.html
 
 angle_charmm.html
 angle_class2.html
 angle_cosine.html
 angle_cosine_delta.html
 angle_cosine_periodic.html
 angle_cosine_shift.html
 angle_cosine_shift_exp.html
 angle_cosine_squared.html
 angle_dipole.html
 angle_fourier.html
 angle_fourier_simple.html
 angle_harmonic.html
 angle_hybrid.html
 angle_none.html
 angle_quartic.html
 angle_sdk.html
 angle_table.html
 angle_zero.html
 
 dihedral_charmm.html
 dihedral_class2.html
 dihedral_cosine_shift_exp.html
 dihedral_fourier.html
 dihedral_harmonic.html
 dihedral_helix.html
 dihedral_hybrid.html
 dihedral_multi_harmonic.html
 dihedral_nharmonic.html
 dihedral_none.html
 dihedral_opls.html
 dihedral_quadratic.html
 dihedral_spherical.html
 dihedral_table.html
 dihedral_zero.html
 
 improper_class2.html
 improper_cossq.html
 improper_cvff.html
 improper_distance.html
 improper_fourier.html
 improper_harmonic.html
 improper_hybrid.html
 improper_none.html
 improper_ring.html
 improper_umbrella.html
 improper_zero.html
 
 USER/atc/man_add_molecule.html
 USER/atc/man_add_species.html
 USER/atc/man_atom_element_map.html
 USER/atc/man_atom_weight.html
 USER/atc/man_atomic_charge.html
 USER/atc/man_boundary.html
 USER/atc/man_boundary_dynamics.html
 USER/atc/man_boundary_faceset.html
 USER/atc/man_boundary_integral.html
 USER/atc/man_consistent_fe_initialization.html
 USER/atc/man_contour_integral.html
 USER/atc/man_control.html
 USER/atc/man_control_momentum.html
 USER/atc/man_control_thermal.html
 USER/atc/man_control_thermal_correction_max_iterations.html
 USER/atc/man_decomposition.html
 USER/atc/man_electron_integration.html
 USER/atc/man_equilibrium_start.html
 USER/atc/man_extrinsic_exchange.html
 USER/atc/man_fe_md_boundary.html
 USER/atc/man_fem_mesh.html
 USER/atc/man_filter_scale.html
 USER/atc/man_filter_type.html
 USER/atc/man_fix_atc.html
 USER/atc/man_fix_flux.html
 USER/atc/man_fix_nodes.html
 USER/atc/man_hardy_computes.html
 USER/atc/man_hardy_fields.html
 USER/atc/man_hardy_gradients.html
 USER/atc/man_hardy_kernel.html
 USER/atc/man_hardy_on_the_fly.html
 USER/atc/man_hardy_rates.html
 USER/atc/man_initial.html
 USER/atc/man_internal_atom_integrate.html
 USER/atc/man_internal_element_set.html
 USER/atc/man_internal_quadrature.html
 USER/atc/man_kernel_function.html
 USER/atc/man_localized_lambda.html
 USER/atc/man_lumped_lambda_solve.html
 USER/atc/man_mask_direction.html
 USER/atc/man_mass_matrix.html
 USER/atc/man_material.html
 USER/atc/man_mesh_add_to_nodeset.html
 USER/atc/man_mesh_create.html
 USER/atc/man_mesh_create_elementset.html
 USER/atc/man_mesh_create_faceset_box.html
 USER/atc/man_mesh_create_faceset_plane.html
 USER/atc/man_mesh_create_nodeset.html
 USER/atc/man_mesh_delete_elements.html
 USER/atc/man_mesh_nodeset_to_elementset.html
 USER/atc/man_mesh_output.html
 USER/atc/man_mesh_quadrature.html
 USER/atc/man_mesh_read.html
 USER/atc/man_mesh_write.html
 USER/atc/man_momentum_time_integration.html
 USER/atc/man_output.html
 USER/atc/man_output_elementset.html
 USER/atc/man_output_nodeset.html
 USER/atc/man_pair_interactions.html
 USER/atc/man_poisson_solver.html
 USER/atc/man_read_restart.html
 USER/atc/man_remove_molecule.html
 USER/atc/man_remove_source.html
 USER/atc/man_remove_species.html
 USER/atc/man_reset_atomic_reference_positions.html
 USER/atc/man_reset_time.html
 USER/atc/man_sample_frequency.html
 USER/atc/man_set.html
 USER/atc/man_source.html
 USER/atc/man_source_integration.html
 USER/atc/man_temperature_definition.html
 USER/atc/man_thermal_time_integration.html
 USER/atc/man_time_filter.html
 USER/atc/man_track_displacement.html
 USER/atc/man_unfix_flux.html
 USER/atc/man_unfix_nodes.html
 USER/atc/man_write_atom_weights.html
 USER/atc/man_write_restart.html
diff --git a/src/fix_move.cpp b/src/fix_move.cpp
index 3c0797738..d3c4b12e3 100644
--- a/src/fix_move.cpp
+++ b/src/fix_move.cpp
@@ -1,1249 +1,1251 @@
 /* ----------------------------------------------------------------------
    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.
 ------------------------------------------------------------------------- */
 
 #include <string.h>
 #include <stdlib.h>
 #include <math.h>
 #include "fix_move.h"
 #include "atom.h"
 #include "group.h"
 #include "update.h"
 #include "modify.h"
 #include "force.h"
 #include "domain.h"
 #include "lattice.h"
 #include "comm.h"
 #include "respa.h"
 #include "input.h"
 #include "variable.h"
 #include "atom_vec_ellipsoid.h"
 #include "atom_vec_line.h"
 #include "atom_vec_tri.h"
 #include "atom_vec_body.h"
 #include "math_const.h"
 #include "math_extra.h"
 #include "memory.h"
 #include "error.h"
 
 using namespace LAMMPS_NS;
 using namespace FixConst;
 using namespace MathConst;
 
 enum{LINEAR,WIGGLE,ROTATE,VARIABLE};
 enum{EQUAL,ATOM};
 
 #define INERTIA 0.2          // moment of inertia prefactor for ellipsoid
 
 /* ---------------------------------------------------------------------- */
 
 FixMove::FixMove(LAMMPS *lmp, int narg, char **arg) :
   Fix(lmp, narg, arg),
-  xvarstr(NULL), yvarstr(NULL), zvarstr(NULL), vxvarstr(NULL), vyvarstr(NULL), vzvarstr(NULL),
-  xoriginal(NULL), toriginal(NULL), qoriginal(NULL), displace(NULL), velocity(NULL)
+  xvarstr(NULL), yvarstr(NULL), zvarstr(NULL), vxvarstr(NULL), 
+  vyvarstr(NULL), vzvarstr(NULL),
+  xoriginal(NULL), toriginal(NULL), qoriginal(NULL), 
+  displace(NULL), velocity(NULL)
 {
   if (narg < 4) error->all(FLERR,"Illegal fix move command");
 
   restart_global = 1;
   restart_peratom = 1;
   peratom_flag = 1;
   size_peratom_cols = 3;
   peratom_freq = 1;
   time_integrate = 1;
   create_attribute = 1;
   displaceflag = 0;
   velocityflag = 0;
   maxatom = 0;
 
   // parse args
 
   int iarg;
   
   if (strcmp(arg[3],"linear") == 0) {
     if (narg < 7) error->all(FLERR,"Illegal fix move command");
     iarg = 7;
     mstyle = LINEAR;
     if (strcmp(arg[4],"NULL") == 0) vxflag = 0;
     else {
       vxflag = 1;
       vx = force->numeric(FLERR,arg[4]);
     }
     if (strcmp(arg[5],"NULL") == 0) vyflag = 0;
     else {
       vyflag = 1;
       vy = force->numeric(FLERR,arg[5]);
     }
     if (strcmp(arg[6],"NULL") == 0) vzflag = 0;
     else {
       vzflag = 1;
       vz = force->numeric(FLERR,arg[6]);
     }
 
   } else if (strcmp(arg[3],"wiggle") == 0) {
     if (narg < 8) error->all(FLERR,"Illegal fix move command");
     iarg = 8;
     mstyle = WIGGLE;
     if (strcmp(arg[4],"NULL") == 0) axflag = 0;
     else {
       axflag = 1;
       ax = force->numeric(FLERR,arg[4]);
     }
     if (strcmp(arg[5],"NULL") == 0) ayflag = 0;
     else {
       ayflag = 1;
       ay = force->numeric(FLERR,arg[5]);
     }
     if (strcmp(arg[6],"NULL") == 0) azflag = 0;
     else {
       azflag = 1;
       az = force->numeric(FLERR,arg[6]);
     }
     period = force->numeric(FLERR,arg[7]);
     if (period <= 0.0) error->all(FLERR,"Illegal fix move command");
 
   } else if (strcmp(arg[3],"rotate") == 0) {
     if (narg < 11) error->all(FLERR,"Illegal fix move command");
     iarg = 11;
     mstyle = ROTATE;
     point[0] = force->numeric(FLERR,arg[4]);
     point[1] = force->numeric(FLERR,arg[5]);
     point[2] = force->numeric(FLERR,arg[6]);
     axis[0] = force->numeric(FLERR,arg[7]);
     axis[1] = force->numeric(FLERR,arg[8]);
     axis[2] = force->numeric(FLERR,arg[9]);
     period = force->numeric(FLERR,arg[10]);
     if (period <= 0.0) error->all(FLERR,"Illegal fix move command");
 
   } else if (strcmp(arg[3],"variable") == 0) {
     if (narg < 10) error->all(FLERR,"Illegal fix move command");
     iarg = 10;
     mstyle = VARIABLE;
     if (strcmp(arg[4],"NULL") == 0) xvarstr = NULL;
     else if (strstr(arg[4],"v_") == arg[4]) {
       int n = strlen(&arg[4][2]) + 1;
       xvarstr = new char[n];
       strcpy(xvarstr,&arg[4][2]);
     } else error->all(FLERR,"Illegal fix move command");
     if (strcmp(arg[5],"NULL") == 0) yvarstr = NULL;
     else if (strstr(arg[5],"v_") == arg[5]) {
       int n = strlen(&arg[5][2]) + 1;
       yvarstr = new char[n];
       strcpy(yvarstr,&arg[5][2]);
     } else error->all(FLERR,"Illegal fix move command");
     if (strcmp(arg[6],"NULL") == 0) zvarstr = NULL;
     else if (strstr(arg[6],"v_") == arg[6]) {
       int n = strlen(&arg[6][2]) + 1;
       zvarstr = new char[n];
       strcpy(zvarstr,&arg[6][2]);
     } else error->all(FLERR,"Illegal fix move command");
     if (strcmp(arg[7],"NULL") == 0) vxvarstr = NULL;
     else if (strstr(arg[7],"v_") == arg[7]) {
       int n = strlen(&arg[7][2]) + 1;
       vxvarstr = new char[n];
       strcpy(vxvarstr,&arg[7][2]);
     } else error->all(FLERR,"Illegal fix move command");
     if (strcmp(arg[8],"NULL") == 0) vyvarstr = NULL;
     else if (strstr(arg[8],"v_") == arg[8]) {
       int n = strlen(&arg[8][2]) + 1;
       vyvarstr = new char[n];
       strcpy(vyvarstr,&arg[8][2]);
     } else error->all(FLERR,"Illegal fix move command");
     if (strcmp(arg[9],"NULL") == 0) vzvarstr = NULL;
     else if (strstr(arg[9],"v_") == arg[9]) {
       int n = strlen(&arg[9][2]) + 1;
       vzvarstr = new char[n];
       strcpy(vzvarstr,&arg[9][2]);
     } else error->all(FLERR,"Illegal fix move command");
 
   } else error->all(FLERR,"Illegal fix move command");
 
   // optional args
 
   int scaleflag = 1;
 
   while (iarg < narg) {
     if (strcmp(arg[iarg],"units") == 0) {
       if (iarg+2 > narg) error->all(FLERR,"Illegal fix move command");
       if (strcmp(arg[iarg+1],"box") == 0) scaleflag = 0;
       else if (strcmp(arg[iarg+1],"lattice") == 0) scaleflag = 1;
       else error->all(FLERR,"Illegal fix move command");
       iarg += 2;
     } else error->all(FLERR,"Illegal fix move command");
   }
 
   // error checks and warnings
 
   if (domain->dimension == 2) {
     if (mstyle == LINEAR && vzflag && vz != 0.0)
       error->all(FLERR,"Fix move cannot set linear z motion for 2d problem");
     if (mstyle == WIGGLE && azflag && az != 0.0)
       error->all(FLERR,"Fix move cannot set wiggle z motion for 2d problem");
     if (mstyle == ROTATE && (axis[0] != 0.0 || axis[1] != 0.0))
       error->all(FLERR,
                  "Fix move cannot rotate around non z-axis for 2d problem");
     if (mstyle == VARIABLE && (zvarstr || vzvarstr))
       error->all(FLERR,
                  "Fix move cannot define z or vz variable for 2d problem");
   }
 
   // setup scaling and apply scaling factors to velocity & amplitude
 
   if ((mstyle == LINEAR || mstyle == WIGGLE || mstyle == ROTATE) &&
       scaleflag) {
     double xscale = domain->lattice->xlattice;
     double yscale = domain->lattice->ylattice;
     double zscale = domain->lattice->zlattice;
 
     if (mstyle == LINEAR) {
       if (vxflag) vx *= xscale;
       if (vyflag) vy *= yscale;
       if (vzflag) vz *= zscale;
     } else if (mstyle == WIGGLE) {
       if (axflag) ax *= xscale;
       if (ayflag) ay *= yscale;
       if (azflag) az *= zscale;
     } else if (mstyle == ROTATE) {
       point[0] *= xscale;
       point[1] *= yscale;
       point[2] *= zscale;
     }
   }
 
   // set omega_rotate from period
 
   if (mstyle == WIGGLE || mstyle == ROTATE) omega_rotate = MY_2PI / period;
 
   // runit = unit vector along rotation axis
 
   if (mstyle == ROTATE) {
     double len = sqrt(axis[0]*axis[0] + axis[1]*axis[1] + axis[2]*axis[2]);
     if (len == 0.0)
       error->all(FLERR,"Zero length rotation vector with fix move");
     runit[0] = axis[0]/len;
     runit[1] = axis[1]/len;
     runit[2] = axis[2]/len;
   }
 
   // set flags for extra attributes particles may store
   // relevant extra attributes = omega, angmom, theta, quat
   
   omega_flag = atom->omega_flag;
   angmom_flag = atom->angmom_flag;
 
   radius_flag = atom->radius_flag;
   ellipsoid_flag = atom->ellipsoid_flag;
   line_flag = atom->line_flag;
   tri_flag = atom->tri_flag;
   body_flag = atom->body_flag;
 
   theta_flag = quat_flag = 0;
   if (line_flag) theta_flag = 1;
   if (ellipsoid_flag || tri_flag || body_flag) quat_flag = 1;
 
   extra_flag = 0;
   if (omega_flag || angmom_flag || theta_flag || quat_flag) extra_flag = 1;
     
   // perform initial allocation of atom-based array
   // register with Atom class
 
   grow_arrays(atom->nmax);
   atom->add_callback(0);
   atom->add_callback(1);
 
   displace = velocity = NULL;
 
   // AtomVec pointers to retrieve per-atom storage of extra quantities
 
   avec_ellipsoid = (AtomVecEllipsoid *) atom->style_match("ellipsoid");
   avec_line = (AtomVecLine *) atom->style_match("line");
   avec_tri = (AtomVecTri *) atom->style_match("tri");
   avec_body = (AtomVecBody *) atom->style_match("body");
 
   // xoriginal = initial unwrapped positions of atoms
   // toriginal = initial theta of lines
   // qoriginal = initial quat of extended particles
 
   double **x = atom->x;
   imageint *image = atom->image;
   int *ellipsoid = atom->ellipsoid;
   int *line = atom->line;
   int *tri = atom->tri;
   int *body = atom->body;
   int *mask = atom->mask;
   int nlocal = atom->nlocal;
 
   for (int i = 0; i < nlocal; i++) {
     if (mask[i] & groupbit) domain->unmap(x[i],image[i],xoriginal[i]);
     else xoriginal[i][0] = xoriginal[i][1] = xoriginal[i][2] = 0.0;
   }
 
   if (theta_flag) {
     for (int i = 0; i < nlocal; i++) {
       if ((mask[i] & groupbit) && line[i] >= 0)
         toriginal[i] = avec_line->bonus[line[i]].theta;
       else toriginal[i] = 0.0;
     }
   }
 
   if (quat_flag) {
     double *quat;
     for (int i = 0; i < nlocal; i++) {
       quat = NULL;
       if (mask[i] & groupbit) {
         if (ellipsoid_flag && ellipsoid[i] >= 0)
           quat = avec_ellipsoid->bonus[ellipsoid[i]].quat;
         else if (tri_flag && tri[i] >= 0)
           quat = avec_tri->bonus[tri[i]].quat;
         else if (body_flag && body[i] >= 0)
           quat = avec_body->bonus[body[i]].quat;
       }
       if (quat) {
         qoriginal[i][0] = quat[0];
         qoriginal[i][1] = quat[1];
         qoriginal[i][2] = quat[2];
         qoriginal[i][3] = quat[3];
       } else qoriginal[i][0] = qoriginal[i][1] = 
                qoriginal[i][2] = qoriginal[i][3] = 0.0;
     }
   }
 
   // nrestart = size of per-atom restart data
   // nrestart = 1 + xorig + torig + qorig
 
   nrestart = 4;
   if (theta_flag) nrestart++;
   if (quat_flag) nrestart += 4;
 
   // time origin for movement = current timestep
 
   time_origin = update->ntimestep;
 }
 
 /* ---------------------------------------------------------------------- */
 
 FixMove::~FixMove()
 {
   // unregister callbacks to this fix from Atom class
 
   atom->delete_callback(id,0);
   atom->delete_callback(id,1);
 
   // delete locally stored arrays
 
   memory->destroy(xoriginal);
   memory->destroy(toriginal);
   memory->destroy(qoriginal);
   memory->destroy(displace);
   memory->destroy(velocity);
 
   delete [] xvarstr;
   delete [] yvarstr;
   delete [] zvarstr;
   delete [] vxvarstr;
   delete [] vyvarstr;
   delete [] vzvarstr;
 }
 
 /* ---------------------------------------------------------------------- */
 
 int FixMove::setmask()
 {
   int mask = 0;
   mask |= INITIAL_INTEGRATE;
   mask |= INITIAL_INTEGRATE_RESPA;
   mask |= FINAL_INTEGRATE;
   mask |= FINAL_INTEGRATE_RESPA;
   return mask;
 }
 
 /* ---------------------------------------------------------------------- */
 
 void FixMove::init()
 {
   dt = update->dt;
   dtv = update->dt;
   dtf = 0.5 * update->dt * force->ftm2v;
 
   // set indices and style of all variables
 
   displaceflag = velocityflag = 0;
 
   if (mstyle == VARIABLE) {
     if (xvarstr) {
       xvar = input->variable->find(xvarstr);
       if (xvar < 0) error->all(FLERR,
                                "Variable name for fix move does not exist");
       if (input->variable->equalstyle(xvar)) xvarstyle = EQUAL;
       else if (input->variable->atomstyle(xvar)) xvarstyle = ATOM;
       else error->all(FLERR,"Variable for fix move is invalid style");
     }
     if (yvarstr) {
       yvar = input->variable->find(yvarstr);
       if (yvar < 0) error->all(FLERR,
                                "Variable name for fix move does not exist");
       if (input->variable->equalstyle(yvar)) yvarstyle = EQUAL;
       else if (input->variable->atomstyle(yvar)) yvarstyle = ATOM;
       else error->all(FLERR,"Variable for fix move is invalid style");
     }
     if (zvarstr) {
       zvar = input->variable->find(zvarstr);
       if (zvar < 0) error->all(FLERR,
                                "Variable name for fix move does not exist");
       if (input->variable->equalstyle(zvar)) zvarstyle = EQUAL;
       else if (input->variable->atomstyle(zvar)) zvarstyle = ATOM;
       else error->all(FLERR,"Variable for fix move is invalid style");
     }
     if (vxvarstr) {
       vxvar = input->variable->find(vxvarstr);
       if (vxvar < 0) error->all(FLERR,
                                 "Variable name for fix move does not exist");
       if (input->variable->equalstyle(vxvar)) vxvarstyle = EQUAL;
       else if (input->variable->atomstyle(vxvar)) vxvarstyle = ATOM;
       else error->all(FLERR,"Variable for fix move is invalid style");
     }
     if (vyvarstr) {
       vyvar = input->variable->find(vyvarstr);
       if (vyvar < 0) error->all(FLERR,
                                 "Variable name for fix move does not exist");
       if (input->variable->equalstyle(vyvar)) vyvarstyle = EQUAL;
       else if (input->variable->atomstyle(vyvar)) vyvarstyle = ATOM;
       else error->all(FLERR,"Variable for fix move is invalid style");
     }
     if (vzvarstr) {
       vzvar = input->variable->find(vzvarstr);
       if (vzvar < 0) error->all(FLERR,
                                 "Variable name for fix move does not exist");
       if (input->variable->equalstyle(vzvar)) vzvarstyle = EQUAL;
       else if (input->variable->atomstyle(vzvar)) vzvarstyle = ATOM;
       else error->all(FLERR,"Variable for fix move is invalid style");
     }
 
     if (xvarstr && xvarstyle == ATOM) displaceflag = 1;
     if (yvarstr && yvarstyle == ATOM) displaceflag = 1;
     if (zvarstr && zvarstyle == ATOM) displaceflag = 1;
     if (vxvarstr && vxvarstyle == ATOM) velocityflag = 1;
     if (vyvarstr && vyvarstyle == ATOM) velocityflag = 1;
     if (vzvarstr && vzvarstyle == ATOM) velocityflag = 1;
   }
 
   maxatom = atom->nmax;
   memory->destroy(displace);
   memory->destroy(velocity);
   if (displaceflag) memory->create(displace,maxatom,3,"move:displace");
   else displace = NULL;
   if (velocityflag) memory->create(velocity,maxatom,3,"move:velocity");
   else velocity = NULL;
 
   if (strstr(update->integrate_style,"respa"))
     nlevels_respa = ((Respa *) update->integrate)->nlevels;
 }
 
 /* ----------------------------------------------------------------------
    set x,v of particles
 ------------------------------------------------------------------------- */
 
 void FixMove::initial_integrate(int vflag)
 {
   int flag;
   double ddotr,dx,dy,dz;
   double dtfm,theta_new;
   double xold[3],a[3],b[3],c[3],d[3],disp[3],w[3],ex[3],ey[3],ez[3];
   double inertia_ellipsoid[3],qrotate[4];
   double *quat,*inertia,*shape;
 
   double delta = (update->ntimestep - time_origin) * dt;
 
   double **x = atom->x;
   double **v = atom->v;
   double **f = atom->f;
   double **omega = atom->omega;
   double **angmom = atom->angmom;
   double *radius = atom->radius;
   double *rmass = atom->rmass;
   double *mass = atom->mass;
   int *type = atom->type;
   int *ellipsoid = atom->ellipsoid;
   int *line = atom->line;
   int *tri = atom->tri;
   int *body = atom->body;
   int *mask = atom->mask;
 
   int nlocal = atom->nlocal;
 
   // for linear: X = X0 + V*dt
 
   if (mstyle == LINEAR) {
     for (int i = 0; i < nlocal; i++) {
       if (mask[i] & groupbit) {
         xold[0] = x[i][0];
         xold[1] = x[i][1];
         xold[2] = x[i][2];
 
         if (vxflag) {
           v[i][0] = vx;
           x[i][0] = xoriginal[i][0] + vx*delta;
         } else if (rmass) {
           dtfm = dtf / rmass[i];
           v[i][0] += dtfm * f[i][0];
           x[i][0] += dtv * v[i][0];
         } else {
           dtfm = dtf / mass[type[i]];
           v[i][0] += dtfm * f[i][0];
           x[i][0] += dtv * v[i][0];
         }
 
         if (vyflag) {
           v[i][1] = vy;
           x[i][1] = xoriginal[i][1] + vy*delta;
         } else if (rmass) {
           dtfm = dtf / rmass[i];
           v[i][1] += dtfm * f[i][1];
           x[i][1] += dtv * v[i][1];
         } else {
           dtfm = dtf / mass[type[i]];
           v[i][1] += dtfm * f[i][1];
           x[i][1] += dtv * v[i][1];
         }
 
         if (vzflag) {
           v[i][2] = vz;
           x[i][2] = xoriginal[i][2] + vz*delta;
         } else if (rmass) {
           dtfm = dtf / rmass[i];
           v[i][2] += dtfm * f[i][2];
           x[i][2] += dtv * v[i][2];
         } else {
           dtfm = dtf / mass[type[i]];
           v[i][2] += dtfm * f[i][2];
           x[i][2] += dtv * v[i][2];
         }
 
         domain->remap_near(x[i],xold);
       }
     }
 
   // for wiggle: X = X0 + A sin(w*dt)
 
   } else if (mstyle == WIGGLE) {
     double arg = omega_rotate * delta;
     double sine = sin(arg);
     double cosine = cos(arg);
 
     for (int i = 0; i < nlocal; i++) {
       if (mask[i] & groupbit) {
         xold[0] = x[i][0];
         xold[1] = x[i][1];
         xold[2] = x[i][2];
 
         if (axflag) {
           v[i][0] = ax*omega_rotate*cosine;
           x[i][0] = xoriginal[i][0] + ax*sine;
         } else if (rmass) {
           dtfm = dtf / rmass[i];
           v[i][0] += dtfm * f[i][0];
           x[i][0] += dtv * v[i][0];
         } else {
           dtfm = dtf / mass[type[i]];
           v[i][0] += dtfm * f[i][0];
           x[i][0] += dtv * v[i][0];
         }
 
         if (ayflag) {
           v[i][1] = ay*omega_rotate*cosine;
           x[i][1] = xoriginal[i][1] + ay*sine;
         } else if (rmass) {
           dtfm = dtf / rmass[i];
           v[i][1] += dtfm * f[i][1];
           x[i][1] += dtv * v[i][1];
         } else {
           dtfm = dtf / mass[type[i]];
           v[i][1] += dtfm * f[i][1];
           x[i][1] += dtv * v[i][1];
         }
 
         if (azflag) {
           v[i][2] = az*omega_rotate*cosine;
           x[i][2] = xoriginal[i][2] + az*sine;
         } else if (rmass) {
           dtfm = dtf / rmass[i];
           v[i][2] += dtfm * f[i][2];
           x[i][2] += dtv * v[i][2];
         } else {
           dtfm = dtf / mass[type[i]];
           v[i][2] += dtfm * f[i][2];
           x[i][2] += dtv * v[i][2];
         }
 
         domain->remap_near(x[i],xold);
       }
     }
 
   // for rotate by right-hand rule around omega:
   // P = point = vector = point of rotation
   // R = vector = axis of rotation
   // w = omega of rotation (from period)
   // X0 = xoriginal = initial coord of atom
   // R0 = runit = unit vector for R
   // D = X0 - P = vector from P to X0
   // C = (D dot R0) R0 = projection of atom coord onto R line
   // A = D - C = vector from R line to X0
   // B = R0 cross A = vector perp to A in plane of rotation
   // A,B define plane of circular rotation around R line
   // X = P + C + A cos(w*dt) + B sin(w*dt)
   // V = w R0 cross (A cos(w*dt) + B sin(w*dt))
 
   } else if (mstyle == ROTATE) {
     double arg = omega_rotate * delta;
     double cosine = cos(arg);
     double sine = sin(arg);
 
     double qcosine = cos(0.5*arg);
     double qsine = sin(0.5*arg);
     qrotate[0] = qcosine;
     qrotate[1] = runit[0]*qsine;
     qrotate[2] = runit[1]*qsine;
     qrotate[3] = runit[2]*qsine;
 
     for (int i = 0; i < nlocal; i++) {
       if (mask[i] & groupbit) {
         xold[0] = x[i][0];
         xold[1] = x[i][1];
         xold[2] = x[i][2];
 
         d[0] = xoriginal[i][0] - point[0];
         d[1] = xoriginal[i][1] - point[1];
         d[2] = xoriginal[i][2] - point[2];
         ddotr = d[0]*runit[0] + d[1]*runit[1] + d[2]*runit[2];
         c[0] = ddotr*runit[0];
         c[1] = ddotr*runit[1];
         c[2] = ddotr*runit[2];
         a[0] = d[0] - c[0];
         a[1] = d[1] - c[1];
         a[2] = d[2] - c[2];
         b[0] = runit[1]*a[2] - runit[2]*a[1];
         b[1] = runit[2]*a[0] - runit[0]*a[2];
         b[2] = runit[0]*a[1] - runit[1]*a[0];
         disp[0] = a[0]*cosine  + b[0]*sine;
         disp[1] = a[1]*cosine  + b[1]*sine;
         disp[2] = a[2]*cosine  + b[2]*sine;
 
         x[i][0] = point[0] + c[0] + disp[0];
         x[i][1] = point[1] + c[1] + disp[1];
         x[i][2] = point[2] + c[2] + disp[2];
         v[i][0] = omega_rotate * (runit[1]*disp[2] - runit[2]*disp[1]);
         v[i][1] = omega_rotate * (runit[2]*disp[0] - runit[0]*disp[2]);
         v[i][2] = omega_rotate * (runit[0]*disp[1] - runit[1]*disp[0]);
 
         // set any extra attributes affected by rotation
 
         if (extra_flag) {
 
           // omega for spheres, lines, tris
 
           if (omega_flag) {
             flag = 0;
             if (radius_flag && radius[i] > 0.0) flag = 1;
             if (line_flag && line[i] >= 0.0) flag = 1;
             if (tri_flag && tri[i] >= 0.0) flag = 1;
             if (flag) {
               omega[i][0] = omega_rotate*runit[0];
               omega[i][1] = omega_rotate*runit[1];
               omega[i][2] = omega_rotate*runit[2];
             }
           }
 
           // angmom for ellipsoids, tris, and bodies
           
           if (angmom_flag) {
             quat = inertia = NULL;
             if (ellipsoid_flag && ellipsoid[i] >= 0) {
               quat = avec_ellipsoid->bonus[ellipsoid[i]].quat;
               shape = avec_ellipsoid->bonus[ellipsoid[i]].shape;
               inertia_ellipsoid[0] = 
                 INERTIA*rmass[i] * (shape[1]*shape[1]+shape[2]*shape[2]);
               inertia_ellipsoid[1] = 
                 INERTIA*rmass[i] * (shape[0]*shape[0]+shape[2]*shape[2]);
               inertia_ellipsoid[2] = 
                 INERTIA*rmass[i] * (shape[0]*shape[0]+shape[1]*shape[1]);
               inertia = inertia_ellipsoid;
             } else if (tri_flag && tri[i] >= 0) {
               quat = avec_tri->bonus[tri[i]].quat;
               inertia = avec_tri->bonus[tri[i]].inertia;
             } else if (body_flag && body[i] >= 0) {
               quat = avec_body->bonus[body[i]].quat;
               inertia = avec_body->bonus[body[i]].inertia;
             }
             if (quat) {
               w[0] = omega_rotate*runit[0];
               w[1] = omega_rotate*runit[1];
               w[2] = omega_rotate*runit[2];
               MathExtra::q_to_exyz(quat,ex,ey,ez);
               MathExtra::omega_to_angmom(w,ex,ey,ez,inertia,angmom[i]);
             }
           }
 
           // theta for lines
           
           if (theta_flag && line[i] >= 0.0) {
             theta_new = fmod(toriginal[i]+arg,MY_2PI);
             avec_line->bonus[atom->line[i]].theta = theta_new;
           }
           
           // quats for ellipsoids, tris, and bodies
           
           if (quat_flag) {
             quat = NULL;
             if (ellipsoid_flag && ellipsoid[i] >= 0)
               quat = avec_ellipsoid->bonus[ellipsoid[i]].quat;
             else if (tri_flag && tri[i] >= 0)
               quat = avec_tri->bonus[tri[i]].quat;
             else if (body_flag && body[i] >= 0)
               quat = avec_body->bonus[body[i]].quat;
             if (quat) MathExtra::quatquat(qrotate,qoriginal[i],quat);
           }
         }
 
         domain->remap_near(x[i],xold);
       }
     }
 
   // for variable: compute x,v from variables
   // NOTE: also allow for changes to extra attributes?
   //       omega, angmom, theta, quat
   //       only necessary if prescribed motion involves rotation
 
   } else if (mstyle == VARIABLE) {
 
     // reallocate displace and velocity arrays as necessary
 
     if ((displaceflag || velocityflag) && atom->nmax > maxatom) {
       maxatom = atom->nmax;
       if (displaceflag) {
         memory->destroy(displace);
         memory->create(displace,maxatom,3,"move:displace");
       }
       if (velocityflag) {
         memory->destroy(velocity);
         memory->create(velocity,maxatom,3,"move:velocity");
       }
     }
 
     // pre-compute variable values, wrap with clear/add
 
     modify->clearstep_compute();
 
     if (xvarstr) {
       if (xvarstyle == EQUAL) dx = input->variable->compute_equal(xvar);
       else input->variable->compute_atom(xvar,igroup,&displace[0][0],3,0);
     }
     if (yvarstr) {
       if (yvarstyle == EQUAL) dy = input->variable->compute_equal(yvar);
       else input->variable->compute_atom(yvar,igroup,&displace[0][1],3,0);
     }
     if (zvarstr) {
       if (zvarstyle == EQUAL) dz = input->variable->compute_equal(zvar);
       else input->variable->compute_atom(zvar,igroup,&displace[0][2],3,0);
     }
     if (vxvarstr) {
       if (vxvarstyle == EQUAL) vx = input->variable->compute_equal(vxvar);
       else input->variable->compute_atom(vxvar,igroup,&velocity[0][0],3,0);
     }
     if (vyvarstr) {
       if (vyvarstyle == EQUAL) vy = input->variable->compute_equal(vyvar);
       else input->variable->compute_atom(vyvar,igroup,&velocity[0][1],3,0);
     }
     if (vzvarstr) {
       if (vzvarstyle == EQUAL) vz = input->variable->compute_equal(vzvar);
       else input->variable->compute_atom(vzvar,igroup,&velocity[0][2],3,0);
     }
 
     modify->addstep_compute(update->ntimestep + 1);
 
     // update x,v
 
     for (int i = 0; i < nlocal; i++) {
       if (mask[i] & groupbit) {
         xold[0] = x[i][0];
         xold[1] = x[i][1];
         xold[2] = x[i][2];
 
         if (xvarstr && vxvarstr) {
           if (vxvarstyle == EQUAL) v[i][0] = vx;
           else v[i][0] = velocity[i][0];
           if (xvarstyle == EQUAL) x[i][0] = xoriginal[i][0] + dx;
           else x[i][0] = xoriginal[i][0] + displace[i][0];
         } else if (xvarstr) {
           if (xvarstyle == EQUAL) x[i][0] = xoriginal[i][0] + dx;
           else x[i][0] = xoriginal[i][0] + displace[i][0];
         } else if (vxvarstr) {
           if (vxvarstyle == EQUAL) v[i][0] = vx;
           else v[i][0] = velocity[i][0];
           if (rmass) {
             dtfm = dtf / rmass[i];
             x[i][0] += dtv * v[i][0];
           } else {
             dtfm = dtf / mass[type[i]];
             x[i][0] += dtv * v[i][0];
           }
         } else {
           if (rmass) {
             dtfm = dtf / rmass[i];
             v[i][0] += dtfm * f[i][0];
             x[i][0] += dtv * v[i][0];
           } else {
             dtfm = dtf / mass[type[i]];
             v[i][0] += dtfm * f[i][0];
             x[i][0] += dtv * v[i][0];
           }
         }
 
         if (yvarstr && vyvarstr) {
           if (vyvarstyle == EQUAL) v[i][1] = vy;
           else v[i][1] = velocity[i][1];
           if (yvarstyle == EQUAL) x[i][1] = xoriginal[i][1] + dy;
           else x[i][1] = xoriginal[i][1] + displace[i][1];
         } else if (yvarstr) {
           if (yvarstyle == EQUAL) x[i][1] = xoriginal[i][1] + dy;
           else x[i][1] = xoriginal[i][1] + displace[i][1];
         } else if (vyvarstr) {
           if (vyvarstyle == EQUAL) v[i][1] = vy;
           else v[i][1] = velocity[i][1];
           if (rmass) {
             dtfm = dtf / rmass[i];
             x[i][1] += dtv * v[i][1];
           } else {
             dtfm = dtf / mass[type[i]];
             x[i][1] += dtv * v[i][1];
           }
         } else {
           if (rmass) {
             dtfm = dtf / rmass[i];
             v[i][1] += dtfm * f[i][1];
             x[i][1] += dtv * v[i][1];
           } else {
             dtfm = dtf / mass[type[i]];
             v[i][1] += dtfm * f[i][1];
             x[i][1] += dtv * v[i][1];
           }
         }
 
         if (zvarstr && vzvarstr) {
           if (vzvarstyle == EQUAL) v[i][2] = vz;
           else v[i][2] = velocity[i][2];
           if (zvarstyle == EQUAL) x[i][2] = xoriginal[i][2] + dz;
           else x[i][2] = xoriginal[i][2] + displace[i][2];
         } else if (zvarstr) {
           if (zvarstyle == EQUAL) x[i][2] = xoriginal[i][2] + dz;
           else x[i][2] = xoriginal[i][2] + displace[i][2];
         } else if (vzvarstr) {
           if (vzvarstyle == EQUAL) v[i][2] = vz;
           else v[i][2] = velocity[i][2];
           if (rmass) {
             dtfm = dtf / rmass[i];
             x[i][2] += dtv * v[i][2];
           } else {
             dtfm = dtf / mass[type[i]];
             x[i][2] += dtv * v[i][2];
           }
         } else {
           if (rmass) {
             dtfm = dtf / rmass[i];
             v[i][2] += dtfm * f[i][2];
             x[i][2] += dtv * v[i][2];
           } else {
             dtfm = dtf / mass[type[i]];
             v[i][2] += dtfm * f[i][2];
             x[i][2] += dtv * v[i][2];
           }
         }
 
         domain->remap_near(x[i],xold);
       }
     }
   }
 }
 
 /* ----------------------------------------------------------------------
    final NVE of particles with NULL components
 ------------------------------------------------------------------------- */
 
 void FixMove::final_integrate()
 {
   double dtfm;
 
   int xflag = 1;
   if (mstyle == LINEAR && vxflag) xflag = 0;
   else if (mstyle == WIGGLE && axflag) xflag = 0;
   else if (mstyle == ROTATE) xflag = 0;
   else if (mstyle == VARIABLE && (xvarstr || vxvarstr)) xflag = 0;
 
   int yflag = 1;
   if (mstyle == LINEAR && vyflag) yflag = 0;
   else if (mstyle == WIGGLE && ayflag) yflag = 0;
   else if (mstyle == ROTATE) yflag = 0;
   else if (mstyle == VARIABLE && (yvarstr || vyvarstr)) yflag = 0;
 
   int zflag = 1;
   if (mstyle == LINEAR && vzflag) zflag = 0;
   else if (mstyle == WIGGLE && azflag) zflag = 0;
   else if (mstyle == ROTATE) zflag = 0;
   else if (mstyle == VARIABLE && (zvarstr || vzvarstr)) zflag = 0;
 
   double **v = atom->v;
   double **f = atom->f;
   double *rmass = atom->rmass;
   double *mass = atom->mass;
   int *type = atom->type;
   int *mask = atom->mask;
   int nlocal = atom->nlocal;
 
   for (int i = 0; i < nlocal; i++) {
     if (mask[i] & groupbit) {
       if (xflag) {
         if (rmass) {
           dtfm = dtf / rmass[i];
           v[i][0] += dtfm * f[i][0];
         } else {
           dtfm = dtf / mass[type[i]];
           v[i][0] += dtfm * f[i][0];
         }
       }
 
       if (yflag) {
         if (rmass) {
           dtfm = dtf / rmass[i];
           v[i][1] += dtfm * f[i][1];
         } else {
           dtfm = dtf / mass[type[i]];
           v[i][1] += dtfm * f[i][1];
         }
       }
 
       if (zflag) {
         if (rmass) {
           dtfm = dtf / rmass[i];
           v[i][2] += dtfm * f[i][2];
         } else {
           dtfm = dtf / mass[type[i]];
           v[i][2] += dtfm * f[i][2];
         }
       }
     }
   }
 }
 
 /* ---------------------------------------------------------------------- */
 
 void FixMove::initial_integrate_respa(int vflag, int ilevel, int iloop)
 {
   // outermost level - update v and x
   // all other levels - nothing
 
   if (ilevel == nlevels_respa-1) initial_integrate(vflag);
 }
 
 /* ---------------------------------------------------------------------- */
 
 void FixMove::final_integrate_respa(int ilevel, int iloop)
 {
   if (ilevel == nlevels_respa-1) final_integrate();
 }
 
 /* ----------------------------------------------------------------------
    memory usage of local atom-based array
 ------------------------------------------------------------------------- */
 
 double FixMove::memory_usage()
 {
   double bytes = atom->nmax*3 * sizeof(double);
   if (theta_flag) bytes += atom->nmax * sizeof(double);
   if (quat_flag) bytes += atom->nmax*4 * sizeof(double);
   if (displaceflag) bytes += atom->nmax*3 * sizeof(double);
   if (velocityflag) bytes += atom->nmax*3 * sizeof(double);
   return bytes;
 }
 
 /* ----------------------------------------------------------------------
    pack entire state of Fix into one write
 ------------------------------------------------------------------------- */
 
 void FixMove::write_restart(FILE *fp)
 {
   int n = 0;
   double list[1];
   list[n++] = time_origin;
 
   if (comm->me == 0) {
     int size = n * sizeof(double);
     fwrite(&size,sizeof(int),1,fp);
     fwrite(list,sizeof(double),n,fp);
   }
 }
 
 /* ----------------------------------------------------------------------
    use state info from restart file to restart the Fix
 ------------------------------------------------------------------------- */
 
 void FixMove::restart(char *buf)
 {
   int n = 0;
   double *list = (double *) buf;
 
   time_origin = static_cast<int> (list[n++]);
 }
 
 /* ----------------------------------------------------------------------
    allocate atom-based array
 ------------------------------------------------------------------------- */
 
 void FixMove::grow_arrays(int nmax)
 {
   memory->grow(xoriginal,nmax,3,"move:xoriginal");
   if (theta_flag) memory->grow(toriginal,nmax,"move:toriginal");
   if (quat_flag) memory->grow(qoriginal,nmax,4,"move:qoriginal");
   array_atom = xoriginal;
 }
 
 /* ----------------------------------------------------------------------
    copy values within local atom-based array
 ------------------------------------------------------------------------- */
 
 void FixMove::copy_arrays(int i, int j, int delflag)
 {
   xoriginal[j][0] = xoriginal[i][0];
   xoriginal[j][1] = xoriginal[i][1];
   xoriginal[j][2] = xoriginal[i][2];
   if (theta_flag) toriginal[j] = toriginal[i];
   if (quat_flag) {
     qoriginal[j][0] = qoriginal[i][0];
     qoriginal[j][1] = qoriginal[i][1];
     qoriginal[j][2] = qoriginal[i][2];
     qoriginal[j][3] = qoriginal[i][3];
   }
 }
 
 /* ----------------------------------------------------------------------
    initialize one atom's array values, called when atom is created
 ------------------------------------------------------------------------- */
 
 void FixMove::set_arrays(int i)
 {
   double theta;
   double *quat;
 
   double **x = atom->x;
   imageint *image = atom->image;
   int *ellipsoid = atom->ellipsoid;
   int *line = atom->line;
   int *tri = atom->tri;
   int *body = atom->body;
   int *mask = atom->mask;
 
   // particle not in group
 
   if (!(mask[i] & groupbit)) {
     xoriginal[i][0] = xoriginal[i][1] = xoriginal[i][2] = 0.0;
     return;
   }
 
   // current time still equal fix creation time
 
   if (update->ntimestep == time_origin) {
     domain->unmap(x[i],image[i],xoriginal[i]);
     return;
   }
 
   // backup particle to time_origin
 
   if (mstyle == VARIABLE)
     error->all(FLERR,"Cannot add atoms to fix move variable");
 
   domain->unmap(x[i],image[i],xoriginal[i]);
   double delta = (update->ntimestep - time_origin) * update->dt;
 
   if (mstyle == LINEAR) {
     if (vxflag) xoriginal[i][0] -= vx * delta;
     if (vyflag) xoriginal[i][1] -= vy * delta;
     if (vzflag) xoriginal[i][2] -= vz * delta;
   } else if (mstyle == WIGGLE) {
     double arg = omega_rotate * delta;
     double sine = sin(arg);
     if (axflag) xoriginal[i][0] -= ax*sine;
     if (ayflag) xoriginal[i][1] -= ay*sine;
     if (azflag) xoriginal[i][2] -= az*sine;
   } else if (mstyle == ROTATE) {
     double a[3],b[3],c[3],d[3],disp[3],ddotr;
     double arg = - omega_rotate * delta;
     double sine = sin(arg);
     double cosine = cos(arg);
     d[0] = x[i][0] - point[0];
     d[1] = x[i][1] - point[1];
     d[2] = x[i][2] - point[2];
     ddotr = d[0]*runit[0] + d[1]*runit[1] + d[2]*runit[2];
     c[0] = ddotr*runit[0];
     c[1] = ddotr*runit[1];
     c[2] = ddotr*runit[2];
 
     a[0] = d[0] - c[0];
     a[1] = d[1] - c[1];
     a[2] = d[2] - c[2];
     b[0] = runit[1]*a[2] - runit[2]*a[1];
     b[1] = runit[2]*a[0] - runit[0]*a[2];
     b[2] = runit[0]*a[1] - runit[1]*a[0];
     disp[0] = a[0]*cosine  + b[0]*sine;
     disp[1] = a[1]*cosine  + b[1]*sine;
     disp[2] = a[2]*cosine  + b[2]*sine;
 
     xoriginal[i][0] = point[0] + c[0] + disp[0];
     xoriginal[i][1] = point[1] + c[1] + disp[1];
     xoriginal[i][2] = point[2] + c[2] + disp[2];
     
     // set theta and quat extra attributes affected by rotation
 
     if (extra_flag) {
 
       // theta for lines
       
       if (theta_flag && line[i] >= 0.0) {
         theta = avec_line->bonus[atom->line[i]].theta;
         toriginal[i] = theta - 0.0;  // NOTE: edit this line
       }
       
       // quats for ellipsoids, tris, and bodies
       
       if (quat_flag) {
         quat = NULL;
         if (ellipsoid_flag && ellipsoid[i] >= 0)
           quat = avec_ellipsoid->bonus[ellipsoid[i]].quat;
         else if (tri_flag && tri[i] >= 0)
           quat = avec_tri->bonus[tri[i]].quat;
         else if (body_flag && body[i] >= 0)
           quat = avec_body->bonus[body[i]].quat;
         if (quat) {
           // qoriginal = f(quat,-delta);   // NOTE: edit this line
         }
       }
     }
   }
 }
 
 /* ----------------------------------------------------------------------
    pack values in local atom-based array for exchange with another proc
 ------------------------------------------------------------------------- */
 
 int FixMove::pack_exchange(int i, double *buf)
 {
   int n = 0;
   buf[n++] = xoriginal[i][0];
   buf[n++] = xoriginal[i][1];
   buf[n++] = xoriginal[i][2];
   if (theta_flag) buf[n++] = toriginal[i];
   if (quat_flag) {
     buf[n++] = qoriginal[i][0];
     buf[n++] = qoriginal[i][1];
     buf[n++] = qoriginal[i][2];
     buf[n++] = qoriginal[i][3];
   }
   return n;
 }
 
 /* ----------------------------------------------------------------------
    unpack values in local atom-based array from exchange with another proc
 ------------------------------------------------------------------------- */
 
 int FixMove::unpack_exchange(int nlocal, double *buf)
 {
   int n = 0;
   xoriginal[nlocal][0] = buf[n++];
   xoriginal[nlocal][1] = buf[n++];
   xoriginal[nlocal][2] = buf[n++];
   if (theta_flag) toriginal[nlocal] = buf[n++];
   if (quat_flag) {
     qoriginal[nlocal][0] = buf[n++];
     qoriginal[nlocal][1] = buf[n++];
     qoriginal[nlocal][2] = buf[n++];
     qoriginal[nlocal][3] = buf[n++];
   }
   return n;
 }
 
 /* ----------------------------------------------------------------------
    pack values in local atom-based arrays for restart file
 ------------------------------------------------------------------------- */
 
 int FixMove::pack_restart(int i, double *buf)
 {
   int n = 1;
   buf[n++] = xoriginal[i][0];
   buf[n++] = xoriginal[i][1];
   buf[n++] = xoriginal[i][2];
   if (theta_flag) buf[n++] = toriginal[i];
   if (quat_flag) {
     buf[n++] = qoriginal[i][0];
     buf[n++] = qoriginal[i][1];
     buf[n++] = qoriginal[i][2];
     buf[n++] = qoriginal[i][3];
   }
   buf[0] = n;
   return n;
 }
 
 /* ----------------------------------------------------------------------
    unpack values from atom->extra array to restart the fix
 ------------------------------------------------------------------------- */
 
 void FixMove::unpack_restart(int nlocal, int nth)
 {
   double **extra = atom->extra;
 
   // skip to Nth set of extra values
 
   int m = 0;
   for (int i = 0; i < nth; i++) m += static_cast<int> (extra[nlocal][m]);
   m++;
 
   xoriginal[nlocal][0] = extra[nlocal][m++];
   xoriginal[nlocal][1] = extra[nlocal][m++];
   xoriginal[nlocal][2] = extra[nlocal][m++];
   if (theta_flag) toriginal[nlocal] = extra[nlocal][m++];
   if (quat_flag) {
     qoriginal[nlocal][0] = extra[nlocal][m++];
     qoriginal[nlocal][1] = extra[nlocal][m++];
     qoriginal[nlocal][2] = extra[nlocal][m++];
     qoriginal[nlocal][3] = extra[nlocal][m++];
   }
 }
 
 /* ----------------------------------------------------------------------
    maxsize of any atom's restart data
 ------------------------------------------------------------------------- */
 
 int FixMove::maxsize_restart()
 {
   return nrestart;
 }
 
 /* ----------------------------------------------------------------------
    size of atom nlocal's restart data
 ------------------------------------------------------------------------- */
 
 int FixMove::size_restart(int nlocal)
 {
   return nrestart;
 }
 
 /* ---------------------------------------------------------------------- */
 
 void FixMove::reset_dt()
 {
   error->all(FLERR,"Resetting timestep size is not allowed with fix move");
 }
diff --git a/src/fix_nh.cpp b/src/fix_nh.cpp
index e2ce97db4..01ba81e55 100644
--- a/src/fix_nh.cpp
+++ b/src/fix_nh.cpp
@@ -1,2336 +1,2337 @@
 /* ----------------------------------------------------------------------
    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: Mark Stevens (SNL), Aidan Thompson (SNL)
 ------------------------------------------------------------------------- */
 
 #include <string.h>
 #include <stdlib.h>
 #include <math.h>
 #include "fix_nh.h"
 #include "math_extra.h"
 #include "atom.h"
 #include "force.h"
 #include "group.h"
 #include "comm.h"
 #include "neighbor.h"
 #include "irregular.h"
 #include "modify.h"
 #include "fix_deform.h"
 #include "compute.h"
 #include "kspace.h"
 #include "update.h"
 #include "respa.h"
 #include "domain.h"
 #include "memory.h"
 #include "error.h"
 
 using namespace LAMMPS_NS;
 using namespace FixConst;
 
 #define DELTAFLIP 0.1
 #define TILTMAX 1.5
 
 enum{NOBIAS,BIAS};
 enum{NONE,XYZ,XY,YZ,XZ};
 enum{ISO,ANISO,TRICLINIC};
 
 /* ----------------------------------------------------------------------
    NVT,NPH,NPT integrators for improved Nose-Hoover equations of motion
  ---------------------------------------------------------------------- */
 
 FixNH::FixNH(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg),
 tstat_flag(0), pstat_flag(0),
 rfix(NULL), id_dilate(NULL), irregular(NULL), id_temp(NULL), id_press(NULL),
 tcomputeflag(0), pcomputeflag(0), eta(NULL), eta_dot(NULL), eta_dotdot(NULL),
-eta_mass(NULL), etap(NULL), etap_dot(NULL), etap_dotdot(NULL), etap_mass(NULL), mpchain(0)
+eta_mass(NULL), etap(NULL), etap_dot(NULL), etap_dotdot(NULL), 
+etap_mass(NULL), mpchain(0)
 {
   if (narg < 4) error->all(FLERR,"Illegal fix nvt/npt/nph command");
 
   restart_global = 1;
   dynamic_group_allow = 1;
   time_integrate = 1;
   scalar_flag = 1;
   vector_flag = 1;
   global_freq = 1;
   extscalar = 1;
   extvector = 0;
 
   // default values
 
   pcouple = NONE;
   drag = 0.0;
   allremap = 1;
   id_dilate = NULL;
   mtchain = mpchain = 3;
   nc_tchain = nc_pchain = 1;
   mtk_flag = 1;
   deviatoric_flag = 0;
   nreset_h0 = 0;
   eta_mass_flag = 1;
   omega_mass_flag = 0;
   etap_mass_flag = 0;
   flipflag = 1;
   dipole_flag = 0;
   dlm_flag = 0;
 
   tcomputeflag = 0;
   pcomputeflag = 0;
   id_temp = NULL;
   id_press = NULL;
 
   // turn on tilt factor scaling, whenever applicable
 
   dimension = domain->dimension;
 
   scaleyz = scalexz = scalexy = 0;
   if (domain->yperiodic && domain->xy != 0.0) scalexy = 1;
   if (domain->zperiodic && dimension == 3) {
     if (domain->yz != 0.0) scaleyz = 1;
     if (domain->xz != 0.0) scalexz = 1;
   }
 
   // set fixed-point to default = center of cell
 
   fixedpoint[0] = 0.5*(domain->boxlo[0]+domain->boxhi[0]);
   fixedpoint[1] = 0.5*(domain->boxlo[1]+domain->boxhi[1]);
   fixedpoint[2] = 0.5*(domain->boxlo[2]+domain->boxhi[2]);
 
   // used by FixNVTSllod to preserve non-default value
 
   mtchain_default_flag = 1;
 
   tstat_flag = 0;
   double t_period = 0.0;
 
   double p_period[6];
   for (int i = 0; i < 6; i++) {
     p_start[i] = p_stop[i] = p_period[i] = p_target[i] = 0.0;
     p_flag[i] = 0;
   }
 
   // process keywords
 
   int iarg = 3;
 
   while (iarg < narg) {
     if (strcmp(arg[iarg],"temp") == 0) {
       if (iarg+4 > narg) error->all(FLERR,"Illegal fix nvt/npt/nph command");
       tstat_flag = 1;
       t_start = force->numeric(FLERR,arg[iarg+1]);
       t_target = t_start;
       t_stop = force->numeric(FLERR,arg[iarg+2]);
       t_period = force->numeric(FLERR,arg[iarg+3]);
       if (t_start < 0.0 || t_stop <= 0.0)
         error->all(FLERR,
                    "Target temperature for fix nvt/npt/nph cannot be 0.0");
       iarg += 4;
 
     } else if (strcmp(arg[iarg],"iso") == 0) {
       if (iarg+4 > narg) error->all(FLERR,"Illegal fix nvt/npt/nph command");
       pcouple = XYZ;
       p_start[0] = p_start[1] = p_start[2] = force->numeric(FLERR,arg[iarg+1]);
       p_stop[0] = p_stop[1] = p_stop[2] = force->numeric(FLERR,arg[iarg+2]);
       p_period[0] = p_period[1] = p_period[2] =
         force->numeric(FLERR,arg[iarg+3]);
       p_flag[0] = p_flag[1] = p_flag[2] = 1;
       if (dimension == 2) {
         p_start[2] = p_stop[2] = p_period[2] = 0.0;
         p_flag[2] = 0;
       }
       iarg += 4;
     } else if (strcmp(arg[iarg],"aniso") == 0) {
       if (iarg+4 > narg) error->all(FLERR,"Illegal fix nvt/npt/nph command");
       pcouple = NONE;
       p_start[0] = p_start[1] = p_start[2] = force->numeric(FLERR,arg[iarg+1]);
       p_stop[0] = p_stop[1] = p_stop[2] = force->numeric(FLERR,arg[iarg+2]);
       p_period[0] = p_period[1] = p_period[2] =
         force->numeric(FLERR,arg[iarg+3]);
       p_flag[0] = p_flag[1] = p_flag[2] = 1;
       if (dimension == 2) {
         p_start[2] = p_stop[2] = p_period[2] = 0.0;
         p_flag[2] = 0;
       }
       iarg += 4;
     } else if (strcmp(arg[iarg],"tri") == 0) {
       if (iarg+4 > narg) error->all(FLERR,"Illegal fix nvt/npt/nph command");
       pcouple = NONE;
       scalexy = scalexz = scaleyz = 0;
       p_start[0] = p_start[1] = p_start[2] = force->numeric(FLERR,arg[iarg+1]);
       p_stop[0] = p_stop[1] = p_stop[2] = force->numeric(FLERR,arg[iarg+2]);
       p_period[0] = p_period[1] = p_period[2] =
         force->numeric(FLERR,arg[iarg+3]);
       p_flag[0] = p_flag[1] = p_flag[2] = 1;
       p_start[3] = p_start[4] = p_start[5] = 0.0;
       p_stop[3] = p_stop[4] = p_stop[5] = 0.0;
       p_period[3] = p_period[4] = p_period[5] =
         force->numeric(FLERR,arg[iarg+3]);
       p_flag[3] = p_flag[4] = p_flag[5] = 1;
       if (dimension == 2) {
         p_start[2] = p_stop[2] = p_period[2] = 0.0;
         p_flag[2] = 0;
         p_start[3] = p_stop[3] = p_period[3] = 0.0;
         p_flag[3] = 0;
         p_start[4] = p_stop[4] = p_period[4] = 0.0;
         p_flag[4] = 0;
       }
       iarg += 4;
     } else if (strcmp(arg[iarg],"x") == 0) {
       if (iarg+4 > narg) error->all(FLERR,"Illegal fix nvt/npt/nph command");
       p_start[0] = force->numeric(FLERR,arg[iarg+1]);
       p_stop[0] = force->numeric(FLERR,arg[iarg+2]);
       p_period[0] = force->numeric(FLERR,arg[iarg+3]);
       p_flag[0] = 1;
       deviatoric_flag = 1;
       iarg += 4;
     } else if (strcmp(arg[iarg],"y") == 0) {
       if (iarg+4 > narg) error->all(FLERR,"Illegal fix nvt/npt/nph command");
       p_start[1] = force->numeric(FLERR,arg[iarg+1]);
       p_stop[1] = force->numeric(FLERR,arg[iarg+2]);
       p_period[1] = force->numeric(FLERR,arg[iarg+3]);
       p_flag[1] = 1;
       deviatoric_flag = 1;
       iarg += 4;
     } else if (strcmp(arg[iarg],"z") == 0) {
       if (iarg+4 > narg) error->all(FLERR,"Illegal fix nvt/npt/nph command");
       p_start[2] = force->numeric(FLERR,arg[iarg+1]);
       p_stop[2] = force->numeric(FLERR,arg[iarg+2]);
       p_period[2] = force->numeric(FLERR,arg[iarg+3]);
       p_flag[2] = 1;
       deviatoric_flag = 1;
       iarg += 4;
       if (dimension == 2)
         error->all(FLERR,"Invalid fix nvt/npt/nph command for a 2d simulation");
 
     } else if (strcmp(arg[iarg],"yz") == 0) {
       if (iarg+4 > narg) error->all(FLERR,"Illegal fix nvt/npt/nph command");
       p_start[3] = force->numeric(FLERR,arg[iarg+1]);
       p_stop[3] = force->numeric(FLERR,arg[iarg+2]);
       p_period[3] = force->numeric(FLERR,arg[iarg+3]);
       p_flag[3] = 1;
       deviatoric_flag = 1;
       scaleyz = 0;
       iarg += 4;
       if (dimension == 2)
         error->all(FLERR,"Invalid fix nvt/npt/nph command for a 2d simulation");
     } else if (strcmp(arg[iarg],"xz") == 0) {
       if (iarg+4 > narg) error->all(FLERR,"Illegal fix nvt/npt/nph command");
       p_start[4] = force->numeric(FLERR,arg[iarg+1]);
       p_stop[4] = force->numeric(FLERR,arg[iarg+2]);
       p_period[4] = force->numeric(FLERR,arg[iarg+3]);
       p_flag[4] = 1;
       deviatoric_flag = 1;
       scalexz = 0;
       iarg += 4;
       if (dimension == 2)
         error->all(FLERR,"Invalid fix nvt/npt/nph command for a 2d simulation");
     } else if (strcmp(arg[iarg],"xy") == 0) {
       if (iarg+4 > narg) error->all(FLERR,"Illegal fix nvt/npt/nph command");
       p_start[5] = force->numeric(FLERR,arg[iarg+1]);
       p_stop[5] = force->numeric(FLERR,arg[iarg+2]);
       p_period[5] = force->numeric(FLERR,arg[iarg+3]);
       p_flag[5] = 1;
       deviatoric_flag = 1;
       scalexy = 0;
       iarg += 4;
 
     } else if (strcmp(arg[iarg],"couple") == 0) {
       if (iarg+2 > narg) error->all(FLERR,"Illegal fix nvt/npt/nph command");
       if (strcmp(arg[iarg+1],"xyz") == 0) pcouple = XYZ;
       else if (strcmp(arg[iarg+1],"xy") == 0) pcouple = XY;
       else if (strcmp(arg[iarg+1],"yz") == 0) pcouple = YZ;
       else if (strcmp(arg[iarg+1],"xz") == 0) pcouple = XZ;
       else if (strcmp(arg[iarg+1],"none") == 0) pcouple = NONE;
       else error->all(FLERR,"Illegal fix nvt/npt/nph command");
       iarg += 2;
 
     } else if (strcmp(arg[iarg],"drag") == 0) {
       if (iarg+2 > narg) error->all(FLERR,"Illegal fix nvt/npt/nph command");
       drag = force->numeric(FLERR,arg[iarg+1]);
       if (drag < 0.0) error->all(FLERR,"Illegal fix nvt/npt/nph command");
       iarg += 2;
     } else if (strcmp(arg[iarg],"dilate") == 0) {
       if (iarg+2 > narg) error->all(FLERR,"Illegal fix nvt/npt/nph command");
       if (strcmp(arg[iarg+1],"all") == 0) allremap = 1;
       else {
         allremap = 0;
         delete [] id_dilate;
         int n = strlen(arg[iarg+1]) + 1;
         id_dilate = new char[n];
         strcpy(id_dilate,arg[iarg+1]);
         int idilate = group->find(id_dilate);
         if (idilate == -1)
           error->all(FLERR,"Fix nvt/npt/nph dilate group ID does not exist");
       }
       iarg += 2;
 
     } else if (strcmp(arg[iarg],"tchain") == 0) {
       if (iarg+2 > narg) error->all(FLERR,"Illegal fix nvt/npt/nph command");
       mtchain = force->inumeric(FLERR,arg[iarg+1]);
       // used by FixNVTSllod to preserve non-default value
       mtchain_default_flag = 0;
       if (mtchain < 1) error->all(FLERR,"Illegal fix nvt/npt/nph command");
       iarg += 2;
     } else if (strcmp(arg[iarg],"pchain") == 0) {
       if (iarg+2 > narg) error->all(FLERR,"Illegal fix nvt/npt/nph command");
       mpchain = force->inumeric(FLERR,arg[iarg+1]);
       if (mpchain < 0) error->all(FLERR,"Illegal fix nvt/npt/nph command");
       iarg += 2;
     } else if (strcmp(arg[iarg],"mtk") == 0) {
       if (iarg+2 > narg) error->all(FLERR,"Illegal fix nvt/npt/nph command");
       if (strcmp(arg[iarg+1],"yes") == 0) mtk_flag = 1;
       else if (strcmp(arg[iarg+1],"no") == 0) mtk_flag = 0;
       else error->all(FLERR,"Illegal fix nvt/npt/nph command");
       iarg += 2;
     } else if (strcmp(arg[iarg],"tloop") == 0) {
       if (iarg+2 > narg) error->all(FLERR,"Illegal fix nvt/npt/nph command");
       nc_tchain = force->inumeric(FLERR,arg[iarg+1]);
       if (nc_tchain < 0) error->all(FLERR,"Illegal fix nvt/npt/nph command");
       iarg += 2;
     } else if (strcmp(arg[iarg],"ploop") == 0) {
       if (iarg+2 > narg) error->all(FLERR,"Illegal fix nvt/npt/nph command");
       nc_pchain = force->inumeric(FLERR,arg[iarg+1]);
       if (nc_pchain < 0) error->all(FLERR,"Illegal fix nvt/npt/nph command");
       iarg += 2;
     } else if (strcmp(arg[iarg],"nreset") == 0) {
       if (iarg+2 > narg) error->all(FLERR,"Illegal fix nvt/npt/nph command");
       nreset_h0 = force->inumeric(FLERR,arg[iarg+1]);
       if (nreset_h0 < 0) error->all(FLERR,"Illegal fix nvt/npt/nph command");
       iarg += 2;
     } else if (strcmp(arg[iarg],"scalexy") == 0) {
       if (iarg+2 > narg) error->all(FLERR,"Illegal fix nvt/npt/nph command");
       if (strcmp(arg[iarg+1],"yes") == 0) scalexy = 1;
       else if (strcmp(arg[iarg+1],"no") == 0) scalexy = 0;
       else error->all(FLERR,"Illegal fix nvt/npt/nph command");
       iarg += 2;
     } else if (strcmp(arg[iarg],"scalexz") == 0) {
       if (iarg+2 > narg) error->all(FLERR,"Illegal fix nvt/npt/nph command");
       if (strcmp(arg[iarg+1],"yes") == 0) scalexz = 1;
       else if (strcmp(arg[iarg+1],"no") == 0) scalexz = 0;
       else error->all(FLERR,"Illegal fix nvt/npt/nph command");
       iarg += 2;
     } else if (strcmp(arg[iarg],"scaleyz") == 0) {
       if (iarg+2 > narg) error->all(FLERR,"Illegal fix nvt/npt/nph command");
       if (strcmp(arg[iarg+1],"yes") == 0) scaleyz = 1;
       else if (strcmp(arg[iarg+1],"no") == 0) scaleyz = 0;
       else error->all(FLERR,"Illegal fix nvt/npt/nph command");
       iarg += 2;
     } else if (strcmp(arg[iarg],"flip") == 0) {
       if (iarg+2 > narg) error->all(FLERR,"Illegal fix nvt/npt/nph command");
       if (strcmp(arg[iarg+1],"yes") == 0) flipflag = 1;
       else if (strcmp(arg[iarg+1],"no") == 0) flipflag = 0;
       else error->all(FLERR,"Illegal fix nvt/npt/nph command");
       iarg += 2;
     } else if (strcmp(arg[iarg],"update") == 0) {
       if (iarg+2 > narg) error->all(FLERR,"Illegal fix nvt/npt/nph command");
       if (strcmp(arg[iarg+1],"dipole") == 0) dipole_flag = 1;
       else if (strcmp(arg[iarg+1],"dipole/dlm") == 0) {
         dipole_flag = 1;
         dlm_flag = 1;
       } else error->all(FLERR,"Illegal fix nvt/npt/nph command");
       iarg += 2;
     } else if (strcmp(arg[iarg],"fixedpoint") == 0) {
       if (iarg+4 > narg) error->all(FLERR,"Illegal fix nvt/npt/nph command");
       fixedpoint[0] = force->numeric(FLERR,arg[iarg+1]);
       fixedpoint[1] = force->numeric(FLERR,arg[iarg+2]);
       fixedpoint[2] = force->numeric(FLERR,arg[iarg+3]);
       iarg += 4;
     } else error->all(FLERR,"Illegal fix nvt/npt/nph command");
   }
 
   // error checks
 
   if (dimension == 2 && (p_flag[2] || p_flag[3] || p_flag[4]))
     error->all(FLERR,"Invalid fix nvt/npt/nph command for a 2d simulation");
   if (dimension == 2 && (pcouple == YZ || pcouple == XZ))
     error->all(FLERR,"Invalid fix nvt/npt/nph command for a 2d simulation");
   if (dimension == 2 && (scalexz == 1 || scaleyz == 1 ))
     error->all(FLERR,"Invalid fix nvt/npt/nph command for a 2d simulation");
 
   if (pcouple == XYZ && (p_flag[0] == 0 || p_flag[1] == 0))
     error->all(FLERR,"Invalid fix nvt/npt/nph command pressure settings");
   if (pcouple == XYZ && dimension == 3 && p_flag[2] == 0)
     error->all(FLERR,"Invalid fix nvt/npt/nph command pressure settings");
   if (pcouple == XY && (p_flag[0] == 0 || p_flag[1] == 0))
     error->all(FLERR,"Invalid fix nvt/npt/nph command pressure settings");
   if (pcouple == YZ && (p_flag[1] == 0 || p_flag[2] == 0))
     error->all(FLERR,"Invalid fix nvt/npt/nph command pressure settings");
   if (pcouple == XZ && (p_flag[0] == 0 || p_flag[2] == 0))
     error->all(FLERR,"Invalid fix nvt/npt/nph command pressure settings");
 
   // require periodicity in tensile dimension
 
   if (p_flag[0] && domain->xperiodic == 0)
     error->all(FLERR,"Cannot use fix nvt/npt/nph on a non-periodic dimension");
   if (p_flag[1] && domain->yperiodic == 0)
     error->all(FLERR,"Cannot use fix nvt/npt/nph on a non-periodic dimension");
   if (p_flag[2] && domain->zperiodic == 0)
     error->all(FLERR,"Cannot use fix nvt/npt/nph on a non-periodic dimension");
 
   // require periodicity in 2nd dim of off-diagonal tilt component
 
   if (p_flag[3] && domain->zperiodic == 0)
     error->all(FLERR,
                "Cannot use fix nvt/npt/nph on a 2nd non-periodic dimension");
   if (p_flag[4] && domain->zperiodic == 0)
     error->all(FLERR,
                "Cannot use fix nvt/npt/nph on a 2nd non-periodic dimension");
   if (p_flag[5] && domain->yperiodic == 0)
     error->all(FLERR,
                "Cannot use fix nvt/npt/nph on a 2nd non-periodic dimension");
 
   if (scaleyz == 1 && domain->zperiodic == 0)
     error->all(FLERR,"Cannot use fix nvt/npt/nph "
                "with yz scaling when z is non-periodic dimension");
   if (scalexz == 1 && domain->zperiodic == 0)
     error->all(FLERR,"Cannot use fix nvt/npt/nph "
                "with xz scaling when z is non-periodic dimension");
   if (scalexy == 1 && domain->yperiodic == 0)
     error->all(FLERR,"Cannot use fix nvt/npt/nph "
                "with xy scaling when y is non-periodic dimension");
 
   if (p_flag[3] && scaleyz == 1)
     error->all(FLERR,"Cannot use fix nvt/npt/nph with "
                "both yz dynamics and yz scaling");
   if (p_flag[4] && scalexz == 1)
     error->all(FLERR,"Cannot use fix nvt/npt/nph with "
                "both xz dynamics and xz scaling");
   if (p_flag[5] && scalexy == 1)
     error->all(FLERR,"Cannot use fix nvt/npt/nph with "
                "both xy dynamics and xy scaling");
 
   if (!domain->triclinic && (p_flag[3] || p_flag[4] || p_flag[5]))
     error->all(FLERR,"Can not specify Pxy/Pxz/Pyz in "
                "fix nvt/npt/nph with non-triclinic box");
 
   if (pcouple == XYZ && dimension == 3 &&
       (p_start[0] != p_start[1] || p_start[0] != p_start[2] ||
        p_stop[0] != p_stop[1] || p_stop[0] != p_stop[2] ||
        p_period[0] != p_period[1] || p_period[0] != p_period[2]))
     error->all(FLERR,"Invalid fix nvt/npt/nph pressure settings");
   if (pcouple == XYZ && dimension == 2 &&
       (p_start[0] != p_start[1] || p_stop[0] != p_stop[1] ||
        p_period[0] != p_period[1]))
     error->all(FLERR,"Invalid fix nvt/npt/nph pressure settings");
   if (pcouple == XY &&
       (p_start[0] != p_start[1] || p_stop[0] != p_stop[1] ||
        p_period[0] != p_period[1]))
     error->all(FLERR,"Invalid fix nvt/npt/nph pressure settings");
   if (pcouple == YZ &&
       (p_start[1] != p_start[2] || p_stop[1] != p_stop[2] ||
        p_period[1] != p_period[2]))
     error->all(FLERR,"Invalid fix nvt/npt/nph pressure settings");
   if (pcouple == XZ &&
       (p_start[0] != p_start[2] || p_stop[0] != p_stop[2] ||
        p_period[0] != p_period[2]))
     error->all(FLERR,"Invalid fix nvt/npt/nph pressure settings");
 
   if (dipole_flag) {
     if (!atom->sphere_flag)
       error->all(FLERR,"Using update dipole flag requires atom style sphere");
     if (!atom->mu_flag)
       error->all(FLERR,"Using update dipole flag requires atom attribute mu");
   }
   
   if ((tstat_flag && t_period <= 0.0) ||
       (p_flag[0] && p_period[0] <= 0.0) ||
       (p_flag[1] && p_period[1] <= 0.0) ||
       (p_flag[2] && p_period[2] <= 0.0) ||
       (p_flag[3] && p_period[3] <= 0.0) ||
       (p_flag[4] && p_period[4] <= 0.0) ||
       (p_flag[5] && p_period[5] <= 0.0))
     error->all(FLERR,"Fix nvt/npt/nph damping parameters must be > 0.0");
 
   // set pstat_flag and box change and restart_pbc variables
 
   pre_exchange_flag = 0;
   pstat_flag = 0;
   pstyle = ISO;
 
   for (int i = 0; i < 6; i++)
     if (p_flag[i]) pstat_flag = 1;
 
   if (pstat_flag) {
     if (p_flag[0] || p_flag[1] || p_flag[2]) box_change_size = 1;
     if (p_flag[3] || p_flag[4] || p_flag[5]) box_change_shape = 1;
     no_change_box = 1;
     if (allremap == 0) restart_pbc = 1;
 
     // pstyle = TRICLINIC if any off-diagonal term is controlled -> 6 dof
     // else pstyle = ISO if XYZ coupling or XY coupling in 2d -> 1 dof
     // else pstyle = ANISO -> 3 dof
 
     if (p_flag[3] || p_flag[4] || p_flag[5]) pstyle = TRICLINIC;
     else if (pcouple == XYZ || (dimension == 2 && pcouple == XY)) pstyle = ISO;
     else pstyle = ANISO;
 
     // pre_exchange only required if flips can occur due to shape changes
 
     if (flipflag && (p_flag[3] || p_flag[4] || p_flag[5])) pre_exchange_flag = 1;
     if (flipflag && (domain->yz != 0.0 || domain->xz != 0.0 || domain->xy != 0.0))
       pre_exchange_flag = 1;
   }
 
   // convert input periods to frequencies
 
   t_freq = 0.0;
   p_freq[0] = p_freq[1] = p_freq[2] = p_freq[3] = p_freq[4] = p_freq[5] = 0.0;
 
   if (tstat_flag) t_freq = 1.0 / t_period;
   if (p_flag[0]) p_freq[0] = 1.0 / p_period[0];
   if (p_flag[1]) p_freq[1] = 1.0 / p_period[1];
   if (p_flag[2]) p_freq[2] = 1.0 / p_period[2];
   if (p_flag[3]) p_freq[3] = 1.0 / p_period[3];
   if (p_flag[4]) p_freq[4] = 1.0 / p_period[4];
   if (p_flag[5]) p_freq[5] = 1.0 / p_period[5];
 
   // Nose/Hoover temp and pressure init
 
   size_vector = 0;
 
   if (tstat_flag) {
     int ich;
     eta = new double[mtchain];
 
     // add one extra dummy thermostat, set to zero
 
     eta_dot = new double[mtchain+1];
     eta_dot[mtchain] = 0.0;
     eta_dotdot = new double[mtchain];
     for (ich = 0; ich < mtchain; ich++) {
       eta[ich] = eta_dot[ich] = eta_dotdot[ich] = 0.0;
     }
     eta_mass = new double[mtchain];
     size_vector += 2*2*mtchain;
   }
 
   if (pstat_flag) {
     omega[0] = omega[1] = omega[2] = 0.0;
     omega_dot[0] = omega_dot[1] = omega_dot[2] = 0.0;
     omega_mass[0] = omega_mass[1] = omega_mass[2] = 0.0;
     omega[3] = omega[4] = omega[5] = 0.0;
     omega_dot[3] = omega_dot[4] = omega_dot[5] = 0.0;
     omega_mass[3] = omega_mass[4] = omega_mass[5] = 0.0;
     if (pstyle == ISO) size_vector += 2*2*1;
     else if (pstyle == ANISO) size_vector += 2*2*3;
     else if (pstyle == TRICLINIC) size_vector += 2*2*6;
 
     if (mpchain) {
       int ich;
       etap = new double[mpchain];
 
       // add one extra dummy thermostat, set to zero
 
       etap_dot = new double[mpchain+1];
       etap_dot[mpchain] = 0.0;
       etap_dotdot = new double[mpchain];
       for (ich = 0; ich < mpchain; ich++) {
         etap[ich] = etap_dot[ich] =
           etap_dotdot[ich] = 0.0;
       }
       etap_mass = new double[mpchain];
       size_vector += 2*2*mpchain;
     }
 
     if (deviatoric_flag) size_vector += 1;
   }
 
   nrigid = 0;
   rfix = NULL;
 
   if (pre_exchange_flag) irregular = new Irregular(lmp);
   else irregular = NULL;
 
   // initialize vol0,t0 to zero to signal uninitialized
   // values then assigned in init(), if necessary
 
   vol0 = t0 = 0.0;
 }
 
 /* ---------------------------------------------------------------------- */
 
 FixNH::~FixNH()
 {
   if (copymode) return;
 
   delete [] id_dilate;
   delete [] rfix;
 
   delete irregular;
 
   // delete temperature and pressure if fix created them
 
   if (tcomputeflag) modify->delete_compute(id_temp);
   delete [] id_temp;
 
   if (tstat_flag) {
     delete [] eta;
     delete [] eta_dot;
     delete [] eta_dotdot;
     delete [] eta_mass;
   }
 
   if (pstat_flag) {
     if (pcomputeflag) modify->delete_compute(id_press);
     delete [] id_press;
     if (mpchain) {
       delete [] etap;
       delete [] etap_dot;
       delete [] etap_dotdot;
       delete [] etap_mass;
     }
   }
 }
 
 /* ---------------------------------------------------------------------- */
 
 int FixNH::setmask()
 {
   int mask = 0;
   mask |= INITIAL_INTEGRATE;
   mask |= FINAL_INTEGRATE;
   mask |= THERMO_ENERGY;
   mask |= INITIAL_INTEGRATE_RESPA;
   mask |= FINAL_INTEGRATE_RESPA;
   if (pre_exchange_flag) mask |= PRE_EXCHANGE;
   return mask;
 }
 
 /* ---------------------------------------------------------------------- */
 
 void FixNH::init()
 {
   // recheck that dilate group has not been deleted
 
   if (allremap == 0) {
     int idilate = group->find(id_dilate);
     if (idilate == -1)
       error->all(FLERR,"Fix nvt/npt/nph dilate group ID does not exist");
     dilate_group_bit = group->bitmask[idilate];
   }
 
   // ensure no conflict with fix deform
 
   if (pstat_flag)
     for (int i = 0; i < modify->nfix; i++)
       if (strcmp(modify->fix[i]->style,"deform") == 0) {
         int *dimflag = ((FixDeform *) modify->fix[i])->dimflag;
         if ((p_flag[0] && dimflag[0]) || (p_flag[1] && dimflag[1]) ||
             (p_flag[2] && dimflag[2]) || (p_flag[3] && dimflag[3]) ||
             (p_flag[4] && dimflag[4]) || (p_flag[5] && dimflag[5]))
           error->all(FLERR,"Cannot use fix npt and fix deform on "
                      "same component of stress tensor");
       }
 
   // set temperature and pressure ptrs
 
   int icompute = modify->find_compute(id_temp);
   if (icompute < 0)
     error->all(FLERR,"Temperature ID for fix nvt/npt does not exist");
   temperature = modify->compute[icompute];
 
   if (temperature->tempbias) which = BIAS;
   else which = NOBIAS;
 
   if (pstat_flag) {
     icompute = modify->find_compute(id_press);
     if (icompute < 0)
       error->all(FLERR,"Pressure ID for fix npt/nph does not exist");
     pressure = modify->compute[icompute];
   }
 
   // set timesteps and frequencies
 
   dtv = update->dt;
   dtf = 0.5 * update->dt * force->ftm2v;
   dthalf = 0.5 * update->dt;
   dt4 = 0.25 * update->dt;
   dt8 = 0.125 * update->dt;
   dto = dthalf;
 
   p_freq_max = 0.0;
   if (pstat_flag) {
     p_freq_max = MAX(p_freq[0],p_freq[1]);
     p_freq_max = MAX(p_freq_max,p_freq[2]);
     if (pstyle == TRICLINIC) {
       p_freq_max = MAX(p_freq_max,p_freq[3]);
       p_freq_max = MAX(p_freq_max,p_freq[4]);
       p_freq_max = MAX(p_freq_max,p_freq[5]);
     }
     pdrag_factor = 1.0 - (update->dt * p_freq_max * drag / nc_pchain);
   }
 
   if (tstat_flag)
     tdrag_factor = 1.0 - (update->dt * t_freq * drag / nc_tchain);
 
   // tally the number of dimensions that are barostatted
   // set initial volume and reference cell, if not already done
 
   if (pstat_flag) {
     pdim = p_flag[0] + p_flag[1] + p_flag[2];
     if (vol0 == 0.0) {
       if (dimension == 3) vol0 = domain->xprd * domain->yprd * domain->zprd;
       else vol0 = domain->xprd * domain->yprd;
       h0_inv[0] = domain->h_inv[0];
       h0_inv[1] = domain->h_inv[1];
       h0_inv[2] = domain->h_inv[2];
       h0_inv[3] = domain->h_inv[3];
       h0_inv[4] = domain->h_inv[4];
       h0_inv[5] = domain->h_inv[5];
     }
   }
 
   boltz = force->boltz;
   nktv2p = force->nktv2p;
 
   if (force->kspace) kspace_flag = 1;
   else kspace_flag = 0;
 
   if (strstr(update->integrate_style,"respa")) {
     nlevels_respa = ((Respa *) update->integrate)->nlevels;
     step_respa = ((Respa *) update->integrate)->step;
     dto = 0.5*step_respa[0];
   }
 
   // detect if any rigid fixes exist so rigid bodies move when box is remapped
   // rfix[] = indices to each fix rigid
 
   delete [] rfix;
   nrigid = 0;
   rfix = NULL;
 
   for (int i = 0; i < modify->nfix; i++)
     if (modify->fix[i]->rigid_flag) nrigid++;
   if (nrigid) {
     rfix = new int[nrigid];
     nrigid = 0;
     for (int i = 0; i < modify->nfix; i++)
       if (modify->fix[i]->rigid_flag) rfix[nrigid++] = i;
   }
 }
 
 /* ----------------------------------------------------------------------
    compute T,P before integrator starts
 ------------------------------------------------------------------------- */
 
 void FixNH::setup(int vflag)
 {
   // tdof needed by compute_temp_target()
 
   t_current = temperature->compute_scalar();
   tdof = temperature->dof;
   
   // t_target is needed by NVT and NPT in compute_scalar()
   // If no thermostat or using fix nphug,
   // t_target must be defined by other means.
 
   if (tstat_flag && strstr(style,"nphug") == NULL) {
     compute_temp_target();
   } else if (pstat_flag) {
 
     // t0 = reference temperature for masses
     // cannot be done in init() b/c temperature cannot be called there
     // is b/c Modify::init() inits computes after fixes due to dof dependence
     // guesstimate a unit-dependent t0 if actual T = 0.0
     // if it was read in from a restart file, leave it be
 
     if (t0 == 0.0) {
       t0 = temperature->compute_scalar();
       if (t0 == 0.0) {
         if (strcmp(update->unit_style,"lj") == 0) t0 = 1.0;
         else t0 = 300.0;
       }
     }
     t_target = t0;
   }
 
   if (pstat_flag) compute_press_target();
 
   if (pstat_flag) {
     if (pstyle == ISO) pressure->compute_scalar();
     else pressure->compute_vector();
     couple();
     pressure->addstep(update->ntimestep+1);
   }
 
   // masses and initial forces on thermostat variables
 
   if (tstat_flag) {
     eta_mass[0] = tdof * boltz * t_target / (t_freq*t_freq);
     for (int ich = 1; ich < mtchain; ich++)
       eta_mass[ich] = boltz * t_target / (t_freq*t_freq);
     for (int ich = 1; ich < mtchain; ich++) {
       eta_dotdot[ich] = (eta_mass[ich-1]*eta_dot[ich-1]*eta_dot[ich-1] -
                          boltz * t_target) / eta_mass[ich];
     }
   }
 
   // masses and initial forces on barostat variables
 
   if (pstat_flag) {
     double kt = boltz * t_target;
     double nkt = atom->natoms * kt;
 
     for (int i = 0; i < 3; i++)
       if (p_flag[i])
         omega_mass[i] = nkt/(p_freq[i]*p_freq[i]);
 
     if (pstyle == TRICLINIC) {
       for (int i = 3; i < 6; i++)
         if (p_flag[i]) omega_mass[i] = nkt/(p_freq[i]*p_freq[i]);
     }
 
   // masses and initial forces on barostat thermostat variables
 
     if (mpchain) {
       etap_mass[0] = boltz * t_target / (p_freq_max*p_freq_max);
       for (int ich = 1; ich < mpchain; ich++)
         etap_mass[ich] = boltz * t_target / (p_freq_max*p_freq_max);
       for (int ich = 1; ich < mpchain; ich++)
         etap_dotdot[ich] =
           (etap_mass[ich-1]*etap_dot[ich-1]*etap_dot[ich-1] -
            boltz * t_target) / etap_mass[ich];
     }
 
   }
 }
 
 /* ----------------------------------------------------------------------
    1st half of Verlet update
 ------------------------------------------------------------------------- */
 
 void FixNH::initial_integrate(int vflag)
 {
   // update eta_press_dot
 
   if (pstat_flag && mpchain) nhc_press_integrate();
 
   // update eta_dot
 
   if (tstat_flag) {
     compute_temp_target();
     nhc_temp_integrate();
   }
 
   // need to recompute pressure to account for change in KE
   // t_current is up-to-date, but compute_temperature is not
   // compute appropriately coupled elements of mvv_current
 
   if (pstat_flag) {
     if (pstyle == ISO) {
       temperature->compute_scalar();
       pressure->compute_scalar();
     } else {
       temperature->compute_vector();
       pressure->compute_vector();
     }
     couple();
     pressure->addstep(update->ntimestep+1);
   }
 
   if (pstat_flag) {
     compute_press_target();
     nh_omega_dot();
     nh_v_press();
   }
 
   nve_v();
 
   // remap simulation box by 1/2 step
 
   if (pstat_flag) remap();
 
   nve_x();
 
   // remap simulation box by 1/2 step
   // redo KSpace coeffs since volume has changed
 
   if (pstat_flag) {
     remap();
     if (kspace_flag) force->kspace->setup();
   }
 }
 
 /* ----------------------------------------------------------------------
    2nd half of Verlet update
 ------------------------------------------------------------------------- */
 
 void FixNH::final_integrate()
 {
   nve_v();
 
   // re-compute temp before nh_v_press()
   // only needed for temperature computes with BIAS on reneighboring steps:
   //   b/c some biases store per-atom values (e.g. temp/profile)
   //   per-atom values are invalid if reneigh/comm occurred
   //     since temp->compute() in initial_integrate()
 
   if (which == BIAS && neighbor->ago == 0)
     t_current = temperature->compute_scalar();
 
   if (pstat_flag) nh_v_press();
 
   // compute new T,P after velocities rescaled by nh_v_press()
   // compute appropriately coupled elements of mvv_current
 
   t_current = temperature->compute_scalar();
   tdof = temperature->dof;
 
   if (pstat_flag) {
     if (pstyle == ISO) pressure->compute_scalar();
     else pressure->compute_vector();
     couple();
     pressure->addstep(update->ntimestep+1);
   }
 
   if (pstat_flag) nh_omega_dot();
 
   // update eta_dot
   // update eta_press_dot
 
   if (tstat_flag) nhc_temp_integrate();
   if (pstat_flag && mpchain) nhc_press_integrate();
 }
 
 /* ---------------------------------------------------------------------- */
 
 void FixNH::initial_integrate_respa(int vflag, int ilevel, int iloop)
 {
   // set timesteps by level
 
   dtv = step_respa[ilevel];
   dtf = 0.5 * step_respa[ilevel] * force->ftm2v;
   dthalf = 0.5 * step_respa[ilevel];
 
   // outermost level - update eta_dot and omega_dot, apply to v
   // all other levels - NVE update of v
   // x,v updates only performed for atoms in group
 
   if (ilevel == nlevels_respa-1) {
 
     // update eta_press_dot
 
     if (pstat_flag && mpchain) nhc_press_integrate();
 
     // update eta_dot
 
     if (tstat_flag) {
       compute_temp_target();
       nhc_temp_integrate();
     }
 
     // recompute pressure to account for change in KE
     // t_current is up-to-date, but compute_temperature is not
     // compute appropriately coupled elements of mvv_current
 
     if (pstat_flag) {
       if (pstyle == ISO) {
         temperature->compute_scalar();
         pressure->compute_scalar();
       } else {
         temperature->compute_vector();
         pressure->compute_vector();
       }
       couple();
       pressure->addstep(update->ntimestep+1);
     }
 
     if (pstat_flag) {
       compute_press_target();
       nh_omega_dot();
       nh_v_press();
     }
 
     nve_v();
 
   } else nve_v();
 
   // innermost level - also update x only for atoms in group
   // if barostat, perform 1/2 step remap before and after
 
   if (ilevel == 0) {
     if (pstat_flag) remap();
     nve_x();
     if (pstat_flag) remap();
   }
 
   // if barostat, redo KSpace coeffs at outermost level,
   // since volume has changed
 
   if (ilevel == nlevels_respa-1 && kspace_flag && pstat_flag)
     force->kspace->setup();
 }
 
 /* ---------------------------------------------------------------------- */
 
 void FixNH::final_integrate_respa(int ilevel, int iloop)
 {
   // set timesteps by level
 
   dtf = 0.5 * step_respa[ilevel] * force->ftm2v;
   dthalf = 0.5 * step_respa[ilevel];
 
   // outermost level - update eta_dot and omega_dot, apply via final_integrate
   // all other levels - NVE update of v
 
   if (ilevel == nlevels_respa-1) final_integrate();
   else nve_v();
 }
 
 /* ---------------------------------------------------------------------- */
 
 void FixNH::couple()
 {
   double *tensor = pressure->vector;
 
   if (pstyle == ISO)
     p_current[0] = p_current[1] = p_current[2] = pressure->scalar;
   else if (pcouple == XYZ) {
     double ave = 1.0/3.0 * (tensor[0] + tensor[1] + tensor[2]);
     p_current[0] = p_current[1] = p_current[2] = ave;
   } else if (pcouple == XY) {
     double ave = 0.5 * (tensor[0] + tensor[1]);
     p_current[0] = p_current[1] = ave;
     p_current[2] = tensor[2];
   } else if (pcouple == YZ) {
     double ave = 0.5 * (tensor[1] + tensor[2]);
     p_current[1] = p_current[2] = ave;
     p_current[0] = tensor[0];
   } else if (pcouple == XZ) {
     double ave = 0.5 * (tensor[0] + tensor[2]);
     p_current[0] = p_current[2] = ave;
     p_current[1] = tensor[1];
   } else {
     p_current[0] = tensor[0];
     p_current[1] = tensor[1];
     p_current[2] = tensor[2];
   }
 
   if (!ISFINITE(p_current[0]) || !ISFINITE(p_current[1]) || !ISFINITE(p_current[2]))
     error->all(FLERR,"Non-numeric pressure - simulation unstable");
 
   // switch order from xy-xz-yz to Voigt
 
   if (pstyle == TRICLINIC) {
     p_current[3] = tensor[5];
     p_current[4] = tensor[4];
     p_current[5] = tensor[3];
 
     if (!ISFINITE(p_current[3]) || !ISFINITE(p_current[4]) || !ISFINITE(p_current[5]))
       error->all(FLERR,"Non-numeric pressure - simulation unstable");
   }
 }
 
 /* ----------------------------------------------------------------------
    change box size
    remap all atoms or dilate group atoms depending on allremap flag
    if rigid bodies exist, scale rigid body centers-of-mass
 ------------------------------------------------------------------------- */
 
 void FixNH::remap()
 {
   int i;
   double oldlo,oldhi;
   double expfac;
 
   double **x = atom->x;
   int *mask = atom->mask;
   int nlocal = atom->nlocal;
   double *h = domain->h;
 
   // omega is not used, except for book-keeping
 
   for (int i = 0; i < 6; i++) omega[i] += dto*omega_dot[i];
 
   // convert pertinent atoms and rigid bodies to lamda coords
 
   if (allremap) domain->x2lamda(nlocal);
   else {
     for (i = 0; i < nlocal; i++)
       if (mask[i] & dilate_group_bit)
         domain->x2lamda(x[i],x[i]);
   }
 
   if (nrigid)
     for (i = 0; i < nrigid; i++)
       modify->fix[rfix[i]]->deform(0);
 
   // reset global and local box to new size/shape
 
   // this operation corresponds to applying the
   // translate and scale operations
   // corresponding to the solution of the following ODE:
   //
   // h_dot = omega_dot * h
   //
   // where h_dot, omega_dot and h are all upper-triangular
   // 3x3 tensors. In Voigt notation, the elements of the
   // RHS product tensor are:
   // h_dot = [0*0, 1*1, 2*2, 1*3+3*2, 0*4+5*3+4*2, 0*5+5*1]
   //
   // Ordering of operations preserves time symmetry.
 
   double dto2 = dto/2.0;
   double dto4 = dto/4.0;
   double dto8 = dto/8.0;
 
   // off-diagonal components, first half
 
   if (pstyle == TRICLINIC) {
 
     if (p_flag[4]) {
       expfac = exp(dto8*omega_dot[0]);
       h[4] *= expfac;
       h[4] += dto4*(omega_dot[5]*h[3]+omega_dot[4]*h[2]);
       h[4] *= expfac;
     }
 
     if (p_flag[3]) {
       expfac = exp(dto4*omega_dot[1]);
       h[3] *= expfac;
       h[3] += dto2*(omega_dot[3]*h[2]);
       h[3] *= expfac;
     }
 
     if (p_flag[5]) {
       expfac = exp(dto4*omega_dot[0]);
       h[5] *= expfac;
       h[5] += dto2*(omega_dot[5]*h[1]);
       h[5] *= expfac;
     }
 
     if (p_flag[4]) {
       expfac = exp(dto8*omega_dot[0]);
       h[4] *= expfac;
       h[4] += dto4*(omega_dot[5]*h[3]+omega_dot[4]*h[2]);
       h[4] *= expfac;
     }
   }
 
   // scale diagonal components
   // scale tilt factors with cell, if set
 
   if (p_flag[0]) {
     oldlo = domain->boxlo[0];
     oldhi = domain->boxhi[0];
     expfac = exp(dto*omega_dot[0]);
     domain->boxlo[0] = (oldlo-fixedpoint[0])*expfac + fixedpoint[0];
     domain->boxhi[0] = (oldhi-fixedpoint[0])*expfac + fixedpoint[0];
   }
 
   if (p_flag[1]) {
     oldlo = domain->boxlo[1];
     oldhi = domain->boxhi[1];
     expfac = exp(dto*omega_dot[1]);
     domain->boxlo[1] = (oldlo-fixedpoint[1])*expfac + fixedpoint[1];
     domain->boxhi[1] = (oldhi-fixedpoint[1])*expfac + fixedpoint[1];
     if (scalexy) h[5] *= expfac;
   }
 
   if (p_flag[2]) {
     oldlo = domain->boxlo[2];
     oldhi = domain->boxhi[2];
     expfac = exp(dto*omega_dot[2]);
     domain->boxlo[2] = (oldlo-fixedpoint[2])*expfac + fixedpoint[2];
     domain->boxhi[2] = (oldhi-fixedpoint[2])*expfac + fixedpoint[2];
     if (scalexz) h[4] *= expfac;
     if (scaleyz) h[3] *= expfac;
   }
 
   // off-diagonal components, second half
 
   if (pstyle == TRICLINIC) {
 
     if (p_flag[4]) {
       expfac = exp(dto8*omega_dot[0]);
       h[4] *= expfac;
       h[4] += dto4*(omega_dot[5]*h[3]+omega_dot[4]*h[2]);
       h[4] *= expfac;
     }
 
     if (p_flag[3]) {
       expfac = exp(dto4*omega_dot[1]);
       h[3] *= expfac;
       h[3] += dto2*(omega_dot[3]*h[2]);
       h[3] *= expfac;
     }
 
     if (p_flag[5]) {
       expfac = exp(dto4*omega_dot[0]);
       h[5] *= expfac;
       h[5] += dto2*(omega_dot[5]*h[1]);
       h[5] *= expfac;
     }
 
     if (p_flag[4]) {
       expfac = exp(dto8*omega_dot[0]);
       h[4] *= expfac;
       h[4] += dto4*(omega_dot[5]*h[3]+omega_dot[4]*h[2]);
       h[4] *= expfac;
     }
 
   }
 
   domain->yz = h[3];
   domain->xz = h[4];
   domain->xy = h[5];
 
   // tilt factor to cell length ratio can not exceed TILTMAX in one step
 
   if (domain->yz < -TILTMAX*domain->yprd ||
       domain->yz > TILTMAX*domain->yprd ||
       domain->xz < -TILTMAX*domain->xprd ||
       domain->xz > TILTMAX*domain->xprd ||
       domain->xy < -TILTMAX*domain->xprd ||
       domain->xy > TILTMAX*domain->xprd)
     error->all(FLERR,"Fix npt/nph has tilted box too far in one step - "
                "periodic cell is too far from equilibrium state");
 
   domain->set_global_box();
   domain->set_local_box();
 
   // convert pertinent atoms and rigid bodies back to box coords
 
   if (allremap) domain->lamda2x(nlocal);
   else {
     for (i = 0; i < nlocal; i++)
       if (mask[i] & dilate_group_bit)
         domain->lamda2x(x[i],x[i]);
   }
 
   if (nrigid)
     for (i = 0; i < nrigid; i++)
       modify->fix[rfix[i]]->deform(1);
 }
 
 /* ----------------------------------------------------------------------
    pack entire state of Fix into one write
 ------------------------------------------------------------------------- */
 
 void FixNH::write_restart(FILE *fp)
 {
   int nsize = size_restart_global();
 
   double *list;
   memory->create(list,nsize,"nh:list");
 
   pack_restart_data(list);
 
   if (comm->me == 0) {
     int size = nsize * sizeof(double);
     fwrite(&size,sizeof(int),1,fp);
     fwrite(list,sizeof(double),nsize,fp);
   }
 
   memory->destroy(list);
 }
 
 /* ----------------------------------------------------------------------
     calculate the number of data to be packed
 ------------------------------------------------------------------------- */
 
 int FixNH::size_restart_global()
 {
   int nsize = 2;
   if (tstat_flag) nsize += 1 + 2*mtchain;
   if (pstat_flag) {
     nsize += 16 + 2*mpchain;
     if (deviatoric_flag) nsize += 6;
   }
 
   return nsize;
 }
 
 /* ----------------------------------------------------------------------
    pack restart data
 ------------------------------------------------------------------------- */
 
 int FixNH::pack_restart_data(double *list)
 {
   int n = 0;
 
   list[n++] = tstat_flag;
   if (tstat_flag) {
     list[n++] = mtchain;
     for (int ich = 0; ich < mtchain; ich++)
       list[n++] = eta[ich];
     for (int ich = 0; ich < mtchain; ich++)
       list[n++] = eta_dot[ich];
   }
 
   list[n++] = pstat_flag;
   if (pstat_flag) {
     list[n++] = omega[0];
     list[n++] = omega[1];
     list[n++] = omega[2];
     list[n++] = omega[3];
     list[n++] = omega[4];
     list[n++] = omega[5];
     list[n++] = omega_dot[0];
     list[n++] = omega_dot[1];
     list[n++] = omega_dot[2];
     list[n++] = omega_dot[3];
     list[n++] = omega_dot[4];
     list[n++] = omega_dot[5];
     list[n++] = vol0;
     list[n++] = t0;
     list[n++] = mpchain;
     if (mpchain) {
       for (int ich = 0; ich < mpchain; ich++)
         list[n++] = etap[ich];
       for (int ich = 0; ich < mpchain; ich++)
         list[n++] = etap_dot[ich];
     }
 
     list[n++] = deviatoric_flag;
     if (deviatoric_flag) {
       list[n++] = h0_inv[0];
       list[n++] = h0_inv[1];
       list[n++] = h0_inv[2];
       list[n++] = h0_inv[3];
       list[n++] = h0_inv[4];
       list[n++] = h0_inv[5];
     }
   }
 
   return n;
 }
 
 /* ----------------------------------------------------------------------
    use state info from restart file to restart the Fix
 ------------------------------------------------------------------------- */
 
 void FixNH::restart(char *buf)
 {
   int n = 0;
   double *list = (double *) buf;
   int flag = static_cast<int> (list[n++]);
   if (flag) {
     int m = static_cast<int> (list[n++]);
     if (tstat_flag && m == mtchain) {
       for (int ich = 0; ich < mtchain; ich++)
         eta[ich] = list[n++];
       for (int ich = 0; ich < mtchain; ich++)
         eta_dot[ich] = list[n++];
     } else n += 2*m;
   }
   flag = static_cast<int> (list[n++]);
   if (flag) {
     omega[0] = list[n++];
     omega[1] = list[n++];
     omega[2] = list[n++];
     omega[3] = list[n++];
     omega[4] = list[n++];
     omega[5] = list[n++];
     omega_dot[0] = list[n++];
     omega_dot[1] = list[n++];
     omega_dot[2] = list[n++];
     omega_dot[3] = list[n++];
     omega_dot[4] = list[n++];
     omega_dot[5] = list[n++];
     vol0 = list[n++];
     t0 = list[n++];
     int m = static_cast<int> (list[n++]);
     if (pstat_flag && m == mpchain) {
       for (int ich = 0; ich < mpchain; ich++)
         etap[ich] = list[n++];
       for (int ich = 0; ich < mpchain; ich++)
         etap_dot[ich] = list[n++];
     } else n+=2*m;
     flag = static_cast<int> (list[n++]);
     if (flag) {
       h0_inv[0] = list[n++];
       h0_inv[1] = list[n++];
       h0_inv[2] = list[n++];
       h0_inv[3] = list[n++];
       h0_inv[4] = list[n++];
       h0_inv[5] = list[n++];
     }
   }
 }
 
 /* ---------------------------------------------------------------------- */
 
 int FixNH::modify_param(int narg, char **arg)
 {
   if (strcmp(arg[0],"temp") == 0) {
     if (narg < 2) error->all(FLERR,"Illegal fix_modify command");
     if (tcomputeflag) {
       modify->delete_compute(id_temp);
       tcomputeflag = 0;
     }
     delete [] id_temp;
     int n = strlen(arg[1]) + 1;
     id_temp = new char[n];
     strcpy(id_temp,arg[1]);
 
     int icompute = modify->find_compute(arg[1]);
     if (icompute < 0)
       error->all(FLERR,"Could not find fix_modify temperature ID");
     temperature = modify->compute[icompute];
 
     if (temperature->tempflag == 0)
       error->all(FLERR,
                  "Fix_modify temperature ID does not compute temperature");
     if (temperature->igroup != 0 && comm->me == 0)
       error->warning(FLERR,"Temperature for fix modify is not for group all");
 
     // reset id_temp of pressure to new temperature ID
 
     if (pstat_flag) {
       icompute = modify->find_compute(id_press);
       if (icompute < 0)
         error->all(FLERR,"Pressure ID for fix modify does not exist");
       modify->compute[icompute]->reset_extra_compute_fix(id_temp);
     }
 
     return 2;
 
   } else if (strcmp(arg[0],"press") == 0) {
     if (narg < 2) error->all(FLERR,"Illegal fix_modify command");
     if (!pstat_flag) error->all(FLERR,"Illegal fix_modify command");
     if (pcomputeflag) {
       modify->delete_compute(id_press);
       pcomputeflag = 0;
     }
     delete [] id_press;
     int n = strlen(arg[1]) + 1;
     id_press = new char[n];
     strcpy(id_press,arg[1]);
 
     int icompute = modify->find_compute(arg[1]);
     if (icompute < 0) error->all(FLERR,"Could not find fix_modify pressure ID");
     pressure = modify->compute[icompute];
 
     if (pressure->pressflag == 0)
       error->all(FLERR,"Fix_modify pressure ID does not compute pressure");
     return 2;
   }
 
   return 0;
 }
 
 /* ---------------------------------------------------------------------- */
 
 double FixNH::compute_scalar()
 {
   int i;
   double volume;
   double energy;
   double kt = boltz * t_target;
   double lkt_press = kt;
   int ich;
   if (dimension == 3) volume = domain->xprd * domain->yprd * domain->zprd;
   else volume = domain->xprd * domain->yprd;
 
   energy = 0.0;
 
   // thermostat chain energy is equivalent to Eq. (2) in
   // Martyna, Tuckerman, Tobias, Klein, Mol Phys, 87, 1117
   // Sum(0.5*p_eta_k^2/Q_k,k=1,M) + L*k*T*eta_1 + Sum(k*T*eta_k,k=2,M),
   // where L = tdof
   //       M = mtchain
   //       p_eta_k = Q_k*eta_dot[k-1]
   //       Q_1 = L*k*T/t_freq^2
   //       Q_k = k*T/t_freq^2, k > 1
 
   if (tstat_flag) {
     energy += ke_target * eta[0] + 0.5*eta_mass[0]*eta_dot[0]*eta_dot[0];
     for (ich = 1; ich < mtchain; ich++)
       energy += kt * eta[ich] + 0.5*eta_mass[ich]*eta_dot[ich]*eta_dot[ich];
   }
 
   // barostat energy is equivalent to Eq. (8) in
   // Martyna, Tuckerman, Tobias, Klein, Mol Phys, 87, 1117
   // Sum(0.5*p_omega^2/W + P*V),
   // where N = natoms
   //       p_omega = W*omega_dot
   //       W = N*k*T/p_freq^2
   //       sum is over barostatted dimensions
 
   if (pstat_flag) {
     for (i = 0; i < 3; i++)
       if (p_flag[i])
         energy += 0.5*omega_dot[i]*omega_dot[i]*omega_mass[i] +
           p_hydro*(volume-vol0) / (pdim*nktv2p);
 
     if (pstyle == TRICLINIC) {
       for (i = 3; i < 6; i++)
         if (p_flag[i])
           energy += 0.5*omega_dot[i]*omega_dot[i]*omega_mass[i];
     }
 
     // extra contributions from thermostat chain for barostat
 
     if (mpchain) {
       energy += lkt_press * etap[0] + 0.5*etap_mass[0]*etap_dot[0]*etap_dot[0];
       for (ich = 1; ich < mpchain; ich++)
         energy += kt * etap[ich] +
           0.5*etap_mass[ich]*etap_dot[ich]*etap_dot[ich];
     }
 
     // extra contribution from strain energy
 
     if (deviatoric_flag) energy += compute_strain_energy();
   }
 
   return energy;
 }
 
 /* ----------------------------------------------------------------------
    return a single element of the following vectors, in this order:
       eta[tchain], eta_dot[tchain], omega[ndof], omega_dot[ndof]
       etap[pchain], etap_dot[pchain], PE_eta[tchain], KE_eta_dot[tchain]
       PE_omega[ndof], KE_omega_dot[ndof], PE_etap[pchain], KE_etap_dot[pchain]
       PE_strain[1]
   if no thermostat exists, related quantities are omitted from the list
   if no barostat exists, related quantities are omitted from the list
   ndof = 1,3,6 degrees of freedom for pstyle = ISO,ANISO,TRI
 ------------------------------------------------------------------------- */
 
 double FixNH::compute_vector(int n)
 {
   int ilen;
 
   if (tstat_flag) {
     ilen = mtchain;
     if (n < ilen) return eta[n];
     n -= ilen;
     ilen = mtchain;
     if (n < ilen) return eta_dot[n];
     n -= ilen;
   }
 
   if (pstat_flag) {
     if (pstyle == ISO) {
       ilen = 1;
       if (n < ilen) return omega[n];
       n -= ilen;
     } else if (pstyle == ANISO) {
       ilen = 3;
       if (n < ilen) return omega[n];
       n -= ilen;
     } else {
       ilen = 6;
       if (n < ilen) return omega[n];
       n -= ilen;
     }
 
     if (pstyle == ISO) {
       ilen = 1;
       if (n < ilen) return omega_dot[n];
       n -= ilen;
     } else if (pstyle == ANISO) {
       ilen = 3;
       if (n < ilen) return omega_dot[n];
       n -= ilen;
     } else {
       ilen = 6;
       if (n < ilen) return omega_dot[n];
       n -= ilen;
     }
 
     if (mpchain) {
       ilen = mpchain;
       if (n < ilen) return etap[n];
       n -= ilen;
       ilen = mpchain;
       if (n < ilen) return etap_dot[n];
       n -= ilen;
     }
   }
 
   double volume;
   double kt = boltz * t_target;
   double lkt_press = kt;
   int ich;
   if (dimension == 3) volume = domain->xprd * domain->yprd * domain->zprd;
   else volume = domain->xprd * domain->yprd;
 
   if (tstat_flag) {
     ilen = mtchain;
     if (n < ilen) {
       ich = n;
       if (ich == 0)
         return ke_target * eta[0];
       else
         return kt * eta[ich];
     }
     n -= ilen;
     ilen = mtchain;
     if (n < ilen) {
       ich = n;
       if (ich == 0)
         return 0.5*eta_mass[0]*eta_dot[0]*eta_dot[0];
       else
         return 0.5*eta_mass[ich]*eta_dot[ich]*eta_dot[ich];
     }
     n -= ilen;
   }
 
   if (pstat_flag) {
     if (pstyle == ISO) {
       ilen = 1;
       if (n < ilen)
         return p_hydro*(volume-vol0) / nktv2p;
       n -= ilen;
     } else if (pstyle == ANISO) {
       ilen = 3;
       if (n < ilen) {
         if (p_flag[n])
           return p_hydro*(volume-vol0) / (pdim*nktv2p);
         else
           return 0.0;
       }
       n -= ilen;
     } else {
       ilen = 6;
       if (n < ilen) {
         if (n > 2) return 0.0;
         else if (p_flag[n])
           return p_hydro*(volume-vol0) / (pdim*nktv2p);
         else
           return 0.0;
       }
       n -= ilen;
     }
 
     if (pstyle == ISO) {
       ilen = 1;
       if (n < ilen)
         return pdim*0.5*omega_dot[n]*omega_dot[n]*omega_mass[n];
       n -= ilen;
     } else if (pstyle == ANISO) {
       ilen = 3;
       if (n < ilen) {
         if (p_flag[n])
           return 0.5*omega_dot[n]*omega_dot[n]*omega_mass[n];
         else return 0.0;
       }
       n -= ilen;
     } else {
       ilen = 6;
       if (n < ilen) {
         if (p_flag[n])
           return 0.5*omega_dot[n]*omega_dot[n]*omega_mass[n];
         else return 0.0;
       }
       n -= ilen;
     }
 
     if (mpchain) {
       ilen = mpchain;
       if (n < ilen) {
         ich = n;
         if (ich == 0) return lkt_press * etap[0];
         else return kt * etap[ich];
       }
       n -= ilen;
       ilen = mpchain;
       if (n < ilen) {
         ich = n;
         if (ich == 0)
           return 0.5*etap_mass[0]*etap_dot[0]*etap_dot[0];
         else
           return 0.5*etap_mass[ich]*etap_dot[ich]*etap_dot[ich];
       }
       n -= ilen;
     }
 
     if (deviatoric_flag) {
       ilen = 1;
       if (n < ilen)
         return compute_strain_energy();
       n -= ilen;
     }
   }
 
   return 0.0;
 }
 
 /* ---------------------------------------------------------------------- */
 
 void FixNH::reset_target(double t_new)
 {
   t_target = t_start = t_stop = t_new;
 }
 
 /* ---------------------------------------------------------------------- */
 
 void FixNH::reset_dt()
 {
   dtv = update->dt;
   dtf = 0.5 * update->dt * force->ftm2v;
   dthalf = 0.5 * update->dt;
   dt4 = 0.25 * update->dt;
   dt8 = 0.125 * update->dt;
   dto = dthalf;
 
   // If using respa, then remap is performed in innermost level
 
   if (strstr(update->integrate_style,"respa"))
     dto = 0.5*step_respa[0];
 
   if (pstat_flag)
     pdrag_factor = 1.0 - (update->dt * p_freq_max * drag / nc_pchain);
 
   if (tstat_flag)
     tdrag_factor = 1.0 - (update->dt * t_freq * drag / nc_tchain);
 }
 
 /* ----------------------------------------------------------------------
    extract thermostat properties
 ------------------------------------------------------------------------- */
 
 void *FixNH::extract(const char *str, int &dim)
 {
   dim=0;
   if (strcmp(str,"t_target") == 0) {
     return &t_target;
   } else if (strcmp(str,"mtchain") == 0) {
     return &mtchain;
   }
   dim=1;
   if (strcmp(str,"eta") == 0) {
     return &eta;
   }
   return NULL;
 }
 
 /* ----------------------------------------------------------------------
    perform half-step update of chain thermostat variables
 ------------------------------------------------------------------------- */
 
 void FixNH::nhc_temp_integrate()
 {
   int ich;
   double expfac;
   double kecurrent = tdof * boltz * t_current;
 
   // Update masses, to preserve initial freq, if flag set
 
   if (eta_mass_flag) {
     eta_mass[0] = tdof * boltz * t_target / (t_freq*t_freq);
     for (int ich = 1; ich < mtchain; ich++)
       eta_mass[ich] = boltz * t_target / (t_freq*t_freq);
   }
 
   if (eta_mass[0] > 0.0)
     eta_dotdot[0] = (kecurrent - ke_target)/eta_mass[0];
   else eta_dotdot[0] = 0.0;
 
   double ncfac = 1.0/nc_tchain;
   for (int iloop = 0; iloop < nc_tchain; iloop++) {
 
     for (ich = mtchain-1; ich > 0; ich--) {
       expfac = exp(-ncfac*dt8*eta_dot[ich+1]);
       eta_dot[ich] *= expfac;
       eta_dot[ich] += eta_dotdot[ich] * ncfac*dt4;
       eta_dot[ich] *= tdrag_factor;
       eta_dot[ich] *= expfac;
     }
 
     expfac = exp(-ncfac*dt8*eta_dot[1]);
     eta_dot[0] *= expfac;
     eta_dot[0] += eta_dotdot[0] * ncfac*dt4;
     eta_dot[0] *= tdrag_factor;
     eta_dot[0] *= expfac;
 
     factor_eta = exp(-ncfac*dthalf*eta_dot[0]);
     nh_v_temp();
 
     // rescale temperature due to velocity scaling
     // should not be necessary to explicitly recompute the temperature
 
     t_current *= factor_eta*factor_eta;
     kecurrent = tdof * boltz * t_current;
 
     if (eta_mass[0] > 0.0)
       eta_dotdot[0] = (kecurrent - ke_target)/eta_mass[0];
     else eta_dotdot[0] = 0.0;
 
     for (ich = 0; ich < mtchain; ich++)
       eta[ich] += ncfac*dthalf*eta_dot[ich];
 
     eta_dot[0] *= expfac;
     eta_dot[0] += eta_dotdot[0] * ncfac*dt4;
     eta_dot[0] *= expfac;
 
     for (ich = 1; ich < mtchain; ich++) {
       expfac = exp(-ncfac*dt8*eta_dot[ich+1]);
       eta_dot[ich] *= expfac;
       eta_dotdot[ich] = (eta_mass[ich-1]*eta_dot[ich-1]*eta_dot[ich-1]
                          - boltz * t_target)/eta_mass[ich];
       eta_dot[ich] += eta_dotdot[ich] * ncfac*dt4;
       eta_dot[ich] *= expfac;
     }
   }
 }
 
 /* ----------------------------------------------------------------------
    perform half-step update of chain thermostat variables for barostat
    scale barostat velocities
 ------------------------------------------------------------------------- */
 
 void FixNH::nhc_press_integrate()
 {
   int ich,i;
   double expfac,factor_etap,kecurrent;
   double kt = boltz * t_target;
   double lkt_press = kt;
 
   // Update masses, to preserve initial freq, if flag set
 
   if (omega_mass_flag) {
     double nkt = atom->natoms * kt;
     for (int i = 0; i < 3; i++)
       if (p_flag[i])
         omega_mass[i] = nkt/(p_freq[i]*p_freq[i]);
 
     if (pstyle == TRICLINIC) {
       for (int i = 3; i < 6; i++)
         if (p_flag[i]) omega_mass[i] = nkt/(p_freq[i]*p_freq[i]);
     }
   }
 
   if (etap_mass_flag) {
     if (mpchain) {
       etap_mass[0] = boltz * t_target / (p_freq_max*p_freq_max);
       for (int ich = 1; ich < mpchain; ich++)
         etap_mass[ich] = boltz * t_target / (p_freq_max*p_freq_max);
       for (int ich = 1; ich < mpchain; ich++)
         etap_dotdot[ich] =
           (etap_mass[ich-1]*etap_dot[ich-1]*etap_dot[ich-1] -
            boltz * t_target) / etap_mass[ich];
     }
   }
 
   kecurrent = 0.0;
   for (i = 0; i < 3; i++)
     if (p_flag[i]) kecurrent += omega_mass[i]*omega_dot[i]*omega_dot[i];
 
   if (pstyle == TRICLINIC) {
     for (i = 3; i < 6; i++)
       if (p_flag[i]) kecurrent += omega_mass[i]*omega_dot[i]*omega_dot[i];
   }
 
   etap_dotdot[0] = (kecurrent - lkt_press)/etap_mass[0];
 
   double ncfac = 1.0/nc_pchain;
   for (int iloop = 0; iloop < nc_pchain; iloop++) {
 
     for (ich = mpchain-1; ich > 0; ich--) {
       expfac = exp(-ncfac*dt8*etap_dot[ich+1]);
       etap_dot[ich] *= expfac;
       etap_dot[ich] += etap_dotdot[ich] * ncfac*dt4;
       etap_dot[ich] *= pdrag_factor;
       etap_dot[ich] *= expfac;
     }
 
     expfac = exp(-ncfac*dt8*etap_dot[1]);
     etap_dot[0] *= expfac;
     etap_dot[0] += etap_dotdot[0] * ncfac*dt4;
     etap_dot[0] *= pdrag_factor;
     etap_dot[0] *= expfac;
 
     for (ich = 0; ich < mpchain; ich++)
       etap[ich] += ncfac*dthalf*etap_dot[ich];
 
     factor_etap = exp(-ncfac*dthalf*etap_dot[0]);
     for (i = 0; i < 3; i++)
       if (p_flag[i]) omega_dot[i] *= factor_etap;
 
     if (pstyle == TRICLINIC) {
       for (i = 3; i < 6; i++)
         if (p_flag[i]) omega_dot[i] *= factor_etap;
     }
 
     kecurrent = 0.0;
     for (i = 0; i < 3; i++)
       if (p_flag[i]) kecurrent += omega_mass[i]*omega_dot[i]*omega_dot[i];
 
     if (pstyle == TRICLINIC) {
       for (i = 3; i < 6; i++)
         if (p_flag[i]) kecurrent += omega_mass[i]*omega_dot[i]*omega_dot[i];
     }
 
     etap_dotdot[0] = (kecurrent - lkt_press)/etap_mass[0];
 
     etap_dot[0] *= expfac;
     etap_dot[0] += etap_dotdot[0] * ncfac*dt4;
     etap_dot[0] *= expfac;
 
     for (ich = 1; ich < mpchain; ich++) {
       expfac = exp(-ncfac*dt8*etap_dot[ich+1]);
       etap_dot[ich] *= expfac;
       etap_dotdot[ich] =
         (etap_mass[ich-1]*etap_dot[ich-1]*etap_dot[ich-1] - boltz*t_target) /
         etap_mass[ich];
       etap_dot[ich] += etap_dotdot[ich] * ncfac*dt4;
       etap_dot[ich] *= expfac;
     }
   }
 }
 
 /* ----------------------------------------------------------------------
    perform half-step barostat scaling of velocities
 -----------------------------------------------------------------------*/
 
 void FixNH::nh_v_press()
 {
   double factor[3];
   double **v = atom->v;
   int *mask = atom->mask;
   int nlocal = atom->nlocal;
   if (igroup == atom->firstgroup) nlocal = atom->nfirst;
 
   factor[0] = exp(-dt4*(omega_dot[0]+mtk_term2));
   factor[1] = exp(-dt4*(omega_dot[1]+mtk_term2));
   factor[2] = exp(-dt4*(omega_dot[2]+mtk_term2));
 
   if (which == NOBIAS) {
     for (int i = 0; i < nlocal; i++) {
       if (mask[i] & groupbit) {
         v[i][0] *= factor[0];
         v[i][1] *= factor[1];
         v[i][2] *= factor[2];
         if (pstyle == TRICLINIC) {
           v[i][0] += -dthalf*(v[i][1]*omega_dot[5] + v[i][2]*omega_dot[4]);
           v[i][1] += -dthalf*v[i][2]*omega_dot[3];
         }
         v[i][0] *= factor[0];
         v[i][1] *= factor[1];
         v[i][2] *= factor[2];
       }
     }
   } else if (which == BIAS) {
     for (int i = 0; i < nlocal; i++) {
       if (mask[i] & groupbit) {
         temperature->remove_bias(i,v[i]);
         v[i][0] *= factor[0];
         v[i][1] *= factor[1];
         v[i][2] *= factor[2];
         if (pstyle == TRICLINIC) {
           v[i][0] += -dthalf*(v[i][1]*omega_dot[5] + v[i][2]*omega_dot[4]);
           v[i][1] += -dthalf*v[i][2]*omega_dot[3];
         }
         v[i][0] *= factor[0];
         v[i][1] *= factor[1];
         v[i][2] *= factor[2];
         temperature->restore_bias(i,v[i]);
       }
     }
   }
 }
 
 /* ----------------------------------------------------------------------
    perform half-step update of velocities
 -----------------------------------------------------------------------*/
 
 void FixNH::nve_v()
 {
   double dtfm;
   double **v = atom->v;
   double **f = atom->f;
   double *rmass = atom->rmass;
   double *mass = atom->mass;
   int *type = atom->type;
   int *mask = atom->mask;
   int nlocal = atom->nlocal;
   if (igroup == atom->firstgroup) nlocal = atom->nfirst;
 
   if (rmass) {
     for (int i = 0; i < nlocal; i++) {
       if (mask[i] & groupbit) {
         dtfm = dtf / rmass[i];
         v[i][0] += dtfm*f[i][0];
         v[i][1] += dtfm*f[i][1];
         v[i][2] += dtfm*f[i][2];
       }
     }
   } else {
     for (int i = 0; i < nlocal; i++) {
       if (mask[i] & groupbit) {
         dtfm = dtf / mass[type[i]];
         v[i][0] += dtfm*f[i][0];
         v[i][1] += dtfm*f[i][1];
         v[i][2] += dtfm*f[i][2];
       }
     }
   }
 }
 
 /* ----------------------------------------------------------------------
    perform full-step update of positions
 -----------------------------------------------------------------------*/
 
 void FixNH::nve_x()
 {
   double **x = atom->x;
   double **v = atom->v;
   int *mask = atom->mask;
   int nlocal = atom->nlocal;
   if (igroup == atom->firstgroup) nlocal = atom->nfirst;
 
   // x update by full step only for atoms in group
 
   for (int i = 0; i < nlocal; i++) {
     if (mask[i] & groupbit) {
       x[i][0] += dtv * v[i][0];
       x[i][1] += dtv * v[i][1];
       x[i][2] += dtv * v[i][2];
     }
   }
 }
 
 /* ----------------------------------------------------------------------
    perform half-step thermostat scaling of velocities
 -----------------------------------------------------------------------*/
 
 void FixNH::nh_v_temp()
 {
   double **v = atom->v;
   int *mask = atom->mask;
   int nlocal = atom->nlocal;
   if (igroup == atom->firstgroup) nlocal = atom->nfirst;
 
   if (which == NOBIAS) {
     for (int i = 0; i < nlocal; i++) {
       if (mask[i] & groupbit) {
         v[i][0] *= factor_eta;
         v[i][1] *= factor_eta;
         v[i][2] *= factor_eta;
       }
     }
   } else if (which == BIAS) {
     for (int i = 0; i < nlocal; i++) {
       if (mask[i] & groupbit) {
         temperature->remove_bias(i,v[i]);
         v[i][0] *= factor_eta;
         v[i][1] *= factor_eta;
         v[i][2] *= factor_eta;
         temperature->restore_bias(i,v[i]);
       }
     }
   }
 }
 
 /* ----------------------------------------------------------------------
    compute sigma tensor
    needed whenever p_target or h0_inv changes
 -----------------------------------------------------------------------*/
 
 void FixNH::compute_sigma()
 {
   // if nreset_h0 > 0, reset vol0 and h0_inv
   // every nreset_h0 timesteps
 
   if (nreset_h0 > 0) {
     int delta = update->ntimestep - update->beginstep;
     if (delta % nreset_h0 == 0) {
       if (dimension == 3) vol0 = domain->xprd * domain->yprd * domain->zprd;
       else vol0 = domain->xprd * domain->yprd;
       h0_inv[0] = domain->h_inv[0];
       h0_inv[1] = domain->h_inv[1];
       h0_inv[2] = domain->h_inv[2];
       h0_inv[3] = domain->h_inv[3];
       h0_inv[4] = domain->h_inv[4];
       h0_inv[5] = domain->h_inv[5];
     }
   }
 
   // generate upper-triangular half of
   // sigma = vol0*h0inv*(p_target-p_hydro)*h0inv^t
   // units of sigma are are PV/L^2 e.g. atm.A
   //
   // [ 0 5 4 ]   [ 0 5 4 ] [ 0 5 4 ] [ 0 - - ]
   // [ 5 1 3 ] = [ - 1 3 ] [ 5 1 3 ] [ 5 1 - ]
   // [ 4 3 2 ]   [ - - 2 ] [ 4 3 2 ] [ 4 3 2 ]
 
   sigma[0] =
     vol0*(h0_inv[0]*((p_target[0]-p_hydro)*h0_inv[0] +
                      p_target[5]*h0_inv[5]+p_target[4]*h0_inv[4]) +
           h0_inv[5]*(p_target[5]*h0_inv[0] +
                      (p_target[1]-p_hydro)*h0_inv[5]+p_target[3]*h0_inv[4]) +
           h0_inv[4]*(p_target[4]*h0_inv[0]+p_target[3]*h0_inv[5] +
                      (p_target[2]-p_hydro)*h0_inv[4]));
   sigma[1] =
     vol0*(h0_inv[1]*((p_target[1]-p_hydro)*h0_inv[1] +
                      p_target[3]*h0_inv[3]) +
           h0_inv[3]*(p_target[3]*h0_inv[1] +
                      (p_target[2]-p_hydro)*h0_inv[3]));
   sigma[2] =
     vol0*(h0_inv[2]*((p_target[2]-p_hydro)*h0_inv[2]));
   sigma[3] =
     vol0*(h0_inv[1]*(p_target[3]*h0_inv[2]) +
           h0_inv[3]*((p_target[2]-p_hydro)*h0_inv[2]));
   sigma[4] =
     vol0*(h0_inv[0]*(p_target[4]*h0_inv[2]) +
           h0_inv[5]*(p_target[3]*h0_inv[2]) +
           h0_inv[4]*((p_target[2]-p_hydro)*h0_inv[2]));
   sigma[5] =
     vol0*(h0_inv[0]*(p_target[5]*h0_inv[1]+p_target[4]*h0_inv[3]) +
           h0_inv[5]*((p_target[1]-p_hydro)*h0_inv[1]+p_target[3]*h0_inv[3]) +
           h0_inv[4]*(p_target[3]*h0_inv[1]+(p_target[2]-p_hydro)*h0_inv[3]));
 }
 
 /* ----------------------------------------------------------------------
    compute strain energy
 -----------------------------------------------------------------------*/
 
 double FixNH::compute_strain_energy()
 {
   // compute strain energy = 0.5*Tr(sigma*h*h^t) in energy units
 
   double* h = domain->h;
   double d0,d1,d2;
 
   d0 =
     sigma[0]*(h[0]*h[0]+h[5]*h[5]+h[4]*h[4]) +
     sigma[5]*(          h[1]*h[5]+h[3]*h[4]) +
     sigma[4]*(                    h[2]*h[4]);
   d1 =
     sigma[5]*(          h[5]*h[1]+h[4]*h[3]) +
     sigma[1]*(          h[1]*h[1]+h[3]*h[3]) +
     sigma[3]*(                    h[2]*h[3]);
   d2 =
     sigma[4]*(                    h[4]*h[2]) +
     sigma[3]*(                    h[3]*h[2]) +
     sigma[2]*(                    h[2]*h[2]);
 
   double energy = 0.5*(d0+d1+d2)/nktv2p;
   return energy;
 }
 
 /* ----------------------------------------------------------------------
    compute deviatoric barostat force = h*sigma*h^t
 -----------------------------------------------------------------------*/
 
 void FixNH::compute_deviatoric()
 {
   // generate upper-triangular part of h*sigma*h^t
   // units of fdev are are PV, e.g. atm*A^3
   // [ 0 5 4 ]   [ 0 5 4 ] [ 0 5 4 ] [ 0 - - ]
   // [ 5 1 3 ] = [ - 1 3 ] [ 5 1 3 ] [ 5 1 - ]
   // [ 4 3 2 ]   [ - - 2 ] [ 4 3 2 ] [ 4 3 2 ]
 
   double* h = domain->h;
 
   fdev[0] =
     h[0]*(sigma[0]*h[0]+sigma[5]*h[5]+sigma[4]*h[4]) +
     h[5]*(sigma[5]*h[0]+sigma[1]*h[5]+sigma[3]*h[4]) +
     h[4]*(sigma[4]*h[0]+sigma[3]*h[5]+sigma[2]*h[4]);
   fdev[1] =
     h[1]*(              sigma[1]*h[1]+sigma[3]*h[3]) +
     h[3]*(              sigma[3]*h[1]+sigma[2]*h[3]);
   fdev[2] =
     h[2]*(                            sigma[2]*h[2]);
   fdev[3] =
     h[1]*(                            sigma[3]*h[2]) +
     h[3]*(                            sigma[2]*h[2]);
   fdev[4] =
     h[0]*(                            sigma[4]*h[2]) +
     h[5]*(                            sigma[3]*h[2]) +
     h[4]*(                            sigma[2]*h[2]);
   fdev[5] =
     h[0]*(              sigma[5]*h[1]+sigma[4]*h[3]) +
     h[5]*(              sigma[1]*h[1]+sigma[3]*h[3]) +
     h[4]*(              sigma[3]*h[1]+sigma[2]*h[3]);
 }
 
 /* ----------------------------------------------------------------------
    compute target temperature and kinetic energy
 -----------------------------------------------------------------------*/
 
 void FixNH::compute_temp_target()
 {
   double delta = update->ntimestep - update->beginstep;
   if (delta != 0.0) delta /= update->endstep - update->beginstep;
 
   t_target = t_start + delta * (t_stop-t_start);
   ke_target = tdof * boltz * t_target;
 }
 
 /* ----------------------------------------------------------------------
    compute hydrostatic target pressure
 -----------------------------------------------------------------------*/
 
 void FixNH::compute_press_target()
 {
   double delta = update->ntimestep - update->beginstep;
   if (delta != 0.0) delta /= update->endstep - update->beginstep;
 
   p_hydro = 0.0;
   for (int i = 0; i < 3; i++)
     if (p_flag[i]) {
       p_target[i] = p_start[i] + delta * (p_stop[i]-p_start[i]);
       p_hydro += p_target[i];
     }
   if (pdim > 0) p_hydro /= pdim;
 
   if (pstyle == TRICLINIC)
     for (int i = 3; i < 6; i++)
       p_target[i] = p_start[i] + delta * (p_stop[i]-p_start[i]);
 
   // if deviatoric, recompute sigma each time p_target changes
 
   if (deviatoric_flag) compute_sigma();
 }
 
 /* ----------------------------------------------------------------------
    update omega_dot, omega
 -----------------------------------------------------------------------*/
 
 void FixNH::nh_omega_dot()
 {
   double f_omega,volume;
 
   if (dimension == 3) volume = domain->xprd*domain->yprd*domain->zprd;
   else volume = domain->xprd*domain->yprd;
 
   if (deviatoric_flag) compute_deviatoric();
 
   mtk_term1 = 0.0;
   if (mtk_flag) {
     if (pstyle == ISO) {
       mtk_term1 = tdof * boltz * t_current;
       mtk_term1 /= pdim * atom->natoms;
     } else {
       double *mvv_current = temperature->vector;
       for (int i = 0; i < 3; i++)
         if (p_flag[i])
           mtk_term1 += mvv_current[i];
       mtk_term1 /= pdim * atom->natoms;
     }
   }
 
   for (int i = 0; i < 3; i++)
     if (p_flag[i]) {
       f_omega = (p_current[i]-p_hydro)*volume /
         (omega_mass[i] * nktv2p) + mtk_term1 / omega_mass[i];
       if (deviatoric_flag) f_omega -= fdev[i]/(omega_mass[i] * nktv2p);
       omega_dot[i] += f_omega*dthalf;
       omega_dot[i] *= pdrag_factor;
     }
 
   mtk_term2 = 0.0;
   if (mtk_flag) {
     for (int i = 0; i < 3; i++)
       if (p_flag[i])
         mtk_term2 += omega_dot[i];
     if (pdim > 0) mtk_term2 /= pdim * atom->natoms;
   }
 
   if (pstyle == TRICLINIC) {
     for (int i = 3; i < 6; i++) {
       if (p_flag[i]) {
         f_omega = p_current[i]*volume/(omega_mass[i] * nktv2p);
         if (deviatoric_flag)
           f_omega -= fdev[i]/(omega_mass[i] * nktv2p);
         omega_dot[i] += f_omega*dthalf;
         omega_dot[i] *= pdrag_factor;
       }
     }
   }
 }
 
 /* ----------------------------------------------------------------------
   if any tilt ratios exceed limits, set flip = 1 and compute new tilt values
   do not flip in x or y if non-periodic (can tilt but not flip)
     this is b/c the box length would be changed (dramatically) by flip
   if yz tilt exceeded, adjust C vector by one B vector
   if xz tilt exceeded, adjust C vector by one A vector
   if xy tilt exceeded, adjust B vector by one A vector
   check yz first since it may change xz, then xz check comes after
   if any flip occurs, create new box in domain
   image_flip() adjusts image flags due to box shape change induced by flip
   remap() puts atoms outside the new box back into the new box
   perform irregular on atoms in lamda coords to migrate atoms to new procs
   important that image_flip comes before remap, since remap may change
     image flags to new values, making eqs in doc of Domain:image_flip incorrect
 ------------------------------------------------------------------------- */
 
 void FixNH::pre_exchange()
 {
   double xprd = domain->xprd;
   double yprd = domain->yprd;
 
   // flip is only triggered when tilt exceeds 0.5 by DELTAFLIP
   // this avoids immediate re-flipping due to tilt oscillations
 
   double xtiltmax = (0.5+DELTAFLIP)*xprd;
   double ytiltmax = (0.5+DELTAFLIP)*yprd;
 
   int flipxy,flipxz,flipyz;
   flipxy = flipxz = flipyz = 0;
 
   if (domain->yperiodic) {
     if (domain->yz < -ytiltmax) {
       domain->yz += yprd;
       domain->xz += domain->xy;
       flipyz = 1;
     } else if (domain->yz >= ytiltmax) {
       domain->yz -= yprd;
       domain->xz -= domain->xy;
       flipyz = -1;
     }
   }
 
   if (domain->xperiodic) {
     if (domain->xz < -xtiltmax) {
       domain->xz += xprd;
       flipxz = 1;
     } else if (domain->xz >= xtiltmax) {
       domain->xz -= xprd;
       flipxz = -1;
     }
     if (domain->xy < -xtiltmax) {
       domain->xy += xprd;
       flipxy = 1;
     } else if (domain->xy >= xtiltmax) {
       domain->xy -= xprd;
       flipxy = -1;
     }
   }
 
   int flip = 0;
   if (flipxy || flipxz || flipyz) flip = 1;
 
   if (flip) {
     domain->set_global_box();
     domain->set_local_box();
 
     domain->image_flip(flipxy,flipxz,flipyz);
 
     double **x = atom->x;
     imageint *image = atom->image;
     int nlocal = atom->nlocal;
     for (int i = 0; i < nlocal; i++) domain->remap(x[i],image[i]);
 
     domain->x2lamda(atom->nlocal);
     irregular->migrate_atoms();
     domain->lamda2x(atom->nlocal);
   }
 }
 
 /* ----------------------------------------------------------------------
    memory usage of Irregular
 ------------------------------------------------------------------------- */
 
 double FixNH::memory_usage()
 {
   double bytes = 0.0;
   if (irregular) bytes += irregular->memory_usage();
   return bytes;
 }
diff --git a/src/version.h b/src/version.h
index c55689c66..c6df1e417 100644
--- a/src/version.h
+++ b/src/version.h
@@ -1 +1 @@
-#define LAMMPS_VERSION "20 Sep 2016"
+#define LAMMPS_VERSION "22 Sep 2016"