From f96520ffc65c8f624dc37424ffcd18978dc7e7db Mon Sep 17 00:00:00 2001
From: Oliver Henrich <ohenrich@users.noreply.github.com>
Date: Fri, 15 Nov 2019 17:49:17 +0000
Subject: [PATCH] Added unique base pairing

---
 examples/USER/cgdna/README                    |  17 +
 .../oxDNA2/unique_bp/data.duplex4.4type       | 130 +++
 .../oxDNA2/unique_bp/data.duplex4.8type       | 130 +++
 .../oxDNA2/unique_bp/generate_unique.py       | 828 ++++++++++++++++++
 .../oxDNA2/unique_bp/in.duplex4.4type         |  94 ++
 .../oxDNA2/unique_bp/in.duplex4.8type         |  94 ++
 .../unique_bp/log.07Aug19.duplex4.4type.g++1  | 232 +++++
 .../unique_bp/log.07Aug19.duplex4.4type.g++4  | 232 +++++
 .../unique_bp/log.07Aug19.duplex4.8type.g++1  | 274 ++++++
 .../unique_bp/log.07Aug19.duplex4.8type.g++4  | 274 ++++++
 10 files changed, 2305 insertions(+)
 create mode 100644 examples/USER/cgdna/examples/oxDNA2/unique_bp/data.duplex4.4type
 create mode 100644 examples/USER/cgdna/examples/oxDNA2/unique_bp/data.duplex4.8type
 create mode 100644 examples/USER/cgdna/examples/oxDNA2/unique_bp/generate_unique.py
 create mode 100644 examples/USER/cgdna/examples/oxDNA2/unique_bp/in.duplex4.4type
 create mode 100644 examples/USER/cgdna/examples/oxDNA2/unique_bp/in.duplex4.8type
 create mode 100644 examples/USER/cgdna/examples/oxDNA2/unique_bp/log.07Aug19.duplex4.4type.g++1
 create mode 100644 examples/USER/cgdna/examples/oxDNA2/unique_bp/log.07Aug19.duplex4.4type.g++4
 create mode 100644 examples/USER/cgdna/examples/oxDNA2/unique_bp/log.07Aug19.duplex4.8type.g++1
 create mode 100644 examples/USER/cgdna/examples/oxDNA2/unique_bp/log.07Aug19.duplex4.8type.g++4

diff --git a/examples/USER/cgdna/README b/examples/USER/cgdna/README
index a03490b630..7cbf76fd5a 100644
--- a/examples/USER/cgdna/README
+++ b/examples/USER/cgdna/README
@@ -55,6 +55,23 @@ are sequence-dependent (cf. keyword 'seqdep' in according pair styles).
 
 /******************************************************************************/
 
+/examples/oxDNA2/unique_bp:
+
+This example uses atom types 1-8 to model a 13 base pair duplex.
+The nucleotide types are assigned as follows:
+A = 1,5; C = 2,6; G = 3,7; T = 4,8
+When a large number of atom types is used, this feature allows
+quasi-unique base pairing between two individual nucleotides.
+
+The topology is
+A   C   G   T   A   C   G   T   A   C   G   T   A
+1 - 2 - 3 - 4 - 5 - 6 - 7 - 8 - 1 - 2 - 7 - 8 - 1
+|   |   |   |   |   |   |   |   |   |   |   |   |
+4 - 3 - 2 - 1 - 8 - 7 - 6 - 5 - 4 - 3 - 6 - 5 - 4
+T   G   C   A   T   G   C   A   T   G   C   A   T
+
+/******************************************************************************/
+
 /examples/oxRNA2/duplex4
 
 This is the duplex2 run with the oxRNA2 force field instead of the oxDNA or 
