add neighbor list kind output to screen

Merge pull request #349 from akohlmey/collected-small-fixes
collected fixes and improvements
2017-01-20 13:24:09 -07:00 · 2017-01-20 13:19:43 -07:00 · 2017-01-20 13:15:16 -07:00 · 2017-01-20 13:14:59 -07:00 · 2017-01-20 14:30:46 -05:00 · 2017-01-19 20:47:10 -05:00
77 changed files with 13631 additions and 227 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="20 Jan 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
+20 Jan 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,
@ -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_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/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_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/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_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/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/KOKKOS/domain_kokkos.cpp
+++ b/src/KOKKOS/domain_kokkos.cpp
@ -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);
 }

 /* ----------------------------------------------------------------------
--- 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/verlet_kokkos.cpp
+++ b/src/KOKKOS/verlet_kokkos.cpp
@ -170,7 +170,7 @@ void VerletKokkos::setup()

  modify->setup(vflag);
  output->setup();
-  lmp->kokkos->auto_sync = 0;
+  lmp->kokkos->auto_sync = 1;
  update->setupflag = 1;
 }

--- a/src/MC/fix_atom_swap.cpp
+++ b/src/MC/fix_atom_swap.cpp
@ -57,7 +57,10 @@ using namespace MathConst;

 FixAtomSwap::FixAtomSwap(LAMMPS *lmp, int narg, char **arg) :
  Fix(lmp, narg, arg),
-  idregion(NULL), type_list(NULL), mu(NULL), qtype(NULL), sqrt_mass_ratio(NULL), local_swap_iatom_list(NULL), local_swap_jatom_list(NULL), local_swap_atom_list(NULL), random_equal(NULL), random_unequal(NULL), c_pe(NULL)
+  idregion(NULL), type_list(NULL), mu(NULL), qtype(NULL), 
+  sqrt_mass_ratio(NULL), local_swap_iatom_list(NULL), 
+  local_swap_jatom_list(NULL), local_swap_atom_list(NULL), 
+  random_equal(NULL), random_unequal(NULL), c_pe(NULL)
 {
  if (narg < 10) error->all(FLERR,"Illegal fix atom/swap command");

--- a/src/Makefile
+++ b/src/Makefile
@ -48,7 +48,7 @@ PACKAGE = asphere body class2 colloid compress coreshell dipole gpu \
 	  mpiio mscg opt peri poems \
 	  python qeq reax replica rigid shock snap srd voronoi

-PACKUSER = user-atc user-awpmd user-cg-cmm user-colvars \
+PACKUSER = user-atc user-awpmd user-cg-cmm user-cgdna user-colvars \
 	   user-diffraction user-dpd user-drude user-eff user-fep user-h5md \
 	   user-intel user-lb user-manifold user-mgpt user-misc user-molfile \
 	   user-nc-dump user-omp user-phonon user-qmmm user-qtb \
--- a/src/USER-CGDNA/Install.sh
+++ b/src/USER-CGDNA/Install.sh
@ -0,0 +1,47 @@
+# Install/unInstall package files in LAMMPS
+# mode = 0/1/2 for uninstall/install/update
+
+mode=$1
+
+# arg1 = file, arg2 = file it depends on
+
+# enforce using portable C locale
+LC_ALL=C
+export LC_ALL
+
+action () {
+  if (test $mode = 0) then
+    rm -f ../$1
+  elif (! cmp -s $1 ../$1) then
+    if (test -z "$2" || test -e ../$2) then
+      cp $1 ..
+      if (test $mode = 2) then
+        echo "  updating src/$1"
+      fi
+    fi
+  elif (test -n "$2") then
+    if (test ! -e ../$2) then
+      rm -f ../$1
+    fi
+  fi
+}
+
+# list of files with dependcies
+
+action  bond_oxdna_fene.cpp bond_fene.h
+action  bond_oxdna_fene.h bond_fene.h
+action  fix_nve_dotc_langevin.cpp atom_vec_ellipsoid.h
+action  fix_nve_dotc_langevin.h atom_vec_ellipsoid.h
+action  fix_nve_dot.cpp atom_vec_ellipsoid.h
+action  fix_nve_dot.h atom_vec_ellipsoid.h
+action  mf_oxdna.h atom_vec_ellipsoid.h
+action  pair_oxdna_coaxstk.cpp atom_vec_ellipsoid.h
+action  pair_oxdna_coaxstk.h atom_vec_ellipsoid.h
+action  pair_oxdna_excv.cpp atom_vec_ellipsoid.h
+action  pair_oxdna_excv.h atom_vec_ellipsoid.h
+action  pair_oxdna_hbond.cpp atom_vec_ellipsoid.h
+action  pair_oxdna_hbond.h atom_vec_ellipsoid.h
+action  pair_oxdna_stk.cpp atom_vec_ellipsoid.h
+action  pair_oxdna_stk.h atom_vec_ellipsoid.h
+action  pair_oxdna_xstk.cpp atom_vec_ellipsoid.h
+action  pair_oxdna_xstk.h atom_vec_ellipsoid.h
--- a/src/USER-CGDNA/README
+++ b/src/USER-CGDNA/README
@ -0,0 +1,67 @@
+This package contains a LAMMPS implementation of coarse-grained
+models of DNA, which can be used to model sequence-specific
+DNA strands.
+
+See the doc pages and [1,2] for the individual bond and pair styles. 
+The packages contains also a new Langevin-type rigid-body integrator,
+which has also its own doc page and is explained in [3].
+
+[1] T. Ouldridge, A. Louis, J. Doye, "Structural, mechanical, 
+and thermodynamic properties of a coarse-grained DNA model",
+J. Chem. Phys. 134, 085101 (2011).
+
+[2] T.E. Ouldridge, Coarse-grained modelling of DNA and DNA 
+self-assembly, DPhil. University of Oxford (2011).
+
+[3] R. Davidchack, T. Ouldridge, M. Tretyakov, "New Langevin and 
+gradient thermostats for rigid body dynamics", J. Chem. Phys. 142, 
+144114 (2015).
+
+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 as
+well as DNA duplexes and arrays of duplexes can be found in
+/examples/USER/cgdna/util/.  A technical report with more information
+on the model, the structure of the input and data file, the setup tool
+and the performance of the LAMMPS-implementation of oxDNA can be found
+in /doc/src/PDF/USER-CGDNA-overview.pdf.
+
+IMPORTANT NOTE: This package can only be used if LAMMPS is compiled
+with the MOLECULE and ASPHERE packages.  These should be included in
+the LAMMPS build by typing "make yes-asphere yes-molecule" prior to
+the usual compilation (see the "Including/excluding packages" section
+of the LAMMPS manual).
+
+The creator of this package is:
+
+Dr Oliver Henrich
+University of Edinburgh, UK
+ohenrich@ph.ed.ac.uk
+o.henrich@epcc.ed.ac.uk
+
+--------------------------------------------------------------------------
+
+** Bond styles provided by this package:
+
+bond_oxdna_fene.cpp:  backbone connectivity, a modified FENE potential
+
+** Pair styles provided by this package:
+
+pair_oxdna_excv.cpp:  excluded volume interaction between the nucleotides
+
+pair_oxdna_stk.cpp:  stacking interaction between consecutive nucleotides
+                     on the same strand
+
+pair_oxdna_hbond.cpp:  hydrogen-bonding interaction between complementary
+                       nucleotides on different strands, e.g. A-T and C-G
+
+pair_oxdna_xstk.cpp:  cross-stacking interaction between nucleotides
+
+pair_oxdna_coaxstk.cpp:  coaxial stacking interaction between nucleotides
+
+** Fixes provided by this package:
+
+fix_nve_dotc_langevin.cpp:  fix for Langevin-type rigid body integrator "C"
+                            in above Ref. [3] 
+
+fix_nve_dot.cpp:  NVE-type rigid body integrator without noise
--- a/src/USER-CGDNA/bond_oxdna_fene.cpp
+++ b/src/USER-CGDNA/bond_oxdna_fene.cpp
@ -0,0 +1,335 @@
+/* ----------------------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   http://lammps.sandia.gov, Sandia National Laboratories
+   Steve Plimpton, sjplimp@sandia.gov
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+/* ----------------------------------------------------------------------
+   Contributing author: Oliver Henrich (EPCC, University of Edinburgh)
+------------------------------------------------------------------------- */
+
+#include <math.h>
+#include <stdlib.h>
+#include "bond_oxdna_fene.h"
+#include "atom.h"
+#include "neighbor.h"
+#include "domain.h"
+#include "comm.h"
+#include "update.h"
+#include "force.h"
+#include "memory.h"
+#include "error.h"
+#include "atom_vec_ellipsoid.h"
+#include "math_extra.h"
+
+using namespace LAMMPS_NS;
+
+/* ---------------------------------------------------------------------- */
+
+BondOxdnaFene::BondOxdnaFene(LAMMPS *lmp) : Bond(lmp)
+{
+
+}
+
+/* ---------------------------------------------------------------------- */
+
+BondOxdnaFene::~BondOxdnaFene()
+{
+  if (allocated) {
+
+    memory->destroy(setflag);
+    memory->destroy(k);
+    memory->destroy(Delta);
+    memory->destroy(r0);
+
+  }
+}
+
+/* ----------------------------------------------------------------------
+   compute function for oxDNA FENE-bond interaction
+   s=sugar-phosphate backbone site, b=base site, st=stacking site
+------------------------------------------------------------------------- */
+void BondOxdnaFene::compute(int eflag, int vflag)
+{
+  int a,b,in,type;
+  double delf[3],delta[3],deltb[3]; // force, torque increment;;
+  double delr[3],ebond,fbond;
+  double rsq,Deltasq,rlogarg;
+  double r,rr0,rr0sq;
+  // distances COM-backbone site
+  double d_cs=-0.24;
+  // vectors COM-backbone site in lab frame
+  double ra_cs[3],rb_cs[3];
+
+  double *qa,ax[3],ay[3],az[3];
+  double *qb,bx[3],by[3],bz[3];
+
+  double **x = atom->x;
+  double **f = atom->f;
+  double **torque = atom->torque;
+
+  AtomVecEllipsoid *avec = (AtomVecEllipsoid *) atom->style_match("ellipsoid");
+  AtomVecEllipsoid::Bonus *bonus = avec->bonus;
+
+  int **bondlist = neighbor->bondlist;
+  int nbondlist = neighbor->nbondlist;
+  int nlocal = atom->nlocal;
+  int newton_bond = force->newton_bond;
+
+  ebond = 0.0;
+  if (eflag || vflag) ev_setup(eflag,vflag);
+  else evflag = 0;
+
+  // loop over FENE bonds
+
+  for (in = 0; in < nbondlist; in++) {
+
+    a = bondlist[in][1];
+    b = bondlist[in][0];
+    type = bondlist[in][2];
+
+    qa=bonus[a].quat;
+    MathExtra::q_to_exyz(qa,ax,ay,az);
+    qb=bonus[b].quat;
+    MathExtra::q_to_exyz(qb,bx,by,bz);
+
+    // vector COM-backbone site a and b
+    ra_cs[0] = d_cs*ax[0];
+    ra_cs[1] = d_cs*ax[1];
+    ra_cs[2] = d_cs*ax[2];
+    rb_cs[0] = d_cs*bx[0];
+    rb_cs[1] = d_cs*bx[1];
+    rb_cs[2] = d_cs*bx[2];
+
+    // vector backbone site b to a
+    delr[0] = x[a][0] + ra_cs[0] - x[b][0] - rb_cs[0];
+    delr[1] = x[a][1] + ra_cs[1] - x[b][1] - rb_cs[1];
+    delr[2] = x[a][2] + ra_cs[2] - x[b][2] - rb_cs[2];
+    rsq = delr[0]*delr[0] + delr[1]*delr[1] + delr[2]*delr[2];
+    r = sqrt(rsq);
+
+    rr0 = r - r0[type];
+    rr0sq = rr0*rr0;
+    Deltasq = Delta[type] * Delta[type];
+    rlogarg = 1.0 - rr0sq/Deltasq;
+
+    // if r -> Delta, then rlogarg < 0.0 which is an error
+    // issue a warning and reset rlogarg = epsilon
+    // if r > 2*Delta something serious is wrong, abort
+
+    if (rlogarg < 0.1) {
+      char str[128];
+      sprintf(str,"FENE bond too long: " BIGINT_FORMAT " "
+              TAGINT_FORMAT " " TAGINT_FORMAT " %g",
+              update->ntimestep,atom->tag[a],atom->tag[b],r);
+      error->warning(FLERR,str,0);
+      if (rlogarg <= -3.0) error->one(FLERR,"Bad FENE bond");
+    }
+
+    fbond = -k[type]*rr0/rlogarg/Deltasq/r;
+    delf[0] = delr[0]*fbond;
+    delf[1] = delr[1]*fbond;
+    delf[2] = delr[2]*fbond;
+
+    // energy
+
+    if (eflag) {
+      ebond = -0.5 * k[type]*log(rlogarg);
+    }
+
+    // apply force and torque to each of 2 atoms
+
+    if (newton_bond || a < nlocal) {
+
+      f[a][0] += delf[0];
+      f[a][1] += delf[1];
+      f[a][2] += delf[2];
+
+      MathExtra::cross3(ra_cs,delf,delta);
+
+      torque[a][0] += delta[0];
+      torque[a][1] += delta[1];
+      torque[a][2] += delta[2];
+
+    }
+
+    if (newton_bond || b < nlocal) {
+
+      f[b][0] -= delf[0];
+      f[b][1] -= delf[1];
+      f[b][2] -= delf[2];
+
+      MathExtra::cross3(rb_cs,delf,deltb);
+
+      torque[b][0] -= deltb[0];
+      torque[b][1] -= deltb[1];
+      torque[b][2] -= deltb[2];
+
+    }
+
+    // increment energy and virial
+    if (evflag) ev_tally(a,b,nlocal,newton_bond,ebond,fbond,delr[0],delr[1],delr[2]);
+
+  }
+
+}
+
+/* ---------------------------------------------------------------------- */
+
+void BondOxdnaFene::allocate()
+{
+  allocated = 1;
+  int n = atom->nbondtypes;
+
+  memory->create(k,n+1,"bond:k");
+  memory->create(Delta,n+1,"bond:Delta");
+  memory->create(r0,n+1,"bond:r0");
+  memory->create(setflag,n+1,"bond:setflag");
+
+  for (int i = 1; i <= n; i++) setflag[i] = 0;
+
+}
+
+/* ----------------------------------------------------------------------
+   set coeffs for one type
+------------------------------------------------------------------------- */
+
+void BondOxdnaFene::coeff(int narg, char **arg)
+{
+  if (narg != 4) error->all(FLERR,"Incorrect args for bond coefficients in oxdna_fene");
+  if (!allocated) allocate();
+
+  int ilo,ihi;
+  force->bounds(FLERR,arg[0],atom->nbondtypes,ilo,ihi);
+
+  double k_one = force->numeric(FLERR,arg[1]);
+  double Delta_one = force->numeric(FLERR,arg[2]);
+  double r0_one = force->numeric(FLERR,arg[3]);
+
+  int count = 0;
+
+  for (int i = ilo; i <= ihi; i++) {
+    k[i] = k_one;
+    Delta[i] = Delta_one;
+    r0[i] = r0_one;
+    setflag[i] = 1;
+    count++;
+  }
+
+  if (count == 0) error->all(FLERR,"Incorrect args for bond coefficients in oxdna_fene");
+
+}
+
+/* ----------------------------------------------------------------------
+   set special_bond settings and check if valid
+------------------------------------------------------------------------- */
+
+void BondOxdnaFene::init_style()
+{
+  /* special bonds have to be lj = 0 1 1 and coul = 1 1 1 to exclude
+     the ss excluded volume interaction between nearest neighbours   */
+
+  force->special_lj[1] = 0.0;
+  force->special_lj[2] = 1.0;
+  force->special_lj[3] = 1.0;
+  force->special_coul[1] = 1.0;
+  force->special_coul[2] = 1.0;
+  force->special_coul[3] = 1.0;
+
+  fprintf(screen,"Finding 1-2 1-3 1-4 neighbors ...\n"
+                 " Special bond factors lj:   %-10g %-10g %-10g\n"
+                 " Special bond factors coul: %-10g %-10g %-10g\n",
+                 force->special_lj[1],force->special_lj[2],force->special_lj[3],
+                 force->special_coul[1],force->special_coul[2],force->special_coul[3]);
+
+  if (force->special_lj[1] != 0.0 || force->special_lj[2] != 1.0 || force->special_lj[3] != 1.0 ||
+      force->special_coul[1] != 1.0 || force->special_coul[2] != 1.0 || force->special_coul[3] != 1.0)
+  {
+    if (comm->me == 0)
+      error->warning(FLERR,"Use special bonds lj = 0,1,1 and coul = 1,1,1 with bond style oxdna_fene");
+  }
+
+}
+
+/* ---------------------------------------------------------------------- */
+
+double BondOxdnaFene::equilibrium_distance(int i)
+{
+  return r0[i];
+}
+
+/* ----------------------------------------------------------------------
+   proc 0 writes to restart file
+------------------------------------------------------------------------- */
+
+void BondOxdnaFene::write_restart(FILE *fp)
+{
+  fwrite(&k[1],sizeof(double),atom->nbondtypes,fp);
+  fwrite(&Delta[1],sizeof(double),atom->nbondtypes,fp);
+  fwrite(&r0[1],sizeof(double),atom->nbondtypes,fp);
+}
+
+/* ----------------------------------------------------------------------
+   proc 0 reads from restart file, bcasts
+------------------------------------------------------------------------- */
+
+void BondOxdnaFene::read_restart(FILE *fp)
+{
+  allocate();
+
+  if (comm->me == 0) {
+    fread(&k[1],sizeof(double),atom->nbondtypes,fp);
+    fread(&Delta[1],sizeof(double),atom->nbondtypes,fp);
+    fread(&r0[1],sizeof(double),atom->nbondtypes,fp);
+  }
+  MPI_Bcast(&k[1],atom->nbondtypes,MPI_DOUBLE,0,world);
+  MPI_Bcast(&Delta[1],atom->nbondtypes,MPI_DOUBLE,0,world);
+  MPI_Bcast(&r0[1],atom->nbondtypes,MPI_DOUBLE,0,world);
+
+  for (int i = 1; i <= atom->nbondtypes; i++) setflag[i] = 1;
+}
+
+/* ----------------------------------------------------------------------
+   proc 0 writes to data file
+------------------------------------------------------------------------- */
+
+void BondOxdnaFene::write_data(FILE *fp)
+{
+  for (int i = 1; i <= atom->nbondtypes; i++)
+    fprintf(fp,"%d %g %g %g\n",i,k[i],r0[i],Delta[i]);
+}
+
+/* ---------------------------------------------------------------------- */
+
+double BondOxdnaFene::single(int type, double rsq, int i, int j,
+                        double &fforce)
+{
+  double r = sqrt(rsq);
+  double rr0 = r - r0[type];
+  double rr0sq = rr0*rr0;
+  double Deltasq = Delta[type] * Delta[type];
+  double rlogarg = 1.0 - rr0sq/Deltasq;
+
+  // if r -> Delta, then rlogarg < 0.0 which is an error
+  // issue a warning and reset rlogarg = epsilon
+  // if r > 2*Delta something serious is wrong, abort
+
+  if (rlogarg < 0.1) {
+    char str[128];
+    sprintf(str,"FENE bond too long: " BIGINT_FORMAT " %g",
+            update->ntimestep,sqrt(rsq));
+    error->warning(FLERR,str,0);
+    if (rlogarg <= -3.0) error->one(FLERR,"Bad FENE bond");
+  }
+
+  double eng = -0.5 * k[type]*log(rlogarg);
+  fforce = -k[type]*rr0/rlogarg/Deltasq/r;
+
+  return eng;
+}
--- a/src/USER-CGDNA/bond_oxdna_fene.h
+++ b/src/USER-CGDNA/bond_oxdna_fene.h
@ -0,0 +1,79 @@
+/* -*- c++ -*- ----------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   http://lammps.sandia.gov, Sandia National Laboratories
+   Steve Plimpton, sjplimp@sandia.gov
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+/* ----------------------------------------------------------------------
+   Contributing author: Oliver Henrich (EPCC, University of Edinburgh)
+------------------------------------------------------------------------- */
+
+#ifdef BOND_CLASS
+
+BondStyle(oxdna_fene,BondOxdnaFene)
+
+#else
+
+#ifndef LMP_BOND_OXDNA_FENE_H
+#define LMP_BOND_OXDNA_FENE_H
+
+#include "bond.h"
+
+namespace LAMMPS_NS {
+
+class BondOxdnaFene : public Bond {
+ public:
+  BondOxdnaFene(class LAMMPS *);
+  virtual ~BondOxdnaFene();
+  virtual void compute(int, int);
+  void coeff(int, char **);
+  void init_style();
+  double equilibrium_distance(int);
+  void write_restart(FILE *);
+  void read_restart(FILE *);
+  void write_data(FILE *);
+  double single(int, double, int, int, double &);
+
+ protected:
+  double *k,*Delta,*r0; // FENE
+
+  void allocate();
+};
+
+}
+
+#endif
+#endif
+
+/* ERROR/WARNING messages:
+
+W: FENE bond too long: %ld %d %d %g
+
+A FENE bond has stretched dangerously far.  It's interaction strength
+will be truncated to attempt to prevent the bond from blowing up.
+
+E: Bad FENE bond
+
+Two atoms in a FENE bond have become so far apart that the bond cannot
+be computed.
+
+E: Incorrect args for bond coefficients
+
+Self-explanatory.  Check the input script or data file.
+
+W: Use special bonds = 0,1,1 with bond style oxdna
+
+Most FENE models need this setting for the special_bonds command.
+
+W: FENE bond too long: %ld %g
+
+A FENE bond has stretched dangerously far.  It's interaction strength
+will be truncated to attempt to prevent the bond from blowing up.
+
+*/
--- a/src/USER-CGDNA/fix_nve_dot.cpp
+++ b/src/USER-CGDNA/fix_nve_dot.cpp
@ -0,0 +1,211 @@
+/* ----------------------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   http://lammps.sandia.gov, Sandia National Laboratories
+   Steve Plimpton, sjplimp@sandia.gov
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+/* ----------------------------------------------------------------------
+   Contributing author: Oliver Henrich (EPCC, University of Edinburgh)
+------------------------------------------------------------------------- */
+
+#include <math.h>
+#include <stdio.h>
+#include <string.h>
+#include "fix_nve_dot.h"
+#include "math_extra.h"
+#include "atom.h"
+#include "atom_vec_ellipsoid.h"
+#include "force.h"
+#include "update.h"
+#include "memory.h"
+#include "error.h"
+
+using namespace LAMMPS_NS;
+using namespace FixConst;
+using namespace MathExtra;
+
+#define INERTIA 0.2          // moment of inertia prefactor for ellipsoid
+
+/* ---------------------------------------------------------------------- */
+
+FixNVEDot::FixNVEDot(LAMMPS *lmp, int narg, char **arg) :
+  FixNVE(lmp, narg, arg) {}
+
+/* ---------------------------------------------------------------------- */
+
+void FixNVEDot::init()
+{
+  avec = (AtomVecEllipsoid *) atom->style_match("ellipsoid");
+  if (!avec)
+    error->all(FLERR,"Compute nve/dot requires atom style ellipsoid");
+
+  // check that all particles are finite-size ellipsoids
+  // no point particles allowed, spherical is OK
+
+  int *ellipsoid = atom->ellipsoid;
+  int *mask = atom->mask;
+  int nlocal = atom->nlocal;
+
+  for (int i = 0; i < nlocal; i++)
+    if (mask[i] & groupbit)
+      if (ellipsoid[i] < 0)
+        error->one(FLERR,"Fix nve/dot requires extended particles");
+
+  FixNVE::init();
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixNVEDot::initial_integrate(int vflag)
+{
+  double *shape,*quat;
+  double fquat[4],conjqm[4],inertia[3];
+
+  AtomVecEllipsoid::Bonus *bonus = avec->bonus;
+  int *ellipsoid = atom->ellipsoid;
+  double **x = atom->x;
+  double **v = atom->v;
+  double **f = atom->f;
+  double **angmom = atom->angmom;
+  double **torque = atom->torque;
+  double *rmass = atom->rmass;
+  int *mask = atom->mask;
+  int nlocal = atom->nlocal;
+  if (igroup == atom->firstgroup) nlocal = atom->nfirst;
+
+  // set timestep here since dt may have changed or come via rRESPA
+
+  dt = update->dt;
+  dthlf = 0.5 * dt;
+
+  for (int i = 0; i < nlocal; i++)
+    if (mask[i] & groupbit) {
+
+      dthlfm = dthlf / rmass[i];
+      quat = bonus[ellipsoid[i]].quat;
+      shape = bonus[ellipsoid[i]].shape;
+
+      // update momentum by 1/2 step
+      v[i][0] += dthlfm * f[i][0];
+      v[i][1] += dthlfm * f[i][1];
+      v[i][2] += dthlfm * f[i][2];
+
+      // update position by full step
+      x[i][0] += dt * v[i][0];
+      x[i][1] += dt * v[i][1];
+      x[i][2] += dt * v[i][2];
+
+      // convert angular momentum and torque in space frame into
+      // quaternion 4-momentum and 1/2 of 4-torque in body frame
+      vec3_to_vec4(quat,angmom[i],conjqm);
+      conjqm[0] *= 2.0;
+      conjqm[1] *= 2.0;
+      conjqm[2] *= 2.0;
+      conjqm[3] *= 2.0;
+      vec3_to_vec4(quat,torque[i],fquat);
+
+      // update quaternion 4-momentum by 1/2 step
+      conjqm[0] += dt * fquat[0];
+      conjqm[1] += dt * fquat[1];
+      conjqm[2] += dt * fquat[2];
+      conjqm[3] += dt * fquat[3];
+
+      // principal moments of inertia
+      inertia[0] = INERTIA*rmass[i] * (shape[1]*shape[1]+shape[2]*shape[2]);
+      inertia[1] = INERTIA*rmass[i] * (shape[0]*shape[0]+shape[2]*shape[2]);
+      inertia[2] = INERTIA*rmass[i] * (shape[0]*shape[0]+shape[1]*shape[1]);
+
+      // rotate quaternion and quaternion 4-momentum by full step
+      no_squish_rotate(3,conjqm,quat,inertia,dthlf);
+      no_squish_rotate(2,conjqm,quat,inertia,dthlf);
+      no_squish_rotate(1,conjqm,quat,inertia,dt);
+      no_squish_rotate(2,conjqm,quat,inertia,dthlf);
+      no_squish_rotate(3,conjqm,quat,inertia,dthlf);
+
+      qnormalize(quat);
+
+      // convert quaternion 4-momentum in body frame back to angular momentum in space frame
+      vec4_to_vec3(quat,conjqm,angmom[i]);
+
+      angmom[i][0] *= 0.5;
+      angmom[i][1] *= 0.5;
+      angmom[i][2] *= 0.5;
+
+    }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixNVEDot::final_integrate()
+{
+
+  double *quat;
+  double fquat[4],conjqm[4];
+  double conjqm_dot_quat;
+
+  AtomVecEllipsoid::Bonus *bonus = avec->bonus;
+  int *ellipsoid = atom->ellipsoid;
+  double **v = atom->v;
+  double **f = atom->f;
+  double **angmom = atom->angmom;
+  double **torque = atom->torque;
+  double *rmass = atom->rmass;
+  int *mask = atom->mask;
+  int nlocal = atom->nlocal;
+  if (igroup == atom->firstgroup) nlocal = atom->nfirst;
+
+  // set timestep here since dt may have changed or come via rRESPA
+
+  dt = update->dt;
+  dthlf = 0.5 * dt;
+
+  for (int i = 0; i < nlocal; i++)
+    if (mask[i] & groupbit) {
+
+      dthlfm = dthlf / rmass[i];
+      quat = bonus[ellipsoid[i]].quat;
+
+      // update momentum
+      v[i][0] += dthlfm * f[i][0];
+      v[i][1] += dthlfm * f[i][1];
+      v[i][2] += dthlfm * f[i][2];
+
+      // convert angular momentum and torque in space frame into
+      // quaternion 4-momentum and 1/2 of 4-torque in body frame
+      vec3_to_vec4(quat,angmom[i],conjqm);
+      conjqm[0] *= 2.0;
+      conjqm[1] *= 2.0;
+      conjqm[2] *= 2.0;
+      conjqm[3] *= 2.0;
+      vec3_to_vec4(quat,torque[i],fquat);
+
+      // update quaternion 4-momentum by 1/2 step
+      conjqm[0] += dt * fquat[0];
+      conjqm[1] += dt * fquat[1];
+      conjqm[2] += dt * fquat[2];
+      conjqm[3] += dt * fquat[3];
+
+      // subtract component parallel to quaternion for improved numerical accuracy
+      conjqm_dot_quat = conjqm[0]*quat[0] + conjqm[1]*quat[1] + conjqm[2]*quat[2] + conjqm[3]*quat[3];
+
+      conjqm[0] -= conjqm_dot_quat * quat[0];
+      conjqm[1] -= conjqm_dot_quat * quat[1];
+      conjqm[2] -= conjqm_dot_quat * quat[2];
+      conjqm[3] -= conjqm_dot_quat * quat[3];
+
+      // convert quaternion 4-momentum in body frame back to angular momentum in space frame
+      vec4_to_vec3(quat,conjqm,angmom[i]);
+
+      angmom[i][0] *= 0.5;
+      angmom[i][1] *= 0.5;
+      angmom[i][2] *= 0.5;
+
+    }
+}
--- a/src/USER-CGDNA/fix_nve_dot.h
+++ b/src/USER-CGDNA/fix_nve_dot.h
@ -0,0 +1,68 @@
+/* -*- c++ -*- ----------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   http://lammps.sandia.gov, Sandia National Laboratories
+   Steve Plimpton, sjplimp@sandia.gov
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#ifdef FIX_CLASS
+
+FixStyle(nve/dot,FixNVEDot)
+
+#else
+
+#ifndef LMP_FIX_NVE_DOT_H
+#define LMP_FIX_NVE_DOT_H
+
+#include "fix_nve.h"
+
+namespace LAMMPS_NS {
+
+class FixNVEDot : public FixNVE {
+ public:
+  FixNVEDot(class LAMMPS *, int, char **);
+  void init();
+  void initial_integrate(int);
+  void final_integrate();
+
+ private:
+  double dt,dthlf,dthlfm;
+  class AtomVecEllipsoid *avec;
+  // conversion from 3-vector in space frame to 4-vector in body frame
+  inline void vec3_to_vec4(const double * q, const double * v3, double * v4)
+  {
+    v4[0] = -q[1]*v3[0] - q[2]*v3[1] - q[3]*v3[2];
+    v4[1] =  q[0]*v3[0] + q[3]*v3[1] - q[2]*v3[2];
+    v4[2] = -q[3]*v3[0] + q[0]*v3[1] + q[1]*v3[2];
+    v4[3] =  q[2]*v3[0] - q[1]*v3[1] + q[0]*v3[2];
+  }
+  // conversion from 4-vector in body frame to 3-vector in space frame
+  inline void vec4_to_vec3(const double * q, const double * v4, double * v3)
+  {
+    v3[0] = -q[1]*v4[0] + q[0]*v4[1] - q[3]*v4[2] + q[2]*v4[3];
+    v3[1] = -q[2]*v4[0] + q[3]*v4[1] + q[0]*v4[2] - q[1]*v4[3];
+    v3[2] = -q[3]*v4[0] - q[2]*v4[1] + q[1]*v4[2] + q[0]*v4[3];
+  }
+};
+
+}
+#endif
+#endif
+
+/* ERROR/WARNING messages:
+
+E: Compute nve/dot requires atom style ellipsoid
+
+Self-explanatory.
+
+E: Fix nve/dot requires extended particles
+
+This fix can only be used for particles with a shape setting.
+
+*/
--- a/src/USER-CGDNA/fix_nve_dotc_langevin.cpp
+++ b/src/USER-CGDNA/fix_nve_dotc_langevin.cpp
@ -0,0 +1,329 @@
+/* ----------------------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   http://lammps.sandia.gov, Sandia National Laboratories
+   Steve Plimpton, sjplimp@sandia.gov
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+/* ----------------------------------------------------------------------
+   Contributing author: Oliver Henrich (EPCC, University of Edinburgh)
+------------------------------------------------------------------------- */
+
+#include <math.h>
+#include <stdio.h>
+#include <string.h>
+#include "fix_nve_dotc_langevin.h"
+#include "math_extra.h"
+#include "atom.h"
+#include "atom_vec_ellipsoid.h"
+#include "force.h"
+#include "update.h"
+#include "comm.h"
+#include "random_mars.h"
+#include "memory.h"
+#include "error.h"
+
+using namespace LAMMPS_NS;
+using namespace FixConst;
+using namespace MathExtra;
+
+#define INERTIA 0.2          // moment of inertia prefactor for ellipsoid
+
+/* ---------------------------------------------------------------------- */
+
+FixNVEDotcLangevin::FixNVEDotcLangevin(LAMMPS *lmp, int narg, char **arg) :
+  FixNVE(lmp, narg, arg)
+{
+  if (narg != 9) error->all(FLERR,"Illegal fix nve/dotc/langevin command");
+
+  t_start = force->numeric(FLERR,arg[3]);
+  t_target = t_start;
+  t_stop = force->numeric(FLERR,arg[4]);
+  t_period = force->numeric(FLERR,arg[5]);
+  if (t_period <= 0.0) error->all(FLERR,"Fix nve/dotc/langevin period must be > 0.0");
+  gamma = 1.0/t_period;
+  seed = force->inumeric(FLERR,arg[6]);
+  if (seed <= 0) error->all(FLERR,"Illegal fix nve/dotc/langevin command");
+
+  if (strcmp(arg[7],"angmom") == 0) {
+    if (9 > narg) error->all(FLERR,"Illegal fix nve/dotc/langevin command");
+    if (strcmp(arg[8],"no") == 0) {
+      ascale = 0.0;
+      Gamma = 0.0;
+    }
+    else {
+      ascale = force->numeric(FLERR,arg[8]);
+      Gamma = gamma * ascale;
+    }
+
+  }
+
+  // initialize Marsaglia RNG with processor-unique seed
+
+  random = new RanMars(lmp,seed + comm->me);
+
+}
+
+/* ---------------------------------------------------------------------- */
+
+FixNVEDotcLangevin::~FixNVEDotcLangevin()
+{
+
+  delete random;
+
+}
+
+
+/* ---------------------------------------------------------------------- */
+
+void FixNVEDotcLangevin::init()
+{
+
+  int *ellipsoid = atom->ellipsoid;
+  int *mask = atom->mask;
+  int nlocal = atom->nlocal;
+
+  avec = (AtomVecEllipsoid *) atom->style_match("ellipsoid");
+
+  if (!avec)
+    error->all(FLERR,"Fix nve/dotc/langevin requires atom style ellipsoid");
+
+  for (int i = 0; i < nlocal; i++)
+    if (mask[i] & groupbit)
+      if (ellipsoid[i] < 0)
+        error->one(FLERR,"Fix nve/dotc/langevin requires extended particles");
+
+
+  // set prefactor
+  gfactor1 = exp(-gamma*update->dt);
+
+  // set square root of temperature
+  compute_target();
+
+  FixNVE::init();
+}
+
+/* ----------------------------------------------------------------------
+   set current t_target and t_sqrt
+------------------------------------------------------------------------- */
+
+void FixNVEDotcLangevin::compute_target()
+{
+  double delta = update->ntimestep - update->beginstep;
+  if (delta != 0.0) delta /= update->endstep - update->beginstep;
+
+  // Only homogeneous temperature supported
+  t_target = t_start + delta * (t_stop-t_start);
+  tsqrt = sqrt(t_target);
+
+}
+
+
+/* ---------------------------------------------------------------------- */
+
+void FixNVEDotcLangevin::initial_integrate(int vflag)
+{
+  double *shape,*quat;
+  double fquat[4],conjqm[4],inertia[3];
+  double slq_conjqm[3];
+
+  AtomVecEllipsoid::Bonus *bonus = avec->bonus;
+  int *ellipsoid = atom->ellipsoid;
+  double **x = atom->x;
+  double **v = atom->v;
+  double **f = atom->f;
+  double **angmom = atom->angmom;
+  double **torque = atom->torque;
+  double *rmass = atom->rmass;
+  int *mask = atom->mask;
+  int nlocal = atom->nlocal;
+
+  if (igroup == atom->firstgroup) nlocal = atom->nfirst;
+
+  // set timestep here since dt may have changed or come via rRESPA
+
+  dt = update->dt;
+  dthlf = 0.5 * dt;
+  dtqrt = 0.25 * dt;
+
+  // set square root of temperature
+  compute_target();
+
+  for (int i = 0; i < nlocal; i++)
+    if (mask[i] & groupbit) {
+
+      dthlfm = dthlf / rmass[i];
+      quat = bonus[ellipsoid[i]].quat;
+      shape = bonus[ellipsoid[i]].shape;
+
+      // update momentum by 1/2 step
+      v[i][0] += dthlfm * f[i][0];
+      v[i][1] += dthlfm * f[i][1];
+      v[i][2] += dthlfm * f[i][2];
+
+      // update position by 1/2 step
+      x[i][0] += dthlf * v[i][0];
+      x[i][1] += dthlf * v[i][1];
+      x[i][2] += dthlf * v[i][2];
+
+      // convert angular momentum and torque in space frame into
+      // quaternion 4-momentum and 1/2 of 4-torque in body frame
+      vec3_to_vec4(quat,angmom[i],conjqm);
+      conjqm[0] *= 2.0;
+      conjqm[1] *= 2.0;
+      conjqm[2] *= 2.0;
+      conjqm[3] *= 2.0;
+      vec3_to_vec4(quat,torque[i],fquat);
+
+      // update quaternion 4-momentum by 1/2 step
+      conjqm[0] += dt * fquat[0];
+      conjqm[1] += dt * fquat[1];
+      conjqm[2] += dt * fquat[2];
+      conjqm[3] += dt * fquat[3];
+
+      // principal moments of inertia
+      inertia[0] = INERTIA*rmass[i] * (shape[1]*shape[1]+shape[2]*shape[2]);
+      inertia[1] = INERTIA*rmass[i] * (shape[0]*shape[0]+shape[2]*shape[2]);
+      inertia[2] = INERTIA*rmass[i] * (shape[0]*shape[0]+shape[1]*shape[1]);
+
+      M = inertia[0]*inertia[1]*inertia[2];
+      M /= inertia[1]*inertia[2]+inertia[0]*inertia[2]+inertia[0]*inertia[1];
+
+      // set prefactors
+      // factors 12 and 48 reflect the variance of the uniform distribution:
+      // var = 1/12*(b-a)^2
+      gfactor2 = sqrt(12.0*(1.0-gfactor1*gfactor1)/rmass[i])*tsqrt;
+
+      gfactor3[0] = exp(-Gamma*M*dt/inertia[0]);
+      gfactor3[1] = exp(-Gamma*M*dt/inertia[1]);
+      gfactor3[2] = exp(-Gamma*M*dt/inertia[2]);
+
+      gfactor4[0] = sqrt(48.0*inertia[0]*(1.0-gfactor3[0]*gfactor3[0]))*tsqrt;
+      gfactor4[1] = sqrt(48.0*inertia[1]*(1.0-gfactor3[1]*gfactor3[1]))*tsqrt;
+      gfactor4[2] = sqrt(48.0*inertia[2]*(1.0-gfactor3[2]*gfactor3[2]))*tsqrt;
+
+      // rotate quaternion and quaternion 4-momentum by 1/2 step
+      no_squish_rotate(3,conjqm,quat,inertia,dtqrt);
+      no_squish_rotate(2,conjqm,quat,inertia,dtqrt);
+      no_squish_rotate(1,conjqm,quat,inertia,dthlf);
+      no_squish_rotate(2,conjqm,quat,inertia,dtqrt);
+      no_squish_rotate(3,conjqm,quat,inertia,dtqrt);
+
+      // apply stochastic force to velocities
+      v[i][0] = v[i][0] * gfactor1 + gfactor2 * (random->uniform()-0.5);
+      v[i][1] = v[i][1] * gfactor1 + gfactor2 * (random->uniform()-0.5);
+      v[i][2] = v[i][2] * gfactor1 + gfactor2 * (random->uniform()-0.5);
+
+      // update position by 1/2 step
+      x[i][0] += dthlf * v[i][0];
+      x[i][1] += dthlf * v[i][1];
+      x[i][2] += dthlf * v[i][2];
+
+      // apply stochastic force to quaternion 4-momentum
+      slq_conjqm[0] = -quat[1]*conjqm[0] + quat[0]*conjqm[1] + quat[3]*conjqm[2] - quat[2]*conjqm[3];
+      slq_conjqm[1] = -quat[2]*conjqm[0] - quat[3]*conjqm[1] + quat[0]*conjqm[2] + quat[1]*conjqm[3];
+      slq_conjqm[2] = -quat[3]*conjqm[0] + quat[2]*conjqm[1] - quat[1]*conjqm[2] + quat[0]*conjqm[3];
+
+      gfactor5[0] = gfactor3[0] * slq_conjqm[0] + gfactor4[0] * (random->uniform()-0.5);
+      gfactor5[1] = gfactor3[1] * slq_conjqm[1] + gfactor4[1] * (random->uniform()-0.5);
+      gfactor5[2] = gfactor3[2] * slq_conjqm[2] + gfactor4[2] * (random->uniform()-0.5);
+
+      conjqm[0] = -quat[1] * gfactor5[0] - quat[2] * gfactor5[1] - quat[3] * gfactor5[2];
+      conjqm[1] =  quat[0] * gfactor5[0] - quat[3] * gfactor5[1] + quat[2] * gfactor5[2];
+      conjqm[2] =  quat[3] * gfactor5[0] + quat[0] * gfactor5[1] - quat[1] * gfactor5[2];
+      conjqm[3] = -quat[2] * gfactor5[0] + quat[1] * gfactor5[1] + quat[0] * gfactor5[2];
+
+      // rotate quaternion and quaternion 4-momentum by 1/2 step
+      no_squish_rotate(3,conjqm,quat,inertia,dtqrt);
+      no_squish_rotate(2,conjqm,quat,inertia,dtqrt);
+      no_squish_rotate(1,conjqm,quat,inertia,dthlf);
+      no_squish_rotate(2,conjqm,quat,inertia,dtqrt);
+      no_squish_rotate(3,conjqm,quat,inertia,dtqrt);
+      qnormalize(quat);
+
+      // convert quaternion 4-momentum in body frame back to angular momentum in space frame
+      vec4_to_vec3(quat,conjqm,angmom[i]);
+
+      angmom[i][0] *= 0.5;
+      angmom[i][1] *= 0.5;
+      angmom[i][2] *= 0.5;
+
+    }
+
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixNVEDotcLangevin::final_integrate()
+{
+
+  double *quat;
+  double fquat[4],conjqm[4];
+  double conjqm_dot_quat;
+
+  AtomVecEllipsoid::Bonus *bonus = avec->bonus;
+  int *ellipsoid = atom->ellipsoid;
+  double **v = atom->v;
+  double **f = atom->f;
+  double **angmom = atom->angmom;
+  double **torque = atom->torque;
+  double *rmass = atom->rmass;
+  int *mask = atom->mask;
+  int nlocal = atom->nlocal;
+  if (igroup == atom->firstgroup) nlocal = atom->nfirst;
+
+  // set timestep here since dt may have changed or come via rRESPA
+
+  dt = update->dt;
+  dthlf = 0.5 * dt;
+
+  for (int i = 0; i < nlocal; i++)
+    if (mask[i] & groupbit) {
+
+      dthlfm = dthlf / rmass[i];
+      quat = bonus[ellipsoid[i]].quat;
+
+      // update momentum by 1/2 step
+      v[i][0] += dthlfm * f[i][0];
+      v[i][1] += dthlfm * f[i][1];
+      v[i][2] += dthlfm * f[i][2];
+
+      // convert angular momentum and torque in space frame into
+      // quaternion 4-momentum and 1/2 of 4-torque in body frame
+      vec3_to_vec4(quat,angmom[i],conjqm);
+      conjqm[0] *= 2.0;
+      conjqm[1] *= 2.0;
+      conjqm[2] *= 2.0;
+      conjqm[3] *= 2.0;
+      vec3_to_vec4(quat,torque[i],fquat);
+
+      // update quaternion 4-momentum by 1/2 step
+      conjqm[0] += dt * fquat[0];
+      conjqm[1] += dt * fquat[1];
+      conjqm[2] += dt * fquat[2];
+      conjqm[3] += dt * fquat[3];
+
+      // subtract component parallel to quaternion for improved numerical accuracy
+
+      conjqm_dot_quat = conjqm[0]*quat[0] + conjqm[1]*quat[1] + conjqm[2]*quat[2] + conjqm[3]*quat[3];
+
+      conjqm[0] -= conjqm_dot_quat * quat[0];
+      conjqm[1] -= conjqm_dot_quat * quat[1];
+      conjqm[2] -= conjqm_dot_quat * quat[2];
+      conjqm[3] -= conjqm_dot_quat * quat[3];
+
+      // convert quaternion 4-momentum in body frame back to angular momentum in space frame
+      vec4_to_vec3(quat,conjqm,angmom[i]);
+
+      angmom[i][0] *= 0.5;
+      angmom[i][1] *= 0.5;
+      angmom[i][2] *= 0.5;
+
+    }
+}
--- a/src/USER-CGDNA/fix_nve_dotc_langevin.h
+++ b/src/USER-CGDNA/fix_nve_dotc_langevin.h
@ -0,0 +1,77 @@
+/* -*- c++ -*- ----------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   http://lammps.sandia.gov, Sandia National Laboratories
+   Steve Plimpton, sjplimp@sandia.gov
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#ifdef FIX_CLASS
+
+FixStyle(nve/dotc/langevin,FixNVEDotcLangevin)
+
+#else
+
+#ifndef LMP_FIX_NVE_DOTC_LANGEVIN_H
+#define LMP_FIX_NVE_DOTC_LANGEVIN_H
+
+#include "fix_nve.h"
+
+namespace LAMMPS_NS {
+
+class FixNVEDotcLangevin : public FixNVE {
+ public:
+  FixNVEDotcLangevin(class LAMMPS *, int, char **);
+  virtual ~FixNVEDotcLangevin();
+  void init();
+  void initial_integrate(int);
+  void final_integrate();
+
+ private:
+  double dt,dthlf,dthlfm,dtqrt;
+  // conversion from 3-vector in space frame to 4-vector in body frame
+  inline void vec3_to_vec4(const double * q, const double * v3, double * v4)
+  {
+    v4[0] = -q[1]*v3[0] - q[2]*v3[1] - q[3]*v3[2];
+    v4[1] =  q[0]*v3[0] + q[3]*v3[1] - q[2]*v3[2];
+    v4[2] = -q[3]*v3[0] + q[0]*v3[1] + q[1]*v3[2];
+    v4[3] =  q[2]*v3[0] - q[1]*v3[1] + q[0]*v3[2];
+  }
+  // conversion from 4-vector in body frame to 3-vector in space frame
+  inline void vec4_to_vec3(const double * q, const double * v4, double * v3)
+  {
+    v3[0] = -q[1]*v4[0] + q[0]*v4[1] - q[3]*v4[2] + q[2]*v4[3];
+    v3[1] = -q[2]*v4[0] + q[3]*v4[1] + q[0]*v4[2] - q[1]*v4[3];
+    v3[2] = -q[3]*v4[0] - q[2]*v4[1] + q[1]*v4[2] + q[0]*v4[3];
+  }
+ protected:
+  int seed;
+  class AtomVecEllipsoid *avec;
+  double t_start,t_stop,t_period,t_target,tsqrt;
+  double gamma,Gamma,ascale;
+  double M,gfactor1,gfactor2;
+  double gfactor3[3],gfactor4[3],gfactor5[3];
+  class RanMars *random;
+  void compute_target();
+};
+
+}
+#endif
+#endif
+
+/* ERROR/WARNING messages:
+
+E: Compute nve/dotc/langevin requires atom style ellipsoid
+
+Self-explanatory.
+
+E: Fix nve/dotc/langevin requires extended particles
+
+This fix can only be used for particles with a shape setting.
+
+*/
--- a/src/USER-CGDNA/mf_oxdna.h
+++ b/src/USER-CGDNA/mf_oxdna.h
@ -0,0 +1,263 @@
+/* ----------------------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   http://lammps.sandia.gov, Sandia National Laboratories
+   Steve Plimpton, sjplimp@sandia.gov
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+/* ----------------------------------------------------------------------
+   Contributing author: Oliver Henrich (EPCC, University of Edinburgh)
+------------------------------------------------------------------------- */
+
+#ifndef MF_OXDNA_H
+#define MF_OXDNA_H
+
+#include <stdio.h>
+#include "math_extra.h"
+
+namespace MFOxdna {
+
+  inline double F1(double, double, double, double, double, double, double, double, double, double, double);
+  inline double DF1(double, double, double, double, double, double, double, double, double, double);
+  inline double F2(double, double, double, double, double, double, double, double, double, double);
+  inline double DF2(double, double, double, double, double, double, double, double, double);
+  inline double F3(double, double, double, double, double, double, double, double &);
+  inline double F4(double, double, double, double, double, double);
+  inline double DF4(double, double, double, double, double, double);
+  inline double F5(double, double, double, double, double);
+  inline double DF5(double, double, double, double, double);
+  inline double is_3pto5p(const double *, const double *);
+
+}
+
+/* ----------------------------------------------------------------------
+   f1 modulation factor
+   ------------------------------------------------------------------------- */
+inline double MFOxdna::F1(double r, double eps, double a, double cut_0,
+                          double cut_lc, double cut_hc, double cut_lo,
+                          double cut_hi, double b_lo, double b_hi,
+                          double shift)
+{
+
+  if (r > cut_hc) {
+    return 0.0;
+  }
+  else if (r > cut_hi) {
+    return eps * b_hi * (r-cut_hc) * (r-cut_hc);
+  }
+  else if (r > cut_lo) {
+    double tmp = 1 - exp(-(r-cut_0) * a);
+    return eps * tmp * tmp - shift;
+  }
+  else if (r > cut_lc) {
+    return eps * b_lo * (r-cut_lc) * (r-cut_lc);
+  }
+  else {
+    return 0.0;
+  }
+
+}
+
+/* ----------------------------------------------------------------------
+   derivative of f1 modulation factor
+   ------------------------------------------------------------------------- */
+inline double MFOxdna::DF1(double r, double eps, double a, double cut_0,
+                           double cut_lc, double cut_hc, double cut_lo,
+                           double cut_hi, double b_lo, double b_hi)
+{
+
+  if (r > cut_hc) {
+    return 0.0;
+  }
+  else if (r > cut_hi) {
+    return 2 * eps * b_hi * (1 - cut_hc / r);
+  }
+  else if (r > cut_lo) {
+    double tmp = exp(-(r-cut_0) * a);
+    return 2 * eps * (1 - tmp) * tmp * a / r;
+  }
+  else if (r > cut_lc) {
+    return 2 * eps * b_lo * (1 - cut_lc / r);
+  }
+  else {
+    return 0.0;
+  }
+
+}
+
+/* ----------------------------------------------------------------------
+   f2 modulation factor
+   ------------------------------------------------------------------------- */
+inline double MFOxdna::F2(double r, double k, double cut_0, double cut_lc,
+                          double cut_hc, double cut_lo, double cut_hi,
+                          double b_lo, double b_hi, double cut_c)
+{
+
+  if(r < cut_lc || r > cut_hc){
+    return 0;
+  }
+  else if(r < cut_lo){
+    return k * b_lo * (cut_lc - r)*(cut_lc-r);
+  }
+  else if(r < cut_hi){
+    return k * 0.5 * ((r - cut_0)*(r-cut_0) - (cut_0 - cut_c)*(cut_0 - cut_c));
+  }
+  else{
+    return k * b_hi * (cut_hc - r) * (cut_hc - r);
+  }
+
+}
+
+/* ----------------------------------------------------------------------
+   derivative of f2 modulation factor
+   ------------------------------------------------------------------------- */
+inline double MFOxdna::DF2(double r, double k, double cut_0, double cut_lc,
+                           double cut_hc, double cut_lo, double cut_hi,
+                           double b_lo, double b_hi)
+{
+  if(r < cut_lc || r > cut_hc){
+    return 0;
+  }
+  else if(r < cut_lo){
+    return 2*k * b_lo * (r - cut_lc);
+  }
+  else if(r < cut_hi){
+    return k * (r - cut_0);
+  }
+  else{
+    return 2*k * b_hi * (r - cut_hc);
+  }
+
+}
+
+/* ----------------------------------------------------------------------
+   f3 modulation factor, force and energy calculation
+   ------------------------------------------------------------------------- */
+inline double MFOxdna::F3(double rsq, double cutsq_ast, double cut_c,
+                          double lj1, double lj2, double eps, double b,
+                          double & fpair)
+{
+  double evdwl = 0.0;
+
+  if (rsq < cutsq_ast) {
+    double r2inv = 1.0/rsq;
+    double r6inv = r2inv*r2inv*r2inv;
+    fpair = r2inv*r6inv*(12*lj1*r6inv - 6*lj2);
+    evdwl = r6inv*(lj1*r6inv-lj2);
+  }
+  else {
+    double r = sqrt(rsq);
+    double rinv = 1.0/r;
+    fpair = 2*eps*b*(cut_c*rinv - 1);
+    evdwl = eps*b*(cut_c-r)*(cut_c-r);
+  }
+  return evdwl;
+}
+
+/* ----------------------------------------------------------------------
+   f4 modulation factor
+   ------------------------------------------------------------------------- */
+inline double MFOxdna::F4(double theta, double a, double theta_0,
+                          double dtheta_ast, double b, double dtheta_c)
+{
+  double dtheta = theta-theta_0;
+
+  if (fabs(dtheta) > dtheta_c) {
+    return 0.0;
+  }
+  else if (dtheta > dtheta_ast) {
+    return b * (dtheta-dtheta_c)*(dtheta-dtheta_c);
+  }
+  else if(dtheta > -dtheta_ast) {
+    return 1 - a * dtheta*dtheta;
+  }
+  else {
+    return b * (dtheta+dtheta_c)*(dtheta+dtheta_c);
+  }
+
+}
+
+/* ----------------------------------------------------------------------
+   derivative of f4 modulation factor
+
+   NOTE: We handle the sin(theta) factor from the partial derivative
+   of d(cos(theta))/dtheta externally. The reason for this is
+   because the sign of DF4 depends on the sign of theta in the
+   function call. It is also more efficient to store sin(theta).
+   ------------------------------------------------------------------------- */
+inline double MFOxdna::DF4(double theta, double a, double theta_0,
+                           double dtheta_ast, double b, double dtheta_c)
+{
+  double dtheta = theta-theta_0;
+
+  if (fabs(dtheta) > dtheta_c) {
+    return 0.0;
+  }
+  else if (dtheta > dtheta_ast) {
+    return 2*b* (dtheta-dtheta_c);
+  }
+  else if (dtheta > -dtheta_ast) {
+    return -2*a * dtheta;
+  }
+  else {
+    return 2*b* (dtheta+dtheta_c);
+  }
+
+}
+
+/* ----------------------------------------------------------------------
+   f5 modulation factor
+   ------------------------------------------------------------------------- */
+inline double MFOxdna::F5(double x, double a, double x_ast,
+                          double b, double x_c)
+{
+
+  if (x >= 0) {
+    return 1.0;
+  }
+  else if (x > x_ast) {
+    return 1 - a * x * x;
+  }
+  else if (x > x_c) {
+    return b * (x-x_c) * (x-x_c);
+  }
+  else {
+    return 0.0;
+  }
+}
+
+/* ----------------------------------------------------------------------
+   derivative of f5 modulation factor
+   ------------------------------------------------------------------------- */
+inline double MFOxdna::DF5(double x, double a, double x_ast,
+                           double b, double x_c)
+{
+  if(x >= 0) {
+    return 0.0;
+  }
+  else if (x > x_ast) {
+    return -2 * a * x;
+  }
+  else if(x > x_c) {
+    return 2 * b * (x-x_c);
+  }
+  else {
+    return 0.0;
+  }
+  return 0;
+}
+
+/* ----------------------------------------------------------------------
+   test for directionality by projecting base normal n onto delr,
+   returns 1 if nucleotide a to nucleotide b is 3' to 5', otherwise -1
+   ------------------------------------------------------------------------- */
+inline double MFOxdna::is_3pto5p(const double * delr, const double * n)
+{
+  return copysign(1.0,MathExtra::dot3(delr,n));
+}
+#endif
--- a/src/USER-CGDNA/pair_oxdna_coaxstk.cpp
+++ b/src/USER-CGDNA/pair_oxdna_coaxstk.cpp
--- a/src/USER-CGDNA/pair_oxdna_coaxstk.h
+++ b/src/USER-CGDNA/pair_oxdna_coaxstk.h
@ -0,0 +1,89 @@
+/* ----------------------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   http://lammps.sandia.gov, Sandia National Laboratories
+   Steve Plimpton, sjplimp@sandia.gov
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+/* ----------------------------------------------------------------------
+   Contributing author: Oliver Henrich (EPCC, University of Edinburgh)
+------------------------------------------------------------------------- */
+
+#ifdef PAIR_CLASS
+
+PairStyle(oxdna_coaxstk,PairOxdnaCoaxstk)
+
+#else
+
+#ifndef LMP_PAIR_OXDNA_COAXSTK_H
+#define LMP_PAIR_OXDNA_COAXSTK_H
+
+#include "pair.h"
+
+namespace LAMMPS_NS {
+
+class PairOxdnaCoaxstk : public Pair {
+ public:
+  PairOxdnaCoaxstk(class LAMMPS *);
+  virtual ~PairOxdnaCoaxstk();
+  virtual void compute(int, int);
+  void settings(int, char **);
+  void coeff(int, char **);
+  void init_style();
+  void init_list(int, class NeighList *);
+  double init_one(int, int);
+  void write_restart(FILE *);
+  void read_restart(FILE *);
+  void write_restart_settings(FILE *);
+  void read_restart_settings(FILE *);
+  void write_data(FILE *);
+  void write_data_all(FILE *);
+  void *extract(const char *, int &);
+
+ protected:
+  // coaxial stacking interaction
+  double **k_cxst, **cut_cxst_0, **cut_cxst_c, **cut_cxst_lo, **cut_cxst_hi;
+  double **cut_cxst_lc, **cut_cxst_hc, **b_cxst_lo, **b_cxst_hi;
+  double **cutsq_cxst_hc;
+
+  double **a_cxst1, **theta_cxst1_0, **dtheta_cxst1_ast;
+  double **b_cxst1, **dtheta_cxst1_c;
+
+  double **a_cxst4, **theta_cxst4_0, **dtheta_cxst4_ast;
+  double **b_cxst4, **dtheta_cxst4_c;
+
+  double **a_cxst5, **theta_cxst5_0, **dtheta_cxst5_ast;
+  double **b_cxst5, **dtheta_cxst5_c;
+
+  double **a_cxst6, **theta_cxst6_0, **dtheta_cxst6_ast;
+  double **b_cxst6, **dtheta_cxst6_c;
+
+  double **a_cxst3p, **cosphi_cxst3p_ast, **b_cxst3p, **cosphi_cxst3p_c;
+  double **a_cxst4p, **cosphi_cxst4p_ast, **b_cxst4p, **cosphi_cxst4p_c;
+
+  virtual void allocate();
+};
+
+}
+
+#endif
+#endif
+
+/* ERROR/WARNING messages:
+
+E: Illegal ... command
+
+Self-explanatory.  Check the input script syntax and compare to the
+documentation for the command.  You can use -echo screen as a
+command-line option when running LAMMPS to see the offending line.
+
+E: Incorrect args for pair coefficients
+
+Self-explanatory.  Check the input script or data file.
+
+*/
--- a/src/USER-CGDNA/pair_oxdna_excv.cpp
+++ b/src/USER-CGDNA/pair_oxdna_excv.cpp
@ -0,0 +1,840 @@
+/* ----------------------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   http://lammps.sandia.gov, Sandia National Laboratories
+   Steve Plimpton, sjplimp@sandia.gov
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+/* ----------------------------------------------------------------------
+   Contributing author: Oliver Henrich (EPCC, University of Edinburgh)
+------------------------------------------------------------------------- */
+
+#include <math.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include "pair_oxdna_excv.h"
+#include "mf_oxdna.h"
+#include "atom.h"
+#include "comm.h"
+#include "force.h"
+#include "neighbor.h"
+#include "neigh_list.h"
+#include "neigh_request.h"
+#include "update.h"
+#include "integrate.h"
+#include "math_const.h"
+#include "memory.h"
+#include "error.h"
+#include "atom_vec_ellipsoid.h"
+#include "math_extra.h"
+
+using namespace LAMMPS_NS;
+using namespace MathConst;
+using namespace MFOxdna;
+
+/* ---------------------------------------------------------------------- */
+
+PairOxdnaExcv::PairOxdnaExcv(LAMMPS *lmp) : Pair(lmp)
+{
+  single_enable = 0;
+  writedata = 1;
+}
+
+/* ---------------------------------------------------------------------- */
+
+PairOxdnaExcv::~PairOxdnaExcv()
+{
+  if (allocated) {
+
+    memory->destroy(setflag);
+    memory->destroy(cutsq);
+
+    memory->destroy(epsilon_ss);
+    memory->destroy(sigma_ss);
+    memory->destroy(cut_ss_ast);
+    memory->destroy(b_ss);
+    memory->destroy(cut_ss_c);
+    memory->destroy(lj1_ss);
+    memory->destroy(lj2_ss);
+    memory->destroy(cutsq_ss_ast);
+    memory->destroy(cutsq_ss_c);
+
+    memory->destroy(epsilon_sb);
+    memory->destroy(sigma_sb);
+    memory->destroy(cut_sb_ast);
+    memory->destroy(b_sb);
+    memory->destroy(cut_sb_c);
+    memory->destroy(lj1_sb);
+    memory->destroy(lj2_sb);
+    memory->destroy(cutsq_sb_ast);
+    memory->destroy(cutsq_sb_c);
+
+    memory->destroy(epsilon_bb);
+    memory->destroy(sigma_bb);
+    memory->destroy(cut_bb_ast);
+    memory->destroy(b_bb);
+    memory->destroy(cut_bb_c);
+    memory->destroy(lj1_bb);
+    memory->destroy(lj2_bb);
+    memory->destroy(cutsq_bb_ast);
+    memory->destroy(cutsq_bb_c);
+
+  }
+}
+
+/* ----------------------------------------------------------------------
+   compute function for oxDNA pair interactions
+   s=sugar-phosphate backbone site, b=base site, st=stacking site
+------------------------------------------------------------------------- */
+
+void PairOxdnaExcv::compute(int eflag, int vflag)
+{
+
+  double delf[3],delta[3],deltb[3]; // force, torque increment;
+  double evdwl,fpair,factor_lj;
+  double rtmp_s[3],rtmp_b[3];
+  double delr_ss[3],rsq_ss,delr_sb[3],rsq_sb;
+  double delr_bs[3],rsq_bs,delr_bb[3],rsq_bb;
+
+  // distances COM-backbone site, COM-base site
+  double d_cs=-0.24, d_cb=0.56;
+  // vectors COM-backbone site, COM-base site in lab frame
+  double ra_cs[3],ra_cb[3];
+  double rb_cs[3],rb_cb[3];
+
+  // quaternions and Cartesian unit vectors in lab frame
+  double *qa,ax[3],ay[3],az[3];
+  double *qb,bx[3],by[3],bz[3];
+  double *special_lj = force->special_lj;
+
+  double **x = atom->x;
+  double **f = atom->f;
+  double **torque = atom->torque;
+  int *type = atom->type;
+
+  int nlocal = atom->nlocal;
+  int newton_pair = force->newton_pair;
+  int *alist,*blist,*numneigh,**firstneigh;
+
+  AtomVecEllipsoid *avec = (AtomVecEllipsoid *) atom->style_match("ellipsoid");
+  AtomVecEllipsoid::Bonus *bonus = avec->bonus;
+
+  int a,b,ia,ib,anum,bnum,atype,btype;
+
+  evdwl = 0.0;
+  if (eflag || vflag) ev_setup(eflag,vflag);
+  else evflag = vflag_fdotr = 0;
+
+  anum = list->inum;
+  alist = list->ilist;
+  numneigh = list->numneigh;
+  firstneigh = list->firstneigh;
+
+  // loop over pair interaction neighbours of my atoms
+
+  for (ia = 0; ia < anum; ia++) {
+
+    a = alist[ia];
+    atype = type[a];
+
+    qa=bonus[a].quat;
+    MathExtra::q_to_exyz(qa,ax,ay,az);
+
+    // position of backbone site a
+    ra_cs[0] = d_cs*ax[0];
+    ra_cs[1] = d_cs*ax[1];
+    ra_cs[2] = d_cs*ax[2];
+    rtmp_s[0] = x[a][0] + ra_cs[0];
+    rtmp_s[1] = x[a][1] + ra_cs[1];
+    rtmp_s[2] = x[a][2] + ra_cs[2];
+
+    // position of base site a
+    ra_cb[0] = d_cb*ax[0];
+    ra_cb[1] = d_cb*ax[1];
+    ra_cb[2] = d_cb*ax[2];
+    rtmp_b[0] = x[a][0] + ra_cb[0];
+    rtmp_b[1] = x[a][1] + ra_cb[1];
+    rtmp_b[2] = x[a][2] + ra_cb[2];
+
+    blist = firstneigh[a];
+    bnum = numneigh[a];
+
+    for (ib = 0; ib < bnum; ib++) {
+
+      b = blist[ib];
+      factor_lj = special_lj[sbmask(b)]; // = 0 for nearest neighbours
+      b &= NEIGHMASK;
+
+      btype = type[b];
+
+      qb=bonus[b].quat;
+      MathExtra::q_to_exyz(qb,bx,by,bz);
+
+      rb_cs[0] = d_cs*bx[0];
+      rb_cs[1] = d_cs*bx[1];
+      rb_cs[2] = d_cs*bx[2];
+      rb_cb[0] = d_cb*bx[0];
+      rb_cb[1] = d_cb*bx[1];
+      rb_cb[2] = d_cb*bx[2];
+
+      // vector backbone site b to a
+      delr_ss[0] = rtmp_s[0] - (x[b][0] + rb_cs[0]);
+      delr_ss[1] = rtmp_s[1] - (x[b][1] + rb_cs[1]);
+      delr_ss[2] = rtmp_s[2] - (x[b][2] + rb_cs[2]);
+      rsq_ss = delr_ss[0]*delr_ss[0] + delr_ss[1]*delr_ss[1] + delr_ss[2]*delr_ss[2];
+
+      // vector base site b to backbone site a
+      delr_sb[0] =  rtmp_s[0] - (x[b][0] + rb_cb[0]);
+      delr_sb[1] =  rtmp_s[1] - (x[b][1] + rb_cb[1]);
+      delr_sb[2] =  rtmp_s[2] - (x[b][2] + rb_cb[2]);
+      rsq_sb = delr_sb[0]*delr_sb[0] + delr_sb[1]*delr_sb[1] + delr_sb[2]*delr_sb[2];
+
+      // vector backbone site b to base site a
+      delr_bs[0] = rtmp_b[0] - (x[b][0] + rb_cs[0]);
+      delr_bs[1] = rtmp_b[1] - (x[b][1] + rb_cs[1]);
+      delr_bs[2] = rtmp_b[2] - (x[b][2] + rb_cs[2]);
+      rsq_bs = delr_bs[0]*delr_bs[0] + delr_bs[1]*delr_bs[1] + delr_bs[2]*delr_bs[2];
+
+      // vector base site b to a
+      delr_bb[0] = rtmp_b[0] - (x[b][0] + rb_cb[0]);
+      delr_bb[1] = rtmp_b[1] - (x[b][1] + rb_cb[1]);
+      delr_bb[2] = rtmp_b[2] - (x[b][2] + rb_cb[2]);
+      rsq_bb = delr_bb[0]*delr_bb[0] + delr_bb[1]*delr_bb[1] + delr_bb[2]*delr_bb[2];
+
+      // excluded volume interaction
+
+      // backbone-backbone
+      if (rsq_ss < cutsq_ss_c[atype][btype]) {
+
+        evdwl = F3(rsq_ss,cutsq_ss_ast[atype][btype],cut_ss_c[atype][btype],lj1_ss[atype][btype],
+                        lj2_ss[atype][btype],epsilon_ss[atype][btype],b_ss[atype][btype],fpair);
+
+        // knock out nearest-neighbour interaction between ss
+        fpair *= factor_lj;
+        evdwl *= factor_lj;
+
+        // increment energy and virial
+        if (evflag) ev_tally(a,b,nlocal,newton_pair,
+                evdwl,0.0,fpair,delr_ss[0],delr_ss[1],delr_ss[2]);
+
+        delf[0] = delr_ss[0]*fpair;
+        delf[1] = delr_ss[1]*fpair;
+        delf[2] = delr_ss[2]*fpair;
+
+        f[a][0] += delf[0];
+        f[a][1] += delf[1];
+        f[a][2] += delf[2];
+
+        MathExtra::cross3(ra_cs,delf,delta);
+
+        torque[a][0] += delta[0];
+        torque[a][1] += delta[1];
+        torque[a][2] += delta[2];
+
+        if (newton_pair || b < nlocal) {
+
+          f[b][0] -= delf[0];
+          f[b][1] -= delf[1];
+          f[b][2] -= delf[2];
+
+          MathExtra::cross3(rb_cs,delf,deltb);
+
+          torque[b][0] -= deltb[0];
+          torque[b][1] -= deltb[1];
+          torque[b][2] -= deltb[2];
+
+        }
+
+      }
+
+
+      // backbone-base
+      if (rsq_sb < cutsq_sb_c[atype][btype]) {
+
+        evdwl = F3(rsq_sb,cutsq_sb_ast[atype][btype],cut_sb_c[atype][btype],lj1_sb[atype][btype],
+                        lj2_sb[atype][btype],epsilon_sb[atype][btype],b_sb[atype][btype],fpair);
+
+        // increment energy and virial
+        if (evflag) ev_tally(a,b,nlocal,newton_pair,
+                evdwl,0.0,fpair,delr_sb[0],delr_sb[1],delr_sb[2]);
+
+        delf[0] = delr_sb[0]*fpair;
+        delf[1] = delr_sb[1]*fpair;
+        delf[2] = delr_sb[2]*fpair;
+
+        f[a][0] += delf[0];
+        f[a][1] += delf[1];
+        f[a][2] += delf[2];
+
+        MathExtra::cross3(ra_cs,delf,delta);
+
+        torque[a][0] += delta[0];
+        torque[a][1] += delta[1];
+        torque[a][2] += delta[2];
+
+        if (newton_pair || b < nlocal) {
+
+          f[b][0] -= delf[0];
+          f[b][1] -= delf[1];
+          f[b][2] -= delf[2];
+
+          MathExtra::cross3(rb_cb,delf,deltb);
+
+          torque[b][0] -= deltb[0];
+          torque[b][1] -= deltb[1];
+          torque[b][2] -= deltb[2];
+
+        }
+
+      }
+
+      // base-backbone
+      if (rsq_bs < cutsq_sb_c[atype][btype]) {
+
+        evdwl = F3(rsq_bs,cutsq_sb_ast[atype][btype],cut_sb_c[atype][btype],lj1_sb[atype][btype],
+                        lj2_sb[atype][btype],epsilon_sb[atype][btype],b_sb[atype][btype],fpair);
+
+        // increment energy and virial
+        if (evflag) ev_tally(a,b,nlocal,newton_pair,
+                evdwl,0.0,fpair,delr_bs[0],delr_bs[1],delr_bs[2]);
+
+        delf[0] = delr_bs[0]*fpair;
+        delf[1] = delr_bs[1]*fpair;
+        delf[2] = delr_bs[2]*fpair;
+
+        f[a][0] += delf[0];
+        f[a][1] += delf[1];
+        f[a][2] += delf[2];
+
+        MathExtra::cross3(ra_cb,delf,delta);
+
+        torque[a][0] += delta[0];
+        torque[a][1] += delta[1];
+        torque[a][2] += delta[2];
+
+        if (newton_pair || b < nlocal) {
+
+          f[b][0] -= delf[0];
+          f[b][1] -= delf[1];
+          f[b][2] -= delf[2];
+
+          MathExtra::cross3(rb_cs,delf,deltb);
+
+          torque[b][0] -= deltb[0];
+          torque[b][1] -= deltb[1];
+          torque[b][2] -= deltb[2];
+
+        }
+
+      }
+
+      // base-base
+      if (rsq_bb < cutsq_bb_c[atype][btype]) {
+
+        evdwl = F3(rsq_bb,cutsq_bb_ast[atype][btype],cut_bb_c[atype][btype],lj1_bb[atype][btype],
+                        lj2_bb[atype][btype],epsilon_bb[atype][btype],b_bb[atype][btype],fpair);
+
+        // increment energy and virial
+        if (evflag) ev_tally(a,b,nlocal,newton_pair,
+                evdwl,0.0,fpair,delr_bb[0],delr_bb[1],delr_bb[2]);
+
+        delf[0] = delr_bb[0]*fpair;
+        delf[1] = delr_bb[1]*fpair;
+        delf[2] = delr_bb[2]*fpair;
+
+        f[a][0] += delf[0];
+        f[a][1] += delf[1];
+        f[a][2] += delf[2];
+
+        MathExtra::cross3(ra_cb,delf,delta);
+
+        torque[a][0] += delta[0];
+        torque[a][1] += delta[1];
+        torque[a][2] += delta[2];
+
+        if (newton_pair || b < nlocal) {
+
+          f[b][0] -= delf[0];
+          f[b][1] -= delf[1];
+          f[b][2] -= delf[2];
+
+          MathExtra::cross3(rb_cb,delf,deltb);
+
+          torque[b][0] -= deltb[0];
+          torque[b][1] -= deltb[1];
+          torque[b][2] -= deltb[2];
+
+        }
+
+      }
+      // end excluded volume interaction
+
+    }
+  }
+
+  if (vflag_fdotr) virial_fdotr_compute();
+
+}
+
+/* ----------------------------------------------------------------------
+   allocate all arrays
+------------------------------------------------------------------------- */
+
+void PairOxdnaExcv::allocate()
+{
+  allocated = 1;
+  int n = atom->ntypes;
+
+  memory->create(setflag,n+1,n+1,"pair:setflag");
+  for (int i = 1; i <= n; i++)
+    for (int j = i; j <= n; j++)
+      setflag[i][j] = 0;
+
+  memory->create(cutsq,n+1,n+1,"pair:cutsq");
+
+  memory->create(epsilon_ss,n+1,n+1,"pair:epsilon_ss");
+  memory->create(sigma_ss,n+1,n+1,"pair:sigma_ss");
+  memory->create(cut_ss_ast,n+1,n+1,"pair:cut_ss_ast");
+  memory->create(b_ss,n+1,n+1,"pair:b_ss");
+  memory->create(cut_ss_c,n+1,n+1,"pair:cut_ss_c");
+  memory->create(lj1_ss,n+1,n+1,"pair:lj1_ss");
+  memory->create(lj2_ss,n+1,n+1,"pair:lj2_ss");
+  memory->create(cutsq_ss_ast,n+1,n+1,"pair:cutsq_ss_ast");
+  memory->create(cutsq_ss_c,n+1,n+1,"pair:cutsq_ss_c");
+
+  memory->create(epsilon_sb,n+1,n+1,"pair:epsilon_sb");
+  memory->create(sigma_sb,n+1,n+1,"pair:sigma_sb");
+  memory->create(cut_sb_ast,n+1,n+1,"pair:cut_sb_ast");
+  memory->create(b_sb,n+1,n+1,"pair:b_sb");
+  memory->create(cut_sb_c,n+1,n+1,"pair:cut_sb_c");
+  memory->create(lj1_sb,n+1,n+1,"pair:lj1_sb");
+  memory->create(lj2_sb,n+1,n+1,"pair:lj2_sb");
+  memory->create(cutsq_sb_ast,n+1,n+1,"pair:cutsq_sb_ast");
+  memory->create(cutsq_sb_c,n+1,n+1,"pair:cutsq_sb_c");
+
+  memory->create(epsilon_bb,n+1,n+1,"pair:epsilon_bb");
+  memory->create(sigma_bb,n+1,n+1,"pair:sigma_bb");
+  memory->create(cut_bb_ast,n+1,n+1,"pair:cut_bb_ast");
+  memory->create(b_bb,n+1,n+1,"pair:b_bb");
+  memory->create(cut_bb_c,n+1,n+1,"pair:cut_bb_c");
+  memory->create(lj1_bb,n+1,n+1,"pair:lj1_bb");
+  memory->create(lj2_bb,n+1,n+1,"pair:lj2_bb");
+  memory->create(cutsq_bb_ast,n+1,n+1,"pair:cutsq_bb_ast");
+  memory->create(cutsq_bb_c,n+1,n+1,"pair:cutsq_bb_c");
+
+}
+
+/* ----------------------------------------------------------------------
+   global settings
+------------------------------------------------------------------------- */
+
+void PairOxdnaExcv::settings(int narg, char **arg)
+{
+  if (narg != 0) error->all(FLERR,"Illegal pair_style command");
+
+}
+
+/* ----------------------------------------------------------------------
+   set coeffs for one or more type pairs
+------------------------------------------------------------------------- */
+
+void PairOxdnaExcv::coeff(int narg, char **arg)
+{
+  int count;
+
+  if (narg != 11) error->all(FLERR,"Incorrect args for pair coefficients in oxdna_excv");
+  if (!allocated) allocate();
+
+  int ilo,ihi,jlo,jhi;
+  force->bounds(FLERR,arg[0],atom->ntypes,ilo,ihi);
+  force->bounds(FLERR,arg[1],atom->ntypes,jlo,jhi);
+
+  count = 0;
+
+  double epsilon_ss_one, sigma_ss_one;
+  double cut_ss_ast_one, cut_ss_c_one, b_ss_one;
+
+  double epsilon_sb_one, sigma_sb_one;
+  double cut_sb_ast_one, cut_sb_c_one, b_sb_one;
+
+  double epsilon_bb_one, sigma_bb_one;
+  double cut_bb_ast_one, cut_bb_c_one, b_bb_one;
+
+  // Excluded volume interaction
+  // LJ parameters
+  epsilon_ss_one = force->numeric(FLERR,arg[2]);
+  sigma_ss_one = force->numeric(FLERR,arg[3]);
+  cut_ss_ast_one = force->numeric(FLERR,arg[4]);
+
+  // smoothing - determined through continuity and differentiability
+  b_ss_one = 4.0/sigma_ss_one
+      *(6.0*pow(sigma_ss_one/cut_ss_ast_one,7)-12.0*pow(sigma_ss_one/cut_ss_ast_one,13))
+      *4.0/sigma_ss_one*(6.0*pow(sigma_ss_one/cut_ss_ast_one,7)-12.0*pow(sigma_ss_one/cut_ss_ast_one,13))
+      /4.0/(4.0*(pow(sigma_ss_one/cut_ss_ast_one,12)-pow(sigma_ss_one/cut_ss_ast_one,6)));
+
+  cut_ss_c_one = cut_ss_ast_one
+      - 2.0*4.0*(pow(sigma_ss_one/cut_ss_ast_one,12)-pow(sigma_ss_one/cut_ss_ast_one,6))
+      /(4.0/sigma_ss_one*(6.0*pow(sigma_ss_one/cut_ss_ast_one,7)-12.0*pow(sigma_ss_one/cut_ss_ast_one,13)));
+
+  for (int i = ilo; i <= ihi; i++) {
+    for (int j = MAX(jlo,i); j <= jhi; j++) {
+      epsilon_ss[i][j] = epsilon_ss_one;
+      sigma_ss[i][j] = sigma_ss_one;
+      cut_ss_ast[i][j] = cut_ss_ast_one;
+      b_ss[i][j] = b_ss_one;
+      cut_ss_c[i][j] = cut_ss_c_one;
+      setflag[i][j] = 1;
+      count++;
+    }
+  }
+
+  if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients in oxdna_excv");
+
+  count = 0;
+
+  // LJ parameters
+  epsilon_sb_one = force->numeric(FLERR,arg[5]);
+  sigma_sb_one = force->numeric(FLERR,arg[6]);
+  cut_sb_ast_one = force->numeric(FLERR,arg[7]);
+
+  // smoothing - determined through continuity and differentiability
+  b_sb_one = 4.0/sigma_sb_one
+      *(6.0*pow(sigma_sb_one/cut_sb_ast_one,7)-12.0*pow(sigma_sb_one/cut_sb_ast_one,13))
+      *4.0/sigma_sb_one*(6.0*pow(sigma_sb_one/cut_sb_ast_one,7)-12.0*pow(sigma_sb_one/cut_sb_ast_one,13))
+      /4.0/(4.0*(pow(sigma_sb_one/cut_sb_ast_one,12)-pow(sigma_sb_one/cut_sb_ast_one,6)));
+
+  cut_sb_c_one = cut_sb_ast_one
+      - 2.0*4.0*(pow(sigma_sb_one/cut_sb_ast_one,12)-pow(sigma_sb_one/cut_sb_ast_one,6))
+      /(4.0/sigma_sb_one*(6.0*pow(sigma_sb_one/cut_sb_ast_one,7)-12.0*pow(sigma_sb_one/cut_sb_ast_one,13)));
+
+  for (int i = ilo; i <= ihi; i++) {
+    for (int j = MAX(jlo,i); j <= jhi; j++) {
+      epsilon_sb[i][j] = epsilon_sb_one;
+      sigma_sb[i][j] = sigma_sb_one;
+      cut_sb_ast[i][j] = cut_sb_ast_one;
+      b_sb[i][j] = b_sb_one;
+      cut_sb_c[i][j] = cut_sb_c_one;
+      setflag[i][j] = 1;
+      count++;
+    }
+  }
+
+  if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients in oxdna_excv");
+
+  count = 0;
+
+  // LJ parameters
+  epsilon_bb_one = force->numeric(FLERR,arg[8]);
+  sigma_bb_one = force->numeric(FLERR,arg[9]);
+  cut_bb_ast_one = force->numeric(FLERR,arg[10]);
+
+  // smoothing - determined through continuity and differentiability
+  b_bb_one = 4.0/sigma_bb_one
+      *(6.0*pow(sigma_bb_one/cut_bb_ast_one,7)-12.0*pow(sigma_bb_one/cut_bb_ast_one,13))
+      *4.0/sigma_bb_one*(6.0*pow(sigma_bb_one/cut_bb_ast_one,7)-12.0*pow(sigma_bb_one/cut_bb_ast_one,13))
+      /4.0/(4.0*(pow(sigma_bb_one/cut_bb_ast_one,12)-pow(sigma_bb_one/cut_bb_ast_one,6)));
+
+  cut_bb_c_one = cut_bb_ast_one
+      - 2.0*4.0*(pow(sigma_bb_one/cut_bb_ast_one,12)-pow(sigma_bb_one/cut_bb_ast_one,6))
+      /(4.0/sigma_bb_one*(6.0*pow(sigma_bb_one/cut_bb_ast_one,7)-12.0*pow(sigma_bb_one/cut_bb_ast_one,13)));
+
+  for (int i = ilo; i <= ihi; i++) {
+    for (int j = MAX(jlo,i); j <= jhi; j++) {
+      epsilon_bb[i][j] = epsilon_bb_one;
+      sigma_bb[i][j] = sigma_bb_one;
+      cut_bb_ast[i][j] = cut_bb_ast_one;
+      b_bb[i][j] = b_bb_one;
+      cut_bb_c[i][j] = cut_bb_c_one;
+      setflag[i][j] = 1;
+      count++;
+    }
+  }
+
+  if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients in oxdna_excv");
+
+}
+
+/* ----------------------------------------------------------------------
+   init specific to this pair style
+------------------------------------------------------------------------- */
+
+void PairOxdnaExcv::init_style()
+{
+  int irequest;
+
+  // request regular neighbor lists
+
+  irequest = neighbor->request(this,instance_me);
+
+}
+
+/* ----------------------------------------------------------------------
+   neighbor callback to inform pair style of neighbor list to use regular
+------------------------------------------------------------------------- */
+
+void PairOxdnaExcv::init_list(int id, NeighList *ptr)
+{
+  if (id == 0) list = ptr;
+  if (id  > 0) error->all(FLERR,"Respa not supported");
+
+}
+
+
+/* ----------------------------------------------------------------------
+   init for one type pair i,j and corresponding j,i
+------------------------------------------------------------------------- */
+
+double PairOxdnaExcv::init_one(int i, int j)
+{
+
+  if (setflag[i][j] == 0) {
+    error->all(FLERR,"Coefficient mixing not defined in oxDNA");
+  }
+  if (offset_flag) {
+    error->all(FLERR,"Offset not supported in oxDNA");
+  }
+
+  epsilon_ss[j][i] = epsilon_ss[i][j];
+  sigma_ss[j][i] = sigma_ss[i][j];
+  cut_ss_ast[j][i] = cut_ss_ast[i][j];
+  cut_ss_c[j][i] = cut_ss_c[i][j];
+  b_ss[j][i] = b_ss[i][j];
+
+  epsilon_sb[j][i] = epsilon_sb[i][j];
+  sigma_sb[j][i] = sigma_sb[i][j];
+  cut_sb_ast[j][i] = cut_sb_ast[i][j];
+  cut_sb_c[j][i] = cut_sb_c[i][j];
+  b_sb[j][i] = b_sb[i][j];
+
+  epsilon_bb[j][i] = epsilon_bb[i][j];
+  sigma_bb[j][i] = sigma_bb[i][j];
+  cut_bb_ast[j][i] = cut_bb_ast[i][j];
+  cut_bb_c[j][i] = cut_bb_c[i][j];
+  b_bb[j][i] = b_bb[i][j];
+
+  // excluded volume auxiliary parameters
+
+  lj1_ss[i][j] = 4.0 * epsilon_ss[i][j] * pow(sigma_ss[i][j],12.0);
+  lj2_ss[i][j] = 4.0 * epsilon_ss[i][j] * pow(sigma_ss[i][j],6.0);
+
+  lj1_sb[i][j] = 4.0 * epsilon_sb[i][j] * pow(sigma_sb[i][j],12.0);
+  lj2_sb[i][j] = 4.0 * epsilon_sb[i][j] * pow(sigma_sb[i][j],6.0);
+
+  lj1_bb[i][j] = 4.0 * epsilon_bb[i][j] * pow(sigma_bb[i][j],12.0);
+  lj2_bb[i][j] = 4.0 * epsilon_bb[i][j] * pow(sigma_bb[i][j],6.0);
+
+  lj1_ss[j][i] = lj1_ss[i][j];
+  lj2_ss[j][i] = lj2_ss[i][j];
+
+  lj1_sb[j][i] = lj1_sb[i][j];
+  lj2_sb[j][i] = lj2_sb[i][j];
+
+  lj1_bb[j][i] = lj1_bb[i][j];
+  lj2_bb[j][i] = lj2_bb[i][j];
+
+  cutsq_ss_ast[i][j] = cut_ss_ast[i][j]*cut_ss_ast[i][j];
+  cutsq_ss_c[i][j]  = cut_ss_c[i][j]*cut_ss_c[i][j];
+
+  cutsq_sb_ast[i][j] = cut_sb_ast[i][j]*cut_sb_ast[i][j];
+  cutsq_sb_c[i][j]  = cut_sb_c[i][j]*cut_sb_c[i][j];
+
+  cutsq_bb_ast[i][j] = cut_bb_ast[i][j]*cut_bb_ast[i][j];
+  cutsq_bb_c[i][j]  = cut_bb_c[i][j]*cut_bb_c[i][j];
+
+  cutsq_ss_ast[j][i] = cutsq_ss_ast[i][j];
+  cutsq_ss_c[j][i]  = cutsq_ss_c[i][j];
+
+  cutsq_sb_ast[j][i] = cutsq_sb_ast[i][j];
+  cutsq_sb_c[j][i]  = cutsq_sb_c[i][j];
+
+  cutsq_bb_ast[j][i] = cutsq_bb_ast[i][j];
+  cutsq_bb_c[j][i]  = cutsq_bb_c[i][j];
+
+  // set the master list distance cutoff
+  return cut_ss_ast[i][j];
+
+}
+
+/* ----------------------------------------------------------------------
+   proc 0 writes to restart file
+------------------------------------------------------------------------- */
+
+void PairOxdnaExcv::write_restart(FILE *fp)
+{
+  write_restart_settings(fp);
+
+  int i,j;
+  for (i = 1; i <= atom->ntypes; i++)
+    for (j = i; j <= atom->ntypes; j++) {
+      fwrite(&setflag[i][j],sizeof(int),1,fp);
+      if (setflag[i][j]) {
+
+        fwrite(&epsilon_ss[i][j],sizeof(double),1,fp);
+        fwrite(&sigma_ss[i][j],sizeof(double),1,fp);
+        fwrite(&cut_ss_ast[i][j],sizeof(double),1,fp);
+        fwrite(&b_ss[i][j],sizeof(double),1,fp);
+        fwrite(&cut_ss_c[i][j],sizeof(double),1,fp);
+        fwrite(&epsilon_sb[i][j],sizeof(double),1,fp);
+        fwrite(&sigma_sb[i][j],sizeof(double),1,fp);
+        fwrite(&cut_sb_ast[i][j],sizeof(double),1,fp);
+        fwrite(&b_sb[i][j],sizeof(double),1,fp);
+        fwrite(&cut_sb_c[i][j],sizeof(double),1,fp);
+        fwrite(&epsilon_bb[i][j],sizeof(double),1,fp);
+        fwrite(&sigma_bb[i][j],sizeof(double),1,fp);
+        fwrite(&cut_bb_ast[i][j],sizeof(double),1,fp);
+        fwrite(&b_bb[i][j],sizeof(double),1,fp);
+        fwrite(&cut_bb_c[i][j],sizeof(double),1,fp);
+
+    }
+  }
+}
+
+/* ----------------------------------------------------------------------
+   proc 0 reads from restart file, bcasts
+------------------------------------------------------------------------- */
+
+void PairOxdnaExcv::read_restart(FILE *fp)
+{
+  read_restart_settings(fp);
+  allocate();
+
+  int i,j;
+  int me = comm->me;
+  for (i = 1; i <= atom->ntypes; i++)
+    for (j = i; j <= atom->ntypes; j++) {
+      if (me == 0) fread(&setflag[i][j],sizeof(int),1,fp);
+      MPI_Bcast(&setflag[i][j],1,MPI_INT,0,world);
+      if (setflag[i][j]) {
+        if (me == 0) {
+
+          fread(&epsilon_ss[i][j],sizeof(double),1,fp);
+          fread(&sigma_ss[i][j],sizeof(double),1,fp);
+          fread(&cut_ss_ast[i][j],sizeof(double),1,fp);
+          fread(&b_ss[i][j],sizeof(double),1,fp);
+          fread(&cut_ss_c[i][j],sizeof(double),1,fp);
+          fread(&epsilon_sb[i][j],sizeof(double),1,fp);
+          fread(&sigma_sb[i][j],sizeof(double),1,fp);
+          fread(&cut_sb_ast[i][j],sizeof(double),1,fp);
+          fread(&b_sb[i][j],sizeof(double),1,fp);
+          fread(&cut_sb_c[i][j],sizeof(double),1,fp);
+          fread(&epsilon_bb[i][j],sizeof(double),1,fp);
+          fread(&sigma_bb[i][j],sizeof(double),1,fp);
+          fread(&cut_bb_ast[i][j],sizeof(double),1,fp);
+          fread(&b_bb[i][j],sizeof(double),1,fp);
+          fread(&cut_bb_c[i][j],sizeof(double),1,fp);
+
+         }
+
+        MPI_Bcast(&epsilon_ss[i][j],1,MPI_DOUBLE,0,world);
+        MPI_Bcast(&sigma_ss[i][j],1,MPI_DOUBLE,0,world);
+        MPI_Bcast(&cut_ss_ast[i][j],1,MPI_DOUBLE,0,world);
+        MPI_Bcast(&b_ss[i][j],1,MPI_DOUBLE,0,world);
+        MPI_Bcast(&cut_ss_c[i][j],1,MPI_DOUBLE,0,world);
+        MPI_Bcast(&epsilon_sb[i][j],1,MPI_DOUBLE,0,world);
+        MPI_Bcast(&sigma_sb[i][j],1,MPI_DOUBLE,0,world);
+        MPI_Bcast(&cut_sb_ast[i][j],1,MPI_DOUBLE,0,world);
+        MPI_Bcast(&b_sb[i][j],1,MPI_DOUBLE,0,world);
+        MPI_Bcast(&cut_sb_c[i][j],1,MPI_DOUBLE,0,world);
+        MPI_Bcast(&epsilon_bb[i][j],1,MPI_DOUBLE,0,world);
+        MPI_Bcast(&sigma_bb[i][j],1,MPI_DOUBLE,0,world);
+        MPI_Bcast(&cut_bb_ast[i][j],1,MPI_DOUBLE,0,world);
+        MPI_Bcast(&b_bb[i][j],1,MPI_DOUBLE,0,world);
+        MPI_Bcast(&cut_bb_c[i][j],1,MPI_DOUBLE,0,world);
+
+      }
+    }
+}
+
+/* ----------------------------------------------------------------------
+   proc 0 writes to restart file
+------------------------------------------------------------------------- */
+
+void PairOxdnaExcv::write_restart_settings(FILE *fp)
+{
+  fwrite(&offset_flag,sizeof(int),1,fp);
+  fwrite(&mix_flag,sizeof(int),1,fp);
+  fwrite(&tail_flag,sizeof(int),1,fp);
+}
+
+/* ----------------------------------------------------------------------
+   proc 0 reads from restart file, bcasts
+------------------------------------------------------------------------- */
+
+void PairOxdnaExcv::read_restart_settings(FILE *fp)
+{
+  int me = comm->me;
+  if (me == 0) {
+    fread(&offset_flag,sizeof(int),1,fp);
+    fread(&mix_flag,sizeof(int),1,fp);
+    fread(&tail_flag,sizeof(int),1,fp);
+  }
+  MPI_Bcast(&offset_flag,1,MPI_INT,0,world);
+  MPI_Bcast(&mix_flag,1,MPI_INT,0,world);
+  MPI_Bcast(&tail_flag,1,MPI_INT,0,world);
+}
+
+/* ----------------------------------------------------------------------
+   proc 0 writes to data file
+------------------------------------------------------------------------- */
+
+void PairOxdnaExcv::write_data(FILE *fp)
+{
+  for (int i = 1; i <= atom->ntypes; i++)
+    fprintf(fp,"%d\
+         %g %g %g %g %g\
+         %g %g %g %g %g\
+         %g %g %g %g %g\
+         \n",i,
+        epsilon_ss[i][i],sigma_ss[i][i],cut_ss_ast[i][i],b_ss[i][i],cut_ss_c[i][i],
+        epsilon_sb[i][i],sigma_sb[i][i],cut_sb_ast[i][i],b_sb[i][i],cut_sb_c[i][i],
+        epsilon_bb[i][i],sigma_bb[i][i],cut_bb_ast[i][i],b_bb[i][i],cut_bb_c[i][i]);
+}
+
+/* ----------------------------------------------------------------------
+   proc 0 writes all pairs to data file
+------------------------------------------------------------------------- */
+
+void PairOxdnaExcv::write_data_all(FILE *fp)
+{
+  for (int i = 1; i <= atom->ntypes; i++)
+    for (int j = i; j <= atom->ntypes; j++)
+      fprintf(fp,"%d %d\
+         %g %g %g %g %g\
+         %g %g %g %g %g\
+         %g %g %g %g %g\
+         \n",i,j,
+        epsilon_ss[i][j],sigma_ss[i][j],cut_ss_ast[i][j],b_ss[i][j],cut_ss_c[i][j],
+        epsilon_sb[i][j],sigma_sb[i][j],cut_sb_ast[i][j],b_sb[i][j],cut_sb_c[i][j],
+        epsilon_bb[i][j],sigma_bb[i][j],cut_bb_ast[i][j],b_bb[i][j],cut_bb_c[i][j]);
+}
+
+/* ---------------------------------------------------------------------- */
+
+void *PairOxdnaExcv::extract(const char *str, int &dim)
+{
+  dim = 2;
+
+  if (strcmp(str,"epsilon_ss") == 0) return (void *) epsilon_ss;
+  if (strcmp(str,"sigma_ss") == 0) return (void *) sigma_ss;
+  if (strcmp(str,"cut_ss_ast") == 0) return (void *) cut_ss_ast;
+  if (strcmp(str,"b_ss") == 0) return (void *) b_ss;
+  if (strcmp(str,"cut_ss_c") == 0) return (void *) cut_ss_c;
+  if (strcmp(str,"epsilon_sb") == 0) return (void *) epsilon_sb;
+  if (strcmp(str,"sigma_sb") == 0) return (void *) sigma_sb;
+  if (strcmp(str,"cut_sb_ast") == 0) return (void *) cut_sb_ast;
+  if (strcmp(str,"b_sb") == 0) return (void *) b_sb;
+  if (strcmp(str,"cut_sb_c") == 0) return (void *) cut_sb_c;
+  if (strcmp(str,"epsilon_bb") == 0) return (void *) epsilon_bb;
+  if (strcmp(str,"sigma_bb") == 0) return (void *) sigma_bb;
+  if (strcmp(str,"cut_bb_ast") == 0) return (void *) cut_bb_ast;
+  if (strcmp(str,"b_bb") == 0) return (void *) b_bb;
+  if (strcmp(str,"cut_bb_c") == 0) return (void *) cut_bb_c;
+
+  return NULL;
+}
--- a/src/USER-CGDNA/pair_oxdna_excv.h
+++ b/src/USER-CGDNA/pair_oxdna_excv.h
@ -0,0 +1,80 @@
+/* ----------------------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   http://lammps.sandia.gov, Sandia National Laboratories
+   Steve Plimpton, sjplimp@sandia.gov
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+/* ----------------------------------------------------------------------
+   Contributing author: Oliver Henrich (EPCC, University of Edinburgh)
+------------------------------------------------------------------------- */
+
+#ifdef PAIR_CLASS
+
+PairStyle(oxdna_excv,PairOxdnaExcv)
+
+#else
+
+#ifndef LMP_PAIR_OXDNA_EXCV_H
+#define LMP_PAIR_OXDNA_EXCV_H
+
+#include "pair.h"
+
+namespace LAMMPS_NS {
+
+class PairOxdnaExcv : public Pair {
+ public:
+  PairOxdnaExcv(class LAMMPS *);
+  virtual ~PairOxdnaExcv();
+  virtual void compute(int, int);
+  void settings(int, char **);
+  void coeff(int, char **);
+  void init_style();
+  void init_list(int, class NeighList *);
+  double init_one(int, int);
+  void write_restart(FILE *);
+  void read_restart(FILE *);
+  void write_restart_settings(FILE *);
+  void read_restart_settings(FILE *);
+  void write_data(FILE *);
+  void write_data_all(FILE *);
+  void *extract(const char *, int &);
+
+ protected:
+  // s=sugar-phosphate backbone site, b=base site, st=stacking site
+
+  // excluded volume interaction
+  double **epsilon_ss, **sigma_ss, **cut_ss_ast, **cutsq_ss_ast;
+  double **lj1_ss, **lj2_ss, **b_ss, **cut_ss_c, **cutsq_ss_c;
+  double **epsilon_sb, **sigma_sb, **cut_sb_ast, **cutsq_sb_ast;
+  double **lj1_sb, **lj2_sb, **b_sb, **cut_sb_c, **cutsq_sb_c;
+  double **epsilon_bb, **sigma_bb, **cut_bb_ast, **cutsq_bb_ast;
+  double **lj1_bb, **lj2_bb, **b_bb, **cut_bb_c, **cutsq_bb_c;
+
+  virtual void allocate();
+
+};
+
+}
+
+#endif
+#endif
+
+/* ERROR/WARNING messages:
+
+E: Illegal ... command
+
+Self-explanatory.  Check the input script syntax and compare to the
+documentation for the command.  You can use -echo screen as a
+command-line option when running LAMMPS to see the offending line.
+
+E: Incorrect args for pair coefficients
+
+Self-explanatory.  Check the input script or data file.
+
+*/
--- a/src/USER-CGDNA/pair_oxdna_hbond.cpp
+++ b/src/USER-CGDNA/pair_oxdna_hbond.cpp
--- a/src/USER-CGDNA/pair_oxdna_hbond.h
+++ b/src/USER-CGDNA/pair_oxdna_hbond.h
@ -0,0 +1,92 @@
+/* ----------------------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   http://lammps.sandia.gov, Sandia National Laboratories
+   Steve Plimpton, sjplimp@sandia.gov
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+/* ----------------------------------------------------------------------
+   Contributing author: Oliver Henrich (EPCC, University of Edinburgh)
+------------------------------------------------------------------------- */
+
+#ifdef PAIR_CLASS
+
+PairStyle(oxdna_hbond,PairOxdnaHbond)
+
+#else
+
+#ifndef LMP_PAIR_OXDNA_HBOND_H
+#define LMP_PAIR_OXDNA_HBOND_H
+
+#include "pair.h"
+
+namespace LAMMPS_NS {
+
+class PairOxdnaHbond : public Pair {
+ public:
+  PairOxdnaHbond(class LAMMPS *);
+  virtual ~PairOxdnaHbond();
+  virtual void compute(int, int);
+  void settings(int, char **);
+  void coeff(int, char **);
+  void init_style();
+  void init_list(int, class NeighList *);
+  double init_one(int, int);
+  void write_restart(FILE *);
+  void read_restart(FILE *);
+  void write_restart_settings(FILE *);
+  void read_restart_settings(FILE *);
+  void write_data(FILE *);
+  void write_data_all(FILE *);
+  void *extract(const char *, int &);
+
+ protected:
+  // h-bonding interaction
+  double **epsilon_hb, **a_hb, **cut_hb_0, **cut_hb_c, **cut_hb_lo, **cut_hb_hi;
+  double **cut_hb_lc, **cut_hb_hc, **b_hb_lo, **b_hb_hi, **shift_hb;
+  double **cutsq_hb_hc;
+
+  double **a_hb1, **theta_hb1_0, **dtheta_hb1_ast;
+  double **b_hb1, **dtheta_hb1_c;
+
+  double **a_hb2, **theta_hb2_0, **dtheta_hb2_ast;
+  double **b_hb2, **dtheta_hb2_c;
+
+  double **a_hb3, **theta_hb3_0, **dtheta_hb3_ast;
+  double **b_hb3, **dtheta_hb3_c;
+
+  double **a_hb4, **theta_hb4_0, **dtheta_hb4_ast;
+  double **b_hb4, **dtheta_hb4_c;
+
+  double **a_hb7, **theta_hb7_0, **dtheta_hb7_ast;
+  double **b_hb7, **dtheta_hb7_c;
+
+  double **a_hb8, **theta_hb8_0, **dtheta_hb8_ast;
+  double **b_hb8, **dtheta_hb8_c;
+
+  virtual void allocate();
+};
+
+}
+
+#endif
+#endif
+
+/* ERROR/WARNING messages:
+
+E: Illegal ... command
+
+Self-explanatory.  Check the input script syntax and compare to the
+documentation for the command.  You can use -echo screen as a
+command-line option when running LAMMPS to see the offending line.
+
+E: Incorrect args for pair coefficients
+
+Self-explanatory.  Check the input script or data file.
+
+*/
--- a/src/USER-CGDNA/pair_oxdna_stk.cpp
+++ b/src/USER-CGDNA/pair_oxdna_stk.cpp
--- a/src/USER-CGDNA/pair_oxdna_stk.h
+++ b/src/USER-CGDNA/pair_oxdna_stk.h
@ -0,0 +1,83 @@
+/* ----------------------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   http://lammps.sandia.gov, Sandia National Laboratories
+   Steve Plimpton, sjplimp@sandia.gov
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+/* ----------------------------------------------------------------------
+   Contributing author: Oliver Henrich (EPCC, University of Edinburgh)
+------------------------------------------------------------------------- */
+
+#ifdef PAIR_CLASS
+
+PairStyle(oxdna_stk,PairOxdnaStk)
+
+#else
+
+#ifndef LMP_PAIR_OXDNA_STK_H
+#define LMP_PAIR_OXDNA_STK_H
+
+#include "pair.h"
+
+namespace LAMMPS_NS {
+
+class PairOxdnaStk : public Pair {
+ public:
+  PairOxdnaStk(class LAMMPS *);
+  virtual ~PairOxdnaStk();
+  virtual void compute(int, int);
+  void settings(int, char **);
+  void coeff(int, char **);
+  void init_style();
+  void init_list(int, class NeighList *);
+  double init_one(int, int);
+  void write_restart(FILE *);
+  void read_restart(FILE *);
+  void write_restart_settings(FILE *);
+  void read_restart_settings(FILE *);
+  void write_data(FILE *);
+  void write_data_all(FILE *);
+  void *extract(const char *, int &);
+
+ protected:
+  // stacking interaction
+  double **epsilon_st, **a_st, **cut_st_0, **cut_st_c;
+  double **cut_st_lo, **cut_st_hi;
+  double **cut_st_lc, **cut_st_hc, **b_st_lo, **b_st_hi, **shift_st;
+  double **cutsq_st_hc;
+  double **a_st4, **theta_st4_0, **dtheta_st4_ast;
+  double **b_st4, **dtheta_st4_c;
+  double **a_st5, **theta_st5_0, **dtheta_st5_ast;
+  double **b_st5, **dtheta_st5_c;
+  double **a_st6, **theta_st6_0, **dtheta_st6_ast;
+  double **b_st6, **dtheta_st6_c;
+  double **a_st1, **cosphi_st1_ast, **b_st1, **cosphi_st1_c;
+  double **a_st2, **cosphi_st2_ast, **b_st2, **cosphi_st2_c;
+
+  virtual void allocate();
+};
+
+}
+
+#endif
+#endif
+
+/* ERROR/WARNING messages:
+
+E: Illegal ... command
+
+Self-explanatory.  Check the input script syntax and compare to the
+documentation for the command.  You can use -echo screen as a
+command-line option when running LAMMPS to see the offending line.
+
+E: Incorrect args for pair coefficients
+
+Self-explanatory.  Check the input script or data file.
+
+*/
--- a/src/USER-CGDNA/pair_oxdna_xstk.cpp
+++ b/src/USER-CGDNA/pair_oxdna_xstk.cpp
--- a/src/USER-CGDNA/pair_oxdna_xstk.h
+++ b/src/USER-CGDNA/pair_oxdna_xstk.h
@ -0,0 +1,92 @@
+/* ----------------------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   http://lammps.sandia.gov, Sandia National Laboratories
+   Steve Plimpton, sjplimp@sandia.gov
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+/* ----------------------------------------------------------------------
+   Contributing author: Oliver Henrich (EPCC, University of Edinburgh)
+------------------------------------------------------------------------- */
+
+#ifdef PAIR_CLASS
+
+PairStyle(oxdna_xstk,PairOxdnaXstk)
+
+#else
+
+#ifndef LMP_PAIR_OXDNA_XSTK_H
+#define LMP_PAIR_OXDNA_XSTK_H
+
+#include "pair.h"
+
+namespace LAMMPS_NS {
+
+class PairOxdnaXstk : public Pair {
+ public:
+  PairOxdnaXstk(class LAMMPS *);
+  virtual ~PairOxdnaXstk();
+  virtual void compute(int, int);
+  void settings(int, char **);
+  void coeff(int, char **);
+  void init_style();
+  void init_list(int, class NeighList *);
+  double init_one(int, int);
+  void write_restart(FILE *);
+  void read_restart(FILE *);
+  void write_restart_settings(FILE *);
+  void read_restart_settings(FILE *);
+  void write_data(FILE *);
+  void write_data_all(FILE *);
+  void *extract(const char *, int &);
+
+ protected:
+  // cross-stacking interaction
+  double **k_xst, **cut_xst_0, **cut_xst_c, **cut_xst_lo, **cut_xst_hi;
+  double **cut_xst_lc, **cut_xst_hc, **b_xst_lo, **b_xst_hi;
+  double **cutsq_xst_hc;
+
+  double **a_xst1, **theta_xst1_0, **dtheta_xst1_ast;
+  double **b_xst1, **dtheta_xst1_c;
+
+  double **a_xst2, **theta_xst2_0, **dtheta_xst2_ast;
+  double **b_xst2, **dtheta_xst2_c;
+
+  double **a_xst3, **theta_xst3_0, **dtheta_xst3_ast;
+  double **b_xst3, **dtheta_xst3_c;
+
+  double **a_xst4, **theta_xst4_0, **dtheta_xst4_ast;
+  double **b_xst4, **dtheta_xst4_c;
+
+  double **a_xst7, **theta_xst7_0, **dtheta_xst7_ast;
+  double **b_xst7, **dtheta_xst7_c;
+
+  double **a_xst8, **theta_xst8_0, **dtheta_xst8_ast;
+  double **b_xst8, **dtheta_xst8_c;
+
+  virtual void allocate();
+};
+
+}
+
+#endif
+#endif
+
+/* ERROR/WARNING messages:
+
+E: Illegal ... command
+
+Self-explanatory.  Check the input script syntax and compare to the
+documentation for the command.  You can use -echo screen as a
+command-line option when running LAMMPS to see the offending line.
+
+E: Incorrect args for pair coefficients
+
+Self-explanatory.  Check the input script or data file.
+
+*/
--- a/src/USER-DPD/nbin_ssa.cpp
+++ b/src/USER-DPD/nbin_ssa.cpp
@ -18,6 +18,7 @@

 #include "nbin_ssa.h"
 #include "atom.h"
+#include "update.h"
 #include "group.h"
 #include "memory.h"
 #include "error.h"
@ -59,6 +60,8 @@ void NBinSSA::bin_atoms()
  int *mask = atom->mask;
  int *ssaAIR = atom->ssaAIR;

+  last_bin = update->ntimestep;
+
  for (i = 0; i < mbins; i++) {
    gbinhead_ssa[i] = -1;
    binhead_ssa[i] = -1;
@ -126,4 +129,3 @@ bigint NBinSSA::memory_usage()
  }
  return bytes;
 }
-
--- a/src/USER-DPD/npair_half_bin_newton_ssa.cpp
+++ b/src/USER-DPD/npair_half_bin_newton_ssa.cpp
@ -251,8 +251,8 @@ void NPairHalfBinNewtonSSA::build(NeighList *list)

 static int cmp_ssaAIR(const void *iptr, const void *jptr)
 {
-  int i = *((int *) iptr);
-  int j = *((int *) jptr);
+  int i = NEIGHMASK & *((int *) iptr);
+  int j = NEIGHMASK & *((int *) jptr);
  if (ssaAIRptr[i] < ssaAIRptr[j]) return -1;
  if (ssaAIRptr[i] > ssaAIRptr[j]) return 1;
  return 0;
--- a/src/USER-DPD/npair_halffull_newton_ssa.cpp
+++ b/src/USER-DPD/npair_halffull_newton_ssa.cpp
@ -123,8 +123,8 @@ void NPairHalffullNewtonSSA::build(NeighList *list)

 static int cmp_ssaAIR(const void *iptr, const void *jptr)
 {
-  int i = *((int *) iptr);
-  int j = *((int *) jptr);
+  int i = NEIGHMASK & *((int *) iptr);
+  int j = NEIGHMASK & *((int *) jptr);
  if (ssaAIRptr[i] < ssaAIRptr[j]) return -1;
  if (ssaAIRptr[i] > ssaAIRptr[j]) return 1;
  return 0;
--- a/src/USER-DPD/pair_dpd_fdt.cpp
+++ b/src/USER-DPD/pair_dpd_fdt.cpp
@ -43,6 +43,8 @@ using namespace LAMMPS_NS;
 PairDPDfdt::PairDPDfdt(LAMMPS *lmp) : Pair(lmp)
 {
  random = NULL;
+  splitFDT_flag = false;
+  a0_is_zero = false;
 }

 /* ---------------------------------------------------------------------- */
@ -94,7 +96,7 @@ void PairDPDfdt::compute(int eflag, int vflag)
  // loop over neighbors of my atoms

  if (splitFDT_flag) {
-    for (ii = 0; ii < inum; ii++) {
+    if (!a0_is_zero) for (ii = 0; ii < inum; ii++) {
      i = ilist[ii];
      xtmp = x[i][0];
      ytmp = x[i][1];
@ -287,6 +289,8 @@ void PairDPDfdt::coeff(int narg, char **arg)
  double sigma_one = force->numeric(FLERR,arg[3]);
  double cut_one = cut_global;

+  a0_is_zero = (a0_one == 0.0); // Typical use with SSA is to set a0 to zero
+
  if (narg == 5) cut_one = force->numeric(FLERR,arg[4]);

  int count = 0;
@ -371,6 +375,7 @@ void PairDPDfdt::read_restart(FILE *fp)

  allocate();

+  a0_is_zero = true; // start with assumption that a0 is zero
  int i,j;
  int me = comm->me;
  for (i = 1; i <= atom->ntypes; i++)
@ -386,6 +391,7 @@ void PairDPDfdt::read_restart(FILE *fp)
        MPI_Bcast(&a0[i][j],1,MPI_DOUBLE,0,world);
        MPI_Bcast(&sigma[i][j],1,MPI_DOUBLE,0,world);
        MPI_Bcast(&cut[i][j],1,MPI_DOUBLE,0,world);
+        a0_is_zero = a0_is_zero && (a0[i][j] == 0.0); // verify the zero assumption
      }
    }
 }
--- a/src/USER-DPD/pair_dpd_fdt.h
+++ b/src/USER-DPD/pair_dpd_fdt.h
@ -50,6 +50,7 @@ class PairDPDfdt : public Pair {
  double cut_global;
  int seed;
  bool splitFDT_flag;
+  bool a0_is_zero;

  void allocate();

--- a/src/USER-DPD/pair_dpd_fdt_energy.cpp
+++ b/src/USER-DPD/pair_dpd_fdt_energy.cpp
@ -46,6 +46,7 @@ PairDPDfdtEnergy::PairDPDfdtEnergy(LAMMPS *lmp) : Pair(lmp)
  duCond = NULL;
  duMech = NULL;
  splitFDT_flag = false;
+  a0_is_zero = false;

  comm_reverse = 2;
 }
@ -110,7 +111,7 @@ void PairDPDfdtEnergy::compute(int eflag, int vflag)
  // loop over neighbors of my atoms

  if (splitFDT_flag) {
-    for (ii = 0; ii < inum; ii++) {
+    if (!a0_is_zero) for (ii = 0; ii < inum; ii++) {
      i = ilist[ii];
      xtmp = x[i][0];
      ytmp = x[i][1];
@ -373,6 +374,8 @@ void PairDPDfdtEnergy::coeff(int narg, char **arg)
  double cut_one = cut_global;
  double kappa_one;

+  a0_is_zero = (a0_one == 0.0); // Typical use with SSA is to set a0 to zero
+
  kappa_one = force->numeric(FLERR,arg[4]);
  if (narg == 6) cut_one = force->numeric(FLERR,arg[5]);

@ -466,6 +469,7 @@ void PairDPDfdtEnergy::read_restart(FILE *fp)

  allocate();

+  a0_is_zero = true; // start with assumption that a0 is zero
  int i,j;
  int me = comm->me;
  for (i = 1; i <= atom->ntypes; i++)
@ -483,6 +487,7 @@ void PairDPDfdtEnergy::read_restart(FILE *fp)
        MPI_Bcast(&sigma[i][j],1,MPI_DOUBLE,0,world);
        MPI_Bcast(&kappa[i][j],1,MPI_DOUBLE,0,world);
        MPI_Bcast(&cut[i][j],1,MPI_DOUBLE,0,world);
+        a0_is_zero = a0_is_zero && (a0[i][j] == 0.0); // verify the zero assumption
      }
    }
 }
--- a/src/USER-DPD/pair_dpd_fdt_energy.h
+++ b/src/USER-DPD/pair_dpd_fdt_energy.h
@ -52,6 +52,7 @@ class PairDPDfdtEnergy : public Pair {
  double cut_global;
  int seed;
  bool splitFDT_flag;
+  bool a0_is_zero;

  void allocate();

--- a/src/USER-INTEL/nbin_intel.cpp
+++ b/src/USER-INTEL/nbin_intel.cpp
@ -86,7 +86,6 @@ void NBinIntel::bin_atoms_setup(int nall)
         nocopy(binhead:length(maxbin+1) alloc_if(1) free_if(0))
    }
    #endif
-    last_bin_memory = update->ntimestep;
  }

  // bins = per-atom vector
@ -127,11 +126,7 @@ void NBinIntel::bin_atoms_setup(int nall)
      _fix->get_single_buffers()->set_bininfo(_atombin,_binpacked);
    else
      _fix->get_double_buffers()->set_bininfo(_atombin,_binpacked);
-
-    last_bin_memory = update->ntimestep;
  }
-
-  last_bin = update->ntimestep;
 }

 /* ----------------------------------------------------------------------
@ -140,6 +135,8 @@ void NBinIntel::bin_atoms_setup(int nall)

 void NBinIntel::bin_atoms()
 {
+  last_bin = update->ntimestep;
+
  if (_precision_mode == FixIntel::PREC_MODE_MIXED)
    bin_atoms(_fix->get_mixed_buffers());
  else if (_precision_mode == FixIntel::PREC_MODE_SINGLE)
--- a/src/USER-INTEL/npair_intel.cpp
+++ b/src/USER-INTEL/npair_intel.cpp
@ -46,9 +46,7 @@ NPairIntel::~NPairIntel() {
  #endif
 }

-/* ----------------------------------------------------------------------
-   copy needed info from NStencil class to this build class
------------------------------------------------------------------------- */
+/* ---------------------------------------------------------------------- */

 #ifdef _LMP_INTEL_OFFLOAD
 void NPairIntel::grow_stencil()
--- a/src/compute_coord_atom.cpp
+++ b/src/compute_coord_atom.cpp
@ -81,11 +81,11 @@ ComputeCoordAtom::ComputeCoordAtom(LAMMPS *lmp, int narg, char **arg) :
    if (iorientorder < 0)
      error->all(FLERR,"Could not find compute coord/atom compute ID");
    if (strcmp(modify->compute[iorientorder]->style,"orientorder/atom") != 0)
-      error->all(FLERR,"Compute coord/atom compute ID does not compute orientorder/atom");
+      error->all(FLERR,"Compute coord/atom compute ID is not orientorder/atom");

    threshold = force->numeric(FLERR,arg[5]);
    if (threshold <= -1.0 || threshold >= 1.0)
-      error->all(FLERR,"Compute coord/atom threshold value must lie between -1 and 1");
+      error->all(FLERR,"Compute coord/atom threshold not between -1 and 1");

    ncol = 1;
    typelo = new int[ncol];
@ -126,7 +126,7 @@ void ComputeCoordAtom::init()
    comm_forward = 2*(2*l+1);
    if (c_orientorder->iqlcomp < 0)
      error->all(FLERR,"Compute coord/atom requires components "
-                 "option in compute orientorder/atom be defined");
+                 "option in compute orientorder/atom");
  }

  if (force->pair == NULL)
