sync 13Feb17 patch back to GH

Merge pull request #372 from akohlmey/fft-cleanup
simplify FFT3d code by removing support for outdated FFT libraries
2017-02-13 09:05:03 -07:00 · 2017-02-13 08:53:13 -07:00 · 2017-02-13 08:40:54 -07:00 · 2017-02-10 09:25:32 -07:00 · 2017-02-10 09:23:51 -07:00 · 2017-02-10 09:22:22 -07:00
219 changed files with 15659 additions and 1535 deletions
--- a/doc/src/Eqs/bond_oxdna_fene.jpg
+++ b/doc/src/Eqs/bond_oxdna_fene.jpg
--- a/doc/src/Eqs/bond_oxdna_fene.tex
+++ b/doc/src/Eqs/bond_oxdna_fene.tex
@ -0,0 +1,10 @@
+\documentclass[12pt]{article}
+\pagestyle{empty}
+
+\begin{document}
+
+$$ 
+  E = - \frac{\epsilon}{2} \ln \left[ 1 - \left(\frac{r-r0}{\Delta}\right)^2\right]
+$$
+
+\end{document}
--- a/doc/src/Manual.txt
+++ b/doc/src/Manual.txt
@ -1,7 +1,7 @@
 <!-- HTML_ONLY -->
 <HEAD>
 <TITLE>LAMMPS Users Manual</TITLE>
-<META NAME="docnumber" CONTENT="17 Jan 2017 version">
+<META NAME="docnumber" CONTENT="13 Feb 2017 version">
 <META NAME="author" CONTENT="http://lammps.sandia.gov - Sandia National Laboratories">
 <META NAME="copyright" CONTENT="Copyright (2003) Sandia Corporation.  This software and manual is distributed under the GNU General Public License.">
 </HEAD>
@ -21,7 +21,7 @@
 <H1></H1>

 LAMMPS Documentation :c,h3
-17 Jan 2017 version :c,h4
+13 Feb 2017 version :c,h4

 Version info: :h4

--- a/doc/src/PDF/USER-CGDNA-overview.pdf
+++ b/doc/src/PDF/USER-CGDNA-overview.pdf
--- a/doc/src/Section_commands.txt
+++ b/doc/src/Section_commands.txt
@ -702,6 +702,8 @@ package"_Section_start.html#start_3.
 "meso"_fix_meso.html,
 "manifoldforce"_fix_manifoldforce.html,
 "meso/stationary"_fix_meso_stationary.html,
+"nve/dot"_fix_nve_dot.html,
+"nve/dotc/langevin"_fix_nve_dotc_langevin.html,
 "nve/manifold/rattle"_fix_nve_manifold_rattle.html,
 "nvk"_fix_nvk.html,
 "nvt/manifold/rattle"_fix_nvt_manifold_rattle.html,
@ -967,7 +969,7 @@ KOKKOS, o = USER-OMP, t = OPT.
 "lubricateU/poly"_pair_lubricateU.html,
 "meam"_pair_meam.html,
 "mie/cut (o)"_pair_mie.html,
-"morse (got)"_pair_morse.html,
+"morse (gkot)"_pair_morse.html,
 "nb3b/harmonic (o)"_pair_nb3b_harmonic.html,
 "nm/cut (o)"_pair_nm.html,
 "nm/cut/coul/cut (o)"_pair_nm.html,
@ -1035,6 +1037,11 @@ package"_Section_start.html#start_3.
 "morse/soft"_pair_morse.html,
 "multi/lucy"_pair_multi_lucy.html,
 "multi/lucy/rx"_pair_multi_lucy_rx.html,
+"oxdna/coaxstk"_pair_oxdna.html,
+"oxdna/excv"_pair_oxdna.html,
+"oxdna/hbond"_pair_oxdna.html,
+"oxdna/stk"_pair_oxdna.html,
+"oxdna/xstk"_pair_oxdna.html,
 "quip"_pair_quip.html,
 "reax/c (k)"_pair_reax_c.html,
 "smd/hertz"_pair_smd_hertz.html,
@ -1083,7 +1090,8 @@ if "LAMMPS is built with the appropriate
 package"_Section_start.html#start_3.

 "harmonic/shift (o)"_bond_harmonic_shift.html,
-"harmonic/shift/cut (o)"_bond_harmonic_shift_cut.html :tb(c=4,ea=c)
+"harmonic/shift/cut (o)"_bond_harmonic_shift_cut.html,
+"oxdna/fene"_bond_oxdna_fene.html :tb(c=4,ea=c)

 :line

--- a/doc/src/Section_errors.txt
+++ b/doc/src/Section_errors.txt
@ -22,7 +22,7 @@ either conceptually, or as printed out by the program.

 12.1 Common problems :link(err_1),h4

-If two LAMMPS runs do not produce the same answer on different
+If two LAMMPS runs do not produce the exact same answer on different
 machines or different numbers of processors, this is typically not a
 bug.  In theory you should get identical answers on any number of
 processors and on any machine.  In practice, numerical round-off can
@ -80,12 +80,24 @@ order.  If you mess this up, LAMMPS will often flag the error, but it
 may also simply read a bogus argument and assign a value that is
 valid, but not what you wanted.  E.g. trying to read the string "abc"
 as an integer value of 0.  Careful reading of the associated doc page
-for the command should allow you to fix these problems.  Note that
-some commands allow for variables to be specified in place of numeric
-constants so that the value can be evaluated and change over the
-course of a run.  This is typically done with the syntax {v_name} for
-a parameter, where name is the name of the variable.  This is only
-allowed if the command documentation says it is.
+for the command should allow you to fix these problems. In most cases,
+where LAMMPS expects to read a number, either integer or floating point,
+it performs a stringent test on whether the provided input actually
+is an integer or floating-point number, respectively, and reject the
+input with an error message (for instance, when an integer is required,
+but a floating-point number 1.0 is provided):
+
+ERROR: Expected integer parameter in input script or data file :pre
+
+Some commands allow for using variable references in place of numeric
+constants so that the value can be evaluated and may change over the
+course of a run.  This is typically done with the syntax {v_name} for a
+parameter, where name is the name of the variable. On the other hand,
+immediate variable expansion with the syntax ${name} is performed while
+reading the input and before parsing commands,
+
+NOTE: Using a variable reference (i.e. {v_name}) is only allowed if
+the documentation of the corresponding command explicitly says it is.

 Generally, LAMMPS will print a message to the screen and logfile and
 exit gracefully when it encounters a fatal error.  Sometimes it will
--- a/doc/src/Section_packages.txt
+++ b/doc/src/Section_packages.txt
@ -84,7 +84,7 @@ Package, Description, Author(s), Doc page, Example, Library
 "PERI"_#PERI, Peridynamics models, Mike Parks (Sandia), "pair_style peri"_pair_peri.html, peri, -
 "POEMS"_#POEMS, coupled rigid body motion, Rudra Mukherjee (JPL), "fix poems"_fix_poems.html, rigid, lib/poems
 "PYTHON"_#PYTHON, embed Python code in an input script, -, "python"_python.html, python, lib/python
-"REAX"_#REAX, ReaxFF potential, Aidan Thompson (Sandia), "pair_style reax"_pair_reax.html, reax,  lib/reax
+"REAX"_#REAX, ReaxFF potential, Aidan Thompson (Sandia), "pair_style reax"_pair_reax.html, reax, lib/reax
 "REPLICA"_#REPLICA, multi-replica methods, -, "Section 6.6.5"_Section_howto.html#howto_5, tad, -
 "RIGID"_#RIGID, rigid bodies, -, "fix rigid"_fix_rigid.html, rigid, -
 "SHOCK"_#SHOCK, shock loading methods, -, "fix msst"_fix_msst.html, -, -
@ -1140,6 +1140,7 @@ Package, Description, Author(s), Doc page, Example, Pic/movie, Library
 "USER-ATC"_#USER-ATC, atom-to-continuum coupling, Jones & Templeton & Zimmerman (1), "fix atc"_fix_atc.html, USER/atc, "atc"_atc, lib/atc
 "USER-AWPMD"_#USER-AWPMD, wave-packet MD, Ilya Valuev (JIHT), "pair_style awpmd/cut"_pair_awpmd.html, USER/awpmd, -, lib/awpmd
 "USER-CG-CMM"_#USER-CG-CMM, coarse-graining model, Axel Kohlmeyer (Temple U), "pair_style lj/sdk"_pair_sdk.html, USER/cg-cmm, "cg"_cg, -
+"USER-CGDNA"_#USER-CGDNA, coarse-grained DNA force fields, Oliver Henrich (U Edinburgh), src/USER-CGDNA/README, USER/cgdna, -, -
 "USER-COLVARS"_#USER-COLVARS, collective variables, Fiorin & Henin & Kohlmeyer (2), "fix colvars"_fix_colvars.html, USER/colvars, "colvars"_colvars, lib/colvars
 "USER-DIFFRACTION"_#USER-DIFFRACTION, virutal x-ray and electron diffraction, Shawn Coleman (ARL),"compute xrd"_compute_xrd.html, USER/diffraction, -, -
 "USER-DPD"_#USER-DPD, reactive dissipative particle dynamics (DPD), Larentzos & Mattox & Brennan (5), src/USER-DPD/README, USER/dpd, -, -
@ -1284,6 +1285,31 @@ him directly if you have questions.

 :line

+USER-CGDNA package :link(USER-CGDNA),h5
+
+Contents: The CGDNA package implements coarse-grained force fields for
+single- and double-stranded DNA. This is at the moment mainly the
+oxDNA model, developed by Doye, Louis and Ouldridge at the University
+of Oxford.  The package also contains Langevin-type rigid-body
+integrators with improved stability.
+
+See these doc pages to get started:
+
+"bond_style oxdna_fene"_bond_oxdna_fene.html
+"pair_style oxdna_excv"_pair_oxdna_excv.html
+"fix nve/dotc/langevin"_fix_nve_dotc_langevin.html :ul
+
+Supporting info: /src/USER-CGDNA/README, "bond_style
+oxdna_fene"_bond_oxdna_fene.html, "pair_style
+oxdna_excv"_pair_oxdna_excv.html, "fix
+nve/dotc/langevin"_fix_nve_dotc_langevin.html
+
+Author: Oliver Henrich at the University of Edinburgh, UK (o.henrich
+at epcc.ed.ac.uk or ohenrich at ph.ed.ac.uk).  Contact him directly if
+you have any questions.
+
+:line
+
 USER-COLVARS package :link(USER-COLVARS),h5

 Contents: COLVARS stands for collective variables which can be used to
--- a/doc/src/Section_start.txt
+++ b/doc/src/Section_start.txt
@ -413,7 +413,7 @@ uses (for performing 1d FFTs) when running the particle-particle
 particle-mesh (PPPM) option for long-range Coulombics via the
 "kspace_style"_kspace_style.html command.

-LAMMPS supports various open-source or vendor-supplied FFT libraries
+LAMMPS supports common open-source or vendor-supplied FFT libraries
 for this purpose.  If you leave these 3 variables blank, LAMMPS will
 use the open-source "KISS FFT library"_http://kissfft.sf.net, which is
 included in the LAMMPS distribution.  This library is portable to all
@ -423,10 +423,9 @@ package in your build, you can also leave the 3 variables blank.

 Otherwise, select which kinds of FFTs to use as part of the FFT_INC
 setting by a switch of the form -DFFT_XXX.  Recommended values for XXX
-are: MKL, SCSL, FFTW2, and FFTW3.  Legacy options are: INTEL, SGI,
-ACML, and T3E.  For backward compatability, using -DFFT_FFTW will use
-the FFTW2 library.  Using -DFFT_NONE will use the KISS library
-described above.
+are: MKL or FFTW3.  FFTW2 and NONE are supported as legacy options.
+Selecting -DFFT_FFTW will use the FFTW3 library and -DFFT_NONE will
+use the KISS library described above.

 You may also need to set the FFT_INC, FFT_PATH, and FFT_LIB variables,
 so the compiler and linker can find the needed FFT header and library
--- a/doc/src/bond_oxdna_fene.txt
+++ b/doc/src/bond_oxdna_fene.txt
@ -0,0 +1,70 @@
+"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
+
+bond_style oxdna_fene command :h3
+
+[Syntax:]
+
+bond_style oxdna_fene :pre
+
+[Examples:]
+
+bond_style oxdna_fene
+bond_coeff * 2.0 0.25 0.7525 :pre
+
+[Description:]
+
+The {oxdna_fene} bond style uses the potential
+
+:c,image(Eqs/bond_oxdna_fene.jpg)
+
+to define a modified finite extensible nonlinear elastic (FENE) potential
+"(Ouldridge)"_#oxdna_fene to model the connectivity of the phosphate backbone
+in the oxDNA force field for coarse-grained modelling of DNA. 
+
+The following coefficients must be defined for the bond type via the
+"bond_coeff"_bond_coeff.html command as given in the above example, or in
+the data file or restart files read by the "read_data"_read_data.html
+or "read_restart"_read_restart.html commands:
+
+epsilon (energy)
+Delta (distance)
+r0 (distance) :ul
+
+NOTE: This bond style has to be used together with the corresponding oxDNA pair styles
+for excluded volume interaction {oxdna_excv}, stacking {oxdna_stk}, cross-stacking {oxdna_xstk}
+and coaxial stacking interaction {oxdna_coaxstk} as well as hydrogen-bonding interaction {oxdna_hbond} (see also documentation of 
+"pair_style oxdna_excv"_pair_oxdna_excv.html). The coefficients 
+in the above example have to be kept fixed and cannot be changed without reparametrizing the entire model.
+
+Example input and data files can be found in /examples/USER/cgdna/examples/duplex1/ and /duplex2/.
+A simple python setup tool which creates single straight or helical DNA strands,
+DNA duplexes or arrays of DNA duplexes can be found in /examples/USER/cgdna/util/.
+A technical report with more information on the model, the structure of the input file,
+the setup tool and the performance of the LAMMPS-implementation of oxDNA
+can be found "here"_PDF/USER-CGDNA-overview.pdf.
+
+:line
+
+[Restrictions:]
+
+This bond style can only be used if LAMMPS was built with the
+USER-CGDNA package and the MOLECULE and ASPHERE package.  See the "Making
+LAMMPS"_Section_start.html#start_3 section for more info on packages.
+
+
+[Related commands:]
+
+"pair_style oxdna_excv"_pair_oxdna_excv.html, "fix nve/dotc/langevin"_fix_nve_dotc_langevin.html, "bond_coeff"_bond_coeff.html 
+
+[Default:] none
+
+:line
+
+:link(oxdna_fene)
+[(Ouldridge)] T.E. Ouldridge, A.A. Louis, J.P.K. Doye, J. Chem. Phys. 134, 085101 (2011).
--- a/doc/src/bonds.txt
+++ b/doc/src/bonds.txt
@ -15,6 +15,7 @@ Bond Styles :h1
   bond_morse
   bond_none
   bond_nonlinear
+   bond_oxdna_fene
   bond_quartic
   bond_table
   bond_zero
--- a/doc/src/compute_coord_atom.txt
+++ b/doc/src/compute_coord_atom.txt
@ -12,19 +12,16 @@ compute coord/atom command :h3

 compute ID group-ID coord/atom cstyle args ... :pre

-ID, group-ID are documented in "compute"_compute.html command
-coord/atom = style name of this compute command
-one cstyle must be appended :ul
-
-cstyle = {cutoff} or {orientorder}
-
-{cutoff} args = cutoff typeN
-  cutoff = distance within which to count coordination neighbors (distance units)
-  typeN = atom type for Nth coordination count (see asterisk form below) :pre
-
-{orientorder} args = orientorderID threshold
-  orientorderID = ID of a previously defined orientorder/atom compute
-  threshold = minimum value of the scalar product between two 'connected' atoms (see text for explanation) :pre
+ID, group-ID are documented in "compute"_compute.html command :ulb,l
+coord/atom = style name of this compute command :l
+cstyle = {cutoff} or {orientorder} :l
+  {cutoff} args = cutoff typeN
+    cutoff = distance within which to count coordination neighbors (distance units)
+    typeN = atom type for Nth coordination count (see asterisk form below)
+  {orientorder} args = orientorderID threshold
+    orientorderID = ID of an orientorder/atom compute
+    threshold = minimum value of the product of two "connected" atoms :pre
+:ule

 [Examples:]

@ -35,21 +32,21 @@ compute 1 all coord/atom orientorder 2 0.5 :pre

 [Description:]

-This compute performs generic calculations between neighboring atoms. So far,
-there are two cstyles implemented: {cutoff} and {orientorder}.
-The {cutoff} cstyle calculates one or more coordination numbers
-for each atom in a group.
+This compute performs calculations between neighboring atoms to
+determine a coordination value.  The specific calculation and the
+meaning of the resulting value depend on the {cstyle} keyword used.

-A coordination number is defined as the number of neighbor atoms with
-specified atom type(s) that are within the specified cutoff distance
-from the central atom.  Atoms not in the group are included in a
-coordination number of atoms in the group.
+The {cutoff} cstyle calculates one or more traditional coordination
+numbers for each atom.  A coordination number is defined as the number
+of neighbor atoms with specified atom type(s) that are within the
+specified cutoff distance from the central atom.  Atoms not in the
+specified group are included in the coordination number tally.

-The {typeN} keywords allow you to specify which atom types contribute
-to each coordination number.  One coordination number is computed for
-each of the {typeN} keywords listed.  If no {typeN} keywords are
-listed, a single coordination number is calculated, which includes
-atoms of all types (same as the "*" format, see below).
+The {typeN} keywords allow specification of which atom types
+contribute to each coordination number.  One coordination number is
+computed for each of the {typeN} keywords listed.  If no {typeN}
+keywords are listed, a single coordination number is calculated, which
+includes atoms of all types (same as the "*" format, see below).

 The {typeN} keywords can be specified in one of two ways.  An explicit
 numeric value can be used, as in the 2nd example above.  Or a
@ -61,16 +58,27 @@ from 1 to N.  A leading asterisk means all types from 1 to n
 (inclusive).  A middle asterisk means all types from m to n
 (inclusive).

