diff --git a/src/pair_smatb.cpp b/src/pair_smatb.cpp
index 5e2fd5839f..32e4e91a8d 100644
--- a/src/pair_smatb.cpp
+++ b/src/pair_smatb.cpp
@@ -1,58 +1,58 @@
+/* ----------------------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   https://www.lammps.org/, 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.
+------------------------------------------------------------------------- */
+
+/* ----------------------------------------------------------------------
+This pair style is written by Daniele Rapetti (iximiel@gmail.com)
+------------------------------------------------------------------------- */
+#include "pair_smatb.h"
+#include "atom.h"
+#include "comm.h"
+#include "error.h"
+#include "force.h"
+#include "memory.h"
+#include "neigh_list.h"
+#include "neighbor.h"
 #include <math.h>
 #include <stdio.h>
 #include <stdlib.h>
 #include <string.h>
-#include "pair_smatb.h"
-#include "atom.h"
-#include "force.h"
-#include "comm.h"
-#include "neighbor.h"
-#include "neigh_list.h"
-#include "memory.h"
-#include "error.h"
-
-//#define DR_DEBUG
-#ifdef DR_DEBUG
-#include <iostream>
-#include <iomanip>
-#endif //DR_DEBUG
 
 using namespace LAMMPS_NS;
 
 #define MAXLINE 1024
 