@ -169,9 +169,6 @@ void ComputeCoordAtom::compute_peratom()

  invoked_peratom = update->ntimestep;

-//  printf("Number of degrees %i components degree %i",nqlist,l);
-//  printf("Particle \t %i \t Norm \t %g \n",0,norm[0][0]);
-
  // grow coordination array if necessary

  if (atom->nmax > nmax) {
@ -195,8 +192,6 @@ void ComputeCoordAtom::compute_peratom()
    }
    nqlist = c_orientorder->nqlist;
    int ltmp = l;
-//    l = c_orientorder->qlcomp;
-    if (ltmp != l) error->all(FLERR,"Debug error, ltmp != l\n");
    normv = c_orientorder->array_atom;
    comm->forward_comm_compute(this);
  }
@ -317,7 +312,7 @@ void ComputeCoordAtom::compute_peratom()
 /* ---------------------------------------------------------------------- */

 int ComputeCoordAtom::pack_forward_comm(int n, int *list, double *buf,
-                                  int pbc_flag, int *pbc)
+                                        int pbc_flag, int *pbc)
 {
  int i,m=0,j;
  for (i = 0; i < n; ++i) {
@ -340,7 +335,6 @@ void ComputeCoordAtom::unpack_forward_comm(int n, int first, double *buf)
      normv[i][j] = buf[m++];
    }
  }