diff --git a/examples/USER/cgdna/examples/oxDNA2/unique_bp/data.duplex4.4type b/examples/USER/cgdna/examples/oxDNA2/unique_bp/data.duplex4.4type
new file mode 100644
index 0000000000..a8412eef07
--- /dev/null
+++ b/examples/USER/cgdna/examples/oxDNA2/unique_bp/data.duplex4.4type
@@ -0,0 +1,130 @@
+# LAMMPS data file
+26 atoms
+26 ellipsoids
+24 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-ID, type, position, molecule-ID, ellipsoid flag, density
+Atoms
+
+1 1 -6.000000000000001e-01  0.000000000000000e+00  0.000000000000000e+00 1 1 1
+2 2 -4.957432645895970e-01 -3.379920348381733e-01  3.897628551303122e-01 1 1 1
+3 3 -2.192046146198370e-01 -5.585242491865227e-01  7.795257102606244e-01 1 1 1
+4 4  1.335125603737887e-01 -5.849567473090943e-01  1.169288565390937e+00 1 1 1
+5 1  4.398311230978960e-01 -4.081036426625517e-01  1.559051420521249e+00 1 1 1
+6 2  5.932984957350773e-01 -8.942535970570469e-02  1.948814275651561e+00 1 1 1
+7 3  5.405813207414517e-01  2.603302434705350e-01  2.338577130781873e+00 1 1 1
+8 4  3.000000000000002e-01  5.196152422706634e-01  2.728339985912185e+00 1 1 1
+9 1 -4.483805615185452e-02  5.983222783087083e-01  3.118102841042497e+00 1 1 1
+10 2 -3.740938811152403e-01  4.690988894808181e-01  3.507865696172809e+00 1 1 1
+11 3 -5.733436834716847e-01  1.768531046465427e-01  3.897628551303121e+00 1 1 1
+12 4 -5.733436834716849e-01 -1.768531046465427e-01  4.287391406433434e+00 1 1 1
+13 1 -3.740938811152403e-01 -4.690988894808182e-01  4.677154261563746e+00 1 1 1
+14 4  3.740938811152403e-01  4.690988894808182e-01  4.677154261563746e+00 2 1 1
+15 1  5.733436834716849e-01  1.768531046465427e-01  4.287391406433434e+00 2 1 1
+16 2  5.733436834716849e-01 -1.768531046465426e-01  3.897628551303122e+00 2 1 1
+17 3  3.740938811152403e-01 -4.690988894808181e-01  3.507865696172810e+00 2 1 1
+18 4  4.483805615185462e-02 -5.983222783087085e-01  3.118102841042498e+00 2 1 1
+19 1 -3.000000000000003e-01 -5.196152422706636e-01  2.728339985912186e+00 2 1 1
+20 2 -5.405813207414519e-01 -2.603302434705351e-01  2.338577130781874e+00 2 1 1
+21 3 -5.932984957350776e-01  8.942535970570474e-02  1.948814275651561e+00 2 1 1
+22 4 -4.398311230978962e-01  4.081036426625520e-01  1.559051420521249e+00 2 1 1
+23 1 -1.335125603737888e-01  5.849567473090947e-01  1.169288565390937e+00 2 1 1
+24 2  2.192046146198373e-01  5.585242491865231e-01  7.795257102606246e-01 2 1 1
+25 3  4.957432645895974e-01  3.379920348381736e-01  3.897628551303123e-01 2 1 1
+26 4  6.000000000000006e-01  0.000000000000000e+00  1.110223024625157e-16 2 1 1
+
+# Atom-ID, translational, rotational velocity
+Velocities
+
+1  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00
+2  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00
+3  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00
+4  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00
+5  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00
+6  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00
+7  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00
+8  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00
+9  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00
+10  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00
+11  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00
+12  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00
+13  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00
+14  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00
+15  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00
+16  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00
+17  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00
+18  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00
+19  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00
+20  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00
+21  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00
+22  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00
+23  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00
+24  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00
+25  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00
+26  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00
+
+# Atom-ID, shape, quaternion
+Ellipsoids
+
+1  1.173984503142341e+00  1.173984503142341e+00  1.173984503142341e+00  1.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00
+2  1.173984503142341e+00  1.173984503142341e+00  1.173984503142341e+00  9.555728057861408e-01  0.000000000000000e+00  0.000000000000000e+00  2.947551744109042e-01
+3  1.173984503142341e+00  1.173984503142341e+00  1.173984503142341e+00  8.262387743159949e-01  0.000000000000000e+00  0.000000000000000e+00  5.633200580636221e-01
+4  1.173984503142341e+00  1.173984503142341e+00  1.173984503142341e+00  6.234898018587335e-01  0.000000000000000e+00  0.000000000000000e+00  7.818314824680299e-01
+5  1.173984503142341e+00  1.173984503142341e+00  1.173984503142341e+00  3.653410243663950e-01  0.000000000000000e+00  0.000000000000000e+00  9.308737486442042e-01
+6  1.173984503142341e+00  1.173984503142341e+00  1.173984503142341e+00  7.473009358642424e-02  0.000000000000000e+00  0.000000000000000e+00  9.972037971811802e-01
+7  1.173984503142341e+00  1.173984503142341e+00  1.173984503142341e+00 -2.225209339563144e-01  0.000000000000000e+00  0.000000000000000e+00  9.749279121818237e-01
+8  1.173984503142341e+00  1.173984503142341e+00  1.173984503142341e+00 -5.000000000000001e-01  0.000000000000000e+00  0.000000000000000e+00  8.660254037844387e-01
+9  1.173984503142341e+00  1.173984503142341e+00  1.173984503142341e+00  7.330518718298263e-01 -0.000000000000000e+00  0.000000000000000e+00 -6.801727377709196e-01
+10  1.173984503142341e+00  1.173984503142341e+00  1.173984503142341e+00  9.009688679024190e-01 -0.000000000000000e+00  0.000000000000000e+00 -4.338837391175581e-01
+11  1.173984503142341e+00  1.173984503142341e+00  1.173984503142341e+00  9.888308262251286e-01 -0.000000000000000e+00  0.000000000000000e+00 -1.490422661761745e-01
+12  1.173984503142341e+00  1.173984503142341e+00  1.173984503142341e+00  9.888308262251286e-01  0.000000000000000e+00  0.000000000000000e+00  1.490422661761745e-01
+13  1.173984503142341e+00  1.173984503142341e+00  1.173984503142341e+00  9.009688679024190e-01  0.000000000000000e+00  0.000000000000000e+00  4.338837391175582e-01
+14  1.173984503142341e+00  1.173984503142341e+00  1.173984503142341e+00  0.000000000000000e+00 -4.338837391175582e-01  9.009688679024190e-01  0.000000000000000e+00
+15  1.173984503142341e+00  1.173984503142341e+00  1.173984503142341e+00 -0.000000000000000e+00 -1.490422661761746e-01  9.888308262251286e-01  0.000000000000000e+00
+16  1.173984503142341e+00  1.173984503142341e+00  1.173984503142341e+00 -0.000000000000000e+00  1.490422661761745e-01  9.888308262251286e-01 -0.000000000000000e+00
+17  1.173984503142341e+00  1.173984503142341e+00  1.173984503142341e+00 -0.000000000000000e+00  4.338837391175581e-01  9.009688679024190e-01 -0.000000000000000e+00
+18  1.173984503142341e+00  1.173984503142341e+00  1.173984503142341e+00  0.000000000000000e+00  6.801727377709192e-01  7.330518718298267e-01  0.000000000000000e+00
+19  1.173984503142341e+00  1.173984503142341e+00  1.173984503142341e+00  0.000000000000000e+00  8.660254037844386e-01  5.000000000000001e-01  0.000000000000000e+00
+20  1.173984503142341e+00  1.173984503142341e+00  1.173984503142341e+00  0.000000000000000e+00  9.749279121818235e-01  2.225209339563145e-01  0.000000000000000e+00
+21  1.173984503142341e+00  1.173984503142341e+00  1.173984503142341e+00  0.000000000000000e+00  9.972037971811801e-01 -7.473009358642428e-02  0.000000000000000e+00
+22  1.173984503142341e+00  1.173984503142341e+00  1.173984503142341e+00  0.000000000000000e+00  9.308737486442041e-01 -3.653410243663952e-01  0.000000000000000e+00
+23  1.173984503142341e+00  1.173984503142341e+00  1.173984503142341e+00  0.000000000000000e+00  7.818314824680296e-01 -6.234898018587339e-01  0.000000000000000e+00
+24  1.173984503142341e+00  1.173984503142341e+00  1.173984503142341e+00  0.000000000000000e+00  5.633200580636215e-01 -8.262387743159952e-01  0.000000000000000e+00
+25  1.173984503142341e+00  1.173984503142341e+00  1.173984503142341e+00 -0.000000000000000e+00 -2.947551744109044e-01  9.555728057861407e-01  0.000000000000000e+00
+26  1.173984503142341e+00  1.173984503142341e+00  1.173984503142341e+00 -0.000000000000000e+00  0.000000000000000e+00  1.000000000000000e+00 -0.000000000000000e+00
+
+# Bond topology
+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  8  9
+9  1  9  10
+10  1  10  11
+11  1  11  12
+12  1  12  13
+13  1  14  15
+14  1  15  16
+15  1  16  17
+16  1  17  18
+17  1  18  19
+18  1  19  20
+19  1  20  21
+20  1  21  22
+21  1  22  23
+22  1  23  24
+23  1  24  25
+24  1  25  26
diff --git a/examples/USER/cgdna/examples/oxDNA2/unique_bp/data.duplex4.8type b/examples/USER/cgdna/examples/oxDNA2/unique_bp/data.duplex4.8type
new file mode 100644
index 0000000000..b4d622ac46
--- /dev/null
+++ b/examples/USER/cgdna/examples/oxDNA2/unique_bp/data.duplex4.8type
@@ -0,0 +1,130 @@
+# LAMMPS data file
+26 atoms
+26 ellipsoids
+24 bonds
+
+8 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-ID, type, position, molecule-ID, ellipsoid flag, density
+Atoms
+
+1 1 -6.000000000000001e-01  0.000000000000000e+00  0.000000000000000e+00 1 1 1
+2 2 -4.957432645895970e-01 -3.379920348381733e-01  3.897628551303122e-01 1 1 1
+3 3 -2.192046146198370e-01 -5.585242491865227e-01  7.795257102606244e-01 1 1 1
+4 4  1.335125603737887e-01 -5.849567473090943e-01  1.169288565390937e+00 1 1 1
+5 5  4.398311230978960e-01 -4.081036426625517e-01  1.559051420521249e+00 1 1 1
+6 6  5.932984957350773e-01 -8.942535970570469e-02  1.948814275651561e+00 1 1 1
+7 7  5.405813207414517e-01  2.603302434705350e-01  2.338577130781873e+00 1 1 1
+8 8  3.000000000000002e-01  5.196152422706634e-01  2.728339985912185e+00 1 1 1
+9 1 -4.483805615185452e-02  5.983222783087083e-01  3.118102841042497e+00 1 1 1
+10 2 -3.740938811152403e-01  4.690988894808181e-01  3.507865696172809e+00 1 1 1
+11 7 -5.733436834716847e-01  1.768531046465427e-01  3.897628551303121e+00 1 1 1
+12 8 -5.733436834716849e-01 -1.768531046465427e-01  4.287391406433434e+00 1 1 1
+13 1 -3.740938811152403e-01 -4.690988894808182e-01  4.677154261563746e+00 1 1 1
+14 4  3.740938811152403e-01  4.690988894808182e-01  4.677154261563746e+00 2 1 1
+15 5  5.733436834716849e-01  1.768531046465427e-01  4.287391406433434e+00 2 1 1
+16 6  5.733436834716849e-01 -1.768531046465426e-01  3.897628551303122e+00 2 1 1
+17 3  3.740938811152403e-01 -4.690988894808181e-01  3.507865696172810e+00 2 1 1
+18 4  4.483805615185462e-02 -5.983222783087085e-01  3.118102841042498e+00 2 1 1
+19 5 -3.000000000000003e-01 -5.196152422706636e-01  2.728339985912186e+00 2 1 1
+20 6 -5.405813207414519e-01 -2.603302434705351e-01  2.338577130781874e+00 2 1 1
+21 7 -5.932984957350776e-01  8.942535970570474e-02  1.948814275651561e+00 2 1 1
+22 8 -4.398311230978962e-01  4.081036426625520e-01  1.559051420521249e+00 2 1 1
+23 1 -1.335125603737888e-01  5.849567473090947e-01  1.169288565390937e+00 2 1 1
+24 2  2.192046146198373e-01  5.585242491865231e-01  7.795257102606246e-01 2 1 1
+25 3  4.957432645895974e-01  3.379920348381736e-01  3.897628551303123e-01 2 1 1
+26 4  6.000000000000006e-01  0.000000000000000e+00  1.110223024625157e-16 2 1 1
+
+# Atom-ID, translational, rotational velocity
+Velocities
+
+1  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00
+2  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00
+3  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00
+4  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00
+5  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00
+6  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00
+7  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00
+8  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00
+9  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00
+10  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00
+11  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00
+12  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00
+13  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00
+14  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00
+15  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00
+16  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00
+17  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00
+18  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00
+19  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00
+20  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00
+21  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00
+22  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00
+23  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00
+24  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00
+25  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00
+26  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00
+
+# Atom-ID, shape, quaternion
+Ellipsoids
+
+1  1.173984503142341e+00  1.173984503142341e+00  1.173984503142341e+00  1.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00  0.000000000000000e+00
+2  1.173984503142341e+00  1.173984503142341e+00  1.173984503142341e+00  9.555728057861408e-01  0.000000000000000e+00  0.000000000000000e+00  2.947551744109042e-01
+3  1.173984503142341e+00  1.173984503142341e+00  1.173984503142341e+00  8.262387743159949e-01  0.000000000000000e+00  0.000000000000000e+00  5.633200580636221e-01
+4  1.173984503142341e+00  1.173984503142341e+00  1.173984503142341e+00  6.234898018587335e-01  0.000000000000000e+00  0.000000000000000e+00  7.818314824680299e-01
+5  1.173984503142341e+00  1.173984503142341e+00  1.173984503142341e+00  3.653410243663950e-01  0.000000000000000e+00  0.000000000000000e+00  9.308737486442042e-01
+6  1.173984503142341e+00  1.173984503142341e+00  1.173984503142341e+00  7.473009358642424e-02  0.000000000000000e+00  0.000000000000000e+00  9.972037971811802e-01
+7  1.173984503142341e+00  1.173984503142341e+00  1.173984503142341e+00 -2.225209339563144e-01  0.000000000000000e+00  0.000000000000000e+00  9.749279121818237e-01
+8  1.173984503142341e+00  1.173984503142341e+00  1.173984503142341e+00 -5.000000000000001e-01  0.000000000000000e+00  0.000000000000000e+00  8.660254037844387e-01
+9  1.173984503142341e+00  1.173984503142341e+00  1.173984503142341e+00  7.330518718298263e-01 -0.000000000000000e+00  0.000000000000000e+00 -6.801727377709196e-01
+10  1.173984503142341e+00  1.173984503142341e+00  1.173984503142341e+00  9.009688679024190e-01 -0.000000000000000e+00  0.000000000000000e+00 -4.338837391175581e-01
+11  1.173984503142341e+00  1.173984503142341e+00  1.173984503142341e+00  9.888308262251286e-01 -0.000000000000000e+00  0.000000000000000e+00 -1.490422661761745e-01
+12  1.173984503142341e+00  1.173984503142341e+00  1.173984503142341e+00  9.888308262251286e-01  0.000000000000000e+00  0.000000000000000e+00  1.490422661761745e-01
+13  1.173984503142341e+00  1.173984503142341e+00  1.173984503142341e+00  9.009688679024190e-01  0.000000000000000e+00  0.000000000000000e+00  4.338837391175582e-01
+14  1.173984503142341e+00  1.173984503142341e+00  1.173984503142341e+00  0.000000000000000e+00 -4.338837391175582e-01  9.009688679024190e-01  0.000000000000000e+00
+15  1.173984503142341e+00  1.173984503142341e+00  1.173984503142341e+00 -0.000000000000000e+00 -1.490422661761746e-01  9.888308262251286e-01  0.000000000000000e+00
+16  1.173984503142341e+00  1.173984503142341e+00  1.173984503142341e+00 -0.000000000000000e+00  1.490422661761745e-01  9.888308262251286e-01 -0.000000000000000e+00
+17  1.173984503142341e+00  1.173984503142341e+00  1.173984503142341e+00 -0.000000000000000e+00  4.338837391175581e-01  9.009688679024190e-01 -0.000000000000000e+00
+18  1.173984503142341e+00  1.173984503142341e+00  1.173984503142341e+00  0.000000000000000e+00  6.801727377709192e-01  7.330518718298267e-01  0.000000000000000e+00
+19  1.173984503142341e+00  1.173984503142341e+00  1.173984503142341e+00  0.000000000000000e+00  8.660254037844386e-01  5.000000000000001e-01  0.000000000000000e+00
+20  1.173984503142341e+00  1.173984503142341e+00  1.173984503142341e+00  0.000000000000000e+00  9.749279121818235e-01  2.225209339563145e-01  0.000000000000000e+00
+21  1.173984503142341e+00  1.173984503142341e+00  1.173984503142341e+00  0.000000000000000e+00  9.972037971811801e-01 -7.473009358642428e-02  0.000000000000000e+00
+22  1.173984503142341e+00  1.173984503142341e+00  1.173984503142341e+00  0.000000000000000e+00  9.308737486442041e-01 -3.653410243663952e-01  0.000000000000000e+00
+23  1.173984503142341e+00  1.173984503142341e+00  1.173984503142341e+00  0.000000000000000e+00  7.818314824680296e-01 -6.234898018587339e-01  0.000000000000000e+00
+24  1.173984503142341e+00  1.173984503142341e+00  1.173984503142341e+00  0.000000000000000e+00  5.633200580636215e-01 -8.262387743159952e-01  0.000000000000000e+00
+25  1.173984503142341e+00  1.173984503142341e+00  1.173984503142341e+00 -0.000000000000000e+00 -2.947551744109044e-01  9.555728057861407e-01  0.000000000000000e+00
+26  1.173984503142341e+00  1.173984503142341e+00  1.173984503142341e+00 -0.000000000000000e+00  0.000000000000000e+00  1.000000000000000e+00 -0.000000000000000e+00
+
+# Bond topology
+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  8  9
+9  1  9  10
+10  1  10  11
+11  1  11  12
+12  1  12  13
+13  1  14  15
+14  1  15  16
+15  1  16  17
+16  1  17  18
+17  1  18  19
+18  1  19  20
+19  1  20  21
+20  1  21  22
+21  1  22  23
+22  1  23  24
+23  1  24  25
+24  1  25  26
diff --git a/examples/USER/cgdna/examples/oxDNA2/unique_bp/generate_unique.py b/examples/USER/cgdna/examples/oxDNA2/unique_bp/generate_unique.py
new file mode 100644
index 0000000000..e5141bc47a
--- /dev/null
+++ b/examples/USER/cgdna/examples/oxDNA2/unique_bp/generate_unique.py
@@ -0,0 +1,828 @@
+#!/usr/bin/env python
+"""
+/* ----------------------------------------------------------------------
+   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 (University of Strathclyde, Glasgow)
+------------------------------------------------------------------------- */
+"""
+"""
+Creates unique base-pairings to avoid asymmetrical H-bonds.  
+
+Modified to create the bead wall setup.
+N_BEADS is the number of beads along one direction, the final system will have N_BEADS^2 beads in the wall. N_BEADS should be set to be odd number.
+"""
+
+
+#Define number of base-pairs per turn for B-form DNA
+N = 10.5
+#Define distance between the big bead and the centre of mass of the last base-pair
+BEAD_OFFSET = 2.0
+WALL_PARTICLE_SIZE = 2.0
+N_BEADS = 11
+
+#Number of unique base type groups (1-4) ACGT counts as one group
+N_BASE_TYPES = 20
+
+
+"""
+Import basic modules
+"""
+import sys, os, timeit
+
+from timeit import default_timer as timer
+start_time = timer()
+"""
+Try to import numpy; if failed, import a local version mynumpy 
+which needs to be provided
+"""
+try:
+    import numpy as np
+except:
+    print("numpy not found. Exiting.", file=sys.stderr)
+    sys.exit(1)
+
+"""
+Check that the required arguments (box offset and size in simulation units 
+and the sequence file were provided
+"""
+try:
+    box_offset = float(sys.argv[1])
+    box_length = float(sys.argv[2])
+    infile = sys.argv[3]
+    if len(sys.argv) == 4:
+        topo = 'strand'
+        lk = 0
+    elif len(sys.argv) == 5:
+        topo = 'strand'
+        lk = int(sys.argv[4])
+           
+except:
+    print("Usage: %s <%s> <%s> <%s> <%s> " % (sys.argv[0], \
+	"box offset", "box length", "file with sequences", "[Lk]"), file=sys.stderr)
+    sys.exit(1)
+box = np.array ([box_length, box_length, box_length])
+
+"""
+Try to open the file and fail gracefully if file cannot be opened
+"""
+try:
+    inp = open (infile, 'r')
+    inp.close()
+except:
+    print("Could not open file '%s' for reading. \
+					      Aborting." % infile, file=sys.stderr)
+    sys.exit(2)
+
+# return parts of a string
+def partition(s, d):
+    if d in s:
+        sp = s.split(d, 1)
+        return sp[0], d, sp[1]
+    else:
+        return s, "", ""
+
+"""
+Define the model constants
+"""
+# set model constants
+PI = np.pi
+POS_BASE =  0.4
+POS_BACK = -0.4
+EXCL_RC1 =  0.711879214356
+EXCL_RC2 =  0.335388426126
+EXCL_RC3 =  0.52329943261
+
+"""
+Define auxillary variables for the construction of a helix
+"""
+# center of the double strand
+COM_CENTRE_DS = POS_BASE + 0.2
+
+# ideal rise between two consecutive nucleotides on the
+# same strand which are to be base paired in a duplex
+BASE_BASE = 0.3897628551303122
+
+# cutoff distance for overlap check
+RC2 = 16
+
+# squares of the excluded volume distances for overlap check
+RC2_BACK = EXCL_RC1**2
+RC2_BASE = EXCL_RC2**2
+RC2_BACK_BASE = EXCL_RC3**2
+
+# enumeration to translate from letters to numbers and vice versa
+number_to_base = {1 : 'A', 2 : 'C', 3 : 'G', 4 : 'T'}
+base_to_number = {'A' : 1, 'a' : 1, 'C' : 2, 'c' : 2,
+                  'G' : 3, 'g' : 3, 'T' : 4, 't' : 4}
+
+# auxillary arrays
+positions = []
+a1s = []
+a3s = []
+quaternions = []
+
+newpositions = []
+newa1s = []
+newa3s = []
+
+basetype  = []
+strandnum = []
+
+bonds = []
+
+""" 
+Convert local body frame to quaternion DOF
+"""
+def exyz_to_quat (mya1, mya3):
+
+    mya2 = np.cross(mya3, mya1)
+    myquat = [1,0,0,0]
+
+    q0sq = 0.25 * (mya1[0] + mya2[1] + mya3[2] + 1.0)
+    q1sq = q0sq - 0.5 * (mya2[1] + mya3[2])
+    q2sq = q0sq - 0.5 * (mya1[0] + mya3[2])
+    q3sq = q0sq - 0.5 * (mya1[0] + mya2[1])
+
+    # some component must be greater than 1/4 since they sum to 1
+    # compute other components from it
+
+    if q0sq >= 0.25:
+        myquat[0] = np.sqrt(q0sq)
+        myquat[1] = (mya2[2] - mya3[1]) / (4.0*myquat[0])
+        myquat[2] = (mya3[0] - mya1[2]) / (4.0*myquat[0])
+        myquat[3] = (mya1[1] - mya2[0]) / (4.0*myquat[0])
+    elif q1sq >= 0.25:
+        myquat[1] = np.sqrt(q1sq)
+        myquat[0] = (mya2[2] - mya3[1]) / (4.0*myquat[1])
+        myquat[2] = (mya2[0] + mya1[1]) / (4.0*myquat[1])
+        myquat[3] = (mya1[2] + mya3[0]) / (4.0*myquat[1])
+    elif q2sq >= 0.25:
+        myquat[2] = np.sqrt(q2sq)
+        myquat[0] = (mya3[0] - mya1[2]) / (4.0*myquat[2])
+        myquat[1] = (mya2[0] + mya1[1]) / (4.0*myquat[2])
+        myquat[3] = (mya3[1] + mya2[2]) / (4.0*myquat[2])
+    elif q3sq >= 0.25:
+        myquat[3] = np.sqrt(q3sq)
+        myquat[0] = (mya1[1] - mya2[0]) / (4.0*myquat[3])
+        myquat[1] = (mya3[0] + mya1[2]) / (4.0*myquat[3])
+        myquat[2] = (mya3[1] + mya2[2]) / (4.0*myquat[3])
+
+    norm = 1.0/np.sqrt(myquat[0]*myquat[0] + myquat[1]*myquat[1] + \
+			  myquat[2]*myquat[2] + myquat[3]*myquat[3])
+    myquat[0] *= norm
+    myquat[1] *= norm
+    myquat[2] *= norm
+    myquat[3] *= norm
+
+    return np.array([myquat[0],myquat[1],myquat[2],myquat[3]])
+
+"""
+Adds a strand to the system by appending it to the array of previous strands
+"""
+def add_strands (mynewpositions, mynewa1s, mynewa3s):
+    overlap = False
+	
+    # This is a simple check for each of the particles where for previously 
+    # placed particles i we check whether it overlaps with any of the 
+    # newly created particles j
+
+    print("## Checking for overlaps", file=sys.stdout)
+
+    for i in range(len(positions)):
+
+        p = positions[i]
+        pa1 = a1s[i]
+
+        for j in range (len(mynewpositions)):
+
+            q = mynewpositions[j]
+            qa1 = mynewa1s[j]
+
+            # skip particles that are anyway too far away
+            dr = p - q
+            dr -= box * np.rint (dr / box)
+            if np.dot(dr, dr) > RC2:
+                continue
+
+            # base site and backbone site of the two particles
+            p_pos_back = p + pa1 * POS_BACK
+            p_pos_base = p + pa1 * POS_BASE
+            q_pos_back = q + qa1 * POS_BACK
+            q_pos_base = q + qa1 * POS_BASE
+
+            # check for no overlap between the two backbone sites
+            dr = p_pos_back - q_pos_back
+            dr -= box * np.rint (dr / box)
+            if np.dot(dr, dr) < RC2_BACK:
+                overlap = True
+
+            # check for no overlap between the two base sites
+            dr = p_pos_base -  q_pos_base
+            dr -= box * np.rint (dr / box)
+            if np.dot(dr, dr) < RC2_BASE:
+                overlap = True
+
+            # check for no overlap between backbone site of particle p 
+            # with base site of particle q
+            dr = p_pos_back - q_pos_base
+            dr -= box * np.rint (dr / box)
+            if np.dot(dr, dr) < RC2_BACK_BASE:
+                overlap = True
+
+            # check for no overlap between base site of particle p and 
+            # backbone site of particle q
+            dr = p_pos_base - q_pos_back
+            dr -= box * np.rint (dr / box)
+            if np.dot(dr, dr) < RC2_BACK_BASE:
+                overlap = True
+
+            # exit if there is an overlap
+            if overlap:
+                return False
+
+    # append to the existing list if no overlap is found
+    if not overlap:
+
+        for p in mynewpositions:
+            positions.append(p)
+        for p in mynewa1s:
+            a1s.append (p)
+        for p in mynewa3s:
+            a3s.append (p)
+	# calculate quaternion from local body frame and append
+        for ia in range(len(mynewpositions)):
+            mynewquaternions = exyz_to_quat(mynewa1s[ia],mynewa3s[ia])
+            quaternions.append(mynewquaternions)
+
+    return True
+
+"""
+Calculate angle of rotation site to site
+"""
+def get_angle(bp):
+   #n, minimal number of bases per turn
+   n = 10.5
+   found = False
+   while found == False:
+       turns = bp/n
+       diff = abs( turns - round(turns))
+       if diff < 0.03:
+           found = True
+           turns = round(turns)+lk
+           angle = (360*turns)/bp
+           angle = round (angle,2)
+           #angle =round( 360/n,2)
+       elif n > 11.5:
+           angle = 35.9
+           found = True
+       else:
+           n += 0.02
+   return angle
+
+
+def get_angle2(bp):
+    turns = bp/N + lk
+    angle = (360*turns)/bp
+    
+    return angle
+
+    
+
+"""
+Returns the rotation matrix defined by an axis and angle
+"""
+def get_rotation_matrix(axis, anglest, nbp=0):
+    # The argument anglest can be either an angle in radiants
+    # (accepted types are float, int or np.float64 or np.float64)
+    # or a tuple [angle, units] where angle is a number and
+    # units is a string. It tells the routine whether to use degrees,
+    # radiants (the default) or base pairs turns.
+    if not isinstance (anglest, (np.float64, np.float32, float, int)):
+        if len(anglest) > 1:
+            if anglest[1] in ["degrees", "deg", "o"]:
+                angle = (np.pi / 180.) * (anglest[0])
+            elif anglest[1] in ["bp"]:
+                if nbp == 0:
+                    angle = int(anglest[0]) * (np.pi / 180.) * (35.9)
+                else:    
+                    ang = get_angle2(nbp)
+                    angle = int(anglest[0]) * (np.pi / 180.) * (ang)
+            else:
+                angle = float(anglest[0])
+        else:
+            angle = float(anglest[0])
+    else:
+        angle = float(anglest) # in degrees (?)
+
+    axis = np.array(axis)
+    axis /= np.sqrt(np.dot(axis, axis))
+
+    ct = np.cos(angle)
+    st = np.sin(angle)
+    olc = 1. - ct
+    x, y, z = axis
+
+    return np.array([[olc*x*x+ct, olc*x*y-st*z, olc*x*z+st*y],
+                    [olc*x*y+st*z, olc*y*y+ct, olc*y*z-st*x],
+                    [olc*x*z-st*y, olc*y*z+st*x, olc*z*z+ct]])
+
+"""
+Generates the position and orientation vectors of a 
+(single or double) strand from a sequence string
+"""
+def generate_strand(bp, sequence=None, start_pos=np.array([0, 0, 0]), \
+	  dir=np.array([0, 0, 1]), perp=False, double=True, rot=0.):
+    # generate empty arrays
+    mynewpositions, mynewa1s, mynewa3s = [], [], []
+
+    # cast the provided start_pos array into a numpy array
+    start_pos = np.array(start_pos, dtype=float)
+
+    # overall direction of the helix
+    dir = np.array(dir, dtype=float)
+    #if sequence == None:
+    #    sequence = np.random.randint(1, 5, bp)
+
+    # the elseif here is most likely redundant 
+    #elif len(sequence) != bp:
+    #    n = bp - len(sequence)
+    #    sequence += np.random.randint(1, 5, n)
+    #    print("sequence is too short, adding %d random bases" % n, file=sys.stderr)
+
+    # normalize direction
+    dir_norm = np.sqrt(np.dot(dir,dir))
+    if dir_norm < 1e-10:
+        print("direction must be a valid vector,\
+                defaulting to (0, 0, 1)", file=sys.stderr)
+        dir = np.array([0, 0, 1])
+    else: dir /= dir_norm
+
+    # find a vector orthogonal to dir to act as helix direction,
+    # if not provided switch off random orientation
+    if perp is None or perp is False:
+        v1 = np.random.random_sample(3)
+        # comment in to suppress randomised base vector
+        v1 = [1,0,0]
+        v1 -= dir * (np.dot(dir, v1))
+        v1 /= np.sqrt(sum(v1*v1))
+    else:
+        v1 = perp;
+
+    # generate rotational matrix representing the overall rotation of the helix
+    R0 = get_rotation_matrix(dir, rot)
+	    
+    # rotation matrix corresponding to one step along the helix
+    R = get_rotation_matrix(dir, [1, "bp"],bp)
+
+    # set the vector a1 (backbone to base) to v1 
+    a1 = v1
+    
+    # apply the global rotation to a1 
+    a1 = np.dot(R0, a1)
+    
+    # set the position of the fist backbone site to start_pos
+    rb = np.array(start_pos)
+	    
+    # set a3 to the direction of the helix
+    a3 = dir
+
+    for i in range(bp):
+    # work out the position of the centre of mass of the nucleotide
+        rcom = rb - COM_CENTRE_DS * a1
+
+        # append to newpositions
+        mynewpositions.append(rcom)
+        mynewa1s.append(a1)
+        mynewa3s.append(a3)
+
+        # if we are not at the end of the helix, we work out a1 and rb for the 
+        # next nucleotide along the helix
+        if i != bp - 1:
+            a1 = np.dot(R, a1)
+            rb += a3 * BASE_BASE
+
+    # if we are working on a double strand, we do a cycle similar 
+    # to the previous one but backwards
+    if double == True:
+        a1 = -a1
+        a3 = -dir
+        R = R.transpose()
+        for i in range(bp):
+            rcom = rb - COM_CENTRE_DS * a1
+            mynewpositions.append (rcom)
+            mynewa1s.append (a1)
+            mynewa3s.append (a3)
+            a1 = np.dot(R, a1)
+            rb += a3 * BASE_BASE
+    
+  
+    #Calculate the positions of the bead wall
+    
+    last_base1 = mynewpositions[int( len(mynewpositions)/2 - 1) ]
+    last_base2 = mynewpositions[int( len(mynewpositions)/2) ]
+    mid_point =  (last_base1 + last_base2) / 2 
+    
+    NN = N_BEADS**2
+    p1 = [mid_point[0] - (N_BEADS-1)*WALL_PARTICLE_SIZE, mid_point[1] - (N_BEADS-1)*WALL_PARTICLE_SIZE, mid_point[2] + BEAD_OFFSET ]
+    for i in range(N_BEADS):
+        for j in range(N_BEADS):
+            position = [ p1[0] + 2*i*WALL_PARTICLE_SIZE, p1[1] + 2*j*WALL_PARTICLE_SIZE, p1[2]]
+            mynewa1s.append([1,0,0])
+            mynewa3s.append([1,0,0])
+            mynewpositions.append(position)
+
+    assert (len (mynewpositions) > 0)
+
+    return [mynewpositions, mynewa1s, mynewa3s]
+
+
+
+"""
+Main function for this script.
+Reads a text file with the following format:
+- Each line contains the sequence for a single strand (A,C,G,T)
+- Lines beginning with the keyword 'DOUBLE' produce double-stranded DNA
+
+Ex: Two ssDNA (single stranded DNA)
+ATATATA
+GCGCGCG
+
+Ex: Two strands, one double stranded, the other single stranded.
+DOUBLE AGGGCT
+CCTGTA
+
+"""
+
+def read_strands(filename):
+    try:
+        infile = open (filename)
+    except:
+        print("Could not open file '%s'. Aborting." % filename, file=sys.stderr)
+        sys.exit(2)
+
+    # This block works out the number of nucleotides and strands by reading 
+    # the number of non-empty lines in the input file and the number of letters,
+    # taking the possible DOUBLE keyword into account.
+    nstrands, nnucl, nbonds = 0, 0, 0
+    lines = infile.readlines()
+    for line in lines:
+        line = line.upper().strip()
+        if len(line) == 0:
+            continue
+        if line[:6] == 'DOUBLE':
+            line = line.split()[1]
+            length = len(line)
+            print("## Found duplex of %i base pairs" % length, file=sys.stdout)
+            nnucl += 2*length
+            nstrands += 2
+            nbonds+= 2*length
+            
+        else:
+            line = line.split()[0]
+            length = len(line)
+            print("## Found single strand of %i bases" % length, file=sys.stdout)
+            nnucl += length
+            nstrands += 1
+            if topo == 'ring':
+                nbonds =+ length
+            else:
+                nbonds += length+1
+    # rewind the sequence input file
+    infile.seek(0)
+
+    print("## nstrands, nnucl = ", nstrands, nnucl, file=sys.stdout)
+
+    # generate the data file in LAMMPS format
+    try:
+        out = open ("data.oxdna", "w")
+    except:
+        print("Could not open data file for writing. Aborting.", file=sys.stderr)
+        sys.exit(2)
+	
+    lines = infile.readlines()
+    nlines = len(lines)
+    i = 1
+    myns = 0
+    noffset = 1
+
+    for line in lines:
+        line = line.upper().strip()
+
+        # skip empty lines
+        if len(line) == 0: 
+            i += 1
+            continue
+
+	# block for duplexes: last argument of the generate function 
+	# is set to 'True'
+        if line[:6] == 'DOUBLE':
+            line = line.split()[1]
+            length = len(line)
+            seq = [(base_to_number[x]) for x in line]
+            seq = np.array(seq,dtype=int)
+            n_a, n_c, n_g, n_t = 0, 0, 0, 0
+            for s in range(seq.size):
+                if seq[s] == 1:
+                    n_a += 1
+                elif seq[s] == 2:
+                    n_c += 1
+                elif seq[s] ==3:
+                    n_g += 1
+                elif seq[s] == 4:
+                    n_t += 1
+            smallest_n_bases = n_c
+            if n_a < n_c:
+                smallest_n_bases = n_a
+            if smallest_n_bases > n_t:
+                smallest_n_bases = n_t
+            if smallest_n_bases > n_g:
+                smallest_n_bases = n_g
+             
+            if smallest_n_bases < N_BASE_TYPES:
+                print('## Not enough occurances of base types in the sequence for ' + str(N_BASE_TYPES))
+                print('## unique base types, switching to ' + str(smallest_n_bases) + ' unique types')
+            else:
+                smallest_n_bases = N_BASE_TYPES
+           
+            a, c, g, t = -3, -2, -1, 0
+            for s in range(seq.size):
+                if seq[s] == 1: 
+                    if a < (smallest_n_bases*4-3):
+                        a += 4
+                    else:
+                        a = 1
+                    seq[s] = a
+                
+                elif seq[s] == 2: 
+                    if c < (smallest_n_bases*4-2):
+                        c += 4
+                    else:
+                        c = 2
+                    seq[s] = c
+
+                elif seq[s] == 3:
+                    if g < (smallest_n_bases*4-1):
+                        g += 4
+                    else:
+                        g = 3
+                    seq[s] = g
+                elif seq[s] == 4:
+                    if t < (smallest_n_bases*4):
+                        t += 4
+                    else:
+                        t = 4
+                    seq[s] = t
+                
+
+
+            myns += 1
+
+            for b in range(length):
+                basetype.append(seq[b])
+                strandnum.append(myns)
+
+            for b in range(length-1):
+                bondpair = [noffset + b, noffset + b + 1]
+                bonds.append(bondpair)
+     
+            
+            noffset += length 
+
+            # create the sequence of the second strand as made of 
+            # complementary bases
+            #seq2 = [5-s for s in seq]
+            seq2 = seq
+            for s in range(seq2.size):
+                if seq2[s]%4 == 1:
+                    seq2[s] += 3
+                elif seq2[s]%4 == 2:
+                    seq2[s] += 1
+                elif seq2[s]%4 == 3:
+                    seq2[s] -= 1
+                elif seq2[s]%4 == 0:
+                    seq2[s] -= 3    
+
+            #seq2.reverse()
+
+            myns += 1
+
+            for b in range(length):
+                basetype.append(seq2[b])
+                strandnum.append(myns)
+
+            for b in range(length-1):
+                bondpair = [noffset + b, noffset + b + 1]
+                bonds.append(bondpair)
+
+                     
+            #create wall bead types
+            bead_type = 4*smallest_n_bases + 1
+            for i in range(N_BEADS**2):
+                
+                basetype.append(bead_type)
+                basetype.append(bead_type)
+                strandnum.append(bead_type)
+                strandnum.append(bead_type)
+            #bonds.append([length, noffset + length])
+            #bonds.append([length+1, noffset + length])
+
+            noffset += length
+
+            print("## Created duplex of %i bases" % (2*length), file=sys.stdout)
+
+            # generate random position of the first nucleotide
+            com = box_offset + np.random.random_sample(3) * box
+            # comment out to randomise
+            com = [0,0,0]
+
+            # generate the random direction of the helix 
+            axis = np.random.random_sample(3)
+            # comment out to randomise
+            axis = [0,0,1]
+            axis /= np.sqrt(np.dot(axis, axis))
+
+            # use the generate function defined above to create 
+            # the position and orientation vector of the strand
+            if topo == 'ring':
+                newpositions, newa1s, newa3s = generate_ring(len(line), \
+                      sequence=seq, dir=axis, start_pos=com, double=True)
+            else:
+                newpositions, newa1s, newa3s = generate_strand(len(line), \
+                    sequence=seq, dir=axis, start_pos=com, double=True)
+
+            # generate a new position for the strand until it does not overlap
+            # with anything already present
+            start = timer()
+            while not add_strands(newpositions, newa1s, newa3s):
+                com = box_offset + np.random.random_sample(3) * box
+                axis = np.random.random_sample(3)
+                axis /= np.sqrt(np.dot(axis, axis))
+                if topo == 'ring':
+                    newpositions, newa1s, newa3s = generate_ring(len(line), \
+                        sequence=seq, dir=axis, start_pos=com, double=True)
+                else:
+                    newpositions, newa1s, newa3s = generate_strand(len(line), \
+                        sequence=seq, dir=axis, start_pos=com, double=True)
+                print("## Trying %i" % i, file=sys.stdout)
+            end = timer()
+            print("## Added duplex of %i bases (line %i/%i) in %.2fs, now at %i/%i" % \
+                          (2*length, i, nlines, end-start, len(positions), nnucl), file=sys.stdout)
+
+	# block for single strands: last argument of the generate function 
+	# is set to 'False'
+        else:
+            length = len(line)
+            seq = [(base_to_number[x]) for x in line]
+
+            myns += 1
+            for b in range(length):
+                basetype.append(seq[b])
+                strandnum.append(myns)
+
+            for b in range(length-1):
+                bondpair = [noffset + b, noffset + b + 1]
+                bonds.append(bondpair)
+            if topo == 'ring':    
+                bondpair = [noffset, noffset + length-1]
+                bonds.append(bondpair)
+            noffset += length
+            
+
+	    # generate random position of the first nucleotide
+            com = box_offset + np.random.random_sample(3) * box
+            # comment out to randomise
+            com = [-30,0,0]
+
+            # generate the random direction of the helix 
+            axis = np.random.random_sample(3)
+            # comment out to randomise
+            axis = [0,0,1]
+            axis /= np.sqrt(np.dot(axis, axis))
+
+            print("## Created single strand of %i bases" % length, file=sys.stdout)
+            if topo == 'ring':    
+                newpositions, newa1s, newa3s = generate_ring(length, \
+		              sequence=seq, dir=axis, start_pos=com, double=False)
+            else:
+                newpositions, newa1s, newa3s = generate_strand(length, \
+		               sequence=seq, dir=axis, start_pos=com, double=False)
+            start = timer()
+            while not add_strands(newpositions, newa1s, newa3s):
+                com = box_offset + np.random.random_sample(3) * box
+                axis = np.random.random_sample(3)
+                axis /= np.sqrt(np.dot(axis, axis))
+                if topo == 'ring':
+                    newpositions, newa1s, newa3s = generate_ring(length, \
+                          sequence=seq, dir=axis, start_pos=com, double=False)
+
+                else:    
+                    newpositions, newa1s, newa3s = generate_strand(length, \
+                          sequence=seq, dir=axis, start_pos=com, double=False)
+                print("## Trying  %i" % (i), file=sys.stdout)
+            end = timer()
+            print("## Added single strand of %i bases (line %i/%i) in %.2fs, now at %i/%i" % \
+				      (length, i, nlines, end-start,len(positions), nnucl), file=sys.stdout)
+
+        i += 1
+
+    # sanity check
+    #if not len(positions) == nnucl:
+    #    print(len(positions), nnucl)
+    #    raise AssertionError
+    nnucl = nnucl + (N_BEADS**2)
+    nbonds -= 4
+	
+    out.write('# LAMMPS data file\n')
+    out.write('%d atoms\n' % nnucl)
+    out.write('%d ellipsoids\n' % nnucl)
+    out.write('%d bonds\n' % nbonds)
+    out.write('\n')
+    out.write('%d atom types\n' %bead_type )
+    out.write('1 bond types\n')
+    out.write('\n')
+    out.write('# System size\n')
+    out.write('%f %f xlo xhi\n' % (box_offset,box_offset+box_length))
+    out.write('%f %f ylo yhi\n' % (box_offset,box_offset+box_length))
+    out.write('%f %f zlo zhi\n' % (0,box_length))
+
+    #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('5 3.1575\n')
+
+    # for each nucleotide print a line under the headers
+    # Atoms, Velocities, Ellipsoids and Bonds
+    out.write('\n')
+    out.write(\
+      '# Atom-ID, type, position, molecule-ID, ellipsoid flag, density\n')
+    out.write('Atoms\n')
+    out.write('\n')
+    
+    for i in range(nnucl):
+        out.write('%d %d %22.15le %22.15le %22.15le %d 1 1\n' \
+                  % (i+1, basetype[i], \
+                     positions[i][0], positions[i][1], positions[i][2], \
+                     strandnum[i]))
+
+    out.write('\n')
+    out.write('# Atom-ID, translational, rotational velocity\n')
+    out.write('Velocities\n')
+    out.write('\n')
+
+    for i in range(nnucl):
+        out.write("%d %22.15le %22.15le %22.15le %22.15le %22.15le %22.15le\n" \
+                  % (i+1,0.0,0.0,0.0,0.0,0.0,0.0))
+
+    out.write('\n')
+    out.write('# Atom-ID, shape, quaternion\n')
+    out.write('Ellipsoids\n')
+    out.write('\n')
+
+    for i in range(nnucl):
+        out.write(\
+    "%d %22.15le %22.15le %22.15le %22.15le %22.15le %22.15le %22.15le\n"  \
+      % (i+1,1.1739845031423408,1.1739845031423408,1.1739845031423408, \
+	quaternions[i][0],quaternions[i][1], quaternions[i][2],quaternions[i][3]))
+ 
+    out.write('\n')
+    out.write('# Bond topology\n')
+    out.write('Bonds\n')
+    out.write('\n')
+
+    for i in range(nbonds):
+        if i < nbonds-2:
+             out.write("%d  %d  %d  %d\n" % (i+1,1,bonds[i][0],bonds[i][1]))
+        #else:
+        
+        #out.write("%d  %d  %d  %d\n" % (i+1,2,bonds[i][0],bonds[i][1]))
+
+    out.close()
+
+    print("## Wrote data to 'data.oxdna'", file=sys.stdout)
+    print("## DONE", file=sys.stdout)
+
+# call the above main() function, which executes the program
+read_strands (infile)
+
+end_time=timer()
+runtime = end_time-start_time
+hours = runtime/3600
+minutes = (runtime-np.rint(hours)*3600)/60
+seconds = (runtime-np.rint(hours)*3600-np.rint(minutes)*60)%60
+print("## Total runtime %ih:%im:%.2fs" % (hours,minutes,seconds), file=sys.stdout)
diff --git a/examples/USER/cgdna/examples/oxDNA2/unique_bp/in.duplex4.4type b/examples/USER/cgdna/examples/oxDNA2/unique_bp/in.duplex4.4type
new file mode 100644
index 0000000000..0570cf042d
--- /dev/null
+++ b/examples/USER/cgdna/examples/oxDNA2/unique_bp/in.duplex4.4type
@@ -0,0 +1,94 @@
+variable number	equal 1
+variable ofreq	equal 10000
+variable efreq	equal 10000
+
+variable ntype equal 4
+
+variable T equal 0.1
+
+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 2.0 bin 
+neigh_modify every 10 delay 0 check yes
+
+read_data data.duplex4.4type
+mass * 3.1575
+
+group all type 1 4
+
+# oxDNA bond interactions - FENE backbone
+bond_style oxdna2/fene 
+bond_coeff * 2.0 0.25 0.7564
+
+# oxDNA pair interactions
+pair_style hybrid/overlay oxdna2/excv oxdna2/stk oxdna2/hbond oxdna2/xstk oxdna2/coaxstk oxdna2/dh 
+pair_coeff * * oxdna2/excv    2.0 0.7 0.675 2.0 0.515 0.5 2.0 0.33 0.32
+pair_coeff * * oxdna2/stk     seqdep ${T} 1.3523 2.6717 6.0 0.4 0.9 0.32 0.75 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 * * oxdna2/hbond   seqdep 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
+
+label loop
+variable base loop ${ntype} 
+  variable basemod equal ${base}%4
+  if "${basemod} == 1" then &
+    "variable comp equal ${base}+3" &
+    "pair_coeff ${base} ${comp} oxdna2/hbond   seqdep 1.0678 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"
+  if "${basemod} == 2" then &
+    "variable comp equal ${base}+1" &
+    "pair_coeff ${base} ${comp} oxdna2/hbond   seqdep 1.0678 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"
+next base
+jump in.duplex4.4type loop
+
+pair_coeff * * oxdna2/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 * * oxdna2/coaxstk 58.5 0.4 0.6 0.22 0.58 2.0 2.891592653589793 0.65 1.3 0 0.8 0.9 0 0.95 0.9 0 0.95 40.0 3.116592653589793
+pair_coeff * * oxdna2/dh      ${T} 0.2 0.815
+
+# Langevin dynamics
+fix 1 all nve/asphere
+fix 2 all langevin ${T} ${T} 25.0 457145 angmom 10
+
+timestep 1e-4
+
+#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 mol type x y z ix iy iz c_quat[1] c_quat[2] c_quat[3] c_quat[4]  vx vy vz 
+dump_modify out sort id
+dump_modify out format line "%d %d %d  %13.6le %13.6le %13.6le %d %d %d  %13.6le %13.6le %13.6le  %13.6le %13.6le %13.6le  %13.6le "
+
+#restart 10000 config0_restart config1_restart
+
+run 100000
+
+#write_restart config.${number}.*
+
+
diff --git a/examples/USER/cgdna/examples/oxDNA2/unique_bp/in.duplex4.8type b/examples/USER/cgdna/examples/oxDNA2/unique_bp/in.duplex4.8type
new file mode 100644
index 0000000000..dae853a113
--- /dev/null
+++ b/examples/USER/cgdna/examples/oxDNA2/unique_bp/in.duplex4.8type
@@ -0,0 +1,94 @@
+variable number	equal 1
+variable ofreq	equal 10000
+variable efreq	equal 10000
+
+variable ntype equal 8
+
+variable T equal 0.1
+
+units lj
+
+dimension 3
+
+newton on
+
+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 10 delay 0 check yes
+
+read_data data.duplex4.8type
+mass * 3.1575
+
+group all type 1 8
+
+# oxDNA bond interactions - FENE backbone
+bond_style oxdna2/fene 
+bond_coeff * 2.0 0.25 0.7564
+
+# oxDNA pair interactions
+pair_style hybrid/overlay oxdna2/excv oxdna2/stk oxdna2/hbond oxdna2/xstk oxdna2/coaxstk oxdna2/dh 
+pair_coeff * * oxdna2/excv    2.0 0.7 0.675 2.0 0.515 0.5 2.0 0.33 0.32
+pair_coeff * * oxdna2/stk     seqdep ${T} 1.3523 2.6717 6.0 0.4 0.9 0.32 0.75 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 * * oxdna2/hbond   seqdep 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
+
+label loop
+variable base loop ${ntype} 
+  variable basemod equal ${base}%4
+  if "${basemod} == 1" then &
+    "variable comp equal ${base}+3" &
+    "pair_coeff ${base} ${comp} oxdna2/hbond   seqdep 1.0678 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"
+  if "${basemod} == 2" then &
+    "variable comp equal ${base}+1" &
+    "pair_coeff ${base} ${comp} oxdna2/hbond   seqdep 1.0678 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"
+next base
+jump in.duplex4.8type loop
+
+pair_coeff * * oxdna2/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 * * oxdna2/coaxstk 58.5 0.4 0.6 0.22 0.58 2.0 2.891592653589793 0.65 1.3 0 0.8 0.9 0 0.95 0.9 0 0.95 40.0 3.116592653589793
+pair_coeff * * oxdna2/dh      ${T} 0.2 0.815
+
+# Langevin dynamics
+fix 1 all nve/asphere
+fix 2 all langevin ${T} ${T} 25.0 457145 angmom 10
+
+timestep 1e-4
+
+#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 mol type x y z ix iy iz c_quat[1] c_quat[2] c_quat[3] c_quat[4]  vx vy vz 
+dump_modify out sort id
+dump_modify out format line "%d %d %d  %13.6le %13.6le %13.6le %d %d %d  %13.6le %13.6le %13.6le  %13.6le %13.6le %13.6le  %13.6le "
+
+#restart 10000 config0_restart config1_restart
+
+run 100000
+
+#write_restart config.${number}.*
+
+
diff --git a/examples/USER/cgdna/examples/oxDNA2/unique_bp/log.07Aug19.duplex4.4type.g++1 b/examples/USER/cgdna/examples/oxDNA2/unique_bp/log.07Aug19.duplex4.4type.g++1
new file mode 100644
index 0000000000..d0195dd8e8
--- /dev/null
+++ b/examples/USER/cgdna/examples/oxDNA2/unique_bp/log.07Aug19.duplex4.4type.g++1
@@ -0,0 +1,232 @@
+LAMMPS (7 Aug 2019)
+variable number	equal 1
+variable ofreq	equal 10000
+variable efreq	equal 10000
+
+variable ntype equal 4
+
+variable T equal 0.1
+
+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 2.0 bin
+neigh_modify every 10 delay 0 check yes
+
+read_data data.duplex4.4type
+  orthogonal box = (-20 -20 -20) to (20 20 20)
+  1 by 1 by 1 MPI processor grid
+  reading atoms ...
+  26 atoms
+  reading velocities ...
+  26 velocities
+  26 ellipsoids
+  scanning bonds ...
+  2 = max bonds/atom
+  reading bonds ...
+  24 bonds
+  2 = max # of 1-2 neighbors
+  2 = max # of 1-3 neighbors
+  4 = max # of 1-4 neighbors
+  6 = max # of special neighbors
+  special bonds CPU = 6.2e-05 secs
+  read_data CPU = 0.001124 secs
+mass * 3.1575
+
+group all type 1 4
+26 atoms in group all
+
+# oxDNA bond interactions - FENE backbone
+bond_style oxdna2/fene
+bond_coeff * 2.0 0.25 0.7564
+
+# oxDNA pair interactions
+pair_style hybrid/overlay oxdna2/excv oxdna2/stk oxdna2/hbond oxdna2/xstk oxdna2/coaxstk oxdna2/dh
+pair_coeff * * oxdna2/excv    2.0 0.7 0.675 2.0 0.515 0.5 2.0 0.33 0.32
+pair_coeff * * oxdna2/stk     seqdep ${T} 1.3523 2.6717 6.0 0.4 0.9 0.32 0.75 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 * * oxdna2/stk     seqdep 0.1 1.3523 2.6717 6.0 0.4 0.9 0.32 0.75 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 * * oxdna2/hbond   seqdep 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
+
+label loop
+variable base loop ${ntype}
+variable base loop 4
+  variable basemod equal ${base}%4
+  variable basemod equal 1%4
+  if "${basemod} == 1" then     "variable comp equal ${base}+3"     "pair_coeff ${base} ${comp} oxdna2/hbond   seqdep 1.0678 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"
+variable comp equal ${base}+3
+variable comp equal 1+3
+pair_coeff ${base} ${comp} oxdna2/hbond   seqdep 1.0678 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 ${comp} oxdna2/hbond   seqdep 1.0678 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 oxdna2/hbond   seqdep 1.0678 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
+  if "${basemod} == 2" then     "variable comp equal ${base}+1"     "pair_coeff ${base} ${comp} oxdna2/hbond   seqdep 1.0678 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"
+next base
+jump in.duplex4.4type loop
+variable base loop ${ntype}
+variable base loop 4
+  variable basemod equal ${base}%4
+  variable basemod equal 2%4
+  if "${basemod} == 1" then     "variable comp equal ${base}+3"     "pair_coeff ${base} ${comp} oxdna2/hbond   seqdep 1.0678 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"
+  if "${basemod} == 2" then     "variable comp equal ${base}+1"     "pair_coeff ${base} ${comp} oxdna2/hbond   seqdep 1.0678 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"
+variable comp equal ${base}+1
+variable comp equal 2+1
+pair_coeff ${base} ${comp} oxdna2/hbond   seqdep 1.0678 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 ${comp} oxdna2/hbond   seqdep 1.0678 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 oxdna2/hbond   seqdep 1.0678 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
+next base
+jump in.duplex4.4type loop
+variable base loop ${ntype}
+variable base loop 4
+  variable basemod equal ${base}%4
+  variable basemod equal 3%4
+  if "${basemod} == 1" then     "variable comp equal ${base}+3"     "pair_coeff ${base} ${comp} oxdna2/hbond   seqdep 1.0678 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"
+  if "${basemod} == 2" then     "variable comp equal ${base}+1"     "pair_coeff ${base} ${comp} oxdna2/hbond   seqdep 1.0678 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"
+next base
+jump in.duplex4.4type loop
+variable base loop ${ntype}
+variable base loop 4
+  variable basemod equal ${base}%4
+  variable basemod equal 4%4
+  if "${basemod} == 1" then     "variable comp equal ${base}+3"     "pair_coeff ${base} ${comp} oxdna2/hbond   seqdep 1.0678 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"
+  if "${basemod} == 2" then     "variable comp equal ${base}+1"     "pair_coeff ${base} ${comp} oxdna2/hbond   seqdep 1.0678 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"
+next base
+jump in.duplex4.4type loop
+
+pair_coeff * * oxdna2/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 * * oxdna2/coaxstk 58.5 0.4 0.6 0.22 0.58 2.0 2.891592653589793 0.65 1.3 0 0.8 0.9 0 0.95 0.9 0 0.95 40.0 3.116592653589793
+pair_coeff * * oxdna2/dh      ${T} 0.2 0.815
+pair_coeff * * oxdna2/dh      0.1 0.2 0.815
+
+# Langevin dynamics
+fix 1 all nve/asphere
+fix 2 all langevin ${T} ${T} 25.0 457145 angmom 10
+fix 2 all langevin 0.1 ${T} 25.0 457145 angmom 10
+fix 2 all langevin 0.1 0.1 25.0 457145 angmom 10
+
+timestep 1e-4
+
+#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
+fix 5 all print 10000 "$(step)  ekin = ${ekin} |  erot = ${erot} | epot = ${epot} | etot = ${etot}" screen yes
+
+dump out all custom ${ofreq} out.${number}.txt id mol type x y z ix iy iz c_quat[1] c_quat[2] c_quat[3] c_quat[4]  vx vy vz
+dump out all custom 10000 out.${number}.txt id mol type x y z ix iy iz c_quat[1] c_quat[2] c_quat[3] c_quat[4]  vx vy vz
+dump out all custom 10000 out.1.txt id mol type x y z ix iy iz c_quat[1] c_quat[2] c_quat[3] c_quat[4]  vx vy vz
+dump_modify out sort id
+dump_modify out format line "%d %d %d  %13.6le %13.6le %13.6le %d %d %d  %13.6le %13.6le %13.6le  %13.6le %13.6le %13.6le  %13.6le "
+
+#restart 10000 config0_restart config1_restart
+
+run 100000
+Neighbor list info ...
+  update every 10 steps, delay 0 steps, check yes
+  max neighbors/atom: 2000, page size: 100000
+  master list distance cutoff = 5.63899
+  ghost atom cutoff = 5.63899
+  binsize = 2.81949, bins = 15 15 15
+  6 neighbor lists, perpetual/occasional/extra = 6 0 0
+  (1) pair oxdna2/excv, perpetual
+      attributes: half, newton off
+      pair build: half/bin/newtoff
+      stencil: half/bin/3d/newtoff
+      bin: standard
+  (2) pair oxdna2/stk, perpetual, copy from (1)
+      attributes: half, newton off
+      pair build: copy
+      stencil: none
+      bin: none
+  (3) pair oxdna2/hbond, perpetual, copy from (1)
+      attributes: half, newton off
+      pair build: copy
+      stencil: none
+      bin: none
+  (4) pair oxdna2/xstk, perpetual, copy from (1)
+      attributes: half, newton off
+      pair build: copy
+      stencil: none
+      bin: none
+  (5) pair oxdna2/coaxstk, perpetual, copy from (1)
+      attributes: half, newton off
+      pair build: copy
+      stencil: none
+      bin: none
+  (6) pair oxdna2/dh, perpetual, copy from (1)
+      attributes: half, newton off
+      pair build: copy
+      stencil: none
+      bin: none
+Per MPI rank memory allocation (min/avg/max) = 3.89 | 3.89 | 3.89 Mbytes
+Step Temp E_pair E_mol TotEng Press 
+       0            0   -1.6384018  0.037304147   -1.6010976 6.7769766e-05 
+10000  ekin = 1.01296379553452 |  erot = 1.38476360245491 | epot = -41.9785360036035 | etot = -39.5808086056141
+20000  ekin = 1.75607695499755 |  erot = 1.86620102096623 | epot = -41.9558496477913 | etot = -38.3335716718275
+30000  ekin = 2.25878673688255 |  erot = 2.60035907744906 | epot = -41.2869827431736 | etot = -36.427836928842
+40000  ekin = 2.11105717434601 |  erot = 3.16611217122795 | epot = -40.2762731426449 | etot = -34.999103797071
+50000  ekin = 3.05379083545185 |  erot = 2.45050275456902 | epot = -40.3032957287999 | etot = -34.799002138779
+60000  ekin = 3.05154113920288 |  erot = 2.52539576708931 | epot = -39.6800099891302 | etot = -34.103073082838
+70000  ekin = 3.23212209366503 |  erot = 3.28914763888578 | epot = -40.1240667487288 | etot = -33.602797016178
+80000  ekin = 2.82623224207568 |  erot = 2.91292948677805 | epot = -39.0983962904507 | etot = -33.3592345615969
+90000  ekin = 3.05995314905566 |  erot = 3.5429982411034 | epot = -39.1646070343971 | etot = -32.561655644238
+100000  ekin = 3.46139546969836 |  erot = 4.11704803295073 | epot = -38.5124679837256 | etot = -30.9340244810765
+  100000  0.092303879   -1.5243833  0.043134544   -1.3481182 6.862374e-06 
+Loop time of 7.99505 on 1 procs for 100000 steps with 26 atoms
+
+Performance: 108066.893 tau/day, 12507.742 timesteps/s
+99.8% CPU use with 1 MPI tasks x no OpenMP threads
+
+MPI task timing breakdown:
+Section |  min time  |  avg time  |  max time  |%varavg| %total
+---------------------------------------------------------------
+Pair    | 7.2072     | 7.2072     | 7.2072     |   0.0 | 90.15
+Bond    | 0.14292    | 0.14292    | 0.14292    |   0.0 |  1.79
+Neigh   | 1.9e-05    | 1.9e-05    | 1.9e-05    |   0.0 |  0.00
+Comm    | 0.015676   | 0.015676   | 0.015676   |   0.0 |  0.20
+Output  | 0.001372   | 0.001372   | 0.001372   |   0.0 |  0.02
+Modify  | 0.61071    | 0.61071    | 0.61071    |   0.0 |  7.64
+Other   |            | 0.01714    |            |       |  0.21
+
+Nlocal:    26 ave 26 max 26 min
+Histogram: 1 0 0 0 0 0 0 0 0 0
+Nghost:    0 ave 0 max 0 min
+Histogram: 1 0 0 0 0 0 0 0 0 0
+Neighs:    325 ave 325 max 325 min
+Histogram: 1 0 0 0 0 0 0 0 0 0
+
+Total # of neighbors = 325
+Ave neighs/atom = 12.5
+Ave special neighs/atom = 5.07692
+Neighbor list builds = 1
+Dangerous builds = 0
+
+#write_restart config.${number}.*
+
+
+Total wall time: 0:00:07
diff --git a/examples/USER/cgdna/examples/oxDNA2/unique_bp/log.07Aug19.duplex4.4type.g++4 b/examples/USER/cgdna/examples/oxDNA2/unique_bp/log.07Aug19.duplex4.4type.g++4
new file mode 100644
index 0000000000..a5647ccd77
--- /dev/null
+++ b/examples/USER/cgdna/examples/oxDNA2/unique_bp/log.07Aug19.duplex4.4type.g++4
@@ -0,0 +1,232 @@
+LAMMPS (7 Aug 2019)
+variable number	equal 1
+variable ofreq	equal 10000
+variable efreq	equal 10000
+
+variable ntype equal 4
+
+variable T equal 0.1
+
+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 2.0 bin
+neigh_modify every 10 delay 0 check yes
+
+read_data data.duplex4.4type
+  orthogonal box = (-20 -20 -20) to (20 20 20)
+  1 by 2 by 2 MPI processor grid
+  reading atoms ...
+  26 atoms
+  reading velocities ...
+  26 velocities
+  26 ellipsoids
+  scanning bonds ...
+  2 = max bonds/atom
+  reading bonds ...
+  24 bonds
+  2 = max # of 1-2 neighbors
+  2 = max # of 1-3 neighbors
+  4 = max # of 1-4 neighbors
+  6 = max # of special neighbors
+  special bonds CPU = 0.000173 secs
+  read_data CPU = 0.00242 secs
+mass * 3.1575
+
+group all type 1 4
+26 atoms in group all
+
+# oxDNA bond interactions - FENE backbone
+bond_style oxdna2/fene
+bond_coeff * 2.0 0.25 0.7564
+
+# oxDNA pair interactions
+pair_style hybrid/overlay oxdna2/excv oxdna2/stk oxdna2/hbond oxdna2/xstk oxdna2/coaxstk oxdna2/dh
+pair_coeff * * oxdna2/excv    2.0 0.7 0.675 2.0 0.515 0.5 2.0 0.33 0.32
+pair_coeff * * oxdna2/stk     seqdep ${T} 1.3523 2.6717 6.0 0.4 0.9 0.32 0.75 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 * * oxdna2/stk     seqdep 0.1 1.3523 2.6717 6.0 0.4 0.9 0.32 0.75 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 * * oxdna2/hbond   seqdep 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
+
+label loop
+variable base loop ${ntype}
+variable base loop 4
+  variable basemod equal ${base}%4
+  variable basemod equal 1%4
+  if "${basemod} == 1" then     "variable comp equal ${base}+3"     "pair_coeff ${base} ${comp} oxdna2/hbond   seqdep 1.0678 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"
+variable comp equal ${base}+3
+variable comp equal 1+3
+pair_coeff ${base} ${comp} oxdna2/hbond   seqdep 1.0678 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 ${comp} oxdna2/hbond   seqdep 1.0678 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 oxdna2/hbond   seqdep 1.0678 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
+  if "${basemod} == 2" then     "variable comp equal ${base}+1"     "pair_coeff ${base} ${comp} oxdna2/hbond   seqdep 1.0678 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"
+next base
+jump in.duplex4.4type loop
+variable base loop ${ntype}
+variable base loop 4
+  variable basemod equal ${base}%4
+  variable basemod equal 2%4
+  if "${basemod} == 1" then     "variable comp equal ${base}+3"     "pair_coeff ${base} ${comp} oxdna2/hbond   seqdep 1.0678 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"
+  if "${basemod} == 2" then     "variable comp equal ${base}+1"     "pair_coeff ${base} ${comp} oxdna2/hbond   seqdep 1.0678 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"
+variable comp equal ${base}+1
+variable comp equal 2+1
+pair_coeff ${base} ${comp} oxdna2/hbond   seqdep 1.0678 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 ${comp} oxdna2/hbond   seqdep 1.0678 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 oxdna2/hbond   seqdep 1.0678 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
+next base
+jump in.duplex4.4type loop
+variable base loop ${ntype}
+variable base loop 4
+  variable basemod equal ${base}%4
+  variable basemod equal 3%4
+  if "${basemod} == 1" then     "variable comp equal ${base}+3"     "pair_coeff ${base} ${comp} oxdna2/hbond   seqdep 1.0678 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"
+  if "${basemod} == 2" then     "variable comp equal ${base}+1"     "pair_coeff ${base} ${comp} oxdna2/hbond   seqdep 1.0678 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"
+next base
+jump in.duplex4.4type loop
+variable base loop ${ntype}
+variable base loop 4
+  variable basemod equal ${base}%4
+  variable basemod equal 4%4
+  if "${basemod} == 1" then     "variable comp equal ${base}+3"     "pair_coeff ${base} ${comp} oxdna2/hbond   seqdep 1.0678 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"
+  if "${basemod} == 2" then     "variable comp equal ${base}+1"     "pair_coeff ${base} ${comp} oxdna2/hbond   seqdep 1.0678 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"
+next base
+jump in.duplex4.4type loop
+
+pair_coeff * * oxdna2/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 * * oxdna2/coaxstk 58.5 0.4 0.6 0.22 0.58 2.0 2.891592653589793 0.65 1.3 0 0.8 0.9 0 0.95 0.9 0 0.95 40.0 3.116592653589793
+pair_coeff * * oxdna2/dh      ${T} 0.2 0.815
+pair_coeff * * oxdna2/dh      0.1 0.2 0.815
+
+# Langevin dynamics
+fix 1 all nve/asphere
+fix 2 all langevin ${T} ${T} 25.0 457145 angmom 10
+fix 2 all langevin 0.1 ${T} 25.0 457145 angmom 10
+fix 2 all langevin 0.1 0.1 25.0 457145 angmom 10
+
+timestep 1e-4
+
+#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
+fix 5 all print 10000 "$(step)  ekin = ${ekin} |  erot = ${erot} | epot = ${epot} | etot = ${etot}" screen yes
+
+dump out all custom ${ofreq} out.${number}.txt id mol type x y z ix iy iz c_quat[1] c_quat[2] c_quat[3] c_quat[4]  vx vy vz
+dump out all custom 10000 out.${number}.txt id mol type x y z ix iy iz c_quat[1] c_quat[2] c_quat[3] c_quat[4]  vx vy vz
+dump out all custom 10000 out.1.txt id mol type x y z ix iy iz c_quat[1] c_quat[2] c_quat[3] c_quat[4]  vx vy vz
+dump_modify out sort id
+dump_modify out format line "%d %d %d  %13.6le %13.6le %13.6le %d %d %d  %13.6le %13.6le %13.6le  %13.6le %13.6le %13.6le  %13.6le "
+
+#restart 10000 config0_restart config1_restart
+
+run 100000
+Neighbor list info ...
+  update every 10 steps, delay 0 steps, check yes
+  max neighbors/atom: 2000, page size: 100000
+  master list distance cutoff = 5.63899
+  ghost atom cutoff = 5.63899
+  binsize = 2.81949, bins = 15 15 15
+  6 neighbor lists, perpetual/occasional/extra = 6 0 0
+  (1) pair oxdna2/excv, perpetual
+      attributes: half, newton off
+      pair build: half/bin/newtoff
+      stencil: half/bin/3d/newtoff
+      bin: standard
+  (2) pair oxdna2/stk, perpetual, copy from (1)
+      attributes: half, newton off
+      pair build: copy
+      stencil: none
+      bin: none
+  (3) pair oxdna2/hbond, perpetual, copy from (1)
+      attributes: half, newton off
+      pair build: copy
+      stencil: none
+      bin: none
+  (4) pair oxdna2/xstk, perpetual, copy from (1)
+      attributes: half, newton off
+      pair build: copy
+      stencil: none
+      bin: none
+  (5) pair oxdna2/coaxstk, perpetual, copy from (1)
+      attributes: half, newton off
+      pair build: copy
+      stencil: none
+      bin: none
+  (6) pair oxdna2/dh, perpetual, copy from (1)
+      attributes: half, newton off
+      pair build: copy
+      stencil: none
+      bin: none
+Per MPI rank memory allocation (min/avg/max) = 9.499 | 9.682 | 9.864 Mbytes
+Step Temp E_pair E_mol TotEng Press 
+       0            0   -1.6384018  0.037304147   -1.6010976 6.7769766e-05 
+10000  ekin = 1.20524984175817 |  erot = 1.2679122962101 | epot = -41.9183873815797 | etot = -39.4452252436114
+20000  ekin = 1.60856540355485 |  erot = 2.64140880508964 | epot = -41.0106168994545 | etot = -36.76064269081
+30000  ekin = 2.821101638199 |  erot = 3.06373823252434 | epot = -40.9348041675027 | etot = -35.0499642967793
+40000  ekin = 2.36296917392486 |  erot = 3.53353555114673 | epot = -39.5591676362755 | etot = -33.6626629112039
+50000  ekin = 2.93383938780818 |  erot = 3.39920211088179 | epot = -40.3602175538497 | etot = -34.0271760551597
+60000  ekin = 2.42811021532829 |  erot = 2.75130045540831 | epot = -39.3178531324554 | etot = -34.1384424617188
+70000  ekin = 2.7282200818863 |  erot = 3.76502037022117 | epot = -40.5673105068195 | etot = -34.074070054712
+80000  ekin = 2.90877851656705 |  erot = 2.43617899356904 | epot = -38.4099618426175 | etot = -33.0650043324814
+90000  ekin = 3.89203816080263 |  erot = 2.67194811393274 | epot = -39.1880466354802 | etot = -32.6240603607448
+100000  ekin = 3.98445325397832 |  erot = 2.77079124049636 | epot = -39.2805505658488 | etot = -32.5253060713741
+  100000   0.10625209   -1.5603519  0.049561514   -1.3575422 -1.7755557e-05 
+Loop time of 7.14827 on 4 procs for 100000 steps with 26 atoms
+
+Performance: 120868.376 tau/day, 13989.395 timesteps/s
+99.7% CPU use with 4 MPI tasks x no OpenMP threads
+
+MPI task timing breakdown:
+Section |  min time  |  avg time  |  max time  |%varavg| %total
+---------------------------------------------------------------
+Pair    | 0.047776   | 3.0726     | 6.3143     | 172.8 | 42.98
+Bond    | 0.00736    | 0.056846   | 0.10864    |  20.7 |  0.80
+Neigh   | 1.9e-05    | 1.9e-05    | 1.9e-05    |   0.0 |  0.00
+Comm    | 0.30925    | 3.451      | 6.3739     | 157.5 | 48.28
+Output  | 0.000896   | 0.0010197  | 0.001301   |   0.5 |  0.01
+Modify  | 0.015512   | 0.18417    | 0.37845    |  39.4 |  2.58
+Other   |            | 0.3826     |            |       |  5.35
+
+Nlocal:    6.5 ave 13 max 0 min
+Histogram: 2 0 0 0 0 0 0 0 0 2
+Nghost:    19.5 ave 26 max 13 min
+Histogram: 2 0 0 0 0 0 0 0 0 2
+Neighs:    123.5 ave 247 max 0 min
+Histogram: 2 0 0 0 0 0 0 0 0 2
+
+Total # of neighbors = 494
+Ave neighs/atom = 19
+Ave special neighs/atom = 5.07692
+Neighbor list builds = 1
+Dangerous builds = 0
+
+#write_restart config.${number}.*
+
+
+Total wall time: 0:00:07
diff --git a/examples/USER/cgdna/examples/oxDNA2/unique_bp/log.07Aug19.duplex4.8type.g++1 b/examples/USER/cgdna/examples/oxDNA2/unique_bp/log.07Aug19.duplex4.8type.g++1
new file mode 100644
index 0000000000..c93943a4c7
--- /dev/null
+++ b/examples/USER/cgdna/examples/oxDNA2/unique_bp/log.07Aug19.duplex4.8type.g++1
@@ -0,0 +1,274 @@
+LAMMPS (7 Aug 2019)
+variable number	equal 1
+variable ofreq	equal 10000
+variable efreq	equal 10000
+
+variable ntype equal 8
+
+variable T equal 0.1
+
+units lj
+
+dimension 3
+
+newton on
+
+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 10 delay 0 check yes
+
+read_data data.duplex4.8type
+  orthogonal box = (-20 -20 -20) to (20 20 20)
+  1 by 1 by 1 MPI processor grid
+  reading atoms ...
+  26 atoms
+  reading velocities ...
+  26 velocities
+  26 ellipsoids
+  scanning bonds ...
+  1 = max bonds/atom
+  reading bonds ...
+  24 bonds
+  2 = max # of 1-2 neighbors
+  2 = max # of 1-3 neighbors
+  4 = max # of 1-4 neighbors
+  6 = max # of special neighbors
+  special bonds CPU = 6.5e-05 secs
+  read_data CPU = 0.001139 secs
+mass * 3.1575
+
+group all type 1 8
+26 atoms in group all
+
+# oxDNA bond interactions - FENE backbone
+bond_style oxdna2/fene
+bond_coeff * 2.0 0.25 0.7564
+
+# oxDNA pair interactions
+pair_style hybrid/overlay oxdna2/excv oxdna2/stk oxdna2/hbond oxdna2/xstk oxdna2/coaxstk oxdna2/dh
+pair_coeff * * oxdna2/excv    2.0 0.7 0.675 2.0 0.515 0.5 2.0 0.33 0.32
+pair_coeff * * oxdna2/stk     seqdep ${T} 1.3523 2.6717 6.0 0.4 0.9 0.32 0.75 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 * * oxdna2/stk     seqdep 0.1 1.3523 2.6717 6.0 0.4 0.9 0.32 0.75 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 * * oxdna2/hbond   seqdep 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
+
+label loop
+variable base loop ${ntype}
+variable base loop 8
+  variable basemod equal ${base}%4
+  variable basemod equal 1%4
+  if "${basemod} == 1" then     "variable comp equal ${base}+3"     "pair_coeff ${base} ${comp} oxdna2/hbond   seqdep 1.0678 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"
+variable comp equal ${base}+3
+variable comp equal 1+3
+pair_coeff ${base} ${comp} oxdna2/hbond   seqdep 1.0678 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 ${comp} oxdna2/hbond   seqdep 1.0678 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 oxdna2/hbond   seqdep 1.0678 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
+  if "${basemod} == 2" then     "variable comp equal ${base}+1"     "pair_coeff ${base} ${comp} oxdna2/hbond   seqdep 1.0678 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"
+next base
+jump in.duplex4.8type loop
+variable base loop ${ntype}
+variable base loop 8
+  variable basemod equal ${base}%4
+  variable basemod equal 2%4
+  if "${basemod} == 1" then     "variable comp equal ${base}+3"     "pair_coeff ${base} ${comp} oxdna2/hbond   seqdep 1.0678 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"
+  if "${basemod} == 2" then     "variable comp equal ${base}+1"     "pair_coeff ${base} ${comp} oxdna2/hbond   seqdep 1.0678 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"
+variable comp equal ${base}+1
+variable comp equal 2+1
+pair_coeff ${base} ${comp} oxdna2/hbond   seqdep 1.0678 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 ${comp} oxdna2/hbond   seqdep 1.0678 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 oxdna2/hbond   seqdep 1.0678 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
+next base
+jump in.duplex4.8type loop
+variable base loop ${ntype}
+variable base loop 8
+  variable basemod equal ${base}%4
+  variable basemod equal 3%4
+  if "${basemod} == 1" then     "variable comp equal ${base}+3"     "pair_coeff ${base} ${comp} oxdna2/hbond   seqdep 1.0678 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"
+  if "${basemod} == 2" then     "variable comp equal ${base}+1"     "pair_coeff ${base} ${comp} oxdna2/hbond   seqdep 1.0678 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"
+next base
+jump in.duplex4.8type loop
+variable base loop ${ntype}
+variable base loop 8
+  variable basemod equal ${base}%4
+  variable basemod equal 4%4
+  if "${basemod} == 1" then     "variable comp equal ${base}+3"     "pair_coeff ${base} ${comp} oxdna2/hbond   seqdep 1.0678 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"
+  if "${basemod} == 2" then     "variable comp equal ${base}+1"     "pair_coeff ${base} ${comp} oxdna2/hbond   seqdep 1.0678 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"
+next base
+jump in.duplex4.8type loop
+variable base loop ${ntype}
+variable base loop 8
+  variable basemod equal ${base}%4
+  variable basemod equal 5%4
+  if "${basemod} == 1" then     "variable comp equal ${base}+3"     "pair_coeff ${base} ${comp} oxdna2/hbond   seqdep 1.0678 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"
+variable comp equal ${base}+3
+variable comp equal 5+3
+pair_coeff ${base} ${comp} oxdna2/hbond   seqdep 1.0678 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 5 ${comp} oxdna2/hbond   seqdep 1.0678 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 5 8 oxdna2/hbond   seqdep 1.0678 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
+  if "${basemod} == 2" then     "variable comp equal ${base}+1"     "pair_coeff ${base} ${comp} oxdna2/hbond   seqdep 1.0678 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"
+next base
+jump in.duplex4.8type loop
+variable base loop ${ntype}
+variable base loop 8
+  variable basemod equal ${base}%4
+  variable basemod equal 6%4
+  if "${basemod} == 1" then     "variable comp equal ${base}+3"     "pair_coeff ${base} ${comp} oxdna2/hbond   seqdep 1.0678 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"
+  if "${basemod} == 2" then     "variable comp equal ${base}+1"     "pair_coeff ${base} ${comp} oxdna2/hbond   seqdep 1.0678 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"
+variable comp equal ${base}+1
+variable comp equal 6+1
+pair_coeff ${base} ${comp} oxdna2/hbond   seqdep 1.0678 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 6 ${comp} oxdna2/hbond   seqdep 1.0678 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 6 7 oxdna2/hbond   seqdep 1.0678 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
+next base
+jump in.duplex4.8type loop
+variable base loop ${ntype}
+variable base loop 8
+  variable basemod equal ${base}%4
+  variable basemod equal 7%4
+  if "${basemod} == 1" then     "variable comp equal ${base}+3"     "pair_coeff ${base} ${comp} oxdna2/hbond   seqdep 1.0678 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"
+  if "${basemod} == 2" then     "variable comp equal ${base}+1"     "pair_coeff ${base} ${comp} oxdna2/hbond   seqdep 1.0678 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"
+next base
+jump in.duplex4.8type loop
+variable base loop ${ntype}
+variable base loop 8
+  variable basemod equal ${base}%4
+  variable basemod equal 8%4
+  if "${basemod} == 1" then     "variable comp equal ${base}+3"     "pair_coeff ${base} ${comp} oxdna2/hbond   seqdep 1.0678 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"
+  if "${basemod} == 2" then     "variable comp equal ${base}+1"     "pair_coeff ${base} ${comp} oxdna2/hbond   seqdep 1.0678 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"
+next base
+jump in.duplex4.8type loop
+
+pair_coeff * * oxdna2/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 * * oxdna2/coaxstk 58.5 0.4 0.6 0.22 0.58 2.0 2.891592653589793 0.65 1.3 0 0.8 0.9 0 0.95 0.9 0 0.95 40.0 3.116592653589793
+pair_coeff * * oxdna2/dh      ${T} 0.2 0.815
+pair_coeff * * oxdna2/dh      0.1 0.2 0.815
+
+# Langevin dynamics
+fix 1 all nve/asphere
+fix 2 all langevin ${T} ${T} 25.0 457145 angmom 10
+fix 2 all langevin 0.1 ${T} 25.0 457145 angmom 10
+fix 2 all langevin 0.1 0.1 25.0 457145 angmom 10
+
+timestep 1e-4
+
+#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
+fix 5 all print 10000 "$(step)  ekin = ${ekin} |  erot = ${erot} | epot = ${epot} | etot = ${etot}" screen yes
+
+dump out all custom ${ofreq} out.${number}.txt id mol type x y z ix iy iz c_quat[1] c_quat[2] c_quat[3] c_quat[4]  vx vy vz
+dump out all custom 10000 out.${number}.txt id mol type x y z ix iy iz c_quat[1] c_quat[2] c_quat[3] c_quat[4]  vx vy vz
+dump out all custom 10000 out.1.txt id mol type x y z ix iy iz c_quat[1] c_quat[2] c_quat[3] c_quat[4]  vx vy vz
+dump_modify out sort id
+dump_modify out format line "%d %d %d  %13.6le %13.6le %13.6le %d %d %d  %13.6le %13.6le %13.6le  %13.6le %13.6le %13.6le  %13.6le "
+
+#restart 10000 config0_restart config1_restart
+
+run 100000
+Neighbor list info ...
+  update every 10 steps, delay 0 steps, check yes
+  max neighbors/atom: 2000, page size: 100000
+  master list distance cutoff = 4.63899
+  ghost atom cutoff = 4.63899
+  binsize = 2.31949, bins = 18 18 18
+  6 neighbor lists, perpetual/occasional/extra = 6 0 0
+  (1) pair oxdna2/excv, perpetual
+      attributes: half, newton on
+      pair build: half/bin/newton
+      stencil: half/bin/3d/newton
+      bin: standard
+  (2) pair oxdna2/stk, perpetual, copy from (1)
+      attributes: half, newton on
+      pair build: copy
+      stencil: none
+      bin: none
+  (3) pair oxdna2/hbond, perpetual, copy from (1)
+      attributes: half, newton on
+      pair build: copy
+      stencil: none
+      bin: none
+  (4) pair oxdna2/xstk, perpetual, copy from (1)
+      attributes: half, newton on
+      pair build: copy
+      stencil: none
+      bin: none
+  (5) pair oxdna2/coaxstk, perpetual, copy from (1)
+      attributes: half, newton on
+      pair build: copy
+      stencil: none
+      bin: none
+  (6) pair oxdna2/dh, perpetual, copy from (1)
+      attributes: half, newton on
+      pair build: copy
+      stencil: none
+      bin: none
+Per MPI rank memory allocation (min/avg/max) = 3.905 | 3.905 | 3.905 Mbytes
+Step Temp E_pair E_mol TotEng Press 
+       0            0   -1.6384018  0.037304147   -1.6010976 6.7769766e-05 
+10000  ekin = 1.0129637955345 |  erot = 1.38476360245492 | epot = -41.9785360036034 | etot = -39.580808605614
+20000  ekin = 1.75607695499755 |  erot = 1.86620102096618 | epot = -41.9558496477911 | etot = -38.3335716718274
+30000  ekin = 2.25878673688259 |  erot = 2.60035907744898 | epot = -41.2869827431736 | etot = -36.427836928842
+40000  ekin = 2.11105717434584 |  erot = 3.16611217122799 | epot = -40.2762731426449 | etot = -34.9991037970711
+50000  ekin = 3.05379083545187 |  erot = 2.45050275456888 | epot = -40.3032957287998 | etot = -34.7990021387791
+60000  ekin = 3.05154113920291 |  erot = 2.5253957670889 | epot = -39.6800099891299 | etot = -34.1030730828381
+70000  ekin = 3.23212209366464 |  erot = 3.28914763888543 | epot = -40.1240667487278 | etot = -33.6027970161777
+80000  ekin = 2.82623224207591 |  erot = 2.91292948677732 | epot = -39.0983962904496 | etot = -33.3592345615963
+90000  ekin = 3.05995314905694 |  erot = 3.54299824110274 | epot = -39.164607034397 | etot = -32.5616556442373
+100000  ekin = 3.46139546969892 |  erot = 4.11704803295143 | epot = -38.512467983727 | etot = -30.9340244810767
+  100000  0.092303879   -1.5243833  0.043134544   -1.3481182 6.862374e-06 
+Loop time of 7.94086 on 1 procs for 100000 steps with 26 atoms
+
+Performance: 108804.336 tau/day, 12593.094 timesteps/s
+99.8% CPU use with 1 MPI tasks x no OpenMP threads
+
+MPI task timing breakdown:
+Section |  min time  |  avg time  |  max time  |%varavg| %total
+---------------------------------------------------------------
+Pair    | 7.1507     | 7.1507     | 7.1507     |   0.0 | 90.05
+Bond    | 0.14487    | 0.14487    | 0.14487    |   0.0 |  1.82
+Neigh   | 6.4e-05    | 6.4e-05    | 6.4e-05    |   0.0 |  0.00
+Comm    | 0.032259   | 0.032259   | 0.032259   |   0.0 |  0.41
+Output  | 0.003484   | 0.003484   | 0.003484   |   0.0 |  0.04
+Modify  | 0.59013    | 0.59013    | 0.59013    |   0.0 |  7.43
+Other   |            | 0.01934    |            |       |  0.24
+
+Nlocal:    26 ave 26 max 26 min
+Histogram: 1 0 0 0 0 0 0 0 0 0
+Nghost:    0 ave 0 max 0 min
+Histogram: 1 0 0 0 0 0 0 0 0 0
+Neighs:    323 ave 323 max 323 min
+Histogram: 1 0 0 0 0 0 0 0 0 0
+
+Total # of neighbors = 323
+Ave neighs/atom = 12.4231
+Ave special neighs/atom = 5.07692
+Neighbor list builds = 4
+Dangerous builds = 0
+
+#write_restart config.${number}.*
+
+
+Total wall time: 0:00:07
diff --git a/examples/USER/cgdna/examples/oxDNA2/unique_bp/log.07Aug19.duplex4.8type.g++4 b/examples/USER/cgdna/examples/oxDNA2/unique_bp/log.07Aug19.duplex4.8type.g++4
new file mode 100644
index 0000000000..5cda079e9c
--- /dev/null
+++ b/examples/USER/cgdna/examples/oxDNA2/unique_bp/log.07Aug19.duplex4.8type.g++4
@@ -0,0 +1,274 @@
+LAMMPS (7 Aug 2019)
+variable number	equal 1
+variable ofreq	equal 10000
+variable efreq	equal 10000
+
+variable ntype equal 8
+
+variable T equal 0.1
+
+units lj
+
+dimension 3
+
+newton on
+
+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 10 delay 0 check yes
+
+read_data data.duplex4.8type
+  orthogonal box = (-20 -20 -20) to (20 20 20)
+  1 by 2 by 2 MPI processor grid
+  reading atoms ...
+  26 atoms
+  reading velocities ...
+  26 velocities
+  26 ellipsoids
+  scanning bonds ...
+  1 = max bonds/atom
+  reading bonds ...
+  24 bonds
+  2 = max # of 1-2 neighbors
+  2 = max # of 1-3 neighbors
+  4 = max # of 1-4 neighbors
+  6 = max # of special neighbors
+  special bonds CPU = 0.000232 secs
+  read_data CPU = 0.002447 secs
+mass * 3.1575
+
+group all type 1 8
+26 atoms in group all
+
+# oxDNA bond interactions - FENE backbone
+bond_style oxdna2/fene
+bond_coeff * 2.0 0.25 0.7564
+
+# oxDNA pair interactions
+pair_style hybrid/overlay oxdna2/excv oxdna2/stk oxdna2/hbond oxdna2/xstk oxdna2/coaxstk oxdna2/dh
+pair_coeff * * oxdna2/excv    2.0 0.7 0.675 2.0 0.515 0.5 2.0 0.33 0.32
+pair_coeff * * oxdna2/stk     seqdep ${T} 1.3523 2.6717 6.0 0.4 0.9 0.32 0.75 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 * * oxdna2/stk     seqdep 0.1 1.3523 2.6717 6.0 0.4 0.9 0.32 0.75 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 * * oxdna2/hbond   seqdep 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
+
+label loop
+variable base loop ${ntype}
+variable base loop 8
+  variable basemod equal ${base}%4
+  variable basemod equal 1%4
+  if "${basemod} == 1" then     "variable comp equal ${base}+3"     "pair_coeff ${base} ${comp} oxdna2/hbond   seqdep 1.0678 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"
+variable comp equal ${base}+3
+variable comp equal 1+3
+pair_coeff ${base} ${comp} oxdna2/hbond   seqdep 1.0678 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 ${comp} oxdna2/hbond   seqdep 1.0678 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 oxdna2/hbond   seqdep 1.0678 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
+  if "${basemod} == 2" then     "variable comp equal ${base}+1"     "pair_coeff ${base} ${comp} oxdna2/hbond   seqdep 1.0678 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"
+next base
+jump in.duplex4.8type loop
+variable base loop ${ntype}
+variable base loop 8
+  variable basemod equal ${base}%4
+  variable basemod equal 2%4
+  if "${basemod} == 1" then     "variable comp equal ${base}+3"     "pair_coeff ${base} ${comp} oxdna2/hbond   seqdep 1.0678 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"
+  if "${basemod} == 2" then     "variable comp equal ${base}+1"     "pair_coeff ${base} ${comp} oxdna2/hbond   seqdep 1.0678 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"
+variable comp equal ${base}+1
+variable comp equal 2+1
+pair_coeff ${base} ${comp} oxdna2/hbond   seqdep 1.0678 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 ${comp} oxdna2/hbond   seqdep 1.0678 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 oxdna2/hbond   seqdep 1.0678 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
+next base
+jump in.duplex4.8type loop
+variable base loop ${ntype}
+variable base loop 8
+  variable basemod equal ${base}%4
+  variable basemod equal 3%4
+  if "${basemod} == 1" then     "variable comp equal ${base}+3"     "pair_coeff ${base} ${comp} oxdna2/hbond   seqdep 1.0678 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"
+  if "${basemod} == 2" then     "variable comp equal ${base}+1"     "pair_coeff ${base} ${comp} oxdna2/hbond   seqdep 1.0678 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"
+next base
+jump in.duplex4.8type loop
+variable base loop ${ntype}
+variable base loop 8
+  variable basemod equal ${base}%4
+  variable basemod equal 4%4
+  if "${basemod} == 1" then     "variable comp equal ${base}+3"     "pair_coeff ${base} ${comp} oxdna2/hbond   seqdep 1.0678 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"
+  if "${basemod} == 2" then     "variable comp equal ${base}+1"     "pair_coeff ${base} ${comp} oxdna2/hbond   seqdep 1.0678 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"
+next base
+jump in.duplex4.8type loop
+variable base loop ${ntype}
+variable base loop 8
+  variable basemod equal ${base}%4
+  variable basemod equal 5%4
+  if "${basemod} == 1" then     "variable comp equal ${base}+3"     "pair_coeff ${base} ${comp} oxdna2/hbond   seqdep 1.0678 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"
+variable comp equal ${base}+3
+variable comp equal 5+3
+pair_coeff ${base} ${comp} oxdna2/hbond   seqdep 1.0678 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 5 ${comp} oxdna2/hbond   seqdep 1.0678 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 5 8 oxdna2/hbond   seqdep 1.0678 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
+  if "${basemod} == 2" then     "variable comp equal ${base}+1"     "pair_coeff ${base} ${comp} oxdna2/hbond   seqdep 1.0678 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"
+next base
+jump in.duplex4.8type loop
+variable base loop ${ntype}
+variable base loop 8
+  variable basemod equal ${base}%4
+  variable basemod equal 6%4
+  if "${basemod} == 1" then     "variable comp equal ${base}+3"     "pair_coeff ${base} ${comp} oxdna2/hbond   seqdep 1.0678 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"
+  if "${basemod} == 2" then     "variable comp equal ${base}+1"     "pair_coeff ${base} ${comp} oxdna2/hbond   seqdep 1.0678 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"
+variable comp equal ${base}+1
+variable comp equal 6+1
+pair_coeff ${base} ${comp} oxdna2/hbond   seqdep 1.0678 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 6 ${comp} oxdna2/hbond   seqdep 1.0678 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 6 7 oxdna2/hbond   seqdep 1.0678 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
+next base
+jump in.duplex4.8type loop
+variable base loop ${ntype}
+variable base loop 8
+  variable basemod equal ${base}%4
+  variable basemod equal 7%4
+  if "${basemod} == 1" then     "variable comp equal ${base}+3"     "pair_coeff ${base} ${comp} oxdna2/hbond   seqdep 1.0678 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"
+  if "${basemod} == 2" then     "variable comp equal ${base}+1"     "pair_coeff ${base} ${comp} oxdna2/hbond   seqdep 1.0678 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"
+next base
+jump in.duplex4.8type loop
+variable base loop ${ntype}
+variable base loop 8
+  variable basemod equal ${base}%4
+  variable basemod equal 8%4
+  if "${basemod} == 1" then     "variable comp equal ${base}+3"     "pair_coeff ${base} ${comp} oxdna2/hbond   seqdep 1.0678 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"
+  if "${basemod} == 2" then     "variable comp equal ${base}+1"     "pair_coeff ${base} ${comp} oxdna2/hbond   seqdep 1.0678 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"
+next base
+jump in.duplex4.8type loop
+
+pair_coeff * * oxdna2/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 * * oxdna2/coaxstk 58.5 0.4 0.6 0.22 0.58 2.0 2.891592653589793 0.65 1.3 0 0.8 0.9 0 0.95 0.9 0 0.95 40.0 3.116592653589793
+pair_coeff * * oxdna2/dh      ${T} 0.2 0.815
+pair_coeff * * oxdna2/dh      0.1 0.2 0.815
+
+# Langevin dynamics
+fix 1 all nve/asphere
+fix 2 all langevin ${T} ${T} 25.0 457145 angmom 10
+fix 2 all langevin 0.1 ${T} 25.0 457145 angmom 10
+fix 2 all langevin 0.1 0.1 25.0 457145 angmom 10
+
+timestep 1e-4
+
+#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
+fix 5 all print 10000 "$(step)  ekin = ${ekin} |  erot = ${erot} | epot = ${epot} | etot = ${etot}" screen yes
+
+dump out all custom ${ofreq} out.${number}.txt id mol type x y z ix iy iz c_quat[1] c_quat[2] c_quat[3] c_quat[4]  vx vy vz
+dump out all custom 10000 out.${number}.txt id mol type x y z ix iy iz c_quat[1] c_quat[2] c_quat[3] c_quat[4]  vx vy vz
+dump out all custom 10000 out.1.txt id mol type x y z ix iy iz c_quat[1] c_quat[2] c_quat[3] c_quat[4]  vx vy vz
+dump_modify out sort id
+dump_modify out format line "%d %d %d  %13.6le %13.6le %13.6le %d %d %d  %13.6le %13.6le %13.6le  %13.6le %13.6le %13.6le  %13.6le "
+
+#restart 10000 config0_restart config1_restart
+
+run 100000
+Neighbor list info ...
+  update every 10 steps, delay 0 steps, check yes
+  max neighbors/atom: 2000, page size: 100000
+  master list distance cutoff = 4.63899
+  ghost atom cutoff = 4.63899
+  binsize = 2.31949, bins = 18 18 18
+  6 neighbor lists, perpetual/occasional/extra = 6 0 0
+  (1) pair oxdna2/excv, perpetual
+      attributes: half, newton on
+      pair build: half/bin/newton
+      stencil: half/bin/3d/newton
+      bin: standard
+  (2) pair oxdna2/stk, perpetual, copy from (1)
+      attributes: half, newton on
+      pair build: copy
+      stencil: none
+      bin: none
+  (3) pair oxdna2/hbond, perpetual, copy from (1)
+      attributes: half, newton on
+      pair build: copy
+      stencil: none
+      bin: none
+  (4) pair oxdna2/xstk, perpetual, copy from (1)
+      attributes: half, newton on
+      pair build: copy
+      stencil: none
+      bin: none
+  (5) pair oxdna2/coaxstk, perpetual, copy from (1)
+      attributes: half, newton on
+      pair build: copy
+      stencil: none
+      bin: none
+  (6) pair oxdna2/dh, perpetual, copy from (1)
+      attributes: half, newton on
+      pair build: copy
+      stencil: none
+      bin: none
+Per MPI rank memory allocation (min/avg/max) = 9.381 | 9.563 | 9.746 Mbytes
+Step Temp E_pair E_mol TotEng Press 
+       0            0   -1.6384018  0.037304147   -1.6010976 6.7769766e-05 
+10000  ekin = 1.20524984175816 |  erot = 1.26791229621008 | epot = -41.9183873815797 | etot = -39.4452252436115
+20000  ekin = 1.6085654035548 |  erot = 2.6414088050897 | epot = -41.0106168994546 | etot = -36.7606426908101
+30000  ekin = 2.82110163819891 |  erot = 3.06373823252419 | epot = -40.9348041675026 | etot = -35.0499642967795
+40000  ekin = 2.35553603700794 |  erot = 3.60405945883761 | epot = -39.5418687902286 | etot = -33.5822732943831
+50000  ekin = 3.06061403125833 |  erot = 3.35548669087246 | epot = -40.3236585928169 | etot = -33.9075578706861
+60000  ekin = 2.7347216991605 |  erot = 4.23926242244084 | epot = -39.5379959923263 | etot = -32.564011870725
+70000  ekin = 3.46562604991873 |  erot = 2.8521307045909 | epot = -38.7237000550356 | etot = -32.4059433005259
+80000  ekin = 3.17815543267806 |  erot = 3.11078099027451 | epot = -39.7525036398366 | etot = -33.463567216884
+90000  ekin = 3.51635117668857 |  erot = 2.58384069017333 | epot = -39.4762533092413 | etot = -33.3760614423795
+100000  ekin = 3.64218174280962 |  erot = 2.88607181210736 | epot = -37.6002449519119 | etot = -31.0719913969949
+  100000  0.097124846   -1.5164679  0.070304678   -1.3060794 -2.9111933e-05 
+Loop time of 5.64462 on 4 procs for 100000 steps with 26 atoms
+
+Performance: 153066.031 tau/day, 17715.976 timesteps/s
+99.0% CPU use with 4 MPI tasks x no OpenMP threads
+
+MPI task timing breakdown:
+Section |  min time  |  avg time  |  max time  |%varavg| %total
+---------------------------------------------------------------
+Pair    | 0.069668   | 2.1601     | 4.2189     | 138.7 | 38.27
+Bond    | 0.008542   | 0.052544   | 0.095929   |  18.7 |  0.93
+Neigh   | 4.8e-05    | 5.6e-05    | 6.4e-05    |   0.0 |  0.00
+Comm    | 0.90402    | 3.1443     | 5.4155     | 124.9 | 55.70
+Output  | 0.001003   | 0.0011317  | 0.001468   |   0.6 |  0.02
+Modify  | 0.021098   | 0.20024    | 0.38154    |  39.2 |  3.55
+Other   |            | 0.0863     |            |       |  1.53
+
+Nlocal:    6.5 ave 13 max 0 min
+Histogram: 2 0 0 0 0 0 0 0 0 2
+Nghost:    17.5 ave 22 max 13 min
+Histogram: 2 0 0 0 0 0 0 0 0 2
+Neighs:    80.25 ave 167 max 0 min
+Histogram: 2 0 0 0 0 0 0 0 0 2
+
+Total # of neighbors = 321
+Ave neighs/atom = 12.3462
+Ave special neighs/atom = 5.07692
+Neighbor list builds = 3
+Dangerous builds = 0
+
+#write_restart config.${number}.*
+
+
+Total wall time: 0:00:05