Corrected virial, altered calculation of stacking strength

doc/src/pair_oxdna.txt

@@ -23,9 +23,11 @@ style1 = {hybrid/overlay oxdna/excv oxdna/stk oxdna/hbond oxdna/xstk oxdna/coaxs
 style2 = {oxdna/excv} or {oxdna/stk} or {oxdna/hbond} or {oxdna/xstk} or {oxdna/coaxstk}
 args = list of arguments for these particular styles :ul
 
-  {oxdna/stk} args = seq T 6.0 0.4 0.9 0.32 0.6 1.3 0 0.8 0.9 0 0.95 0.9 0 0.95 2.0 0.65 2.0 0.65
+  {oxdna/stk} args = seq T xi kappa 6.0 0.4 0.9 0.32 0.6 1.3 0 0.8 0.9 0 0.95 0.9 0 0.95 2.0 0.65 2.0 0.65
     seq = seqav (for average sequence stacking strength) or seqdep (for sequence-dependent stacking strength)
     T = temperature (oxDNA units, 0.1 = 300 K)
+    xi = temperature-independent coefficient in stacking strength
+    kappa = coefficient of linear temperature dependence in stacking strength
   {oxdna/hbond} args = seq eps 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
     seq = seqav (for average sequence base-pairing strength) or seqdep (for sequence-dependent base-pairing strength)
     eps = 1.077 (between base pairs A-T and C-G) or 0 (all other pairs) :pre
@@ -34,7 +36,7 @@ args = list of arguments for these particular styles :ul
 
 pair_style hybrid/overlay oxdna/excv oxdna/stk oxdna/hbond oxdna/xstk oxdna/coaxstk
 pair_coeff * * oxdna/excv    2.0 0.7 0.675 2.0 0.515 0.5 2.0 0.33 0.32
-pair_coeff * * oxdna/stk     seqdep 0.1 6.0 0.4 0.9 0.32 0.6 1.3 0 0.8 0.9 0 0.95 0.9 0 0.95 2.0 0.65 2.0 0.65
+pair_coeff * * oxdna/stk     seqdep 0.1 1.3448 2.6568 6.0 0.4 0.9 0.32 0.6 1.3 0 0.8 0.9 0 0.95 0.9 0 0.95 2.0 0.65 2.0 0.65
 pair_coeff * * oxdna/hbond   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
 pair_coeff 1 4 oxdna/hbond   seqdep 1.077 8.0 0.4 0.75 0.34 0.7 1.5 0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45
 pair_coeff 2 3 oxdna/hbond   seqdep 1.077 8.0 0.4 0.75 0.34 0.7 1.5 0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45
@@ -62,7 +64,7 @@ NOTE: These pair styles have to be used together with the related oxDNA bond sty
 {oxdna/fene} for the connectivity of the phosphate backbone (see also documentation of
 "bond_style oxdna/fene"_bond_oxdna.html). Most of the coefficients
 in the above example have to be kept fixed and cannot be changed without reparameterizing the entire model.
-Exceptions are the first and second coefficient after {oxdna/stk} (seq=seqdep and T=0.1 in the above example)
+Exceptions are the first four coefficients after {oxdna/stk} (seq=seqdep, T=0.1, xi=1.3448 and kappa=2.6568 in the above example)
 and the first coefficient after {oxdna/hbond} (seq=seqdep in the above example).
 When using a Langevin thermostat, e.g. through "fix langevin"_fix_langevin.html
 or "fix nve/dotc/langevin"_fix_nve_dotc_langevin.html

doc/src/pair_oxdna2.txt

@@ -24,10 +24,12 @@ style1 = {hybrid/overlay oxdna2/excv oxdna2/stk oxdna2/hbond oxdna2/xstk oxdna2/
 style2 = {oxdna2/excv} or {oxdna2/stk} or {oxdna2/hbond} or {oxdna2/xstk} or {oxdna2/coaxstk} or {oxdna2/dh}
 args = list of arguments for these particular styles :ul
 
-  {oxdna2/stk} args = seq T 6.0 0.4 0.9 0.32 0.6 1.3 0 0.8 0.9 0 0.95 0.9 0 0.95 2.0 0.65 2.0 0.65
+  {oxdna2/stk} args = seq T xi kappa 6.0 0.4 0.9 0.32 0.6 1.3 0 0.8 0.9 0 0.95 0.9 0 0.95 2.0 0.65 2.0 0.65
     seq = seqav (for average sequence stacking strength) or seqdep (for sequence-dependent stacking strength)
     T = temperature (oxDNA units, 0.1 = 300 K)
-  {oxdna/hbond} args = seq eps 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
+    xi = temperature-independent coefficient in stacking strength 
+    kappa = coefficient of linear temperature dependence in stacking strength 
+  {oxdna2/hbond} args = seq eps 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
     seq = seqav (for average sequence base-pairing strength) or seqdep (for sequence-dependent base-pairing strength)
     eps = 1.0678 (between base pairs A-T and C-G) or 0 (all other pairs)
   {oxdna2/dh} args = T rhos qeff
@@ -39,7 +41,7 @@ args = list of arguments for these particular styles :ul
 
 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 0.1 6.0 0.4 0.9 0.32 0.6 1.3 0 0.8 0.9 0 0.95 0.9 0 0.95 2.0 0.65 2.0 0.65
+pair_coeff * * oxdna2/stk     seqdep 0.1 1.3523 2.6717 6.0 0.4 0.9 0.32 0.6 1.3 0 0.8 0.9 0 0.95 0.9 0 0.95 2.0 0.65 2.0 0.65
 pair_coeff * * 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
 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
 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
@@ -68,8 +70,8 @@ NOTE: These pair styles have to be used together with the related oxDNA2 bond st
 {oxdna2/fene} for the connectivity of the phosphate backbone (see also documentation of
 "bond_style oxdna2/fene"_bond_oxdna.html). Most of the coefficients
 in the above example have to be kept fixed and cannot be changed without reparameterizing the entire model.
-Exceptions are the first and the second coefficient after {oxdna2/stk} (seq=seqdep and T=0.1 in the above example),
+Exceptions are the first four coefficients after {oxdna2/stk} (seq=seqdep, T=0.1, xi=1.3523 and kappa=2.6717 in the above example),
-the first coefficient after {oxdna/hbond} (seq=seqdep in the above example) and the three coefficients
+the first coefficient after {oxdna2/hbond} (seq=seqdep in the above example) and the three coefficients
 after {oxdna2/dh} (T=0.1, rhos=1.0, qeff=0.815 in the above example). When using a Langevin thermostat
 e.g. through "fix langevin"_fix_langevin.html or "fix nve/dotc/langevin"_fix_nve_dotc_langevin.html
 the temperature coefficients have to be matched to the one used in the fix.

src/USER-CGDNA/bond_oxdna2_fene.cpp

@@ -37,8 +37,8 @@ BondOxdna2Fene::~BondOxdna2Fene()
 /* ----------------------------------------------------------------------
     compute vector COM-sugar-phosphate backbone interaction site in oxDNA2
 ------------------------------------------------------------------------- */
-void BondOxdna2Fene::compute_interaction_sites(double e1[3],
+void BondOxdna2Fene::compute_interaction_sites(double e1[3], double e2[3],
-  double e2[3], double r[3])
+  double /*e3*/[3], double r[3])
 {
   double d_cs_x=-0.34, d_cs_y=+0.3408;
 

src/USER-CGDNA/bond_oxdna2_fene.h

@@ -28,7 +28,8 @@ class BondOxdna2Fene : public BondOxdnaFene {
  public:
   BondOxdna2Fene(class LAMMPS *);
   virtual ~BondOxdna2Fene();
-  virtual void compute_interaction_sites(double *, double *, double *);
+  virtual void compute_interaction_sites(double *, double *, double *,
+    double *);
 };
 
 }

src/USER-CGDNA/bond_oxdna_fene.cpp

@@ -55,8 +55,8 @@ BondOxdnaFene::~BondOxdnaFene()
 /* ----------------------------------------------------------------------
     compute vector COM-sugar-phosphate backbone interaction site in oxDNA
 ------------------------------------------------------------------------- */
-void BondOxdnaFene::compute_interaction_sites(double e1[3],
+void BondOxdnaFene::compute_interaction_sites(double e1[3], double /*e2*/[3],
-  double /*e2*/[3], double r[3])
+  double /*e3*/[3], double r[3])
 {
   double d_cs=-0.4;
 
@@ -66,6 +66,90 @@ void BondOxdnaFene::compute_interaction_sites(double e1[3],
 
 }
 
+/* ----------------------------------------------------------------------
+   tally energy and virial into global and per-atom accumulators
+------------------------------------------------------------------------- */
+
+void BondOxdnaFene::ev_tally_xyz(int i, int j, int nlocal, int newton_bond,
+                    double ebond,
+                    double fx, double fy, double fz,
+                    double delx, double dely, double delz)
+{
+  double ebondhalf,v[6];
+
+  if (eflag_either) {
+    if (eflag_global) {
+      if (newton_bond) energy += ebond;
+      else {
+        ebondhalf = 0.5*ebond;
+        if (i < nlocal) energy += ebondhalf;
+        if (j < nlocal) energy += ebondhalf;
+      }
+    }
+    if (eflag_atom) {
+      ebondhalf = 0.5*ebond;
+      if (newton_bond || i < nlocal) eatom[i] += ebondhalf;
+      if (newton_bond || j < nlocal) eatom[j] += ebondhalf;
+    }
+  }
+
+  if (vflag_either) {
+    v[0] = delx*fx;
+    v[1] = dely*fy;
+    v[2] = delz*fz;
+    v[3] = delx*fy;
+    v[4] = delx*fz;
+    v[5] = dely*fz;
+
+    if (vflag_global) {
+      if (newton_bond) {
+        virial[0] += v[0];
+        virial[1] += v[1];
+        virial[2] += v[2];
+        virial[3] += v[3];
+        virial[4] += v[4];
+        virial[5] += v[5];
+      } else {
+        if (i < nlocal) {
+          virial[0] += 0.5*v[0];
+          virial[1] += 0.5*v[1];
+          virial[2] += 0.5*v[2];
+          virial[3] += 0.5*v[3];
+          virial[4] += 0.5*v[4];
+          virial[5] += 0.5*v[5];
+        }
+        if (j < nlocal) {
+          virial[0] += 0.5*v[0];
+          virial[1] += 0.5*v[1];
+          virial[2] += 0.5*v[2];
+          virial[3] += 0.5*v[3];
+          virial[4] += 0.5*v[4];
+          virial[5] += 0.5*v[5];
+        }
+      }
+    }
+
+    if (vflag_atom) {
+      if (newton_bond || i < nlocal) {
+        vatom[i][0] += 0.5*v[0];
+        vatom[i][1] += 0.5*v[1];
+        vatom[i][2] += 0.5*v[2];
+        vatom[i][3] += 0.5*v[3];
+        vatom[i][4] += 0.5*v[4];
+        vatom[i][5] += 0.5*v[5];
+      }
+      if (newton_bond || j < nlocal) {
+        vatom[j][0] += 0.5*v[0];
+        vatom[j][1] += 0.5*v[1];
+        vatom[j][2] += 0.5*v[2];
+        vatom[j][3] += 0.5*v[3];
+        vatom[j][4] += 0.5*v[4];
+        vatom[j][5] += 0.5*v[5];
+      }
+    }
+  }
+}
+
 /* ----------------------------------------------------------------------
    compute function for oxDNA FENE-bond interaction
    s=sugar-phosphate backbone site, b=base site, st=stacking site
@@ -112,8 +196,8 @@ void BondOxdnaFene::compute(int eflag, int vflag)
     MathExtra::q_to_exyz(qb,bx,by,bz);
 
     // vector COM-backbone site a and b
-    compute_interaction_sites(ax,ay,ra_cs);
+    compute_interaction_sites(ax,ay,az,ra_cs);
-    compute_interaction_sites(bx,by,rb_cs);
+    compute_interaction_sites(bx,by,bz,rb_cs);
 
     // vector backbone site b to a
     delr[0] = x[a][0] + ra_cs[0] - x[b][0] - rb_cs[0];
@@ -182,7 +266,11 @@ void BondOxdnaFene::compute(int eflag, int vflag)
     }
 
     // increment energy and virial
-    if (evflag) ev_tally(a,b,nlocal,newton_bond,ebond,fbond,delr[0],delr[1],delr[2]);
+    // NOTE: The virial is calculated on the 'molecular' basis.
+    // (see G. Ciccotti and J.P. Ryckaert, Comp. Phys. Rep. 4, 345-392 (1986))
+
+    if (evflag) ev_tally_xyz(a,b,nlocal,newton_bond,ebond,
+        delf[0],delf[1],delf[2],x[a][0]-x[b][0],x[a][1]-x[b][1],x[a][2]-x[b][2]);
 
   }
 

src/USER-CGDNA/bond_oxdna_fene.h

@@ -28,7 +28,8 @@ class BondOxdnaFene : public Bond {
  public:
   BondOxdnaFene(class LAMMPS *);
   virtual ~BondOxdnaFene();
-  virtual void compute_interaction_sites(double *, double *, double *);
+  virtual void compute_interaction_sites(double *, double *, double *,
+    double *);
   virtual void compute(int, int);
   void coeff(int, char **);
   void init_style();
@@ -42,6 +43,7 @@ class BondOxdnaFene : public Bond {
   double *k,*Delta,*r0; // FENE
 
   void allocate();
+  void ev_tally_xyz(int, int, int, int, double, double, double, double, double, double, double);
 };
 
 }

src/USER-CGDNA/pair_oxdna2_coaxstk.cpp

@@ -383,7 +383,11 @@ void PairOxdna2Coaxstk::compute(int eflag, int vflag)
       }
 
       // increment energy and virial
-      if (evflag) ev_tally(a,b,nlocal,newton_pair,evdwl,0.0,fpair,delr_st[0],delr_st[1],delr_st[2]);
+      // NOTE: The virial is calculated on the 'molecular' basis.
+      // (see G. Ciccotti and J.P. Ryckaert, Comp. Phys. Rep. 4, 345-392 (1986))
+
+      if (evflag) ev_tally_xyz(a,b,nlocal,newton_pair,evdwl,0.0,
+          delf[0],delf[1],delf[2],x[a][0]-x[b][0],x[a][1]-x[b][1],x[a][2]-x[b][2]);
 
       // pure torques not expressible as r x f
 

src/USER-CGDNA/pair_oxdna2_dh.cpp

@@ -223,10 +223,14 @@ void PairOxdna2Dh::compute(int eflag, int vflag)
 
         }
 
-        // increment energy and virial
 
-        if (evflag) ev_tally(a,b,nlocal,newton_pair,
+        // increment energy and virial
-                evdwl,0.0,fpair,delr[0],delr[1],delr[2]);
+        // NOTE: The virial is calculated on the 'molecular' basis.
+        // (see G. Ciccotti and J.P. Ryckaert, Comp. Phys. Rep. 4, 345-392 (1986))
+
+        if (evflag) ev_tally_xyz(a,b,nlocal,newton_pair,evdwl,0.0,
+            delf[0],delf[1],delf[2],x[a][0]-x[b][0],x[a][1]-x[b][1],x[a][2]-x[b][2]);
+
       }
 
     }

src/USER-CGDNA/pair_oxdna2_excv.cpp

@@ -39,8 +39,8 @@ PairOxdna2Excv::~PairOxdna2Excv()
 /* ----------------------------------------------------------------------
     compute vector COM-excluded volume interaction sites in oxDNA2
 ------------------------------------------------------------------------- */
-void PairOxdna2Excv::compute_interaction_sites(double e1[3],
+void PairOxdna2Excv::compute_interaction_sites(double e1[3], double e2[3], 
-  double e2[3], double rs[3], double rb[3])
+    double /*e3*/[3], double rs[3], double rb[3])
 {
   double d_cs_x=-0.34, d_cs_y=+0.3408, d_cb=+0.4;
 

src/USER-CGDNA/pair_oxdna2_excv.h

@@ -28,7 +28,7 @@ class PairOxdna2Excv : public PairOxdnaExcv {
  public:
   PairOxdna2Excv(class LAMMPS *);
   virtual ~PairOxdna2Excv();
-  virtual void compute_interaction_sites(double *,
+  virtual void compute_interaction_sites(double *, double *,
     double *, double *, double *);
 };
 

src/USER-CGDNA/pair_oxdna2_stk.cpp

@@ -1,50 +0,0 @@
-/* ----------------------------------------------------------------------
-   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)
-------------------------------------------------------------------------- */
-
-#include <cmath>
-#include <cstdio>
-#include <cstdlib>
-#include <cstring>
-#include "pair_oxdna2_stk.h"
-
-using namespace LAMMPS_NS;
-
-/* ---------------------------------------------------------------------- */
-
-PairOxdna2Stk::PairOxdna2Stk(LAMMPS *lmp) : PairOxdnaStk(lmp)
-{
-
-}
-
-/* ---------------------------------------------------------------------- */
-
-PairOxdna2Stk::~PairOxdna2Stk()
-{
-
-}
-
-/* ----------------------------------------------------------------------
-   return temperature dependent oxDNA2 stacking strength
-------------------------------------------------------------------------- */
-
-double PairOxdna2Stk::stacking_strength(double T)
-{
-  double eps;
-
-  eps = 1.3523 + 2.6717 * T;
-
-  return eps;
-}

src/USER-CGDNA/pair_oxdna2_stk.h

@@ -1,53 +0,0 @@
-/* -*- c++ -*- ----------------------------------------------------------
-   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
-   http://lammps.sandia.gov, Sandia National Laboratories
-   Steve Plimpton, sjplimp@sandia.gov
-
-   Copyright (2003) Sandia Corporation.  Under the terms of Contract
-   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
-   certain rights in this software.  This software is distributed under
-   the GNU General Public License.
-
-   See the README file in the top-level LAMMPS directory.
-------------------------------------------------------------------------- */
-
-#ifdef PAIR_CLASS
-
-PairStyle(oxdna2/stk,PairOxdna2Stk)
-
-#else
-
-#ifndef LMP_PAIR_OXDNA2_STK_H
-#define LMP_PAIR_OXDNA2_STK_H
-
-#include "pair_oxdna_stk.h"
-
-namespace LAMMPS_NS {
-
-class PairOxdna2Stk : public PairOxdnaStk {
- public:
-  PairOxdna2Stk(class LAMMPS *);
-  virtual ~PairOxdna2Stk();
-
- protected:
-  virtual double stacking_strength(double);
-};
-
-}
-
-#endif
-#endif
-
-/* ERROR/WARNING messages:
-
-E: Illegal ... command
-
-Self-explanatory.  Check the input script syntax and compare to the
-documentation for the command.  You can use -echo screen as a
-command-line option when running LAMMPS to see the offending line.
-
-E: Incorrect args for pair coefficients
-
-Self-explanatory.  Check the input script or data file.
-
-*/

src/USER-CGDNA/pair_oxdna_coaxstk.cpp

@@ -450,7 +450,11 @@ void PairOxdnaCoaxstk::compute(int eflag, int vflag)
       }
 
       // increment energy and virial
-      if (evflag) ev_tally(a,b,nlocal,newton_pair,evdwl,0.0,fpair,delr_st[0],delr_st[1],delr_st[2]);
+      // NOTE: The virial is calculated on the 'molecular' basis.
+      // (see G. Ciccotti and J.P. Ryckaert, Comp. Phys. Rep. 4, 345-392 (1986))
+
+      if (evflag) ev_tally_xyz(a,b,nlocal,newton_pair,evdwl,0.0,
+          delf[0],delf[1],delf[2],x[a][0]-x[b][0],x[a][1]-x[b][1],x[a][2]-x[b][2]);
 
       // pure torques not expressible as r x f
 

src/USER-CGDNA/pair_oxdna_excv.cpp

@@ -91,8 +91,8 @@ PairOxdnaExcv::~PairOxdnaExcv()
 /* ----------------------------------------------------------------------
     compute vector COM-excluded volume interaction sites in oxDNA
 ------------------------------------------------------------------------- */
-void PairOxdnaExcv::compute_interaction_sites(double e1[3],
+void PairOxdnaExcv::compute_interaction_sites(double e1[3], double /*e2*/[3], 
-  double /*e2*/[3], double rs[3], double rb[3])
+    double /*e3*/[3], double rs[3], double rb[3])
 {
   double d_cs=-0.4, d_cb=+0.4;
 
@@ -162,7 +162,7 @@ void PairOxdnaExcv::compute(int eflag, int vflag)
     MathExtra::q_to_exyz(qa,ax,ay,az);
 
     // vector COM - backbone and base site a
-    compute_interaction_sites(ax,ay,ra_cs,ra_cb);
+    compute_interaction_sites(ax,ay,az,ra_cs,ra_cb);
 
     rtmp_s[0] = x[a][0] + ra_cs[0];
     rtmp_s[1] = x[a][1] + ra_cs[1];
@@ -187,7 +187,7 @@ void PairOxdnaExcv::compute(int eflag, int vflag)
       MathExtra::q_to_exyz(qb,bx,by,bz);
 
       // vector COM - backbone and base site b
-      compute_interaction_sites(bx,by,rb_cs,rb_cb);
+      compute_interaction_sites(bx,by,bz,rb_cs,rb_cb);
 
       // vector backbone site b to a
       delr_ss[0] = rtmp_s[0] - (x[b][0] + rb_cs[0]);
@@ -225,14 +225,17 @@ void PairOxdnaExcv::compute(int eflag, int vflag)
         fpair *= factor_lj;
         evdwl *= factor_lj;
 
-        // increment energy and virial
-        if (evflag) ev_tally(a,b,nlocal,newton_pair,
-                evdwl,0.0,fpair,delr_ss[0],delr_ss[1],delr_ss[2]);
-
         delf[0] = delr_ss[0]*fpair;
         delf[1] = delr_ss[1]*fpair;
         delf[2] = delr_ss[2]*fpair;
 
+        // increment energy and virial
+        // NOTE: The virial is calculated on the 'molecular' basis.
+        // (see G. Ciccotti and J.P. Ryckaert, Comp. Phys. Rep. 4, 345-392 (1986))
+
+        if (evflag) ev_tally_xyz(a,b,nlocal,newton_pair,evdwl,0.0,
+            delf[0],delf[1],delf[2],x[a][0]-x[b][0],x[a][1]-x[b][1],x[a][2]-x[b][2]);
+
         f[a][0] += delf[0];
         f[a][1] += delf[1];
         f[a][2] += delf[2];
@@ -259,21 +262,20 @@ void PairOxdnaExcv::compute(int eflag, int vflag)
 
       }
 
-
       // backbone-base
       if (rsq_sb < cutsq_sb_c[atype][btype]) {
 
         evdwl = F3(rsq_sb,cutsq_sb_ast[atype][btype],cut_sb_c[atype][btype],lj1_sb[atype][btype],
                         lj2_sb[atype][btype],epsilon_sb[atype][btype],b_sb[atype][btype],fpair);
 
-        // increment energy and virial
-        if (evflag) ev_tally(a,b,nlocal,newton_pair,
-                evdwl,0.0,fpair,delr_sb[0],delr_sb[1],delr_sb[2]);
-
         delf[0] = delr_sb[0]*fpair;
         delf[1] = delr_sb[1]*fpair;
         delf[2] = delr_sb[2]*fpair;
 
+        // increment energy and virial
+        if (evflag) ev_tally_xyz(a,b,nlocal,newton_pair,evdwl,0.0,
+            delf[0],delf[1],delf[2],x[a][0]-x[b][0],x[a][1]-x[b][1],x[a][2]-x[b][2]);
+
         f[a][0] += delf[0];
         f[a][1] += delf[1];
         f[a][2] += delf[2];
@@ -306,14 +308,14 @@ void PairOxdnaExcv::compute(int eflag, int vflag)
         evdwl = F3(rsq_bs,cutsq_sb_ast[atype][btype],cut_sb_c[atype][btype],lj1_sb[atype][btype],
                         lj2_sb[atype][btype],epsilon_sb[atype][btype],b_sb[atype][btype],fpair);
 
-        // increment energy and virial
-        if (evflag) ev_tally(a,b,nlocal,newton_pair,
-                evdwl,0.0,fpair,delr_bs[0],delr_bs[1],delr_bs[2]);
-
         delf[0] = delr_bs[0]*fpair;
         delf[1] = delr_bs[1]*fpair;
         delf[2] = delr_bs[2]*fpair;
 
+        // increment energy and virial
+        if (evflag) ev_tally_xyz(a,b,nlocal,newton_pair,evdwl,0.0,
+            delf[0],delf[1],delf[2],x[a][0]-x[b][0],x[a][1]-x[b][1],x[a][2]-x[b][2]);
+
         f[a][0] += delf[0];
         f[a][1] += delf[1];
         f[a][2] += delf[2];
@@ -346,14 +348,14 @@ void PairOxdnaExcv::compute(int eflag, int vflag)
         evdwl = F3(rsq_bb,cutsq_bb_ast[atype][btype],cut_bb_c[atype][btype],lj1_bb[atype][btype],
                         lj2_bb[atype][btype],epsilon_bb[atype][btype],b_bb[atype][btype],fpair);
 
-        // increment energy and virial
-        if (evflag) ev_tally(a,b,nlocal,newton_pair,
-                evdwl,0.0,fpair,delr_bb[0],delr_bb[1],delr_bb[2]);
-
         delf[0] = delr_bb[0]*fpair;
         delf[1] = delr_bb[1]*fpair;
         delf[2] = delr_bb[2]*fpair;
 
+        // increment energy and virial
+        if (evflag) ev_tally_xyz(a,b,nlocal,newton_pair,evdwl,0.0,
+            delf[0],delf[1],delf[2],x[a][0]-x[b][0],x[a][1]-x[b][1],x[a][2]-x[b][2]);
+
         f[a][0] += delf[0];
         f[a][1] += delf[1];
         f[a][2] += delf[2];

src/USER-CGDNA/pair_oxdna_excv.h

@@ -28,7 +28,7 @@ class PairOxdnaExcv : public Pair {
  public:
   PairOxdnaExcv(class LAMMPS *);
   virtual ~PairOxdnaExcv();
-  virtual void compute_interaction_sites(double *, double *,
+  virtual void compute_interaction_sites(double *, double *, double *,
     double *, double *);
   virtual void compute(int, int);
   void settings(int, char **);

src/USER-CGDNA/pair_oxdna_hbond.cpp

@@ -410,7 +410,11 @@ void PairOxdnaHbond::compute(int eflag, int vflag)
       }
 
       // increment energy and virial
-      if (evflag) ev_tally(a,b,nlocal,newton_pair,evdwl,0.0,fpair,delr_hb[0],delr_hb[1],delr_hb[2]);
+      // NOTE: The virial is calculated on the 'molecular' basis.
+      // (see G. Ciccotti and J.P. Ryckaert, Comp. Phys. Rep. 4, 345-392 (1986))
+
+      if (evflag) ev_tally_xyz(a,b,nlocal,newton_pair,evdwl,0.0,
+          delf[0],delf[1],delf[2],x[a][0]-x[b][0],x[a][1]-x[b][1],x[a][2]-x[b][2]);
 
       // pure torques not expressible as r x f
 

src/USER-CGDNA/pair_oxdna_stk.cpp

@@ -106,6 +106,94 @@ PairOxdnaStk::~PairOxdnaStk()
   }
 }
 
+/* ----------------------------------------------------------------------
+   tally energy and virial into global and per-atom accumulators
+
+   NOTE: Although this is a pair style interaction, the algorithm below
+   follows the virial incrementation of the bond style. This is because 
+   the bond topology is used in the main compute loop.
+------------------------------------------------------------------------- */
+
+void PairOxdnaStk::ev_tally_xyz(int i, int j, int nlocal, int newton_bond,
+                    double evdwl,
+                    double fx, double fy, double fz,
+                    double delx, double dely, double delz)
+{
+  double evdwlhalf,v[6];
+
+  if (eflag_either) {
+    if (eflag_global) {
+      if (newton_bond) eng_vdwl += evdwl;
+      else {
+        evdwlhalf = 0.5*evdwl;
+        if (i < nlocal) eng_vdwl += evdwlhalf;
+        if (j < nlocal) eng_vdwl += evdwlhalf;
+      }
+    }
+    if (eflag_atom) {
+      evdwlhalf = 0.5*evdwl;
+      if (newton_bond || i < nlocal) eatom[i] += evdwlhalf;
+      if (newton_bond || j < nlocal) eatom[j] += evdwlhalf;
+    }
+  }
+
+  if (vflag_either) {
+    v[0] = delx*fx;
+    v[1] = dely*fy;
+    v[2] = delz*fz;
+    v[3] = delx*fy;
+    v[4] = delx*fz;
+    v[5] = dely*fz;
+
+    if (vflag_global) {
+      if (newton_bond) {
+        virial[0] += v[0];
+        virial[1] += v[1];
+        virial[2] += v[2];
+        virial[3] += v[3];
+        virial[4] += v[4];
+        virial[5] += v[5];
+      } else {
+        if (i < nlocal) {
+          virial[0] += 0.5*v[0];
+          virial[1] += 0.5*v[1];
+          virial[2] += 0.5*v[2];
+          virial[3] += 0.5*v[3];
+          virial[4] += 0.5*v[4];
+          virial[5] += 0.5*v[5];
+        }
+        if (j < nlocal) {
+          virial[0] += 0.5*v[0];
+          virial[1] += 0.5*v[1];
+          virial[2] += 0.5*v[2];
+          virial[3] += 0.5*v[3];
+          virial[4] += 0.5*v[4];
+          virial[5] += 0.5*v[5];
+        }
+      }
+    }
+
+    if (vflag_atom) {
+      if (newton_bond || i < nlocal) {
+        vatom[i][0] += 0.5*v[0];
+        vatom[i][1] += 0.5*v[1];
+        vatom[i][2] += 0.5*v[2];
+        vatom[i][3] += 0.5*v[3];
+        vatom[i][4] += 0.5*v[4];
+        vatom[i][5] += 0.5*v[5];
+      }
+      if (newton_bond || j < nlocal) {
+        vatom[j][0] += 0.5*v[0];
+        vatom[j][1] += 0.5*v[1];
+        vatom[j][2] += 0.5*v[2];
+        vatom[j][3] += 0.5*v[3];
+        vatom[j][4] += 0.5*v[4];
+        vatom[j][5] += 0.5*v[5];
+      }
+    }
+  }
+}
+
 /* ----------------------------------------------------------------------
    compute function for oxDNA pair interactions
    s=sugar-phosphate backbone site, b=base site, st=stacking site
@@ -295,9 +383,6 @@ void PairOxdnaStk::compute(int eflag, int vflag)
     // early rejection criterium
     if (evdwl) {
 
-    // increment energy
-    if (evflag) ev_tally(a,b,nlocal,newton_bond,evdwl,0.0,0.0,0.0,0.0,0.0);
-
     df1 = DF1(r_st, epsilon_st[atype][btype], a_st[atype][btype], cut_st_0[atype][btype],
         cut_st_lc[atype][btype], cut_st_hc[atype][btype], cut_st_lo[atype][btype], cut_st_hi[atype][btype],
         b_st_lo[atype][btype], b_st_hi[atype][btype]);
@@ -366,7 +451,7 @@ void PairOxdnaStk::compute(int eflag, int vflag)
 
     }
 
-    // increment forces, torques and virial
+    // increment forces and torques
 
     if (newton_bond || a < nlocal) {
 
@@ -402,7 +487,12 @@ void PairOxdnaStk::compute(int eflag, int vflag)
 
     }
 
-    if (evflag) ev_tally(a,b,nlocal,newton_bond,0.0,0.0,fpair,delr_st[0],delr_st[1],delr_st[2]);
+    // increment energy and virial
+    // NOTE: The virial is calculated on the 'molecular' basis.
+    // (see G. Ciccotti and J.P. Ryckaert, Comp. Phys. Rep. 4, 345-392 (1986))
+
+    if (evflag) ev_tally_xyz(a,b,nlocal,newton_bond,evdwl,
+        delf[0],delf[1],delf[2],x[a][0]-x[b][0],x[a][1]-x[b][1],x[a][2]-x[b][2]);
 
     // force, torque and virial contribution for forces between backbone sites
 
@@ -444,7 +534,7 @@ void PairOxdnaStk::compute(int eflag, int vflag)
 
     }
 
-    // increment forces, torques and virial
+    // increment forces and torques
 
     if (newton_bond || a < nlocal) {
 
@@ -480,8 +570,9 @@ void PairOxdnaStk::compute(int eflag, int vflag)
 
     }
 
-    if (evflag) ev_tally(a,b,nlocal,newton_bond,0.0,0.0,fpair,delr_ss[0],delr_ss[1],delr_ss[2]);
+    // increment virial only
-
+    if (evflag) ev_tally_xyz(a,b,nlocal,newton_bond,0.0,
+        delf[0],delf[1],delf[2],x[a][0]-x[b][0],x[a][1]-x[b][1],x[a][2]-x[b][2]);
 
     // pure torques not expressible as r x f
 
@@ -656,11 +747,11 @@ void PairOxdnaStk::settings(int narg, char **/*arg*/)
    return temperature dependent oxDNA stacking strength
 ------------------------------------------------------------------------- */
 
-double PairOxdnaStk::stacking_strength(double T)
+double PairOxdnaStk::stacking_strength(double xi_st, double kappa_st, double T)
 {
   double eps;
 
-  eps = 1.3448 + 2.6568 * T;
+  eps = xi_st + kappa_st * T;
 
   return eps;
 }
@@ -673,7 +764,7 @@ void PairOxdnaStk::coeff(int narg, char **arg)
 {
   int count;
 
-  if (narg != 22) error->all(FLERR,"Incorrect args for pair coefficients in oxdna/stk");
+  if (narg != 24) error->all(FLERR,"Incorrect args for pair coefficients in oxdna/stk");
   if (!allocated) allocate();
 
   int ilo,ihi,jlo,jhi;
@@ -683,7 +774,7 @@ void PairOxdnaStk::coeff(int narg, char **arg)
   // stacking interaction
   count = 0;
 
-  double T, epsilon_st_one, a_st_one, b_st_lo_one, b_st_hi_one;
+  double T, epsilon_st_one, xi_st_one, kappa_st_one, a_st_one, b_st_lo_one, b_st_hi_one;
   double cut_st_0_one, cut_st_c_one, cut_st_lo_one, cut_st_hi_one;
   double cut_st_lc_one, cut_st_hc_one, tmp, shift_st_one;
 
@@ -706,27 +797,29 @@ void PairOxdnaStk::coeff(int narg, char **arg)
   if (strcmp(arg[2],"seqdep") == 0) seqdepflag = 1;
 
   T = force->numeric(FLERR,arg[3]);
-  epsilon_st_one = stacking_strength(T);
+  xi_st_one = force->numeric(FLERR,arg[4]);
+  kappa_st_one = force->numeric(FLERR,arg[5]);
+  epsilon_st_one = stacking_strength(xi_st_one, kappa_st_one, T);
 
-  a_st_one = force->numeric(FLERR,arg[4]);
+  a_st_one = force->numeric(FLERR,arg[6]);
-  cut_st_0_one = force->numeric(FLERR,arg[5]);
+  cut_st_0_one = force->numeric(FLERR,arg[7]);
-  cut_st_c_one = force->numeric(FLERR,arg[6]);
+  cut_st_c_one = force->numeric(FLERR,arg[8]);
-  cut_st_lo_one = force->numeric(FLERR,arg[7]);
+  cut_st_lo_one = force->numeric(FLERR,arg[9]);
-  cut_st_hi_one = force->numeric(FLERR,arg[8]);
+  cut_st_hi_one = force->numeric(FLERR,arg[10]);
 
-  a_st4_one = force->numeric(FLERR,arg[9]);
+  a_st4_one = force->numeric(FLERR,arg[11]);
-  theta_st4_0_one = force->numeric(FLERR,arg[10]);
+  theta_st4_0_one = force->numeric(FLERR,arg[12]);
-  dtheta_st4_ast_one = force->numeric(FLERR,arg[11]);
+  dtheta_st4_ast_one = force->numeric(FLERR,arg[13]);
-  a_st5_one = force->numeric(FLERR,arg[12]);
+  a_st5_one = force->numeric(FLERR,arg[14]);
-  theta_st5_0_one = force->numeric(FLERR,arg[13]);
+  theta_st5_0_one = force->numeric(FLERR,arg[15]);
-  dtheta_st5_ast_one = force->numeric(FLERR,arg[14]);
+  dtheta_st5_ast_one = force->numeric(FLERR,arg[16]);
-  a_st6_one = force->numeric(FLERR,arg[15]);
+  a_st6_one = force->numeric(FLERR,arg[17]);
-  theta_st6_0_one = force->numeric(FLERR,arg[16]);
+  theta_st6_0_one = force->numeric(FLERR,arg[18]);
-  dtheta_st6_ast_one = force->numeric(FLERR,arg[17]);
+  dtheta_st6_ast_one = force->numeric(FLERR,arg[19]);
-  a_st1_one = force->numeric(FLERR,arg[18]);
+  a_st1_one = force->numeric(FLERR,arg[20]);
-  cosphi_st1_ast_one = force->numeric(FLERR,arg[19]);
+  cosphi_st1_ast_one = force->numeric(FLERR,arg[21]);
-  a_st2_one = force->numeric(FLERR,arg[20]);
+  a_st2_one = force->numeric(FLERR,arg[22]);
-  cosphi_st2_ast_one = force->numeric(FLERR,arg[21]);
+  cosphi_st2_ast_one = force->numeric(FLERR,arg[23]);
 
   b_st_lo_one = 2*a_st_one*exp(-a_st_one*(cut_st_lo_one-cut_st_0_one))*
         2*a_st_one*exp(-a_st_one*(cut_st_lo_one-cut_st_0_one))*

src/USER-CGDNA/pair_oxdna_stk.h

@@ -14,6 +14,7 @@
 #ifdef PAIR_CLASS
 
 PairStyle(oxdna/stk,PairOxdnaStk)
+PairStyle(oxdna2/stk,PairOxdnaStk)
 
 #else
 
@@ -44,7 +45,7 @@ class PairOxdnaStk : public Pair {
 
  protected:
   // stacking interaction
-  virtual double stacking_strength(double);
+  double stacking_strength(double, double, double);
   double **epsilon_st, **a_st, **cut_st_0, **cut_st_c;
   double **cut_st_lo, **cut_st_hi;
   double **cut_st_lc, **cut_st_hc, **b_st_lo, **b_st_hi, **shift_st;
@@ -61,6 +62,7 @@ class PairOxdnaStk : public Pair {
   int seqdepflag;
 
   virtual void allocate();
+  void ev_tally_xyz(int, int, int, int, double, double, double, double, double, double, double);
 };
 
 }

src/USER-CGDNA/pair_oxdna_xstk.cpp

@@ -422,7 +422,11 @@ void PairOxdnaXstk::compute(int eflag, int vflag)
       }
 
       // increment energy and virial
-      if (evflag) ev_tally(a,b,nlocal,newton_pair,evdwl,0.0,fpair,delr_hb[0],delr_hb[1],delr_hb[2]);
+      // NOTE: The virial is calculated on the 'molecular' basis.
+      // (see G. Ciccotti and J.P. Ryckaert, Comp. Phys. Rep. 4, 345-392 (1986))
+
+      if (evflag) ev_tally_xyz(a,b,nlocal,newton_pair,evdwl,0.0,
+          delf[0],delf[1],delf[2],x[a][0]-x[b][0],x[a][1]-x[b][1],x[a][2]-x[b][2]);
 
       // pure torques not expressible as r x f