-
 }

 /* ----------------------------------------------------------------------
--- a/src/compute_orientorder_atom.cpp
+++ b/src/compute_orientorder_atom.cpp
@ -73,7 +73,8 @@ ComputeOrientOrderAtom::ComputeOrientOrderAtom(LAMMPS *lmp, int narg, char **arg
  int iarg = 3;
  while (iarg < narg) {
    if (strcmp(arg[iarg],"nnn") == 0) {
-      if (iarg+2 > narg) error->all(FLERR,"Illegal compute orientorder/atom command");
+      if (iarg+2 > narg) 
+        error->all(FLERR,"Illegal compute orientorder/atom command");
      if (strcmp(arg[iarg+1],"NULL") == 0) {
        nnn = 0;
      } else {
@ -91,7 +92,8 @@ ComputeOrientOrderAtom::ComputeOrientOrderAtom(LAMMPS *lmp, int narg, char **arg
      memory->destroy(qlist);
      memory->create(qlist,nqlist,"orientorder/atom:qlist");
      iarg += 2;
-      if (iarg+nqlist > narg) error->all(FLERR,"Illegal compute orientorder/atom command");
+      if (iarg+nqlist > narg) 
+        error->all(FLERR,"Illegal compute orientorder/atom command");
      qmax = 0;
      for (int iw = 0; iw < nqlist; iw++) {
        qlist[iw] = force->numeric(FLERR,arg[iarg+iw]);
@ -157,11 +159,12 @@ ComputeOrientOrderAtom::~ComputeOrientOrderAtom()
 void ComputeOrientOrderAtom::init()
 {
  if (force->pair == NULL)
-    error->all(FLERR,"Compute orientorder/atom requires a pair style be defined");
+    error->all(FLERR,"Compute orientorder/atom requires a "
+               "pair style be defined");
  if (cutsq == 0.0) cutsq = force->pair->cutforce * force->pair->cutforce;
  else if (sqrt(cutsq) > force->pair->cutforce)
-    error->all(FLERR,
-               "Compute orientorder/atom cutoff is longer than pairwise cutoff");
+    error->all(FLERR,"Compute orientorder/atom cutoff is "
+               "longer than pairwise cutoff");

  memory->create(qnm_r,qmax,2*qmax+1,"orientorder/atom:qnm_r");
  memory->create(qnm_i,qmax,2*qmax+1,"orientorder/atom:qnm_i");
@ -488,7 +491,8 @@ void ComputeOrientOrderAtom::calc_boop(double **rlist,
   calculate scalar distance
 ------------------------------------------------------------------------- */

