diff --git a/src/CORESHELL/compute_temp_cs.cpp b/src/CORESHELL/compute_temp_cs.cpp
new file mode 100644
index 0000000000..f786b6533b
--- /dev/null
+++ b/src/CORESHELL/compute_temp_cs.cpp
@@ -0,0 +1,503 @@
+/* ----------------------------------------------------------------------
+   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: Hendrik Heenen (Technical University of Munich)
+                        (hendrik.heenen at mytum.com)
+------------------------------------------------------------------------- */
+
+#include "mpi.h"
+#include "stdlib.h"
+#include "string.h"
+#include "math.h"
+#include "compute_temp_cs.h"
+#include "atom.h"
+#include "atom_vec.h"
+#include "domain.h"
+#include "update.h"
+#include "force.h"
+#include "group.h"
+#include "modify.h"
+#include "fix.h"
+#include "fix_store.h"
+#include "comm.h"
+#include "memory.h"
+#include "error.h"
+
+using namespace LAMMPS_NS;
+
+/* ---------------------------------------------------------------------- */
+
+ComputeTempCS::ComputeTempCS(LAMMPS *lmp, int narg, char **arg) :
+  Compute(lmp, narg, arg)
+{
+  if (narg != 5) error->all(FLERR,"Illegal compute temp/cs command");
+
+  if (atom->avec->bonds_allow == 0)
+    error->all(FLERR,"Compute temp/cs used when bonds are not allowed");
+
+  scalar_flag = vector_flag = 1;
+  size_vector = 6;
+  extscalar = 0;
+  extvector = 1;
+  tempflag = 1;
+  tempbias = 1;
+  extarray = 0;
+
+  int *mask = atom->mask;
+  int nlocal = atom->nlocal;
+
+  // find and define groupbits for core and shell groups
+
+  cgroup = group->find(arg[3]);
+  if (cgroup == -1) 
+    error->all(FLERR,"Could not find specified group ID for core particles");
+  groupbit_c = group->bitmask[cgroup];
+
+  sgroup = group->find(arg[4]);
+  if (sgroup == -1) 
+    error->all(FLERR,"Could not find specified group ID for shell particles");
+  groupbit_s = group->bitmask[sgroup];
+
+  // create a new fix STORE style
+  // id = compute-ID + COMPUTE_STORE, fix group = compute group
+
+  int n = strlen(id) + strlen("_COMPUTE_STORE") + 1;
+  id_fix = new char[n];
+  strcpy(id_fix,id);
+  strcat(id_fix,"_COMPUTE_STORE");
+
+  char **newarg = new char*[5];
+  newarg[0] = id_fix;
+  newarg[1] = group->names[igroup];
+  newarg[2] = (char *) "STORE";
+  newarg[3] = (char *) "0";
+  newarg[4] = (char *) "1";
+  modify->add_fix(5,newarg);
+  fix = (FixStore *) modify->fix[modify->nfix-1];
+  delete [] newarg;
+
+  // set fix store values = 0 for now
+  // fill them in via setup() once Comm::borders() has been called
+  // skip if resetting from restart file
+
+  if (fix->restart_reset) {
+    fix->restart_reset = 0;
+    firstflag = 0;
+  } else {
+    double *partner = fix->vstore;
+    int nlocal = atom->nlocal;
+    for (int i = 0; i < nlocal; i++) partner[i] = ubuf(0).d;
+    firstflag = 1;
+  }
+
+  // allocate memory
+
+  vector = new double[6];
+  maxatom = 0;
+  vint = NULL;
+
+  // set comm size needed by this Compute
+
+  comm_reverse = 1;
+}
+
+/* ---------------------------------------------------------------------- */
+
+ComputeTempCS::~ComputeTempCS()
+{
+  // check nfix in case all fixes have already been deleted
+
+  if (modify->nfix) modify->delete_fix(id_fix);
+
+  delete [] id_fix;
+  delete [] vector;
+  memory->destroy(vint);
+}
+
+/* ---------------------------------------------------------------------- */
+
+void ComputeTempCS::init()
+{
+  if (comm->ghost_velocity == 0)
+    error->all(FLERR,"Compute temp/cs requires ghost atoms store velocity");
+}
+
+/* ---------------------------------------------------------------------- */
+
+void ComputeTempCS::setup()
+{
+  if (firstflag) {
+    firstflag = 0;
+
+    // insure # of core atoms = # of shell atoms
+
+    int ncores = group->count(cgroup);
+    nshells = group->count(sgroup);
+    if (ncores != nshells)
+      error->all(FLERR,"Number of core atoms != number of shell atoms");
+
+    // for each C/S pair:
+    // set partner IDs of both atoms if this atom stores bond between them
+    // will set partner IDs for ghost atoms if needed by another proc
+    // nall loop insures all ghost atom partner IDs are set before reverse comm
+
+    int *num_bond = atom->num_bond;
+    tagint **bond_atom = atom->bond_atom;
+    tagint *tag = atom->tag;
+    int *mask = atom->mask;
+    int nlocal = atom->nlocal;
+
+    double *partner = fix->vstore;
+    tagint partnerID;
+
+    int nall = nlocal + atom->nghost;
+    for (int i = nlocal; i < nall; i++) partner[i] = ubuf(0).d;
+
+    int i,j,m,match;
+    for (i = 0; i < nlocal; i++) {
+      if (mask[i] & groupbit_c || mask[i] & groupbit_s) {
+        for (m = 0; m < num_bond[i]; m++) {
+          partnerID = bond_atom[i][m];
+          j = atom->map(partnerID);
+          if (j == -1) error->one(FLERR,"Core/shell partner atom not found");
+          match = 0;
+          if (mask[i] & groupbit_c && mask[j] & groupbit_s) match = 1;
+          if (mask[i] & groupbit_s && mask[j] & groupbit_c) match = 1;
+          if (match) {
+            partner[i] = ubuf(partnerID).d;
+            partner[j] = ubuf(tag[i]).d;
+          }
+        }
+      }
+    }
+
+    // reverse comm to acquire unknown partner IDs from ghost atoms
+    // only needed if newton_bond = on
+
+    if (force->newton_bond) comm->reverse_comm_compute(this);
+
+    // check that all C/S partners were found
+
+    int flag = 0;
+    for (i = 0; i < nlocal; i++) {
+      if (mask[i] & groupbit_c || mask[i] & groupbit_s) {
+        partnerID = (tagint) ubuf(partner[i]).i;
+        if (partnerID == 0) flag = 1;
+      }
+    }
+
+    int flagall;
+    MPI_Allreduce(&flag,&flagall,1,MPI_INT,MPI_SUM,world);
+    if (flagall) error->all(FLERR,"Core/shell partners were not all found");
+  }
+
+  // calculate DOF for temperature
+
+  fix_dof = -1;
+  dof_compute();
+}
+
+/* ---------------------------------------------------------------------- */
+
+void ComputeTempCS::dof_compute()
+{
+  if (fix_dof) adjust_dof_fix();
+  int nper = domain->dimension;
+  double natoms = group->count(igroup);
+  dof = nper * natoms;
+  dof -= nper * nshells;
+  dof -= extra_dof + fix_dof;
+  if (dof > 0) tfactor = force->mvv2e / (dof * force->boltz);
+  else tfactor = 0.0;
+}
+
+/* ---------------------------------------------------------------------- */
+
+double ComputeTempCS::compute_scalar()
+{
+  int i;
+  double vthermal[3];
+
+  invoked_scalar = update->ntimestep;
+
+  vcm_pairs();
+
+  // calculate thermal scalar in respect to atom velocities as center-of-mass 
+  // velocities of its according core/shell pairs
+  
+  double **v = atom->v;
+  int *mask = atom->mask;
+  int *type = atom->type;
+  double *mass = atom->mass;
+  double *rmass = atom->rmass;
+  int nlocal = atom->nlocal;
+
+  double t = 0.0;
+
+  for (int i = 0; i < nlocal; i++){
+    if (mask[i] & groupbit) {
+      vthermal[0] = v[i][0] - vint[i][0];
+      vthermal[1] = v[i][1] - vint[i][1];
+      vthermal[2] = v[i][2] - vint[i][2];
+      if (rmass)
+        t += (vthermal[0]*vthermal[0] + vthermal[1]*vthermal[1] +
+              vthermal[2]*vthermal[2]) * rmass[i];
+      else
+        t += (vthermal[0]*vthermal[0] + vthermal[1]*vthermal[1] +
+              vthermal[2]*vthermal[2]) * mass[type[i]];
+    }
+  }
+
+  MPI_Allreduce(&t,&scalar,1,MPI_DOUBLE,MPI_SUM,world);
+  if (dynamic) dof_compute();
+  scalar *= tfactor;
+  return scalar;
+}
+    
+/* ---------------------------------------------------------------------- */
+
+void ComputeTempCS::compute_vector()
+{
+  double vthermal[3];
+
+  invoked_vector = update->ntimestep;
+
+  vcm_pairs();
+
+  // calculate thermal vector in respect to atom velocities as center-of-mass 
+  // velocities of its according C/S pairs
+
+  double **v = atom->v;
+  int *mask = atom->mask;
+  tagint *molecule = atom->molecule;
+  int *type = atom->type;
+  double *mass = atom->mass;
+  double *rmass = atom->rmass;
+  int nlocal = atom->nlocal;
+  double massone;
+
+  double t[6];
+  for (int i = 0; i < 6; i++) t[i] = 0.0;
+
+  for (int i = 0; i < nlocal; i++){
+    if (mask[i] & groupbit) {
+      vthermal[0] = v[i][0] - vint[i][0];
+      vthermal[1] = v[i][1] - vint[i][1];
+      vthermal[2] = v[i][2] - vint[i][2];
+      if (rmass) massone = rmass[i];
+      else massone = mass[type[i]];
+      t[0] += massone * vthermal[0]*vthermal[0];
+      t[1] += massone * vthermal[1]*vthermal[1];
+      t[2] += massone * vthermal[2]*vthermal[2];
+      t[3] += massone * vthermal[0]*vthermal[1];
+      t[4] += massone * vthermal[0]*vthermal[2];
+      t[5] += massone * vthermal[1]*vthermal[2];
+    }
+  }
+
+  MPI_Allreduce(t,vector,6,MPI_DOUBLE,MPI_SUM,world);
+  for (int i = 0; i < 6; i++) vector[i] *= force->mvv2e;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void ComputeTempCS::vcm_pairs()
+{
+  int i,j;
+  double massone,masstwo;
+  double vcm[3];
+
+  // reallocate vint if necessary
+
+  int nlocal = atom->nlocal;
+
+  if (nlocal > maxatom) {
+    memory->destroy(vint);
+    maxatom = atom->nmax;
+    memory->create(vint,maxatom,3,"temp/cs:vint");
+  }
+  
+  // vcm = COM velocity of each CS pair
+  // vint = internal velocity of each C/S atom, used as bias
+
+  double **v = atom->v;
+  int *mask = atom->mask;
+  int *type = atom->type;
+  double *mass = atom->mass;
+  double *rmass = atom->rmass;
+
+  double *partner = fix->vstore;
+  tagint partnerID;
+
+  for (int i = 0; i < nlocal; i++) {
+    if ((mask[i] & groupbit) && 
+        (mask[i] & groupbit_c || mask[i] & groupbit_s)) {
+      if (rmass) massone = rmass[i];
+      else massone = mass[type[i]];
+      vcm[0] = v[i][0]*massone;
+      vcm[1] = v[i][1]*massone;
+      vcm[2] = v[i][2]*massone;
+
+      partnerID = (tagint) ubuf(partner[i]).i;
+      j = atom->map(partnerID);
+      if (j == -1) error->one(FLERR,"Core/shell partner atom not found");
+
+      if (rmass) masstwo = rmass[j];
+      else masstwo = mass[type[j]];
+      vcm[0] += v[j][0]*masstwo;
+      vcm[1] += v[j][1]*masstwo;
+      vcm[2] += v[j][2]*masstwo;
+      vcm[0] /= (massone + masstwo);
+      vcm[1] /= (massone + masstwo);
+      vcm[2] /= (massone + masstwo);
+
+      vint[i][0] = v[i][0] - vcm[0];
+      vint[i][1] = v[i][1] - vcm[1];
+      vint[i][2] = v[i][2] - vcm[2];
+    } else vint[i][0] = vint[i][1] = vint[i][2] = 0.0;
+  }
+}
+
+/* ----------------------------------------------------------------------
+   remove velocity bias from atom I to leave thermal velocity
+   thermal velocity in this case is COM velocity of C/S pair
+------------------------------------------------------------------------- */
+
+void ComputeTempCS::remove_bias(int i, double *v)
+{
+  v[0] -= vint[i][0];
+  v[1] -= vint[i][1];
+  v[2] -= vint[i][2];
+}
+
+/* ----------------------------------------------------------------------
+   remove velocity bias from all atoms to leave thermal velocity
+   thermal velocity in this case is COM velocity of C/S pair
+------------------------------------------------------------------------- */
+
+void ComputeTempCS::remove_bias_all()
+{
+  double **v = atom->v;
+  int *mask = atom->mask;
+  int nlocal = atom->nlocal;
+
+  for (int i = 0; i < nlocal; i++)
+    if (mask[i] & groupbit) {
+      v[i][0] -= vint[i][0];
+      v[i][1] -= vint[i][1];
+      v[i][2] -= vint[i][2];
+    }
+}
+
+/* ----------------------------------------------------------------------
+   reset thermal velocity of all atoms to be consistent with bias
+   called from velocity command after it creates thermal velocities
+   this resets each atom's velocity to COM velocity of C/S pair
+------------------------------------------------------------------------- */
+
+void ComputeTempCS::reapply_bias_all()
+{
+  double **v = atom->v;
+  int *mask = atom->mask;
+  int nlocal = atom->nlocal;
+
+  // recalculate current COM velocities
+
+  vcm_pairs();
+
+  // zero vint after using ti so that Velocity call to restore_bias_all()
+  // will not further alter the velocities within a C/S pair
+
+  for (int i = 0; i < nlocal; i++)
+    if (mask[i] & groupbit) {
+      v[i][0] -= vint[i][0];
+      v[i][1] -= vint[i][1];
+      v[i][2] -= vint[i][2];
+      vint[i][0] = 0.0;
+      vint[i][1] = 0.0;
+      vint[i][2] = 0.0;
+    }
+}
+
+/* ----------------------------------------------------------------------
+   add back in velocity bias to atom I removed by remove_bias()
+   assume remove_bias() was previously called
+------------------------------------------------------------------------- */
+
+void ComputeTempCS::restore_bias(int i, double *v)
+{
+  v[0] += vint[i][0];
+  v[1] += vint[i][1];
+  v[2] += vint[i][2];
+}
+
+/* ----------------------------------------------------------------------
+   add back in velocity bias to all atoms removed by remove_bias_all()
+   assume remove_bias_all() was previously called
+------------------------------------------------------------------------- */
+
+void ComputeTempCS::restore_bias_all()
+{
+  double **v = atom->v;
+  int *mask = atom->mask;
+  int nlocal = atom->nlocal;
+
+  for (int i = 0; i < nlocal; i++)
+    if (mask[i] & groupbit) {
+      v[i][0] += vint[i][0];
+      v[i][1] += vint[i][1];
+      v[i][2] += vint[i][2];
+    }
+}
+
+/* ---------------------------------------------------------------------- */
+
+int ComputeTempCS::pack_reverse_comm(int n, int first, double *buf)
+{
+  int i,m,last;
+
+  double *partner = fix->vstore;
+
+  m = 0;
+  last = first + n;
+  for (i = first; i < last; i++) buf[m++] = partner[i];
+  return m;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void ComputeTempCS::unpack_reverse_comm(int n, int *list, double *buf)
+{
+  int i,j,m;
+
+  double *partner = fix->vstore;
+  tagint partnerID;
+
+  m = 0;
+  for (i = 0; i < n; i++) {
+    j = list[i];
+    partnerID = (tagint) ubuf(buf[m++]).i;
+    if (partnerID) partner[j] = ubuf(partnerID).d;
+  }
+}
+
+/* ----------------------------------------------------------------------
+   memory usage of local data
+------------------------------------------------------------------------- */
+
+double ComputeTempCS::memory_usage()
+{
+  double bytes = (bigint) maxatom * 3 * sizeof(double);
+  return bytes;
+}
diff --git a/src/CORESHELL/compute_temp_cs.h b/src/CORESHELL/compute_temp_cs.h
new file mode 100644
index 0000000000..7b25e9878d
--- /dev/null
+++ b/src/CORESHELL/compute_temp_cs.h
@@ -0,0 +1,120 @@
+/* ----------------------------------------------------------------------
+   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 COMPUTE_CLASS
+
+ComputeStyle(temp/cs,ComputeTempCS)
+
+#else
+
+#ifndef LMP_COMPUTE_TEMP_CS_H
+#define LMP_COMPUTE_TEMP_CS_H
+
+#include "compute.h"
+
+namespace LAMMPS_NS {
+
+class ComputeTempCS : public Compute {
+ public:
+  ComputeTempCS(class LAMMPS *, int, char **);
+  ~ComputeTempCS();
+  void init();
+  void setup();
+  double compute_scalar();
+  void compute_vector();
+  double memory_usage();
+
+  void remove_bias(int, double *);
+  void remove_bias_all();
+  void reapply_bias_all();
+  void restore_bias(int, double *);
+  void restore_bias_all();
+
+  int pack_reverse_comm(int, int, double *);
+  void unpack_reverse_comm(int, int *, double *);
+
+ private:
+  int groupbit_c,groupbit_s;
+  int nshells;
+  int firstflag;
+  int maxatom;
+  int cgroup,sgroup;
+
+  int fix_dof;
+  double tfactor;
+  double **vint;
+
+  char *id_fix;
+  class FixStore *fix;
+
+  void dof_compute();
+  void vcm_pairs();
+};
+
+}
+
+#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: Option mol of compute temp/vcm requires molecular atom style
+
+Self-explanatory.
+
+E: Option prop of compute temp/vcm requires one set of parameters 
+added by the property/atom fix
+
+Self-explanatory.
+
+E: Fix property/atom vector must contain only intergers to assign 
+sub-ID property
+
+Self-explanatory.
+
+E: Specified sub-ID property does not exist or has not been created 
+by the property/atom fix
+
+Self-explanatory. Usually this means that the specified fix 
+property/atom ID does not match the ID stated in the compute temp/vcm.
+
+E: Molecule count changed in compute com/temp/molecule
+
+Number of molecules must remain constant over time.
+
+E: Sub-ID count changed in compute vcm/temp
+
+Number of Sub-ID groups must remain constant over time.
+
+W: Atom with sub-ID = 0 included in compute group
+
+Self-explanatory. A sub-ID with value 0 will be counted as a normal sub-ID 
+and not left out of by the compute treatment. Therefore a sub-ID of 0 is to 
+be avoided.
+
+E: Too many sub-ID groups for compute
+
+Self-explanatory.
+
+W: More than 2 atoms specified with the same sub-ID, in the case 
+of a core-shell model simulation only core and shell should share the same ID
+
+Self-explanatory.
+
+*/
diff --git a/src/CORESHELL/pair_born_coul_long_cs.cpp b/src/CORESHELL/pair_born_coul_long_cs.cpp
new file mode 100644
index 0000000000..6e87f47e72
--- /dev/null
+++ b/src/CORESHELL/pair_born_coul_long_cs.cpp
@@ -0,0 +1,199 @@
+/* ----------------------------------------------------------------------
+   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: Hendrik Heenen (hendrik.heenen@mytum.de)
+------------------------------------------------------------------------- */
+
+#include "math.h"
+#include "stdio.h"
+#include "stdlib.h"
+#include "string.h"
+#include "pair_born_coul_long_cs.h"
+#include "atom.h"
+#include "comm.h"
+#include "force.h"
+#include "kspace.h"
+#include "neighbor.h"
+#include "neigh_list.h"
+#include "math_const.h"
+#include "memory.h"
+#include "error.h"
+
+using namespace LAMMPS_NS;
+using namespace MathConst;
+
+#define EWALD_F   1.12837917
+#define EWALD_P   0.3275911
+#define A1        0.254829592
+#define A2       -0.284496736
+#define A3        1.421413741
+#define A4       -1.453152027
+#define A5        1.061405429
+
+#define EPSILON 1.0e-20
+#define EPS_EWALD 1.0e-6
+#define EPS_EWALD_SQR 1.0e-12
+
+/* ---------------------------------------------------------------------- */
+
+PairBornCoulLongCS::PairBornCoulLongCS(LAMMPS *lmp) : PairBornCoulLong(lmp)
+{
+  ewaldflag = pppmflag = 1;
+  ftable = NULL;
+  writedata = 1;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void PairBornCoulLongCS::compute(int eflag, int vflag)
+{
+  int i,j,ii,jj,inum,jnum,itable,itype,jtype;
+  double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,evdwl,ecoul,fpair;
+  double fraction,table;
+  double rsq,r2inv,r6inv,forcecoul,forceborn,factor_coul,factor_lj;
+  double grij,expm2,prefactor,t,erfc;
+  double r,rexp;
+  int *ilist,*jlist,*numneigh,**firstneigh;
+
+  evdwl = ecoul = 0.0;
+  if (eflag || vflag) ev_setup(eflag,vflag);
+  else evflag = vflag_fdotr = 0;
+
+  double **x = atom->x;
+  double **f = atom->f;
+  double *q = atom->q;
+  int *type = atom->type;
+  int nlocal = atom->nlocal;
+  double *special_coul = force->special_coul;
+  double *special_lj = force->special_lj;
+  int newton_pair = force->newton_pair;
+  double qqrd2e = force->qqrd2e;
+
+  inum = list->inum;
+  ilist = list->ilist;
+  numneigh = list->numneigh;
+  firstneigh = list->firstneigh;
+
+  // loop over neighbors of my atoms
+
+  for (ii = 0; ii < inum; ii++) {
+    i = ilist[ii];
+    qtmp = q[i];
+    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];
+      factor_lj = special_lj[sbmask(j)];
+      factor_coul = special_coul[sbmask(j)];
+      j &= NEIGHMASK;
+
+      delx = xtmp - x[j][0];
+      dely = ytmp - x[j][1];
+      delz = ztmp - x[j][2];
+      rsq = delx*delx + dely*dely + delz*delz;
+      jtype = type[j];
+
+      if (rsq < cutsq[itype][jtype]) {
+        rsq += EPSILON; // Add Epsilon for case: r = 0; Interaction must be removed by special bond;
+        r2inv = 1.0/rsq;
+        if (rsq < cut_coulsq) {
+          if (!ncoultablebits || rsq <= tabinnersq) {
+            r = sqrt(rsq);
+            prefactor = qqrd2e * qtmp*q[j];
+            if (factor_coul < 1.0) {
+              // When bonded parts are being calculated a minimal distance (EPS_EWALD) 
+              // has to be added to the prefactor and erfc in order to make the 
+              // used approximation functions valid
+              grij = g_ewald * (r+EPS_EWALD);
+              expm2 = exp(-grij*grij);
+              t = 1.0 / (1.0 + EWALD_P*grij);
+              erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2;
+              prefactor /= (r+EPS_EWALD);
+              forcecoul = prefactor * (erfc + EWALD_F*grij*expm2 - (1.0-factor_coul));
+              // Additionally r2inv needs to be accordingly modified since the later
+              // scaling of the overall force shall be consistent
+              r2inv = 1.0/(rsq + EPS_EWALD_SQR);
+            } else {
+              grij = g_ewald * r;
+              expm2 = exp(-grij*grij);
+              t = 1.0 / (1.0 + EWALD_P*grij);
+              erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2;
+              prefactor /= r;
+              forcecoul = prefactor * (erfc + EWALD_F*grij*expm2);
+            }
+          } else {
+            union_int_float_t rsq_lookup;
+            rsq_lookup.f = rsq;
+            itable = rsq_lookup.i & ncoulmask;
+            itable >>= ncoulshiftbits;
+            fraction = (rsq_lookup.f - rtable[itable]) * drtable[itable];
+            table = ftable[itable] + fraction*dftable[itable];
+            forcecoul = qtmp*q[j] * table;
+            if (factor_coul < 1.0) {
+              table = ctable[itable] + fraction*dctable[itable];
+              prefactor = qtmp*q[j] * table;
+              forcecoul -= (1.0-factor_coul)*prefactor;
+            }
+          }
+        } else forcecoul = 0.0;
+
+        if (rsq < cut_ljsq[itype][jtype]) {
+          r = sqrt(rsq);
+          r6inv = r2inv*r2inv*r2inv;
+          rexp = exp((sigma[itype][jtype]-r)*rhoinv[itype][jtype]);
+          forceborn = born1[itype][jtype]*r*rexp - born2[itype][jtype]*r6inv
+            + born3[itype][jtype]*r2inv*r6inv;
+        } else forceborn = 0.0;
+
+        fpair = (forcecoul + factor_lj*forceborn) * r2inv;
+
+        f[i][0] += delx*fpair;
+        f[i][1] += dely*fpair;
+        f[i][2] += delz*fpair;
+        if (newton_pair || j < nlocal) {
+          f[j][0] -= delx*fpair;
+          f[j][1] -= dely*fpair;
+          f[j][2] -= delz*fpair;
+        }
+
+        if (eflag) {
+          if (rsq < cut_coulsq) {
+            if (!ncoultablebits || rsq <= tabinnersq)
+              ecoul = prefactor*erfc;
+            else {
+              table = etable[itable] + fraction*detable[itable];
+              ecoul = qtmp*q[j] * table;
+            }
+            if (factor_coul < 1.0) ecoul -= (1.0-factor_coul)*prefactor;
+          } else ecoul = 0.0;
+          if (rsq < cut_ljsq[itype][jtype]) {
+            evdwl = a[itype][jtype]*rexp - c[itype][jtype]*r6inv
+              + d[itype][jtype]*r6inv*r2inv - offset[itype][jtype];
+            evdwl *= factor_lj;
+          } else evdwl = 0.0;
+        }
+
+        if (evflag) ev_tally(i,j,nlocal,newton_pair,
+                             evdwl,ecoul,fpair,delx,dely,delz);
+      }
+    }
+  }
+
+  if (vflag_fdotr) virial_fdotr_compute();
+}
diff --git a/src/CORESHELL/pair_born_coul_long_cs.h b/src/CORESHELL/pair_born_coul_long_cs.h
new file mode 100644
index 0000000000..1271de38ab
--- /dev/null
+++ b/src/CORESHELL/pair_born_coul_long_cs.h
@@ -0,0 +1,67 @@
+/* ----------------------------------------------------------------------
+   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: Hendrik Heenen (hendrik.heenen@mytum.com)
+------------------------------------------------------------------------- */
+
+#ifdef PAIR_CLASS
+
+PairStyle(born/coul/long/cs,PairBornCoulLongCS)
+
+#else
+
+#ifndef LMP_PAIR_BORN_COUL_LONG_CS_H
+#define LMP_PAIR_BORN_COUL_LONG_CS_H
+
+#include "pair_born_coul_long.h"
+
+namespace LAMMPS_NS {
+
+class PairBornCoulLongCS : public PairBornCoulLong {
+ public:
+  PairBornCoulLongCS(class LAMMPS *);
+  virtual void compute(int, int);
+};
+
+}
+
+#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.
+
+E: All pair coeffs are not set
+
+All pair coefficients must be set in the data file or by the
+pair_coeff command before running a simulation.
+
+E: Pair style born/coul/long requires atom attribute q
+
+An atom style that defines this attribute must be used.
+
+E: Pair style requires a KSpace style
+
+No kspace style is defined.
+
+*/
diff --git a/src/CORESHELL/pair_buck_coul_long_cs.cpp b/src/CORESHELL/pair_buck_coul_long_cs.cpp
new file mode 100644
index 0000000000..7e2820cf42
--- /dev/null
+++ b/src/CORESHELL/pair_buck_coul_long_cs.cpp
@@ -0,0 +1,197 @@
+/* ----------------------------------------------------------------------
+   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: Hendrik Heenen (hendrik.heenen@mytum.de)
+------------------------------------------------------------------------- */
+
+#include "math.h"
+#include "stdio.h"
+#include "stdlib.h"
+#include "string.h"
+#include "pair_buck_coul_long_cs.h"
+#include "atom.h"
+#include "comm.h"
+#include "force.h"
+#include "kspace.h"
+#include "neighbor.h"
+#include "neigh_list.h"
+#include "math_const.h"
+#include "memory.h"
+#include "error.h"
+
+using namespace LAMMPS_NS;
+using namespace MathConst;
+
+#define EWALD_F   1.12837917
+#define EWALD_P   0.3275911
+#define A1        0.254829592
+#define A2       -0.284496736
+#define A3        1.421413741
+#define A4       -1.453152027
+#define A5        1.061405429
+
+#define EPSILON 1.0e-20
+#define EPS_EWALD 1.0e-6
+#define EPS_EWALD_SQR 1.0e-12
+
+/* ---------------------------------------------------------------------- */
+
+PairBuckCoulLongCS::PairBuckCoulLongCS(LAMMPS *lmp) : PairBuckCoulLong(lmp)
+{
+  ewaldflag = pppmflag = 1;
+  writedata = 1;
+  ftable = NULL;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void PairBuckCoulLongCS::compute(int eflag, int vflag)
+{
+  int i,j,ii,jj,inum,jnum,itable,itype,jtype;
+  double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,evdwl,ecoul,fpair;
+  double fraction,table;
+  double rsq,r2inv,r6inv,forcecoul,forcebuck,factor_coul,factor_lj;
+  double grij,expm2,prefactor,t,erfc;
+  double r,rexp;
+  int *ilist,*jlist,*numneigh,**firstneigh;
+
+  evdwl = ecoul = 0.0;
+  if (eflag || vflag) ev_setup(eflag,vflag);
+  else evflag = vflag_fdotr = 0;
+
+  double **x = atom->x;
+  double **f = atom->f;
+  double *q = atom->q;
+  int *type = atom->type;
+  int nlocal = atom->nlocal;
+  double *special_coul = force->special_coul;
+  double *special_lj = force->special_lj;
+  int newton_pair = force->newton_pair;
+  double qqrd2e = force->qqrd2e;
+
+  inum = list->inum;
+  ilist = list->ilist;
+  numneigh = list->numneigh;
+  firstneigh = list->firstneigh;
+
+  // loop over neighbors of my atoms
+
+  for (ii = 0; ii < inum; ii++) {
+    i = ilist[ii];
+    qtmp = q[i];
+    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];
+      factor_lj = special_lj[sbmask(j)];
+      factor_coul = special_coul[sbmask(j)];
+      j &= NEIGHMASK;
+
+      delx = xtmp - x[j][0];
+      dely = ytmp - x[j][1];
+      delz = ztmp - x[j][2];
+      rsq = delx*delx + dely*dely + delz*delz;
+      jtype = type[j];
+
+      if (rsq < cutsq[itype][jtype]) {
+        rsq += EPSILON; // Add Epsilon for case: r = 0; Interaction must be removed by special bond;
+        r2inv = 1.0/rsq;
+        if (rsq < cut_coulsq) {
+          if (!ncoultablebits || rsq <= tabinnersq) {
+            r = sqrt(rsq);
+            prefactor = qqrd2e * qtmp*q[j];
+            if (factor_coul < 1.0) {
+              // When bonded parts are being calculated a minimal distance (EPS_EWALD) 
+              // has to be added to the prefactor and erfc in order to make the 
+              // used approximation functions for the Ewald correction valid
+              grij = g_ewald * (r+EPS_EWALD);
+              expm2 = exp(-grij*grij);
+              t = 1.0 / (1.0 + EWALD_P*grij);
+              erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2;
+              prefactor /= (r+EPS_EWALD);
+              forcecoul = prefactor * (erfc + EWALD_F*grij*expm2 - (1.0-factor_coul));
+              // Additionally r2inv needs to be accordingly modified since the later
+              // scaling of the overall force shall be consistent
+              r2inv = 1.0/(rsq + EPS_EWALD_SQR);
+            } else {
+              grij = g_ewald * r;
+              expm2 = exp(-grij*grij);
+              t = 1.0 / (1.0 + EWALD_P*grij);
+              erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2;
+              prefactor /= r;
+              forcecoul = prefactor * (erfc + EWALD_F*grij*expm2);
+            }
+          } else {
+            union_int_float_t rsq_lookup;
+            rsq_lookup.f = rsq;
+            itable = rsq_lookup.i & ncoulmask;
+            itable >>= ncoulshiftbits;
+            fraction = (rsq_lookup.f - rtable[itable]) * drtable[itable];
+            table = ftable[itable] + fraction*dftable[itable];
+            forcecoul = qtmp*q[j] * table;
+            if (factor_coul < 1.0) {
+              table = ctable[itable] + fraction*dctable[itable];
+              prefactor = qtmp*q[j] * table;
+              forcecoul -= (1.0-factor_coul)*prefactor;
+            }
+          }
+        } else forcecoul = 0.0;
+
+        if (rsq < cut_ljsq[itype][jtype]) {
+          r = sqrt(rsq);
+          r6inv = r2inv*r2inv*r2inv;
+          rexp = exp(-r*rhoinv[itype][jtype]);
+          forcebuck = buck1[itype][jtype]*r*rexp - buck2[itype][jtype]*r6inv;
+        } else forcebuck = 0.0;
+        fpair = (forcecoul + factor_lj*forcebuck) * r2inv;
+
+        f[i][0] += delx*fpair;
+        f[i][1] += dely*fpair;
+        f[i][2] += delz*fpair;
+        if (newton_pair || j < nlocal) {
+          f[j][0] -= delx*fpair;
+          f[j][1] -= dely*fpair;
+          f[j][2] -= delz*fpair;
+        }
+
+        if (eflag) {
+          if (rsq < cut_coulsq) {
+            if (!ncoultablebits || rsq <= tabinnersq)
+              ecoul = prefactor*erfc;
+            else {
+              table = etable[itable] + fraction*detable[itable];
+              ecoul = qtmp*q[j] * table;
+            }
+            if (factor_coul < 1.0) ecoul -= (1.0-factor_coul)*prefactor;
+          } else ecoul = 0.0;
+          if (rsq < cut_ljsq[itype][jtype]) {
+            evdwl = a[itype][jtype]*rexp - c[itype][jtype]*r6inv -
+              offset[itype][jtype];
+            evdwl *= factor_lj;
+          } else evdwl = 0.0;
+        }
+
+        if (evflag) ev_tally(i,j,nlocal,newton_pair,
+                             evdwl,ecoul,fpair,delx,dely,delz);
+      }
+    }
+  }
+
+  if (vflag_fdotr) virial_fdotr_compute();
+}
diff --git a/src/CORESHELL/pair_buck_coul_long_cs.h b/src/CORESHELL/pair_buck_coul_long_cs.h
new file mode 100644
index 0000000000..ea715fef44
--- /dev/null
+++ b/src/CORESHELL/pair_buck_coul_long_cs.h
@@ -0,0 +1,67 @@
+/* ----------------------------------------------------------------------
+   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: Hendrik Heenen (hendrik.heenen@mytum.com)
+------------------------------------------------------------------------- */
+
+#ifdef PAIR_CLASS
+
+PairStyle(buck/coul/long/cs,PairBuckCoulLongCS)
+
+#else
+
+#ifndef LMP_PAIR_BUCK_COUL_LONG_CS_H
+#define LMP_PAIR_BUCK_COUL_LONG_CS_H
+
+#include "pair_buck_coul_long.h"
+
+namespace LAMMPS_NS {
+
+class PairBuckCoulLongCS : public PairBuckCoulLong {
+ public:
+  PairBuckCoulLongCS(class LAMMPS *);
+  virtual void compute(int, int);
+};
+
+}
+
+#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.
+
+E: All pair coeffs are not set
+
+All pair coefficients must be set in the data file or by the
+pair_coeff command before running a simulation.
+
+E: Pair style buck/coul/long requires atom attribute q
+
+The atom style defined does not have these attributes.
+
+E: Pair style requires a KSpace style
+
+No kspace style is defined.
+
+*/