-PairSMATB::PairSMATB(LAMMPS *lmp)
-    : Pair(lmp),
-      nmax(0),
-      on_eb(nullptr),
-      r0(nullptr),
-      p(nullptr),
-      A(nullptr),
-      q(nullptr),
-      QSI(nullptr),
-      cutOffStart(nullptr),
-      cutOffEnd(nullptr),
-      cutOffEnd2(nullptr),
-      a3(nullptr), a4(nullptr), a5(nullptr),
-      x3(nullptr), x4(nullptr), x5(nullptr) {
-  single_enable = 0;           // 1 if single() routine exists
-  restartinfo   = 1;           // 1 if pair style writes restart info
-  respa_enable  = 0;           // 1 if inner/middle/outer rRESPA routines
-  one_coeff     = 0;           // 1 if allows only one coeff * * call
-  manybody_flag = 1;           // 1 if a manybody potential
-  no_virial_fdotr_compute = 0; // 1 if does not invoke virial_fdotr_compute()
-  writedata     = 1;           // 1 if writes coeffs to data file
-  ghostneigh    = 0;           // 1 if pair style needs neighbors of ghosts
+PairSMATB::PairSMATB(LAMMPS *lmp) :
+    Pair(lmp), nmax(0), on_eb(nullptr), r0(nullptr), p(nullptr), A(nullptr), q(nullptr),
+    QSI(nullptr), cutOffStart(nullptr), cutOffEnd(nullptr), cutOffEnd2(nullptr), a3(nullptr),
+    a4(nullptr), a5(nullptr), x3(nullptr), x4(nullptr), x5(nullptr)
+{
+  single_enable = 0;              // 1 if single() routine exists
+  restartinfo = 1;                // 1 if pair style writes restart info
+  respa_enable = 0;               // 1 if inner/middle/outer rRESPA routines
+  one_coeff = 0;                  // 1 if allows only one coeff * * call
+  manybody_flag = 1;              // 1 if a manybody potential
+  no_virial_fdotr_compute = 0;    // 1 if does not invoke virial_fdotr_compute()
+  writedata = 1;                  // 1 if writes coeffs to data file
+  ghostneigh = 0;                 // 1 if pair style needs neighbors of ghosts
 
   // set comm size needed by this Pair
   comm_forward = 1;
   comm_reverse = 1;
 }
 
-PairSMATB::~PairSMATB() {
-  if (copymode) {
-    return;
-  }
+PairSMATB::~PairSMATB()
+{
+  if (copymode) { return; }
   memory->destroy(on_eb);
   if (allocated) {
     memory->destroy(setflag);
@@ -75,17 +75,18 @@ PairSMATB::~PairSMATB() {
   }
 }
 
-void PairSMATB::compute(int eflag, int vflag) {//workhorse routine that computes pairwise interactions
+void PairSMATB::compute(int eflag, int vflag)
+{    //workhorse routine that computes pairwise interactions
   //eflag means compute energy
   //vflag means compute virial
-  int i,j,ii,jj,jnum,itype,jtype;
-  double xtmp,ytmp,ztmp,del[3],fpair;
-  double dijsq,dij;
+  int i, j, ii, jj, jnum, itype, jtype;
+  double xtmp, ytmp, ztmp, del[3], fpair;
+  double dijsq, dij;
   double espo, aexpp, qsiexpq, eb_i, Fb, Fr;
   double polyval, polyval2, polyval3, polyval4, polyval5;
   //sets up the flags for energy caclulations
   if (eflag || vflag) {
-    ev_setup(eflag,vflag);
+    ev_setup(eflag, vflag);
     eng_vdwl = 0;
   } else {
     evflag = vflag_fdotr = eflag_global = eflag_atom = 0;
@@ -96,14 +97,14 @@ void PairSMATB::compute(int eflag, int vflag) {//workhorse routine that computes
 
   if (atom->nmax > nmax) {
     nmax = atom->nmax;
-    memory->grow(on_eb,nmax,"pair_smatb:on_eb");
+    memory->grow(on_eb, nmax, "pair_smatb:on_eb");
   }
 
   double **x = atom->x;
   double **f = atom->f;
-  int *type  = atom->type;
+  int *type = atom->type;
   int nlocal = atom->nlocal;
-  int nall   = nlocal + atom->nghost;
+  int nall = nlocal + atom->nghost;
 
   int newton_pair = force->newton_pair;
 
@@ -122,39 +123,41 @@ void PairSMATB::compute(int eflag, int vflag) {//workhorse routine that computes
 
   //FIRST LOOP: CALCULATES the squared bounding energy and accumulate it in on_eb for each atom
   for (ii = 0; ii < inum; ++ii) {
-    i     = ilist[ii];
-    xtmp  = x[i][0];
-    ytmp  = x[i][1];
-    ztmp  = x[i][2];
+    i = ilist[ii];
+    xtmp = x[i][0];
+    ytmp = x[i][1];
+    ztmp = x[i][2];
     itype = type[i];
     jlist = firstneigh[i];
     jnum = numneigh[i];
 
     for (jj = 0; jj < jnum; ++jj) {
-      j  = jlist[jj];
+      j = jlist[jj];
       j &= NEIGHMASK;
       jtype = type[j];
       del[0] = xtmp - x[j][0];
       del[1] = ytmp - x[j][1];
       del[2] = ztmp - x[j][2];
-      dijsq  = del[0]*del[0] + del[1]*del[1] + del[2]*del[2];
+      dijsq = del[0] * del[0] + del[1] * del[1] + del[2] * del[2];
 
-      if ( dijsq < cutOffEnd2[itype][jtype] ) {
-	dij = sqrt (dijsq);
-	if ( dij < cutOffStart[itype][jtype] ) {
-	  qsiexpq = (QSI[itype][jtype]*QSI[itype][jtype]) * exp(2.0*q[itype][jtype]*(1.0 - dij/r0[itype][jtype]));
-	} else {
-	  polyval  = dij-cutOffEnd[itype][jtype];
-	  polyval3 = polyval*polyval*polyval;
-	  polyval4 = polyval3*polyval;
-	  polyval5 = polyval4*polyval;
-	  qsiexpq  = x5[itype][jtype]*polyval5+x4[itype][jtype]*polyval4+x3[itype][jtype]*polyval3;
-	  qsiexpq  = qsiexpq* qsiexpq;
-	}
-	on_eb[i]+=qsiexpq;
-	//if (newton_pair || j < nlocal) {
-	on_eb[j]+=qsiexpq;
-	//}
+      if (dijsq < cutOffEnd2[itype][jtype]) {
+        dij = sqrt(dijsq);
+        if (dij < cutOffStart[itype][jtype]) {
+          qsiexpq = (QSI[itype][jtype] * QSI[itype][jtype]) *
+              exp(2.0 * q[itype][jtype] * (1.0 - dij / r0[itype][jtype]));
+        } else {
+          polyval = dij - cutOffEnd[itype][jtype];
+          polyval3 = polyval * polyval * polyval;
+          polyval4 = polyval3 * polyval;
+          polyval5 = polyval4 * polyval;
+          qsiexpq = x5[itype][jtype] * polyval5 + x4[itype][jtype] * polyval4 +
+              x3[itype][jtype] * polyval3;
+          qsiexpq = qsiexpq * qsiexpq;
+        }
+        on_eb[i] += qsiexpq;
+        //if (newton_pair || j < nlocal) {
+        on_eb[j] += qsiexpq;
+        //}
       }
     }
   }
@@ -167,22 +170,18 @@ void PairSMATB::compute(int eflag, int vflag) {//workhorse routine that computes
   // the store the reciprocal in on_eb in order to not do it in the SECOND LOOP
 
   for (ii = 0; ii < inum; ++ii) {
-    i     = ilist[ii];
-    if (i < nlocal){
-      eb_i=sqrt(on_eb[i]);
-      if (eb_i!=0.0) {
-        on_eb[i]=1.0/eb_i;
+    i = ilist[ii];
+    if (i < nlocal) {
+      eb_i = sqrt(on_eb[i]);
+      if (eb_i != 0.0) {
+        on_eb[i] = 1.0 / eb_i;
       } else {
         on_eb[i] = 0.0;
       }
       //if needed the bounding energy is accumulated:
       if (eflag_either) {
-        if (eflag_atom) {
-          eatom[i] -= eb_i;
-        }
-        if (eflag_global) {
-          eng_vdwl -= eb_i;
-        }
+        if (eflag_atom) { eatom[i] -= eb_i; }
+        if (eflag_global) { eng_vdwl -= eb_i; }
       }
     }
   }
@@ -207,60 +206,63 @@ void PairSMATB::compute(int eflag, int vflag) {//workhorse routine that computes
       del[1] = ytmp - x[j][1];
       del[2] = ztmp - x[j][2];
 
-      dijsq = del[0]*del[0] + del[1]*del[1] + del[2]*del[2];
-      if ( dijsq < cutOffEnd2[itype][jtype] ) {
-	dij = sqrt (dijsq);
-	if ( dij < cutOffStart[itype][jtype] ) {
-	  espo    = 1.0 - dij/r0[itype][jtype];
-	  aexpp   = exp(p[itype][jtype]*espo)*A[itype][jtype];
-	  Fr      =  (2.0*aexpp)*(p[itype][jtype]/r0[itype][jtype]);
-	  qsiexpq = (QSI[itype][jtype]*QSI[itype][jtype]) * exp(2.0*q[itype][jtype]*espo);
-	  Fb      = -qsiexpq * q[itype][jtype]/r0[itype][jtype];
-	} else {
-	  polyval  = dij-cutOffEnd[itype][jtype];
-	  polyval2 = polyval*polyval;
-	  polyval3 = polyval2*polyval;
-	  polyval4 = polyval3*polyval;
-	  polyval5 = polyval4*polyval;
-	  aexpp    = a5[itype][jtype]*polyval5+a4[itype][jtype]*polyval4+a3[itype][jtype]*polyval3;
-	  Fr       = -2.0*(5.0*a5[itype][jtype]*polyval4+4.0*a4[itype][jtype]*polyval3+3.0*a3[itype][jtype]*polyval2);
-	  qsiexpq  = x5[itype][jtype]*polyval5+x4[itype][jtype]*polyval4+x3[itype][jtype]*polyval3;
-	  Fb       = ((5.0*x5[itype][jtype]*polyval4+4.0*x4[itype][jtype]*polyval3+3.0*x3[itype][jtype]*polyval2))*qsiexpq;
-	}
-	//if needed the repulsive energy is accumulated:
-	if (eflag_either) {
-	  if (eflag_atom) {
-	    eatom[i] += aexpp;
-	    if (newton_pair || j < nlocal) {
-	      eatom[j] += aexpp;
-	    }
-	  }      
-	  if (eflag_global) {
-	    if (newton_pair || j < nlocal) {
-	      eng_vdwl += 2.0*(aexpp);
-	    } else {
-	      eng_vdwl += aexpp;
-	    }
-	  }
-	}
-	//calculates the module of the pair energy between i and j
-	fpair = (Fb*(on_eb[i]+on_eb[j]) + Fr)/dij;
+      dijsq = del[0] * del[0] + del[1] * del[1] + del[2] * del[2];
+      if (dijsq < cutOffEnd2[itype][jtype]) {
+        dij = sqrt(dijsq);
+        if (dij < cutOffStart[itype][jtype]) {
+          espo = 1.0 - dij / r0[itype][jtype];
+          aexpp = exp(p[itype][jtype] * espo) * A[itype][jtype];
+          Fr = (2.0 * aexpp) * (p[itype][jtype] / r0[itype][jtype]);
+          qsiexpq = (QSI[itype][jtype] * QSI[itype][jtype]) * exp(2.0 * q[itype][jtype] * espo);
+          Fb = -qsiexpq * q[itype][jtype] / r0[itype][jtype];
+        } else {
+          polyval = dij - cutOffEnd[itype][jtype];
+          polyval2 = polyval * polyval;
+          polyval3 = polyval2 * polyval;
+          polyval4 = polyval3 * polyval;
+          polyval5 = polyval4 * polyval;
+          aexpp = a5[itype][jtype] * polyval5 + a4[itype][jtype] * polyval4 +
+              a3[itype][jtype] * polyval3;
+          Fr = -2.0 *
+              (5.0 * a5[itype][jtype] * polyval4 + 4.0 * a4[itype][jtype] * polyval3 +
+               3.0 * a3[itype][jtype] * polyval2);
+          qsiexpq = x5[itype][jtype] * polyval5 + x4[itype][jtype] * polyval4 +
+              x3[itype][jtype] * polyval3;
+          Fb = ((5.0 * x5[itype][jtype] * polyval4 + 4.0 * x4[itype][jtype] * polyval3 +
+                 3.0 * x3[itype][jtype] * polyval2)) *
+              qsiexpq;
+        }
+        //if needed the repulsive energy is accumulated:
+        if (eflag_either) {
+          if (eflag_atom) {
+            eatom[i] += aexpp;
+            if (newton_pair || j < nlocal) { eatom[j] += aexpp; }
+          }
+          if (eflag_global) {
+            if (newton_pair || j < nlocal) {
+              eng_vdwl += 2.0 * (aexpp);
+            } else {
+              eng_vdwl += aexpp;
+            }
+          }
+        }
+        //calculates the module of the pair energy between i and j
+        fpair = (Fb * (on_eb[i] + on_eb[j]) + Fr) / dij;
 
-	f[i][0] += del[0]*fpair;
-	f[i][1] += del[1]*fpair;
-	f[i][2] += del[2]*fpair;
-	if (newton_pair || j < nlocal) {
-	  f[j][0] -= del[0]*fpair;
-	  f[j][1] -= del[1]*fpair;
-	  f[j][2] -= del[2]*fpair;
-	}      
-	if (vflag_atom) {
-	  ev_tally(i, j, nlocal, newton_pair,
-		   0.0, 0.0,//Energy is tally'd in the other parts of the potential
-		   fpair, del[0], del[1], del[2]);
-	}
+        f[i][0] += del[0] * fpair;
+        f[i][1] += del[1] * fpair;
+        f[i][2] += del[2] * fpair;
+        if (newton_pair || j < nlocal) {
+          f[j][0] -= del[0] * fpair;
+          f[j][1] -= del[1] * fpair;
+          f[j][2] -= del[2] * fpair;
+        }
+        if (vflag_atom) {
+          ev_tally(i, j, nlocal, newton_pair, 0.0,
+                   0.0,    //Energy is tally'd in the other parts of the potential
+                   fpair, del[0], del[1], del[2]);
+        }
       }
-
     }
   }
   if (vflag_fdotr) virial_fdotr_compute();
@@ -270,43 +272,42 @@ void PairSMATB::compute(int eflag, int vflag) {//workhorse routine that computes
    global settings
    ------------------------------------------------------------------------- */
 
-void PairSMATB::settings(int narg, char **) {//reads the input script line with arguments you define
-  if (narg > 0) error->all(FLERR,"Illegal pair_style command: smatb accepts no options");
+void PairSMATB::settings(int narg, char **)
+{    //reads the input script line with arguments you define
+  if (narg > 0) error->all(FLERR, "Illegal pair_style command: smatb accepts no options");
 }
 
 /* ----------------------------------------------------------------------
    allocate all arrays
    ------------------------------------------------------------------------- */
 
-void PairSMATB::allocate() {
+void PairSMATB::allocate()
+{
   int n = atom->ntypes;
-  int natoms=atom->natoms;
+  int natoms = atom->natoms;
 
-  memory->create(setflag, n+1, n+1, "pair_smatb:setflag");
+  memory->create(setflag, n + 1, n + 1, "pair_smatb:setflag");
   for (int i = 1; i <= n; i++) {
-    for (int j = i; j <= n; j++) {
-      setflag[i][j] = 0;
-    }
+    for (int j = i; j <= n; j++) { setflag[i][j] = 0; }
   }
 
-  memory->create(cutsq, n+1, n+1, "pair_smatb:cutsq");
-
+  memory->create(cutsq, n + 1, n + 1, "pair_smatb:cutsq");
 
   //memory->create is needed to make a false nxn array on a n^2x1 line of data
-  memory->create(r0, n+1, n+1, "pair_smatb:r0");
-  memory->create(p, n+1, n+1, "pair_smatb:p");
-  memory->create(A, n+1, n+1, "pair_smatb:A");
-  memory->create(q, n+1, n+1, "pair_smatb:q");
-  memory->create(QSI, n+1, n+1, "pair_smatb:QSI");
-  memory->create(cutOffStart, n+1, n+1, "pair_smatb:cutOffStart");
-  memory->create(cutOffEnd, n+1, n+1, "pair_smatb:cutOffEnd");
-  memory->create(cutOffEnd2, n+1, n+1, "pair_smatb:cutOffEnd2");
-  memory->create(a3, n+1, n+1, "pair_smatb:a1");
-  memory->create(a4, n+1, n+1, "pair_smatb:a2");
-  memory->create(a5, n+1, n+1, "pair_smatb:a5");
-  memory->create(x3, n+1, n+1, "pair_smatb:x1");
-  memory->create(x4, n+1, n+1, "pair_smatb:x2");
-  memory->create(x5, n+1, n+1, "pair_smatb:x3");
+  memory->create(r0, n + 1, n + 1, "pair_smatb:r0");
+  memory->create(p, n + 1, n + 1, "pair_smatb:p");
+  memory->create(A, n + 1, n + 1, "pair_smatb:A");
+  memory->create(q, n + 1, n + 1, "pair_smatb:q");
+  memory->create(QSI, n + 1, n + 1, "pair_smatb:QSI");
+  memory->create(cutOffStart, n + 1, n + 1, "pair_smatb:cutOffStart");
+  memory->create(cutOffEnd, n + 1, n + 1, "pair_smatb:cutOffEnd");
+  memory->create(cutOffEnd2, n + 1, n + 1, "pair_smatb:cutOffEnd2");
+  memory->create(a3, n + 1, n + 1, "pair_smatb:a1");
+  memory->create(a4, n + 1, n + 1, "pair_smatb:a2");
+  memory->create(a5, n + 1, n + 1, "pair_smatb:a5");
+  memory->create(x3, n + 1, n + 1, "pair_smatb:x1");
+  memory->create(x4, n + 1, n + 1, "pair_smatb:x2");
+  memory->create(x5, n + 1, n + 1, "pair_smatb:x3");
 
   allocated = 1;
 }
@@ -315,28 +316,29 @@ void PairSMATB::allocate() {
    set coeffs for one or more type pairs
    ------------------------------------------------------------------------- */
 
-void PairSMATB::coeff(int narg, char **arg) {//set coefficients for one i,j type pair
-  if (!allocated) {
-    allocate();
-  }
+void PairSMATB::coeff(int narg, char **arg)
+{    //set coefficients for one i,j type pair
+  if (!allocated) { allocate(); }
   if (narg != 9) {
-    error->all(FLERR,"Incorrect args for pair coefficients:\n SMATB needs \"i j r0 p q A QSI CO_start CO_end\"");
+    error->all(
+        FLERR,
+        "Incorrect args for pair coefficients:\n SMATB needs \"i j r0 p q A QSI CO_start CO_end\"");
   }
-  int ilo,ihi,jlo,jhi;
-  utils::bounds(FLERR,arg[0],1,atom->ntypes,ilo,ihi,error);
-  utils::bounds(FLERR,arg[1],1,atom->ntypes,jlo,jhi,error);
+  int ilo, ihi, jlo, jhi;
+  utils::bounds(FLERR, arg[0], 1, atom->ntypes, ilo, ihi, error);
+  utils::bounds(FLERR, arg[1], 1, atom->ntypes, jlo, jhi, error);
   //reading parameters from input
-  double myr0          = utils::numeric(FLERR,arg[2],false,lmp),
-         myp           = utils::numeric(FLERR,arg[3],false,lmp),
-         myq           = utils::numeric(FLERR,arg[4],false,lmp),
-         myA           = utils::numeric(FLERR,arg[5],false,lmp),
-         myQSI         = utils::numeric(FLERR,arg[6],false,lmp),
-         mycutOffStart = utils::numeric(FLERR,arg[7],false,lmp),
-         mycutOffEnd   = utils::numeric(FLERR,arg[8],false,lmp);
+  double myr0 = utils::numeric(FLERR, arg[2], false, lmp),
+         myp = utils::numeric(FLERR, arg[3], false, lmp),
+         myq = utils::numeric(FLERR, arg[4], false, lmp),
+         myA = utils::numeric(FLERR, arg[5], false, lmp),
+         myQSI = utils::numeric(FLERR, arg[6], false, lmp),
+         mycutOffStart = utils::numeric(FLERR, arg[7], false, lmp),
+         mycutOffEnd = utils::numeric(FLERR, arg[8], false, lmp);
   int count = 0;
 
   for (int i = ilo; i <= ihi; i++) {
-    for (int j = MAX(jlo,i); j <= jhi; j++) {
+    for (int j = MAX(jlo, i); j <= jhi; j++) {
       r0[i][j] = myr0;
       p[i][j] = myp;
       A[i][j] = myA;
@@ -351,8 +353,7 @@ void PairSMATB::coeff(int narg, char **arg) {//set coefficients for one i,j type
     }
   }
 
-  if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients");
-
+  if (count == 0) error->all(FLERR, "Incorrect args for pair coefficients");
 }
 
 /* ----------------------------------------------------------------------
@@ -368,41 +369,42 @@ void PairSMATB::coeff(int narg, char **arg) {//set coefficients for one i,j type
    init for one type pair i,j and corresponding j,i
    ------------------------------------------------------------------------- */
 
-double PairSMATB::init_one(int i, int j) {//perform initialization for one i,j type pair
+double PairSMATB::init_one(int i, int j)
+{    //perform initialization for one i,j type pair
   if (setflag[i][j] == 0) {
     ///@todo implement smatb mixing rules
-    cutOffStart[i][j] = MIN(cutOffStart[i][i],cutOffStart[j][j]);
-    cutOffEnd[i][j] = MAX(cutOffEnd[i][i],cutOffEnd[j][j]);
+    cutOffStart[i][j] = MIN(cutOffStart[i][i], cutOffStart[j][j]);
+    cutOffEnd[i][j] = MAX(cutOffEnd[i][i], cutOffEnd[j][j]);
 
-    error->all(FLERR,"All pair coeffs are not set");
+    error->all(FLERR, "All pair coeffs are not set");
   }
   //calculating the polynomial linking to zero
   double es = cutOffEnd[i][j] - cutOffStart[i][j];
-  double es2 = es*es;
-  double es3 = es2*es;
+  double es2 = es * es;
+  double es3 = es2 * es;
 
   //variables for poly for p and A
-  double expp = A[i][j] * exp(p[i][j]*(1.-cutOffStart[i][j]/r0[i][j]));
-  double ap = -1./es3;
-  double bp = p[i][j]/(r0[i][j]*es2);
-  double cp = -(p[i][j]*p[i][j])/(es*r0[i][j]*r0[i][j]);
+  double expp = A[i][j] * exp(p[i][j] * (1. - cutOffStart[i][j] / r0[i][j]));
+  double ap = -1. / es3;
+  double bp = p[i][j] / (r0[i][j] * es2);
+  double cp = -(p[i][j] * p[i][j]) / (es * r0[i][j] * r0[i][j]);
 
-  a5[i][j]= expp * (12.*ap + 6.*bp + cp)/(2.*es2);
-  a4[i][j]= expp * (15.*ap + 7.*bp + cp)/es;
-  a3[i][j]= expp * (20.*ap + 8.*bp + cp)/2.;
+  a5[i][j] = expp * (12. * ap + 6. * bp + cp) / (2. * es2);
+  a4[i][j] = expp * (15. * ap + 7. * bp + cp) / es;
+  a3[i][j] = expp * (20. * ap + 8. * bp + cp) / 2.;
 
   //variables for poly for q and qsi
-  double expq = QSI[i][j]*exp(q[i][j]*(1.-cutOffStart[i][j]/r0[i][j]));
-  double aq = -1/es3;
-  double bq = q[i][j]/(es2*r0[i][j]);
-  double cq = -(q[i][j]*q[i][j])/(es*r0[i][j]*r0[i][j]);
+  double expq = QSI[i][j] * exp(q[i][j] * (1. - cutOffStart[i][j] / r0[i][j]));
+  double aq = -1 / es3;
+  double bq = q[i][j] / (es2 * r0[i][j]);
+  double cq = -(q[i][j] * q[i][j]) / (es * r0[i][j] * r0[i][j]);
+
+  x5[i][j] = expq * (12. * aq + 6. * bq + cq) / (2. * es2);
+  x4[i][j] = expq * (15. * aq + 7. * bq + cq) / es;
+  x3[i][j] = expq * (20. * aq + 8. * bq + cq) / 2.;
 
-  x5[i][j] = expq * (12.*aq + 6.*bq + cq)/(2.*es2);
-  x4[i][j] = expq * (15.*aq + 7.*bq + cq)/es;
-  x3[i][j] = expq * (20.*aq + 8.*bq + cq)/2.;
-  
   cutOffEnd2[i][j] = cutOffEnd[i][j] * cutOffEnd[i][j];
-  if ( i!=j ) {
+  if (i != j) {
     setflag[j][i] = 1;
     cutOffEnd2[j][i] = cutOffEnd2[i][j];
 
@@ -421,26 +423,15 @@ double PairSMATB::init_one(int i, int j) {//perform initialization for one i,j t
     x4[j][i] = x4[i][j];
     x5[j][i] = x5[i][j];
   }
-#ifdef DR_DEBUG
-  std::cout << i << " " << j <<std::endl;
-  std::cout<< "r0\tp\tq\tA\tXi\tCS\tCE\tCE^2\n";
-  std::cout <<r0[j][i]<<"\t"<< p[j][i]<<"\t"<< q[j][i]
-	    <<"\t"<< A[j][i]<<"\t"<< QSI[j][i]<<"\t"
-	    << cutOffStart[j][i]<<"\t"<< cutOffEnd[j][i]<<"\t"<< cutOffEnd2[j][i]<<"\n"
-	    << "Polynomial link\n"
-	    << "a: " << a5[i][j] <<"\t"<< a4[i][j]<<"\t" << a3[i][j]<<"\n"
-	    << "q: " << x5[i][j] <<"\t"<< x4[i][j]<<"\t" << x3[i][j]<<"\n";
-#endif //DR_DEBUG
   //returns cutOffEnd to calculate cutforce and cutsq
   return cutOffEnd[i][j];
 }
 
 /* ---------------------------------------------------------------------- */
 
-int PairSMATB::pack_forward_comm(int n, int *list, double *buf,
-				 int pbc_flag, int *pbc)
+int PairSMATB::pack_forward_comm(int n, int *list, double *buf, int pbc_flag, int *pbc)
 {
-  int i,j,m;
+  int i, j, m;
 
   m = 0;
   for (i = 0; i < n; ++i) {
@@ -452,33 +443,32 @@ int PairSMATB::pack_forward_comm(int n, int *list, double *buf,
 
 /* ---------------------------------------------------------------------- */
 
-void PairSMATB::unpack_forward_comm(int n, int first, double *buf) {
-  int i,m,last;
-  
+void PairSMATB::unpack_forward_comm(int n, int first, double *buf)
+{
+  int i, m, last;
+
   m = 0;
   last = first + n;
-  for (i = first; i < last; ++i) {
-    on_eb[i] = buf[m++];
-  }
+  for (i = first; i < last; ++i) { on_eb[i] = buf[m++]; }
 }
 
 /* ---------------------------------------------------------------------- */
 
-int PairSMATB::pack_reverse_comm(int n, int first, double *buf) {
-  int i,m,last;
-  
+int PairSMATB::pack_reverse_comm(int n, int first, double *buf)
+{
+  int i, m, last;
+
   m = 0;
   last = first + n;
-  for (i = first; i < last; ++i) {
-    buf[m++] = on_eb[i];
-  }
+  for (i = first; i < last; ++i) { buf[m++] = on_eb[i]; }
   return m;
 }
 
 /* ---------------------------------------------------------------------- */
 
-void PairSMATB::unpack_reverse_comm(int n, int *list, double *buf) {
-  int i,j,m;
+void PairSMATB::unpack_reverse_comm(int n, int *list, double *buf)
+{
+  int i, j, m;
   m = 0;
   for (i = 0; i < n; i++) {
     j = list[i];
@@ -487,93 +477,96 @@ void PairSMATB::unpack_reverse_comm(int n, int *list, double *buf) {
 }
 
 //write binary data of this simulation:
-void PairSMATB::write_restart_settings(FILE *fp) {
-  fwrite(&offset_flag,sizeof(int),1,fp);
-  fwrite(&mix_flag,   sizeof(int),1,fp);
-  fwrite(&tail_flag,  sizeof(int),1,fp);
+void PairSMATB::write_restart_settings(FILE *fp)
+{
+  fwrite(&offset_flag, sizeof(int), 1, fp);
+  fwrite(&mix_flag, sizeof(int), 1, fp);
+  fwrite(&tail_flag, sizeof(int), 1, fp);
 }
 
-void PairSMATB::read_restart_settings(FILE *fp) {
+void PairSMATB::read_restart_settings(FILE *fp)
+{
   int me = comm->me;
   size_t result;
   if (me == 0) {
-    result = fread(&offset_flag,sizeof(int),1,fp);
-    result = fread(&mix_flag,   sizeof(int),1,fp);
-    result = fread(&tail_flag,  sizeof(int),1,fp);
+    result = fread(&offset_flag, sizeof(int), 1, fp);
+    result = fread(&mix_flag, sizeof(int), 1, fp);
+    result = fread(&tail_flag, sizeof(int), 1, fp);
   }
-  MPI_Bcast(&offset_flag,  1,MPI_INT,0,world);
-  MPI_Bcast(&mix_flag,     1,MPI_INT,0,world);
-  MPI_Bcast(&tail_flag,    1,MPI_INT,0,world);
+  MPI_Bcast(&offset_flag, 1, MPI_INT, 0, world);
+  MPI_Bcast(&mix_flag, 1, MPI_INT, 0, world);
+  MPI_Bcast(&tail_flag, 1, MPI_INT, 0, world);
 }
 
-void PairSMATB::write_restart(FILE *fp) {
+void PairSMATB::write_restart(FILE *fp)
+{
   write_restart_settings(fp);
   //"I J r0 p q A QSI CO_start CO_end"
-  int i,j;
+  int i, j;
   for (i = 1; i <= atom->ntypes; i++)
     for (j = i; j <= atom->ntypes; j++) {
-      fwrite(&setflag[i][j],        sizeof(int),   1,fp);
+      fwrite(&setflag[i][j], sizeof(int), 1, fp);
       if (setflag[i][j]) {
-        fwrite(&r0[i][j],         sizeof(double),1,fp);
-        fwrite(&p[i][j],          sizeof(double),1,fp);
-        fwrite(&q[i][j],          sizeof(double),1,fp);
-        fwrite(&A[i][j],          sizeof(double),1,fp);
-        fwrite(&QSI[i][j],        sizeof(double),1,fp);
-        fwrite(&cutOffStart[i][j],sizeof(double),1,fp);
-        fwrite(&cutOffEnd[i][j],  sizeof(double),1,fp);
+        fwrite(&r0[i][j], sizeof(double), 1, fp);
+        fwrite(&p[i][j], sizeof(double), 1, fp);
+        fwrite(&q[i][j], sizeof(double), 1, fp);
+        fwrite(&A[i][j], sizeof(double), 1, fp);
+        fwrite(&QSI[i][j], sizeof(double), 1, fp);
+        fwrite(&cutOffStart[i][j], sizeof(double), 1, fp);
+        fwrite(&cutOffEnd[i][j], sizeof(double), 1, fp);
       }
     }
   //maybe we need to save also the values of the various polynomials
 }
 
-void PairSMATB::read_restart(FILE *fp) {
+void PairSMATB::read_restart(FILE *fp)
+{
   read_restart_settings(fp);
   allocate();
   size_t result;
 
-  int i,j;
+  int i, j;
   int me = comm->me;
   for (i = 1; i <= atom->ntypes; i++)
     for (j = i; j <= atom->ntypes; j++) {
-      if (me == 0) {
-        result = fread(&setflag[i][j],sizeof(int),1,fp);
-      }
-      MPI_Bcast(&setflag[i][j],1,MPI_INT,0,world);
+      if (me == 0) { result = fread(&setflag[i][j], sizeof(int), 1, fp); }
+      MPI_Bcast(&setflag[i][j], 1, MPI_INT, 0, world);
       if (setflag[i][j]) {
         if (me == 0) {
-          result = fread(&r0[i][j],         sizeof(double),1,fp);
-          result = fread(&p[i][j],          sizeof(double),1,fp);
-          result = fread(&q[i][j],          sizeof(double),1,fp);
-          result = fread(&A[i][j],          sizeof(double),1,fp);
-          result = fread(&QSI[i][j],        sizeof(double),1,fp);
-          result = fread(&cutOffStart[i][j],sizeof(double),1,fp);
-          result = fread(&cutOffEnd[i][j],  sizeof(double),1,fp);
+          result = fread(&r0[i][j], sizeof(double), 1, fp);
+          result = fread(&p[i][j], sizeof(double), 1, fp);
+          result = fread(&q[i][j], sizeof(double), 1, fp);
+          result = fread(&A[i][j], sizeof(double), 1, fp);
+          result = fread(&QSI[i][j], sizeof(double), 1, fp);
+          result = fread(&cutOffStart[i][j], sizeof(double), 1, fp);
+          result = fread(&cutOffEnd[i][j], sizeof(double), 1, fp);
         }
-        MPI_Bcast(&r0[i][j],         1,MPI_DOUBLE,0,world);
-        MPI_Bcast(&p[i][j],          1,MPI_DOUBLE,0,world);
-        MPI_Bcast(&q[i][j],          1,MPI_DOUBLE,0,world);
-        MPI_Bcast(&A[i][j],          1,MPI_DOUBLE,0,world);
-        MPI_Bcast(&QSI[i][j],        1,MPI_DOUBLE,0,world);
-        MPI_Bcast(&cutOffStart[i][j],1,MPI_DOUBLE,0,world);
-        MPI_Bcast(&cutOffEnd[i][j],  1,MPI_DOUBLE,0,world);
+        MPI_Bcast(&r0[i][j], 1, MPI_DOUBLE, 0, world);
+        MPI_Bcast(&p[i][j], 1, MPI_DOUBLE, 0, world);
+        MPI_Bcast(&q[i][j], 1, MPI_DOUBLE, 0, world);
+        MPI_Bcast(&A[i][j], 1, MPI_DOUBLE, 0, world);
+        MPI_Bcast(&QSI[i][j], 1, MPI_DOUBLE, 0, world);
+        MPI_Bcast(&cutOffStart[i][j], 1, MPI_DOUBLE, 0, world);
+        MPI_Bcast(&cutOffEnd[i][j], 1, MPI_DOUBLE, 0, world);
       }
     }
 }
 
-void PairSMATB::write_data(FILE *fp) {
+void PairSMATB::write_data(FILE *fp)
+{
   //smatb needs I J r0 p q A QSI CO_start CO_end
   for (int i = 1; i <= atom->ntypes; i++) {
-    fprintf(fp,"%d %g %g %g %g %g %g %g\n",
-            i, r0[i][i], p[i][i], q[i][i], A[i][i], QSI[i][i], cutOffStart[i][i], cutOffEnd[i][i]);
+    fprintf(fp, "%d %g %g %g %g %g %g %g\n", i, r0[i][i], p[i][i], q[i][i], A[i][i], QSI[i][i],
+            cutOffStart[i][i], cutOffEnd[i][i]);
   }
 }
 
-void PairSMATB::write_data_all(FILE *fp) {
+void PairSMATB::write_data_all(FILE *fp)
+{
   for (int i = 1; i <= atom->ntypes; i++) {
     for (int j = i; j <= atom->ntypes; j++) {
-      fprintf(fp,"%d %d %g %g %g %g %g %g %g\n",
-              i, j, r0[i][j], p[i][j], q[i][j], A[i][j], QSI[i][j], cutOffStart[i][j], cutOffEnd[i][j]);
+      fprintf(fp, "%d %d %g %g %g %g %g %g %g\n", i, j, r0[i][j], p[i][j], q[i][j], A[i][j],
+              QSI[i][j], cutOffStart[i][j], cutOffEnd[i][j]);
     }
   }
 }
-
diff --git a/src/pair_smatb.h b/src/pair_smatb.h
index 3013a1d40e..8a5fdea91f 100644
--- a/src/pair_smatb.h
+++ b/src/pair_smatb.h
@@ -26,44 +26,45 @@ PairStyle(smatb,PairSMATB)
 
 namespace LAMMPS_NS {
 
-
-  class PairSMATB : public Pair {
-  public:
-    PairSMATB(class LAMMPS *);
-    virtual ~PairSMATB();
-    void compute(int, int);//workhorse routine that computes pairwise interactions
-    /*
+class PairSMATB : public Pair {
+ public:
+  PairSMATB(class LAMMPS *);
+  virtual ~PairSMATB();
+  void compute(int, int);    //workhorse routine that computes pairwise interactions
+  /*
     void compute_inner();
     void compute_middle();
-    void compute_outer(int, int);*/
-    void settings(int, char **);//reads the input script line with arguments you define
-    void coeff(int, char **);//set coefficients for one i,j type pair
-    //void init_style();//initialization specific to this pair style
-    double init_one(int, int);//perform initialization for one i,j type pair
-    //double single(int, int, int, int, double, double, double, double &);//force and energy of a single pairwise interaction between 2 atoms    
+    void compute_outer(int, int);
+    */
+  void settings(int, char **);    //reads the input script line with arguments you define
+  void coeff(int, char **);       //set coefficients for one i,j type pair
+  //void init_style();//initialization specific to this pair style
+  double init_one(int, int);    //perform initialization for one i,j type pair
+  //double single(int, int, int, int, double, double, double, double &);//force and energy of a single pairwise interaction between 2 atoms
 
-    void write_restart(FILE *);
-    void read_restart(FILE *);
-    void write_restart_settings(FILE *);
-    void read_restart_settings(FILE *);
-    void write_data(FILE *);
-    void write_data_all(FILE *);
+  void write_restart(FILE *);
+  void read_restart(FILE *);
+  void write_restart_settings(FILE *);
+  void read_restart_settings(FILE *);
+  void write_data(FILE *);
+  void write_data_all(FILE *);
 
-    virtual int    pack_forward_comm(int, int *, double *, int, int *);
-    virtual void unpack_forward_comm(int, int, double *);
-    virtual int    pack_reverse_comm(int, int, double *);
-    virtual void unpack_reverse_comm(int, int *, double *);
-  protected: 
-    virtual void allocate();
-    int nmax;                    // allocated size of per-atom arrays
-    double *on_eb;               //allocated to store up caclulation values
-    double **r0;                 // interaction radius, user-given
-    double **p, **A, **q, **QSI; // parameters user-given
-    double **cutOffStart, **cutOffEnd;//cut offs, user given
-    double **cutOffEnd2;         //squared cut off end, calculated
-    double **a3, **a4, **a5;     //polynomial for cutoff linking to zero:   Ae^p substitution
-    double **x3, **x4, **x5;     //polynomial for cutoff linking to zero: QSIe^q substitution
-    /* latex form of the potential (R_c is cutOffEnd, \Xi is QSI):
+  virtual int pack_forward_comm(int, int *, double *, int, int *);
+  virtual void unpack_forward_comm(int, int, double *);
+  virtual int pack_reverse_comm(int, int, double *);
+  virtual void unpack_reverse_comm(int, int *, double *);
+
+ protected:
+  virtual void allocate();
+  int nmax;                             // allocated size of per-atom arrays
+  double *on_eb;                        //allocated to store up caclulation values
+  double **r0;                          // interaction radius, user-given
+  double **p, **A, **q, **QSI;          // parameters user-given
+  double **cutOffStart, **cutOffEnd;    //cut offs, user given
+  double **cutOffEnd2;                  //squared cut off end, calculated
+  double **a3, **a4, **a5;              //polynomial for cutoff linking to zero:   Ae^p substitution
+  double **x3, **x4, **x5;              //polynomial for cutoff linking to zero: QSIe^q substitution
+  /* latex form of the potential (R_c is cutOffEnd, \Xi is QSI):
 
        E_i =
        \sum_{j,R_{ij}\leq R_c} A  e^{-p \lrt{\frac{R_{ij}}{R_{0}}-1}}
@@ -71,14 +72,14 @@ namespace LAMMPS_NS {
 
        NB::this form does not have the polynomial link to 0 for the cut off
      */
-  };
+};
 
-}
+}    // namespace LAMMPS_NS
 
 #endif
 #endif
 
-/* ERROR/WARNING messages:
+    /* ERROR/WARNING messages:
 
    E: Illegal ... command
 
@@ -90,9 +91,4 @@ namespace LAMMPS_NS {
 
    Self-explanatory.  Check the input script or data file.
 
-   E: Cannot open EAM potential file %s
-
-   The specified EAM potential file cannot be opened.  Check that the
-   path and name are correct.
-
 */
diff --git a/src/pair_smatb_single.cpp b/src/pair_smatb_single.cpp
index 3430ec1e73..f5b0e340bc 100644
--- a/src/pair_smatb_single.cpp
+++ b/src/pair_smatb_single.cpp
@@ -1,58 +1,57 @@
+/* ----------------------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   https://www.lammps.org/, 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.
+------------------------------------------------------------------------- */
+
+/* ----------------------------------------------------------------------
+This pair style is written by Daniele Rapetti (iximiel@gmail.com)
+------------------------------------------------------------------------- */
+#include "pair_smatb_single.h"
+#include "atom.h"
+#include "comm.h"
+#include "error.h"
+#include "force.h"
+#include "memory.h"
+#include "neigh_list.h"
+#include "neighbor.h"
 #include <math.h>
 #include <stdio.h>
 #include <stdlib.h>
 #include <string.h>
-#include "pair_smatb_single.h"
-#include "atom.h"
-#include "force.h"
-#include "comm.h"
-#include "neighbor.h"
-#include "neigh_list.h"
-#include "memory.h"
-#include "error.h"
-
-//#define DR_DEBUG
-#ifdef DR_DEBUG
-#include <iostream>
-#include <iomanip>
-#endif //DR_DEBUG
 
 using namespace LAMMPS_NS;
 
 #define MAXLINE 1024
 
-PairSMATBSingle::PairSMATBSingle(LAMMPS *lmp)
-    : Pair(lmp),
-      nmax(0),
-      on_eb(nullptr),
-      r0(0),
-      p(0),
-      A(0),
-      q(0),
-      QSI(0),
-      cutOffStart(0),
-      cutOffEnd(0),
-      cutOffEnd2(0),
-      a3(0), a4(0), a5(0),
-      x3(0), x4(0), x5(0) {
-  single_enable = 0;           // 1 if single() routine exists
-  restartinfo   = 1;           // 1 if pair style writes restart info
-  respa_enable  = 0;           // 1 if inner/middle/outer rRESPA routines
-  one_coeff     = 0;           // 1 if allows only one coeff * * call
-  manybody_flag = 1;           // 1 if a manybody potential
-  no_virial_fdotr_compute = 0; // 1 if does not invoke virial_fdotr_compute()
-  writedata     = 1;           // 1 if writes coeffs to data file
-  ghostneigh    = 0;           // 1 if pair style needs neighbors of ghosts
+PairSMATBSingle::PairSMATBSingle(LAMMPS *lmp) :
+    Pair(lmp), nmax(0), on_eb(nullptr), r0(0), p(0), A(0), q(0), QSI(0), cutOffStart(0),
+    cutOffEnd(0), cutOffEnd2(0), a3(0), a4(0), a5(0), x3(0), x4(0), x5(0)
+{
+  single_enable = 0;              // 1 if single() routine exists
+  restartinfo = 1;                // 1 if pair style writes restart info
+  respa_enable = 0;               // 1 if inner/middle/outer rRESPA routines
+  one_coeff = 0;                  // 1 if allows only one coeff * * call
+  manybody_flag = 1;              // 1 if a manybody potential
+  no_virial_fdotr_compute = 0;    // 1 if does not invoke virial_fdotr_compute()
+  writedata = 1;                  // 1 if writes coeffs to data file
+  ghostneigh = 0;                 // 1 if pair style needs neighbors of ghosts
 
   // set comm size needed by this Pair
   comm_forward = 1;
   comm_reverse = 1;
 }
 
-PairSMATBSingle::~PairSMATBSingle() {
-  if (copymode) {
-    return;
-  }
+PairSMATBSingle::~PairSMATBSingle()
+{
+  if (copymode) { return; }
   memory->destroy(on_eb);
   if (allocated) {
     memory->destroy(setflag);
@@ -60,17 +59,18 @@ PairSMATBSingle::~PairSMATBSingle() {
   }
 }
 
-void PairSMATBSingle::compute(int eflag, int vflag) {//workhorse routine that computes pairwise interactions
+void PairSMATBSingle::compute(int eflag, int vflag)
+{    //workhorse routine that computes pairwise interactions
   //eflag means compute energy
   //vflag means compute virial
-  int i,j,ii,jj,jnum;//,itype,jtype;
-  double xtmp,ytmp,ztmp,del[3],fpair;
-  double dijsq,dij;
+  int i, j, ii, jj, jnum;    //,itype,jtype;
+  double xtmp, ytmp, ztmp, del[3], fpair;
+  double dijsq, dij;
   double espo, aexpp, qsiexpq, eb_i, Fb, Fr;
   double polyval, polyval2, polyval3, polyval4, polyval5;
   //sets up the flags for energy caclulations
   if (eflag || vflag) {
-    ev_setup(eflag,vflag);
+    ev_setup(eflag, vflag);
     eng_vdwl = 0;
   } else {
     evflag = vflag_fdotr = eflag_global = eflag_atom = 0;
@@ -81,14 +81,14 @@ void PairSMATBSingle::compute(int eflag, int vflag) {//workhorse routine that co
 
   if (atom->nmax > nmax) {
     nmax = atom->nmax;
-    memory->grow(on_eb,nmax,"pair_smatb:on_eb");
+    memory->grow(on_eb, nmax, "pair_smatb:on_eb");
   }
 
   double **x = atom->x;
   double **f = atom->f;
-  int *type  = atom->type;
+  int *type = atom->type;
   int nlocal = atom->nlocal;
-  int nall   = nlocal + atom->nghost;
+  int nall = nlocal + atom->nghost;
 
   int newton_pair = force->newton_pair;
 
@@ -107,39 +107,39 @@ void PairSMATBSingle::compute(int eflag, int vflag) {//workhorse routine that co
 
   //FIRST LOOP: CALCULATES the squared bounding energy and accumulate it in on_eb for each atom
   for (ii = 0; ii < inum; ++ii) {
-    i     = ilist[ii];
-    xtmp  = x[i][0];
-    ytmp  = x[i][1];
-    ztmp  = x[i][2];
+    i = ilist[ii];
+    xtmp = x[i][0];
+    ytmp = x[i][1];
+    ztmp = x[i][2];
     //itype = type[i];
     jlist = firstneigh[i];
     jnum = numneigh[i];
 
     for (jj = 0; jj < jnum; ++jj) {
-      j  = jlist[jj];
+      j = jlist[jj];
       j &= NEIGHMASK;
       //jtype = type[j];
       del[0] = xtmp - x[j][0];
       del[1] = ytmp - x[j][1];
       del[2] = ztmp - x[j][2];
-      dijsq  = del[0]*del[0] + del[1]*del[1] + del[2]*del[2];
+      dijsq = del[0] * del[0] + del[1] * del[1] + del[2] * del[2];
 
-      if ( dijsq < cutOffEnd2 ) {
-	dij = sqrt (dijsq);
-	if ( dij < cutOffStart ) {
-	  qsiexpq = (QSI*QSI) * exp(2.0*q*(1.0 - dij/r0));
-	} else {
-	  polyval  = dij-cutOffEnd;
-	  polyval3 = polyval*polyval*polyval;
-	  polyval4 = polyval3*polyval;
-	  polyval5 = polyval4*polyval;
-	  qsiexpq  = x5*polyval5+x4*polyval4+x3*polyval3;
-	  qsiexpq  = qsiexpq* qsiexpq;
-	}
-	on_eb[i]+=qsiexpq;
-	//if (newton_pair || j < nlocal) {
-	on_eb[j]+=qsiexpq;
-	//}
+      if (dijsq < cutOffEnd2) {
+        dij = sqrt(dijsq);
+        if (dij < cutOffStart) {
+          qsiexpq = (QSI * QSI) * exp(2.0 * q * (1.0 - dij / r0));
+        } else {
+          polyval = dij - cutOffEnd;
+          polyval3 = polyval * polyval * polyval;
+          polyval4 = polyval3 * polyval;
+          polyval5 = polyval4 * polyval;
+          qsiexpq = x5 * polyval5 + x4 * polyval4 + x3 * polyval3;
+          qsiexpq = qsiexpq * qsiexpq;
+        }
+        on_eb[i] += qsiexpq;
+        //if (newton_pair || j < nlocal) {
+        on_eb[j] += qsiexpq;
+        //}
       }
     }
   }
@@ -152,22 +152,18 @@ void PairSMATBSingle::compute(int eflag, int vflag) {//workhorse routine that co
   // the store the reciprocal in on_eb in order to not do it in the SECOND LOOP
 
   for (ii = 0; ii < inum; ++ii) {
-    i     = ilist[ii];
-    if (i < nlocal){
-      eb_i=sqrt(on_eb[i]);
-      if (eb_i!=0.0) {
-        on_eb[i]=1.0/eb_i;
+    i = ilist[ii];
+    if (i < nlocal) {
+      eb_i = sqrt(on_eb[i]);
+      if (eb_i != 0.0) {
+        on_eb[i] = 1.0 / eb_i;
       } else {
         on_eb[i] = 0.0;
       }
       //if needed the bounding energy is accumulated:
       if (eflag_either) {
-        if (eflag_atom) {
-          eatom[i] -= eb_i;
-        }
-        if (eflag_global) {
-          eng_vdwl -= eb_i;
-        }
+        if (eflag_atom) { eatom[i] -= eb_i; }
+        if (eflag_global) { eng_vdwl -= eb_i; }
       }
     }
   }
@@ -192,60 +188,57 @@ void PairSMATBSingle::compute(int eflag, int vflag) {//workhorse routine that co
       del[1] = ytmp - x[j][1];
       del[2] = ztmp - x[j][2];
 
-      dijsq = del[0]*del[0] + del[1]*del[1] + del[2]*del[2];
-      if ( dijsq < cutOffEnd2 ) {
-	dij = sqrt (dijsq);
-	if ( dij < cutOffStart ) {
-	  espo    = 1.0 - dij/r0;
-	  aexpp   = exp(p*espo)*A;
-	  Fr      =  (2.0*aexpp)*(p/r0);
-	  qsiexpq = (QSI*QSI) * exp(2.0*q*espo);
-	  Fb      = -qsiexpq * q/r0;
-	} else {
-	  polyval  = dij-cutOffEnd;
-	  polyval2 = polyval*polyval;
-	  polyval3 = polyval2*polyval;
-	  polyval4 = polyval3*polyval;
-	  polyval5 = polyval4*polyval;
-	  aexpp    = a5*polyval5+a4*polyval4+a3*polyval3;
-	  Fr       = -2.0*(5.0*a5*polyval4+4.0*a4*polyval3+3.0*a3*polyval2);
-	  qsiexpq  = x5*polyval5+x4*polyval4+x3*polyval3;
-	  Fb       = ((5.0*x5*polyval4+4.0*x4*polyval3+3.0*x3*polyval2))*qsiexpq;
-	}
-	//if needed the repulsive energy is accumulated:
-	if (eflag_either) {
-	  if (eflag_atom) {
-	    eatom[i] += aexpp;
-	    if (newton_pair || j < nlocal) {
-	      eatom[j] += aexpp;
-	    }
-	  }      
-	  if (eflag_global) {
-	    if (newton_pair || j < nlocal) {
-	      eng_vdwl += 2.0*(aexpp);
-	    } else {
-	      eng_vdwl += aexpp;
-	    }
-	  }
-	}
-	//calculates the module of the pair energy between i and j
-	fpair = (Fb*(on_eb[i]+on_eb[j]) + Fr)/dij;
+      dijsq = del[0] * del[0] + del[1] * del[1] + del[2] * del[2];
+      if (dijsq < cutOffEnd2) {
+        dij = sqrt(dijsq);
+        if (dij < cutOffStart) {
+          espo = 1.0 - dij / r0;
+          aexpp = exp(p * espo) * A;
+          Fr = (2.0 * aexpp) * (p / r0);
+          qsiexpq = (QSI * QSI) * exp(2.0 * q * espo);
+          Fb = -qsiexpq * q / r0;
+        } else {
+          polyval = dij - cutOffEnd;
+          polyval2 = polyval * polyval;
+          polyval3 = polyval2 * polyval;
+          polyval4 = polyval3 * polyval;
+          polyval5 = polyval4 * polyval;
+          aexpp = a5 * polyval5 + a4 * polyval4 + a3 * polyval3;
+          Fr = -2.0 * (5.0 * a5 * polyval4 + 4.0 * a4 * polyval3 + 3.0 * a3 * polyval2);
+          qsiexpq = x5 * polyval5 + x4 * polyval4 + x3 * polyval3;
+          Fb = ((5.0 * x5 * polyval4 + 4.0 * x4 * polyval3 + 3.0 * x3 * polyval2)) * qsiexpq;
+        }
+        //if needed the repulsive energy is accumulated:
+        if (eflag_either) {
+          if (eflag_atom) {
+            eatom[i] += aexpp;
+            if (newton_pair || j < nlocal) { eatom[j] += aexpp; }
+          }
+          if (eflag_global) {
+            if (newton_pair || j < nlocal) {
+              eng_vdwl += 2.0 * (aexpp);
+            } else {
+              eng_vdwl += aexpp;
+            }
+          }
+        }
+        //calculates the module of the pair energy between i and j
+        fpair = (Fb * (on_eb[i] + on_eb[j]) + Fr) / dij;
 
-	f[i][0] += del[0]*fpair;
-	f[i][1] += del[1]*fpair;
-	f[i][2] += del[2]*fpair;
-	if (newton_pair || j < nlocal) {
-	  f[j][0] -= del[0]*fpair;
-	  f[j][1] -= del[1]*fpair;
-	  f[j][2] -= del[2]*fpair;
-	}      
-	if (vflag_atom) {
-	  ev_tally(i, j, nlocal, newton_pair,
-		   0.0, 0.0,//Energy is tally'd in the other parts of the potential
-		   fpair, del[0], del[1], del[2]);
-	}
+        f[i][0] += del[0] * fpair;
+        f[i][1] += del[1] * fpair;
+        f[i][2] += del[2] * fpair;
+        if (newton_pair || j < nlocal) {
+          f[j][0] -= del[0] * fpair;
+          f[j][1] -= del[1] * fpair;
+          f[j][2] -= del[2] * fpair;
+        }
+        if (vflag_atom) {
+          ev_tally(i, j, nlocal, newton_pair, 0.0,
+                   0.0,    //Energy is tally'd in the other parts of the potential
+                   fpair, del[0], del[1], del[2]);
+        }
       }
-
     }
   }
   if (vflag_fdotr) virial_fdotr_compute();
@@ -255,26 +248,26 @@ void PairSMATBSingle::compute(int eflag, int vflag) {//workhorse routine that co
    global settings
    ------------------------------------------------------------------------- */
 
-void PairSMATBSingle::settings(int narg, char **) {//reads the input script line with arguments you define
-  if (narg > 0) error->all(FLERR,"Illegal pair_style command: smatb accepts no options");
+void PairSMATBSingle::settings(int narg, char **)
+{    //reads the input script line with arguments you define
+  if (narg > 0) error->all(FLERR, "Illegal pair_style command: smatb accepts no options");
 }
 
 /* ----------------------------------------------------------------------
    allocate all arrays
    ------------------------------------------------------------------------- */
 
-void PairSMATBSingle::allocate() {
+void PairSMATBSingle::allocate()
+{
   int n = atom->ntypes;
-  int natoms=atom->natoms;
+  int natoms = atom->natoms;
 
-  memory->create(setflag, n+1, n+1, "pair_smatb:setflag");
+  memory->create(setflag, n + 1, n + 1, "pair_smatb:setflag");
   for (int i = 1; i <= n; i++) {
-    for (int j = i; j <= n; j++) {
-      setflag[i][j] = 0;
-    }
+    for (int j = i; j <= n; j++) { setflag[i][j] = 0; }
   }
 
-  memory->create(cutsq, n+1, n+1, "pair_smatb:cutsq");
+  memory->create(cutsq, n + 1, n + 1, "pair_smatb:cutsq");
 
   allocated = 1;
 }
@@ -283,28 +276,29 @@ void PairSMATBSingle::allocate() {
    set coeffs for one or more type pairs
    ------------------------------------------------------------------------- */
 
-void PairSMATBSingle::coeff(int narg, char **arg) {//set coefficients for one i,j type pair
-  if (!allocated) {
-    allocate();
-  }
+void PairSMATBSingle::coeff(int narg, char **arg)
+{    //set coefficients for one i,j type pair
+  if (!allocated) { allocate(); }
   if (narg != 9) {
-    error->all(FLERR,"Incorrect args for pair coefficients:\n SMATB needs \"i j r0 p q A QSI CO_start CO_end\"");
+    error->all(
+        FLERR,
+        "Incorrect args for pair coefficients:\n SMATB needs \"i j r0 p q A QSI CO_start CO_end\"");
   }
-  int ilo,ihi,jlo,jhi;
-  utils::bounds(FLERR,arg[0],1,atom->ntypes,ilo,ihi,error);
-  utils::bounds(FLERR,arg[1],1,atom->ntypes,jlo,jhi,error);
+  int ilo, ihi, jlo, jhi;
+  utils::bounds(FLERR, arg[0], 1, atom->ntypes, ilo, ihi, error);
+  utils::bounds(FLERR, arg[1], 1, atom->ntypes, jlo, jhi, error);
   //reading parameters from input
-  double myr0          = utils::numeric(FLERR,arg[2],false,lmp),
-         myp           = utils::numeric(FLERR,arg[3],false,lmp),
-         myq           = utils::numeric(FLERR,arg[4],false,lmp),
-         myA           = utils::numeric(FLERR,arg[5],false,lmp),
-         myQSI         = utils::numeric(FLERR,arg[6],false,lmp),
-         mycutOffStart = utils::numeric(FLERR,arg[7],false,lmp),
-         mycutOffEnd   = utils::numeric(FLERR,arg[8],false,lmp);
+  double myr0 = utils::numeric(FLERR, arg[2], false, lmp),
+         myp = utils::numeric(FLERR, arg[3], false, lmp),
+         myq = utils::numeric(FLERR, arg[4], false, lmp),
+         myA = utils::numeric(FLERR, arg[5], false, lmp),
+         myQSI = utils::numeric(FLERR, arg[6], false, lmp),
+         mycutOffStart = utils::numeric(FLERR, arg[7], false, lmp),
+         mycutOffEnd = utils::numeric(FLERR, arg[8], false, lmp);
   int count = 0;
 
   for (int i = ilo; i <= ihi; i++) {
-    for (int j = MAX(jlo,i); j <= jhi; j++) {
+    for (int j = MAX(jlo, i); j <= jhi; j++) {
       r0 = myr0;
       p = myp;
       A = myA;
@@ -318,8 +312,7 @@ void PairSMATBSingle::coeff(int narg, char **arg) {//set coefficients for one i,
     }
   }
 
-  if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients");
-
+  if (count == 0) error->all(FLERR, "Incorrect args for pair coefficients");
 }
 
 /* ----------------------------------------------------------------------
@@ -335,35 +328,36 @@ void PairSMATBSingle::coeff(int narg, char **arg) {//set coefficients for one i,
    init for one type pair i,j and corresponding j,i
    ------------------------------------------------------------------------- */
 
-double PairSMATBSingle::init_one(int i, int j) {//perform initialization for one i,j type pair
-  if (setflag[i][j] == 0) error->all(FLERR,"All pair coeffs are not set");
+double PairSMATBSingle::init_one(int i, int j)
+{    //perform initialization for one i,j type pair
+  if (setflag[i][j] == 0) error->all(FLERR, "All pair coeffs are not set");
   //calculating the polynomial linking to zero
   double es = cutOffEnd - cutOffStart;
-  double es2 = es*es;
-  double es3 = es2*es;
+  double es2 = es * es;
+  double es3 = es2 * es;
 
   //variables for poly for p and A
-  double expp = A * exp(p*(1.-cutOffStart/r0));
-  double ap = -1./es3;
-  double bp = p/(r0*es2);
-  double cp = -(p*p)/(es*r0*r0);
+  double expp = A * exp(p * (1. - cutOffStart / r0));
+  double ap = -1. / es3;
+  double bp = p / (r0 * es2);
+  double cp = -(p * p) / (es * r0 * r0);
 
-  a5= expp * (12.*ap + 6.*bp + cp)/(2.*es2);
-  a4= expp * (15.*ap + 7.*bp + cp)/es;
-  a3= expp * (20.*ap + 8.*bp + cp)/2.;
+  a5 = expp * (12. * ap + 6. * bp + cp) / (2. * es2);
+  a4 = expp * (15. * ap + 7. * bp + cp) / es;
+  a3 = expp * (20. * ap + 8. * bp + cp) / 2.;
 
   //variables for poly for q and qsi
-  double expq = QSI*exp(q*(1.-cutOffStart/r0));
-  double aq = -1/es3;
-  double bq = q/(es2*r0);
-  double cq = -(q*q)/(es*r0*r0);
+  double expq = QSI * exp(q * (1. - cutOffStart / r0));
+  double aq = -1 / es3;
+  double bq = q / (es2 * r0);
+  double cq = -(q * q) / (es * r0 * r0);
+
+  x5 = expq * (12. * aq + 6. * bq + cq) / (2. * es2);
+  x4 = expq * (15. * aq + 7. * bq + cq) / es;
+  x3 = expq * (20. * aq + 8. * bq + cq) / 2.;
 
-  x5 = expq * (12.*aq + 6.*bq + cq)/(2.*es2);
-  x4 = expq * (15.*aq + 7.*bq + cq)/es;
-  x3 = expq * (20.*aq + 8.*bq + cq)/2.;
-  
   cutOffEnd2 = cutOffEnd * cutOffEnd;
-  if ( i!=j ) {
+  if (i != j) {
     setflag[j][i] = 1;
     cutOffEnd2 = cutOffEnd2;
 
@@ -382,26 +376,15 @@ double PairSMATBSingle::init_one(int i, int j) {//perform initialization for one
     x4 = x4;
     x5 = x5;
   }
-#ifdef DR_DEBUG
-  std::cout << i << " " << j <<std::endl;
-  std::cout<< "r0\tp\tq\tA\tXi\tCS\tCE\n";
-  std::cout <<r0<<"\t"<< p<<"\t"<< q
-	    <<"\t"<< A<<"\t"<< QSI<<"\t"
-	    << cutOffStart<<"\t"<< cutOffEnd<<"\n"
-	    << "Polynomial link\n"
-	    << "a: " << a5 <<"\t"<< a4<<"\t" << a3<<"\n"
-	    << "q: " << x5 <<"\t"<< x4<<"\t" << x3<<"\n";
-#endif //DR_DEBUG
   //returns cutOffEnd to calculate cutforce and cutsq
   return cutOffEnd;
 }
 
 /* ---------------------------------------------------------------------- */
 
-int PairSMATBSingle::pack_forward_comm(int n, int *list, double *buf,
-                                       int pbc_flag, int *pbc)
+int PairSMATBSingle::pack_forward_comm(int n, int *list, double *buf, int pbc_flag, int *pbc)
 {
-  int i,j,m;
+  int i, j, m;
 
   m = 0;
   for (i = 0; i < n; ++i) {
@@ -413,33 +396,32 @@ int PairSMATBSingle::pack_forward_comm(int n, int *list, double *buf,
 
 /* ---------------------------------------------------------------------- */
 
-void PairSMATBSingle::unpack_forward_comm(int n, int first, double *buf) {
-  int i,m,last;
-  
+void PairSMATBSingle::unpack_forward_comm(int n, int first, double *buf)
+{
+  int i, m, last;
+
   m = 0;
   last = first + n;
-  for (i = first; i < last; ++i) {
-    on_eb[i] = buf[m++];
-  }
+  for (i = first; i < last; ++i) { on_eb[i] = buf[m++]; }
 }
 
 /* ---------------------------------------------------------------------- */
 
-int PairSMATBSingle::pack_reverse_comm(int n, int first, double *buf) {
-  int i,m,last;
-  
+int PairSMATBSingle::pack_reverse_comm(int n, int first, double *buf)
+{
+  int i, m, last;
+
   m = 0;
   last = first + n;
-  for (i = first; i < last; ++i) {
-    buf[m++] = on_eb[i];
-  }
+  for (i = first; i < last; ++i) { buf[m++] = on_eb[i]; }
   return m;
 }
 
 /* ---------------------------------------------------------------------- */
 
-void PairSMATBSingle::unpack_reverse_comm(int n, int *list, double *buf) {
-  int i,j,m;
+void PairSMATBSingle::unpack_reverse_comm(int n, int *list, double *buf)
+{
+  int i, j, m;
   m = 0;
   for (i = 0; i < n; i++) {
     j = list[i];
@@ -448,93 +430,94 @@ void PairSMATBSingle::unpack_reverse_comm(int n, int *list, double *buf) {
 }
 
 //write binary data of this simulation:
-void PairSMATBSingle::write_restart_settings(FILE *fp) {
-  fwrite(&offset_flag,sizeof(int),1,fp);
-  fwrite(&mix_flag,   sizeof(int),1,fp);
-  fwrite(&tail_flag,  sizeof(int),1,fp);
+void PairSMATBSingle::write_restart_settings(FILE *fp)
+{
+  fwrite(&offset_flag, sizeof(int), 1, fp);
+  fwrite(&mix_flag, sizeof(int), 1, fp);
+  fwrite(&tail_flag, sizeof(int), 1, fp);
 }
 
-void PairSMATBSingle::read_restart_settings(FILE *fp) {
+void PairSMATBSingle::read_restart_settings(FILE *fp)
+{
   int me = comm->me;
   size_t result;
   if (me == 0) {
-    result = fread(&offset_flag,sizeof(int),1,fp);
-    result = fread(&mix_flag,   sizeof(int),1,fp);
-    result = fread(&tail_flag,  sizeof(int),1,fp);
+    result = fread(&offset_flag, sizeof(int), 1, fp);
+    result = fread(&mix_flag, sizeof(int), 1, fp);
+    result = fread(&tail_flag, sizeof(int), 1, fp);
   }
-  MPI_Bcast(&offset_flag,  1,MPI_INT,0,world);
-  MPI_Bcast(&mix_flag,     1,MPI_INT,0,world);
-  MPI_Bcast(&tail_flag,    1,MPI_INT,0,world);
+  MPI_Bcast(&offset_flag, 1, MPI_INT, 0, world);
+  MPI_Bcast(&mix_flag, 1, MPI_INT, 0, world);
+  MPI_Bcast(&tail_flag, 1, MPI_INT, 0, world);
 }
 
-void PairSMATBSingle::write_restart(FILE *fp) {
+void PairSMATBSingle::write_restart(FILE *fp)
+{
   write_restart_settings(fp);
   //"I J r0 p q A QSI CO_start CO_end"
-  int i,j;
+  int i, j;
   for (i = 1; i <= atom->ntypes; i++)
     for (j = i; j <= atom->ntypes; j++) {
-      fwrite(&setflag[i][j],        sizeof(int),   1,fp);
+      fwrite(&setflag[i][j], sizeof(int), 1, fp);
       if (setflag[i][j]) {
-        fwrite(&r0,         sizeof(double),1,fp);
-        fwrite(&p,          sizeof(double),1,fp);
-        fwrite(&q,          sizeof(double),1,fp);
-        fwrite(&A,          sizeof(double),1,fp);
-        fwrite(&QSI,        sizeof(double),1,fp);
-        fwrite(&cutOffStart,sizeof(double),1,fp);
-        fwrite(&cutOffEnd,  sizeof(double),1,fp);
+        fwrite(&r0, sizeof(double), 1, fp);
+        fwrite(&p, sizeof(double), 1, fp);
+        fwrite(&q, sizeof(double), 1, fp);
+        fwrite(&A, sizeof(double), 1, fp);
+        fwrite(&QSI, sizeof(double), 1, fp);
+        fwrite(&cutOffStart, sizeof(double), 1, fp);
+        fwrite(&cutOffEnd, sizeof(double), 1, fp);
       }
     }
   //maybe we need to save also the values of the various polynomials
 }
 
-void PairSMATBSingle::read_restart(FILE *fp) {
+void PairSMATBSingle::read_restart(FILE *fp)
+{
   read_restart_settings(fp);
   allocate();
   size_t result;
 
-  int i,j;
+  int i, j;
   int me = comm->me;
   for (i = 1; i <= atom->ntypes; i++)
     for (j = i; j <= atom->ntypes; j++) {
-      if (me == 0) {
-        result = fread(&setflag[i][j],sizeof(int),1,fp);
-      }
-      MPI_Bcast(&setflag[i][j],1,MPI_INT,0,world);
+      if (me == 0) { result = fread(&setflag[i][j], sizeof(int), 1, fp); }
+      MPI_Bcast(&setflag[i][j], 1, MPI_INT, 0, world);
       if (setflag[i][j]) {
         if (me == 0) {
-          result = fread(&r0,         sizeof(double),1,fp);
-          result = fread(&p,          sizeof(double),1,fp);
-          result = fread(&q,          sizeof(double),1,fp);
-          result = fread(&A,          sizeof(double),1,fp);
-          result = fread(&QSI,        sizeof(double),1,fp);
-          result = fread(&cutOffStart,sizeof(double),1,fp);
-          result = fread(&cutOffEnd,  sizeof(double),1,fp);
+          result = fread(&r0, sizeof(double), 1, fp);
+          result = fread(&p, sizeof(double), 1, fp);
+          result = fread(&q, sizeof(double), 1, fp);
+          result = fread(&A, sizeof(double), 1, fp);
+          result = fread(&QSI, sizeof(double), 1, fp);
+          result = fread(&cutOffStart, sizeof(double), 1, fp);
+          result = fread(&cutOffEnd, sizeof(double), 1, fp);
         }
-        MPI_Bcast(&r0,         1,MPI_DOUBLE,0,world);
-        MPI_Bcast(&p,          1,MPI_DOUBLE,0,world);
-        MPI_Bcast(&q,          1,MPI_DOUBLE,0,world);
-        MPI_Bcast(&A,          1,MPI_DOUBLE,0,world);
-        MPI_Bcast(&QSI,        1,MPI_DOUBLE,0,world);
-        MPI_Bcast(&cutOffStart,1,MPI_DOUBLE,0,world);
-        MPI_Bcast(&cutOffEnd,  1,MPI_DOUBLE,0,world);
+        MPI_Bcast(&r0, 1, MPI_DOUBLE, 0, world);
+        MPI_Bcast(&p, 1, MPI_DOUBLE, 0, world);
+        MPI_Bcast(&q, 1, MPI_DOUBLE, 0, world);
+        MPI_Bcast(&A, 1, MPI_DOUBLE, 0, world);
+        MPI_Bcast(&QSI, 1, MPI_DOUBLE, 0, world);
+        MPI_Bcast(&cutOffStart, 1, MPI_DOUBLE, 0, world);
+        MPI_Bcast(&cutOffEnd, 1, MPI_DOUBLE, 0, world);
       }
     }
 }
 
-void PairSMATBSingle::write_data(FILE *fp) {
+void PairSMATBSingle::write_data(FILE *fp)
+{
   //smatb needs I J r0 p q A QSI CO_start CO_end
   for (int i = 1; i <= atom->ntypes; i++) {
-    fprintf(fp,"%d %g %g %g %g %g %g %g\n",
-            i, r0, p, q, A, QSI, cutOffStart, cutOffEnd);
+    fprintf(fp, "%d %g %g %g %g %g %g %g\n", i, r0, p, q, A, QSI, cutOffStart, cutOffEnd);
   }
 }
 
-void PairSMATBSingle::write_data_all(FILE *fp) {
+void PairSMATBSingle::write_data_all(FILE *fp)
+{
   for (int i = 1; i <= atom->ntypes; i++) {
     for (int j = i; j <= atom->ntypes; j++) {
-      fprintf(fp,"%d %d %g %g %g %g %g %g %g\n",
-              i, j, r0, p, q, A, QSI, cutOffStart, cutOffEnd);
+      fprintf(fp, "%d %d %g %g %g %g %g %g %g\n", i, j, r0, p, q, A, QSI, cutOffStart, cutOffEnd);
     }
   }
 }
-
diff --git a/src/pair_smatb_single.h b/src/pair_smatb_single.h
index 774578f395..8eb4b2fdf0 100644
--- a/src/pair_smatb_single.h
+++ b/src/pair_smatb_single.h
@@ -10,7 +10,7 @@
 
    See the README file in the top-level LAMMPS directory.
 
-   This style is written by Daniele Rapetti (iximiel@gmail.com)
+   This pair style is written by Daniele Rapetti (iximiel@gmail.com)
    ------------------------------------------------------------------------- */
 
 #ifdef PAIR_CLASS
@@ -26,53 +26,51 @@ PairStyle(smatb/single,PairSMATBSingle)
 
 namespace LAMMPS_NS {
 
-
-  class PairSMATBSingle : public Pair {
-  public:
-    // public variables so USER-ATC package can access them
-
-    PairSMATBSingle(class LAMMPS *);
-    virtual ~PairSMATBSingle();
-    void compute(int, int);//workhorse routine that computes pairwise interactions
-    /*
+class PairSMATBSingle : public Pair {
+ public:
+  PairSMATBSingle(class LAMMPS *);
+  virtual ~PairSMATBSingle();
+  void compute(int, int);    //workhorse routine that computes pairwise interactions
+  /*
     void compute_inner();
     void compute_middle();
-    void compute_outer(int, int);*/
-    void settings(int, char **);//reads the input script line with arguments you define
-    void coeff(int, char **);//set coefficients for one i,j type pair
-    //void init_style();//initialization specific to this pair style
-    double init_one(int, int);//perform initialization for one i,j type pair
-    //double single(int, int, int, int, double, double, double, double &);//force and energy of a single pairwise interaction between 2 atoms
+    void compute_outer(int, int);
+    */
+  void settings(int, char **);    //reads the input script line with arguments you define
+  void coeff(int, char **);       //set coefficients for one i,j type pair
+  //void init_style();//initialization specific to this pair style
+  double init_one(int, int);    //perform initialization for one i,j type pair
 
-    virtual void write_restart(FILE *);
-    virtual void read_restart(FILE *);
-    virtual void write_restart_settings(FILE *);
-    virtual void read_restart_settings(FILE *);
-    virtual void write_data(FILE *);
-    virtual void write_data_all(FILE *);
+  virtual void write_restart(FILE *);
+  virtual void read_restart(FILE *);
+  virtual void write_restart_settings(FILE *);
+  virtual void read_restart_settings(FILE *);
+  virtual void write_data(FILE *);
+  virtual void write_data_all(FILE *);
 
-    virtual int    pack_forward_comm(int, int *, double *, int, int *);
-    virtual void unpack_forward_comm(int, int, double *);
-    virtual int    pack_reverse_comm(int , int , double *);
-    virtual void unpack_reverse_comm(int , int *, double *);
-  protected: 
-    virtual void allocate();
-    int nmax;            // allocated size of per-atom arrays
-    double *on_eb;       //allocated to store up caclulation values
-    double r0;           // interaction radius, user-given
-    double p, A, q, QSI; // parameters user-given
-    double cutOffStart, cutOffEnd;//cut offs, user given
-    double cutOffEnd2;  //squared cut off end, calculated
-    double a3, a4, a5;  //polynomial for cutoff linking to zero:   Ae^p substitution
-    double x3, x4, x5;  //polynomial for cutoff linking to zero: QSIe^q substitution
-  };
+  virtual int pack_forward_comm(int, int *, double *, int, int *);
+  virtual void unpack_forward_comm(int, int, double *);
+  virtual int pack_reverse_comm(int, int, double *);
+  virtual void unpack_reverse_comm(int, int *, double *);
 
-}
+ protected:
+  virtual void allocate();
+  int nmax;                         // allocated size of per-atom arrays
+  double *on_eb;                    //allocated to store up caclulation values
+  double r0;                        // interaction radius, user-given
+  double p, A, q, QSI;              // parameters user-given
+  double cutOffStart, cutOffEnd;    //cut offs, user given
+  double cutOffEnd2;                //squared cut off end, calculated
+  double a3, a4, a5;                //polynomial for cutoff linking to zero:   Ae^p substitution
+  double x3, x4, x5;                //polynomial for cutoff linking to zero: QSIe^q substitution
+};
+
+}    // namespace LAMMPS_NS
 
 #endif
 #endif
 
-/* ERROR/WARNING messages:
+    /* ERROR/WARNING messages:
 
    E: Illegal ... command
 
@@ -83,10 +81,4 @@ namespace LAMMPS_NS {
    E: Incorrect args for pair coefficients
 
    Self-explanatory.  Check the input script or data file.
-
-   E: Cannot open EAM potential file %s
-
-   The specified EAM potential file cannot be opened.  Check that the
-   path and name are correct.
-
 */