-double ComputeOrientOrderAtom::dist(const double r[]) {
+double ComputeOrientOrderAtom::dist(const double r[])
+{
  return sqrt(r[0]*r[0] + r[1]*r[1] + r[2]*r[2]);
 }

@ -497,8 +501,8 @@ double ComputeOrientOrderAtom::dist(const double r[]) {
   Y_l^m (theta, phi) = prefactor(l, m, cos(theta)) * exp(i*m*phi)
 ------------------------------------------------------------------------- */

-double ComputeOrientOrderAtom::
-polar_prefactor(int l, int m, double costheta) {
+double ComputeOrientOrderAtom::polar_prefactor(int l, int m, double costheta) 
+{
  const int mabs = abs(m);

  double prefactor = 1.0;
@ -517,8 +521,8 @@ polar_prefactor(int l, int m, double costheta) {
   associated legendre polynomial
 ------------------------------------------------------------------------- */

-double ComputeOrientOrderAtom::
-associated_legendre(int l, int m, double x) {
+double ComputeOrientOrderAtom::associated_legendre(int l, int m, double x) 
+{
  if (l < m) return 0.0;

  double p(1.0), pm1(0.0), pm2(0.0);
--- a/src/compute_rdf.cpp
+++ b/src/compute_rdf.cpp
@ -40,7 +40,8 @@ using namespace MathConst;
 ComputeRDF::ComputeRDF(LAMMPS *lmp, int narg, char **arg) :
  Compute(lmp, narg, arg),
  rdfpair(NULL), nrdfpair(NULL), ilo(NULL), ihi(NULL), jlo(NULL), jhi(NULL),
-  hist(NULL), histall(NULL), typecount(NULL), icount(NULL), jcount(NULL), duplicates(NULL)
+  hist(NULL), histall(NULL), typecount(NULL), icount(NULL), jcount(NULL), 
+  duplicates(NULL)
 {
  if (narg < 4 || (narg-4) % 2) error->all(FLERR,"Illegal compute rdf command");

--- a/src/nbin.cpp
+++ b/src/nbin.cpp
@ -24,7 +24,7 @@ using namespace LAMMPS_NS;

 NBin::NBin(LAMMPS *lmp) : Pointers(lmp)
 {
-  last_setup = last_bin = last_bin_memory = -1;
+  last_bin = -1;
  maxbin = maxatom = 0;
  binhead = NULL;
  bins = NULL;
@ -71,7 +71,6 @@ void NBin::bin_atoms_setup(int nall)
    maxbin = mbins;
    memory->destroy(binhead);
    memory->create(binhead,maxbin,"neigh:binhead");
-    last_bin_memory = update->ntimestep;
  }

  // bins = per-atom vector
@ -80,10 +79,7 @@ void NBin::bin_atoms_setup(int nall)
    maxatom = nall;
    memory->destroy(bins);
    memory->create(bins,maxatom,"neigh:bins");
-    last_bin_memory = update->ntimestep;
  }
-
-  last_bin = update->ntimestep;
 }

 /* ----------------------------------------------------------------------
--- a/src/nbin.h
+++ b/src/nbin.h
@ -21,9 +21,7 @@ namespace LAMMPS_NS {
 class NBin : protected Pointers {
 public:
  int istyle;                      // 1-N index into binnames
-
-  bigint last_setup,last_bin;      // timesteps for last operations performed
-  bigint last_bin_memory;
+  bigint last_bin;                 // last timestep atoms were binned

  int nbinx,nbiny,nbinz;           // # of global bins
  int mbins;                       // # of local bins and offset on this proc
--- a/src/nbin_standard.cpp
+++ b/src/nbin_standard.cpp
@ -54,8 +54,6 @@ NBinStandard::NBinStandard(LAMMPS *lmp) : NBin(lmp) {}

 void NBinStandard::setup_bins(int style)
 {
-  last_setup = update->ntimestep;
-
  // bbox = size of bbox of entire domain
  // bsubbox lo/hi = bounding box of my subdomain extended by comm->cutghost
  // for triclinic:
@ -197,6 +195,7 @@ void NBinStandard::bin_atoms()
 {
  int i,ibin;

+  last_bin = update->ntimestep;
  for (i = 0; i < mbins; i++) binhead[i] = -1;

  // bin in reverse order so linked list will be in forward order
--- a/src/neighbor.cpp
+++ b/src/neighbor.cpp
@ -695,27 +695,38 @@ void Neighbor::init_pair()
  }
      
  // (C) rule: 
-  // for fix/compute requests, occasional or not does not matter
+  // for fix/compute requests to be copied:
  // 1st check:
+  // occasional or not does not matter
+  // Kokkos host/device flags must match
+  // SSA flag must match
+  // newton setting must match (else list has different neighbors in it)
+  // 2nd check:
  // if request = half and non-skip/copy pair half/respaouter request exists,
  // or if request = full and non-skip/copy pair full request exists,
  // or if request = gran and non-skip/copy pair gran request exists,
  // then morph to copy of the matching parent list
-  // 2nd check: only if no match to 1st check
+  // 3rd check: only if no match to 1st check
  // if request = half and non-skip/copy pair full request exists,
  // then morph to half_from_full of the matching parent list
  // for 1st or 2nd check, parent can be copy list or pair or fix
  
+  int inewton,jnewton;
+
  for (i = 0; i < nrequest; i++) {
    if (!requests[i]->fix && !requests[i]->compute) continue;
    for (j = 0; j < nrequest; j++) {
-      // Kokkos flags must match
      if (requests[i]->kokkos_device != requests[j]->kokkos_device) continue;
      if (requests[i]->kokkos_host != requests[j]->kokkos_host) continue;

      if (requests[i]->ssa != requests[j]->ssa) continue;
-      // newton 2 and newton 0 both are newton off
-      if ((requests[i]->newton & 2) != (requests[j]->newton & 2)) continue;
+
+      // IJ newton = 1 for newton on, 2 for newton off
+      inewton = requests[i]->newton;
+      if (inewton == 0) inewton = force->newton_pair ? 1 : 2; 
+      jnewton = requests[j]->newton;
+      if (jnewton == 0) jnewton = force->newton_pair ? 1 : 2;
+      if (inewton != jnewton) continue;

      if (requests[i]->half && requests[j]->pair && 
          !requests[j]->skip && requests[j]->half && !requests[j]->copy)
@ -899,33 +910,37 @@ void Neighbor::init_pair()
  }

  // reorder plist vector if necessary
-  // relevant for lists that copy/skip/half-full from parent
+  // relevant for lists that are derived from a parent list:
+  //   half-full,copy,skip
  // the child index must appear in plist after the parent index
  // swap two indices within plist when dependency is mis-ordered
-  // done when entire pass thru plist results in no swaps
-  
+  // start double loop check again whenever a swap is made
+  // done when entire double loop test results in no swaps
+
  NeighList *ptr;

  int done = 0;
  while (!done) {
    done = 1;
    for (i = 0; i < npair_perpetual; i++) {
-      ptr = NULL;
-      if (lists[plist[i]]->listfull) ptr = lists[plist[i]]->listfull;
-      if (lists[plist[i]]->listcopy) ptr = lists[plist[i]]->listcopy;
-      // listskip check must be after listfull check
-      if (lists[plist[i]]->listskip) ptr = lists[plist[i]]->listskip;
-      if (ptr == NULL) continue;
-      for (m = 0; m < nrequest; m++)
-        if (ptr == lists[m]) break;
-      for (j = 0; j < npair_perpetual; j++)
-        if (m == plist[j]) break;
-      if (j < i) continue;
-      int tmp = plist[i];     // swap I,J indices
-      plist[i] = plist[j];
-      plist[j] = tmp;
-      done = 0;
-      break;
+      for (k = 0; k < 3; k++) {
+        ptr = NULL;
+        if (k == 0) ptr = lists[plist[i]]->listcopy;
+        if (k == 1) ptr = lists[plist[i]]->listskip;
+        if (k == 2) ptr = lists[plist[i]]->listfull;
+        if (ptr == NULL) continue;
+        for (m = 0; m < nrequest; m++)
+          if (ptr == lists[m]) break;
+        for (j = 0; j < npair_perpetual; j++)
+          if (m == plist[j]) break;
+        if (j < i) continue;
+        int tmp = plist[i];     // swap I,J indices
+        plist[i] = plist[j];
+        plist[j] = tmp;
+        done = 0;
+        break;
+      }
+      if (!done) break;
    }
  }

@ -1085,7 +1100,7 @@ void Neighbor::init_topology()

 void Neighbor::print_pairwise_info()
 {
-  int i,j,m;
+  int i,m;
  char str[128];
  const char *kind;
  FILE *out;
@ -1154,8 +1169,8 @@ void Neighbor::print_pairwise_info()
        else if (requests[i]->respamiddle) kind = "respa/middle";
        else if (requests[i]->respaouter) kind = "respa/outer";
        else if (requests[i]->half_from_full) kind = "half/from/full";
-        if (requests[i]->occasional) fprintf(out,", occasional");
-        else fprintf(out,", perpetual");
+        if (requests[i]->occasional) fprintf(out,", %s, occasional",kind);
+        else fprintf(out,", %s, perpetual",kind);
        if (requests[i]->ghost) fprintf(out,", ghost");
        if (requests[i]->ssa) fprintf(out,", ssa");
        if (requests[i]->omp) fprintf(out,", omp");
@ -1841,7 +1856,7 @@ void Neighbor::build_one(class NeighList *mylist, int preflag)
  // create stencil if hasn't been created since last setup_bins() call

  NStencil *ns = np->ns;
-  if (ns && ns->last_create < last_setup_bins) {
+  if (ns && ns->last_stencil < last_setup_bins) {
    ns->create_setup();
    ns->create();
  }
@ -1874,29 +1889,22 @@ void Neighbor::set(int narg, char **arg)
 /* ----------------------------------------------------------------------
   reset timestamps in all NeignBin, NStencil, NPair classes
   so that neighbor lists will rebuild properly with timestep change
+   ditto for lastcall and last_setup_bins
 ------------------------------------------------------------------------- */

 void Neighbor::reset_timestep(bigint ntimestep)
 {
-  for (int i = 0; i < nbin; i++) {
-    neigh_bin[i]->last_setup = -1;
+  for (int i = 0; i < nbin; i++)
    neigh_bin[i]->last_bin = -1;
-    neigh_bin[i]->last_bin_memory = -1;
-  }
-
-  for (int i = 0; i < nstencil; i++) {
-    neigh_stencil[i]->last_create = -1;
-    neigh_stencil[i]->last_stencil_memory = -1;
-    neigh_stencil[i]->last_copy_bin = -1;
-  }
-
+  for (int i = 0; i < nstencil; i++)
+    neigh_stencil[i]->last_stencil = -1;
  for (int i = 0; i < nlist; i++) {
    if (!neigh_pair[i]) continue;
    neigh_pair[i]->last_build = -1;
-    neigh_pair[i]->last_copy_bin_setup = -1;
-    neigh_pair[i]->last_copy_bin = -1;
-    neigh_pair[i]->last_copy_stencil = -1;
  }
+
+  lastcall = -1;
+  last_setup_bins = -1;
 }

 /* ----------------------------------------------------------------------
--- a/src/npair.cpp
+++ b/src/npair.cpp
@ -28,8 +28,6 @@ NPair::NPair(LAMMPS *lmp)
  : Pointers(lmp), nb(NULL), ns(NULL), bins(NULL), stencil(NULL)
 {
  last_build = -1;
-  last_copy_bin_setup = last_copy_bin = last_copy_stencil = -1;
-
  molecular = atom->molecular;
 }

@ -76,10 +74,10 @@ void NPair::copy_neighbor_info()
 }

 /* ----------------------------------------------------------------------
-   copy bin geometry info from NBin class to this build class
+   copy info from NBin class to this build class
 ------------------------------------------------------------------------- */

-void NPair::copy_bin_setup_info()
+void NPair::copy_bin_info()
 {
  nbinx = nb->nbinx;
  nbiny = nb->nbiny;
@ -95,20 +93,13 @@ void NPair::copy_bin_setup_info()
  bininvx = nb->bininvx;
  bininvy = nb->bininvy;
  bininvz = nb->bininvz;
-}

-/* ----------------------------------------------------------------------
-   copy per-atom and per-bin vectors from NBin class to this build class
------------------------------------------------------------------------- */
-
-void NPair::copy_bin_info()
-{
  bins = nb->bins;
  binhead = nb->binhead;
 }

 /* ----------------------------------------------------------------------
-   copy needed info from NStencil class to this build class
+   copy info from NStencil class to this build class
 ------------------------------------------------------------------------- */

 void NPair::copy_stencil_info()
@ -122,23 +113,15 @@ void NPair::copy_stencil_info()
 }

 /* ----------------------------------------------------------------------
-   copy needed info from NStencil class to this build class
+   copy info from NBin and NStencil classes to this build class
 ------------------------------------------------------------------------- */

 void NPair::build_setup()
 {
-  if (nb && last_copy_bin_setup < nb->last_setup) {
-    copy_bin_setup_info();
-    last_copy_bin_setup = update->ntimestep;
-  }
-  if (nb && ((last_copy_bin < nb->last_bin_memory) || (bins != nb->bins))) {
-    copy_bin_info();
-    last_copy_bin = update->ntimestep;
-  }
-  if (ns && ((last_copy_stencil < ns->last_create) || (stencil != ns->stencil))) {
-    copy_stencil_info();
-    last_copy_stencil = update->ntimestep;
-  }
+  if (nb) copy_bin_info();
+  if (ns) copy_stencil_info();
+
+  // set here, since build_setup() always called before build()

  last_build = update->ntimestep;
 }
--- a/src/npair.h
+++ b/src/npair.h
@ -23,11 +23,7 @@ class NPair : protected Pointers {
  int istyle;                   // 1-N index into pairnames
  class NBin *nb;               // ptr to NBin instance I depend on
  class NStencil *ns;           // ptr to NStencil instance I depend on
-
  bigint last_build;            // last timestep build performed
-  bigint last_copy_bin_setup;   // last timestep I invoked copy_bin_setup_info()
-  bigint last_copy_bin;         // last step I invoked copy_bin_info()
-  bigint last_copy_stencil;     // last step I invoked copy_bin_stencil_info()

  NPair(class LAMMPS *);
  virtual ~NPair() {}
@ -93,7 +89,6 @@ class NPair : protected Pointers {

  // methods for all NPair variants

-  void copy_bin_setup_info();
  virtual void copy_bin_info();
  virtual void copy_stencil_info();

--- a/src/nstencil.cpp
+++ b/src/nstencil.cpp
@ -56,11 +56,9 @@ enum{NSQ,BIN,MULTI};     // also in Neighbor

 NStencil::NStencil(LAMMPS *lmp) : Pointers(lmp)
 {
-  last_create = last_stencil_memory = -1;
-  last_copy_bin = -1;
+  last_stencil = -1;

  xyzflag = 0;
-
  maxstencil = maxstencil_multi = 0;
  stencil = NULL;
  stencilxyz = NULL;
@ -126,10 +124,8 @@ void NStencil::copy_bin_info()

 void NStencil::create_setup()
 {
-  if (nb && last_copy_bin < nb->last_setup) {
-    copy_bin_info();
-    last_copy_bin = update->ntimestep;
-  }
+  if (nb) copy_bin_info();
+  last_stencil = update->ntimestep;

  // sx,sy,sz = max range of stencil in each dim
  // smax = max possible size of entire 3d stencil
@ -157,7 +153,6 @@ void NStencil::create_setup()
        memory->destroy(stencilxyz);
        memory->create(stencilxyz,maxstencil,3,"neighstencil:stencilxyz");
      }
-      last_stencil_memory = update->ntimestep;
    }

  } else {
@ -172,7 +167,6 @@ void NStencil::create_setup()
        stencil_multi[i] = NULL;
        distsq_multi[i] = NULL;
      }
-      last_stencil_memory = update->ntimestep;
    }
    if (smax > maxstencil_multi) {
      maxstencil_multi = smax;
@ -183,12 +177,9 @@ void NStencil::create_setup()
                       "neighstencil:stencil_multi");
        memory->create(distsq_multi[i],maxstencil_multi,
                       "neighstencil:distsq_multi");
-        last_stencil_memory = update->ntimestep;
      }
    }
  }
-
-  last_create = update->ntimestep;
 }

 /* ----------------------------------------------------------------------
--- a/src/nstencil.h
+++ b/src/nstencil.h
@ -22,10 +22,7 @@ class NStencil : protected Pointers {
 public:
  int istyle;                      // 1-N index into binnames
  class NBin *nb;                  // ptr to NBin instance I depend on
-
-  bigint last_create;              // timesteps for last operations performed
-  bigint last_stencil_memory;
-  bigint last_copy_bin;
+  bigint last_stencil;             // last timestep stencil was created

  int nstencil;                    // # of bins in stencil
  int *stencil;                    // list of bin offsets
--- a/src/version.h
+++ b/src/version.h
@ -1 +1 @@
-#define LAMMPS_VERSION "17 Jan 2017"
+#define LAMMPS_VERSION "20 Jan 2017"
Author	SHA1	Message	Date
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