diff --git a/examples/PACKAGES/cgdna/examples/test.sh b/examples/PACKAGES/cgdna/examples/test.sh
index 795be9ef88..152047b94b 100755
--- a/examples/PACKAGES/cgdna/examples/test.sh
+++ b/examples/PACKAGES/cgdna/examples/test.sh
@@ -1,20 +1,24 @@
 #! /bin/bash
 
-DATE='2Jul21'
-LMPDIR=/Users/ohenrich/Work/code/lammps
+DATE='14Dec21'
+TOL=1e-8
 
+LMPDIR=/Users/ohenrich/Work/code/lammps
 SRCDIR=$LMPDIR/src
 EXDIR=$LMPDIR/examples/PACKAGES/cgdna/examples
 
 if [ $# -eq 1 ] && [ $1 = run ]; then
-  echo '# Compiling executable in' $SRCDIR
+  echo '# Compiling executable in' $SRCDIR | tee -a $EXDIR/test.log
 
   cd $SRCDIR
-  make clean-all
-  make -j8 mpi
+  make clean-all | tee -a $EXDIR/test.log
+  make purge | tee -a $EXDIR/test.log
+  make pu | tee -a $EXDIR/test.log
+  make ps | tee -a $EXDIR/test.log
+  make -j8 mpi | tee -a $EXDIR/test.log
 
   ######################################################
-  echo '# Running oxDNA duplex1 test'
+  echo '# Running oxDNA duplex1 test' | tee -a $EXDIR/test.log
   cd $EXDIR/oxDNA/duplex1
   mkdir test
   cd test
@@ -24,20 +28,32 @@ if [ $# -eq 1 ] && [ $1 = run ]; then
 
   mpirun -np 1 ./lmp_mpi < in.duplex1 > /dev/null
   mv log.lammps log.$DATE.duplex1.g++.1
-  grep etot log.$DATE.duplex1.g++.1 > e_test.dat
-  grep etot ../log*1 > e_old.dat
-  ndiff -relerr 1e-8 e_test.dat e_old.dat
+  grep etot log.$DATE.duplex1.g++.1 > e_test.1.dat
+  grep etot ../log*1 > e_old.1.dat
+  ndiff -relerr $TOL e_test.1.dat e_old.1.dat
+  if [ $? -eq 0 ];
+  then 
+      echo "# 1 MPI-task passed" | tee -a $EXDIR/test.log
+  else 
+      echo "# 1 MPI-task unsuccessful" | tee -a $EXDIR/test.log
+  fi
 
   mpirun -np 4 ./lmp_mpi < in.duplex1 > /dev/null
   mv log.lammps log.$DATE.duplex1.g++.4
-  grep etot log.$DATE.duplex1.g++.4 > e_test.dat
-  grep etot ../log*4 > e_old.dat
-  ndiff -relerr 1e-8 e_test.dat e_old.dat
+  grep etot log.$DATE.duplex1.g++.4 > e_test.4.dat
+  grep etot ../log*4 > e_old.4.dat
+  ndiff -relerr $TOL e_test.4.dat e_old.4.dat
+  if [ $? -eq 0 ];
+  then 
+      echo "# 4 MPI-tasks passed" | tee -a $EXDIR/test.log
+  else 
+      echo "# 4 MPI-tasks unsuccessful" | tee -a $EXDIR/test.log
+  fi
 
   ######################################################
 
   ######################################################
-  echo '# Running oxDNA duplex2 test'
+  echo '# Running oxDNA duplex2 test' | tee -a $EXDIR/test.log
   cd $EXDIR/oxDNA/duplex2
   mkdir test
   cd test
@@ -47,20 +63,32 @@ if [ $# -eq 1 ] && [ $1 = run ]; then
 
   mpirun -np 1 ./lmp_mpi < in.duplex2 > /dev/null
   mv log.lammps log.$DATE.duplex2.g++.1
-  grep etot log.$DATE.duplex2.g++.1 > e_test.dat
-  grep etot ../log*1 > e_old.dat
-  ndiff -relerr 1e-8 e_test.dat e_old.dat
+  grep etot log.$DATE.duplex2.g++.1 > e_test.1.dat
+  grep etot ../log*1 > e_old.1.dat
+  ndiff -relerr $TOL e_test.1.dat e_old.1.dat
+  if [ $? -eq 0 ];
+  then 
+      echo "# 1 MPI-task passed" | tee -a $EXDIR/test.log
+  else 
+      echo "# 1 MPI-task unsuccessful" | tee -a $EXDIR/test.log
+  fi
 
   mpirun -np 4 ./lmp_mpi < in.duplex2 > /dev/null
   mv log.lammps log.$DATE.duplex2.g++.4
-  grep etot log.$DATE.duplex2.g++.4 > e_test.dat
-  grep etot ../log*4 > e_old.dat
-  ndiff -relerr 1e-8 e_test.dat e_old.dat
+  grep etot log.$DATE.duplex2.g++.4 > e_test.4.dat
+  grep etot ../log*4 > e_old.4.dat
+  ndiff -relerr $TOL e_test.4.dat e_old.4.dat
+  if [ $? -eq 0 ];
+  then 
+      echo "# 4 MPI-tasks passed" | tee -a $EXDIR/test.log
+  else 
+      echo "# 4 MPI-tasks unsuccessful" | tee -a $EXDIR/test.log
+  fi
 
   ######################################################
 
   ######################################################
-  echo '# Running oxDNA2 duplex1 test'
+  echo '# Running oxDNA2 duplex1 test' | tee -a $EXDIR/test.log
   cd $EXDIR/oxDNA2/duplex1
   mkdir test
   cd test
@@ -70,20 +98,32 @@ if [ $# -eq 1 ] && [ $1 = run ]; then
 
   mpirun -np 1 ./lmp_mpi < in.duplex1 > /dev/null
   mv log.lammps log.$DATE.duplex1.g++.1
-  grep etot log.$DATE.duplex1.g++.1 > e_test.dat
-  grep etot ../log*1 > e_old.dat
-  ndiff -relerr 1e-8 e_test.dat e_old.dat
+  grep etot log.$DATE.duplex1.g++.1 > e_test.1.dat
+  grep etot ../log*1 > e_old.1.dat
+  ndiff -relerr $TOL e_test.1.dat e_old.1.dat
+  if [ $? -eq 0 ];
+  then 
+      echo "# 1 MPI-task passed" | tee -a $EXDIR/test.log
+  else 
+      echo "# 1 MPI-task unsuccessful" | tee -a $EXDIR/test.log
+  fi
 
   mpirun -np 4 ./lmp_mpi < in.duplex1 > /dev/null
   mv log.lammps log.$DATE.duplex1.g++.4
-  grep etot log.$DATE.duplex1.g++.4 > e_test.dat
-  grep etot ../log*4 > e_old.dat
-  ndiff -relerr 1e-8 e_test.dat e_old.dat
+  grep etot log.$DATE.duplex1.g++.4 > e_test.4.dat
+  grep etot ../log*4 > e_old.4.dat
+  ndiff -relerr $TOL e_test.4.dat e_old.4.dat
+  if [ $? -eq 0 ];
+  then 
+      echo "# 4 MPI-tasks passed" | tee -a $EXDIR/test.log
+  else 
+      echo "# 4 MPI-tasks unsuccessful" | tee -a $EXDIR/test.log
+  fi
 
   ######################################################
 
   ######################################################
-  echo '# Running oxDNA2 duplex2 test'
+  echo '# Running oxDNA2 duplex2 test' | tee -a $EXDIR/test.log
   cd $EXDIR/oxDNA2/duplex2
   mkdir test
   cd test
@@ -93,20 +133,32 @@ if [ $# -eq 1 ] && [ $1 = run ]; then
 
   mpirun -np 1 ./lmp_mpi < in.duplex2 > /dev/null
   mv log.lammps log.$DATE.duplex2.g++.1
-  grep etot log.$DATE.duplex2.g++.1 > e_test.dat
-  grep etot ../log*1 > e_old.dat
-  ndiff -relerr 1e-8 e_test.dat e_old.dat
+  grep etot log.$DATE.duplex2.g++.1 > e_test.1.dat
+  grep etot ../log*1 > e_old.1.dat
+  ndiff -relerr $TOL e_test.1.dat e_old.1.dat
+  if [ $? -eq 0 ];
+  then 
+      echo "# 1 MPI-task passed" | tee -a $EXDIR/test.log
+  else 
+      echo "# 1 MPI-task unsuccessful" | tee -a $EXDIR/test.log
+  fi
 
   mpirun -np 4 ./lmp_mpi < in.duplex2 > /dev/null
   mv log.lammps log.$DATE.duplex2.g++.4
-  grep etot log.$DATE.duplex2.g++.4 > e_test.dat
-  grep etot ../log*4 > e_old.dat
-  ndiff -relerr 1e-8 e_test.dat e_old.dat
+  grep etot log.$DATE.duplex2.g++.4 > e_test.4.dat
+  grep etot ../log*4 > e_old.4.dat
+  ndiff -relerr $TOL e_test.4.dat e_old.4.dat
+  if [ $? -eq 0 ];
+  then 
+      echo "# 4 MPI-tasks passed" | tee -a $EXDIR/test.log
+  else 
+      echo "# 4 MPI-tasks unsuccessful" | tee -a $EXDIR/test.log
+  fi
 
   ######################################################
 
   ######################################################
-  echo '# Running oxDNA2 duplex3 test'
+  echo '# Running oxDNA2 duplex3 test' | tee -a $EXDIR/test.log
   cd $EXDIR/oxDNA2/duplex3
   mkdir test
   cd test
@@ -116,20 +168,32 @@ if [ $# -eq 1 ] && [ $1 = run ]; then
 
   mpirun -np 1 ./lmp_mpi < in.duplex3 > /dev/null
   mv log.lammps log.$DATE.duplex3.g++.1
-  grep etot log.$DATE.duplex3.g++.1 > e_test.dat
-  grep etot ../log*1 > e_old.dat
-  ndiff -relerr 1e-8 e_test.dat e_old.dat
+  grep etot log.$DATE.duplex3.g++.1 > e_test.1.dat
+  grep etot ../log*1 > e_old.1.dat
+  ndiff -relerr $TOL e_test.1.dat e_old.1.dat
+  if [ $? -eq 0 ];
+  then 
+      echo "# 1 MPI-task passed" | tee -a $EXDIR/test.log
+  else 
+      echo "# 1 MPI-task unsuccessful" | tee -a $EXDIR/test.log
+  fi
 
   mpirun -np 4 ./lmp_mpi < in.duplex3 > /dev/null
   mv log.lammps log.$DATE.duplex3.g++.4
-  grep etot  log.$DATE.duplex3.g++.4 > e_test.dat
-  grep etot ../log*4 > e_old.dat
-  ndiff -relerr 1e-8 e_test.dat e_old.dat
+  grep etot  log.$DATE.duplex3.g++.4 > e_test.4.dat
+  grep etot ../log*4 > e_old.4.dat
+  ndiff -relerr $TOL e_test.4.dat e_old.4.dat
+  if [ $? -eq 0 ];
+  then 
+      echo "# 4 MPI-tasks passed" | tee -a $EXDIR/test.log
+  else 
+      echo "# 4 MPI-tasks unsuccessful" | tee -a $EXDIR/test.log
+  fi
 
   ######################################################
 
   ######################################################
-  echo '# Running oxDNA2 unique_bp test'
+  echo '# Running oxDNA2 unique_bp test' | tee -a $EXDIR/test.log
   cd $EXDIR/oxDNA2/unique_bp
   mkdir test
   cd test
@@ -141,32 +205,56 @@ if [ $# -eq 1 ] && [ $1 = run ]; then
 
   mpirun -np 1 ./lmp_mpi < in.duplex4.4type > /dev/null
   mv log.lammps log.$DATE.duplex4.4type.g++.1
-  grep etot log.$DATE.duplex4.4type.g++.1 > e_test.dat
-  grep etot ../log*4type*1 > e_old.dat
-  ndiff -relerr 1e-8 e_test.dat e_old.dat
+  grep etot log.$DATE.duplex4.4type.g++.1 > e_test.4type.1.dat
+  grep etot ../log*4type*1 > e_old.4type.1.dat
+  ndiff -relerr $TOL e_test.4type.1.dat e_old.4type.1.dat
+  if [ $? -eq 0 ];
+  then 
+      echo "# 1 MPI-task 4 types passed" | tee -a $EXDIR/test.log
+  else 
+      echo "# 1 MPI-task 4 types unsuccessful" | tee -a $EXDIR/test.log
+  fi
 
   mpirun -np 4 ./lmp_mpi < in.duplex4.4type > /dev/null
   mv log.lammps log.$DATE.duplex4.4type.g++.4
-  grep etot log.$DATE.duplex4.4type.g++.4 > e_test.dat
-  grep etot ../log*4type*4 > e_old.dat
-  ndiff -relerr 1e-8 e_test.dat e_old.dat
+  grep etot log.$DATE.duplex4.4type.g++.4 > e_test.4type.4.dat
+  grep etot ../log*4type*4 > e_old.4type.4.dat
+  ndiff -relerr $TOL e_test.4type.4.dat e_old.4type.4.dat
+  if [ $? -eq 0 ];
+  then 
+      echo "# 4 MPI-tasks 4 types passed" | tee -a $EXDIR/test.log
+  else 
+      echo "# 4 MPI-tasks 4 types unsuccessful" | tee -a $EXDIR/test.log
+  fi
 
   mpirun -np 1 ./lmp_mpi < in.duplex4.8type > /dev/null
   mv log.lammps log.$DATE.duplex4.8type.g++.1
-  grep etot log.$DATE.duplex4.8type.g++.1 > e_test.dat
-  grep etot ../log*8type*1 > e_old.dat
-  ndiff -relerr 1e-8 e_test.dat e_old.dat
+  grep etot log.$DATE.duplex4.8type.g++.1 > e_test.8type.1.dat
+  grep etot ../log*8type*1 > e_old.8type.1.dat
+  ndiff -relerr $TOL e_test.8type.1.dat e_old.8type.1.dat
+  if [ $? -eq 0 ];
+  then 
+      echo "# 1 MPI-task 8 types passed" | tee -a $EXDIR/test.log
+  else 
+      echo "# 1 MPI-task 8 types unsuccessful" | tee -a $EXDIR/test.log
+  fi
 
   mpirun -np 4 ./lmp_mpi < in.duplex4.8type > /dev/null
   mv log.lammps log.$DATE.duplex4.8type.g++.4
-  grep etot log.$DATE.duplex4.8type.g++.4 > e_test.dat
-  grep etot ../log*8type*4 > e_old.dat
-  ndiff -relerr 1e-8 e_test.dat e_old.dat
+  grep etot log.$DATE.duplex4.8type.g++.4 > e_test.8type.4.dat
+  grep etot ../log*8type*4 > e_old.8type.4.dat
+  ndiff -relerr $TOL e_test.8type.4.dat e_old.8type.4.dat
+  if [ $? -eq 0 ];
+  then 
+      echo "# 4 MPI-tasks 8 types passed" | tee -a $EXDIR/test.log
+  else 
+      echo "# 4 MPI-tasks 8 types unsuccessful" | tee -a $EXDIR/test.log
+  fi
 
   ######################################################
 
   ######################################################
-  echo '# Running oxDNA2 dsring test'
+  echo '# Running oxDNA2 dsring test' | tee -a $EXDIR/test.log
   cd $EXDIR/oxDNA2/dsring
   mkdir test
   cd test
@@ -176,20 +264,32 @@ if [ $# -eq 1 ] && [ $1 = run ]; then
 
   mpirun -np 1 ./lmp_mpi < in.dsring > /dev/null
   mv log.lammps log.$DATE.dsring.g++.1
-  grep etot log.$DATE.dsring.g++.1 > e_test.dat
-  grep etot ../log*1 > e_old.dat
-  ndiff -relerr 1e-8 e_test.dat e_old.dat
+  grep etot log.$DATE.dsring.g++.1 > e_test.1.dat
+  grep etot ../log*1 > e_old.1.dat
+  ndiff -relerr $TOL e_test.1.dat e_old.1.dat
+  if [ $? -eq 0 ];
+  then 
+      echo "# 1 MPI-task passed" | tee -a $EXDIR/test.log
+  else 
+      echo "# 1 MPI-task unsuccessful" | tee -a $EXDIR/test.log
+  fi
 
   mpirun -np 4 ./lmp_mpi < in.dsring > /dev/null
   mv log.lammps log.$DATE.dsring.g++.4
-  grep etot log.$DATE.dsring.g++.4 > e_test.dat
-  grep etot ../log*4 > e_old.dat
-  ndiff -relerr 1e-8 e_test.dat e_old.dat
+  grep etot log.$DATE.dsring.g++.4 > e_test.4.dat
+  grep etot ../log*4 > e_old.4.dat
+  ndiff -relerr $TOL e_test.4.dat e_old.4.dat
+  if [ $? -eq 0 ];
+  then 
+      echo "# 4 MPI-tasks passed" | tee -a $EXDIR/test.log
+  else 
+      echo "# 4 MPI-tasks unsuccessful" | tee -a $EXDIR/test.log
+  fi
 
   ######################################################
 
   ######################################################
-  echo '# Running oxRNA2 duplex2 test'
+  echo '# Running oxRNA2 duplex2 test' | tee -a $EXDIR/test.log
   cd $EXDIR/oxRNA2/duplex2
   mkdir test
   cd test
@@ -199,18 +299,30 @@ if [ $# -eq 1 ] && [ $1 = run ]; then
 
   mpirun -np 1 ./lmp_mpi < in.duplex2 > /dev/null
   mv log.lammps log.$DATE.duplex2.g++.1
-  grep etot log.$DATE.duplex2.g++.1 > e_test.dat
-  grep etot ../log*1 > e_old.dat
-  ndiff -relerr 1e-8 e_test.dat e_old.dat
+  grep etot log.$DATE.duplex2.g++.1 > e_test.1.dat
+  grep etot ../log*1 > e_old.1.dat
+  ndiff -relerr $TOL e_test.1.dat e_old.1.dat
+  if [ $? -eq 0 ];
+  then 
+      echo "# 1 MPI-task passed" | tee -a $EXDIR/test.log
+  else 
+      echo "# 1 MPI-task unsuccessful" | tee -a $EXDIR/test.log
+  fi
 
   mpirun -np 4 ./lmp_mpi < in.duplex2 > /dev/null
   mv log.lammps log.$DATE.duplex2.g++.4
-  grep etot log.$DATE.duplex2.g++.4 > e_test.dat
-  grep etot ../log*4 > e_old.dat
-  ndiff -relerr 1e-8 e_test.dat e_old.dat
+  grep etot log.$DATE.duplex2.g++.4 > e_test.4.dat
+  grep etot ../log*4 > e_old.4.dat
+  ndiff -relerr $TOL e_test.4.dat e_old.4.dat
+  if [ $? -eq 0 ];
+  then 
+      echo "# 4 MPI-tasks passed" | tee -a $EXDIR/test.log
+  else 
+      echo "# 4 MPI-tasks unsuccessful" | tee -a $EXDIR/test.log
+  fi
 
   ######################################################
-  echo '# Done'
+  echo '# Done' | tee -a $EXDIR/test.log
 
 elif [ $# -eq 1 ] && [ $1 = clean ]; then
   echo '# Deleting test directories'
@@ -222,6 +334,7 @@ elif [ $# -eq 1 ] && [ $1 = clean ]; then
   rm -rf $EXDIR/oxDNA2/unique_bp/test
   rm -rf $EXDIR/oxDNA2/dsring/test
   rm -rf $EXDIR/oxRNA2/duplex2/test
+  rm -rf $EXDIR/test.log
   echo '# Done'
 
 else
diff --git a/src/CG-DNA/bond_oxdna_fene.cpp b/src/CG-DNA/bond_oxdna_fene.cpp
index 6b720f6a7f..4db0fb4cf1 100644
--- a/src/CG-DNA/bond_oxdna_fene.cpp
+++ b/src/CG-DNA/bond_oxdna_fene.cpp
@@ -26,6 +26,7 @@
 
 #include "atom_vec_ellipsoid.h"
 #include "math_extra.h"
+#include "pair.h"
 
 #include <cmath>
 
@@ -145,16 +146,16 @@ void BondOxdnaFene::ev_tally_xyz(int i, int j, int nlocal, int newton_bond, doub
 ------------------------------------------------------------------------- */
 void BondOxdnaFene::compute(int eflag, int vflag)
 {
-  int a, b, in, type;
-  double delf[3], delta[3], deltb[3];    // force, torque increment;;
-  double delr[3], ebond, fbond;
-  double rsq, Deltasq, rlogarg;
-  double r, rr0, rr0sq;
+  int a,b,in,type;
+  double delf[3],delta[3],deltb[3];    // force, torque increment;;
+  double delr[3],ebond,fbond;
+  double rsq,Deltasq,rlogarg;
+  double r,rr0,rr0sq;
   // vectors COM-backbone site in lab frame
-  double ra_cs[3], rb_cs[3];
-
-  double *qa, ax[3], ay[3], az[3];
-  double *qb, bx[3], by[3], bz[3];
+  double ra_cs[3],rb_cs[3];
+  // Cartesian unit vectors in lab frame
+  double ax[3],ay[3],az[3];
+  double bx[3],by[3],bz[3];
 
   double **x = atom->x;
   double **f = atom->f;
@@ -172,6 +173,12 @@ void BondOxdnaFene::compute(int eflag, int vflag)
   ebond = 0.0;
   ev_init(eflag, vflag);
 
+  // n(x/y/z)_xtrct = extracted local unit vectors in lab frame from oxdna_excv
+  int dim;
+  nx_xtrct = (double **) force->pair->extract("nx",dim);
+  ny_xtrct = (double **) force->pair->extract("ny",dim);
+  nz_xtrct = (double **) force->pair->extract("nz",dim);
+
   // loop over FENE bonds
 
   for (in = 0; in < nbondlist; in++) {
@@ -180,10 +187,24 @@ void BondOxdnaFene::compute(int eflag, int vflag)
     b = bondlist[in][0];
     type = bondlist[in][2];
 
-    qa = bonus[ellipsoid[a]].quat;
-    MathExtra::q_to_exyz(qa, ax, ay, az);
-    qb = bonus[ellipsoid[b]].quat;
-    MathExtra::q_to_exyz(qb, bx, by, bz);
+    ax[0] = nx_xtrct[a][0];
+    ax[1] = nx_xtrct[a][1];
+    ax[2] = nx_xtrct[a][2];
+    ay[0] = ny_xtrct[a][0];
+    ay[1] = ny_xtrct[a][1];
+    ay[2] = ny_xtrct[a][2];
+    az[0] = nz_xtrct[a][0];
+    az[1] = nz_xtrct[a][1];
+    az[2] = nz_xtrct[a][2];
+    bx[0] = nx_xtrct[b][0];
+    bx[1] = nx_xtrct[b][1];
+    bx[2] = nx_xtrct[b][2];
+    by[0] = ny_xtrct[b][0];
+    by[1] = ny_xtrct[b][1];
+    by[2] = ny_xtrct[b][2];
+    bz[0] = nz_xtrct[b][0];
+    bz[1] = nz_xtrct[b][1];
+    bz[2] = nz_xtrct[b][2];
 
     // vector COM-backbone site a and b
     compute_interaction_sites(ax, ay, az, ra_cs);
diff --git a/src/CG-DNA/bond_oxdna_fene.h b/src/CG-DNA/bond_oxdna_fene.h
index d2e4eed183..d2aa84612a 100644
--- a/src/CG-DNA/bond_oxdna_fene.h
+++ b/src/CG-DNA/bond_oxdna_fene.h
@@ -40,6 +40,7 @@ class BondOxdnaFene : public Bond {
 
  protected:
   double *k, *Delta, *r0;    // FENE
+  double **nx_xtrct, **ny_xtrct, **nz_xtrct; // per-atom arrays for local unit vectors
 
   void allocate();
   void ev_tally_xyz(int, int, int, int, double, double, double, double, double, double, double);
diff --git a/src/CG-DNA/pair_oxdna2_coaxstk.cpp b/src/CG-DNA/pair_oxdna2_coaxstk.cpp
index e8ab197be1..8b3609983b 100644
--- a/src/CG-DNA/pair_oxdna2_coaxstk.cpp
+++ b/src/CG-DNA/pair_oxdna2_coaxstk.cpp
@@ -118,9 +118,9 @@ void PairOxdna2Coaxstk::compute(int eflag, int vflag)
   double ra_cs[3],ra_cst[3];
   double rb_cs[3],rb_cst[3];
 
-  // quaternions and Cartesian unit vectors in lab frame
-  double *qa,ax[3],ay[3],az[3];
-  double *qb,bx[3],by[3],bz[3];
+  // Cartesian unit vectors in lab frame
+  double ax[3],ay[3],az[3];
+  double bx[3],by[3],bz[3];
 
   double **x = atom->x;
   double **f = atom->f;
@@ -149,6 +149,11 @@ void PairOxdna2Coaxstk::compute(int eflag, int vflag)
   numneigh = list->numneigh;
   firstneigh = list->firstneigh;
 
+  // n(x/z)_xtrct = extracted local unit vectors from oxdna_excv
+  int dim;
+  nx_xtrct = (double **) force->pair->extract("nx",dim);
+  nz_xtrct = (double **) force->pair->extract("nz",dim);
+
   // loop over pair interaction neighbors of my atoms
 
   for (ia = 0; ia < anum; ia++) {
@@ -156,8 +161,9 @@ void PairOxdna2Coaxstk::compute(int eflag, int vflag)
     a = alist[ia];
     atype = type[a];
 
-    qa=bonus[ellipsoid[a]].quat;
-    MathExtra::q_to_exyz(qa,ax,ay,az);
+    ax[0] = nx_xtrct[a][0];
+    ax[1] = nx_xtrct[a][1];
+    ax[2] = nx_xtrct[a][2];
 
     // vector COM a - stacking site a
     ra_cst[0] = d_cst*ax[0];
@@ -180,8 +186,9 @@ void PairOxdna2Coaxstk::compute(int eflag, int vflag)
 
       btype = type[b];
 
-      qb=bonus[ellipsoid[b]].quat;
-      MathExtra::q_to_exyz(qb,bx,by,bz);
+      bx[0] = nx_xtrct[b][0];
+      bx[1] = nx_xtrct[b][1];
+      bx[2] = nx_xtrct[b][2];
 
       // vector COM b - stacking site b
       rb_cst[0] = d_cst*bx[0];
@@ -232,6 +239,13 @@ void PairOxdna2Coaxstk::compute(int eflag, int vflag)
       // early rejection criterium
       if (f4f6t1) {
 
+      az[0] = nz_xtrct[a][0];
+      az[1] = nz_xtrct[a][1];
+      az[2] = nz_xtrct[a][2];
+      bz[0] = nz_xtrct[b][0];
+      bz[1] = nz_xtrct[b][1];
+      bz[2] = nz_xtrct[b][2];
+
       cost4 = MathExtra::dot3(az,bz);
       if (cost4 >  1.0) cost4 =  1.0;
       if (cost4 < -1.0) cost4 = -1.0;
@@ -306,7 +320,7 @@ void PairOxdna2Coaxstk::compute(int eflag, int vflag)
               DF4(theta6p, a_cxst6[atype][btype], theta_cxst6_0[atype][btype], dtheta_cxst6_ast[atype][btype],
               b_cxst6[atype][btype], dtheta_cxst6_c[atype][btype])*rsint;
 
-     // force, torque and virial contribution for forces between stacking sites
+      // force, torque and virial contribution for forces between stacking sites
 
       fpair = 0.0;
 
diff --git a/src/CG-DNA/pair_oxdna2_coaxstk.h b/src/CG-DNA/pair_oxdna2_coaxstk.h
index 4fe52b6e9b..6bf763a1fe 100644
--- a/src/CG-DNA/pair_oxdna2_coaxstk.h
+++ b/src/CG-DNA/pair_oxdna2_coaxstk.h
@@ -46,20 +46,16 @@ class PairOxdna2Coaxstk : public Pair {
   double **k_cxst, **cut_cxst_0, **cut_cxst_c, **cut_cxst_lo, **cut_cxst_hi;
   double **cut_cxst_lc, **cut_cxst_hc, **b_cxst_lo, **b_cxst_hi;
   double **cutsq_cxst_hc;
-
   double **a_cxst1, **theta_cxst1_0, **dtheta_cxst1_ast;
   double **b_cxst1, **dtheta_cxst1_c;
-
   double **a_cxst4, **theta_cxst4_0, **dtheta_cxst4_ast;
   double **b_cxst4, **dtheta_cxst4_c;
-
   double **a_cxst5, **theta_cxst5_0, **dtheta_cxst5_ast;
   double **b_cxst5, **dtheta_cxst5_c;
-
   double **a_cxst6, **theta_cxst6_0, **dtheta_cxst6_ast;
   double **b_cxst6, **dtheta_cxst6_c;
-
   double **AA_cxst1, **BB_cxst1;
+  double **nx_xtrct, **nz_xtrct; // per-atom arrays for local unit vectors
 
   virtual void allocate();
 };
diff --git a/src/CG-DNA/pair_oxdna2_dh.cpp b/src/CG-DNA/pair_oxdna2_dh.cpp
index 1a587e82e1..cf79fd07cf 100644
--- a/src/CG-DNA/pair_oxdna2_dh.cpp
+++ b/src/CG-DNA/pair_oxdna2_dh.cpp
@@ -87,10 +87,9 @@ void PairOxdna2Dh::compute(int eflag, int vflag)
   // vectors COM-backbone sites in lab frame
   double ra_cs[3],rb_cs[3];
 
-  // quaternions and Cartesian unit vectors in lab frame
-  double *qa,ax[3],ay[3],az[3];
-  double *qb,bx[3],by[3],bz[3];
-  double *special_lj = force->special_lj;
+  // Cartesian unit vectors in lab frame
+  double ax[3],ay[3],az[3];
+  double bx[3],by[3],bz[3];
 
   double **x = atom->x;
   double **f = atom->f;
@@ -100,6 +99,7 @@ void PairOxdna2Dh::compute(int eflag, int vflag)
   int nlocal = atom->nlocal;
   int newton_pair = force->newton_pair;
   int *alist,*blist,*numneigh,**firstneigh;
+  double *special_lj = force->special_lj;
 
   AtomVecEllipsoid *avec = (AtomVecEllipsoid *) atom->style_match("ellipsoid");
   AtomVecEllipsoid::Bonus *bonus = avec->bonus;
@@ -115,6 +115,12 @@ void PairOxdna2Dh::compute(int eflag, int vflag)
   numneigh = list->numneigh;
   firstneigh = list->firstneigh;
 
+  // n(x/y/z)_xtrct = extracted local unit vectors from oxdna_excv
+  int dim;
+  nx_xtrct = (double **) force->pair->extract("nx",dim);
+  ny_xtrct = (double **) force->pair->extract("ny",dim);
+  nz_xtrct = (double **) force->pair->extract("nz",dim);
+
   // loop over pair interaction neighbors of my atoms
 
   for (ia = 0; ia < anum; ia++) {
@@ -122,8 +128,15 @@ void PairOxdna2Dh::compute(int eflag, int vflag)
     a = alist[ia];
     atype = type[a];
 
-    qa=bonus[ellipsoid[a]].quat;
-    MathExtra::q_to_exyz(qa,ax,ay,az);
+    ax[0] = nx_xtrct[a][0];
+    ax[1] = nx_xtrct[a][1];
+    ax[2] = nx_xtrct[a][2];
+    ay[0] = ny_xtrct[a][0];
+    ay[1] = ny_xtrct[a][1];
+    ay[2] = ny_xtrct[a][2];
+    az[0] = nz_xtrct[a][0];
+    az[1] = nz_xtrct[a][1];
+    az[2] = nz_xtrct[a][2];
 
     // vector COM-backbone site a
     compute_interaction_sites(ax,ay,az,ra_cs);
@@ -142,8 +155,15 @@ void PairOxdna2Dh::compute(int eflag, int vflag)
       b &= NEIGHMASK;
       btype = type[b];
 
-      qb=bonus[ellipsoid[b]].quat;
-      MathExtra::q_to_exyz(qb,bx,by,bz);
+      bx[0] = nx_xtrct[b][0];
+      bx[1] = nx_xtrct[b][1];
+      bx[2] = nx_xtrct[b][2];
+      by[0] = ny_xtrct[b][0];
+      by[1] = ny_xtrct[b][1];
+      by[2] = ny_xtrct[b][2];
+      bz[0] = nz_xtrct[b][0];
+      bz[1] = nz_xtrct[b][1];
+      bz[2] = nz_xtrct[b][2];
 
       // vector COM-backbone site b
       compute_interaction_sites(bx,by,bz,rb_cs);
diff --git a/src/CG-DNA/pair_oxdna2_dh.h b/src/CG-DNA/pair_oxdna2_dh.h
index 5c79641935..db37ba0d43 100644
--- a/src/CG-DNA/pair_oxdna2_dh.h
+++ b/src/CG-DNA/pair_oxdna2_dh.h
@@ -45,6 +45,7 @@ class PairOxdna2Dh : public Pair {
  protected:
   double **qeff_dh_pf, **kappa_dh;
   double **b_dh, **cut_dh_ast, **cutsq_dh_ast, **cut_dh_c, **cutsq_dh_c;
+  double **nx_xtrct, **ny_xtrct, **nz_xtrct; // per-atom arrays for local unit vectors
 
   virtual void allocate();
 };
diff --git a/src/CG-DNA/pair_oxdna_coaxstk.cpp b/src/CG-DNA/pair_oxdna_coaxstk.cpp
index 2715b7c5e2..49b25bb81b 100644
--- a/src/CG-DNA/pair_oxdna_coaxstk.cpp
+++ b/src/CG-DNA/pair_oxdna_coaxstk.cpp
@@ -106,7 +106,6 @@ PairOxdnaCoaxstk::~PairOxdnaCoaxstk()
 
 void PairOxdnaCoaxstk::compute(int eflag, int vflag)
 {
-
   double delf[3],delt[3],delta[3],deltb[3]; // force, torque increment;
   double evdwl,fpair,finc,tpair,factor_lj;
   double v1tmp[3],v2tmp[3],v3tmp[3];
@@ -129,10 +128,9 @@ void PairOxdnaCoaxstk::compute(int eflag, int vflag)
   // vectors COM-backbone site, COM-stacking site in lab frame
   double ra_cs[3],ra_cst[3];
   double rb_cs[3],rb_cst[3];
-
-  // quaternions and Cartesian unit vectors in lab frame
-  double *qa,ax[3],ay[3],az[3];
-  double *qb,bx[3],by[3],bz[3];
+  // Cartesian unit vectors in lab frame
+  double ax[3],ay[3],az[3];
+  double bx[3],by[3],bz[3];
 
   double **x = atom->x;
   double **f = atom->f;
@@ -161,6 +159,12 @@ void PairOxdnaCoaxstk::compute(int eflag, int vflag)
   numneigh = list->numneigh;
   firstneigh = list->firstneigh;
 
+  // n(x/y/z)_xtrct = extracted local unit vectors in lab frame from oxdna_excv
+  int dim;
+  nx_xtrct = (double **) force->pair->extract("nx",dim);
+  ny_xtrct = (double **) force->pair->extract("ny",dim);
+  nz_xtrct = (double **) force->pair->extract("nz",dim);
+
   // loop over pair interaction neighbors of my atoms
 
   for (ia = 0; ia < anum; ia++) {
@@ -168,8 +172,10 @@ void PairOxdnaCoaxstk::compute(int eflag, int vflag)
     a = alist[ia];
     atype = type[a];
 
-    qa=bonus[ellipsoid[a]].quat;
-    MathExtra::q_to_exyz(qa,ax,ay,az);
+    ax[0] = nx_xtrct[a][0];
+    ax[1] = nx_xtrct[a][1];
+    ax[2] = nx_xtrct[a][2];
+    // a(y/z) not needed here as oxDNA(1) co-linear
 
     // vector COM a - stacking site a
     ra_cst[0] = d_cst*ax[0];
@@ -192,8 +198,10 @@ void PairOxdnaCoaxstk::compute(int eflag, int vflag)
 
       btype = type[b];
 
-      qb=bonus[ellipsoid[b]].quat;
-      MathExtra::q_to_exyz(qb,bx,by,bz);
+      bx[0] = nx_xtrct[b][0];
+      bx[1] = nx_xtrct[b][1];
+      bx[2] = nx_xtrct[b][2];
+      // b(y/z) not needed here as oxDNA(1) co-linear
 
       // vector COM b - stacking site b
       rb_cst[0] = d_cst*bx[0];
@@ -245,6 +253,13 @@ void PairOxdnaCoaxstk::compute(int eflag, int vflag)
       // early rejection criterium
       if (f4t1) {
 
+      az[0] = nz_xtrct[a][0];
+      az[1] = nz_xtrct[a][1];
+      az[2] = nz_xtrct[a][2];
+      bz[0] = nz_xtrct[b][0];
+      bz[1] = nz_xtrct[b][1];
+      bz[2] = nz_xtrct[b][2];
+
       cost4 = MathExtra::dot3(az,bz);
       if (cost4 >  1.0) cost4 =  1.0;
       if (cost4 < -1.0) cost4 = -1.0;
@@ -380,6 +395,10 @@ void PairOxdnaCoaxstk::compute(int eflag, int vflag)
       // cosphi3 and cosphi4 (=cosphi3) force and virial
       if (cosphi3) {
 
+        ay[0] = ny_xtrct[a][0];
+        ay[1] = ny_xtrct[a][1];
+        ay[2] = ny_xtrct[a][2];
+
         finc   = -f2 * f4t1* f4t4 * f4t5 * f4t6 * 2.0 * f5c3 * df5c3 * factor_lj;
         fpair += finc;
 
diff --git a/src/CG-DNA/pair_oxdna_coaxstk.h b/src/CG-DNA/pair_oxdna_coaxstk.h
index 86538a432f..62ffb7f2c1 100644
--- a/src/CG-DNA/pair_oxdna_coaxstk.h
+++ b/src/CG-DNA/pair_oxdna_coaxstk.h
@@ -47,21 +47,17 @@ class PairOxdnaCoaxstk : public Pair {
   double **k_cxst, **cut_cxst_0, **cut_cxst_c, **cut_cxst_lo, **cut_cxst_hi;
   double **cut_cxst_lc, **cut_cxst_hc, **b_cxst_lo, **b_cxst_hi;
   double **cutsq_cxst_hc;
-
   double **a_cxst1, **theta_cxst1_0, **dtheta_cxst1_ast;
   double **b_cxst1, **dtheta_cxst1_c;
-
   double **a_cxst4, **theta_cxst4_0, **dtheta_cxst4_ast;
   double **b_cxst4, **dtheta_cxst4_c;
-
   double **a_cxst5, **theta_cxst5_0, **dtheta_cxst5_ast;
   double **b_cxst5, **dtheta_cxst5_c;
-
   double **a_cxst6, **theta_cxst6_0, **dtheta_cxst6_ast;
   double **b_cxst6, **dtheta_cxst6_c;
-
   double **a_cxst3p, **cosphi_cxst3p_ast, **b_cxst3p, **cosphi_cxst3p_c;
   double **a_cxst4p, **cosphi_cxst4p_ast, **b_cxst4p, **cosphi_cxst4p_c;
+  double **nx_xtrct, **ny_xtrct, **nz_xtrct; // per-atom arrays for local unit vectors
 
   virtual void allocate();
 };
diff --git a/src/CG-DNA/pair_oxdna_excv.cpp b/src/CG-DNA/pair_oxdna_excv.cpp
index d41976e86c..ae5752028a 100644
--- a/src/CG-DNA/pair_oxdna_excv.cpp
+++ b/src/CG-DNA/pair_oxdna_excv.cpp
@@ -39,6 +39,9 @@ PairOxdnaExcv::PairOxdnaExcv(LAMMPS *lmp) : Pair(lmp)
 {
   single_enable = 0;
   writedata = 1;
+
+  // set comm size needed by this Pair
+  comm_forward = 9;
 }
 
 /* ---------------------------------------------------------------------- */
@@ -47,6 +50,10 @@ PairOxdnaExcv::~PairOxdnaExcv()
 {
   if (allocated) {
 
+    memory->destroy(nx);
+    memory->destroy(ny);
+    memory->destroy(nz);
+
     memory->destroy(setflag);
     memory->destroy(cutsq);
 
@@ -108,7 +115,6 @@ void PairOxdnaExcv::compute_interaction_sites(double e1[3], double /*e2*/[3],
 
 void PairOxdnaExcv::compute(int eflag, int vflag)
 {
-
   double delf[3],delta[3],deltb[3]; // force, torque increment;
   double evdwl,fpair,factor_lj;
   double rtmp_s[3],rtmp_b[3];
@@ -118,10 +124,10 @@ void PairOxdnaExcv::compute(int eflag, int vflag)
   // vectors COM-backbone site, COM-base site in lab frame
   double ra_cs[3],ra_cb[3];
   double rb_cs[3],rb_cb[3];
+  // Cartesian unit vectors in lab frame
+  double ax[3],ay[3],az[3];
+  double bx[3],by[3],bz[3];
 
-  // quaternions and Cartesian unit vectors in lab frame
-  double *qa,ax[3],ay[3],az[3];
-  double *qb,bx[3],by[3],bz[3];
   double *special_lj = force->special_lj;
 
   double **x = atom->x;
@@ -137,7 +143,7 @@ void PairOxdnaExcv::compute(int eflag, int vflag)
   AtomVecEllipsoid::Bonus *bonus = avec->bonus;
   int *ellipsoid = atom->ellipsoid;
 
-  int a,b,ia,ib,anum,bnum,atype,btype;
+  int n,a,b,in,ia,ib,anum,bnum,atype,btype;
 
   evdwl = 0.0;
   ev_init(eflag,vflag);
@@ -147,6 +153,29 @@ void PairOxdnaExcv::compute(int eflag, int vflag)
   numneigh = list->numneigh;
   firstneigh = list->firstneigh;
 
+  // loop over all local atoms, calculation of local reference frame
+  for (in = 0; in < atom->nlocal; in++) {
+
+    int n = alist[in];
+    double *qn,nx_temp[3],ny_temp[3],nz_temp[3]; // quaternion and Cartesian unit vectors in lab frame
+
+    qn=bonus[ellipsoid[n]].quat;
+    MathExtra::q_to_exyz(qn,nx_temp,ny_temp,nz_temp);
+
+    nx[n][0] = nx_temp[0];
+    nx[n][1] = nx_temp[1];
+    nx[n][2] = nx_temp[2];
+    ny[n][0] = ny_temp[0];
+    ny[n][1] = ny_temp[1];
+    ny[n][2] = ny_temp[2];
+    nz[n][0] = nz_temp[0];
+    nz[n][1] = nz_temp[1];
+    nz[n][2] = nz_temp[2];
+
+  }
+
+  comm->forward_comm_pair(this);
+
   // loop over pair interaction neighbors of my atoms
 
   for (ia = 0; ia < anum; ia++) {
@@ -154,8 +183,15 @@ void PairOxdnaExcv::compute(int eflag, int vflag)
     a = alist[ia];
     atype = type[a];
 
-    qa=bonus[ellipsoid[a]].quat;
-    MathExtra::q_to_exyz(qa,ax,ay,az);
+    ax[0] = nx[a][0];
+    ax[1] = nx[a][1];
+    ax[2] = nx[a][2];
+    ay[0] = ny[a][0];
+    ay[1] = ny[a][1];
+    ay[2] = ny[a][2];
+    az[0] = nz[a][0];
+    az[1] = nz[a][1];
+    az[2] = nz[a][2];
 
     // vector COM - backbone and base site a
     compute_interaction_sites(ax,ay,az,ra_cs,ra_cb);
@@ -179,8 +215,15 @@ void PairOxdnaExcv::compute(int eflag, int vflag)
 
       btype = type[b];
 
-      qb=bonus[ellipsoid[b]].quat;
-      MathExtra::q_to_exyz(qb,bx,by,bz);
+      bx[0] = nx[b][0];
+      bx[1] = nx[b][1];
+      bx[2] = nx[b][2];
+      by[0] = ny[b][0];
+      by[1] = ny[b][1];
+      by[2] = ny[b][2];
+      bz[0] = nz[b][0];
+      bz[1] = nz[b][1];
+      bz[2] = nz[b][2];
 
       // vector COM - backbone and base site b
       compute_interaction_sites(bx,by,bz,rb_cs,rb_cb);
@@ -400,6 +443,10 @@ void PairOxdnaExcv::allocate()
     for (int j = i; j <= n; j++)
       setflag[i][j] = 0;
 
+  memory->create(nx,atom->nmax,3,"pair:nx");
+  memory->create(ny,atom->nmax,3,"pair:ny");
+  memory->create(nz,atom->nmax,3,"pair:nz");
+
   memory->create(cutsq,n+1,n+1,"pair:cutsq");
 
   memory->create(epsilon_ss,n+1,n+1,"pair:epsilon_ss");
@@ -806,10 +853,58 @@ void PairOxdnaExcv::write_data_all(FILE *fp)
 
 /* ---------------------------------------------------------------------- */
 
+int PairOxdnaExcv::pack_forward_comm(int n, int *list, double *buf,
+                               int /*pbc_flag*/, int * /*pbc*/)
+{
+  int i,j,m;
+
+  m = 0;
+  for (i = 0; i < n; i++) {
+    j = list[i];
+	buf[m++] = nx[j][0];
+	buf[m++] = nx[j][1];
+	buf[m++] = nx[j][2];
+	buf[m++] = ny[j][0];
+	buf[m++] = ny[j][1];
+	buf[m++] = ny[j][2];
+	buf[m++] = nz[j][0];
+	buf[m++] = nz[j][1];
+	buf[m++] = nz[j][2];
+  }
+  return m;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void PairOxdnaExcv::unpack_forward_comm(int n, int first, double *buf)
+{
+  int i,m,last;
+
+  m = 0;
+  last = first + n;
+  for (i = first; i < last; i++) {
+	nx[i][0] = buf[m++];
+	nx[i][1] = buf[m++];
+	nx[i][2] = buf[m++];
+	ny[i][0] = buf[m++];
+	ny[i][1] = buf[m++];
+	ny[i][2] = buf[m++];
+	nz[i][0] = buf[m++];
+	nz[i][1] = buf[m++];
+	nz[i][2] = buf[m++];
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
 void *PairOxdnaExcv::extract(const char *str, int &dim)
 {
   dim = 2;
 
+  if (strcmp(str,"nx") == 0) return (void *) nx;
+  if (strcmp(str,"ny") == 0) return (void *) ny;
+  if (strcmp(str,"nz") == 0) return (void *) nz;
+
   if (strcmp(str,"epsilon_ss") == 0) return (void *) epsilon_ss;
   if (strcmp(str,"sigma_ss") == 0) return (void *) sigma_ss;
   if (strcmp(str,"cut_ss_ast") == 0) return (void *) cut_ss_ast;
diff --git a/src/CG-DNA/pair_oxdna_excv.h b/src/CG-DNA/pair_oxdna_excv.h
index 2e92985e0c..5cf781feaa 100644
--- a/src/CG-DNA/pair_oxdna_excv.h
+++ b/src/CG-DNA/pair_oxdna_excv.h
@@ -41,6 +41,8 @@ class PairOxdnaExcv : public Pair {
   void write_data(FILE *) override;
   void write_data_all(FILE *) override;
   void *extract(const char *, int &) override;
+  virtual int pack_forward_comm(int, int *, double *, int, int *);
+  virtual void unpack_forward_comm(int, int, double *);
 
  protected:
   // s=sugar-phosphate backbone site, b=base site, st=stacking site
@@ -52,6 +54,7 @@ class PairOxdnaExcv : public Pair {
   double **lj1_sb, **lj2_sb, **b_sb, **cut_sb_c, **cutsq_sb_c;
   double **epsilon_bb, **sigma_bb, **cut_bb_ast, **cutsq_bb_ast;
   double **lj1_bb, **lj2_bb, **b_bb, **cut_bb_c, **cutsq_bb_c;
+  double **nx, **ny, **nz; // per-atom arrays for local unit vectors
 
   virtual void allocate();
 };
diff --git a/src/CG-DNA/pair_oxdna_hbond.cpp b/src/CG-DNA/pair_oxdna_hbond.cpp
index 4bf0bc9c56..3d2e803aa8 100644
--- a/src/CG-DNA/pair_oxdna_hbond.cpp
+++ b/src/CG-DNA/pair_oxdna_hbond.cpp
@@ -148,10 +148,9 @@ void PairOxdnaHbond::compute(int eflag, int vflag)
   double d_chb=+0.4;
   // vectors COM-h-bonding site in lab frame
   double ra_chb[3],rb_chb[3];
-
-  // quaternions and Cartesian unit vectors in lab frame
-  double *qa,ax[3],ay[3],az[3];
-  double *qb,bx[3],by[3],bz[3];
+  // Cartesian unit vectors in lab frame
+  double ax[3],ay[3],az[3];
+  double bx[3],by[3],bz[3];
 
   double **x = atom->x;
   double **f = atom->f;
@@ -180,6 +179,12 @@ void PairOxdnaHbond::compute(int eflag, int vflag)
   numneigh = list->numneigh;
   firstneigh = list->firstneigh;
 
+  // n(x/y/z)_xtrct = extracted local unit vectors from oxdna_excv
+  int dim;
+  nx_xtrct = (double **) force->pair->extract("nx",dim);
+  ny_xtrct = (double **) force->pair->extract("ny",dim);
+  nz_xtrct = (double **) force->pair->extract("nz",dim);
+
   // loop over pair interaction neighbors of my atoms
 
   for (ia = 0; ia < anum; ia++) {
@@ -187,8 +192,9 @@ void PairOxdnaHbond::compute(int eflag, int vflag)
     a = alist[ia];
     atype = type[a];
 
-    qa=bonus[ellipsoid[a]].quat;
-    MathExtra::q_to_exyz(qa,ax,ay,az);
+    ax[0] = nx_xtrct[a][0];
+    ax[1] = nx_xtrct[a][1];
+    ax[2] = nx_xtrct[a][2];
 
     ra_chb[0] = d_chb*ax[0];
     ra_chb[1] = d_chb*ax[1];
@@ -205,8 +211,9 @@ void PairOxdnaHbond::compute(int eflag, int vflag)
 
       btype = type[b];
 
-      qb=bonus[ellipsoid[b]].quat;
-      MathExtra::q_to_exyz(qb,bx,by,bz);
+      bx[0] = nx_xtrct[b][0];
+      bx[1] = nx_xtrct[b][1];
+      bx[2] = nx_xtrct[b][2];
 
       rb_chb[0] = d_chb*bx[0];
       rb_chb[1] = d_chb*bx[1];
@@ -265,6 +272,13 @@ void PairOxdnaHbond::compute(int eflag, int vflag)
       // early rejection criterium
       if (f4t3) {
 
+      az[0] = nz_xtrct[a][0];
+      az[1] = nz_xtrct[a][1];
+      az[2] = nz_xtrct[a][2];
+      bz[0] = nz_xtrct[b][0];
+      bz[1] = nz_xtrct[b][1];
+      bz[2] = nz_xtrct[b][2];
+
       cost4 = MathExtra::dot3(az,bz);
       if (cost4 >  1.0) cost4 =  1.0;
       if (cost4 < -1.0) cost4 = -1.0;
diff --git a/src/CG-DNA/pair_oxdna_hbond.h b/src/CG-DNA/pair_oxdna_hbond.h
index 7c6add7c4a..32120a79cf 100644
--- a/src/CG-DNA/pair_oxdna_hbond.h
+++ b/src/CG-DNA/pair_oxdna_hbond.h
@@ -48,25 +48,19 @@ class PairOxdnaHbond : public Pair {
   double **epsilon_hb, **a_hb, **cut_hb_0, **cut_hb_c, **cut_hb_lo, **cut_hb_hi;
   double **cut_hb_lc, **cut_hb_hc, **b_hb_lo, **b_hb_hi, **shift_hb;
   double **cutsq_hb_hc;
-
   double **a_hb1, **theta_hb1_0, **dtheta_hb1_ast;
   double **b_hb1, **dtheta_hb1_c;
-
   double **a_hb2, **theta_hb2_0, **dtheta_hb2_ast;
   double **b_hb2, **dtheta_hb2_c;
-
   double **a_hb3, **theta_hb3_0, **dtheta_hb3_ast;
   double **b_hb3, **dtheta_hb3_c;
-
   double **a_hb4, **theta_hb4_0, **dtheta_hb4_ast;
   double **b_hb4, **dtheta_hb4_c;
-
   double **a_hb7, **theta_hb7_0, **dtheta_hb7_ast;
   double **b_hb7, **dtheta_hb7_c;
-
   double **a_hb8, **theta_hb8_0, **dtheta_hb8_ast;
   double **b_hb8, **dtheta_hb8_c;
-
+  double **nx_xtrct, **ny_xtrct, **nz_xtrct; // per-atom arrays for local unit vectors
   int seqdepflag;
 
   virtual void allocate();
diff --git a/src/CG-DNA/pair_oxdna_stk.cpp b/src/CG-DNA/pair_oxdna_stk.cpp
index 4980408986..feeb41d78b 100644
--- a/src/CG-DNA/pair_oxdna_stk.cpp
+++ b/src/CG-DNA/pair_oxdna_stk.cpp
@@ -212,7 +212,6 @@ void PairOxdnaStk::ev_tally_xyz(int i, int j, int nlocal, int newton_bond,
 
 void PairOxdnaStk::compute(int eflag, int vflag)
 {
-
   double delf[3],delta[3],deltb[3]; // force, torque increment;
   double evdwl,fpair,finc,tpair;
   double delr_ss[3],delr_ss_norm[3],rsq_ss,r_ss,rinv_ss;
@@ -227,10 +226,9 @@ void PairOxdnaStk::compute(int eflag, int vflag)
   // vectors COM-backbone site, COM-stacking site in lab frame
   double ra_cs[3],ra_cst[3];
   double rb_cs[3],rb_cst[3];
-
-  // quaternions and Cartesian unit vectors in lab frame
-  double *qa,ax[3],ay[3],az[3];
-  double *qb,bx[3],by[3],bz[3];
+  // Cartesian unit vectors in lab frame
+  double ax[3],ay[3],az[3];
+  double bx[3],by[3],bz[3];
 
   double **x = atom->x;
   double **f = atom->f;
@@ -257,6 +255,12 @@ void PairOxdnaStk::compute(int eflag, int vflag)
   evdwl = 0.0;
   ev_init(eflag,vflag);
 
+  // n(x/y/z)_xtrct = extracted local unit vectors from oxdna_excv
+  int dim;
+  nx_xtrct = (double **) force->pair->extract("nx",dim);
+  ny_xtrct = (double **) force->pair->extract("ny",dim);
+  nz_xtrct = (double **) force->pair->extract("nz",dim);
+
   // loop over stacking interaction neighbors using bond topology
 
   for (in = 0; in < nbondlist; in++) {
@@ -275,10 +279,13 @@ void PairOxdnaStk::compute(int eflag, int vflag)
 
     // a now in 3' direction, b in 5' direction
 
-    qa=bonus[ellipsoid[a]].quat;
-    MathExtra::q_to_exyz(qa,ax,ay,az);
-    qb=bonus[ellipsoid[b]].quat;
-    MathExtra::q_to_exyz(qb,bx,by,bz);
+    ax[0] = nx_xtrct[a][0];
+    ax[1] = nx_xtrct[a][1];
+    ax[2] = nx_xtrct[a][2];
+    bx[0] = nx_xtrct[b][0];
+    bx[1] = nx_xtrct[b][1];
+    bx[2] = nx_xtrct[b][2];
+    // (a/b)y/z not needed here as oxDNA(1) co-linear
 
     // vector COM a - stacking site a
     ra_cst[0] = d_cst*ax[0];
@@ -336,6 +343,13 @@ void PairOxdnaStk::compute(int eflag, int vflag)
     // early rejection criterium
     if (f1) {
 
+    az[0] = nz_xtrct[a][0];
+    az[1] = nz_xtrct[a][1];
+    az[2] = nz_xtrct[a][2];
+    bz[0] = nz_xtrct[b][0];
+    bz[1] = nz_xtrct[b][1];
+    bz[2] = nz_xtrct[b][2];
+
     // theta4 angle and correction
     cost4 = MathExtra::dot3(bz,az);
     if (cost4 >  1.0) cost4 =  1.0;
@@ -360,6 +374,13 @@ void PairOxdnaStk::compute(int eflag, int vflag)
     // early rejection criterium
     if (f4t5) {
 
+    ay[0] = ny_xtrct[a][0];
+    ay[1] = ny_xtrct[a][1];
+    ay[2] = ny_xtrct[a][2];
+    by[0] = ny_xtrct[b][0];
+    by[1] = ny_xtrct[b][1];
+    by[2] = ny_xtrct[b][2];
+
     cost6p = MathExtra::dot3(delr_st_norm,az);
     if (cost6p >  1.0) cost6p =  1.0;
     if (cost6p < -1.0) cost6p = -1.0;
diff --git a/src/CG-DNA/pair_oxdna_stk.h b/src/CG-DNA/pair_oxdna_stk.h
index 1e97f8823c..5e9329ea57 100644
--- a/src/CG-DNA/pair_oxdna_stk.h
+++ b/src/CG-DNA/pair_oxdna_stk.h
@@ -59,7 +59,7 @@ class PairOxdnaStk : public Pair {
   double **b_st6, **dtheta_st6_c;
   double **a_st1, **cosphi_st1_ast, **b_st1, **cosphi_st1_c;
   double **a_st2, **cosphi_st2_ast, **b_st2, **cosphi_st2_c;
-
+  double **nx_xtrct, **ny_xtrct, **nz_xtrct; // per-atom arrays for local unit vectors
   int seqdepflag;
 
   virtual void allocate();
diff --git a/src/CG-DNA/pair_oxdna_xstk.cpp b/src/CG-DNA/pair_oxdna_xstk.cpp
index 447483e413..61dc1a8100 100644
--- a/src/CG-DNA/pair_oxdna_xstk.cpp
+++ b/src/CG-DNA/pair_oxdna_xstk.cpp
@@ -111,7 +111,6 @@ PairOxdnaXstk::~PairOxdnaXstk()
 
 void PairOxdnaXstk::compute(int eflag, int vflag)
 {
-
   double delf[3],delta[3],deltb[3]; // force, torque increment;
   double evdwl,fpair,finc,tpair,factor_lj;
   double delr_hb[3],delr_hb_norm[3],rsq_hb,r_hb,rinv_hb;
@@ -126,10 +125,9 @@ void PairOxdnaXstk::compute(int eflag, int vflag)
   double d_chb=+0.4;
   // vectors COM-h-bonding site in lab frame
   double ra_chb[3],rb_chb[3];
-
-  // quaternions and Cartesian unit vectors in lab frame
-  double *qa,ax[3],ay[3],az[3];
-  double *qb,bx[3],by[3],bz[3];
+  // Cartesian unit vectors in lab frame
+  double ax[3],ay[3],az[3];
+  double bx[3],by[3],bz[3];
 
   double **x = atom->x;
   double **f = atom->f;
@@ -158,6 +156,12 @@ void PairOxdnaXstk::compute(int eflag, int vflag)
   numneigh = list->numneigh;
   firstneigh = list->firstneigh;
 
+  // n(x/y/z)_xtrct = extracted local unit vectors from oxdna_excv
+  int dim;
+  nx_xtrct = (double **) force->pair->extract("nx",dim);
+  ny_xtrct = (double **) force->pair->extract("ny",dim);
+  nz_xtrct = (double **) force->pair->extract("nz",dim);
+
   // loop over pair interaction neighbors of my atoms
 
   for (ia = 0; ia < anum; ia++) {
@@ -165,8 +169,10 @@ void PairOxdnaXstk::compute(int eflag, int vflag)
     a = alist[ia];
     atype = type[a];
 
-    qa=bonus[ellipsoid[a]].quat;
-    MathExtra::q_to_exyz(qa,ax,ay,az);
+    ax[0] = nx_xtrct[a][0];
+    ax[1] = nx_xtrct[a][1];
+    ax[2] = nx_xtrct[a][2];
+    // a(y/z) not needed here as oxDNA(1) co-linear
 
     ra_chb[0] = d_chb*ax[0];
     ra_chb[1] = d_chb*ax[1];
@@ -183,8 +189,10 @@ void PairOxdnaXstk::compute(int eflag, int vflag)
 
       btype = type[b];
 
-      qb=bonus[ellipsoid[b]].quat;
-      MathExtra::q_to_exyz(qb,bx,by,bz);
+      bx[0] = nx_xtrct[b][0];
+      bx[1] = nx_xtrct[b][1];
+      bx[2] = nx_xtrct[b][2];
+      // b(y/z) not needed here as oxDNA(1) co-linear
 
       rb_chb[0] = d_chb*bx[0];
       rb_chb[1] = d_chb*bx[1];
@@ -243,6 +251,13 @@ void PairOxdnaXstk::compute(int eflag, int vflag)
       // early rejection criterium
       if (f4t3) {
 
+      az[0] = nz_xtrct[a][0];
+      az[1] = nz_xtrct[a][1];
+      az[2] = nz_xtrct[a][2];
+      bz[0] = nz_xtrct[b][0];
+      bz[1] = nz_xtrct[b][1];
+      bz[2] = nz_xtrct[b][2];
+
       cost4 = MathExtra::dot3(az,bz);
       if (cost4 >  1.0) cost4 =  1.0;
       if (cost4 < -1.0) cost4 = -1.0;
diff --git a/src/CG-DNA/pair_oxdna_xstk.h b/src/CG-DNA/pair_oxdna_xstk.h
index 28fc7fb2f6..8d0a126943 100644
--- a/src/CG-DNA/pair_oxdna_xstk.h
+++ b/src/CG-DNA/pair_oxdna_xstk.h
@@ -47,24 +47,19 @@ class PairOxdnaXstk : public Pair {
   double **k_xst, **cut_xst_0, **cut_xst_c, **cut_xst_lo, **cut_xst_hi;
   double **cut_xst_lc, **cut_xst_hc, **b_xst_lo, **b_xst_hi;
   double **cutsq_xst_hc;
-
   double **a_xst1, **theta_xst1_0, **dtheta_xst1_ast;
   double **b_xst1, **dtheta_xst1_c;
-
   double **a_xst2, **theta_xst2_0, **dtheta_xst2_ast;
   double **b_xst2, **dtheta_xst2_c;
-
   double **a_xst3, **theta_xst3_0, **dtheta_xst3_ast;
   double **b_xst3, **dtheta_xst3_c;
-
   double **a_xst4, **theta_xst4_0, **dtheta_xst4_ast;
   double **b_xst4, **dtheta_xst4_c;
-
   double **a_xst7, **theta_xst7_0, **dtheta_xst7_ast;
   double **b_xst7, **dtheta_xst7_c;
-
   double **a_xst8, **theta_xst8_0, **dtheta_xst8_ast;
   double **b_xst8, **dtheta_xst8_c;
+  double **nx_xtrct, **ny_xtrct, **nz_xtrct; // per-atom arrays for local unit vectors
 
   virtual void allocate();
 };
diff --git a/src/CG-DNA/pair_oxrna2_stk.cpp b/src/CG-DNA/pair_oxrna2_stk.cpp
index 18cd798fc7..c5fa47e25e 100644
--- a/src/CG-DNA/pair_oxrna2_stk.cpp
+++ b/src/CG-DNA/pair_oxrna2_stk.cpp
@@ -245,9 +245,9 @@ void PairOxrna2Stk::compute(int eflag, int vflag)
   double ra_cs[3],ra_cst[3];
   double rb_cs[3],rb_cst[3];
 
-  // quaternions and Cartesian unit vectors in lab frame
-  double *qa,ax[3],ay[3],az[3];
-  double *qb,bx[3],by[3],bz[3];
+  // Cartesian unit vectors in lab frame
+  double ax[3],ay[3],az[3];
+  double bx[3],by[3],bz[3];
 
   double **x = atom->x;
   double **f = atom->f;
@@ -274,6 +274,12 @@ void PairOxrna2Stk::compute(int eflag, int vflag)
   evdwl = 0.0;
   ev_init(eflag,vflag);
 
+  // n(x/y/z)_xtrct = extracted local unit vectors from oxdna_excv
+  int dim;
+  nx_xtrct = (double **) force->pair->extract("nx",dim);
+  ny_xtrct = (double **) force->pair->extract("ny",dim);
+  nz_xtrct = (double **) force->pair->extract("nz",dim);
+
   // loop over stacking interaction neighbors using bond topology
 
   for (in = 0; in < nbondlist; in++) {
@@ -290,12 +296,26 @@ void PairOxrna2Stk::compute(int eflag, int vflag)
 
     }
 
-   // a now in 3' direction, b in 5' direction
+    // a now in 3' direction, b in 5' direction
 
-    qa=bonus[ellipsoid[a]].quat;
-    MathExtra::q_to_exyz(qa,ax,ay,az);
-    qb=bonus[ellipsoid[b]].quat;
-    MathExtra::q_to_exyz(qb,bx,by,bz);
+    ax[0] = nx_xtrct[a][0];
+    ax[1] = nx_xtrct[a][1];
+    ax[2] = nx_xtrct[a][2];
+    ay[0] = ny_xtrct[a][0];
+    ay[1] = ny_xtrct[a][1];
+    ay[2] = ny_xtrct[a][2];
+    az[0] = nz_xtrct[a][0];
+    az[1] = nz_xtrct[a][1];
+    az[2] = nz_xtrct[a][2];
+    bx[0] = nx_xtrct[b][0];
+    bx[1] = nx_xtrct[b][1];
+    bx[2] = nx_xtrct[b][2];
+    by[0] = ny_xtrct[b][0];
+    by[1] = ny_xtrct[b][1];
+    by[2] = ny_xtrct[b][2];
+    bz[0] = nz_xtrct[b][0];
+    bz[1] = nz_xtrct[b][1];
+    bz[2] = nz_xtrct[b][2];
 
     // vector COM a - 5'-stacking site a
     ra_cst[0] = d_cst_x_5p*ax[0] + d_cst_y_5p*ay[0];
diff --git a/src/CG-DNA/pair_oxrna2_stk.h b/src/CG-DNA/pair_oxrna2_stk.h
index b8b2102a7a..8530544fdc 100644
--- a/src/CG-DNA/pair_oxrna2_stk.h
+++ b/src/CG-DNA/pair_oxrna2_stk.h
@@ -60,7 +60,7 @@ class PairOxrna2Stk : public Pair {
   double **b_st10, **dtheta_st10_c;
   double **a_st1, **cosphi_st1_ast, **b_st1, **cosphi_st1_c;
   double **a_st2, **cosphi_st2_ast, **b_st2, **cosphi_st2_c;
-
+  double **nx_xtrct, **ny_xtrct, **nz_xtrct; // per-atom arrays for local unit vectors
   int seqdepflag;
 
   virtual void allocate();
diff --git a/src/CG-DNA/pair_oxrna2_xstk.cpp b/src/CG-DNA/pair_oxrna2_xstk.cpp
index 8b44a9bde1..1f069e6b2f 100644
--- a/src/CG-DNA/pair_oxrna2_xstk.cpp
+++ b/src/CG-DNA/pair_oxrna2_xstk.cpp
@@ -120,9 +120,9 @@ void PairOxrna2Xstk::compute(int eflag, int vflag)
   // vectors COM-h-bonding site in lab frame
   double ra_chb[3],rb_chb[3];
 
-  // quaternions and Cartesian unit vectors in lab frame
-  double *qa,ax[3],ay[3],az[3];
-  double *qb,bx[3],by[3],bz[3];
+  // Cartesian unit vectors in lab frame
+  double ax[3],ay[3],az[3];
+  double bx[3],by[3],bz[3];
 
   double **x = atom->x;
   double **f = atom->f;
@@ -151,6 +151,11 @@ void PairOxrna2Xstk::compute(int eflag, int vflag)
   numneigh = list->numneigh;
   firstneigh = list->firstneigh;
 
+  // n(x/z)_xtrct = extracted local unit vectors from oxdna_excv
+  int dim;
+  nx_xtrct = (double **) force->pair->extract("nx",dim);
+  nz_xtrct = (double **) force->pair->extract("nz",dim);
+
   // loop over pair interaction neighbors of my atoms
 
   for (ia = 0; ia < anum; ia++) {
@@ -158,8 +163,9 @@ void PairOxrna2Xstk::compute(int eflag, int vflag)
     a = alist[ia];
     atype = type[a];
 
-    qa=bonus[ellipsoid[a]].quat;
-    MathExtra::q_to_exyz(qa,ax,ay,az);
+    ax[0] = nx_xtrct[a][0];
+    ax[1] = nx_xtrct[a][1];
+    ax[2] = nx_xtrct[a][2];
 
     ra_chb[0] = d_chb*ax[0];
     ra_chb[1] = d_chb*ax[1];
@@ -176,8 +182,9 @@ void PairOxrna2Xstk::compute(int eflag, int vflag)
 
       btype = type[b];
 
-      qb=bonus[ellipsoid[b]].quat;
-      MathExtra::q_to_exyz(qb,bx,by,bz);
+      bx[0] = nx_xtrct[b][0];
+      bx[1] = nx_xtrct[b][1];
+      bx[2] = nx_xtrct[b][2];
 
       rb_chb[0] = d_chb*bx[0];
       rb_chb[1] = d_chb*bx[1];
@@ -236,6 +243,10 @@ void PairOxrna2Xstk::compute(int eflag, int vflag)
       // early rejection criterium
       if (f4t3) {
 
+      az[0] = nz_xtrct[a][0];
+      az[1] = nz_xtrct[a][1];
+      az[2] = nz_xtrct[a][2];
+
       cost7 = -1.0*MathExtra::dot3(az,delr_hb_norm);
       if (cost7 >  1.0) cost7 =  1.0;
       if (cost7 < -1.0) cost7 = -1.0;
@@ -250,6 +261,10 @@ void PairOxrna2Xstk::compute(int eflag, int vflag)
       // early rejection criterium
       if (f4t7) {
 
+      bz[0] = nz_xtrct[b][0];
+      bz[1] = nz_xtrct[b][1];
+      bz[2] = nz_xtrct[b][2];
+
       cost8 = MathExtra::dot3(bz,delr_hb_norm);
       if (cost8 >  1.0) cost8 =  1.0;
       if (cost8 < -1.0) cost8 = -1.0;
diff --git a/src/CG-DNA/pair_oxrna2_xstk.h b/src/CG-DNA/pair_oxrna2_xstk.h
index 6a3f02ee8d..c73839b4a7 100644
--- a/src/CG-DNA/pair_oxrna2_xstk.h
+++ b/src/CG-DNA/pair_oxrna2_xstk.h
@@ -46,21 +46,17 @@ class PairOxrna2Xstk : public Pair {
   double **k_xst, **cut_xst_0, **cut_xst_c, **cut_xst_lo, **cut_xst_hi;
   double **cut_xst_lc, **cut_xst_hc, **b_xst_lo, **b_xst_hi;
   double **cutsq_xst_hc;
-
   double **a_xst1, **theta_xst1_0, **dtheta_xst1_ast;
   double **b_xst1, **dtheta_xst1_c;
-
   double **a_xst2, **theta_xst2_0, **dtheta_xst2_ast;
   double **b_xst2, **dtheta_xst2_c;
-
   double **a_xst3, **theta_xst3_0, **dtheta_xst3_ast;
   double **b_xst3, **dtheta_xst3_c;
-
   double **a_xst7, **theta_xst7_0, **dtheta_xst7_ast;
   double **b_xst7, **dtheta_xst7_c;
-
   double **a_xst8, **theta_xst8_0, **dtheta_xst8_ast;
   double **b_xst8, **dtheta_xst8_c;
+  double **nx_xtrct, **nz_xtrct; // per-atom arrays for local unit vectors
 
   virtual void allocate();
 };
diff --git a/unittest/formats/test_atom_styles.cpp b/unittest/formats/test_atom_styles.cpp
index ab00efc949..c4c79238aa 100644
--- a/unittest/formats/test_atom_styles.cpp
+++ b/unittest/formats/test_atom_styles.cpp
@@ -4809,6 +4809,446 @@ TEST_F(AtomStyleTest, property_atom)
     EXPECT_NEAR(three[GETIDX(2)], 0.5, EPSILON);
 }
 
+TEST_F(AtomStyleTest, oxdna)
+{
+    if (!LAMMPS::is_installed_pkg("MOLECULE")) GTEST_SKIP();
+    if (!LAMMPS::is_installed_pkg("ASPHERE")) GTEST_SKIP();
+    if (!LAMMPS::is_installed_pkg("CG-DNA")) GTEST_SKIP();
+
+    BEGIN_HIDE_OUTPUT();
+    command("atom_style hybrid bond ellipsoid oxdna");
+    END_HIDE_OUTPUT();
+
+    AtomState expected;
+    expected.atom_style     = "hybrid";
+    expected.molecular      = Atom::MOLECULAR;
+    expected.tag_enable     = 1;
+    expected.molecule_flag  = 1;
+    expected.ellipsoid_flag = 1;
+    expected.rmass_flag     = 1;
+    expected.torque_flag    = 1;
+    expected.angmom_flag    = 1;
+    expected.has_type       = true;
+    expected.has_mask       = true;
+    expected.has_image      = true;
+    expected.has_x          = true;
+    expected.has_v          = true;
+    expected.has_f          = true;
+    expected.has_bonds      = true;
+    expected.has_nspecial   = true;
+    expected.has_special    = true;
+    expected.map_style      = 3;
+
+    ASSERT_ATOM_STATE_EQ(lmp->atom, expected);
+
+    auto hybrid = (AtomVecHybrid *)lmp->atom->avec;
+    ASSERT_THAT(std::string(lmp->atom->atom_style), Eq("hybrid"));
+    ASSERT_EQ(hybrid->nstyles, 3);
+    ASSERT_THAT(std::string(hybrid->keywords[0]), Eq("bond"));
+    ASSERT_THAT(std::string(hybrid->keywords[1]), Eq("ellipsoid"));
+    ASSERT_THAT(std::string(hybrid->keywords[2]), Eq("oxdna"));
+    ASSERT_NE(hybrid->styles[0], nullptr);
+    ASSERT_NE(hybrid->styles[1], nullptr);
+    ASSERT_NE(hybrid->styles[2], nullptr);
+
+    BEGIN_HIDE_OUTPUT();
+    command("units lj");
+    command("dimension 3");
+    command("newton on");
+    command("boundary  p p p");
+    command("atom_modify sort 0 1.0");
+    command("neighbor 2.0 bin");
+    command("neigh_modify every 1 delay 0 check yes");
+    command("region mybox block -20 20 -20 20 -20 20");
+    command("create_box 4 mybox bond/types 1 extra/bond/per/atom 2 extra/special/per/atom 4");
+    command("create_atoms 1 single -0.33741452300167507 -0.43708835412476305 0.6450685042019271");
+    command("create_atoms 2 single -0.32142606102826937 -0.7137743037592722 1.1817366147004618");
+    command("create_atoms 3 single -0.130363628207774 -0.9147144801536078 1.62581312195109");
+    command("create_atoms 4 single 0.16795127962282844 -0.9808507459807022 2.0894908590909003");
+    command("create_atoms 1 single 0.46370423490634166 -0.7803347954883079 2.4251986815515827");
+    command("create_atoms 4 single -0.4462950185476711 0.09062163051035639 2.4668941268777607");
+    command("create_atoms 1 single -0.03377054097560965 0.20979847489755046 2.078208732038921");
+    command("create_atoms 2 single 0.3297325391466579 0.17657587120899895 1.7206328374934152");
+    command("create_atoms 3 single 0.6063699309305985 0.04682595158675571 1.2335049647817748");
+    command("create_atoms 4 single 0.8003979559814726 -0.364393011459011 0.9884025318908612");
+    command("set type 1 mass 3.1575");
+    command("set type 2 mass 3.1575");
+    command("set type 3 mass 3.1575");
+    command("set type 4 mass 3.1575");
+    command("mass 1 3.1575");
+    command("mass 2 3.1575");
+    command("mass 3 3.1575");
+    command("mass 4 3.1575");
+    command("set atom 1 shape 1.173984503142341 1.173984503142341 1.173984503142341");
+    command("set atom 2 shape 1.173984503142341 1.173984503142341 1.173984503142341");
+    command("set atom 3 shape 1.173984503142341 1.173984503142341 1.173984503142341");
+    command("set atom 4 shape 1.173984503142341 1.173984503142341 1.173984503142341");
+    command("set atom 5 shape 1.173984503142341 1.173984503142341 1.173984503142341");
+    command("set atom 6 shape 1.173984503142341 1.173984503142341 1.173984503142341");
+    command("set atom 7 shape 1.173984503142341 1.173984503142341 1.173984503142341");
+    command("set atom 8 shape 1.173984503142341 1.173984503142341 1.173984503142341");
+    command("set atom 9 shape 1.173984503142341 1.173984503142341 1.173984503142341");
+    command("set atom 10 shape 1.173984503142341 1.173984503142341 1.173984503142341");
+    command("set atom 1 quat 0.120438343611269 -0.970540441176996 0.208676441957758 16.990727782866998");
+    command("set atom 2 quat 0.122039415796829 -0.068232256412985 0.990177125658213 40.001729435287870");
+    command("set atom 3 quat 0.052760168698289 0.030943512185297 0.998127679033382 69.627682451632380");
+    command("set atom 4 quat -0.037622918613871 0.030623545471522 0.998822664169035 97.038820280300570");
+    command("set atom 5 quat 0.055056042946138 0.077631917807377 0.995460756369964 137.7813218321917");
+    command("set atom 6 quat 0.931128471673637 -0.355724722922553 -0.080372201291206 166.2836226291888");
+    command("set atom 7 quat 0.753526078198930 -0.648440397941919 0.108275111595674 200.6802564250672");
+    command("set atom 8 quat 0.553942138074214 -0.829580511279186 0.070315595507185 192.0355407659524");
+    command("set atom 9 quat -0.373540155765431 0.913070802138105 -0.163613759548524 171.0789308260751");
+    command("set atom 10 quat 0.027515673832457 0.998248649922676 -0.052369080773879 161.2621224558284");
+    command("bond_style oxdna2/fene");
+    command("bond_coeff * 2.0 0.25 0.7564");
+    command("special_bonds lj 0 1 1");
+    command("create_bonds single/bond 1 1 2");
+    command("create_bonds single/bond 1 2 3");
+    command("create_bonds single/bond 1 3 4");
+    command("create_bonds single/bond 1 4 5");
+    command("create_bonds single/bond 1 6 7");
+    command("create_bonds single/bond 1 7 8");
+    command("create_bonds single/bond 1 8 9");
+    command("create_bonds single/bond 1 9 10");
+    command("pair_style hybrid/overlay oxdna2/excv oxdna2/stk oxdna2/hbond oxdna2/xstk oxdna2/coaxstk oxdna2/dh");
+    command("pair_coeff * * oxdna2/excv 2.0 0.7 0.675 2.0 0.515 0.5 2.0 0.33 0.32");
+    command("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");
+    command("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");
+    command("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");
+    command("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");
+    command("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");
+    command("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");
+    command("pair_coeff * * oxdna2/dh 0.1 0.2 0.815");
+    END_HIDE_OUTPUT();
+
+    ASSERT_THAT(std::string(lmp->atom->atom_style), Eq("hybrid"));
+    ASSERT_NE(lmp->atom->avec, nullptr);
+    ASSERT_EQ(lmp->atom->natoms, 10);
+    ASSERT_EQ(lmp->atom->nbonds, 8);
+    ASSERT_EQ(lmp->atom->nbondtypes, 1);
+    ASSERT_EQ(lmp->atom->nellipsoids, 10);
+    ASSERT_EQ(lmp->atom->nlocal, 10);
+    ASSERT_EQ(lmp->atom->nghost, 0);
+    ASSERT_NE(lmp->atom->nmax, -1);
+    ASSERT_EQ(lmp->atom->tag_enable, 1);
+    ASSERT_EQ(lmp->atom->molecular, Atom::MOLECULAR);
+    ASSERT_EQ(lmp->atom->ntypes, 4);
+    ASSERT_EQ(lmp->atom->nextra_grow, 0);
+    ASSERT_EQ(lmp->atom->nextra_restart, 0);
+    ASSERT_EQ(lmp->atom->nextra_border, 0);
+    ASSERT_EQ(lmp->atom->nextra_grow_max, 0);
+    ASSERT_EQ(lmp->atom->nextra_restart_max, 0);
+    ASSERT_EQ(lmp->atom->nextra_border_max, 0);
+    ASSERT_EQ(lmp->atom->nextra_store, 0);
+    ASSERT_EQ(lmp->atom->extra_grow, nullptr);
+    ASSERT_EQ(lmp->atom->extra_restart, nullptr);
+    ASSERT_EQ(lmp->atom->extra_border, nullptr);
+    ASSERT_EQ(lmp->atom->extra, nullptr);
+    ASSERT_NE(lmp->atom->mass, nullptr);
+    ASSERT_NE(lmp->atom->rmass, nullptr);
+    ASSERT_EQ(lmp->atom->ellipsoid_flag, 1);
+    ASSERT_NE(lmp->atom->ellipsoid, nullptr);
+    ASSERT_NE(lmp->atom->mass_setflag, nullptr);
+    ASSERT_NE(lmp->atom->id5p, nullptr);
+
+    BEGIN_HIDE_OUTPUT();
+    command("write_data test_atom_styles.data nocoeff");
+    command("clear");
+    command("units lj");
+    command("dimension 3");
+    command("newton on");
+    command("boundary  p p p");
+    command("atom_style hybrid bond ellipsoid oxdna");
+    command("atom_modify sort 0 1.0");
+    command("neighbor 2.0 bin");
+    command("neigh_modify every 1 delay 0 check yes");
+    command("read_data test_atom_styles.data");
+    command("set type 1 mass 3.1575");
+    command("set type 2 mass 3.1575");
+    command("set type 3 mass 3.1575");
+    command("set type 4 mass 3.1575");
+    command("mass 1 3.1575");
+    command("mass 2 3.1575");
+    command("mass 3 3.1575");
+    command("mass 4 3.1575");
+    command("bond_style oxdna2/fene");
+    command("bond_coeff * 2.0 0.25 0.7564");
+    command("special_bonds lj 0 1 1");
+    command("pair_style hybrid/overlay oxdna2/excv oxdna2/stk oxdna2/hbond oxdna2/xstk oxdna2/coaxstk oxdna2/dh");
+    command("pair_coeff * * oxdna2/excv 2.0 0.7 0.675 2.0 0.515 0.5 2.0 0.33 0.32");
+    command("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");
+    command("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");
+    command("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");
+    command("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");
+    command("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");
+    command("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");
+    command("pair_coeff * * oxdna2/dh 0.1 0.2 0.815");
+
+    ASSERT_THAT(std::string(lmp->atom->atom_style), Eq("hybrid"));
+    ASSERT_NE(lmp->atom->avec, nullptr);
+    hybrid = (AtomVecHybrid *)lmp->atom->avec;
+
+    ASSERT_EQ(hybrid->nstyles, 3);
+    ASSERT_THAT(std::string(hybrid->keywords[0]), Eq("bond"));
+    ASSERT_THAT(std::string(hybrid->keywords[1]), Eq("ellipsoid"));
+    ASSERT_THAT(std::string(hybrid->keywords[2]), Eq("oxdna"));
+    ASSERT_NE(hybrid->styles[0], nullptr);
+    ASSERT_NE(hybrid->styles[1], nullptr);
+    ASSERT_NE(hybrid->styles[2], nullptr);
+
+    ASSERT_EQ(lmp->atom->natoms, 10);
+    ASSERT_EQ(lmp->atom->nbonds, 8);
+    ASSERT_EQ(lmp->atom->nbondtypes, 1);
+    ASSERT_EQ(lmp->atom->nellipsoids, 10);
+    ASSERT_EQ(lmp->atom->nlocal, 10);
+    ASSERT_EQ(lmp->atom->nghost, 0);
+    ASSERT_NE(lmp->atom->nmax, -1);
+    ASSERT_EQ(lmp->atom->tag_enable, 1);
+    ASSERT_EQ(lmp->atom->molecular, Atom::MOLECULAR);
+    ASSERT_EQ(lmp->atom->ntypes, 4);
+    ASSERT_EQ(lmp->atom->nextra_grow, 0);
+    ASSERT_EQ(lmp->atom->nextra_restart, 0);
+    ASSERT_EQ(lmp->atom->nextra_border, 0);
+    ASSERT_EQ(lmp->atom->nextra_grow_max, 0);
+    ASSERT_EQ(lmp->atom->nextra_restart_max, 0);
+    ASSERT_EQ(lmp->atom->nextra_border_max, 0);
+    ASSERT_EQ(lmp->atom->nextra_store, 0);
+    ASSERT_EQ(lmp->atom->extra_grow, nullptr);
+    ASSERT_EQ(lmp->atom->extra_restart, nullptr);
+    ASSERT_EQ(lmp->atom->extra_border, nullptr);
+    ASSERT_EQ(lmp->atom->extra, nullptr);
+    ASSERT_NE(lmp->atom->mass, nullptr);
+    ASSERT_NE(lmp->atom->rmass, nullptr);
+    ASSERT_EQ(lmp->atom->ellipsoid_flag, 1);
+    ASSERT_NE(lmp->atom->ellipsoid, nullptr);
+    ASSERT_NE(lmp->atom->mass_setflag, nullptr);
+    ASSERT_NE(lmp->atom->id5p, nullptr);
+
+    auto x         = lmp->atom->x;
+    auto v         = lmp->atom->v;
+    auto type      = lmp->atom->type;
+    auto ellipsoid = lmp->atom->ellipsoid;
+    auto rmass     = lmp->atom->rmass;
+
+    auto avec  = (AtomVecEllipsoid *)hybrid->styles[1];
+    auto bonus = avec->bonus;
+
+    EXPECT_NEAR(x[GETIDX(1)][0], -0.33741452300167507, EPSILON);
+    EXPECT_NEAR(x[GETIDX(1)][1], -0.43708835412476305, EPSILON);
+    EXPECT_NEAR(x[GETIDX(1)][2], 0.6450685042019271, EPSILON);
+    EXPECT_NEAR(x[GETIDX(2)][0], -0.32142606102826937, EPSILON);
+    EXPECT_NEAR(x[GETIDX(2)][1], -0.7137743037592722, EPSILON);
+    EXPECT_NEAR(x[GETIDX(2)][2], 1.1817366147004618, EPSILON);
+    EXPECT_NEAR(x[GETIDX(3)][0], -0.130363628207774, EPSILON);
+    EXPECT_NEAR(x[GETIDX(3)][1], -0.9147144801536078, EPSILON);
+    EXPECT_NEAR(x[GETIDX(3)][2], 1.62581312195109, EPSILON);
+    EXPECT_NEAR(x[GETIDX(4)][0], 0.16795127962282844, EPSILON);
+    EXPECT_NEAR(x[GETIDX(4)][1], -0.9808507459807022, EPSILON);
+    EXPECT_NEAR(x[GETIDX(4)][2], 2.0894908590909003, EPSILON);
+    EXPECT_NEAR(x[GETIDX(5)][0], 0.46370423490634166, EPSILON);
+    EXPECT_NEAR(x[GETIDX(5)][1], -0.7803347954883079, EPSILON);
+    EXPECT_NEAR(x[GETIDX(5)][2], 2.4251986815515827, EPSILON);
+    EXPECT_NEAR(x[GETIDX(6)][0], -0.4462950185476711, EPSILON);
+    EXPECT_NEAR(x[GETIDX(6)][1], 0.09062163051035639, EPSILON);
+    EXPECT_NEAR(x[GETIDX(6)][2], 2.4668941268777607, EPSILON);
+    EXPECT_NEAR(x[GETIDX(7)][0], -0.03377054097560965, EPSILON);
+    EXPECT_NEAR(x[GETIDX(7)][1], 0.20979847489755046, EPSILON);
+    EXPECT_NEAR(x[GETIDX(7)][2], 2.078208732038921, EPSILON);
+    EXPECT_NEAR(x[GETIDX(8)][0], 0.3297325391466579, EPSILON);
+    EXPECT_NEAR(x[GETIDX(8)][1], 0.17657587120899895, EPSILON);
+    EXPECT_NEAR(x[GETIDX(8)][2], 1.7206328374934152, EPSILON);
+    EXPECT_NEAR(x[GETIDX(9)][0], 0.6063699309305985, EPSILON);
+    EXPECT_NEAR(x[GETIDX(9)][1], 0.04682595158675571, EPSILON);
+    EXPECT_NEAR(x[GETIDX(9)][2], 1.2335049647817748, EPSILON);
+    EXPECT_NEAR(x[GETIDX(10)][0], 0.8003979559814726, EPSILON);
+    EXPECT_NEAR(x[GETIDX(10)][1], -0.364393011459011, EPSILON);
+    EXPECT_NEAR(x[GETIDX(10)][2], 0.9884025318908612, EPSILON);
+
+    EXPECT_NEAR(v[GETIDX(1)][0], 0.0, EPSILON);
+    EXPECT_NEAR(v[GETIDX(1)][1], 0.0, EPSILON);
+    EXPECT_NEAR(v[GETIDX(1)][2], 0.0, EPSILON);
+    EXPECT_NEAR(v[GETIDX(2)][0], 0.0, EPSILON);
+    EXPECT_NEAR(v[GETIDX(2)][1], 0.0, EPSILON);
+    EXPECT_NEAR(v[GETIDX(2)][2], 0.0, EPSILON);
+    EXPECT_NEAR(v[GETIDX(3)][0], 0.0, EPSILON);
+    EXPECT_NEAR(v[GETIDX(3)][1], 0.0, EPSILON);
+    EXPECT_NEAR(v[GETIDX(3)][2], 0.0, EPSILON);
+    EXPECT_NEAR(v[GETIDX(4)][0], 0.0, EPSILON);
+    EXPECT_NEAR(v[GETIDX(4)][1], 0.0, EPSILON);
+    EXPECT_NEAR(v[GETIDX(4)][2], 0.0, EPSILON);
+    EXPECT_NEAR(v[GETIDX(5)][0], 0.0, EPSILON);
+    EXPECT_NEAR(v[GETIDX(5)][1], 0.0, EPSILON);
+    EXPECT_NEAR(v[GETIDX(5)][2], 0.0, EPSILON);
+    EXPECT_NEAR(v[GETIDX(6)][0], 0.0, EPSILON);
+    EXPECT_NEAR(v[GETIDX(6)][1], 0.0, EPSILON);
+    EXPECT_NEAR(v[GETIDX(6)][2], 0.0, EPSILON);
+    EXPECT_NEAR(v[GETIDX(7)][0], 0.0, EPSILON);
+    EXPECT_NEAR(v[GETIDX(7)][1], 0.0, EPSILON);
+    EXPECT_NEAR(v[GETIDX(7)][2], 0.0, EPSILON);
+    EXPECT_NEAR(v[GETIDX(8)][0], 0.0, EPSILON);
+    EXPECT_NEAR(v[GETIDX(8)][1], 0.0, EPSILON);
+    EXPECT_NEAR(v[GETIDX(8)][2], 0.0, EPSILON);
+    EXPECT_NEAR(v[GETIDX(9)][0], 0.0, EPSILON);
+    EXPECT_NEAR(v[GETIDX(9)][1], 0.0, EPSILON);
+    EXPECT_NEAR(v[GETIDX(9)][2], 0.0, EPSILON);
+    EXPECT_NEAR(v[GETIDX(10)][0], 0.0, EPSILON);
+    EXPECT_NEAR(v[GETIDX(10)][1], 0.0, EPSILON);
+    EXPECT_NEAR(v[GETIDX(10)][2], 0.0, EPSILON);
+
+    ASSERT_EQ(type[GETIDX(1)], 1);
+    ASSERT_EQ(type[GETIDX(2)], 2);
+    ASSERT_EQ(type[GETIDX(3)], 3);
+    ASSERT_EQ(type[GETIDX(4)], 4);
+    ASSERT_EQ(type[GETIDX(5)], 1);
+    ASSERT_EQ(type[GETIDX(6)], 4);
+    ASSERT_EQ(type[GETIDX(7)], 1);
+    ASSERT_EQ(type[GETIDX(8)], 2);
+    ASSERT_EQ(type[GETIDX(9)], 3);
+    ASSERT_EQ(type[GETIDX(10)], 4);
+
+    ASSERT_EQ(ellipsoid[GETIDX(1)], 0);
+    ASSERT_EQ(ellipsoid[GETIDX(2)], 1);
+    ASSERT_EQ(ellipsoid[GETIDX(3)], 2);
+    ASSERT_EQ(ellipsoid[GETIDX(4)], 3);
+    ASSERT_EQ(ellipsoid[GETIDX(5)], 4);
+    ASSERT_EQ(ellipsoid[GETIDX(6)], 5);
+    ASSERT_EQ(ellipsoid[GETIDX(7)], 6);
+    ASSERT_EQ(ellipsoid[GETIDX(8)], 7);
+    ASSERT_EQ(ellipsoid[GETIDX(9)], 8);
+    ASSERT_EQ(ellipsoid[GETIDX(10)], 9);
+
+    EXPECT_NEAR(rmass[GETIDX(1)], 3.1575, EPSILON);
+    EXPECT_NEAR(rmass[GETIDX(2)], 3.1575, EPSILON);
+    EXPECT_NEAR(rmass[GETIDX(3)], 3.1575, EPSILON);
+    EXPECT_NEAR(rmass[GETIDX(4)], 3.1575, EPSILON);
+    EXPECT_NEAR(rmass[GETIDX(5)], 3.1575, EPSILON);
+    EXPECT_NEAR(rmass[GETIDX(6)], 3.1575, EPSILON);
+    EXPECT_NEAR(rmass[GETIDX(7)], 3.1575, EPSILON);
+    EXPECT_NEAR(rmass[GETIDX(8)], 3.1575, EPSILON);
+    EXPECT_NEAR(rmass[GETIDX(9)], 3.1575, EPSILON);
+    EXPECT_NEAR(rmass[GETIDX(10)], 3.1575, EPSILON);
+
+    EXPECT_NEAR(bonus[0].shape[0], 0.5869922515711705, EPSILON);
+    EXPECT_NEAR(bonus[0].shape[1], 0.5869922515711705, EPSILON);
+    EXPECT_NEAR(bonus[0].shape[2], 0.5869922515711705, EPSILON);
+    EXPECT_NEAR(bonus[1].shape[0], 0.5869922515711705, EPSILON);
+    EXPECT_NEAR(bonus[1].shape[1], 0.5869922515711705, EPSILON);
+    EXPECT_NEAR(bonus[1].shape[2], 0.5869922515711705, EPSILON);
+    EXPECT_NEAR(bonus[2].shape[0], 0.5869922515711705, EPSILON);
+    EXPECT_NEAR(bonus[2].shape[1], 0.5869922515711705, EPSILON);
+    EXPECT_NEAR(bonus[2].shape[2], 0.5869922515711705, EPSILON);
+    EXPECT_NEAR(bonus[3].shape[0], 0.5869922515711705, EPSILON);
+    EXPECT_NEAR(bonus[3].shape[1], 0.5869922515711705, EPSILON);
+    EXPECT_NEAR(bonus[3].shape[2], 0.5869922515711705, EPSILON);
+    EXPECT_NEAR(bonus[4].shape[0], 0.5869922515711705, EPSILON);
+    EXPECT_NEAR(bonus[4].shape[1], 0.5869922515711705, EPSILON);
+    EXPECT_NEAR(bonus[4].shape[2], 0.5869922515711705, EPSILON);
+    EXPECT_NEAR(bonus[5].shape[0], 0.5869922515711705, EPSILON);
+    EXPECT_NEAR(bonus[5].shape[1], 0.5869922515711705, EPSILON);
+    EXPECT_NEAR(bonus[5].shape[2], 0.5869922515711705, EPSILON);
+    EXPECT_NEAR(bonus[6].shape[0], 0.5869922515711705, EPSILON);
+    EXPECT_NEAR(bonus[6].shape[1], 0.5869922515711705, EPSILON);
+    EXPECT_NEAR(bonus[6].shape[2], 0.5869922515711705, EPSILON);
+    EXPECT_NEAR(bonus[7].shape[0], 0.5869922515711705, EPSILON);
+    EXPECT_NEAR(bonus[7].shape[1], 0.5869922515711705, EPSILON);
+    EXPECT_NEAR(bonus[7].shape[2], 0.5869922515711705, EPSILON);
+    EXPECT_NEAR(bonus[8].shape[0], 0.5869922515711705, EPSILON);
+    EXPECT_NEAR(bonus[8].shape[1], 0.5869922515711705, EPSILON);
+    EXPECT_NEAR(bonus[8].shape[2], 0.5869922515711705, EPSILON);
+    EXPECT_NEAR(bonus[9].shape[0], 0.5869922515711705, EPSILON);
+    EXPECT_NEAR(bonus[9].shape[1], 0.5869922515711705, EPSILON);
+    EXPECT_NEAR(bonus[9].shape[2], 0.5869922515711705, EPSILON);
+
+    EXPECT_NEAR(bonus[0].quat[0], 0.9890278201757743, EPSILON); 
+    EXPECT_NEAR(bonus[0].quat[1], 0.01779228232037064, EPSILON);
+    EXPECT_NEAR(bonus[0].quat[2], -0.14337734159225404, EPSILON);
+    EXPECT_NEAR(bonus[0].quat[3], 0.030827642240801516, EPSILON);
+    EXPECT_NEAR(bonus[1].quat[0], 0.939687458852748, EPSILON);
+    EXPECT_NEAR(bonus[1].quat[1], 0.04174166924055095, EPSILON);
+    EXPECT_NEAR(bonus[1].quat[2], -0.023337773785056866, EPSILON);
+    EXPECT_NEAR(bonus[1].quat[3], 0.338674565089608, EPSILON);
+    EXPECT_NEAR(bonus[2].quat[0], 0.8210113150655425, EPSILON);
+    EXPECT_NEAR(bonus[2].quat[1], 0.03012140921736572, EPSILON);
+    EXPECT_NEAR(bonus[2].quat[2], 0.017666019956944813, EPSILON);
+    EXPECT_NEAR(bonus[2].quat[3], 0.5698429897612057, EPSILON);
+    EXPECT_NEAR(bonus[3].quat[0], 0.6623662858285051, EPSILON);
+    EXPECT_NEAR(bonus[3].quat[1], -0.028186343967346823, EPSILON);
+    EXPECT_NEAR(bonus[3].quat[2], 0.022942552517501488, EPSILON);
+    EXPECT_NEAR(bonus[3].quat[3], 0.7482981175276918, EPSILON);
+    EXPECT_NEAR(bonus[4].quat[0], 0.3601488726765216, EPSILON);
+    EXPECT_NEAR(bonus[4].quat[1], 0.0513614985821682, EPSILON);
+    EXPECT_NEAR(bonus[4].quat[2], 0.0724224158335286, EPSILON);
+    EXPECT_NEAR(bonus[4].quat[3], 0.9286602067807472, EPSILON);
+    EXPECT_NEAR(bonus[5].quat[0], 0.11941234710084649, EPSILON);
+    EXPECT_NEAR(bonus[5].quat[1], 0.9244660117493703, EPSILON);
+    EXPECT_NEAR(bonus[5].quat[2], -0.35317942248051865, EPSILON);
+    EXPECT_NEAR(bonus[5].quat[3], -0.07979711784524246, EPSILON);
+    EXPECT_NEAR(bonus[6].quat[0], -0.17949125421205164, EPSILON);
+    EXPECT_NEAR(bonus[6].quat[1], 0.7412884899431119, EPSILON);
+    EXPECT_NEAR(bonus[6].quat[2], -0.6379094464220707, EPSILON);
+    EXPECT_NEAR(bonus[6].quat[3], 0.1065166771202199, EPSILON);
+    EXPECT_NEAR(bonus[7].quat[0], -0.10483691088405202, EPSILON);
+    EXPECT_NEAR(bonus[7].quat[1], 0.5508895999584645, EPSILON);
+    EXPECT_NEAR(bonus[7].quat[2], -0.8250090480220789, EPSILON);
+    EXPECT_NEAR(bonus[7].quat[3], 0.06992811634525403, EPSILON);
+    EXPECT_NEAR(bonus[8].quat[0], 0.07777239911646, EPSILON);
+    EXPECT_NEAR(bonus[8].quat[1], -0.3724087549185288, EPSILON);
+    EXPECT_NEAR(bonus[8].quat[2], 0.9103052384821374, EPSILON);
+    EXPECT_NEAR(bonus[8].quat[3], -0.1631181963720798, EPSILON);
+    EXPECT_NEAR(bonus[9].quat[0], 0.16279109707978262, EPSILON);
+    EXPECT_NEAR(bonus[9].quat[1], 0.027148630125149613, EPSILON);
+    EXPECT_NEAR(bonus[9].quat[2], 0.9849325709665359, EPSILON);
+    EXPECT_NEAR(bonus[9].quat[3], -0.0516705065113425, EPSILON);
+
+    auto num_bond  = lmp->atom->num_bond;
+    auto bond_type = lmp->atom->bond_type;
+    auto bond_atom = lmp->atom->bond_atom;
+    auto id5p      = lmp->atom->id5p;
+
+    ASSERT_EQ(num_bond[GETIDX(1)], 1);
+    ASSERT_EQ(num_bond[GETIDX(2)], 1);
+    ASSERT_EQ(num_bond[GETIDX(3)], 1);
+    ASSERT_EQ(num_bond[GETIDX(4)], 1);
+    ASSERT_EQ(num_bond[GETIDX(5)], 0);
+    ASSERT_EQ(num_bond[GETIDX(6)], 1);
+    ASSERT_EQ(num_bond[GETIDX(7)], 1);
+    ASSERT_EQ(num_bond[GETIDX(8)], 1);
+    ASSERT_EQ(num_bond[GETIDX(9)], 1);
+    ASSERT_EQ(num_bond[GETIDX(10)], 0);
+
+    ASSERT_EQ(bond_type[GETIDX(1)][0], 1);
+    ASSERT_EQ(bond_type[GETIDX(2)][0], 1);
+    ASSERT_EQ(bond_type[GETIDX(3)][0], 1);
+    ASSERT_EQ(bond_type[GETIDX(4)][0], 1);
+    ASSERT_EQ(bond_type[GETIDX(6)][0], 1);
+    ASSERT_EQ(bond_type[GETIDX(7)][0], 1);
+    ASSERT_EQ(bond_type[GETIDX(8)][0], 1);
+    ASSERT_EQ(bond_type[GETIDX(9)][0], 1);
+
+    ASSERT_EQ(bond_atom[GETIDX(1)][0], 2);
+    ASSERT_EQ(bond_atom[GETIDX(2)][0], 3);
+    ASSERT_EQ(bond_atom[GETIDX(3)][0], 4);
+    ASSERT_EQ(bond_atom[GETIDX(4)][0], 5);
+    ASSERT_EQ(bond_atom[GETIDX(6)][0], 7);
+    ASSERT_EQ(bond_atom[GETIDX(7)][0], 8);
+    ASSERT_EQ(bond_atom[GETIDX(8)][0], 9);
+    ASSERT_EQ(bond_atom[GETIDX(9)][0], 10);
+
+    ASSERT_EQ(id5p[GETIDX(1)], 2);
+    ASSERT_EQ(id5p[GETIDX(2)], 3);
+    ASSERT_EQ(id5p[GETIDX(3)], 4);
+    ASSERT_EQ(id5p[GETIDX(4)], 5);
+    ASSERT_EQ(id5p[GETIDX(5)], -1);
+    ASSERT_EQ(id5p[GETIDX(6)], 7);
+    ASSERT_EQ(id5p[GETIDX(7)], 8);
+    ASSERT_EQ(id5p[GETIDX(8)], 9);
+    ASSERT_EQ(id5p[GETIDX(9)], 10);
+    ASSERT_EQ(id5p[GETIDX(10)], -1);
+
+    END_HIDE_OUTPUT();
+
+}
+
 } // namespace LAMMPS_NS
 
 int main(int argc, char **argv)