-The {orientorder} cstyle calculates the number of 'connected' atoms j 
-around each atom i. The atom j is connected to i if the scalar product
-({Ybar_lm(i)},{Ybar_lm(j)}) is larger than {threshold}. Thus, this cstyle
-will work only if a "compute orientorder/atom"_compute_orientorder_atom.html
-has been previously defined. This cstyle allows one to apply the 
-ten Wolde's criterion to identify cristal-like atoms in a system 
-(see "ten Wolde et al."_#tenWolde).
+The {orientorder} cstyle calculates the number of "connected" neighbor
+atoms J around each central atom I.  For this {cstyle}, connected is
+defined by the orientational order parameter calculated by the
+"compute orientorder/atom"_compute_orientorder_atom.html command.
+This {cstyle} thus allows one to apply the ten Wolde's criterion to
+identify crystal-like atoms in a system, as discussed in "ten
+Wolde"_#tenWolde.

-The value of all coordination numbers will be 0.0 for atoms not in the
-specified compute group.
+The ID of the previously specified "compute
+orientorder/atom"_compute_orientorder/atom command is specified as
+{orientorderID}.  The compute must invoke its {components} option to
+calculate components of the {Ybar_lm} vector for each atoms, as
+described in its documenation.  Note that orientorder/atom compute
+defines its own criteria for identifying neighboring atoms.  If the
+scalar product ({Ybar_lm(i)},{Ybar_lm(j)}), calculated by the
+orientorder/atom compute is larger than the specified {threshold},
+then I and J are connected, and the coordination value of I is
+incremented by one.
+
+For all {cstyle} settings, all coordination values will be 0.0 for
+atoms not in the specified compute group.

 The neighbor list needed to compute this quantity is constructed each
 time the calculation is performed (i.e. each time a snapshot of atoms
@ -92,21 +100,23 @@ the neighbor list.

 [Output info:]

-If single {type1} keyword is specified (or if none are specified),
-this compute calculates a per-atom vector.  If multiple {typeN}
-keywords are specified, this compute calculates a per-atom array, with
-N columns.  These values can be accessed by any command that uses
-per-atom values from a compute as input.  See "Section
+For {cstyle} cutoff, this compute can calculate a per-atom vector or
+array.  If single {type1} keyword is specified (or if none are
+specified), this compute calculates a per-atom vector.  If multiple
+{typeN} keywords are specified, this compute calculates a per-atom
+array, with N columns.
+
+For {cstyle} orientorder, this compute calculates a per-atom vector.
+
+These values can be accessed by any command that uses per-atom values
+from a compute as input.  See "Section
 6.15"_Section_howto.html#howto_15 for an overview of LAMMPS output
 options.

 The per-atom vector or array values will be a number >= 0.0, as
 explained above.

-[Restrictions:]
-The cstyle {orientorder} can only be used if a 
-"compute orientorder/atom"_compute_orientorder_atom.html command
-was previously defined. Otherwise, an error message will be issued.
+[Restrictions:] none

 [Related commands:]

@ -118,4 +128,5 @@ was previously defined. Otherwise, an error message will be issued.
 :line

 :link(tenWolde)
-[(tenWolde)] P. R. ten Wolde, M. J. Ruiz-Montero, D. Frenkel, J. Chem. Phys. 104, 9932 (1996).
+[(tenWolde)] P. R. ten Wolde, M. J. Ruiz-Montero, D. Frenkel,
+J. Chem. Phys. 104, 9932 (1996).
--- a/doc/src/compute_group_group.txt
+++ b/doc/src/compute_group_group.txt
@ -16,10 +16,11 @@ ID, group-ID are documented in "compute"_compute.html command :ulb,l
 group/group = style name of this compute command :l
 group2-ID = group ID of second (or same) group :l
 zero or more keyword/value pairs may be appended :l
-keyword = {pair} or {kspace} or {boundary} :l
+keyword = {pair} or {kspace} or {boundary} or {molecule} :l
  {pair} value = {yes} or {no}
  {kspace} value = {yes} or {no}
-  {boundary} value = {yes} or {no} :pre
+  {boundary} value = {yes} or {no}
+  {molecule} value = {off} or {inter} or {intra} :pre
 :ule

 [Examples:]
@ -46,6 +47,13 @@ NOTE: The energies computed by the {pair} keyword do not include tail
 corrections, even if they are enabled via the
 "pair_modify"_pair_modify.html command.

+If the {molecule} keyword is set to {inter} or {intra} than an
+additional check is made based on the molecule IDs of the two atoms in
+each pair before including their pairwise interaction energy and
+force.  For the {inter} setting, the two atoms must be in different
+molecules.  For the {intra} setting, the two atoms must be in the same
+molecule.
+
 If the {kspace} keyword is set to {yes}, which is not the default, and
 if a "kspace_style"_kspace_style.html is defined, then the interaction
 energy will include a Kspace component which is the long-range
@ -66,6 +74,10 @@ affect the force calculation and will be zero if one or both of the
 groups are charge neutral.  This energy correction term is the same as
 that included in the regular Ewald and PPPM routines.

+NOTE: The {molecule} setting only affects the group/group
+contributions calculated by the {pair} keyword.  It does not affect
+the group/group contributions calculated by the {kspace} keyword.
+
 This compute does not calculate any bond or angle or dihedral or
 improper interactions between atoms in the two groups.

@ -78,6 +90,22 @@ work (FFTs, Ewald summation) as computing long-range forces for the
 entire system.  Thus it can be costly to invoke this compute too
 frequently.

+NOTE: If you have a bonded system, then the settings of
+"special_bonds"_special_bonds.html command can remove pairwise
+interactions between atoms in the same bond, angle, or dihedral.  This
+is the default setting for the "special_bonds"_special_bonds.html
+command, and means those pairwise interactions do not appear in the
+neighbor list.  Because this compute uses a neighbor list, it also
+means those pairs will not be included in the group/group interaction.
+This does not apply when using long-range coulomb interactions
+({coul/long}, {coul/msm}, {coul/wolf} or similar.  One way to get
+around this would be to set special_bond scaling factors to very tiny
+numbers that are not exactly zero (e.g. 1.0e-50). Another workaround
+is to write a dump file, and use the "rerun"_rerun.html command to
+compute the group/group interactions for snapshots in the dump file.
+The rerun script can use a "special_bonds"_special_bonds.html command
+that includes all pairs in the neighbor list.
+
 If you desire a breakdown of the interactions into a pairwise and
 Kspace component, simply invoke the compute twice with the appropriate
 yes/no settings for the {pair} and {kspace} keywords.  This is no more
@ -119,7 +147,8 @@ The {ewald} and {pppm} styles do.

 [Default:]

-The option defaults are pair = yes, kspace = no, and boundary = yes.
+The option defaults are pair = yes, kspace = no, boundary = yes,
+molecule = off.

 :line

--- a/doc/src/compute_orientorder_atom.txt
+++ b/doc/src/compute_orientorder_atom.txt
@ -19,7 +19,7 @@ keyword = {cutoff} or {nnn} or {degrees} or {components}
  {cutoff} value = distance cutoff
  {nnn} value = number of nearest neighbors
  {degrees} values = nlvalues, l1, l2,...
-  {components} value = l  :pre
+  {components} value = ldegree  :pre

 :ule

@ -64,21 +64,21 @@ specified distance cutoff are used.
 The optional keyword {degrees} defines the list of order parameters to
 be computed.  The first argument {nlvalues} is the number of order
 parameters. This is followed by that number of integers giving the
-degree of each order parameter. Because {Q}2 and all odd-degree
-order parameters are zero for atoms in cubic crystals
-(see "Steinhardt"_#Steinhardt), the default order parameters
-are {Q}4, {Q}6, {Q}8, {Q}10, and {Q}12. For the
-FCC crystal with {nnn}=12, {Q}4 = sqrt(7/3)/8 = 0.19094....
-The numerical values of all order parameters up to {Q}12
-for a range of commonly encountered high-symmetry structures are given
-in Table I of "Mickel et al."_#Mickel.
+degree of each order parameter. Because {Q}2 and all odd-degree order
+parameters are zero for atoms in cubic crystals (see
+"Steinhardt"_#Steinhardt), the default order parameters are {Q}4,
+{Q}6, {Q}8, {Q}10, and {Q}12. For the FCC crystal with {nnn}=12, {Q}4
+= sqrt(7/3)/8 = 0.19094....  The numerical values of all order
+parameters up to {Q}12 for a range of commonly encountered
+high-symmetry structures are given in Table I of "Mickel et
+al."_#Mickel.

-The optional keyword {components} will output the components of
-the normalized complex vector {Ybar_lm} of degree {l}, which must be
+The optional keyword {components} will output the components of the
+normalized complex vector {Ybar_lm} of degree {ldegree}, which must be
 explicitly included in the keyword {degrees}. This option can be used
 in conjunction with "compute coord_atom"_compute_coord_atom.html to
-calculate the ten Wolde's criterion to identify crystal-like particles
-(see "ten Wolde et al."_#tenWolde96).
+calculate the ten Wolde's criterion to identify crystal-like
+particles, as discussed in "ten Wolde"_#tenWolde.

 The value of {Ql} is set to zero for atoms not in the
 specified compute group, as well as for atoms that have less than
@ -104,14 +104,16 @@ the neighbor list.

 [Output info:]

-This compute calculates a per-atom array with {nlvalues} columns, giving the
-{Ql} values for each atom, which are real numbers on the range 0 <= {Ql} <= 1.
+This compute calculates a per-atom array with {nlvalues} columns,
+giving the {Ql} values for each atom, which are real numbers on the
+range 0 <= {Ql} <= 1.

-If the keyword {components} is set, then the real and imaginary parts of each
-component of (normalized) {Ybar_lm} will be added to the output array in the
-following order:
-Re({Ybar_-m}) Im({Ybar_-m}) Re({Ybar_-m+1}) Im({Ybar_-m+1}) ... Re({Ybar_m}) Im({Ybar_m}).
-This way, the per-atom array will have a total of {nlvalues}+2*(2{l}+1) columns.
+If the keyword {components} is set, then the real and imaginary parts
+of each component of (normalized) {Ybar_lm} will be added to the
+output array in the following order: Re({Ybar_-m}) Im({Ybar_-m})
+Re({Ybar_-m+1}) Im({Ybar_-m+1}) ... Re({Ybar_m}) Im({Ybar_m}).  This
+way, the per-atom array will have a total of {nlvalues}+2*(2{l}+1)
+columns.

 These values can be accessed by any command that uses
 per-atom values from a compute as input.  See "Section
@ -122,19 +124,25 @@ options.

 [Related commands:]

-"compute coord/atom"_compute_coord_atom.html, "compute centro/atom"_compute_centro_atom.html, "compute hexorder/atom"_compute_hexorder_atom.html
+"compute coord/atom"_compute_coord_atom.html, "compute
+centro/atom"_compute_centro_atom.html, "compute
+hexorder/atom"_compute_hexorder_atom.html

 [Default:]

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

 :line

 :link(Steinhardt)
-[(Steinhardt)] P. Steinhardt, D. Nelson, and M. Ronchetti, Phys. Rev. B 28, 784 (1983).
+[(Steinhardt)] P. Steinhardt, D. Nelson, and M. Ronchetti,
+Phys. Rev. B 28, 784 (1983).

 :link(Mickel)
-[(Mickel)] W. Mickel, S. C. Kapfer, G. E. Schroeder-Turkand, K. Mecke, J. Chem. Phys. 138, 044501 (2013).
+[(Mickel)] W. Mickel, S. C. Kapfer, G. E. Schroeder-Turkand, K. Mecke,
+J. Chem. Phys. 138, 044501 (2013).

-:link(tenWolde96)
-[(tenWolde)] P. R. ten Wolde, M. J. Ruiz-Montero, D. Frenkel, J. Chem. Phys. 104, 9932 (1996).
+:link(tenWolde)
+[(tenWolde)] P. R. ten Wolde, M. J. Ruiz-Montero, D. Frenkel,
+J. Chem. Phys. 104, 9932 (1996).
--- a/doc/src/compute_rdf.txt
+++ b/doc/src/compute_rdf.txt
@ -41,14 +41,14 @@ NOTE: If you have a bonded system, then the settings of
 interactions between atoms in the same bond, angle, or dihedral.  This
 is the default setting for the "special_bonds"_special_bonds.html
 command, and means those pairwise interactions do not appear in the
-neighbor list.  Because this fix uses the neighbor list, it also means
+neighbor list.  Because this fix uses a neighbor list, it also means
 those pairs will not be included in the RDF. This does not apply when
-using long-range coulomb ({coul/long}, {coul/msm}, {coul/wolf} or
-similar.  One way to get around this would be to set special_bond
-scaling factors to very tiny numbers that are not exactly zero
-(e.g. 1.0e-50). Another workaround is to write a dump file, and use
-the "rerun"_rerun.html command to compute the RDF for snapshots in the
-dump file.  The rerun script can use a
+using long-range coulomb interactions ({coul/long}, {coul/msm},
+{coul/wolf} or similar.  One way to get around this would be to set
+special_bond scaling factors to very tiny numbers that are not exactly
+zero (e.g. 1.0e-50). Another workaround is to write a dump file, and
+use the "rerun"_rerun.html command to compute the RDF for snapshots in
+the dump file.  The rerun script can use a
 "special_bonds"_special_bonds.html command that includes all pairs in
 the neighbor list.

--- a/doc/src/fix_nve_dot.txt
+++ b/doc/src/fix_nve_dot.txt
@ -0,0 +1,61 @@
+"LAMMPS WWW Site"_lws - "LAMMPS Documentation"_ld - "LAMMPS Commands"_lc :c
+
+:link(lws,http://lammps.sandia.gov)
+:link(ld,Manual.html)
+:link(lc,Section_commands.html#comm)
+
+:line
+
+fix nve/dot command :h3
+
+[Syntax:]
+
+fix ID group-ID nve/dot :pre
+
+ID, group-ID are documented in "fix"_fix.html command :ulb,l
+nve/dot = style name of this fix command :l
+:ule
+
+[Examples:]
+
+fix 1 all nve/dot :pre
+
+[Description:]
+
+Apply a rigid-body integrator as described in "(Davidchack)"_#Davidchack
+to a group of atoms, but without Langevin dynamics. 
+This command performs Molecular dynamics (MD)
+via a velocity-Verlet algorithm and an evolution operator that rotates 
+the quaternion degrees of freedom, similar to the scheme outlined in "(Miller)"_#Miller. 
+
+This command is the equivalent of the "fix nve/dotc/langevin"_fix_nve_dotc_langevin.html
+without damping and noise and can be used to determine the stability range 
+in a NVE ensemble prior to using the Langevin-type DOTC-integrator
+(see also "fix nve/dotc/langevin"_fix_nve_dotc_langevin.html).
+The command is equivalent to the "fix nve"_fix_nve.html.
+The particles are always considered to have a finite size.
+
+An example input file can be found in /examples/USER/cgdna/examples/duplex1/.
+A technical report with more information on this integrator can be found
+"here"_PDF/USER-CGDNA-overview.pdf.
+
+:line
+
+[Restrictions:]
+
+These pair styles can only be used if LAMMPS was built with the
+USER-CGDNA package and the MOLECULE and ASPHERE package.  See the "Making
+LAMMPS"_Section_start.html#start_3 section for more info on packages.
+
+[Related commands:]
+
+"fix nve/dotc/langevin"_fix_nve_dotc_langevin.html, "fix nve"_fix_nve.html
+
+[Default:] none
+
+:line
+
+:link(Davidchack)
+[(Davidchack)] R.L Davidchack, T.E. Ouldridge, and M.V. Tretyakov. J. Chem. Phys. 142, 144114 (2015).
+:link(Miller)
+[(Miller)] T. F. Miller III, M. Eleftheriou, P. Pattnaik, A. Ndirango, G. J. Martyna, J. Chem. Phys., 116, 8649-8659 (2002).
--- a/doc/src/fix_nve_dotc_langevin.txt
+++ b/doc/src/fix_nve_dotc_langevin.txt
@ -0,0 +1,134 @@
+"LAMMPS WWW Site"_lws - "LAMMPS Documentation"_ld - "LAMMPS Commands"_lc :c
+
+:link(lws,http://lammps.sandia.gov)
+:link(ld,Manual.html)
+:link(lc,Section_commands.html#comm)
+
+:line
+
+fix nve/dotc/langevin command :h3
+
+[Syntax:]
+
+fix ID group-ID nve/dotc/langevin Tstart Tstop damp seed keyword value :pre
+
+ID, group-ID are documented in "fix"_fix.html command :ulb,l
+nve/dotc/langevin = style name of this fix command :l
+Tstart,Tstop = desired temperature at start/end of run (temperature units) :l
+damp = damping parameter (time units) :l
+seed = random number seed to use for white noise (positive integer) :l
+keyword = {angmom} :l
+  {angmom} value = factor
+    factor = do thermostat rotational degrees of freedom via the angular momentum and apply numeric scale factor as discussed below :pre
+:ule
+
+[Examples:]
+
+fix 1 all nve/dotc/langevin 1.0 1.0 0.03 457145 angmom 10 :pre
+
+[Description:]
+
+Apply a rigid-body Langevin-type integrator of the kind "Langevin C" 
+as described in "(Davidchack)"_#Davidchack
+to a group of atoms, which models an interaction with an implicit background
+solvent.  This command performs Brownian dynamics (BD)
+via a technique that splits the integration into a deterministic Hamiltonian 
+part and the Ornstein-Uhlenbeck process for noise and damping. 
+The quaternion degrees of freedom are updated though an evolution
+operator which performs a rotation in quaternion space, preserves
+the quaternion norm and is akin to "(Miller)"_#Miller.
+
+In terms of syntax this command has been closely modelled on the 
+"fix langevin"_fix_langevin.html and its {angmom} option. But it combines 
+the "fix nve"_fix_nve.html and the "fix langevin"_fix_langevin.html in 
+one single command. The main feature is improved stability 
+over the standard integrator, permitting slightly larger timestep sizes.
+
+NOTE: Unlike the "fix langevin"_fix_langevin.html this command performs
+also time integration of the translational and quaternion degrees of freedom.
+
+The total force on each atom will have the form:
+
+F = Fc + Ff + Fr
+Ff = - (m / damp) v
+Fr is proportional to sqrt(Kb T m / (dt damp)) :pre
+
+Fc is the conservative force computed via the usual inter-particle
+interactions ("pair_style"_pair_style.html,
+"bond_style"_bond_style.html, etc).
+
+The Ff and Fr terms are implicitly taken into account by this fix 
+on a per-particle basis.
+
+Ff is a frictional drag or viscous damping term proportional to the
+particle's velocity.  The proportionality constant for each atom is
+computed as m/damp, where m is the mass of the particle and damp is
+the damping factor specified by the user.
+
+Fr is a force due to solvent atoms at a temperature T randomly bumping
+into the particle.  As derived from the fluctuation/dissipation
+theorem, its magnitude as shown above is proportional to sqrt(Kb T m /
+dt damp), where Kb is the Boltzmann constant, T is the desired
+temperature, m is the mass of the particle, dt is the timestep size,
+and damp is the damping factor.  Random numbers are used to randomize
+the direction and magnitude of this force as described in
+"(Dunweg)"_#Dunweg, where a uniform random number is used (instead of
+a Gaussian random number) for speed.
+
+:line
+
+{Tstart} and {Tstop} have to be constant values, i.e. they cannot 
+be variables.
+
+The {damp} parameter is specified in time units and determines how
+rapidly the temperature is relaxed.  For example, a value of 0.03
+means to relax the temperature in a timespan of (roughly) 0.03 time
+units tau (see the "units"_units.html command).
+The damp factor can be thought of as inversely related to the
+viscosity of the solvent, i.e. a small relaxation time implies a
+hi-viscosity solvent and vice versa.  See the discussion about gamma
+and viscosity in the documentation for the "fix
+viscous"_fix_viscous.html command for more details.
+
+The random # {seed} must be a positive integer. A Marsaglia random
+number generator is used.  Each processor uses the input seed to
+generate its own unique seed and its own stream of random numbers.
+Thus the dynamics of the system will not be identical on two runs on
+different numbers of processors.
+
+The keyword/value option has to be used in the following way:
+
+This fix has to be used together with the {angmom} keyword. The 
+particles are always considered to have a finite size. 
+The keyword {angmom} enables thermostatting of the rotational degrees of 
+freedom in addition to the usual translational degrees of freedom. 
+
+The scale factor after the {angmom} keyword gives the ratio of the rotational to 
+the translational friction coefficient.
+
+An example input file can be found in /examples/USER/cgdna/examples/duplex2/.
+A technical report with more information on this integrator can be found 
+"here"_PDF/USER-CGDNA-overview.pdf.
+
+:line
+
+[Restrictions:]
+
+These pair styles can only be used if LAMMPS was built with the
+USER-CGDNA package and the MOLECULE and ASPHERE package.  See the "Making
+LAMMPS"_Section_start.html#start_3 section for more info on packages.
+
+[Related commands:]
+
+"fix nve"_fix_nve.html, "fix langevin"_fix_langevin.html, "fix nve/dot"_fix_nve_dot.html,  
+
+[Default:] none
+
+:line
+
+:link(Davidchack)
+[(Davidchack)] R.L Davidchack, T.E. Ouldridge, M.V. Tretyakov. J. Chem. Phys. 142, 144114 (2015).
+:link(Miller)
+[(Miller)] T. F. Miller III, M. Eleftheriou, P. Pattnaik, A. Ndirango, G. J. Martyna, J. Chem. Phys., 116, 8649-8659 (2002).
+:link(Dunweg)
+[(Dunweg)] B. Dunweg, W. Paul, Int. J. Mod. Phys. C, 2, 817-27 (1991).
--- a/doc/src/fixes.txt
+++ b/doc/src/fixes.txt
@ -84,6 +84,8 @@ Fixes :h1
   fix_nve_asphere
   fix_nve_asphere_noforce
   fix_nve_body
+   fix_nve_dot
+   fix_nve_dotc_langevin
   fix_nve_eff
   fix_nve_limit
   fix_nve_line
--- a/doc/src/kspace_style.txt
+++ b/doc/src/kspace_style.txt
@ -229,11 +229,16 @@ dramatically in z.  For example, for a triclinic system with all three
 tilt factors set to the maximum limit, the PPPM grid should be
 increased roughly by a factor of 1.5 in the y direction and 2.0 in the
 z direction as compared to the same system using a cubic orthogonal
-simulation cell. One way to ensure the accuracy requirement is being
-met is to run a short simulation at the maximum expected tilt or
-length, note the required grid size, and then use the
+simulation cell.  One way to handle this issue if you have a long
+simulation where the box size changes dramatically, is to break it
+into shorter simulations (multiple "run"_run.html commands).  This
+works because the grid size is re-computed at the beginning of each
+run.  Another way to ensure the descired accuracy requirement is met
+is to run a short simulation at the maximum expected tilt or length,
+note the required grid size, and then use the
 "kspace_modify"_kspace_modify.html {mesh} command to manually set the
-PPPM grid size to this value.
+PPPM grid size to this value for the long run.  The simulation then
+will be "too accurate" for some portion of the run.

 RMS force errors in real space for {ewald} and {pppm} are estimated
 using equation 18 of "(Kolafa)"_#Kolafa, which is also referenced as
@ -285,6 +290,8 @@ LAMMPS"_Section_start.html#start_3 section for more info.
 See "Section 5"_Section_accelerate.html of the manual for
 more instructions on how to use the accelerated styles effectively.

+:line
+
 [Restrictions:]

 Note that the long-range electrostatic solvers in LAMMPS assume conducting
--- a/doc/src/lammps.book
+++ b/doc/src/lammps.book
@ -194,7 +194,6 @@ fix_meso.html
 fix_meso_stationary.html
 fix_momentum.html
 fix_move.html
-fix_mscg.html
 fix_msst.html
 fix_neb.html
 fix_nh.html
@ -210,6 +209,8 @@ fix_nve.html
 fix_nve_asphere.html
 fix_nve_asphere_noforce.html
 fix_nve_body.html
+fix_nve_dot.html
+fix_nve_dotc_langevin.html
 fix_nve_eff.html
 fix_nve_limit.html
 fix_nve_line.html
@ -217,7 +218,6 @@ fix_nve_manifold_rattle.html
 fix_nve_noforce.html
 fix_nve_sphere.html
 fix_nve_tri.html
-fix_nvk.html
 fix_nvt_asphere.html
 fix_nvt_body.html
 fix_nvt_manifold_rattle.html
@ -458,6 +458,7 @@ pair_multi_lucy_rx.html
 pair_nb3b_harmonic.html
 pair_nm.html
 pair_none.html
+pair_oxdna_excv.html
 pair_peri.html
 pair_polymorphic.html
 pair_quip.html
@ -496,6 +497,7 @@ pair_zero.html
 bond_class2.html
 bond_fene.html
 bond_fene_expand.html
+bond_oxdna_fene.html
 bond_harmonic.html
 bond_harmonic_shift.html
 bond_harmonic_shift_cut.html
--- a/doc/src/pair_morse.txt
+++ b/doc/src/pair_morse.txt
@ -13,6 +13,7 @@ pair_style morse/opt command :h3
 pair_style morse/smooth/linear command :h3
 pair_style morse/smooth/linear/omp command :h3
 pair_style morse/soft command :h3
+pair_style morse/kk command :h3

 [Syntax:]

--- a/doc/src/pair_oxdna_excv.txt
+++ b/doc/src/pair_oxdna_excv.txt
@ -0,0 +1,80 @@
+"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
+
+pair_style oxdna_excv command :h3
+pair_style oxdna_stk command :h3
+pair_style oxdna_hbond command :h3
+pair_style oxdna_xstk command :h3
+pair_style oxdna_coaxstk command :h3
+
+[Syntax:]
+
+pair_style style :pre
+
+style = {hybrid/overlay oxdna_excv oxdna_stk oxdna_hbond oxdna_xstk oxdna_coaxstk} :ul
+
+[Examples:]
+
+pair_style hybrid/overlay oxdna_excv oxdna_stk oxdna_hbond oxdna_xstk oxdna_coaxstk
+pair_coeff * * oxdna_excv    2.0 0.7 0.675 2.0 0.515 0.5 2.0 0.33 0.32
+pair_coeff * * oxdna_stk     1.61048 6.0 0.4 0.9 0.32 0.6 1.3 0 0.8 0.9 0 0.95 0.9 0 0.95 2.0 0.65 2.0 0.65
+pair_coeff * * oxdna_hbond   0.0   8.0 0.4 0.75 0.34 0.7 1.5 0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45
+pair_coeff 1 4 oxdna_hbond   1.077 8.0 0.4 0.75 0.34 0.7 1.5 0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45
+pair_coeff 2 3 oxdna_hbond   1.077 8.0 0.4 0.75 0.34 0.7 1.5 0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45
+pair_coeff * * oxdna_xstk    47.5 0.575 0.675 0.495 0.655 2.25 0.791592653589793 0.58 1.7 1.0 0.68 1.7 1.0 0.68 1.5 0 0.65 1.7 0.875 0.68 1.7 0.875 0.68
+pair_coeff * * oxdna_coaxstk 46.0 0.4 0.6 0.22 0.58 2.0 2.541592653589793 0.65 1.3 0 0.8 0.9 0 0.95 0.9 0 0.95 2.0 -0.65 2.0 -0.65 :pre
+
+[Description:]
+
+The {oxdna} pair styles compute the pairwise-additive parts of the oxDNA force field 
+for coarse-grained modelling of DNA. The effective interaction between the nucleotides consists of potentials for the 
+excluded volume interaction {oxdna_excv}, the stacking {oxdna_stk}, cross-stacking {oxdna_xstk}
+and coaxial stacking interaction {oxdna_coaxstk} as well
+as the hydrogen-bonding interaction {oxdna_hbond} between complementary pairs of nucleotides on
+opposite strands.
+
+The exact functional form of the pair styles is rather complex, which manifests itself in the 144 coefficients 
+in the above example. The individual potentials consist of products of modulation factors, 
+which themselves are constructed from a number of more basic potentials 
+(Morse, Lennard-Jones, harmonic angle and distance) as well as quadratic smoothing and modulation terms. 
+We refer to "(Ouldridge-DPhil)"_#Ouldridge-DPhil and "(Ouldridge)"_#Ouldridge
+for a detailed description of the oxDNA force field.
+
+NOTE: These pair styles have to be used together with the related oxDNA bond style
+{oxdna_fene} for the connectivity of the phosphate backbone (see also documentation of
+"bond_style oxdna_fene"_bond_oxdna_fene.html). The coefficients
+in the above example have to be kept fixed and cannot be changed without reparametrizing the entire model.
+
+Example input and data files can be found in /examples/USER/cgdna/examples/duplex1/ and /duplex2/.
+A simple python setup tool which creates single straight or helical DNA strands, 
+DNA duplexes or arrays of DNA duplexes can be found in /examples/USER/cgdna/util/.
+A technical report with more information on the model, the structure of the input file,
+the setup tool and the performance of the LAMMPS-implementation of oxDNA 
+can be found "here"_PDF/USER-CGDNA-overview.pdf.
+
+:line
+
+[Restrictions:]
+
+These pair styles can only be used if LAMMPS was built with the
+USER-CGDNA package and the MOLECULE and ASPHERE package.  See the "Making
+LAMMPS"_Section_start.html#start_3 section for more info on packages.
+
+[Related commands:]
+
+"bond_style oxdna_fene"_bond_oxdna_fene.html, "fix nve/dotc/langevin"_fix_nve_dotc_langevin.html, "pair_coeff"_pair_coeff.html 
+
+[Default:] none
+
+:line
+
+:link(Ouldridge-DPhil)
+[(Ouldrigde-DPhil)] T.E. Ouldridge, Coarse-grained modelling of DNA and DNA self-assembly, DPhil. University of Oxford (2011).
+
+:link(Ouldridge)
+[(Ouldridge)] T.E. Ouldridge, A.A. Louis, J.P.K. Doye, J. Chem. Phys. 134, 085101 (2011).
--- a/doc/src/pairs.txt
+++ b/doc/src/pairs.txt
@ -65,6 +65,7 @@ Pair Styles :h1
   pair_nb3b_harmonic
   pair_nm
   pair_none
+   pair_oxdna_excv
   pair_peri
   pair_polymorphic
   pair_quip
--- a/examples/USER/cgdna/README
+++ b/examples/USER/cgdna/README
@ -0,0 +1,28 @@
+This directory contains example data and input files 
+and utility scripts for the oxDNA coarse-grained model 
+for DNA.
+
+/examples/duplex1:
+Input, data and log files for a DNA duplex (double-stranded DNA) 
+consisiting of 5 base pairs. The duplex contains two strands with 
+complementary base pairs. The topology is
+
+A - A - A - A - A
+|   |   |   |   |
+T - T - T - T - T     
+
+/examples/duplex2:
+Input, data and log files for a nicked DNA duplex (double-stranded DNA) 
+consisiting of 8 base pairs. The duplex contains strands with 
+complementary base pairs, but the backbone on one side is not continuous: 
+two individual strands on one side form a duplex with a longer single 
+strand on the other side. The topology is
+
+A - A - A - A - A - A - A - A
+|   |   |   |   |   |   |   |
+T - T - T   T - T - T - T - T
+
+/util:
+This directory contains a simple python setup tool which creates 
+single straight or helical DNA strands, DNA duplexes or arrays of DNA 
+duplexes.
--- a/examples/USER/cgdna/examples/duplex1/data.duplex1
+++ b/examples/USER/cgdna/examples/duplex1/data.duplex1
@ -0,0 +1,74 @@
+# LAMMPS data file
+10 atoms
+10 ellipsoids
+8 bonds
+
+4 atom types
+1 bond types
+
+# System size
+-20.000000 20.000000 xlo xhi
+-20.000000 20.000000 ylo yhi
+-20.000000 20.000000 zlo zhi
+
+# Atom masses for each atom type
+Masses
+
+1 3.1575
+2 3.1575
+3 3.1575
+4 3.1575
+
+# Atom-ID, type, position, molecule-ID, ellipsoid flag, density
+Atoms
+
+1 1 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 1 1 1
+2 1 1.3274493266864451e-01 -4.2912827978022683e-01 3.7506163469402809e-01 1 1 1
+3 1 4.8460810659772807e-01 -7.0834970533509178e-01 7.5012326938805618e-01 1 1 1
+4 1 9.3267359196674593e-01 -7.4012419946742802e-01 1.1251849040820843e+00 1 1 1
+5 1 1.3204192238113461e+00 -5.1335201721887447e-01 1.5002465387761124e+00 1 1 1
+6 4 1.9958077618865377e-01 5.1335201721887447e-01 1.5002465387761124e+00 1 1 1
+7 4 5.8732640803325409e-01 7.4012419946742802e-01 1.1251849040820843e+00 1 1 1
+8 4 1.0353918934022719e+00 7.0834970533509178e-01 7.5012326938805618e-01 1 1 1
+9 4 1.3872550673313555e+00 4.2912827978022683e-01 3.7506163469402809e-01 1 1 1
+10 4 1.5200000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 1 1 1
+
+# Atom-ID, translational, rotational velocity
+Velocities
+
+1 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00
+2 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00
+3 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00
+4 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00
+5 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00
+6 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00
+7 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00
+8 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00
+9 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00
+10 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00
+
+# Atom-ID, shape, quaternion
+Ellipsoids
+
+1 1.1739845031423408e+00 1.1739845031423408e+00 1.1739845031423408e+00 1.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00
+2 1.1739845031423408e+00 1.1739845031423408e+00 1.1739845031423408e+00 9.5533648912560598e-01 0.0000000000000000e+00 0.0000000000000000e+00 2.9552020666133955e-01
+3 1.1739845031423408e+00 1.1739845031423408e+00 1.1739845031423408e+00 8.2533561490967822e-01 0.0000000000000000e+00 0.0000000000000000e+00 5.6464247339503526e-01
+4 1.1739845031423408e+00 1.1739845031423408e+00 1.1739845031423408e+00 6.2160996827066439e-01 0.0000000000000000e+00 0.0000000000000000e+00 7.8332690962748319e-01
+5 1.1739845031423408e+00 1.1739845031423408e+00 1.1739845031423408e+00 3.6235775447667351e-01 0.0000000000000000e+00 0.0000000000000000e+00 9.3203908596722607e-01
+6 1.1739845031423408e+00 1.1739845031423408e+00 1.1739845031423408e+00 0.0000000000000000e+00 9.3203908596722607e-01 -3.6235775447667351e-01 0.0000000000000000e+00
+7 1.1739845031423408e+00 1.1739845031423408e+00 1.1739845031423408e+00 0.0000000000000000e+00 7.8332690962748319e-01 -6.2160996827066439e-01 0.0000000000000000e+00
+8 1.1739845031423408e+00 1.1739845031423408e+00 1.1739845031423408e+00 0.0000000000000000e+00 5.6464247339503526e-01 -8.2533561490967822e-01 0.0000000000000000e+00
+9 1.1739845031423408e+00 1.1739845031423408e+00 1.1739845031423408e+00 0.0000000000000000e+00 2.9552020666133955e-01 -9.5533648912560598e-01 0.0000000000000000e+00
+10 1.1739845031423408e+00 1.1739845031423408e+00 1.1739845031423408e+00 0.0000000000000000e+00 0.0000000000000000e+00 -1.0000000000000000e+00 0.0000000000000000e+00
+
+# Bond topology
+Bonds
+
+1 1 1 2
+2 1 2 3
+3 1 3 4
+4 1 4 5
+5 1 6 7
+6 1 7 8
+7 1 8 9
+8 1 9 10
--- a/examples/USER/cgdna/examples/duplex1/input.duplex1
+++ b/examples/USER/cgdna/examples/duplex1/input.duplex1
@ -0,0 +1,75 @@
+variable number	equal 1
+variable ofreq	equal 1000
+variable efreq	equal 1000
+
+units lj
+
+dimension 3
+
+newton off
+
+boundary  p p p
+
+atom_style hybrid bond ellipsoid
+atom_modify sort 0 1.0
+
+# Pair interactions require lists of neighbours to be calculated
+neighbor 1.0 bin
+neigh_modify every 1 delay 0 check yes
+
+read_data data.duplex1
+
+set atom * mass 3.1575
+
+group all type 1 4
+
+# oxDNA bond interactions - FENE backbone
+bond_style oxdna_fene
+bond_coeff * 2.0 0.25 0.7525
+
+# oxDNA pair interactions
+pair_style hybrid/overlay oxdna_excv oxdna_stk oxdna_hbond oxdna_xstk oxdna_coaxstk
+pair_coeff * * oxdna_excv   2.0 0.7 0.675 2.0 0.515 0.5 2.0 0.33 0.32
+pair_coeff * * oxdna_stk    1.61048 6.0 0.4 0.9 0.32 0.6 1.3 0 0.8 0.9 0 0.95 0.9 0 0.95 2.0 0.65 2.0 0.65
+pair_coeff * * oxdna_hbond  0.0   8.0 0.4 0.75 0.34 0.7 1.5 0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45
+pair_coeff 1 4 oxdna_hbond  1.077 8.0 0.4 0.75 0.34 0.7 1.5 0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45
+pair_coeff 2 3 oxdna_hbond  1.077 8.0 0.4 0.75 0.34 0.7 1.5 0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45
+pair_coeff * * oxdna_xstk   47.5 0.575 0.675 0.495 0.655 2.25 0.791592653589793 0.58 1.7 1.0 0.68 1.7 1.0 0.68 1.5 0 0.65 1.7 0.875 0.68 1.7 0.875 0.68
+pair_coeff * * oxdna_coaxstk 46.0 0.4 0.6 0.22 0.58 2.0 2.541592653589793 0.65 1.3 0 0.8 0.9 0 0.95 0.9 0 0.95 2.0 -0.65 2.0 -0.65
+
+# NVE ensemble
+#fix 1 all   nve/dotc/langevin 0.1 0.1 0.03 457145 angmom 10
+fix 1 all   nve/dot
+
+timestep 1e-5
+
+#comm_style tiled
+#fix 3 all balance 10000 1.1 rcb
+
+#compute mol all chunk/atom molecule
+#compute mychunk all vcm/chunk mol
+#fix 4 all ave/time 10000 1 10000 c_mychunk[1] c_mychunk[2] c_mychunk[3] file vcm.txt mode vector
+
+dump pos all xyz ${ofreq} traj.${number}.xyz
+
+compute quat all property/atom quatw quati quatj quatk
+dump quat all custom ${ofreq} quat.${number}.txt id c_quat[1] c_quat[2] c_quat[3] c_quat[4]
+dump_modify quat sort id
+dump_modify quat format line "%d  %13.6le  %13.6le  %13.6le  %13.6le"
+
+compute erot all erotate/asphere
+compute ekin all ke
+compute epot all pe
+variable erot equal c_erot
+variable ekin equal c_ekin
+variable epot equal c_epot
+variable etot equal c_erot+c_ekin+c_epot
+fix 5 all print ${efreq} "$(step)  ekin = ${ekin} |  erot = ${erot} | epot = ${epot} | etot = ${etot}" screen yes
+
+dump out all custom ${ofreq} out.${number}.txt id x y z vx vy vz fx fy fz tqx tqy tqz
+dump_modify out sort id
+dump_modify out format line "%d   %13.6le %13.6le %13.6le  %13.6le %13.6le %13.6le  %13.6le %13.6le %13.6le  %13.6le %13.6le %13.6le"
+
+run 1000000
+
+#write_restart config.${number}.*
--- a/examples/USER/cgdna/examples/duplex1/log.9Jan17.duplex1.g++.1
+++ b/examples/USER/cgdna/examples/duplex1/log.9Jan17.duplex1.g++.1
--- a/examples/USER/cgdna/examples/duplex1/log.9Jan17.duplex1.g++.4
+++ b/examples/USER/cgdna/examples/duplex1/log.9Jan17.duplex1.g++.4
--- a/examples/USER/cgdna/examples/duplex2/data.duplex2
+++ b/examples/USER/cgdna/examples/duplex2/data.duplex2
@ -0,0 +1,97 @@
+# LAMMPS data file
+16 atoms
+16 ellipsoids
+13 bonds
+
+4 atom types
+1 bond types
+
+# System size
+-20.0 20.0 xlo xhi
+-20.0 20.0 ylo yhi
+-20.0 20.0 zlo zhi
+
+# Atom masses for each atom type
+Masses
+
+1 3.1575
+2 3.1575
+3 3.1575
+4 3.1575
+
+# Atom-ID, type, position, molecule-ID, ellipsoid flag, density
+Atoms
+
+1  1  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  1 1 1
+2  1  1.327449326686445e-01 -4.291282797802268e-01  3.750616346940281e-01  1 1 1
+3  1  4.846081065977281e-01 -7.083497053350921e-01  7.501232693880562e-01  1 1 1
+4  1  9.326735919667459e-01 -7.401241994674285e-01  1.125184904082084e+00  1 1 1
+5  1  1.320419223811347e+00 -5.133520172188747e-01  1.500246538776112e+00  1 1 1
+6  1  1.512394297416339e+00 -1.072512061254991e-01  1.875308173470140e+00  1 1 1
+7  1  1.441536396413952e+00  3.363155369040876e-01  2.250369808164169e+00  1 1 1
+8  1  1.132598224218932e+00  6.623975870343269e-01  2.625431442858197e+00  1 1 1
+9  4  5.873264080332541e-01  7.401241994674285e-01  1.125184904082084e+00  1 1 1
+10 4  1.035391893402272e+00  7.083497053350921e-01  7.501232693880562e-01  1 1 1
+11 4  1.387255067331356e+00  4.291282797802267e-01  3.750616346940281e-01  1 1 1
+12 4  1.520000000000000e+00  1.260981291332700e-33  0.000000000000000e+00  1 1 1
+13 4  3.874017757810680e-01 -6.623975870343268e-01  2.625431442858197e+00  1 1 1
+14 4  7.846360358604798e-02 -3.363155369040874e-01  2.250369808164169e+00  1 1 1
+15 4  7.605702583661333e-03  1.072512061254995e-01  1.875308173470140e+00  1 1 1
+16 4  1.995807761886533e-01  5.133520172188748e-01  1.500246538776112e+00  1 1 1
+
+# Atom-ID, translational, rotational velocity
+Velocities
+
+1  0.0  0.0  0.0  0.0  0.0  0.0 
+2  0.0  0.0  0.0  0.0  0.0  0.0 
+3  0.0  0.0  0.0  0.0  0.0  0.0 
+4  0.0  0.0  0.0  0.0  0.0  0.0 
+5  0.0  0.0  0.0  0.0  0.0  0.0 
+6  0.0  0.0  0.0  0.0  0.0  0.0 
+7  0.0  0.0  0.0  0.0  0.0  0.0 
+8  0.0  0.0  0.0  0.0  0.0  0.0 
+9  0.0  0.0  0.0  0.0  0.0  0.0 
+10 0.0  0.0  0.0  0.0  0.0  0.0 
+11 0.0  0.0  0.0  0.0  0.0  0.0 
+12 0.0  0.0  0.0  0.0  0.0  0.0 
+13 0.0  0.0  0.0  0.0  0.0  0.0 
+14 0.0  0.0  0.0  0.0  0.0  0.0 
+15 0.0  0.0  0.0  0.0  0.0  0.0 
+16 0.0  0.0  0.0  0.0  0.0  0.0 
+
+# Atom-ID, shape, quaternion
+Ellipsoids
+
+1   1.1739845031423408 1.1739845031423408 1.1739845031423408  1.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00
+2   1.1739845031423408 1.1739845031423408 1.1739845031423408  9.553364891256060e-01  0.000000000000000e+00  0.000000000000000e+00  2.955202066613395e-01 
+3   1.1739845031423408 1.1739845031423408 1.1739845031423408  8.253356149096783e-01  0.000000000000000e+00  0.000000000000000e+00  5.646424733950354e-01 
+4   1.1739845031423408 1.1739845031423408 1.1739845031423408  6.216099682706646e-01  0.000000000000000e+00  0.000000000000000e+00  7.833269096274833e-01 
+5   1.1739845031423408 1.1739845031423408 1.1739845031423408  3.623577544766736e-01  0.000000000000000e+00  0.000000000000000e+00  9.320390859672263e-01 
+6   1.1739845031423408 1.1739845031423408 1.1739845031423408  7.073720166770291e-02  0.000000000000000e+00  0.000000000000000e+00  9.974949866040544e-01 
+7   1.1739845031423408 1.1739845031423408 1.1739845031423408 -2.272020946930869e-01 -0.000000000000000e+00  0.000000000000000e+00  9.738476308781953e-01 
+8   1.1739845031423408 1.1739845031423408 1.1739845031423408 -5.048461045998575e-01 -0.000000000000000e+00  0.000000000000000e+00  8.632093666488738e-01 
+9   1.1739845031423408 1.1739845031423408 1.1739845031423408  4.796493962806427e-17  7.833269096274833e-01 -6.216099682706646e-01  3.806263289803786e-17 
+10  1.1739845031423408 1.1739845031423408 1.1739845031423408  5.707093416549944e-17  5.646424733950354e-01 -8.253356149096784e-01  2.218801320830406e-17 
+11  1.1739845031423408 1.1739845031423408 1.1739845031423408  6.107895212550935e-17  2.955202066613394e-01 -9.553364891256061e-01  4.331404380149668e-18 
+12  1.1739845031423408 1.1739845031423408 1.1739845031423408  5.963096920061075e-17  0.000000000000000e+00 -1.000000000000000e+00 -1.391211590127312e-17 
+13  1.1739845031423408 1.1739845031423408 1.1739845031423408  5.285632939302787e-17  8.632093666488739e-01  5.048461045998572e-01 -3.091290830301125e-17 
+14  1.1739845031423408 1.1739845031423408 1.1739845031423408  4.136019110019290e-17  9.738476308781953e-01  2.272020946930868e-01 -4.515234267244800e-17 
+15  1.1739845031423408 1.1739845031423408 1.1739845031423408  2.616947011741696e-17  9.974949866040544e-01 -7.073720166770313e-02 -5.535845274597425e-17 
+16  1.1739845031423408 1.1739845031423408 1.1739845031423408  8.641108308308281e-18  9.320390859672264e-01 -3.623577544766736e-01 -6.061955710708163e-17 
+
+# Bond-ID, type, atom pairs
+Bonds
+
+1	1	1	2
+2	1	2	3
+3	1	3	4
+4	1	4	5
+5	1	5	6
+6	1	6	7
+7	1	7	8
+8	1	13	14
+9	1	14	15
+10	1	15	16
+11	1	9	10
+12	1	10	11
+13	1	11	12
--- a/examples/USER/cgdna/examples/duplex2/input.duplex2
+++ b/examples/USER/cgdna/examples/duplex2/input.duplex2
@ -0,0 +1,75 @@
+variable number	equal 2
+variable ofreq	equal 1000
+variable efreq	equal 1000
+
+units lj
+
+dimension 3
+
+newton off
+
+boundary  p p p
+
+atom_style hybrid bond ellipsoid
+atom_modify sort 0 1.0
+
+# Pair interactions require lists of neighbours to be calculated
+neighbor 1.0 bin
+neigh_modify every 1 delay 0 check yes
+
+read_data data.duplex2
+
+set atom * mass 3.1575
+
+group all type 1 4
+
+# oxDNA bond interactions - FENE backbone
+bond_style oxdna_fene
+bond_coeff * 2.0 0.25 0.7525
+
+# oxDNA pair interactions
+pair_style hybrid/overlay oxdna_excv oxdna_stk oxdna_hbond oxdna_xstk oxdna_coaxstk
+pair_coeff * * oxdna_excv   2.0 0.7 0.675 2.0 0.515 0.5 2.0 0.33 0.32
+pair_coeff * * oxdna_stk    1.61048 6.0 0.4 0.9 0.32 0.6 1.3 0 0.8 0.9 0 0.95 0.9 0 0.95 2.0 0.65 2.0 0.65   
+pair_coeff * * oxdna_hbond  0.0   8.0 0.4 0.75 0.34 0.7 1.5 0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45
+pair_coeff 1 4 oxdna_hbond  1.077 8.0 0.4 0.75 0.34 0.7 1.5 0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45
+pair_coeff 2 3 oxdna_hbond  1.077 8.0 0.4 0.75 0.34 0.7 1.5 0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45
+pair_coeff * * oxdna_xstk   47.5 0.575 0.675 0.495 0.655 2.25 0.791592653589793 0.58 1.7 1.0 0.68 1.7 1.0 0.68 1.5 0 0.65 1.7 0.875 0.68 1.7 0.875 0.68 
+pair_coeff * * oxdna_coaxstk 46.0 0.4 0.6 0.22 0.58 2.0 2.541592653589793 0.65 1.3 0 0.8 0.9 0 0.95 0.9 0 0.95 2.0 -0.65 2.0 -0.65
+
+# NVE ensemble
+fix 1 all   nve/dotc/langevin 0.1 0.1 0.03 457145 angmom 10
+#fix 1 all   nve/dot
+
+timestep 1e-5 
+
+#comm_style tiled
+#fix 3 all balance 10000 1.1 rcb
+
+#compute mol all chunk/atom molecule
+#compute mychunk all vcm/chunk mol
+#fix 4 all ave/time 10000 1 10000 c_mychunk[1] c_mychunk[2] c_mychunk[3] file vcm.txt mode vector
+
+dump pos all xyz ${ofreq} traj.${number}.xyz
+
+compute quat all property/atom quatw quati quatj quatk
+dump quat all custom ${ofreq} quat.${number}.txt id c_quat[1] c_quat[2] c_quat[3] c_quat[4]
+dump_modify quat sort id
+dump_modify quat format line "%d  %13.6le  %13.6le  %13.6le  %13.6le"
+
+compute erot all erotate/asphere
+compute ekin all ke
+compute epot all pe
+variable erot equal c_erot
+variable ekin equal c_ekin
+variable epot equal c_epot
+variable etot equal c_erot+c_ekin+c_epot
+fix 5 all print ${efreq} "$(step)  ekin = ${ekin} |  erot = ${erot} | epot = ${epot} | etot = ${etot}" screen yes
+
+dump out all custom ${ofreq} out.${number}.txt id x y z vx vy vz fx fy fz tqx tqy tqz
+dump_modify out sort id
+dump_modify out format line "%d   %13.6le %13.6le %13.6le  %13.6le %13.6le %13.6le  %13.6le %13.6le %13.6le  %13.6le %13.6le %13.6le"
+
+run 1000000
+
+#write_restart config.${number}.*
--- a/examples/USER/cgdna/examples/duplex2/log.9Jan17.duplex2.g++.1
+++ b/examples/USER/cgdna/examples/duplex2/log.9Jan17.duplex2.g++.1
--- a/examples/USER/cgdna/examples/duplex2/log.9Jan17.duplex2.g++.4
+++ b/examples/USER/cgdna/examples/duplex2/log.9Jan17.duplex2.g++.4
--- a/examples/USER/cgdna/util/generate_input.py
+++ b/examples/USER/cgdna/util/generate_input.py
@ -0,0 +1,388 @@
+# Setup tool for oxDNA input in LAMMPS format.
+
+import math,numpy as np,sys,os
+
+# system size
+lxmin = -115.0
+lxmax = +115.0
+lymin = -115.0
+lymax = +115.0
+lzmin = -115.0
+lzmax = +115.0
+
+# rise in z-direction
+r0 = 0.7
+
+# definition of single untwisted strand
+def single():
+
+  strand = inp[1].split(':')
+
+  com_start=strand[0].split(',')
+
+  posx=float(com_start[0])
+  posy=float(com_start[1])
+  posz=float(com_start[2])
+  risex=0
+  risey=0
+  risez=r0
+
+  strandstart=len(nucleotide)+1
+
+  for letter in strand[2]:
+    temp=[]
+
+    temp.append(nt2num[letter])
+    temp.append([posx,posy,posz])
+    vel=[0,0,0,0,0,0]
+    temp.append(vel)
+    temp.append(shape)
+
+    quat=[1,0,0,0]
+    temp.append(quat)
+
+    posx=posx+risex
+    posy=posy+risey
+    posz=posz+risez
+
+    if (len(nucleotide)+1 > strandstart):
+      topology.append([1,len(nucleotide),len(nucleotide)+1])
+
+    nucleotide.append(temp)
+
+  return
+
+# definition of single twisted strand
+def single_helix():
+
+  strand = inp[1].split(':')
+
+  com_start=strand[0].split(',')
+  twist=float(strand[1])
+
+  posx = float(com_start[0])
+  posy = float(com_start[1])
+  posz = float(com_start[2])
+  risex=0
+  risey=0
+  risez=math.sqrt(r0**2-4.0*math.sin(0.5*twist)**2) 
+
+  dcomh=0.76
+  axisx=dcomh + posx
+  axisy=posy
+
+  strandstart=len(nucleotide)+1
+  quat=[1,0,0,0]
+
+  qrot0=math.cos(0.5*twist)
+  qrot1=0
+  qrot2=0
+  qrot3=math.sin(0.5*twist)
+
+  for letter in strand[2]:
+    temp=[]
+
+    temp.append(nt2num[letter])
+    temp.append([posx,posy,posz])
+    vel=[0,0,0,0,0,0]
+    temp.append(vel)
+    temp.append(shape)
+
+    temp.append(quat)
+
+    quat0 = quat[0]*qrot0 - quat[1]*qrot1 - quat[2]*qrot2 - quat[3]*qrot3 
+    quat1 = quat[0]*qrot1 + quat[1]*qrot0 + quat[2]*qrot3 - quat[3]*qrot2 
+    quat2 = quat[0]*qrot2 + quat[2]*qrot0 + quat[3]*qrot1 - quat[1]*qrot3 
+    quat3 = quat[0]*qrot3 + quat[3]*qrot0 + quat[1]*qrot2 + quat[2]*qrot1 
+
+    quat = [quat0,quat1,quat2,quat3]
+
+    posx=axisx - dcomh*(quat[0]**2+quat[1]**2-quat[2]**2-quat[3]**2)
+    posy=axisy - dcomh*(2*(quat[1]*quat[2]+quat[0]*quat[3]))
+    posz=posz+risez
+
+    if (len(nucleotide)+1 > strandstart):
+      topology.append([1,len(nucleotide),len(nucleotide)+1])
+
+    nucleotide.append(temp)
+
+  return
+
+# definition of twisted duplex  
+def duplex():
+
+  strand = inp[1].split(':')
+
+  com_start=strand[0].split(',')
+  twist=float(strand[1])
+
+  compstrand=[]
+  comptopo=[]
+
+  posx1 = float(com_start[0])
+  posy1 = float(com_start[1])
+  posz1 = float(com_start[2])
+
+  risex=0
+  risey=0
+  risez=math.sqrt(r0**2-4.0*math.sin(0.5*twist)**2) 
+
+  dcomh=0.76
+  axisx=dcomh + posx1
+  axisy=posy1
+
+  posx2 = axisx + dcomh  
+  posy2 = posy1
+  posz2 = posz1
+
+  strandstart=len(nucleotide)+1
+
+  quat1=[1,0,0,0]
+  quat2=[0,0,-1,0]
+
+  qrot0=math.cos(0.5*twist)
+  qrot1=0
+  qrot2=0
+  qrot3=math.sin(0.5*twist)
+
+  for letter in strand[2]:
+    temp1=[]
+    temp2=[]
+
+    temp1.append(nt2num[letter])
+    temp2.append(compnt2num[letter])
+
+    temp1.append([posx1,posy1,posz1])
+    temp2.append([posx2,posy2,posz2])
+
+    vel=[0,0,0,0,0,0]
+    temp1.append(vel)
+    temp2.append(vel)
+
+    temp1.append(shape)
+    temp2.append(shape)
+
+    temp1.append(quat1)
+    temp2.append(quat2)
+
+    quat1_0 = quat1[0]*qrot0 - quat1[1]*qrot1 - quat1[2]*qrot2 - quat1[3]*qrot3 
+    quat1_1 = quat1[0]*qrot1 + quat1[1]*qrot0 + quat1[2]*qrot3 - quat1[3]*qrot2 
+    quat1_2 = quat1[0]*qrot2 + quat1[2]*qrot0 + quat1[3]*qrot1 - quat1[1]*qrot3 
+    quat1_3 = quat1[0]*qrot3 + quat1[3]*qrot0 + quat1[1]*qrot2 + quat1[2]*qrot1 
+
+    quat1 = [quat1_0,quat1_1,quat1_2,quat1_3]
+
+    posx1=axisx - dcomh*(quat1[0]**2+quat1[1]**2-quat1[2]**2-quat1[3]**2)
+    posy1=axisy - dcomh*(2*(quat1[1]*quat1[2]+quat1[0]*quat1[3]))
+    posz1=posz1+risez
+
+    quat2_0 = quat2[0]*qrot0 - quat2[1]*qrot1 - quat2[2]*qrot2 + quat2[3]*qrot3 
+    quat2_1 = quat2[0]*qrot1 + quat2[1]*qrot0 - quat2[2]*qrot3 - quat2[3]*qrot2 
+    quat2_2 = quat2[0]*qrot2 + quat2[2]*qrot0 + quat2[3]*qrot1 + quat2[1]*qrot3 
+    quat2_3 =-quat2[0]*qrot3 + quat2[3]*qrot0 + quat2[1]*qrot2 + quat2[2]*qrot1 
+
+    quat2 = [quat2_0,quat2_1,quat2_2,quat2_3]
+
+    posx2=axisx + dcomh*(quat1[0]**2+quat1[1]**2-quat1[2]**2-quat1[3]**2)
+    posy2=axisy + dcomh*(2*(quat1[1]*quat1[2]+quat1[0]*quat1[3]))
+    posz2=posz1
+
+    if (len(nucleotide)+1 > strandstart):
+      topology.append([1,len(nucleotide),len(nucleotide)+1])
+      comptopo.append([1,len(nucleotide)+len(strand[2]),len(nucleotide)+len(strand[2])+1])
+
+    nucleotide.append(temp1)
+    compstrand.append(temp2)
+
+  for ib in range(len(compstrand)):
+    nucleotide.append(compstrand[len(compstrand)-1-ib])
+
+  for ib in range(len(comptopo)):
+    topology.append(comptopo[ib])
+
+  return
+
+# definition of array of duplexes  
+def duplex_array():
+
+  strand = inp[1].split(':')
+  number=strand[0].split(',')
+  posz1_0 = float(strand[1])
+  twist=float(strand[2])
+
+  nx = int(number[0])
+  ny = int(number[1])
+
+  dx = (lxmax-lxmin)/nx
+  dy = (lymax-lymin)/ny
+
+  risex=0
+  risey=0
+  risez=math.sqrt(r0**2-4.0*math.sin(0.5*twist)**2) 
+  dcomh=0.76
+
+  for ix in range(nx):
+
+    axisx=lxmin + dx/2 + ix * dx
+
+    for iy in range(ny):
+
+      axisy=lymin + dy/2 + iy * dy
+
+      compstrand=[]
+      comptopo=[]
+
+      posx1 = axisx - dcomh
+      posy1 = axisy
+      posz1 = posz1_0
+
+      posx2 = axisx + dcomh  
+      posy2 = posy1
+      posz2 = posz1
+
+      strandstart=len(nucleotide)+1
+      quat1=[1,0,0,0]
+      quat2=[0,0,-1,0]
+
+      qrot0=math.cos(0.5*twist)
+      qrot1=0
+      qrot2=0
+      qrot3=math.sin(0.5*twist)
+
+      for letter in strand[3]:
+	temp1=[]
+	temp2=[]
+
+	temp1.append(nt2num[letter])
+	temp2.append(compnt2num[letter])
+
+	temp1.append([posx1,posy1,posz1])
+	temp2.append([posx2,posy2,posz2])
+
+	vel=[0,0,0,0,0,0]
+	temp1.append(vel)
+	temp2.append(vel)
+
+	temp1.append(shape)
+	temp2.append(shape)
+
+	temp1.append(quat1)
+	temp2.append(quat2)
+
+	quat1_0 = quat1[0]*qrot0 - quat1[1]*qrot1 - quat1[2]*qrot2 - quat1[3]*qrot3 
+	quat1_1 = quat1[0]*qrot1 + quat1[1]*qrot0 + quat1[2]*qrot3 - quat1[3]*qrot2 
+	quat1_2 = quat1[0]*qrot2 + quat1[2]*qrot0 + quat1[3]*qrot1 - quat1[1]*qrot3 
+	quat1_3 = quat1[0]*qrot3 + quat1[3]*qrot0 + quat1[1]*qrot2 + quat1[2]*qrot1 
+
+	quat1 = [quat1_0,quat1_1,quat1_2,quat1_3]
+
+	posx1=axisx - dcomh*(quat1[0]**2+quat1[1]**2-quat1[2]**2-quat1[3]**2)
+	posy1=axisy - dcomh*(2*(quat1[1]*quat1[2]+quat1[0]*quat1[3]))
+	posz1=posz1+risez
+
+	quat2_0 = quat2[0]*qrot0 - quat2[1]*qrot1 - quat2[2]*qrot2 + quat2[3]*qrot3 
+	quat2_1 = quat2[0]*qrot1 + quat2[1]*qrot0 - quat2[2]*qrot3 - quat2[3]*qrot2 
+	quat2_2 = quat2[0]*qrot2 + quat2[2]*qrot0 + quat2[3]*qrot1 + quat2[1]*qrot3 
+	quat2_3 =-quat2[0]*qrot3 + quat2[3]*qrot0 + quat2[1]*qrot2 + quat2[2]*qrot1 
+
+	quat2 = [quat2_0,quat2_1,quat2_2,quat2_3]
+
+	posx2=axisx + dcomh*(quat1[0]**2+quat1[1]**2-quat1[2]**2-quat1[3]**2)
+	posy2=axisy + dcomh*(2*(quat1[1]*quat1[2]+quat1[0]*quat1[3]))
+	posz2=posz1
+
+	if (len(nucleotide)+1 > strandstart):
+	  topology.append([1,len(nucleotide),len(nucleotide)+1])
+	  comptopo.append([1,len(nucleotide)+len(strand[3]),len(nucleotide)+len(strand[3])+1])
+
+	nucleotide.append(temp1)
+	compstrand.append(temp2)
+
+      for ib in range(len(compstrand)):
+	nucleotide.append(compstrand[len(compstrand)-1-ib])
+
+      for ib in range(len(comptopo)):
+	topology.append(comptopo[ib])
+
+  return
+
+# main part
+nt2num = {'A':1, 'C':2, 'G':3, 'T':4}
+compnt2num = {'T':1, 'G':2, 'C':3, 'A':4}
+shape = [1.1739845031423408,1.1739845031423408,1.1739845031423408]
+
+nucleotide=[]
+topology=[]
+
+seqfile = open(sys.argv[1],'r')
+
+# process sequence file line by line
+for line in seqfile:
+
+  inp = line.split()
+  if inp[0] == 'single':
+    single()
+  if inp[0] == 'single_helix':
+    single_helix()
+  if inp[0] == 'duplex':
+    duplex()
+  if inp[0] == 'duplex_array':
+    duplex_array()
+
+# output atom data in LAMMPS format
+out = open(sys.argv[2],'w')
+
+out.write('# LAMMPS data file\n')
+out.write('%d atoms\n' % len(nucleotide))
+out.write('%d ellipsoids\n' % len(nucleotide))
+out.write('%d bonds\n' % len(topology))
+out.write('\n')
+out.write('4 atom types\n')
+out.write('1 bond types\n')
+out.write('\n')
+out.write('# System size\n')
+out.write('%f %f xlo xhi\n' % (lxmin,lxmax))
+out.write('%f %f ylo yhi\n' % (lymin,lymax))
+out.write('%f %f zlo zhi\n' % (lzmin,lzmax))
+out.write('\n')
+out.write('Masses\n')
+out.write('\n')
+out.write('1 3.1575\n')
+out.write('2 3.1575\n')
+out.write('3 3.1575\n')
+out.write('4 3.1575\n')
+
+out.write('\n')
+out.write('# Atom-ID, type, position, molecule-ID, ellipsoid flag, density\n')
+out.write('Atoms\n')
+out.write('\n')
+for ib in range(len(nucleotide)):
+  out.write("%d %d %22.16le %22.16le %22.16le 1 1 1\n" % (ib+1,nucleotide[ib][0],nucleotide[ib][1][0],nucleotide[ib][1][1],nucleotide[ib][1][2]))
+
+out.write('\n')
+out.write('# Atom-ID, translational, rotational velocity\n')
+out.write('Velocities\n')
+out.write('\n')
+for ib in range(len(nucleotide)):
+  out.write("%d %22.16le %22.16le %22.16le %22.16le %22.16le %22.16le\n" % (ib+1,nucleotide[ib][2][0],nucleotide[ib][2][1],nucleotide[ib][2][2],nucleotide[ib][2][3],nucleotide[ib][2][4],nucleotide[ib][2][5]))
+
+out.write('\n')
+out.write('# Atom-ID, shape, quaternion\n')
+out.write('Ellipsoids\n')
+out.write('\n')
+for ib in range(len(nucleotide)):
+  out.write("%d %22.16le %22.16le %22.16le %22.16le %22.16le %22.16le %22.16le\n" % (ib+1,nucleotide[ib][3][0],nucleotide[ib][3][1],nucleotide[ib][3][2],nucleotide[ib][4][0],nucleotide[ib][4][1],nucleotide[ib][4][2],nucleotide[ib][4][3]))
+
+out.write('\n')
+out.write('# Bond topology\n')
+out.write('Bonds\n')
+out.write('\n')
+for ib in range(len(topology)):
+  out.write("%d %d %d %d\n" % (ib+1,topology[ib][0],topology[ib][1],topology[ib][2]))
+
+out.close() 
+
+seqfile.close()
+sys.exit(0)
+
+
--- a/examples/USER/cgdna/util/input.ref
+++ b/examples/USER/cgdna/util/input.ref
@ -0,0 +1,77 @@
+variable number	equal 8
+variable ofreq	equal 1000
+variable efreq	equal 1000
+
+units lj
+
+dimension 3
+
+newton off
+
+processors 1 1 1
+
+boundary  p p p
+
+atom_style hybrid bond ellipsoid
+atom_modify sort 0 1.0
+
+# Pair interactions require lists of neighbours to be calculated
+neighbor 1.0 bin
+neigh_modify every 1 delay 0 check yes
+
+read_data data.duplex2
+
+set atom * mass 3.1575
+
+group all type 1 4
+
+# oxDNA bond interactions - FENE backbone
+bond_style oxdna_fene
+bond_coeff * 2.0 0.25 0.7525
+
+# oxDNA pair interactions
+pair_style hybrid/overlay oxdna_excv oxdna_stk oxdna_hbond oxdna_xstk oxdna_coaxstk
+pair_coeff * * oxdna_excv   2.0 0.7 0.675 2.0 0.515 0.5 2.0 0.33 0.32
+pair_coeff * * oxdna_stk    1.61048 6.0 0.4 0.9 0.32 0.6 1.3 0 0.8 0.9 0 0.95 0.9 0 0.95 2.0 0.65 2.0 0.65   
+pair_coeff * * oxdna_hbond  0.0   8.0 0.4 0.75 0.34 0.7 1.5 0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45
+pair_coeff 1 4 oxdna_hbond  1.077 8.0 0.4 0.75 0.34 0.7 1.5 0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45
+pair_coeff 2 3 oxdna_hbond  1.077 8.0 0.4 0.75 0.34 0.7 1.5 0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45
+pair_coeff * * oxdna_xstk   47.5 0.575 0.675 0.495 0.655 2.25 0.791592653589793 0.58 1.7 1.0 0.68 1.7 1.0 0.68 1.5 0 0.65 1.7 0.875 0.68 1.7 0.875 0.68 
+pair_coeff * * oxdna_coaxstk 46.0 0.4 0.6 0.22 0.58 2.0 2.541592653589793 0.65 1.3 0 0.8 0.9 0 0.95 0.9 0 0.95 2.0 -0.65 2.0 -0.65
+
+# NVE ensemble
+#fix 1 all   nve/dotc/langevin 0.1 0.1 0.03 457145 angmom 10
+fix 1 all   nve/dot
+
+timestep 1e-5 
+
+#comm_style tiled
+#fix 3 all balance 10000 1.1 rcb
+
+#compute mol all chunk/atom molecule
+#compute mychunk all vcm/chunk mol
+#fix 4 all ave/time 10000 1 10000 c_mychunk[1] c_mychunk[2] c_mychunk[3] file vcm.txt mode vector
+
+#dump pos all xyz ${ofreq} traj.${number}.xyz
+
+#compute quat all property/atom quatw quati quatj quatk
+#dump quat all custom ${ofreq} quat.${number}.txt id c_quat[1] c_quat[2] c_quat[3] c_quat[4]
+#dump_modify quat sort id
+#dump_modify quat format line "%d  %13.6le  %13.6le  %13.6le  %13.6le"
+
+compute erot all erotate/asphere
+compute ekin all ke
+compute epot all pe
+variable erot equal c_erot
+variable ekin equal c_ekin
+variable epot equal c_epot
+variable etot equal c_erot+c_ekin+c_epot
+fix 5 all print ${efreq} "$(step)  ekin = ${ekin} |  erot = ${erot} | epot = ${epot} | etot = ${etot}" screen yes
+
+dump out all custom ${ofreq} out.${number}.txt id x y z vx vy vz fx fy fz tqx tqy tqz
+dump_modify out sort id
+dump_modify out format line "%d   %13.6le %13.6le %13.6le  %13.6le %13.6le %13.6le  %13.6le %13.6le %13.6le  %13.6le %13.6le %13.6le"
+
+run 1000000
+
+#write_restart config.${number}.*
--- a/examples/USER/cgdna/util/sequence.txt
+++ b/examples/USER/cgdna/util/sequence.txt
@ -0,0 +1,4 @@
+single 0,0,0:0.6:AAAAA
+single_helix 0,0,0:0.6:AAAAA
+duplex 0,0,0:0.6:AAAAA
+duplex_array 10,10:-112.0:0.6:AAAAA
--- a/examples/VISCOSITY/in.einstein.2d
+++ b/examples/VISCOSITY/in.einstein.2d
@ -54,7 +54,8 @@ reset_timestep  0
 variable         pxy equal pxy
 variable         pxx equal pxx-press

-fix              avstress all ave/time $s $p $d v_pxy v_pxx ave one file einstein.dat
+fix              avstress all ave/time $s $p $d v_pxy v_pxx ave one &
+                 file profile.einstein.2d

 # Diagonal components of SS are larger by factor 2-2/d,
 # which is 4/3 for d=3, but 1 for d=2.
--- a/examples/coreshell/log.9Nov16.coreshell.dsf.g++.1
+++ b/examples/coreshell/log.9Nov16.coreshell.dsf.g++.1
--- a/examples/coreshell/log.9Nov16.coreshell.dsf.g++.4
+++ b/examples/coreshell/log.9Nov16.coreshell.dsf.g++.4
@ -0,0 +1,189 @@
+LAMMPS (26 Jan 2017)
+# Testsystem for core-shell model compared to Mitchel and Finchham
+# Hendrik Heenen, June 2014
+
+# ------------------------ INITIALIZATION ----------------------------
+
+units 		metal
+dimension	3
+boundary	p	p	p
+atom_style	full
+
+# ----------------------- ATOM DEFINITION ----------------------------
+
+fix csinfo all property/atom i_CSID
+read_data data.coreshell fix csinfo NULL CS-Info
+  orthogonal box = (0 0 0) to (24.096 24.096 24.096)
+  1 by 2 by 2 MPI processor grid
+  reading atoms ...
+  432 atoms
+  scanning bonds ...
+  1 = max bonds/atom
+  reading bonds ...
+  216 bonds
+  1 = max # of 1-2 neighbors
+  0 = max # of 1-3 neighbors
+  0 = max # of 1-4 neighbors
+  1 = max # of special neighbors
+
+group cores type 1 2
+216 atoms in group cores
+group shells type 3 4
+216 atoms in group shells
+
+neighbor 2.0 bin
+comm_modify vel yes
+
+# ------------------------ FORCE FIELDS ------------------------------
+
+pair_style   born/coul/dsf/cs 0.1 20.0 20.0    # A, rho, sigma=0, C, D
+pair_coeff   * *      0.0 1.000   0.00  0.00   0.00
+pair_coeff   3 3    487.0 0.23768 0.00  1.05   0.50 #Na-Na
+pair_coeff   3 4 145134.0 0.23768 0.00  6.99   8.70 #Na-Cl
+pair_coeff   4 4 405774.0 0.23768 0.00 72.40 145.40 #Cl-Cl
+
+bond_style harmonic
+bond_coeff 1 63.014 0.0
+bond_coeff 2 25.724 0.0
+
+# ------------------------ Equilibration Run -------------------------------
+
+reset_timestep 0
+
+thermo 50
+thermo_style custom step etotal pe ke temp press              epair evdwl ecoul elong ebond fnorm fmax vol
+
+compute CSequ all temp/cs cores shells
+
+# output via chunk method
+
+#compute prop all property/atom i_CSID
+#compute cs_chunk all chunk/atom c_prop
+#compute cstherm all temp/chunk cs_chunk temp internal com yes cdof 3.0
+#fix ave_chunk all ave/time 100 1 100 c_cstherm file chunk.dump mode vector
+
+thermo_modify temp CSequ
+
+# velocity bias option
+
+velocity all create 1427 134 dist gaussian mom yes rot no bias yes temp CSequ
+Neighbor list info ...
+  update every 1 steps, delay 10 steps, check yes
+  max neighbors/atom: 2000, page size: 100000
+  master list distance cutoff = 22
+  ghost atom cutoff = 22
+  binsize = 11, bins = 3 3 3
+  1 neighbor lists, perpetual/occasional/extra = 1 0 0
+  (1) pair born/coul/dsf/cs, half, perpetual
+      pair build: half/bin/newton
+      stencil: half/bin/3d/newton
+      bin: standard
+velocity all scale 1427 temp CSequ
+
+fix thermoberendsen all temp/berendsen 1427 1427 0.4
+fix nve all nve
+fix_modify thermoberendsen temp CSequ
+
+# 2 fmsec timestep
+
+timestep 0.002
+run 500
+Memory usage per processor = 6.8559 Mbytes
+Step TotEng PotEng KinEng Temp Press E_pair E_vdwl E_coul E_long E_bond Fnorm Fmax Volume 
+       0   -635.80596   -675.46362    39.657659         1427   -21302.622   -675.46362    1.6320365   -677.09565            0            0 1.5814015e-14 3.2317898e-15      13990.5 
+      50   -634.07021   -666.11867    32.048452    1153.1982    -4560.945   -668.28236    37.756542    -706.0389            0     2.163691    13.802484     3.022372      13990.5 
+     100   -631.97128   -662.02544    30.054164    1081.4378    -3497.564   -664.61825    39.275003   -703.89325            0    2.5928078    13.956833    2.5417699      13990.5 
+     150   -630.14953   -663.04215    32.892622    1183.5739    -88.43828   -665.63444    46.239965   -711.87441            0    2.5922927    14.667898    2.4964255      13990.5 
+     200   -628.52878    -663.9795     35.45072    1275.6219   -1755.9004   -666.73564    41.758052   -708.49369            0    2.7561421    14.230743    3.0924004      13990.5 
+     250   -627.27102     -662.025    34.753978    1250.5511   -1234.0918   -665.13519    43.170874   -708.30606            0    3.1101887    14.221086     1.941354      13990.5 
+     300    -626.5495   -663.74287    37.193368    1338.3275   -2049.3444   -666.45574    40.476148   -706.93188            0    2.7128711    13.330425    1.7756755      13990.5 
+     350   -625.87313   -665.21855    39.345421    1415.7647   -1543.1723   -667.90872    41.577366   -709.48609            0    2.6901682    13.541311     1.854662      13990.5 
+     400   -625.09344   -661.26404      36.1706    1301.5253   -729.96729   -664.10334    43.468765   -707.57211            0    2.8392963    13.663555    1.9067551      13990.5 
+     450   -624.46214   -660.01362    35.551477    1279.2474   -1617.7158   -663.06571    41.644856   -704.71057            0    3.0520921    14.527005    1.7280213      13990.5 
+     500   -623.49246    -659.2527     35.76024    1286.7593   -935.99238   -662.32953    43.038808   -705.36834            0    3.0768302    14.099593    1.9831106      13990.5 
+Loop time of 4.09864 on 4 procs for 500 steps with 432 atoms
+
+Performance: 21.080 ns/day, 1.139 hours/ns, 121.992 timesteps/s
+99.7% CPU use with 4 MPI tasks x no OpenMP threads
+
+MPI task timing breakdown:
+Section |  min time  |  avg time  |  max time  |%varavg| %total
+---------------------------------------------------------------
+Pair    | 3.3804     | 3.568      | 3.8354     |   8.9 | 87.05
+Bond    | 0.00074339 | 0.00079519 | 0.00087976 |   0.0 |  0.02
+Neigh   | 0.045851   | 0.046084   | 0.046361   |   0.1 |  1.12
+Comm    | 0.20413    | 0.47123    | 0.65875    |  24.3 | 11.50
+Output  | 0.00044298 | 0.00046057 | 0.00051165 |   0.0 |  0.01
+Modify  | 0.0064909  | 0.0067219  | 0.0069766  |   0.2 |  0.16
+Other   |            | 0.005345   |            |       |  0.13
+
+Nlocal:    108 ave 114 max 105 min
+Histogram: 1 1 1 0 0 0 0 0 0 1
+Nghost:    6527 ave 6599 max 6472 min
+Histogram: 1 0 1 0 1 0 0 0 0 1
+Neighs:    74388.2 ave 75855 max 73680 min
+Histogram: 1 2 0 0 0 0 0 0 0 1
+
+Total # of neighbors = 297553
+Ave neighs/atom = 688.78
+Ave special neighs/atom = 1
+Neighbor list builds = 20
+Dangerous builds = 0
+
+unfix thermoberendsen
+
+# ------------------------ Dynamic Run -------------------------------
+
+run 1000
+Memory usage per processor = 6.85787 Mbytes
+Step TotEng PotEng KinEng Temp Press E_pair E_vdwl E_coul E_long E_bond Fnorm Fmax Volume 
+     500   -623.49319    -659.2527    35.759511    1286.7331   -936.04802   -662.32953    43.038808   -705.36834            0    3.0768302    14.099593    1.9831106      13990.5 
+     550   -623.44059   -663.57938    40.138795    1444.3127   -935.73484    -666.2789    42.563337   -708.84224            0    2.6995167    13.918509    2.3189805      13990.5 
+     600    -623.4703   -660.01592    36.545618    1315.0196    1327.3492   -663.08845    47.985462   -711.07391            0    3.0725254    15.192713    2.4098428      13990.5 
+     650   -623.46796   -661.56776    38.099807    1370.9439    457.82439   -664.81976    45.495622   -710.31538            0    3.2519966    15.026057    1.8500226      13990.5 
+     700   -623.50158    -659.5131    36.011523    1295.8012   -460.03772    -663.1078    43.938203     -707.046            0    3.5946908    14.660979    2.4825518      13990.5 
+     750   -623.44787   -661.93353    38.485658    1384.8279    97.429626    -664.9551    45.083146   -710.03825            0    3.0215753     15.10043    2.3433897      13990.5 
+     800   -623.48215   -659.50655    36.024402    1296.2647    1097.3866   -662.61124    47.251998   -709.86324            0    3.1046914    14.556382    2.0543766      13990.5 
+     850   -623.45868   -661.13782    37.679134    1355.8068   -1802.1624   -664.41257     40.70845   -705.12102            0    3.2747525    14.691444    2.2054332      13990.5 
+     900   -623.43556   -663.59137    40.155815    1444.9251    534.99197   -666.71877    45.601619   -712.32039            0     3.127395    14.741411    2.5807895      13990.5 
+     950   -623.51318   -661.57916     38.06598    1369.7267   -678.12625   -664.37535    43.207862   -707.58322            0    2.7961988    14.430307    2.3936105      13990.5 
+    1000   -623.47287   -661.22274    37.749874    1358.3523     634.7979   -664.42973    46.373361   -710.80309            0    3.2069879    15.891192    2.4042765      13990.5 
+    1050   -623.48133   -661.52868    38.047347    1369.0562   -583.15228    -664.6098    43.618772   -708.22857            0     3.081116    14.806856    2.3447613      13990.5 
+    1100   -623.47867   -661.83761    38.358946    1380.2685    -868.9779    -664.8826     42.84846   -707.73106            0     3.044983     14.69567     2.399143      13990.5 
+    1150   -623.44713   -661.21299    37.765857    1358.9274    405.14554   -664.09567    45.578739    -709.6744            0    2.8826753    15.437367    3.1381305      13990.5 
+    1200   -623.46549   -660.91706    37.451568    1347.6183    699.78996    -664.0883     46.36297   -710.45127            0    3.1712473    15.109665    1.8891886      13990.5 
+    1250   -623.49296    -658.2218    34.728838    1249.6464    1061.0154   -661.29052    47.668699   -708.95922            0    3.0687228    14.901367    2.3964137      13990.5 
+    1300   -623.49837   -660.91022    37.411844    1346.1889    226.99512   -664.35989    45.352287   -709.71217            0    3.4496704    15.161542    2.2137993      13990.5 
+    1350   -623.46718   -658.80365    35.336469    1271.5108    1039.6469   -662.16908    47.565671   -709.73475            0    3.3654314    15.892516    2.7888426      13990.5 
+    1400   -623.47124   -661.45375    37.982513    1366.7233   -379.56023    -664.6321    43.788306   -708.42041            0    3.1783497    14.251126    1.7415409      13990.5 
+    1450   -623.46671   -660.17518    36.708464    1320.8792   -374.37056   -662.92706    44.083648   -707.01071            0    2.7518803    15.210167    1.9984277      13990.5 
+    1500   -623.50515   -659.06488    35.559725    1279.5442    260.37822   -662.39548    45.779764   -708.17524            0    3.3306005    14.682396    2.4201107      13990.5 
+Loop time of 8.26746 on 4 procs for 1000 steps with 432 atoms
+
+Performance: 20.901 ns/day, 1.148 hours/ns, 120.956 timesteps/s
+99.7% CPU use with 4 MPI tasks x no OpenMP threads
+
+MPI task timing breakdown:
+Section |  min time  |  avg time  |  max time  |%varavg| %total
+---------------------------------------------------------------
+Pair    | 6.706      | 7.1568     | 7.6597     |  12.7 | 86.57
+Bond    | 0.0014617  | 0.0015531  | 0.0016506  |   0.2 |  0.02
+Neigh   | 0.10511    | 0.10522    | 0.10532    |   0.0 |  1.27
+Comm    | 0.48547    | 0.98841    | 1.4393     |  34.0 | 11.96
+Output  | 0.0012085  | 0.0012462  | 0.0013196  |   0.1 |  0.02
+Modify  | 0.0021446  | 0.0021989  | 0.0022545  |   0.1 |  0.03
+Other   |            | 0.01204    |            |       |  0.15
+
+Nlocal:    108 ave 114 max 94 min
+Histogram: 1 0 0 0 0 0 0 0 1 2
+Nghost:    6512.25 ave 6586 max 6456 min
+Histogram: 1 0 0 2 0 0 0 0 0 1
+Neighs:    74248.2 ave 77441 max 65858 min
+Histogram: 1 0 0 0 0 0 0 0 0 3
+
+Total # of neighbors = 296993
+Ave neighs/atom = 687.484
+Ave special neighs/atom = 1
+Neighbor list builds = 46
+Dangerous builds = 0
+Total wall time: 0:00:12
--- a/potentials/Ge.tersoff
+++ b/potentials/Ge.tersoff
@ -0,0 +1,17 @@
+# DATE: 2016-12-21 CONTRIBUTOR: Sayyed Jalil Mahdizadh, saja.mahdizadeh@gmail.com CITATION: Mahdizadeh, Akhlamadi, J. Mol. Graph. Model. 72, 1-5 (2017)  
+
+# Tersoff parameters for various elements and mixtures
+# multiple entries can be added to this file, LAMMPS reads the ones it needs
+# these entries are in LAMMPS "metal" units:
+#   A,B = eV; lambda1,lambda2,lambda3 = 1/Angstroms; R,D = Angstroms
+#   other quantities are unitless
+
+# This is the Ge reparameterization from the following paper:
+# Optimized  Tersoff  empirical  potential for germanene
+# Mahdizadeh, Akhlamadi, J. Mol. Graph. Model. 72, 1-5 (2017)
+
+# format of a single entry (one or more lines):
+#   element 1, element 2, element 3, 
+#   m, gamma, lambda3, c, d, costheta0, n, beta, lambda2, B, R, D, lambda1, A
+
+Ge     Ge      Ge        3.0 1.0 0.0 1.0643e5 15.2 -0.35 0.75627 5.017e-7  1.71  430.0    2.95    0.15   2.4451  1760.1
--- a/potentials/WL.meam
+++ b/potentials/WL.meam
@ -0,0 +1,13 @@
+# DATE: 2017-01-25 CONTRIBUTOR: Aidan Thompson, athomps@sandia.gov,  CITATION: Lee, Baskes, Kim, Cho. Phys. Rev. B, 64, 184102 (2001) 
+rc = 3.8
+delr = 0.1
+augt1 = 0
+erose_form = 2
+zbl(1,1) = 0
+nn2(1,1) = 1
+Ec(1,1) = 8.66
+re(1,1) = 2.74
+attrac(1,1) = 0
+repuls(1,1) = 0
+Cmin(1,1,1) = 0.49
+Cmax(1,1,1) = 2.8
--- a/src/.gitignore
+++ b/src/.gitignore
@ -154,6 +154,8 @@
 /bond_morse.h
 /bond_nonlinear.cpp
 /bond_nonlinear.h
+/bond_oxdna_fene.cpp
+/bond_oxdna_fene.h
 /bond_quartic.cpp
 /bond_quartic.h
 /bond_table.cpp
@ -390,6 +392,10 @@
 /fix_nve_asphere.h
 /fix_nve_asphere_noforce.cpp
 /fix_nve_asphere_noforce.h
+/fix_nve_dot.cpp
+/fix_nve_dot.h
+/fix_nve_dotc_langevin.cpp
+/fix_nve_dotc_langevin.h
 /fix_nh_body.cpp
 /fix_nh_body.h
 /fix_nph_body.cpp
@ -751,6 +757,9 @@
 /pair_nm_cut_coul_cut.h
 /pair_nm_cut_coul_long.cpp
 /pair_nm_cut_coul_long.h
+/pair_oxdna_*.cpp
+/pair_oxdna_*.h
+/mf_oxdna.h
 /pair_peri_eps.cpp
 /pair_peri_eps.h
 /pair_peri_lps.cpp
--- a/src/Depend.sh
+++ b/src/Depend.sh
@ -48,6 +48,7 @@ depend () {
 if (test $1 = "ASPHERE") then
  depend GPU
  depend USER-OMP
+  depend USER-CGDNA
  depend USER-INTEL
 fi

@ -96,6 +97,7 @@ if (test $1 = "MOLECULE") then
  depend USER-MISC
  depend USER-OMP
  depend USER-FEP
+  depend USER-CGDNA
  depend USER-INTEL
 fi

--- a/src/GRANULAR/fix_wall_gran_region.cpp
+++ b/src/GRANULAR/fix_wall_gran_region.cpp
@ -278,7 +278,7 @@ void FixWallGranRegion::update_contacts(int i, int nc)
  iold = 0;
  while (iold < ncontact[i]) {
    for (m = 0; m < nc; m++)
-      if (region->contact[m].iwall = walls[i][iold]) break;
+      if (region->contact[m].iwall == walls[i][iold]) break;
    if (m < nc) {
      ilast = ncontact[i]-1;
      for (j = 0; j < sheardim; j++)
--- a/src/GRANULAR/pair_gran_hooke_history.cpp
+++ b/src/GRANULAR/pair_gran_hooke_history.cpp
@ -799,3 +799,14 @@ double PairGranHookeHistory::memory_usage()
  double bytes = nmax * sizeof(double);
  return bytes;
 }
+
+/* ----------------------------------------------------------------------
+   return ptr to FixShearHistory class
+   called by Neighbor when setting up neighbor lists
+------------------------------------------------------------------------- */
+
+void *PairGranHookeHistory::extract(const char *str, int &dim)
+{
+  if (strcmp(str,"history") == 0) return (void *) fix_history;
+  return NULL;
+}
--- a/src/GRANULAR/pair_gran_hooke_history.h
+++ b/src/GRANULAR/pair_gran_hooke_history.h
@ -43,6 +43,7 @@ class PairGranHookeHistory : public Pair {
  int pack_forward_comm(int, int *, double *, int, int *);
  void unpack_forward_comm(int, int, double *);
  double memory_usage();
+  void *extract(const char *, int &);

 protected:
  double kn,kt,gamman,gammat,xmu;
--- a/src/KOKKOS/Install.sh
+++ b/src/KOKKOS/Install.sh
@ -173,6 +173,8 @@ action pair_lj_gromacs_kokkos.cpp
 action pair_lj_gromacs_kokkos.h
 action pair_lj_sdk_kokkos.cpp pair_lj_sdk.cpp
 action pair_lj_sdk_kokkos.h pair_lj_sdk.h
+action pair_morse_kokkos.cpp
+action pair_morse_kokkos.h
 action pair_reax_c_kokkos.cpp pair_reax_c.cpp
 action pair_reax_c_kokkos.h pair_reax_c.h
 action pair_sw_kokkos.cpp pair_sw.cpp
--- a/src/KOKKOS/angle_charmm_kokkos.cpp
+++ b/src/KOKKOS/angle_charmm_kokkos.cpp
@ -70,23 +70,20 @@ void AngleCharmmKokkos<DeviceType>::compute(int eflag_in, int vflag_in)
  // reallocate per-atom arrays if necessary

  if (eflag_atom) {
-    if(k_eatom.dimension_0()<maxeatom) {
+    //if(k_eatom.dimension_0()<maxeatom) { // won't work without adding zero functor
      memory->destroy_kokkos(k_eatom,eatom);
      memory->create_kokkos(k_eatom,eatom,maxeatom,"improper:eatom");
-      d_eatom = k_eatom.d_view;
-    }
+      d_eatom = k_eatom.template view<DeviceType>();
+    //}
  }
  if (vflag_atom) {
-    if(k_vatom.dimension_0()<maxvatom) {
+    //if(k_vatom.dimension_0()<maxvatom) { // won't work without adding zero functor
      memory->destroy_kokkos(k_vatom,vatom);
      memory->create_kokkos(k_vatom,vatom,maxvatom,6,"improper:vatom");
-      d_vatom = k_vatom.d_view;
-    }
+      d_vatom = k_vatom.template view<DeviceType>();
+    //}
  }

-  if (eflag || vflag) atomKK->modified(execution_space,datamask_modify);
-  else atomKK->modified(execution_space,F_MASK);
-
  x = atomKK->k_x.view<DeviceType>();
  f = atomKK->k_f.view<DeviceType>();
  neighborKK->k_anglelist.template sync<DeviceType>();
--- a/src/KOKKOS/atom_vec_angle_kokkos.cpp
+++ b/src/KOKKOS/atom_vec_angle_kokkos.cpp
@ -957,10 +957,10 @@ struct AtomVecAngleKokkos_UnpackBorder {
      _x(i+_first,0) = _buf(i,0);
      _x(i+_first,1) = _buf(i,1);
      _x(i+_first,2) = _buf(i,2);
-      _tag(i+_first) = static_cast<int> (_buf(i,3));
+      _tag(i+_first) = static_cast<tagint> (_buf(i,3));
      _type(i+_first) = static_cast<int>  (_buf(i,4));
      _mask(i+_first) = static_cast<int>  (_buf(i,5));
-      _molecule(i+_first) = static_cast<int> (_buf(i,6));
+      _molecule(i+_first) = static_cast<tagint> (_buf(i,6));

  }
 };
--- a/src/KOKKOS/atom_vec_atomic_kokkos.cpp
+++ b/src/KOKKOS/atom_vec_atomic_kokkos.cpp
@ -836,7 +836,7 @@ struct AtomVecAtomicKokkos_UnpackBorder {
      _x(i+_first,0) = _buf(i,0);
      _x(i+_first,1) = _buf(i,1);
      _x(i+_first,2) = _buf(i,2);
-      _tag(i+_first) = static_cast<int> (_buf(i,3));
+      _tag(i+_first) = static_cast<tagint> (_buf(i,3));
      _type(i+_first) = static_cast<int>  (_buf(i,4));
      _mask(i+_first) = static_cast<int>  (_buf(i,5));
 //      printf("%i %i %lf %lf %lf %i BORDER\n",_tag(i+_first),i+_first,_x(i+_first,0),_x(i+_first,1),_x(i+_first,2),_type(i+_first));
--- a/src/KOKKOS/atom_vec_bond_kokkos.cpp
+++ b/src/KOKKOS/atom_vec_bond_kokkos.cpp
@ -905,10 +905,10 @@ struct AtomVecBondKokkos_UnpackBorder {
      _x(i+_first,0) = _buf(i,0);
      _x(i+_first,1) = _buf(i,1);
      _x(i+_first,2) = _buf(i,2);
-      _tag(i+_first) = static_cast<int> (_buf(i,3));
+      _tag(i+_first) = static_cast<tagint> (_buf(i,3));
      _type(i+_first) = static_cast<int>  (_buf(i,4));
      _mask(i+_first) = static_cast<int>  (_buf(i,5));
-      _molecule(i+_first) = static_cast<int> (_buf(i,6));
+      _molecule(i+_first) = static_cast<tagint> (_buf(i,6));

  }
 };
--- a/src/KOKKOS/atom_vec_charge_kokkos.cpp
+++ b/src/KOKKOS/atom_vec_charge_kokkos.cpp
@ -872,7 +872,7 @@ struct AtomVecChargeKokkos_UnpackBorder {
      _x(i+_first,0) = _buf(i,0);
      _x(i+_first,1) = _buf(i,1);
      _x(i+_first,2) = _buf(i,2);
-      _tag(i+_first) = static_cast<int> (_buf(i,3));
+      _tag(i+_first) = static_cast<tagint> (_buf(i,3));
      _type(i+_first) = static_cast<int>  (_buf(i,4));
      _mask(i+_first) = static_cast<int>  (_buf(i,5));
      _q(i+_first) = _buf(i,6);
--- a/src/KOKKOS/atom_vec_molecular_kokkos.cpp
+++ b/src/KOKKOS/atom_vec_molecular_kokkos.cpp
@ -1029,10 +1029,10 @@ struct AtomVecMolecularKokkos_UnpackBorder {
      _x(i+_first,0) = _buf(i,0);
      _x(i+_first,1) = _buf(i,1);
      _x(i+_first,2) = _buf(i,2);
-      _tag(i+_first) = static_cast<int> (_buf(i,3));
+      _tag(i+_first) = static_cast<tagint> (_buf(i,3));
      _type(i+_first) = static_cast<int>  (_buf(i,4));
      _mask(i+_first) = static_cast<int>  (_buf(i,5));
-      _molecule(i+_first) = static_cast<int> (_buf(i,6));
+      _molecule(i+_first) = static_cast<tagint> (_buf(i,6));

  }
 };
--- a/src/KOKKOS/bond_fene_kokkos.cpp
+++ b/src/KOKKOS/bond_fene_kokkos.cpp
@ -77,21 +77,18 @@ void BondFENEKokkos<DeviceType>::compute(int eflag_in, int vflag_in)
  if (eflag_atom) {
    memory->destroy_kokkos(k_eatom,eatom);
    memory->create_kokkos(k_eatom,eatom,maxeatom,"bond:eatom");
-    d_eatom = k_eatom.d_view;
+    d_eatom = k_eatom.view<DeviceType>();
  }
  if (vflag_atom) {
    memory->destroy_kokkos(k_vatom,vatom);
    memory->create_kokkos(k_vatom,vatom,maxvatom,6,"bond:vatom");
-    d_vatom = k_vatom.d_view;
+    d_vatom = k_vatom.view<DeviceType>();
  }

-  atomKK->sync(execution_space,datamask_read);
  k_k.template sync<DeviceType>();
  k_r0.template sync<DeviceType>();
  k_epsilon.template sync<DeviceType>();
  k_sigma.template sync<DeviceType>();
-  if (eflag || vflag) atomKK->modified(execution_space,datamask_modify);
-  else atomKK->modified(execution_space,F_MASK);

  x = atomKK->k_x.view<DeviceType>();
  f = atomKK->k_f.view<DeviceType>();
--- a/src/KOKKOS/bond_fene_kokkos.h
+++ b/src/KOKKOS/bond_fene_kokkos.h
@ -66,8 +66,8 @@ class BondFENEKokkos : public BondFENE {

  DAT::tdual_efloat_1d k_eatom;
  DAT::tdual_virial_array k_vatom;
-  DAT::t_efloat_1d d_eatom;
-  DAT::t_virial_array d_vatom;
+  typename ArrayTypes<DeviceType>::t_efloat_1d d_eatom;
+  typename ArrayTypes<DeviceType>::t_virial_array d_vatom;

  DAT::tdual_int_scalar k_warning_flag;
  typename AT::t_int_scalar d_warning_flag;
--- a/src/KOKKOS/bond_harmonic_kokkos.cpp
+++ b/src/KOKKOS/bond_harmonic_kokkos.cpp
@ -67,18 +67,18 @@ void BondHarmonicKokkos<DeviceType>::compute(int eflag_in, int vflag_in)
  // reallocate per-atom arrays if necessary

  if (eflag_atom) {
-    if(k_eatom.dimension_0()<maxeatom) {
+    //if(k_eatom.dimension_0()<maxeatom) { // won't work without adding zero functor
      memory->destroy_kokkos(k_eatom,eatom);
      memory->create_kokkos(k_eatom,eatom,maxeatom,"improper:eatom");
-      d_eatom = k_eatom.d_view;
-    }
+      d_eatom = k_eatom.template view<DeviceType>();
+    //}
  }
  if (vflag_atom) {
-    if(k_vatom.dimension_0()<maxvatom) {
+    //if(k_vatom.dimension_0()<maxvatom) { // won't work without adding zero functor
      memory->destroy_kokkos(k_vatom,vatom);
      memory->create_kokkos(k_vatom,vatom,maxvatom,6,"improper:vatom");
-      d_vatom = k_vatom.d_view;
-    }
+      d_vatom = k_vatom.template view<DeviceType>();
+    //}
  }

 //  if (eflag || vflag) atomKK->modified(execution_space,datamask_modify);
--- a/src/KOKKOS/dihedral_charmm_kokkos.cpp
+++ b/src/KOKKOS/dihedral_charmm_kokkos.cpp
@ -80,29 +80,24 @@ void DihedralCharmmKokkos<DeviceType>::compute(int eflag_in, int vflag_in)
  // reallocate per-atom arrays if necessary

  if (eflag_atom) {
-    if(k_eatom.dimension_0()<maxeatom) {
+    //if(k_eatom.dimension_0()<maxeatom) { // won't work without adding zero functor
      memory->destroy_kokkos(k_eatom,eatom);
      memory->create_kokkos(k_eatom,eatom,maxeatom,"dihedral:eatom");
-      d_eatom = k_eatom.d_view;
+      d_eatom = k_eatom.template view<DeviceType>();
      k_eatom_pair = Kokkos::DualView<E_FLOAT*,Kokkos::LayoutRight,DeviceType>("dihedral:eatom_pair",maxeatom);
-      d_eatom_pair = k_eatom.d_view;
-    }
+      d_eatom_pair = k_eatom.template view<DeviceType>();
+    //}
  }
  if (vflag_atom) {
-    if(k_vatom.dimension_0()<maxvatom) {
+    //if(k_vatom.dimension_0()<maxvatom) { // won't work without adding zero functor
      memory->destroy_kokkos(k_vatom,vatom);
      memory->create_kokkos(k_vatom,vatom,maxvatom,6,"dihedral:vatom");
-      d_vatom = k_vatom.d_view;
+      d_vatom = k_vatom.template view<DeviceType>();
      k_vatom_pair = Kokkos::DualView<F_FLOAT*[6],Kokkos::LayoutRight,DeviceType>("dihedral:vatom_pair",maxvatom);
-      d_vatom_pair = k_vatom.d_view;
-    }
+      d_vatom_pair = k_vatom.template view<DeviceType>();
+    //}
  }

-
-  //atomKK->sync(execution_space,datamask_read);
-  if (eflag || vflag) atomKK->modified(execution_space,datamask_modify);
-  else atomKK->modified(execution_space,F_MASK);
-
  x = atomKK->k_x.view<DeviceType>();
  f = atomKK->k_f.view<DeviceType>();
  q = atomKK->k_q.view<DeviceType>();
--- a/src/KOKKOS/dihedral_opls_kokkos.cpp
+++ b/src/KOKKOS/dihedral_opls_kokkos.cpp
@ -77,21 +77,18 @@ void DihedralOPLSKokkos<DeviceType>::compute(int eflag_in, int vflag_in)
  if (eflag_atom) {
    memory->destroy_kokkos(k_eatom,eatom);
    memory->create_kokkos(k_eatom,eatom,maxeatom,"dihedral:eatom");
-    d_eatom = k_eatom.d_view;
+    d_eatom = k_eatom.view<DeviceType>();
  }
  if (vflag_atom) {
    memory->destroy_kokkos(k_vatom,vatom);
    memory->create_kokkos(k_vatom,vatom,maxvatom,6,"dihedral:vatom");
-    d_vatom = k_vatom.d_view;
+    d_vatom = k_vatom.view<DeviceType>();
  }

-  atomKK->sync(execution_space,datamask_read);
  k_k1.template sync<DeviceType>();
  k_k2.template sync<DeviceType>();
  k_k3.template sync<DeviceType>();
  k_k4.template sync<DeviceType>();
-  if (eflag || vflag) atomKK->modified(execution_space,datamask_modify);
-  else atomKK->modified(execution_space,F_MASK);

  x = atomKK->k_x.view<DeviceType>();
  f = atomKK->k_f.view<DeviceType>();
--- a/src/KOKKOS/dihedral_opls_kokkos.h
+++ b/src/KOKKOS/dihedral_opls_kokkos.h
@ -68,8 +68,8 @@ class DihedralOPLSKokkos : public DihedralOPLS {

  DAT::tdual_efloat_1d k_eatom;
  DAT::tdual_virial_array k_vatom;
-  DAT::t_efloat_1d d_eatom;
-  DAT::t_virial_array d_vatom;
+  typename ArrayTypes<DeviceType>::t_efloat_1d d_eatom;
+  typename ArrayTypes<DeviceType>::t_virial_array d_vatom;

  int nlocal,newton_bond;
  int eflag,vflag;
--- a/src/KOKKOS/domain_kokkos.cpp
+++ b/src/KOKKOS/domain_kokkos.cpp
@ -250,7 +250,7 @@ struct DomainPBCFunctor {
        x(i,0) += period[0];
        if (DEFORM_VREMAP && (mask[i] & deform_groupbit)) v(i,0) += h_rate[0];
        imageint idim = image[i] & IMGMASK;
-        const int otherdims = image[i] ^ idim;
+        const imageint otherdims = image[i] ^ idim;
        idim--;
        idim &= IMGMASK;
        image[i] = otherdims | idim;
@ -260,7 +260,7 @@ struct DomainPBCFunctor {
        x(i,0) = MAX(x(i,0),lo[0]);
        if (DEFORM_VREMAP && (mask[i] & deform_groupbit)) v(i,0) -= h_rate[0];
        imageint idim = image[i] & IMGMASK;
-        const int otherdims = image[i] ^ idim;
+        const imageint otherdims = image[i] ^ idim;
        idim++;
        idim &= IMGMASK;
        image[i] = otherdims | idim;
@ -354,7 +354,6 @@ void DomainKokkos::pbc()
  }

  atomKK->sync(Device,X_MASK|V_MASK|MASK_MASK|IMAGE_MASK);
-  atomKK->modified(Device,X_MASK|V_MASK|IMAGE_MASK);

  if (xperiodic || yperiodic || zperiodic) {
    if (deform_vremap) {
@ -385,8 +384,9 @@ void DomainKokkos::pbc()
      Kokkos::parallel_for(nlocal,f);
    }
  }
-
  LMPDeviceType::fence();
+
+  atomKK->modified(Device,X_MASK|V_MASK|IMAGE_MASK);
 }

 /* ----------------------------------------------------------------------
@ -402,7 +402,6 @@ void DomainKokkos::pbc()
 void DomainKokkos::remap_all()
 {
  atomKK->sync(Device,X_MASK | IMAGE_MASK);
-  atomKK->modified(Device,X_MASK | IMAGE_MASK);

  x = atomKK->k_x.view<LMPDeviceType>();
  image = atomKK->k_image.view<LMPDeviceType>();
@ -428,6 +427,8 @@ void DomainKokkos::remap_all()
  LMPDeviceType::fence();
  copymode = 0;

+  atomKK->modified(Device,X_MASK | IMAGE_MASK);
+
  if (triclinic) lamda2x(nlocal);
 }

@ -521,7 +522,6 @@ void DomainKokkos::image_flip(int m_in, int n_in, int p_in)
  p_flip = p_in;

  atomKK->sync(Device,IMAGE_MASK);
-  atomKK->modified(Device,IMAGE_MASK);

  image = atomKK->k_image.view<LMPDeviceType>();
  int nlocal = atomKK->nlocal;
@ -530,6 +530,8 @@ void DomainKokkos::image_flip(int m_in, int n_in, int p_in)
  Kokkos::parallel_for(Kokkos::RangePolicy<LMPDeviceType, TagDomain_image_flip>(0,nlocal),*this);
  LMPDeviceType::fence();
  copymode = 0;
+
+  atomKK->modified(Device,IMAGE_MASK);
 }

 KOKKOS_INLINE_FUNCTION
@ -554,7 +556,6 @@ void DomainKokkos::operator()(TagDomain_image_flip, const int &i) const {
 void DomainKokkos::lamda2x(int n)
 {
  atomKK->sync(Device,X_MASK);
-  atomKK->modified(Device,X_MASK);

  x = atomKK->k_x.view<LMPDeviceType>();

@ -562,6 +563,8 @@ void DomainKokkos::lamda2x(int n)
  Kokkos::parallel_for(Kokkos::RangePolicy<LMPDeviceType, TagDomain_lamda2x>(0,n),*this);
  LMPDeviceType::fence();
  copymode = 0;
+
+  atomKK->modified(Device,X_MASK);
 }

 KOKKOS_INLINE_FUNCTION
@ -579,7 +582,6 @@ void DomainKokkos::operator()(TagDomain_lamda2x, const int &i) const {
 void DomainKokkos::x2lamda(int n)
 {
  atomKK->sync(Device,X_MASK);
-  atomKK->modified(Device,X_MASK);

  x = atomKK->k_x.view<LMPDeviceType>();

@ -587,6 +589,8 @@ void DomainKokkos::x2lamda(int n)
  Kokkos::parallel_for(Kokkos::RangePolicy<LMPDeviceType, TagDomain_x2lamda>(0,n),*this);
  LMPDeviceType::fence();
  copymode = 0;
+
+  atomKK->modified(Device,X_MASK);
 }

 KOKKOS_INLINE_FUNCTION
--- a/src/KOKKOS/fix_langevin_kokkos.cpp
+++ b/src/KOKKOS/fix_langevin_kokkos.cpp
@ -177,11 +177,12 @@ void FixLangevinKokkos<DeviceType>::post_force(int vflag)

  // account for bias velocity
  if(tbiasflag == BIAS){
+    atomKK->sync(temperature->execution_space,temperature->datamask_read);
    temperature->compute_scalar();
    temperature->remove_bias_all(); // modifies velocities
    // if temeprature compute is kokkosized host-devcie comm won't be needed
-    atomKK->modified(Host,V_MASK);
-    atomKK->sync(execution_space,V_MASK);
+    atomKK->modified(temperature->execution_space,temperature->datamask_modify);
+    atomKK->sync(execution_space,temperature->datamask_modify);
  }

  // compute langevin force in parallel on the device
@ -508,8 +509,10 @@ void FixLangevinKokkos<DeviceType>::post_force(int vflag)
  DeviceType::fence();

  if(tbiasflag == BIAS){
+    atomKK->sync(temperature->execution_space,temperature->datamask_read);
    temperature->restore_bias_all(); // modifies velocities
-    atomKK->modified(Host,V_MASK);
+    atomKK->modified(temperature->execution_space,temperature->datamask_modify);
+    atomKK->sync(execution_space,temperature->datamask_modify);
  }

  // set modify flags for the views modified in post_force functor
--- a/src/KOKKOS/fix_qeq_reax_kokkos.cpp
+++ b/src/KOKKOS/fix_qeq_reax_kokkos.cpp
@ -387,7 +387,7 @@ KOKKOS_INLINE_FUNCTION
 void FixQEqReaxKokkos<DeviceType>::compute_h_item(int ii, int &m_fill, const bool &final) const
 {
  const int i = d_ilist[ii];
-  int j,jj,jtag,jtype,flag;
+  int j,jj,jtype,flag;

  if (mask[i] & groupbit) {

@ -395,7 +395,7 @@ void FixQEqReaxKokkos<DeviceType>::compute_h_item(int ii, int &m_fill, const boo
    const X_FLOAT ytmp = x(i,1);
    const X_FLOAT ztmp = x(i,2);
    const int itype = type(i);
-    const int itag = tag(i);
+    const tagint itag = tag(i);
    const int jnum = d_numneigh[i];
    if (final)
      d_firstnbr[i] = m_fill;
@ -403,7 +403,6 @@ void FixQEqReaxKokkos<DeviceType>::compute_h_item(int ii, int &m_fill, const boo
    for (jj = 0; jj < jnum; jj++) {
      j = d_neighbors(i,jj);
      j &= NEIGHMASK;
-
      jtype = type(j);

      const X_FLOAT delx = x(j,0) - xtmp;
@ -411,10 +410,11 @@ void FixQEqReaxKokkos<DeviceType>::compute_h_item(int ii, int &m_fill, const boo
      const X_FLOAT delz = x(j,2) - ztmp;

      if (neighflag != FULL) {
+        const tagint jtag = tag(j);
        flag = 0;
        if (j < nlocal) flag = 1;
-        else if (tag[i] < tag[j]) flag = 1;
-        else if (tag[i] == tag[j]) {
+        else if (itag < jtag) flag = 1;
+        else if (itag == jtag) {
          if (delz > SMALL) flag = 1;
          else if (fabs(delz) < SMALL) {
            if (dely > SMALL) flag = 1;
--- a/src/KOKKOS/fix_qeq_reax_kokkos.h
+++ b/src/KOKKOS/fix_qeq_reax_kokkos.h
@ -159,7 +159,8 @@ class FixQEqReaxKokkos : public FixQEqReax {
  //typename ArrayTypes<DeviceType>::t_float_1d_randomread mass, q;
  typename ArrayTypes<DeviceType>::t_float_1d_randomread mass;
  typename ArrayTypes<DeviceType>::t_float_1d q;
-  typename ArrayTypes<DeviceType>::t_int_1d type, tag, mask;
+  typename ArrayTypes<DeviceType>::t_int_1d type, mask;
+  typename ArrayTypes<DeviceType>::t_tagint_1d tag;

  DAT::tdual_float_1d k_q;
  typename AT::t_float_1d d_q;
--- a/src/KOKKOS/fix_reaxc_bonds_kokkos.cpp
+++ b/src/KOKKOS/fix_reaxc_bonds_kokkos.cpp
@ -76,7 +76,6 @@ void FixReaxCBondsKokkos::init()
 void FixReaxCBondsKokkos::Output_ReaxC_Bonds(bigint ntimestep, FILE *fp)

 {
-  int i, j;
  int nbuf_local;
  int nlocal_max, numbonds, numbonds_max;
  double *buf;
--- a/src/KOKKOS/improper_harmonic_kokkos.cpp
+++ b/src/KOKKOS/improper_harmonic_kokkos.cpp
@ -77,18 +77,18 @@ void ImproperHarmonicKokkos<DeviceType>::compute(int eflag_in, int vflag_in)
  // reallocate per-atom arrays if necessary

  if (eflag_atom) {
-    if(k_eatom.dimension_0()<maxeatom) {
+    //if(k_eatom.dimension_0()<maxeatom) { // won't work without adding zero functor
      memory->destroy_kokkos(k_eatom,eatom);
      memory->create_kokkos(k_eatom,eatom,maxeatom,"improper:eatom");
-      d_eatom = k_eatom.d_view;
-    }
+      d_eatom = k_eatom.template view<DeviceType>();
+    //}
  }
  if (vflag_atom) {
-    if(k_vatom.dimension_0()<maxvatom) {
+    //if(k_vatom.dimension_0()<maxvatom) { // won't work without adding zero functor
      memory->destroy_kokkos(k_vatom,vatom);
      memory->create_kokkos(k_vatom,vatom,maxvatom,6,"improper:vatom");
-      d_vatom = k_vatom.d_view;
-    }
+      d_vatom = k_vatom.template view<DeviceType>();
+    //}
  }

  //atomKK->sync(execution_space,datamask_read);
--- a/src/KOKKOS/math_special_kokkos.h
+++ b/src/KOKKOS/math_special_kokkos.h
@ -42,12 +42,11 @@ namespace MathSpecialKokkos {
  {
    x *= x;
    x *= 1.4426950408889634074; // log_2(e)
-#if defined(__BYTE_ORDER__)
-#if __BYTE_ORDER__ == __ORDER_LITTLE_ENDIAN__
+#if defined(__BYTE_ORDER__) &&  __BYTE_ORDER__ == __ORDER_LITTLE_ENDIAN__
    return (x < 1023.0) ? exp2_x86(-x) : 0.0;
-#endif
-#endif
+#else
    return (x < 1023.0) ? exp2(-x) : 0.0;
+#endif
  }

  // x**2, use instead of pow(x,2.0)
--- a/src/KOKKOS/modify_kokkos.cpp
+++ b/src/KOKKOS/modify_kokkos.cpp
@ -360,9 +360,7 @@ void ModifyKokkos::post_run()
  for (int i = 0; i < nfix; i++) {
    atomKK->sync(fix[i]->execution_space,
                 fix[i]->datamask_read);
-    if (!fix[i]->kokkosable) lmp->kokkos->auto_sync = 1;
    fix[i]->post_run();
-    lmp->kokkos->auto_sync = 0;
    atomKK->modified(fix[i]->execution_space,
                     fix[i]->datamask_modify);
  }
--- a/src/KOKKOS/nbin_kokkos.cpp
+++ b/src/KOKKOS/nbin_kokkos.cpp
@ -74,10 +74,7 @@ void NBinKokkos<DeviceType>::bin_atoms_setup(int nall)

    k_bincount = DAT::tdual_int_1d("Neighbor::d_bincount",mbins);
    bincount = k_bincount.view<DeviceType>();
-    last_bin_memory = update->ntimestep;
  }
-
-  last_bin = update->ntimestep;
 }

 /* ----------------------------------------------------------------------
@ -87,6 +84,8 @@ void NBinKokkos<DeviceType>::bin_atoms_setup(int nall)
 template<class DeviceType>
 void NBinKokkos<DeviceType>::bin_atoms()
 {
+  last_bin = update->ntimestep;
+
  h_resize() = 1;

  while(h_resize() > 0) {
@ -116,7 +115,6 @@ void NBinKokkos<DeviceType>::bin_atoms()
      k_bins = DAT::tdual_int_2d("bins", mbins, atoms_per_bin);
      bins = k_bins.view<DeviceType>();
      c_bins = bins;
-      last_bin_memory = update->ntimestep;
    }
  }
 }
--- a/src/KOKKOS/neighbor_kokkos.cpp
+++ b/src/KOKKOS/neighbor_kokkos.cpp
@ -147,9 +147,9 @@ void NeighborKokkos::init_ex_mol_bit_kokkos()
 int NeighborKokkos::check_distance()
 {
  if (device_flag)
-    check_distance_kokkos<LMPDeviceType>();
+    return check_distance_kokkos<LMPDeviceType>();
  else
-    check_distance_kokkos<LMPHostType>();
+    return check_distance_kokkos<LMPHostType>();
 }

 template<class DeviceType>
@ -157,7 +157,7 @@ int NeighborKokkos::check_distance_kokkos()
 {
  typedef DeviceType device_type;

-  double delx,dely,delz,rsq;
+  double delx,dely,delz;
  double delta,delta1,delta2;

  if (boxcheck) {
--- a/src/KOKKOS/npair_kokkos.cpp
+++ b/src/KOKKOS/npair_kokkos.cpp
@ -24,8 +24,8 @@ namespace LAMMPS_NS {

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

-template<class DeviceType, int HALF_NEIGH, int GHOST>
-NPairKokkos<DeviceType,HALF_NEIGH,GHOST>::NPairKokkos(LAMMPS *lmp) : NPair(lmp) {
+template<class DeviceType, int HALF_NEIGH, int GHOST, int TRI>
+NPairKokkos<DeviceType,HALF_NEIGH,GHOST,TRI>::NPairKokkos(LAMMPS *lmp) : NPair(lmp) {

 }

@ -33,8 +33,8 @@ NPairKokkos<DeviceType,HALF_NEIGH,GHOST>::NPairKokkos(LAMMPS *lmp) : NPair(lmp)
   copy needed info from Neighbor class to this build class
   ------------------------------------------------------------------------- */

-template<class DeviceType, int HALF_NEIGH, int GHOST>
-void NPairKokkos<DeviceType,HALF_NEIGH,GHOST>::copy_neighbor_info()
+template<class DeviceType, int HALF_NEIGH, int GHOST, int TRI>
+void NPairKokkos<DeviceType,HALF_NEIGH,GHOST,TRI>::copy_neighbor_info()
 {
  NPair::copy_neighbor_info();

@ -62,8 +62,8 @@ void NPairKokkos<DeviceType,HALF_NEIGH,GHOST>::copy_neighbor_info()
 copy per-atom and per-bin vectors from NBin class to this build class
 ------------------------------------------------------------------------- */

-template<class DeviceType, int HALF_NEIGH, int GHOST>
-void NPairKokkos<DeviceType,HALF_NEIGH,GHOST>::copy_bin_info()
+template<class DeviceType, int HALF_NEIGH, int GHOST, int TRI>
+void NPairKokkos<DeviceType,HALF_NEIGH,GHOST,TRI>::copy_bin_info()
 {
  NPair::copy_bin_info();

@ -78,8 +78,8 @@ void NPairKokkos<DeviceType,HALF_NEIGH,GHOST>::copy_bin_info()
 copy needed info from NStencil class to this build class
 ------------------------------------------------------------------------- */

-template<class DeviceType, int HALF_NEIGH, int GHOST>
-void NPairKokkos<DeviceType,HALF_NEIGH,GHOST>::copy_stencil_info()
+template<class DeviceType, int HALF_NEIGH, int GHOST, int TRI>
+void NPairKokkos<DeviceType,HALF_NEIGH,GHOST,TRI>::copy_stencil_info()
 {
  NPair::copy_stencil_info();

@ -106,8 +106,8 @@ void NPairKokkos<DeviceType,HALF_NEIGH,GHOST>::copy_stencil_info()

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

-template<class DeviceType, int HALF_NEIGH, int GHOST>
-void NPairKokkos<DeviceType,HALF_NEIGH,GHOST>::build(NeighList *list_)
+template<class DeviceType, int HALF_NEIGH, int GHOST, int TRI>
+void NPairKokkos<DeviceType,HALF_NEIGH,GHOST,TRI>::build(NeighList *list_)
 {
  NeighListKokkos<DeviceType>* list = (NeighListKokkos<DeviceType>*) list_;
  const int nlocal = includegroup?atom->nfirst:atom->nlocal;
@ -196,14 +196,14 @@ if (GHOST) {
  Kokkos::parallel_for(nall, f);
 } else {
  if (newton_pair) {
-    NPairKokkosBuildFunctor<DeviceType,HALF_NEIGH,1> f(data,atoms_per_bin * 5 * sizeof(X_FLOAT) * factor);
+    NPairKokkosBuildFunctor<DeviceType,TRI?0:HALF_NEIGH,1,TRI> f(data,atoms_per_bin * 5 * sizeof(X_FLOAT) * factor);
 #ifdef KOKKOS_HAVE_CUDA
    Kokkos::parallel_for(config, f);
 #else
    Kokkos::parallel_for(nall, f);
 #endif
  } else {
-    NPairKokkosBuildFunctor<DeviceType,HALF_NEIGH,0> f(data,atoms_per_bin * 5 * sizeof(X_FLOAT) * factor);
+    NPairKokkosBuildFunctor<DeviceType,HALF_NEIGH,0,0> f(data,atoms_per_bin * 5 * sizeof(X_FLOAT) * factor);
 #ifdef KOKKOS_HAVE_CUDA
    Kokkos::parallel_for(config, f);
 #else
@ -293,7 +293,7 @@ int NeighborKokkosExecute<DeviceType>::exclusion(const int &i,const int &j,

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

-template<class DeviceType> template<int HalfNeigh,int Newton>
+template<class DeviceType> template<int HalfNeigh,int Newton,int Tri>
 void NeighborKokkosExecute<DeviceType>::
   build_Item(const int &i) const
 {
@ -365,7 +365,7 @@ void NeighborKokkosExecute<DeviceType>::
    const int jbin = ibin + stencil[k];

    // get subview of jbin
-    if(HalfNeigh&&(ibin==jbin)) continue;
+    if(HalfNeigh && (ibin==jbin)) continue;
    //const ArrayTypes<DeviceType>::t_int_1d_const_um =Kokkos::subview<t_int_1d_const_um>(bins,jbin,ALL);
      for(int m = 0; m < c_bincount(jbin); m++) {

@ -374,6 +374,16 @@ void NeighborKokkosExecute<DeviceType>::

        if(HalfNeigh && !Newton && (j < i)) continue;
        if(!HalfNeigh && j==i) continue;
+        if(Tri) {
+          if (x(j,2) < ztmp) continue;
+          if (x(j,2) == ztmp) {
+            if (x(j,1) < ytmp) continue;
+            if (x(j,1) == ytmp) {
+              if (x(j,0) < xtmp) continue;
+              if (x(j,0) == xtmp && j <= i) continue;
+            }
+          }
+        }
        if(exclude && exclusion(i,j,itype,jtype)) continue;

        const X_FLOAT delx = xtmp - x(j, 0);
@ -428,7 +438,7 @@ extern __shared__ X_FLOAT sharedmem[];

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

-template<class DeviceType> template<int HalfNeigh,int Newton>
+template<class DeviceType> template<int HalfNeigh,int Newton,int Tri>
 __device__ inline
 void NeighborKokkosExecute<DeviceType>::build_ItemCuda(typename Kokkos::TeamPolicy<DeviceType>::member_type dev) const
 {
@ -491,6 +501,16 @@ void NeighborKokkosExecute<DeviceType>::build_ItemCuda(typename Kokkos::TeamPoli
            ((j >= nlocal) && ((x(j, 2) < ztmp) || (x(j, 2) == ztmp && x(j, 1) < ytmp) ||
              (x(j, 2) == ztmp && x(j, 1)  == ytmp && x(j, 0) < xtmp)))))
        ) continue;
+        if(Tri) {
+          if (x(j,2) < ztmp) continue;
+          if (x(j,2) == ztmp) {
+            if (x(j,1) < ytmp) continue;
+            if (x(j,1) == ytmp) {
+              if (x(j,0) < xtmp) continue;
+              if (x(j,0) == xtmp && j <= i) continue;
+            }
+          }
+        }
      if(exclude && exclusion(i,j,itype,jtype)) continue;
      const X_FLOAT delx = xtmp - other_x[m];
      const X_FLOAT dely = ytmp - other_x[m + atoms_per_bin];
@ -558,6 +578,16 @@ void NeighborKokkosExecute<DeviceType>::build_ItemCuda(typename Kokkos::TeamPoli
        //if(HalfNeigh && (j < i))  continue;
        if(HalfNeigh && !Newton && (j < i)) continue;
        if(!HalfNeigh && j==i) continue;
+        if(Tri) {
+          if (x(j,2) < ztmp) continue;
+          if (x(j,2) == ztmp) {
+            if (x(j,1) < ytmp) continue;
+            if (x(j,1) == ytmp) {
+              if (x(j,0) < xtmp) continue;
+              if (x(j,0) == xtmp && j <= i) continue;
+            }
+          }
+        }
        if(exclude && exclusion(i,j,itype,jtype)) continue;

        const X_FLOAT delx = xtmp - other_x[m];
@ -736,14 +766,16 @@ void NeighborKokkosExecute<DeviceType>::
 }

 namespace LAMMPS_NS {
-template class NPairKokkos<LMPDeviceType,0,0>;
-template class NPairKokkos<LMPDeviceType,0,1>;
-template class NPairKokkos<LMPDeviceType,1,0>;
-template class NPairKokkos<LMPDeviceType,1,1>;
+template class NPairKokkos<LMPDeviceType,0,0,0>;
+template class NPairKokkos<LMPDeviceType,0,1,0>;
+template class NPairKokkos<LMPDeviceType,1,0,0>;
+template class NPairKokkos<LMPDeviceType,1,1,0>;
+template class NPairKokkos<LMPDeviceType,1,0,1>;
 #ifdef KOKKOS_HAVE_CUDA
-template class NPairKokkos<LMPHostType,0,0>;
-template class NPairKokkos<LMPHostType,0,1>;
-template class NPairKokkos<LMPHostType,1,0>;
-template class NPairKokkos<LMPHostType,1,1>;
+template class NPairKokkos<LMPHostType,0,0,0>;
+template class NPairKokkos<LMPHostType,0,1,0>;
+template class NPairKokkos<LMPHostType,1,0,0>;
+template class NPairKokkos<LMPHostType,1,1,0>;
+template class NPairKokkos<LMPHostType,1,0,1>;
 #endif
 }
--- a/src/KOKKOS/npair_kokkos.h
+++ b/src/KOKKOS/npair_kokkos.h
@ -13,42 +13,52 @@

 #ifdef NPAIR_CLASS

-typedef NPairKokkos<LMPHostType,0,0> NPairKokkosFullBinHost;
+typedef NPairKokkos<LMPHostType,0,0,0> NPairKokkosFullBinHost;
 NPairStyle(full/bin/kk/host,
           NPairKokkosFullBinHost,
           NP_FULL | NP_BIN | NP_KOKKOS_HOST | NP_NEWTON | NP_NEWTOFF | NP_ORTHO | NP_TRI)

-typedef NPairKokkos<LMPDeviceType,0,0> NPairKokkosFullBinDevice;
+typedef NPairKokkos<LMPDeviceType,0,0,0> NPairKokkosFullBinDevice;
 NPairStyle(full/bin/kk/device,
           NPairKokkosFullBinDevice,
           NP_FULL | NP_BIN | NP_KOKKOS_DEVICE | NP_NEWTON | NP_NEWTOFF | NP_ORTHO | NP_TRI)

-typedef NPairKokkos<LMPHostType,0,1> NPairKokkosFullBinGhostHost;
+typedef NPairKokkos<LMPHostType,0,1,0> NPairKokkosFullBinGhostHost;
 NPairStyle(full/bin/ghost/kk/host,
           NPairKokkosFullBinGhostHost,
           NP_FULL | NP_BIN | NP_KOKKOS_HOST | NP_NEWTON | NP_NEWTOFF | NP_GHOST | NP_ORTHO | NP_TRI)

-typedef NPairKokkos<LMPDeviceType,0,1> NPairKokkosFullBinGhostDevice;
+typedef NPairKokkos<LMPDeviceType,0,1,0> NPairKokkosFullBinGhostDevice;
 NPairStyle(full/bin/ghost/kk/device,
           NPairKokkosFullBinGhostDevice,
           NP_FULL | NP_BIN | NP_KOKKOS_DEVICE | NP_NEWTON | NP_NEWTOFF | NP_GHOST | NP_ORTHO | NP_TRI)

-typedef NPairKokkos<LMPHostType,1,0> NPairKokkosHalfBinHost;
+typedef NPairKokkos<LMPHostType,1,0,0> NPairKokkosHalfBinHost;
 NPairStyle(half/bin/kk/host,
           NPairKokkosHalfBinHost,
-           NP_HALF | NP_BIN | NP_KOKKOS_HOST | NP_NEWTON | NP_NEWTOFF | NP_ORTHO | NP_TRI)
+           NP_HALF | NP_BIN | NP_KOKKOS_HOST | NP_NEWTON | NP_NEWTOFF | NP_ORTHO)

-typedef NPairKokkos<LMPDeviceType,1,0> NPairKokkosHalfBinDevice;
+typedef NPairKokkos<LMPDeviceType,1,0,0> NPairKokkosHalfBinDevice;
 NPairStyle(half/bin/kk/device,
           NPairKokkosHalfBinDevice,
-           NP_HALF | NP_BIN | NP_KOKKOS_DEVICE | NP_NEWTON | NP_NEWTOFF | NP_ORTHO | NP_TRI)
+           NP_HALF | NP_BIN | NP_KOKKOS_DEVICE | NP_NEWTON | NP_NEWTOFF | NP_ORTHO)

-typedef NPairKokkos<LMPHostType,1,1> NPairKokkosHalfBinGhostHost;
+typedef NPairKokkos<LMPHostType,1,0,1> NPairKokkosHalfBinHostTri;
+NPairStyle(half/bin/kk/host,
+           NPairKokkosHalfBinHostTri,
+           NP_HALF | NP_BIN | NP_KOKKOS_HOST | NP_NEWTON | NP_NEWTOFF | NP_TRI)
+
+typedef NPairKokkos<LMPDeviceType,1,0,1> NPairKokkosHalfBinDeviceTri;
+NPairStyle(half/bin/kk/device,
+           NPairKokkosHalfBinDeviceTri,
+           NP_HALF | NP_BIN | NP_KOKKOS_DEVICE | NP_NEWTON | NP_NEWTOFF | NP_TRI)
+
+typedef NPairKokkos<LMPHostType,1,1,0> NPairKokkosHalfBinGhostHost;
 NPairStyle(half/bin/ghost/kk/host,
           NPairKokkosHalfBinGhostHost,
           NP_HALF | NP_BIN | NP_KOKKOS_HOST | NP_NEWTON | NP_NEWTOFF | NP_GHOST | NP_ORTHO | NP_TRI)

-typedef NPairKokkos<LMPDeviceType,1,1> NPairKokkosHalfBinGhostDevice;
+typedef NPairKokkos<LMPDeviceType,1,1,0> NPairKokkosHalfBinGhostDevice;
 NPairStyle(half/bin/ghost/kk/device,
           NPairKokkosHalfBinGhostDevice,
           NP_HALF | NP_BIN | NP_KOKKOS_DEVICE | NP_NEWTON | NP_NEWTOFF | NP_GHOST | NP_ORTHO | NP_TRI)
@ -63,7 +73,7 @@ NPairStyle(half/bin/ghost/kk/device,

 namespace LAMMPS_NS {

-template<class DeviceType, int HALF_NEIGH, int GHOST>
+template<class DeviceType, int HALF_NEIGH, int GHOST, int TRI>
 class NPairKokkos : public NPair {
 public:
  NPairKokkos(class LAMMPS *);
@ -147,7 +157,8 @@ class NeighborKokkosExecute
  // data from Atom class

  const typename AT::t_x_array_randomread x;
-  const typename AT::t_int_1d_const type,mask,molecule;
+  const typename AT::t_int_1d_const type,mask;
+  const typename AT::t_tagint_1d_const molecule;
  const typename AT::t_tagint_1d_const tag;
  const typename AT::t_tagint_2d_const special;
  const typename AT::t_int_2d_const nspecial;
@ -184,7 +195,7 @@ class NeighborKokkosExecute
                        const typename AT::t_x_array_randomread &_x,
                        const typename AT::t_int_1d_const &_type,
                        const typename AT::t_int_1d_const &_mask,
-                        const typename AT::t_int_1d_const &_molecule,
+                        const typename AT::t_tagint_1d_const &_molecule,
                        const typename AT::t_tagint_1d_const &_tag,
                        const typename AT::t_tagint_2d_const &_special,
                        const typename AT::t_int_2d_const &_nspecial,
@ -252,7 +263,7 @@ class NeighborKokkosExecute

  ~NeighborKokkosExecute() {neigh_list.clean_copy();};

-  template<int HalfNeigh, int Newton>
+  template<int HalfNeigh, int Newton, int Tri>
  KOKKOS_FUNCTION
  void build_Item(const int &i) const;

@ -261,7 +272,7 @@ class NeighborKokkosExecute
  void build_Item_Ghost(const int &i) const;

 #ifdef KOKKOS_HAVE_CUDA
-  template<int HalfNeigh, int Newton>
+  template<int HalfNeigh, int Newton, int Tri>
  __device__ inline
  void build_ItemCuda(typename Kokkos::TeamPolicy<DeviceType>::member_type dev) const;
 #endif
@ -353,7 +364,7 @@ class NeighborKokkosExecute

 };

-template<class DeviceType,int HALF_NEIGH,int GHOST_NEWTON>
+template<class DeviceType, int HALF_NEIGH, int GHOST_NEWTON, int TRI>
 struct NPairKokkosBuildFunctor {
  typedef DeviceType device_type;

@ -366,20 +377,20 @@ struct NPairKokkosBuildFunctor {

  KOKKOS_INLINE_FUNCTION
  void operator() (const int & i) const {
-    c.template build_Item<HALF_NEIGH,GHOST_NEWTON>(i);
+    c.template build_Item<HALF_NEIGH,GHOST_NEWTON,TRI>(i);
  }
 #ifdef KOKKOS_HAVE_CUDA
  __device__ inline
  
  void operator() (typename Kokkos::TeamPolicy<DeviceType>::member_type dev) const {
-    c.template build_ItemCuda<HALF_NEIGH,GHOST_NEWTON>(dev);
+    c.template build_ItemCuda<HALF_NEIGH,GHOST_NEWTON,TRI>(dev);
  }
  size_t shmem_size(const int team_size) const { (void) team_size; return sharedsize; }
 #endif
 };

-template<int HALF_NEIGH,int GHOST_NEWTON>
-struct NPairKokkosBuildFunctor<LMPHostType,HALF_NEIGH,GHOST_NEWTON> {
+template<int HALF_NEIGH, int GHOST_NEWTON, int TRI>
+struct NPairKokkosBuildFunctor<LMPHostType,HALF_NEIGH,GHOST_NEWTON,TRI> {
  typedef LMPHostType device_type;

  const NeighborKokkosExecute<LMPHostType> c;
@ -391,7 +402,7 @@ struct NPairKokkosBuildFunctor<LMPHostType,HALF_NEIGH,GHOST_NEWTON> {

  KOKKOS_INLINE_FUNCTION
  void operator() (const int & i) const {
-    c.template build_Item<HALF_NEIGH,GHOST_NEWTON>(i);
+    c.template build_Item<HALF_NEIGH,GHOST_NEWTON,TRI>(i);
  }

  void operator() (typename Kokkos::TeamPolicy<LMPHostType>::member_type dev) const {}
--- a/src/KOKKOS/pair_buck_coul_cut_kokkos.cpp
+++ b/src/KOKKOS/pair_buck_coul_cut_kokkos.cpp
@ -154,7 +154,7 @@ void PairBuckCoulCutKokkos<DeviceType>::compute(int eflag_in, int vflag_in)
    virial[5] += ev.v[5];
  }

-  if (vflag_fdotr) virial_fdotr_compute();
+  if (vflag_fdotr) pair_virial_fdotr_compute(this);

  if (eflag_atom) {
    k_eatom.template modify<DeviceType>();
--- a/src/KOKKOS/pair_buck_coul_cut_kokkos.h
+++ b/src/KOKKOS/pair_buck_coul_cut_kokkos.h
@ -34,6 +34,7 @@ class PairBuckCoulCutKokkos : public PairBuckCoulCut {
  enum {EnabledNeighFlags=FULL|HALFTHREAD|HALF};
  enum {COUL_FLAG=1};
  typedef DeviceType device_type;
+  typedef ArrayTypes<DeviceType> AT;
  PairBuckCoulCutKokkos(class LAMMPS *);
  ~PairBuckCoulCutKokkos();

@ -44,7 +45,9 @@ class PairBuckCoulCutKokkos : public PairBuckCoulCut {
  double init_one(int, int);

  struct params_buck_coul{
+    KOKKOS_INLINE_FUNCTION
    params_buck_coul(){cut_ljsq=0;cut_coulsq=0;a=0;c=0;rhoinv=0;buck1=0;buck2=0;offset=0;};
+    KOKKOS_INLINE_FUNCTION
    params_buck_coul(int i){cut_ljsq=0;cut_coulsq=0;a=0;c=0;rhoinv=0;buck1=0;buck2=0;offset=0;};
    F_FLOAT cut_ljsq,cut_coulsq,a,c,rhoinv,buck1,buck2,offset;
  };
@ -81,25 +84,25 @@ class PairBuckCoulCutKokkos : public PairBuckCoulCut {
  F_FLOAT m_cutsq[MAX_TYPES_STACKPARAMS+1][MAX_TYPES_STACKPARAMS+1];
  F_FLOAT m_cut_ljsq[MAX_TYPES_STACKPARAMS+1][MAX_TYPES_STACKPARAMS+1];
  F_FLOAT m_cut_coulsq[MAX_TYPES_STACKPARAMS+1][MAX_TYPES_STACKPARAMS+1];
-  typename ArrayTypes<DeviceType>::t_x_array_randomread x;
-  typename ArrayTypes<DeviceType>::t_x_array c_x;
-  typename ArrayTypes<DeviceType>::t_f_array f;
-  typename ArrayTypes<DeviceType>::t_int_1d_randomread type;
-  typename ArrayTypes<DeviceType>::t_float_1d_randomread q;
+  typename AT::t_x_array_randomread x;
+  typename AT::t_x_array c_x;
+  typename AT::t_f_array f;
+  typename AT::t_int_1d_randomread type;
+  typename AT::t_float_1d_randomread q;

  DAT::tdual_efloat_1d k_eatom;
  DAT::tdual_virial_array k_vatom;
-  typename ArrayTypes<DeviceType>::t_efloat_1d d_eatom;
-  typename ArrayTypes<DeviceType>::t_virial_array d_vatom;
+  typename AT::t_efloat_1d d_eatom;
+  typename AT::t_virial_array d_vatom;

  int newton_pair;

-  typename ArrayTypes<DeviceType>::tdual_ffloat_2d k_cutsq;
-  typename ArrayTypes<DeviceType>::t_ffloat_2d d_cutsq;
-  typename ArrayTypes<DeviceType>::tdual_ffloat_2d k_cut_ljsq;
-  typename ArrayTypes<DeviceType>::t_ffloat_2d d_cut_ljsq;
-  typename ArrayTypes<DeviceType>::tdual_ffloat_2d k_cut_coulsq;
-  typename ArrayTypes<DeviceType>::t_ffloat_2d d_cut_coulsq;
+  typename AT::tdual_ffloat_2d k_cutsq;
+  typename AT::t_ffloat_2d d_cutsq;
+  typename AT::tdual_ffloat_2d k_cut_ljsq;
+  typename AT::t_ffloat_2d d_cut_ljsq;
+  typename AT::tdual_ffloat_2d k_cut_coulsq;
+  typename AT::t_ffloat_2d d_cut_coulsq;


  int neighflag;
@ -121,6 +124,7 @@ class PairBuckCoulCutKokkos : public PairBuckCoulCut {
  friend EV_FLOAT pair_compute_neighlist<PairBuckCoulCutKokkos,HALFTHREAD,void>(PairBuckCoulCutKokkos*,NeighListKokkos<DeviceType>*);
  friend EV_FLOAT pair_compute<PairBuckCoulCutKokkos,void>(PairBuckCoulCutKokkos*,
                                                            NeighListKokkos<DeviceType>*);
+  friend void pair_virial_fdotr_compute<PairBuckCoulCutKokkos>(PairBuckCoulCutKokkos*);

 };

--- a/src/KOKKOS/pair_buck_coul_long_kokkos.cpp
+++ b/src/KOKKOS/pair_buck_coul_long_kokkos.cpp
@ -114,6 +114,19 @@ void PairBuckCoulLongKokkos<DeviceType>::compute(int eflag_in, int vflag_in)
  if (eflag || vflag) ev_setup(eflag,vflag);
  else evflag = vflag_fdotr = 0;

+  // reallocate per-atom arrays if necessary
+
+  if (eflag_atom) {
+    memory->destroy_kokkos(k_eatom,eatom);
+    memory->create_kokkos(k_eatom,eatom,maxeatom,"pair:eatom");
+    d_eatom = k_eatom.view<DeviceType>();
+  }
+  if (vflag_atom) {
+    memory->destroy_kokkos(k_vatom,vatom);
+    memory->create_kokkos(k_vatom,vatom,maxvatom,6,"pair:vatom");
+    d_vatom = k_vatom.view<DeviceType>();
+  }
+
  atomKK->sync(execution_space,datamask_read);
  k_cutsq.template sync<DeviceType>();
  k_cut_ljsq.template sync<DeviceType>();
--- a/src/KOKKOS/pair_buck_coul_long_kokkos.h
+++ b/src/KOKKOS/pair_buck_coul_long_kokkos.h
@ -34,6 +34,7 @@ class PairBuckCoulLongKokkos : public PairBuckCoulLong {
  enum {EnabledNeighFlags=FULL|HALFTHREAD|HALF};
  enum {COUL_FLAG=1};
  typedef DeviceType device_type;
+  typedef ArrayTypes<DeviceType> AT;
  PairBuckCoulLongKokkos(class LAMMPS *);
  ~PairBuckCoulLongKokkos();

@ -45,7 +46,9 @@ class PairBuckCoulLongKokkos : public PairBuckCoulLong {
  double init_one(int, int);

  struct params_buck_coul{
+    KOKKOS_INLINE_FUNCTION
    params_buck_coul(){cut_ljsq=0;cut_coulsq=0;a=0;c=0;rhoinv=0;buck1=0;buck2=0;offset=0;};
+    KOKKOS_INLINE_FUNCTION
    params_buck_coul(int i){cut_ljsq=0;cut_coulsq=0;a=0;c=0;rhoinv=0;buck1=0;buck2=0;offset=0;};
    F_FLOAT cut_ljsq,cut_coulsq,a,c,rhoinv,buck1,buck2,offset;
  };
@ -82,27 +85,27 @@ class PairBuckCoulLongKokkos : public PairBuckCoulLong {
  F_FLOAT m_cutsq[MAX_TYPES_STACKPARAMS+1][MAX_TYPES_STACKPARAMS+1];
  F_FLOAT m_cut_ljsq[MAX_TYPES_STACKPARAMS+1][MAX_TYPES_STACKPARAMS+1];
  F_FLOAT m_cut_coulsq[MAX_TYPES_STACKPARAMS+1][MAX_TYPES_STACKPARAMS+1];
-  typename ArrayTypes<DeviceType>::t_x_array_randomread x;
-  typename ArrayTypes<DeviceType>::t_x_array c_x;
-  typename ArrayTypes<DeviceType>::t_f_array f;
-  typename ArrayTypes<DeviceType>::t_int_1d_randomread type;
-  typename ArrayTypes<DeviceType>::t_float_1d_randomread q;
+  typename AT::t_x_array_randomread x;
+  typename AT::t_x_array c_x;
+  typename AT::t_f_array f;
+  typename AT::t_int_1d_randomread type;
+  typename AT::t_float_1d_randomread q;

  DAT::tdual_efloat_1d k_eatom;
  DAT::tdual_virial_array k_vatom;
-  typename ArrayTypes<DeviceType>::t_efloat_1d d_eatom;
-  typename ArrayTypes<DeviceType>::t_virial_array d_vatom;
+  typename AT::t_efloat_1d d_eatom;
+  typename AT::t_virial_array d_vatom;

  int newton_pair;

-  typename ArrayTypes<DeviceType>::tdual_ffloat_2d k_cutsq;
-  typename ArrayTypes<DeviceType>::t_ffloat_2d d_cutsq;
-  typename ArrayTypes<DeviceType>::tdual_ffloat_2d k_cut_ljsq;
-  typename ArrayTypes<DeviceType>::t_ffloat_2d d_cut_ljsq;
-  typename ArrayTypes<DeviceType>::tdual_ffloat_2d k_cut_coulsq;
-  typename ArrayTypes<DeviceType>::t_ffloat_2d d_cut_coulsq;
+  typename AT::tdual_ffloat_2d k_cutsq;
+  typename AT::t_ffloat_2d d_cutsq;
+  typename AT::tdual_ffloat_2d k_cut_ljsq;
+  typename AT::t_ffloat_2d d_cut_ljsq;
+  typename AT::tdual_ffloat_2d k_cut_coulsq;
+  typename AT::t_ffloat_2d d_cut_coulsq;

-  typename ArrayTypes<DeviceType>::t_ffloat_1d_randomread
+  typename AT::t_ffloat_1d_randomread
    d_rtable, d_drtable, d_ftable, d_dftable,
    d_ctable, d_dctable, d_etable, d_detable;

--- a/src/KOKKOS/pair_buck_kokkos.h
+++ b/src/KOKKOS/pair_buck_kokkos.h
@ -34,6 +34,7 @@ class PairBuckKokkos : public PairBuck {
  enum {EnabledNeighFlags=FULL|HALFTHREAD|HALF|N2};
  enum {COUL_FLAG=0};
  typedef DeviceType device_type;
+  typedef ArrayTypes<DeviceType> AT;
  PairBuckKokkos(class LAMMPS *);
  ~PairBuckKokkos();

@ -43,7 +44,9 @@ class PairBuckKokkos : public PairBuck {
  double init_one(int, int);

  struct params_buck{
+    KOKKOS_INLINE_FUNCTION
    params_buck(){cutsq=0;a=0;c=0;rhoinv=0;buck1=0;buck2=0;offset=0;};
+    KOKKOS_INLINE_FUNCTION
    params_buck(int i){cutsq=0;a=0;c=0;rhoinv=0;buck1=0;buck2=0;offset=0;};
    F_FLOAT cutsq,a,c,rhoinv,buck1,buck2,offset;
  };
@ -70,22 +73,22 @@ class PairBuckKokkos : public PairBuck {
  typename Kokkos::DualView<params_buck**,Kokkos::LayoutRight,DeviceType>::t_dev_const_um params;
  params_buck m_params[MAX_TYPES_STACKPARAMS+1][MAX_TYPES_STACKPARAMS+1];  // hardwired to space for 12 atom types
  F_FLOAT m_cutsq[MAX_TYPES_STACKPARAMS+1][MAX_TYPES_STACKPARAMS+1];
-  typename ArrayTypes<DeviceType>::t_x_array_randomread x;
-  typename ArrayTypes<DeviceType>::t_x_array c_x;
-  typename ArrayTypes<DeviceType>::t_f_array f;
-  typename ArrayTypes<DeviceType>::t_int_1d_randomread type;
-  typename ArrayTypes<DeviceType>::t_tagint_1d tag;
+  typename AT::t_x_array_randomread x;
+  typename AT::t_x_array c_x;
+  typename AT::t_f_array f;
+  typename AT::t_int_1d_randomread type;
+  typename AT::t_tagint_1d tag;

  DAT::tdual_efloat_1d k_eatom;
  DAT::tdual_virial_array k_vatom;
-  typename ArrayTypes<DeviceType>::t_efloat_1d d_eatom;
-  typename ArrayTypes<DeviceType>::t_virial_array d_vatom;
+  typename AT::t_efloat_1d d_eatom;
+  typename AT::t_virial_array d_vatom;

  int newton_pair;
  double special_lj[4];

-  typename ArrayTypes<DeviceType>::tdual_ffloat_2d k_cutsq;
-  typename ArrayTypes<DeviceType>::t_ffloat_2d d_cutsq;
+  typename AT::tdual_ffloat_2d k_cutsq;
+  typename AT::t_ffloat_2d d_cutsq;


  int neighflag;
--- a/src/KOKKOS/pair_coul_cut_kokkos.cpp
+++ b/src/KOKKOS/pair_coul_cut_kokkos.cpp
@ -83,6 +83,19 @@ void PairCoulCutKokkos<DeviceType>::compute(int eflag_in, int vflag_in)
  if (eflag || vflag) ev_setup(eflag,vflag);
  else evflag = vflag_fdotr = 0;

+  // reallocate per-atom arrays if necessary
+
+  if (eflag_atom) {
+    memory->destroy_kokkos(k_eatom,eatom);
+    memory->create_kokkos(k_eatom,eatom,maxeatom,"pair:eatom");
+    d_eatom = k_eatom.view<DeviceType>();
+  }
+  if (vflag_atom) {
+    memory->destroy_kokkos(k_vatom,vatom);
+    memory->create_kokkos(k_vatom,vatom,maxvatom,6,"pair:vatom");
+    d_vatom = k_vatom.view<DeviceType>();
+  }
+
  atomKK->sync(execution_space,datamask_read);
  k_cutsq.template sync<DeviceType>();
  k_cut_ljsq.template sync<DeviceType>();
@ -124,6 +137,16 @@ void PairCoulCutKokkos<DeviceType>::compute(int eflag_in, int vflag_in)
    virial[5] += ev.v[5];
  }

+  if (eflag_atom) {
+    k_eatom.template modify<DeviceType>();
+    k_eatom.template sync<LMPHostType>();
+  }
+
+  if (vflag_atom) {
+    k_vatom.template modify<DeviceType>();
+    k_vatom.template sync<LMPHostType>();
+  }
+
  if (vflag_fdotr) pair_virial_fdotr_compute(this);
 }

--- a/src/KOKKOS/pair_coul_cut_kokkos.h
+++ b/src/KOKKOS/pair_coul_cut_kokkos.h
@ -34,6 +34,7 @@ class PairCoulCutKokkos : public PairCoulCut {
  enum {EnabledNeighFlags=FULL|HALFTHREAD|HALF};
  enum {COUL_FLAG=1};
  typedef DeviceType device_type;
+  typedef ArrayTypes<DeviceType> AT;
  PairCoulCutKokkos(class LAMMPS *);
  ~PairCoulCutKokkos();

@ -87,22 +88,25 @@ class PairCoulCutKokkos : public PairCoulCut {
  F_FLOAT m_cutsq[MAX_TYPES_STACKPARAMS+1][MAX_TYPES_STACKPARAMS+1];
  F_FLOAT m_cut_ljsq[MAX_TYPES_STACKPARAMS+1][MAX_TYPES_STACKPARAMS+1];
  F_FLOAT m_cut_coulsq[MAX_TYPES_STACKPARAMS+1][MAX_TYPES_STACKPARAMS+1];
-  typename ArrayTypes<DeviceType>::t_x_array_randomread x;
-  typename ArrayTypes<DeviceType>::t_x_array c_x;
-  typename ArrayTypes<DeviceType>::t_f_array f;
-  typename ArrayTypes<DeviceType>::t_float_1d_randomread q;
-  typename ArrayTypes<DeviceType>::t_int_1d_randomread type;
-  typename ArrayTypes<DeviceType>::t_efloat_1d d_eatom;
-  typename ArrayTypes<DeviceType>::t_virial_array d_vatom;
+  typename AT::t_x_array_randomread x;
+  typename AT::t_x_array c_x;
+  typename AT::t_f_array f;
+  typename AT::t_float_1d_randomread q;
+  typename AT::t_int_1d_randomread type;
+
+  DAT::tdual_efloat_1d k_eatom;
+  DAT::tdual_virial_array k_vatom;
+  typename AT::t_efloat_1d d_eatom;
+  typename AT::t_virial_array d_vatom;

  int newton_pair;

-  typename ArrayTypes<DeviceType>::tdual_ffloat_2d k_cutsq;
-  typename ArrayTypes<DeviceType>::t_ffloat_2d d_cutsq;
-  typename ArrayTypes<DeviceType>::tdual_ffloat_2d k_cut_ljsq;
-  typename ArrayTypes<DeviceType>::t_ffloat_2d d_cut_ljsq;
-  typename ArrayTypes<DeviceType>::tdual_ffloat_2d k_cut_coulsq;
-  typename ArrayTypes<DeviceType>::t_ffloat_2d d_cut_coulsq;
+  typename AT::tdual_ffloat_2d k_cutsq;
+  typename AT::t_ffloat_2d d_cutsq;
+  typename AT::tdual_ffloat_2d k_cut_ljsq;
+  typename AT::t_ffloat_2d d_cut_ljsq;
+  typename AT::tdual_ffloat_2d k_cut_coulsq;
+  typename AT::t_ffloat_2d d_cut_coulsq;


  int neighflag;
--- a/src/KOKKOS/pair_coul_debye_kokkos.cpp
+++ b/src/KOKKOS/pair_coul_debye_kokkos.cpp
@ -90,6 +90,19 @@ void PairCoulDebyeKokkos<DeviceType>::compute(int eflag_in, int vflag_in)
  if (eflag || vflag) ev_setup(eflag,vflag);
  else evflag = vflag_fdotr = 0;

+  // reallocate per-atom arrays if necessary
+
+  if (eflag_atom) {
+    memory->destroy_kokkos(k_eatom,eatom);
+    memory->create_kokkos(k_eatom,eatom,maxeatom,"pair:eatom");
+    d_eatom = k_eatom.view<DeviceType>();
+  }
+  if (vflag_atom) {
+    memory->destroy_kokkos(k_vatom,vatom);
+    memory->create_kokkos(k_vatom,vatom,maxvatom,6,"pair:vatom");
+    d_vatom = k_vatom.view<DeviceType>();
+  }
+
  atomKK->sync(execution_space,datamask_read);
  k_cutsq.template sync<DeviceType>();
  k_cut_ljsq.template sync<DeviceType>();
@ -136,7 +149,17 @@ void PairCoulDebyeKokkos<DeviceType>::compute(int eflag_in, int vflag_in)
    virial[5] += ev.v[5];
  }

-  if (vflag_fdotr) virial_fdotr_compute();
+  if (eflag_atom) {
+    k_eatom.template modify<DeviceType>();
+    k_eatom.template sync<LMPHostType>();
+  }
+
+  if (vflag_atom) {
+    k_vatom.template modify<DeviceType>();
+    k_vatom.template sync<LMPHostType>();
+  }
+
+  if (vflag_fdotr) pair_virial_fdotr_compute(this);

  copymode = 0;
 }
--- a/src/KOKKOS/pair_coul_debye_kokkos.h
+++ b/src/KOKKOS/pair_coul_debye_kokkos.h
@ -34,6 +34,7 @@ class PairCoulDebyeKokkos : public PairCoulDebye {
  enum {EnabledNeighFlags=FULL|HALFTHREAD|HALF};
  enum {COUL_FLAG=1};
  typedef DeviceType device_type;
+  typedef ArrayTypes<DeviceType> AT;
  PairCoulDebyeKokkos(class LAMMPS *);
  ~PairCoulDebyeKokkos();

@ -44,7 +45,9 @@ class PairCoulDebyeKokkos : public PairCoulDebye {
  double init_one(int, int);

  struct params_coul{
+    KOKKOS_INLINE_FUNCTION
    params_coul(){cutsq=0,scale=0;};
+    KOKKOS_INLINE_FUNCTION
    params_coul(int i){cutsq=0,scale=0;};
    F_FLOAT cutsq, scale;
  };
@ -85,22 +88,25 @@ class PairCoulDebyeKokkos : public PairCoulDebye {
  F_FLOAT m_cutsq[MAX_TYPES_STACKPARAMS+1][MAX_TYPES_STACKPARAMS+1];
  F_FLOAT m_cut_ljsq[MAX_TYPES_STACKPARAMS+1][MAX_TYPES_STACKPARAMS+1];
  F_FLOAT m_cut_coulsq[MAX_TYPES_STACKPARAMS+1][MAX_TYPES_STACKPARAMS+1];
-  typename ArrayTypes<DeviceType>::t_x_array_randomread x;
-  typename ArrayTypes<DeviceType>::t_x_array c_x;
-  typename ArrayTypes<DeviceType>::t_f_array f;
-  typename ArrayTypes<DeviceType>::t_int_1d_randomread type;
-  typename ArrayTypes<DeviceType>::t_float_1d_randomread q;
-  typename ArrayTypes<DeviceType>::t_efloat_1d d_eatom;
-  typename ArrayTypes<DeviceType>::t_virial_array d_vatom;
+  typename AT::t_x_array_randomread x;
+  typename AT::t_x_array c_x;
+  typename AT::t_f_array f;
+  typename AT::t_int_1d_randomread type;
+  typename AT::t_float_1d_randomread q;
+
+  DAT::tdual_efloat_1d k_eatom;
+  DAT::tdual_virial_array k_vatom;
+  typename AT::t_efloat_1d d_eatom;
+  typename AT::t_virial_array d_vatom;

  int newton_pair;

-  typename ArrayTypes<DeviceType>::tdual_ffloat_2d k_cutsq;
-  typename ArrayTypes<DeviceType>::t_ffloat_2d d_cutsq;
-  typename ArrayTypes<DeviceType>::tdual_ffloat_2d k_cut_ljsq;
-  typename ArrayTypes<DeviceType>::t_ffloat_2d d_cut_ljsq;
-  typename ArrayTypes<DeviceType>::tdual_ffloat_2d k_cut_coulsq;
-  typename ArrayTypes<DeviceType>::t_ffloat_2d d_cut_coulsq;
+  typename AT::tdual_ffloat_2d k_cutsq;
+  typename AT::t_ffloat_2d d_cutsq;
+  typename AT::tdual_ffloat_2d k_cut_ljsq;
+  typename AT::t_ffloat_2d d_cut_ljsq;
+  typename AT::tdual_ffloat_2d k_cut_coulsq;
+  typename AT::t_ffloat_2d d_cut_coulsq;


  int neighflag;
@ -122,6 +128,7 @@ class PairCoulDebyeKokkos : public PairCoulDebye {
  friend EV_FLOAT pair_compute_neighlist<PairCoulDebyeKokkos,HALFTHREAD,void>(PairCoulDebyeKokkos*,NeighListKokkos<DeviceType>*);
  friend EV_FLOAT pair_compute<PairCoulDebyeKokkos,void>(PairCoulDebyeKokkos*,
                                                            NeighListKokkos<DeviceType>*);
+  friend void pair_virial_fdotr_compute<PairCoulDebyeKokkos>(PairCoulDebyeKokkos*);

 };

--- a/src/KOKKOS/pair_coul_dsf_kokkos.cpp
+++ b/src/KOKKOS/pair_coul_dsf_kokkos.cpp
@ -88,12 +88,12 @@ void PairCoulDSFKokkos<DeviceType>::compute(int eflag_in, int vflag_in)
  if (eflag_atom) {
    memory->destroy_kokkos(k_eatom,eatom);
    memory->create_kokkos(k_eatom,eatom,maxeatom,"pair:eatom");
-    d_eatom = k_eatom.d_view;
+    d_eatom = k_eatom.view<DeviceType>();
  }
  if (vflag_atom) {
    memory->destroy_kokkos(k_vatom,vatom);
    memory->create_kokkos(k_vatom,vatom,maxvatom,6,"pair:vatom");
-    d_vatom = k_vatom.d_view;
+    d_vatom = k_vatom.view<DeviceType>();
  }

  atomKK->sync(execution_space,datamask_read);
@ -183,8 +183,6 @@ void PairCoulDSFKokkos<DeviceType>::compute(int eflag_in, int vflag_in)
    virial[5] += ev.v[5];
  }

-  if (vflag_fdotr) pair_virial_fdotr_compute(this);
-
  if (eflag_atom) {
    k_eatom.template modify<DeviceType>();
    k_eatom.template sync<LMPHostType>();
@ -195,6 +193,8 @@ void PairCoulDSFKokkos<DeviceType>::compute(int eflag_in, int vflag_in)
    k_vatom.template sync<LMPHostType>();
  }

+  if (vflag_fdotr) pair_virial_fdotr_compute(this);
+
  copymode = 0;
 }

--- a/src/KOKKOS/pair_coul_dsf_kokkos.h
+++ b/src/KOKKOS/pair_coul_dsf_kokkos.h
@ -37,6 +37,7 @@ class PairCoulDSFKokkos : public PairCoulDSF {
  enum {EnabledNeighFlags=FULL|HALFTHREAD|HALF};
  enum {COUL_FLAG=1};
  typedef DeviceType device_type;
+  typedef ArrayTypes<DeviceType> AT;
  typedef EV_FLOAT value_type;
  PairCoulDSFKokkos(class LAMMPS *);
  ~PairCoulDSFKokkos();
@ -61,16 +62,16 @@ class PairCoulDSFKokkos : public PairCoulDSF {
  KOKKOS_INLINE_FUNCTION
  int sbmask(const int& j) const;

- protected:
+  typename AT::t_x_array_randomread x;
+  typename AT::t_f_array f;
+  typename AT::t_float_1d_randomread q;

-  typename ArrayTypes<DeviceType>::t_x_array_randomread x;
-  typename ArrayTypes<DeviceType>::t_f_array f;
-  typename ArrayTypes<DeviceType>::t_float_1d_randomread q;
+ protected:

  DAT::tdual_efloat_1d k_eatom;
  DAT::tdual_virial_array k_vatom;
-  DAT::t_efloat_1d d_eatom;
-  DAT::t_virial_array d_vatom;
+  typename AT::t_efloat_1d d_eatom;
+  typename AT::t_virial_array d_vatom;


  int neighflag,newton_pair;
@ -79,9 +80,9 @@ class PairCoulDSFKokkos : public PairCoulDSF {
  double special_coul[4];
  double qqrd2e;

-  typename ArrayTypes<DeviceType>::t_neighbors_2d d_neighbors;
-  typename ArrayTypes<DeviceType>::t_int_1d_randomread d_ilist;
-  typename ArrayTypes<DeviceType>::t_int_1d_randomread d_numneigh;
+  typename AT::t_neighbors_2d d_neighbors;
+  typename AT::t_int_1d_randomread d_ilist;
+  typename AT::t_int_1d_randomread d_numneigh;
  //NeighListKokkos<DeviceType> k_list;

  friend void pair_virial_fdotr_compute<PairCoulDSFKokkos>(PairCoulDSFKokkos*);
--- a/src/KOKKOS/pair_coul_long_kokkos.cpp
+++ b/src/KOKKOS/pair_coul_long_kokkos.cpp
@ -107,6 +107,19 @@ void PairCoulLongKokkos<DeviceType>::compute(int eflag_in, int vflag_in)
  if (eflag || vflag) ev_setup(eflag,vflag);
  else evflag = vflag_fdotr = 0;

+  // reallocate per-atom arrays if necessary
+
+  if (eflag_atom) {
+    memory->destroy_kokkos(k_eatom,eatom);
+    memory->create_kokkos(k_eatom,eatom,maxeatom,"pair:eatom");
+    d_eatom = k_eatom.view<DeviceType>();
+  }
+  if (vflag_atom) {
+    memory->destroy_kokkos(k_vatom,vatom);
+    memory->create_kokkos(k_vatom,vatom,maxvatom,6,"pair:vatom");
+    d_vatom = k_vatom.view<DeviceType>();
+  }
+
  atomKK->sync(execution_space,datamask_read);
  k_cutsq.template sync<DeviceType>();
  k_cut_coulsq.template sync<DeviceType>();
@ -158,6 +171,16 @@ void PairCoulLongKokkos<DeviceType>::compute(int eflag_in, int vflag_in)
    virial[5] += ev.v[5];
  }

+  if (eflag_atom) {
+    k_eatom.template modify<DeviceType>();
+    k_eatom.template sync<LMPHostType>();
+  }
+
+  if (vflag_atom) {
+    k_vatom.template modify<DeviceType>();
+    k_vatom.template sync<LMPHostType>();
+  }
+
  if (vflag_fdotr) pair_virial_fdotr_compute(this);

  copymode = 0;
--- a/src/KOKKOS/pair_coul_long_kokkos.h
+++ b/src/KOKKOS/pair_coul_long_kokkos.h
@ -34,6 +34,7 @@ class PairCoulLongKokkos : public PairCoulLong {
  enum {EnabledNeighFlags=FULL|HALFTHREAD|HALF};
  enum {COUL_FLAG=1};
  typedef DeviceType device_type;
+  typedef ArrayTypes<DeviceType> AT;
  PairCoulLongKokkos(class LAMMPS *);
  ~PairCoulLongKokkos();

@ -45,7 +46,9 @@ class PairCoulLongKokkos : public PairCoulLong {
  double init_one(int, int);

  struct params_coul{
+    KOKKOS_INLINE_FUNCTION
    params_coul(){cut_coulsq=0;};
+    KOKKOS_INLINE_FUNCTION
    params_coul(int i){cut_coulsq=0;};
    F_FLOAT cut_coulsq;
  };
@ -86,27 +89,27 @@ class PairCoulLongKokkos : public PairCoulLong {
  F_FLOAT m_cutsq[MAX_TYPES_STACKPARAMS+1][MAX_TYPES_STACKPARAMS+1];
  F_FLOAT m_cut_ljsq[MAX_TYPES_STACKPARAMS+1][MAX_TYPES_STACKPARAMS+1];
  F_FLOAT m_cut_coulsq[MAX_TYPES_STACKPARAMS+1][MAX_TYPES_STACKPARAMS+1];
-  typename ArrayTypes<DeviceType>::t_x_array_randomread x;
-  typename ArrayTypes<DeviceType>::t_x_array c_x;
-  typename ArrayTypes<DeviceType>::t_f_array f;
-  typename ArrayTypes<DeviceType>::t_int_1d_randomread type;
-  typename ArrayTypes<DeviceType>::t_float_1d_randomread q;
+  typename AT::t_x_array_randomread x;
+  typename AT::t_x_array c_x;
+  typename AT::t_f_array f;
+  typename AT::t_int_1d_randomread type;
+  typename AT::t_float_1d_randomread q;

  DAT::tdual_efloat_1d k_eatom;
  DAT::tdual_virial_array k_vatom;
-  typename ArrayTypes<DeviceType>::t_efloat_1d d_eatom;
-  typename ArrayTypes<DeviceType>::t_virial_array d_vatom;
+  typename AT::t_efloat_1d d_eatom;
+  typename AT::t_virial_array d_vatom;

  int newton_pair;

-  typename ArrayTypes<DeviceType>::tdual_ffloat_2d k_cutsq;
-  typename ArrayTypes<DeviceType>::t_ffloat_2d d_cutsq;
-  typename ArrayTypes<DeviceType>::tdual_ffloat_2d k_cut_ljsq;
-  typename ArrayTypes<DeviceType>::t_ffloat_2d d_cut_ljsq;
-  typename ArrayTypes<DeviceType>::tdual_ffloat_2d k_cut_coulsq;
-  typename ArrayTypes<DeviceType>::t_ffloat_2d d_cut_coulsq;
+  typename AT::tdual_ffloat_2d k_cutsq;
+  typename AT::t_ffloat_2d d_cutsq;
+  typename AT::tdual_ffloat_2d k_cut_ljsq;
+  typename AT::t_ffloat_2d d_cut_ljsq;
+  typename AT::tdual_ffloat_2d k_cut_coulsq;
+  typename AT::t_ffloat_2d d_cut_coulsq;

-  typename ArrayTypes<DeviceType>::t_ffloat_1d_randomread
+  typename AT::t_ffloat_1d_randomread
    d_rtable, d_drtable, d_ftable, d_dftable,
    d_ctable, d_dctable, d_etable, d_detable;

--- a/src/KOKKOS/pair_coul_wolf_kokkos.cpp
+++ b/src/KOKKOS/pair_coul_wolf_kokkos.cpp
@ -83,12 +83,12 @@ void PairCoulWolfKokkos<DeviceType>::compute(int eflag_in, int vflag_in)
  if (eflag_atom) {
    memory->destroy_kokkos(k_eatom,eatom);
    memory->create_kokkos(k_eatom,eatom,maxeatom,"pair:eatom");
-    d_eatom = k_eatom.d_view;
+    d_eatom = k_eatom.view<DeviceType>();
  }
  if (vflag_atom) {
    memory->destroy_kokkos(k_vatom,vatom);
    memory->create_kokkos(k_vatom,vatom,maxvatom,6,"pair:vatom");
-    d_vatom = k_vatom.d_view;
+    d_vatom = k_vatom.view<DeviceType>();
  }

  atomKK->sync(execution_space,datamask_read);
@ -184,8 +184,6 @@ void PairCoulWolfKokkos<DeviceType>::compute(int eflag_in, int vflag_in)
    virial[5] += ev.v[5];
  }

-  if (vflag_fdotr) pair_virial_fdotr_compute(this);
-
  if (eflag_atom) {
    k_eatom.template modify<DeviceType>();
    k_eatom.template sync<LMPHostType>();
@ -196,6 +194,8 @@ void PairCoulWolfKokkos<DeviceType>::compute(int eflag_in, int vflag_in)
    k_vatom.template sync<LMPHostType>();
  }

+  if (vflag_fdotr) pair_virial_fdotr_compute(this);
+
  copymode = 0;
 }

--- a/src/KOKKOS/pair_coul_wolf_kokkos.h
+++ b/src/KOKKOS/pair_coul_wolf_kokkos.h
@ -37,6 +37,7 @@ class PairCoulWolfKokkos : public PairCoulWolf {
  enum {EnabledNeighFlags=FULL|HALFTHREAD|HALF};
  enum {COUL_FLAG=1};
  typedef DeviceType device_type;
+  typedef ArrayTypes<DeviceType> AT;
  typedef EV_FLOAT value_type;
  PairCoulWolfKokkos(class LAMMPS *);
  ~PairCoulWolfKokkos();
@ -63,14 +64,14 @@ class PairCoulWolfKokkos : public PairCoulWolf {

 protected:

-  typename ArrayTypes<DeviceType>::t_x_array_randomread x;
-  typename ArrayTypes<DeviceType>::t_f_array f;
-  typename ArrayTypes<DeviceType>::t_float_1d_randomread q;
+  typename AT::t_x_array_randomread x;
+  typename AT::t_f_array f;
+  typename AT::t_float_1d_randomread q;

  DAT::tdual_efloat_1d k_eatom;
  DAT::tdual_virial_array k_vatom;
-  DAT::t_efloat_1d d_eatom;
-  DAT::t_virial_array d_vatom;
+  typename ArrayTypes<DeviceType>::t_efloat_1d d_eatom;
+  typename ArrayTypes<DeviceType>::t_virial_array d_vatom;


  int neighflag,newton_pair;
@ -81,9 +82,9 @@ class PairCoulWolfKokkos : public PairCoulWolf {
  double special_coul[4];
  double qqrd2e;

-  typename ArrayTypes<DeviceType>::t_neighbors_2d d_neighbors;
-  typename ArrayTypes<DeviceType>::t_int_1d_randomread d_ilist;
-  typename ArrayTypes<DeviceType>::t_int_1d_randomread d_numneigh;
+  typename AT::t_neighbors_2d d_neighbors;
+  typename AT::t_int_1d_randomread d_ilist;
+  typename AT::t_int_1d_randomread d_numneigh;
  //NeighListKokkos<DeviceType> k_list;

  friend void pair_virial_fdotr_compute<PairCoulWolfKokkos>(PairCoulWolfKokkos*);
--- a/src/KOKKOS/pair_eam_alloy_kokkos.cpp
+++ b/src/KOKKOS/pair_eam_alloy_kokkos.cpp
@ -82,12 +82,12 @@ void PairEAMAlloyKokkos<DeviceType>::compute(int eflag_in, int vflag_in)
  if (eflag_atom) {
    memory->destroy_kokkos(k_eatom,eatom);
    memory->create_kokkos(k_eatom,eatom,maxeatom,"pair:eatom");
-    d_eatom = k_eatom.d_view;
+    d_eatom = k_eatom.view<DeviceType>();
  }
  if (vflag_atom) {
    memory->destroy_kokkos(k_vatom,vatom);
    memory->create_kokkos(k_vatom,vatom,maxvatom,6,"pair:vatom");
-    d_vatom = k_vatom.d_view;
+    d_vatom = k_vatom.view<DeviceType>();
  }

  atomKK->sync(execution_space,datamask_read);
@ -246,8 +246,6 @@ void PairEAMAlloyKokkos<DeviceType>::compute(int eflag_in, int vflag_in)
    virial[5] += ev.v[5];
  }

-  if (vflag_fdotr) pair_virial_fdotr_compute(this);
-
  if (eflag_atom) {
    k_eatom.template modify<DeviceType>();
    k_eatom.template sync<LMPHostType>();
@ -258,6 +256,8 @@ void PairEAMAlloyKokkos<DeviceType>::compute(int eflag_in, int vflag_in)
    k_vatom.template sync<LMPHostType>();
  }

+  if (vflag_fdotr) pair_virial_fdotr_compute(this);
+
  copymode = 0;
 }

--- a/src/KOKKOS/pair_eam_alloy_kokkos.h
+++ b/src/KOKKOS/pair_eam_alloy_kokkos.h
@ -125,8 +125,8 @@ class PairEAMAlloyKokkos : public PairEAM {

  DAT::tdual_efloat_1d k_eatom;
  DAT::tdual_virial_array k_vatom;
-  DAT::t_efloat_1d d_eatom;
-  DAT::t_virial_array d_vatom;
+  typename ArrayTypes<DeviceType>::t_efloat_1d d_eatom;
+  typename ArrayTypes<DeviceType>::t_virial_array d_vatom;

  DAT::tdual_ffloat_1d k_rho;
  DAT::tdual_ffloat_1d k_fp;
@ -154,9 +154,9 @@ class PairEAMAlloyKokkos : public PairEAM {
  void interpolate(int, double, double *, t_host_ffloat_2d_n7, int);
  void read_file(char *);

-  typename ArrayTypes<DeviceType>::t_neighbors_2d d_neighbors;
-  typename ArrayTypes<DeviceType>::t_int_1d_randomread d_ilist;
-  typename ArrayTypes<DeviceType>::t_int_1d_randomread d_numneigh;
+  typename AT::t_neighbors_2d d_neighbors;
+  typename AT::t_int_1d_randomread d_ilist;
+  typename AT::t_int_1d_randomread d_numneigh;
  //NeighListKokkos<DeviceType> k_list;

  int iswap;
--- a/src/KOKKOS/pair_eam_fs_kokkos.cpp
+++ b/src/KOKKOS/pair_eam_fs_kokkos.cpp
@ -82,12 +82,12 @@ void PairEAMFSKokkos<DeviceType>::compute(int eflag_in, int vflag_in)
  if (eflag_atom) {
    memory->destroy_kokkos(k_eatom,eatom);
    memory->create_kokkos(k_eatom,eatom,maxeatom,"pair:eatom");
-    d_eatom = k_eatom.d_view;
+    d_eatom = k_eatom.view<DeviceType>();
  }
  if (vflag_atom) {
    memory->destroy_kokkos(k_vatom,vatom);
    memory->create_kokkos(k_vatom,vatom,maxvatom,6,"pair:vatom");
-    d_vatom = k_vatom.d_view;
+    d_vatom = k_vatom.view<DeviceType>();
  }

  atomKK->sync(execution_space,datamask_read);
--- a/src/KOKKOS/pair_eam_fs_kokkos.h
+++ b/src/KOKKOS/pair_eam_fs_kokkos.h
@ -125,8 +125,8 @@ class PairEAMFSKokkos : public PairEAM {

  DAT::tdual_efloat_1d k_eatom;
  DAT::tdual_virial_array k_vatom;
-  DAT::t_efloat_1d d_eatom;
-  DAT::t_virial_array d_vatom;
+  typename ArrayTypes<DeviceType>::t_efloat_1d d_eatom;
+  typename ArrayTypes<DeviceType>::t_virial_array d_vatom;

  DAT::tdual_ffloat_1d k_rho;
  DAT::tdual_ffloat_1d k_fp;
@ -154,9 +154,9 @@ class PairEAMFSKokkos : public PairEAM {
  void interpolate(int, double, double *, t_host_ffloat_2d_n7, int);
  void read_file(char *);

-  typename ArrayTypes<DeviceType>::t_neighbors_2d d_neighbors;
-  typename ArrayTypes<DeviceType>::t_int_1d_randomread d_ilist;
-  typename ArrayTypes<DeviceType>::t_int_1d_randomread d_numneigh;
+  typename AT::t_neighbors_2d d_neighbors;
+  typename AT::t_int_1d_randomread d_ilist;
+  typename AT::t_int_1d_randomread d_numneigh;
  //NeighListKokkos<DeviceType> k_list;

  int iswap;
--- a/src/KOKKOS/pair_eam_kokkos.cpp
+++ b/src/KOKKOS/pair_eam_kokkos.cpp
@ -77,12 +77,12 @@ void PairEAMKokkos<DeviceType>::compute(int eflag_in, int vflag_in)
  if (eflag_atom) {
    memory->destroy_kokkos(k_eatom,eatom);
    memory->create_kokkos(k_eatom,eatom,maxeatom,"pair:eatom");
-    d_eatom = k_eatom.d_view;
+    d_eatom = k_eatom.view<DeviceType>();
  }
  if (vflag_atom) {
    memory->destroy_kokkos(k_vatom,vatom);
    memory->create_kokkos(k_vatom,vatom,maxvatom,6,"pair:vatom");
-    d_vatom = k_vatom.d_view;
+    d_vatom = k_vatom.view<DeviceType>();
  }

  atomKK->sync(execution_space,datamask_read);
--- a/src/KOKKOS/pair_eam_kokkos.h
+++ b/src/KOKKOS/pair_eam_kokkos.h
@ -122,8 +122,8 @@ class PairEAMKokkos : public PairEAM {

  DAT::tdual_efloat_1d k_eatom;
  DAT::tdual_virial_array k_vatom;
-  DAT::t_efloat_1d d_eatom;
-  DAT::t_virial_array d_vatom;
+  typename ArrayTypes<DeviceType>::t_efloat_1d d_eatom;
+  typename ArrayTypes<DeviceType>::t_virial_array d_vatom;

  DAT::tdual_ffloat_1d k_rho;
  DAT::tdual_ffloat_1d k_fp;
@ -149,9 +149,9 @@ class PairEAMKokkos : public PairEAM {
  virtual void file2array();
  void array2spline();

-  typename ArrayTypes<DeviceType>::t_neighbors_2d d_neighbors;
-  typename ArrayTypes<DeviceType>::t_int_1d_randomread d_ilist;
-  typename ArrayTypes<DeviceType>::t_int_1d_randomread d_numneigh;
+  typename AT::t_neighbors_2d d_neighbors;
+  typename AT::t_int_1d_randomread d_ilist;
+  typename AT::t_int_1d_randomread d_numneigh;
  //NeighListKokkos<DeviceType> k_list;

  int iswap;
--- a/src/KOKKOS/pair_lj_charmm_coul_charmm_implicit_kokkos.cpp
+++ b/src/KOKKOS/pair_lj_charmm_coul_charmm_implicit_kokkos.cpp
@ -115,6 +115,19 @@ void PairLJCharmmCoulCharmmImplicitKokkos<DeviceType>::compute(int eflag_in, int
  if (eflag || vflag) ev_setup(eflag,vflag);
  else evflag = vflag_fdotr = 0;

+  // reallocate per-atom arrays if necessary
+
+  if (eflag_atom) {
+    memory->destroy_kokkos(k_eatom,eatom);
+    memory->create_kokkos(k_eatom,eatom,maxeatom,"pair:eatom");
+    d_eatom = k_eatom.view<DeviceType>();
+  }
+  if (vflag_atom) {
+    memory->destroy_kokkos(k_vatom,vatom);
+    memory->create_kokkos(k_vatom,vatom,maxvatom,6,"pair:vatom");
+    d_vatom = k_vatom.view<DeviceType>();
+  }
+
  atomKK->sync(execution_space,datamask_read);
  k_cutsq.template sync<DeviceType>();
  k_cut_ljsq.template sync<DeviceType>();
@ -167,6 +180,16 @@ void PairLJCharmmCoulCharmmImplicitKokkos<DeviceType>::compute(int eflag_in, int
    virial[5] += ev.v[5];
  }

+  if (eflag_atom) {
+    k_eatom.template modify<DeviceType>();
+    k_eatom.template sync<LMPHostType>();
+  }
+
+  if (vflag_atom) {
+    k_vatom.template modify<DeviceType>();
+    k_vatom.template sync<LMPHostType>();
+  }
+
  if (vflag_fdotr) pair_virial_fdotr_compute(this);

  copymode = 0;
--- a/src/KOKKOS/pair_lj_charmm_coul_charmm_implicit_kokkos.h
+++ b/src/KOKKOS/pair_lj_charmm_coul_charmm_implicit_kokkos.h
@ -34,6 +34,7 @@ class PairLJCharmmCoulCharmmImplicitKokkos : public PairLJCharmmCoulCharmmImplic
  enum {EnabledNeighFlags=FULL|HALFTHREAD|HALF};
  enum {COUL_FLAG=1};
  typedef DeviceType device_type;
+  typedef ArrayTypes<DeviceType> AT;
  PairLJCharmmCoulCharmmImplicitKokkos(class LAMMPS *);
  ~PairLJCharmmCoulCharmmImplicitKokkos();

@ -77,27 +78,27 @@ class PairLJCharmmCoulCharmmImplicitKokkos : public PairLJCharmmCoulCharmmImplic
  F_FLOAT m_cutsq[MAX_TYPES_STACKPARAMS+1][MAX_TYPES_STACKPARAMS+1];
  F_FLOAT m_cut_ljsq[MAX_TYPES_STACKPARAMS+1][MAX_TYPES_STACKPARAMS+1];
  F_FLOAT m_cut_coulsq[MAX_TYPES_STACKPARAMS+1][MAX_TYPES_STACKPARAMS+1];
-  typename ArrayTypes<DeviceType>::t_x_array_randomread x;
-  typename ArrayTypes<DeviceType>::t_x_array c_x;
-  typename ArrayTypes<DeviceType>::t_f_array f;
-  typename ArrayTypes<DeviceType>::t_int_1d_randomread type;
-  typename ArrayTypes<DeviceType>::t_float_1d_randomread q;
+  typename AT::t_x_array_randomread x;
+  typename AT::t_x_array c_x;
+  typename AT::t_f_array f;
+  typename AT::t_int_1d_randomread type;
+  typename AT::t_float_1d_randomread q;

  DAT::tdual_efloat_1d k_eatom;
  DAT::tdual_virial_array k_vatom;
-  typename ArrayTypes<DeviceType>::t_efloat_1d d_eatom;
-  typename ArrayTypes<DeviceType>::t_virial_array d_vatom;
+  typename AT::t_efloat_1d d_eatom;
+  typename AT::t_virial_array d_vatom;

  int newton_pair;

-  typename ArrayTypes<DeviceType>::tdual_ffloat_2d k_cutsq;
-  typename ArrayTypes<DeviceType>::t_ffloat_2d d_cutsq;
-  typename ArrayTypes<DeviceType>::tdual_ffloat_2d k_cut_ljsq;
-  typename ArrayTypes<DeviceType>::t_ffloat_2d d_cut_ljsq;
-  typename ArrayTypes<DeviceType>::tdual_ffloat_2d k_cut_coulsq;
-  typename ArrayTypes<DeviceType>::t_ffloat_2d d_cut_coulsq;
+  typename AT::tdual_ffloat_2d k_cutsq;
+  typename AT::t_ffloat_2d d_cutsq;
+  typename AT::tdual_ffloat_2d k_cut_ljsq;
+  typename AT::t_ffloat_2d d_cut_ljsq;
+  typename AT::tdual_ffloat_2d k_cut_coulsq;
+  typename AT::t_ffloat_2d d_cut_coulsq;

-  typename ArrayTypes<DeviceType>::t_ffloat_1d_randomread
+  typename AT::t_ffloat_1d_randomread
    d_rtable, d_drtable, d_ftable, d_dftable,
    d_ctable, d_dctable, d_etable, d_detable;

--- a/src/KOKKOS/pair_lj_charmm_coul_charmm_kokkos.cpp
+++ b/src/KOKKOS/pair_lj_charmm_coul_charmm_kokkos.cpp
@ -115,6 +115,19 @@ void PairLJCharmmCoulCharmmKokkos<DeviceType>::compute(int eflag_in, int vflag_i
  if (eflag || vflag) ev_setup(eflag,vflag);
  else evflag = vflag_fdotr = 0;

+  // reallocate per-atom arrays if necessary
+
+  if (eflag_atom) {
+    memory->destroy_kokkos(k_eatom,eatom);
+    memory->create_kokkos(k_eatom,eatom,maxeatom,"pair:eatom");
+    d_eatom = k_eatom.view<DeviceType>();
+  }
+  if (vflag_atom) {
+    memory->destroy_kokkos(k_vatom,vatom);
+    memory->create_kokkos(k_vatom,vatom,maxvatom,6,"pair:vatom");
+    d_vatom = k_vatom.view<DeviceType>();
+  }
+
  atomKK->sync(execution_space,datamask_read);
  k_cutsq.template sync<DeviceType>();
  k_cut_ljsq.template sync<DeviceType>();
@ -167,6 +180,16 @@ void PairLJCharmmCoulCharmmKokkos<DeviceType>::compute(int eflag_in, int vflag_i
    virial[5] += ev.v[5];
  }

+  if (eflag_atom) {
+    k_eatom.template modify<DeviceType>();
+    k_eatom.template sync<LMPHostType>();
+  }
+
+  if (vflag_atom) {
+    k_vatom.template modify<DeviceType>();
+    k_vatom.template sync<LMPHostType>();
+  }
+
  if (vflag_fdotr) pair_virial_fdotr_compute(this);

  copymode = 0;
--- a/src/KOKKOS/pair_lj_charmm_coul_charmm_kokkos.h
+++ b/src/KOKKOS/pair_lj_charmm_coul_charmm_kokkos.h
@ -34,6 +34,7 @@ class PairLJCharmmCoulCharmmKokkos : public PairLJCharmmCoulCharmm {
  enum {EnabledNeighFlags=FULL|HALFTHREAD|HALF};
  enum {COUL_FLAG=1};
  typedef DeviceType device_type;
+  typedef ArrayTypes<DeviceType> AT;
  PairLJCharmmCoulCharmmKokkos(class LAMMPS *);
  ~PairLJCharmmCoulCharmmKokkos();

@ -77,27 +78,27 @@ class PairLJCharmmCoulCharmmKokkos : public PairLJCharmmCoulCharmm {
  F_FLOAT m_cutsq[MAX_TYPES_STACKPARAMS+1][MAX_TYPES_STACKPARAMS+1];
  F_FLOAT m_cut_ljsq[MAX_TYPES_STACKPARAMS+1][MAX_TYPES_STACKPARAMS+1];
  F_FLOAT m_cut_coulsq[MAX_TYPES_STACKPARAMS+1][MAX_TYPES_STACKPARAMS+1];
-  typename ArrayTypes<DeviceType>::t_x_array_randomread x;
-  typename ArrayTypes<DeviceType>::t_x_array c_x;
-  typename ArrayTypes<DeviceType>::t_f_array f;
-  typename ArrayTypes<DeviceType>::t_int_1d_randomread type;
-  typename ArrayTypes<DeviceType>::t_float_1d_randomread q;
+  typename AT::t_x_array_randomread x;
+  typename AT::t_x_array c_x;
+  typename AT::t_f_array f;
+  typename AT::t_int_1d_randomread type;
+  typename AT::t_float_1d_randomread q;

  DAT::tdual_efloat_1d k_eatom;
  DAT::tdual_virial_array k_vatom;
-  typename ArrayTypes<DeviceType>::t_efloat_1d d_eatom;
-  typename ArrayTypes<DeviceType>::t_virial_array d_vatom;
+  typename AT::t_efloat_1d d_eatom;
+  typename AT::t_virial_array d_vatom;

  int newton_pair;

-  typename ArrayTypes<DeviceType>::tdual_ffloat_2d k_cutsq;
-  typename ArrayTypes<DeviceType>::t_ffloat_2d d_cutsq;
-  typename ArrayTypes<DeviceType>::tdual_ffloat_2d k_cut_ljsq;
-  typename ArrayTypes<DeviceType>::t_ffloat_2d d_cut_ljsq;
-  typename ArrayTypes<DeviceType>::tdual_ffloat_2d k_cut_coulsq;
-  typename ArrayTypes<DeviceType>::t_ffloat_2d d_cut_coulsq;
+  typename AT::tdual_ffloat_2d k_cutsq;
+  typename AT::t_ffloat_2d d_cutsq;
+  typename AT::tdual_ffloat_2d k_cut_ljsq;
+  typename AT::t_ffloat_2d d_cut_ljsq;
+  typename AT::tdual_ffloat_2d k_cut_coulsq;
+  typename AT::t_ffloat_2d d_cut_coulsq;

-  typename ArrayTypes<DeviceType>::t_ffloat_1d_randomread
+  typename AT::t_ffloat_1d_randomread
    d_rtable, d_drtable, d_ftable, d_dftable,
    d_ctable, d_dctable, d_etable, d_detable;

--- a/src/KOKKOS/pair_lj_charmm_coul_long_kokkos.cpp
+++ b/src/KOKKOS/pair_lj_charmm_coul_long_kokkos.cpp
@ -180,8 +180,6 @@ void PairLJCharmmCoulLongKokkos<DeviceType>::compute(int eflag_in, int vflag_in)
    virial[5] += ev.v[5];
  }

-  if (vflag_fdotr) pair_virial_fdotr_compute(this);
-
  if (eflag_atom) {
    k_eatom.template modify<DeviceType>();
    k_eatom.template sync<LMPHostType>();
@ -192,6 +190,8 @@ void PairLJCharmmCoulLongKokkos<DeviceType>::compute(int eflag_in, int vflag_in)
    k_vatom.template sync<LMPHostType>();
  }

+  if (vflag_fdotr) pair_virial_fdotr_compute(this);
+
  copymode = 0;
 }

--- a/src/KOKKOS/pair_lj_charmm_coul_long_kokkos.h
+++ b/src/KOKKOS/pair_lj_charmm_coul_long_kokkos.h
@ -34,6 +34,7 @@ class PairLJCharmmCoulLongKokkos : public PairLJCharmmCoulLong {
  enum {EnabledNeighFlags=FULL|HALFTHREAD|HALF};
  enum {COUL_FLAG=1};
  typedef DeviceType device_type;
+  typedef ArrayTypes<DeviceType> AT;
  PairLJCharmmCoulLongKokkos(class LAMMPS *);
  ~PairLJCharmmCoulLongKokkos();

@ -75,27 +76,27 @@ class PairLJCharmmCoulLongKokkos : public PairLJCharmmCoulLong {
  F_FLOAT m_cutsq[MAX_TYPES_STACKPARAMS+1][MAX_TYPES_STACKPARAMS+1];
  F_FLOAT m_cut_ljsq[MAX_TYPES_STACKPARAMS+1][MAX_TYPES_STACKPARAMS+1];
  F_FLOAT m_cut_coulsq[MAX_TYPES_STACKPARAMS+1][MAX_TYPES_STACKPARAMS+1];
-  typename ArrayTypes<DeviceType>::t_x_array_randomread x;
-  typename ArrayTypes<DeviceType>::t_x_array c_x;
-  typename ArrayTypes<DeviceType>::t_f_array f;
-  typename ArrayTypes<DeviceType>::t_int_1d_randomread type;
-  typename ArrayTypes<DeviceType>::t_float_1d_randomread q;
+  typename AT::t_x_array_randomread x;
+  typename AT::t_x_array c_x;
+  typename AT::t_f_array f;
+  typename AT::t_int_1d_randomread type;
+  typename AT::t_float_1d_randomread q;

  DAT::tdual_efloat_1d k_eatom;
  DAT::tdual_virial_array k_vatom;
-  typename ArrayTypes<DeviceType>::t_efloat_1d d_eatom;
-  typename ArrayTypes<DeviceType>::t_virial_array d_vatom;
+  typename AT::t_efloat_1d d_eatom;
+  typename AT::t_virial_array d_vatom;

  int newton_pair;

-  typename ArrayTypes<DeviceType>::tdual_ffloat_2d k_cutsq;
-  typename ArrayTypes<DeviceType>::t_ffloat_2d d_cutsq;
-  typename ArrayTypes<DeviceType>::tdual_ffloat_2d k_cut_ljsq;
-  typename ArrayTypes<DeviceType>::t_ffloat_2d d_cut_ljsq;
-  typename ArrayTypes<DeviceType>::tdual_ffloat_2d k_cut_coulsq;
-  typename ArrayTypes<DeviceType>::t_ffloat_2d d_cut_coulsq;
+  typename AT::tdual_ffloat_2d k_cutsq;
+  typename AT::t_ffloat_2d d_cutsq;
+  typename AT::tdual_ffloat_2d k_cut_ljsq;
+  typename AT::t_ffloat_2d d_cut_ljsq;
+  typename AT::tdual_ffloat_2d k_cut_coulsq;
+  typename AT::t_ffloat_2d d_cut_coulsq;

-  typename ArrayTypes<DeviceType>::t_ffloat_1d_randomread
+  typename AT::t_ffloat_1d_randomread
    d_rtable, d_drtable, d_ftable, d_dftable,
    d_ctable, d_dctable, d_etable, d_detable;

--- a/src/KOKKOS/pair_lj_class2_coul_cut_kokkos.cpp
+++ b/src/KOKKOS/pair_lj_class2_coul_cut_kokkos.cpp
@ -92,6 +92,19 @@ void PairLJClass2CoulCutKokkos<DeviceType>::compute(int eflag_in, int vflag_in)
  if (eflag || vflag) ev_setup(eflag,vflag);
  else evflag = vflag_fdotr = 0;

+  // reallocate per-atom arrays if necessary
+
+  if (eflag_atom) {
+    memory->destroy_kokkos(k_eatom,eatom);
+    memory->create_kokkos(k_eatom,eatom,maxeatom,"pair:eatom");
+    d_eatom = k_eatom.view<DeviceType>();
+  }
+  if (vflag_atom) {
+    memory->destroy_kokkos(k_vatom,vatom);
+    memory->create_kokkos(k_vatom,vatom,maxvatom,6,"pair:vatom");
+    d_vatom = k_vatom.view<DeviceType>();
+  }
+
  atomKK->sync(execution_space,datamask_read);
  k_cutsq.template sync<DeviceType>();
  k_cut_ljsq.template sync<DeviceType>();
@ -138,7 +151,17 @@ void PairLJClass2CoulCutKokkos<DeviceType>::compute(int eflag_in, int vflag_in)
    virial[5] += ev.v[5];
  }

-  if (vflag_fdotr) virial_fdotr_compute();
+  if (eflag_atom) {
+    k_eatom.template modify<DeviceType>();
+    k_eatom.template sync<LMPHostType>();
+  }
+
+  if (vflag_atom) {
+    k_vatom.template modify<DeviceType>();
+    k_vatom.template sync<LMPHostType>();
+  }
+
+  if (vflag_fdotr) pair_virial_fdotr_compute(this);

  copymode = 0;
 }
--- a/src/KOKKOS/pair_lj_class2_coul_cut_kokkos.h
+++ b/src/KOKKOS/pair_lj_class2_coul_cut_kokkos.h
@ -34,6 +34,7 @@ class PairLJClass2CoulCutKokkos : public PairLJClass2CoulCut {
  enum {EnabledNeighFlags=FULL|HALFTHREAD|HALF};
  enum {COUL_FLAG=1};
  typedef DeviceType device_type;
+  typedef ArrayTypes<DeviceType> AT;
  PairLJClass2CoulCutKokkos(class LAMMPS *);
  ~PairLJClass2CoulCutKokkos();

@ -76,22 +77,25 @@ class PairLJClass2CoulCutKokkos : public PairLJClass2CoulCut {
  F_FLOAT m_cutsq[MAX_TYPES_STACKPARAMS+1][MAX_TYPES_STACKPARAMS+1];
  F_FLOAT m_cut_ljsq[MAX_TYPES_STACKPARAMS+1][MAX_TYPES_STACKPARAMS+1];
  F_FLOAT m_cut_coulsq[MAX_TYPES_STACKPARAMS+1][MAX_TYPES_STACKPARAMS+1];
-  typename ArrayTypes<DeviceType>::t_x_array_randomread x;
-  typename ArrayTypes<DeviceType>::t_x_array c_x;
-  typename ArrayTypes<DeviceType>::t_f_array f;
-  typename ArrayTypes<DeviceType>::t_int_1d_randomread type;
-  typename ArrayTypes<DeviceType>::t_float_1d_randomread q;
-  typename ArrayTypes<DeviceType>::t_efloat_1d d_eatom;
-  typename ArrayTypes<DeviceType>::t_virial_array d_vatom;
+  typename AT::t_x_array_randomread x;
+  typename AT::t_x_array c_x;
+  typename AT::t_f_array f;
+  typename AT::t_int_1d_randomread type;
+  typename AT::t_float_1d_randomread q;
+
+  DAT::tdual_efloat_1d k_eatom;
+  DAT::tdual_virial_array k_vatom;
+  typename AT::t_efloat_1d d_eatom;
+  typename AT::t_virial_array d_vatom;

  int newton_pair;

-  typename ArrayTypes<DeviceType>::tdual_ffloat_2d k_cutsq;
-  typename ArrayTypes<DeviceType>::t_ffloat_2d d_cutsq;
-  typename ArrayTypes<DeviceType>::tdual_ffloat_2d k_cut_ljsq;
-  typename ArrayTypes<DeviceType>::t_ffloat_2d d_cut_ljsq;
-  typename ArrayTypes<DeviceType>::tdual_ffloat_2d k_cut_coulsq;
-  typename ArrayTypes<DeviceType>::t_ffloat_2d d_cut_coulsq;
+  typename AT::tdual_ffloat_2d k_cutsq;
+  typename AT::t_ffloat_2d d_cutsq;
+  typename AT::tdual_ffloat_2d k_cut_ljsq;
+  typename AT::t_ffloat_2d d_cut_ljsq;
+  typename AT::tdual_ffloat_2d k_cut_coulsq;
+  typename AT::t_ffloat_2d d_cut_coulsq;

  int neighflag;
  int nlocal,nall,eflag,vflag;
@ -112,6 +116,7 @@ class PairLJClass2CoulCutKokkos : public PairLJClass2CoulCut {
  friend EV_FLOAT pair_compute_neighlist<PairLJClass2CoulCutKokkos,HALFTHREAD,void>(PairLJClass2CoulCutKokkos*,NeighListKokkos<DeviceType>*);
  friend EV_FLOAT pair_compute<PairLJClass2CoulCutKokkos,void>(PairLJClass2CoulCutKokkos*,
                                                            NeighListKokkos<DeviceType>*);
+  friend void pair_virial_fdotr_compute<PairLJClass2CoulCutKokkos>(PairLJClass2CoulCutKokkos*);

 };

--- a/src/KOKKOS/pair_lj_class2_coul_long_kokkos.cpp
+++ b/src/KOKKOS/pair_lj_class2_coul_long_kokkos.cpp
@ -100,6 +100,19 @@ void PairLJClass2CoulLongKokkos<DeviceType>::compute(int eflag_in, int vflag_in)
  if (eflag || vflag) ev_setup(eflag,vflag);
  else evflag = vflag_fdotr = 0;

+  // reallocate per-atom arrays if necessary
+
+  if (eflag_atom) {
+    memory->destroy_kokkos(k_eatom,eatom);
+    memory->create_kokkos(k_eatom,eatom,maxeatom,"pair:eatom");
+    d_eatom = k_eatom.view<DeviceType>();
+  }
+  if (vflag_atom) {
+    memory->destroy_kokkos(k_vatom,vatom);
+    memory->create_kokkos(k_vatom,vatom,maxvatom,6,"pair:vatom");
+    d_vatom = k_vatom.view<DeviceType>();
+  }
+
  atomKK->sync(execution_space,datamask_read);
  k_cutsq.template sync<DeviceType>();
  k_cut_ljsq.template sync<DeviceType>();
@ -152,6 +165,16 @@ void PairLJClass2CoulLongKokkos<DeviceType>::compute(int eflag_in, int vflag_in)
    virial[5] += ev.v[5];
  }

+  if (eflag_atom) {
+    k_eatom.template modify<DeviceType>();
+    k_eatom.template sync<LMPHostType>();
+  }
+
+  if (vflag_atom) {
+    k_vatom.template modify<DeviceType>();
+    k_vatom.template sync<LMPHostType>();
+  }
+
  if (vflag_fdotr) pair_virial_fdotr_compute(this);

  copymode = 0;
--- a/src/KOKKOS/pair_lj_class2_coul_long_kokkos.h
+++ b/src/KOKKOS/pair_lj_class2_coul_long_kokkos.h
@ -34,6 +34,7 @@ class PairLJClass2CoulLongKokkos : public PairLJClass2CoulLong {
  enum {EnabledNeighFlags=FULL|HALFTHREAD|HALF};
  enum {COUL_FLAG=1};
  typedef DeviceType device_type;
+  typedef ArrayTypes<DeviceType> AT;
  PairLJClass2CoulLongKokkos(class LAMMPS *);
  ~PairLJClass2CoulLongKokkos();

@ -76,24 +77,27 @@ class PairLJClass2CoulLongKokkos : public PairLJClass2CoulLong {
  F_FLOAT m_cutsq[MAX_TYPES_STACKPARAMS+1][MAX_TYPES_STACKPARAMS+1];
  F_FLOAT m_cut_ljsq[MAX_TYPES_STACKPARAMS+1][MAX_TYPES_STACKPARAMS+1];
  F_FLOAT m_cut_coulsq[MAX_TYPES_STACKPARAMS+1][MAX_TYPES_STACKPARAMS+1];
-  typename ArrayTypes<DeviceType>::t_x_array_randomread x;
-  typename ArrayTypes<DeviceType>::t_x_array c_x;
-  typename ArrayTypes<DeviceType>::t_f_array f;
-  typename ArrayTypes<DeviceType>::t_int_1d_randomread type;
-  typename ArrayTypes<DeviceType>::t_float_1d_randomread q;
-  typename ArrayTypes<DeviceType>::t_efloat_1d d_eatom;
-  typename ArrayTypes<DeviceType>::t_virial_array d_vatom;
+  typename AT::t_x_array_randomread x;
+  typename AT::t_x_array c_x;
+  typename AT::t_f_array f;
+  typename AT::t_int_1d_randomread type;
+  typename AT::t_float_1d_randomread q;
+
+  DAT::tdual_efloat_1d k_eatom;
+  DAT::tdual_virial_array k_vatom;
+  typename AT::t_efloat_1d d_eatom;
+  typename AT::t_virial_array d_vatom;

  int newton_pair;

-  typename ArrayTypes<DeviceType>::tdual_ffloat_2d k_cutsq;
-  typename ArrayTypes<DeviceType>::t_ffloat_2d d_cutsq;
-  typename ArrayTypes<DeviceType>::tdual_ffloat_2d k_cut_ljsq;
-  typename ArrayTypes<DeviceType>::t_ffloat_2d d_cut_ljsq;
-  typename ArrayTypes<DeviceType>::tdual_ffloat_2d k_cut_coulsq;
-  typename ArrayTypes<DeviceType>::t_ffloat_2d d_cut_coulsq;
+  typename AT::tdual_ffloat_2d k_cutsq;
+  typename AT::t_ffloat_2d d_cutsq;
+  typename AT::tdual_ffloat_2d k_cut_ljsq;
+  typename AT::t_ffloat_2d d_cut_ljsq;
+  typename AT::tdual_ffloat_2d k_cut_coulsq;
+  typename AT::t_ffloat_2d d_cut_coulsq;

-  typename ArrayTypes<DeviceType>::t_ffloat_1d_randomread
+  typename AT::t_ffloat_1d_randomread
    d_rtable, d_drtable, d_ftable, d_dftable,
    d_ctable, d_dctable, d_etable, d_detable;

--- a/Show More
+++ b/Show More
Author	SHA1	Message	Date
Steve Plimpton	cb982f2f28	sync 13Feb17 patch back to GH	2017-02-13 09:05:03 -07:00
sjplimp	4843296d4e	Merge pull request #372 from akohlmey/fft-cleanup simplify FFT3d code by removing support for outdated FFT libraries	2017-02-13 08:53:13 -07:00
Steve Plimpton	2bdda8f6c0	patch 12Feb17 - change int to tagint for compute group/group	2017-02-13 08:40:54 -07:00
Steve Plimpton	0068ef5616	added molecule option to compute group/group command	2017-02-10 09:25:32 -07:00
sjplimp	02b0e6cc55	Merge pull request #375 from akohlmey/small-updates-and-fixes Small updates and fixes	2017-02-10 09:23:51 -07:00
sjplimp	fbb24c2406	Merge pull request #374 from agiliopadua/master Updated polarizer.py in USER-DRUDE to use coul/long/cs	2017-02-10 09:22:22 -07:00
Axel Kohlmeyer	a5f830c40c	fix typo (cherry picked from commit `6410797697`)	2017-02-08 14:33:45 -05:00
Axel Kohlmeyer	8c074a363a	Merge branch 'master' into small-updates-and-fixes	2017-02-08 14:32:44 -05:00
Agilio Padua	27aca14094	Updated polarizer.py to use coul/long/cs	2017-02-04 15:02:08 +01:00
Axel Kohlmeyer	191453e1c7	Merge branch 'master' into fft-cleanup	2017-02-03 16:53:10 -05:00
sjplimp	207adc3968	Merge pull request #373 from stanmoore1/kk_more_bugfixes Fixing Kokkos per-atom energy/virial issues	2017-02-03 14:45:31 -07:00
sjplimp	84c517159d	Merge pull request #368 from Pakketeretet2/kokkos_morse Kokkos morse	2017-02-03 14:45:07 -07:00
sjplimp	6ca377436f	Merge pull request #366 from rbberger/kokkos_lammps_bigbig_fix Fix data type of molecule array in npair_kokkos.h	2017-02-03 14:43:33 -07:00
sjplimp	dc34a32602	Merge pull request #362 from ibaned/warnings2 fix Kokkos+kspace warnings	2017-02-03 14:43:12 -07:00
Stan Moore	067119f6c6	Adding missing friend statement to pair_lj_class2_coul_cut_kokkos	2017-02-02 15:21:30 -07:00
Stan Moore	1834a5e46c	Fixing more Kokkos per-atom and fdotr issues	2017-02-02 15:21:21 -07:00
Stan Moore	6a4918b39a	Fixing typo in pair_buck_coul_cut_kokkos	2017-02-02 15:21:05 -07:00
Stan Moore	5da0d39392	Fixing fdotr in pair_buck_coul_cut_kokkos	2017-02-02 13:35:51 -07:00
Stan Moore	6f92429602	Fixing per-atom ev issue	2017-02-02 13:34:27 -07:00
Stan Moore	38e0e4bb69	Add missing typedef in Kokkos pair styles	2017-02-02 13:24:05 -07:00
Stan Moore	daf9f95381	Fixing Kokkos per-atom e/v issue	2017-02-02 13:09:52 -07:00
Axel Kohlmeyer	6595fde0a1	explain in more detail the handling of error checking for numerical inputs	2017-02-02 11:58:12 -05:00
Richard Berger	6bcec9c61d	Merge pull request #2 from stanmoore1/kk_tag_bugfixes Fixing tagint and imageint issues in Kokkos package	2017-02-02 08:57:21 -05:00
Axel Kohlmeyer	9d1991bf84	remove support for obsolete legacy FFT libraries and point -DFFT_FFTW to FFTW3	2017-02-02 08:10:23 -05:00
Stefan Paquay	0a87b7443a	Updated contributing authors and docs	2017-02-02 13:42:47 +01:00
Stan Moore	7ee45ec5f3	Fixing tagint and imageint issues in Kokkos package	2017-02-01 11:52:27 -07:00
Stefan Paquay	d4c9e2500b	Ported Morse to KOKKOS	2017-02-01 17:45:21 +01:00
Stefan Paquay	6232073d3b	Removed traces of pair morse/kk	2017-02-01 17:39:37 +01:00
Stefan Paquay	ed59193d13	Removed traces of pair morse/kk	2017-02-01 17:39:06 +01:00
Richard Berger	67bed8e853	Merge pull request #1 from akohlmey/tagint-issue Fix additional tagint issue in fix qeq/reax/kk	2017-01-31 18:34:35 -05:00
Axel Kohlmeyer	bcb1d94b9a	silence compiler warning about dead code	2017-01-31 18:28:04 -05:00
Axel Kohlmeyer	fbe30b5683	correct issue with compiling for -DLAMMPS_BIGBIG in fix qeq/reax/kk	2017-01-31 18:13:44 -05:00
Axel Kohlmeyer	9ef55fedf7	Merge branch 'kokkos_lammps_bigbig_fix' of https://github.com/rbberger/lammps into tagint-issue	2017-01-31 17:23:51 -05:00
sjplimp	997142a4c1	Merge pull request #364 from stanmoore1/kk_triclinic_neighlist Add triclinic neighbor list support to Kokkos	2017-01-30 07:27:02 -07:00
sjplimp	033b07fdb7	Merge pull request #363 from ibaned/obey-datamask Fix GPU sync bugs	2017-01-30 07:26:49 -07:00
Richard Berger	51a0b6b445	Fix data type of molecule array in npair_kokkos.h This showed up when trying to compile with -DLAMMPS_BIGBIG. Fixes issue #365	2017-01-28 07:49:08 -05:00
Stan Moore	59f4a77dd5	Whitespace change to npair_kokkos	2017-01-27 15:17:39 -07:00
Stan Moore	579cc6d7aa	More tweaks to npair_kokkos for triclinic	2017-01-27 15:13:37 -07:00
Stan Moore	5afd3e995b	Adding support to npair_kokkos for triclinic-newton-on neighborlists	2017-01-27 14:18:01 -07:00
Dan Ibanez	2a6f5e651c	more preference of datamask over custom sync see commit `09fc8b0` for details on why	2017-01-27 09:35:55 -07:00
Dan Ibanez	09fc8b0bd7	kspace & dihedral can't do their own sync/modify because the verlet_kokkos system has a "clever" optimization which will alter the datamasks before calling sync/modify, so the datamask framework must be strictly obeyed for GPU correctness. (the optimization is to concurrently compute forces on the host and GPU, and add them up at the end of an iteration. calling your own sync will overwrite the partial GPU forces with the partial host forces).	2017-01-27 08:39:55 -07:00
Dan Ibanez	e5d0bde783	pppm_kokkos: remove useless statement	2017-01-27 08:35:37 -07:00
Dan Ibanez	9daf7fb650	pppm_kokkos: don't shadow member variables	2017-01-27 08:35:37 -07:00
Dan Ibanez	b5d622c6a3	pppm_kokkos: remove unused variables	2017-01-27 08:35:37 -07:00
Dan Ibanez	2023fa28e0	consistent #ifdefs for fft3d variable (2) this variable is only used when FFTW3 is enabled, so its declaration and initialization should be protected under the same conditions to avoid compiler warnings	2017-01-27 08:35:37 -07:00
Dan Ibanez	5b29515849	fft3d: use C++ loop declarations the variable (offset) is only used in a subset of numerous scenarios with #ifdef, it seems better just to have each loop declare it as needed. (avoids compiler warnings)	2017-01-27 08:35:37 -07:00
Dan Ibanez	5b18421dd2	fft3d : remove unused variables	2017-01-27 08:35:37 -07:00
Dan Ibanez	cf95ea0709	fft3d: only declare variables when used avoids compiler warnings	2017-01-27 08:35:36 -07:00
Dan Ibanez	6a74a81da0	consistent #ifdefs for fft3d variable this variable is only used when FFTW3 is enabled, so its declaration and initialization should be protected under the same conditions to avoid compiler warnings	2017-01-27 08:35:36 -07:00
Dan Ibanez	f0a4ed615d	add missing KOKKOS_INLINE_FUNCTION for params	2017-01-27 08:35:36 -07:00
Dan Ibanez	cfe818a175	remove unused variables from fix_cmap	2017-01-27 08:35:36 -07:00
Steve Plimpton	f8506fee23	sync GHub with SVN	2017-01-26 14:06:43 -07:00
sjplimp	18e5584311	Merge pull request #354 from stanmoore1/kokkos_bugfixes Kokkos bugfixes	2017-01-26 13:51:47 -07:00
sjplimp	851f80464f	Merge pull request #361 from akohlmey/user-omp-fix-per-atom-data fix USER-OMP bug on per-atom data with hybrid styles	2017-01-26 13:50:13 -07:00
sjplimp	5971d4c994	Merge pull request #358 from ibaned/warnings warning fixes (Kokkos+CUDA)	2017-01-26 13:49:56 -07:00
sjplimp	868d95f0a5	Merge pull request #352 from akohlmey/fix-skip-with-ghost-issue Fix skip with ghost issue	2017-01-26 13:47:12 -07:00
sjplimp	a5ff35435a	Merge pull request #351 from timattox/USER-DPD_pair_exp6_bugfix USER-DPD: Possible uninitialized variable in pair_exp6_rx.cpp bugfix.	2017-01-26 13:45:45 -07:00
Axel Kohlmeyer	8b7bd9d88e	fix bug where per atom data for USER-OMP was reducing the wrong arrays with hybrid styles	2017-01-26 14:59:10 -05:00
Dan Ibanez	672bbbe494	add more missing KOKKOS_INLINE_FUNCTION attributes	2017-01-25 16:03:11 -07:00
Dan Ibanez	03c9c46533	add missing KOKKOS_INLINE_FUNCTION attributes	2017-01-25 15:49:05 -07:00
Dan Ibanez	e992bfe510	remove unused variable	2017-01-25 15:40:52 -07:00
Dan Ibanez	053ee54a27	remove unused variable	2017-01-25 15:38:41 -07:00
Dan Ibanez	1074c6734b	add missing return keywords	2017-01-25 15:37:27 -07:00
Dan Ibanez	60b48c9d66	add missing KOKKOS_INLINE_FUNCTION attributes this structure gets put inside a DualView, so these members need to be able to execute on the GPU	2017-01-25 15:36:24 -07:00
Dan Ibanez	3d40b51708	remove unused variable	2017-01-25 15:24:52 -07:00
Dan Ibanez	effbe18c46	fix domain boundary indexing the compiler pointed out that boundary[2][2] doesn't exist. If I understand this correctly, those checks should be against boundary[*][0].	2017-01-25 15:24:01 -07:00
Dan Ibanez	6328beb7d7	fix double-return warning this #ifdef adds a return statement for little endian machines, but leaves the old one, which the compiler comlains is unreachable. this commit combines the conditionals so we can use #else	2017-01-25 15:22:42 -07:00
Stan Moore	26c8d3d98f	Fixing GPU memory issue in fix_langevin_kokkos	2017-01-25 12:53:55 -07:00
Stan Moore	73177d650d	Fixing GPU memory issue in domain_kokkos	2017-01-25 11:18:03 -07:00
Axel Kohlmeyer	b5cb74bd33	skip list build is compatible with NP_GHOST	2017-01-23 19:21:48 -05:00
Axel Kohlmeyer	31976d1dee	skip list definition was missing NP_HALFFULL flag	2017-01-23 19:20:05 -05:00
Christopher Stone	c8260af37c	Possible uninitialized variable in USER-DPD/pair_exp6_rx.cpp bugfix. Added explicit initialization (to zero) for several variables inside the inner j-loop to avoid using them uninitialized or from prior iterations within rmOldij_12 == 0.	2017-01-23 13:34:51 -05:00
Steve Plimpton	caea8973a3	add neighbor list kind output to screen	2017-01-20 13:24:09 -07:00
sjplimp	aa0ad9b483	Merge pull request #349 from akohlmey/collected-small-fixes collected fixes and improvements	2017-01-20 13:19:43 -07:00
sjplimp	5d0e4e1ba9	Merge pull request #346 from stanmoore1/kokkos_fixes Kokkos fixes	2017-01-20 13:15:16 -07:00
sjplimp	f8d3c4c740	Merge pull request #345 from timattox/USER-DPD_another_zero_compute USER-DPD another zero compute optimization	2017-01-20 13:14:59 -07:00
Axel Kohlmeyer	e6996121d1	remove dead code	2017-01-20 14:30:46 -05:00
Axel Kohlmeyer	fbfb1df5eb	fix typo causing wrong neighbor list copy selections	2017-01-19 20:47:10 -05:00
Steve Plimpton	9a299875da	simplified neighbor list copying to avoid possible same-timestep re-build issues	2017-01-19 17:01:15 -07:00
Stan Moore	fc94f1bd18	Fixing GPU memory issues in Kokkos	2017-01-19 12:14:25 -07:00
Stan Moore	5ce8e2fbae	Fixing GPU memory issue in modify_kokkos, need to cherry pick back to Master	2017-01-19 12:13:48 -07:00
Tim Mattox	f6cd98636b	USER-DPD: Also apply "check if a0 is zero" optimization to pair_dpd_fdt This relates to commit `4eb08a5822` that was applied to pair_dpd_fdt_energy	2017-01-18 16:17:11 -05:00
Tim Mattox	05cafb716f	USER-DPD: cleanup initialization of splitFDT_flag in pair_dpd_fdt.cpp	2017-01-18 15:51:50 -05:00
sjplimp	3af4b3c28c	Merge pull request #337 from ohenrich/user-cgdna Added source code and documentation for USER-CGDNA	2017-01-18 11:31:35 -07:00
sjplimp	7fc0970587	Merge pull request #344 from timattox/USER-DPD_zero_compute USER-DPD: Skip a0*stuff computations, if a0 was set to zero in pair_coeff	2017-01-18 11:31:14 -07:00
sjplimp	93262b52b4	Merge pull request #343 from timattox/USER-DPD_bugfix_molecule USER-DPD: bugfix for a segfault when using MOLECULE and DPD together.	2017-01-18 11:30:58 -07:00
Tim Mattox	4eb08a5822	USER-DPD: Skip a0*stuff computations, if a0 was set to zero in pair_coeff. This saves around 10% of the runtime for many of our tests using SSA.	2017-01-17 15:55:39 -05:00
Tim Mattox	01609f55e2	USER-DPD: bugfix for a segfault when using MOLECULE and DPD together.	2017-01-17 12:47:59 -05:00
Axel Kohlmeyer	8360e70f4e	update USER-CGDNA examples to follow LAMMPS style	2017-01-13 18:56:45 -05:00
Axel Kohlmeyer	b988b29413	remove dead code	2017-01-13 18:43:35 -05:00
Axel Kohlmeyer	5d48bfdcab	USER-CGDNA whitespace cleanup: expand tabs and remove trailing whitespace	2017-01-13 18:40:34 -05:00
Axel Kohlmeyer	fe8caa8a56	apply some LAMMPS formatting style conventions for include files	2017-01-13 18:33:32 -05:00
Axel Kohlmeyer	afaacc6173	add USER-CGDNA package with dependencies into the build system	2017-01-13 18:32:32 -05:00
Axel Kohlmeyer	374abea0f0	some minor documentation integration tweaks for USER-CGDNA package	2017-01-13 18:09:45 -05:00
Oliver Henrich	96259ea2d2	Added source code and documentation for USER-CGDNA	2017-01-13 13:36:54 +00:00