integrate pair style tersoff/mod/c contributed by Ganga P Purja Pun (GMU)

This includes docs, added testing and inclusion of USER-OMP support.

doc/src/Eqs/pair_tersoff_mod_c.tex

@@ -0,0 +1,10 @@
+\documentclass[12pt]{article}
+\pagestyle{empty}
+
+\begin{document}
+
+\begin{eqnarray*}
+  V_{ij}  & =  & f_C(r_{ij}) \left[ f_R(r_{ij}) + b_{ij} f_A(r_{ij} + c_0) \right]
+\end{eqnarray*}
+
+\end{document}

doc/src/pair_tersoff_mod.txt

@@ -7,32 +7,43 @@
 :line
 
 pair_style tersoff/mod command :h3
+pair_style tersoff/mod/c command :h3
 pair_style tersoff/mod/gpu command :h3
 pair_style tersoff/mod/kk command :h3
 pair_style tersoff/mod/omp command :h3
+pair_style tersoff/mod/c/omp command :h3
 
 [Syntax:]
 
 pair_style tersoff/mod :pre
 
+pair_style tersoff/mod/c :pre
+
 [Examples:]
 
 pair_style tersoff/mod
 pair_coeff * * Si.tersoff.mod Si Si :pre
 
+pair_style tersoff/mod/c
+pair_coeff * * Si.tersoff.modc Si Si :pre
+
 [Description:]
 
-The {tersoff/mod} style computes a bond-order type interatomic
-potential "(Kumagai)"_#Kumagai based on a 3-body Tersoff potential
-"(Tersoff_1)"_#Tersoff_1, "(Tersoff_2)"_#Tersoff_2 with modified
-cutoff function and angular-dependent term, giving the energy E of a
-system of atoms as
+The {tersoff/mod} and {tersoff/mod/c} styles computes a bond-order type
+interatomic potential "(Kumagai)"_#Kumagai based on a 3-body Tersoff
+potential "(Tersoff_1)"_#Tersoff_1, "(Tersoff_2)"_#Tersoff_2 with
+modified cutoff function and angular-dependent term, giving the energy
+E of a system of atoms as
 
 :c,image(Eqs/pair_tersoff_mod.jpg)
 
 where f_R is a two-body term and f_A includes three-body interactions.
 The summations in the formula are over all neighbors J and K of atom I
 within a cutoff distance = R + D.
+The {tersoff/mod/c} style differs from {tersoff/mod} only in the
+formulation of the V_ij term, where it contains an additional c0 term.
+:c,image(Eqs/pair_tersoff_mod_c.jpg)
+
 
 The modified cutoff function f_C proposed by "(Murty)"_#Murty and
 having a continuous second-order differential is employed. The
@@ -69,10 +80,11 @@ are placeholders for atom types that will be used with other
 potentials.
 
 Tersoff/MOD file in the {potentials} directory of the LAMMPS
-distribution have a ".tersoff.mod" suffix.  Lines that are not blank
-or comments (starting with #) define parameters for a triplet of
-elements.  The parameters in a single entry correspond to coefficients
-in the formula above:
+distribution have a ".tersoff.mod" suffix. Potential files for the
+{tersoff/mod/c} style have the suffix ".tersoff.modc". Lines that are
+not blank or comments (starting with #) define parameters for a triplet
+of elements.  The parameters in a single entry correspond to
+coefficients in the formulae above:
 
 element 1 (the center atom in a 3-body interaction)
 element 2 (the atom bonded to the center atom)
@@ -93,13 +105,15 @@ c1
 c2
 c3
 c4
-c5 :ul
+c5
+c0 (energy units, tersoff/mod/c only):ul
 
 The n, eta, lambda2, B, lambda1, and A parameters are only used for
 two-body interactions.  The beta, alpha, c1, c2, c3, c4, c5, h
 parameters are only used for three-body interactions. The R and D
-parameters are used for both two-body and three-body interactions. The
-non-annotated parameters are unitless.
+parameters are used for both two-body and three-body interactions.
+The c0 term applies to {tersoff/mod/c} only. The non-annotated
+parameters are unitless.
 
 The Tersoff/MOD potential file must contain entries for all the elements
 listed in the pair_coeff command.  It can also contain entries for

examples/threebody/in.threebody

@@ -114,3 +114,35 @@ neighbor        1.0 bin
 neigh_modify    every 1 delay 10 check yes
 run             100
 
+# Test Tersoff/Mod model for Si
+
+clear
+read_restart	restart.equil
+
+pair_style      tersoff/mod
+pair_coeff 	* * Si.tersoff.mod Si Si Si Si Si Si Si Si
+
+thermo          10
+fix             1 all nvt temp $t $t 0.1
+fix_modify 	1 energy yes
+timestep        1.0e-3
+neighbor        1.0 bin
+neigh_modify    every 1 delay 10 check yes
+run             100
+
+# Test Tersoff/Mod/C model for Si
+
+clear
+read_restart	restart.equil
+
+pair_style      tersoff/mod/c
+pair_coeff 	* * Si.tersoff.modc Si Si Si Si Si Si Si Si
+
+thermo          10
+fix             1 all nvt temp $t $t 0.1
+fix_modify 	1 energy yes
+timestep        1.0e-3
+neighbor        1.0 bin
+neigh_modify    every 1 delay 10 check yes
+run             100
+

src/MANYBODY/pair_tersoff.h

@@ -49,6 +49,7 @@ class PairTersoff : public Pair {
     double ZBLcut,ZBLexpscale;
     double c5,ca1,ca4;           // added for TersoffMOD
     double powern_del;
+    double c0;                   // added for TersoffMODC
   };
 
   Param *params;                // parameter set for an I-J-K interaction

src/MANYBODY/pair_tersoff_mod.h

@@ -30,7 +30,7 @@ class PairTersoffMOD : public PairTersoff {
   ~PairTersoffMOD() {}
 
  protected:
-  void read_file(char *);
+  virtual void read_file(char *);
   virtual void setup_params();
   double zeta(Param *, double, double, double *, double *);
 

src/MANYBODY/pair_tersoff_mod_c.cpp

@@ -0,0 +1,196 @@
+/* ----------------------------------------------------------------------
+   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: Ganga P Purja Pun (George Mason University, Fairfax)
+------------------------------------------------------------------------- */
+
+#include <math.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include "pair_tersoff_mod_c.h"
+#include "atom.h"
+#include "neighbor.h"
+#include "neigh_list.h"
+#include "neigh_request.h"
+#include "force.h"
+#include "comm.h"
+#include "memory.h"
+#include "error.h"
+
+#include "math_const.h"
+
+using namespace LAMMPS_NS;
+using namespace MathConst;
+
+#define MAXLINE 1024
+#define DELTA 4
+
+/* ---------------------------------------------------------------------- */
+
+void PairTersoffMODC::read_file(char *file)
+{
+  int params_per_line = 21;
+  char **words = new char*[params_per_line+1];
+
+  memory->sfree(params);
+  params = NULL;
+  nparams = maxparam = 0;
+
+  // open file on proc 0
+
+  FILE *fp;
+  if (comm->me == 0) {
+    fp = force->open_potential(file);
+    if (fp == NULL) {
+      char str[128];
+      sprintf(str,"Cannot open Tersoff potential file %s",file);
+      error->one(FLERR,str);
+    }
+  }
+
+  // read each line out of file, skipping blank lines or leading '#'
+  // store line of params if all 3 element tags are in element list
+
+  int n,nwords,ielement,jelement,kelement;
+  char line[MAXLINE],*ptr;
+  int eof = 0;
+
+  while (1) {
+    if (comm->me == 0) {
+      ptr = fgets(line,MAXLINE,fp);
+      if (ptr == NULL) {
+	    eof = 1;
+	    fclose(fp);
+      } else n = strlen(line) + 1;
+    }
+    MPI_Bcast(&eof,1,MPI_INT,0,world);
+    if (eof) break;
+    MPI_Bcast(&n,1,MPI_INT,0,world);
+    MPI_Bcast(line,n,MPI_CHAR,0,world);
+
+    // strip comment, skip line if blank
+
+    if ((ptr = strchr(line,'#'))) *ptr = '\0';
+    nwords = atom->count_words(line);
+    if (nwords == 0) continue;
+
+    // concatenate additional lines until have params_per_line words
+
+    while (nwords < params_per_line) {
+      n = strlen(line);
+      if (comm->me == 0) {
+        ptr = fgets(&line[n],MAXLINE-n,fp);
+        if (ptr == NULL) {
+	      eof = 1;
+	      fclose(fp);
+        } else n = strlen(line) + 1;
+      }
+      MPI_Bcast(&eof,1,MPI_INT,0,world);
+      if (eof) break;
+      MPI_Bcast(&n,1,MPI_INT,0,world);
+      MPI_Bcast(line,n,MPI_CHAR,0,world);
+      if ((ptr = strchr(line,'#'))) *ptr = '\0';
+      nwords = atom->count_words(line);
+    }
+
+    if (nwords != params_per_line)
+      error->all(FLERR,"Incorrect format in Tersoff potential file");
+
+    // words = ptrs to all words in line
+
+    nwords = 0;
+    words[nwords++] = strtok(line," \t\n\r\f");
+    while ((words[nwords++] = strtok(NULL," \t\n\r\f"))) continue;
+
+    // ielement,jelement,kelement = 1st args
+    // if all 3 args are in element list, then parse this line
+    // else skip to next line
+
+    for (ielement = 0; ielement < nelements; ielement++)
+      if (strcmp(words[0],elements[ielement]) == 0) break;
+    if (ielement == nelements) continue;
+    for (jelement = 0; jelement < nelements; jelement++)
+      if (strcmp(words[1],elements[jelement]) == 0) break;
+    if (jelement == nelements) continue;
+    for (kelement = 0; kelement < nelements; kelement++)
+      if (strcmp(words[2],elements[kelement]) == 0) break;
+    if (kelement == nelements) continue;
+
+    // load up parameter settings and error check their values
+
+    if (nparams == maxparam) {
+      maxparam += DELTA;
+      params = (Param *) memory->srealloc(params,maxparam*sizeof(Param),
+					  "pair:params");
+    }
+
+    params[nparams].ielement = ielement;
+    params[nparams].jelement = jelement;
+    params[nparams].kelement = kelement;
+    params[nparams].powerm = atof(words[3]);
+    params[nparams].lam3 = atof(words[4]);
+    params[nparams].h = atof(words[5]);
+    params[nparams].powern = atof(words[6]);
+    params[nparams].beta = atof(words[7]);
+    params[nparams].lam2 = atof(words[8]);
+    params[nparams].bigb = atof(words[9]);
+    params[nparams].bigr = atof(words[10]);
+    params[nparams].bigd = atof(words[11]);
+    params[nparams].lam1 = atof(words[12]);
+    params[nparams].biga = atof(words[13]);
+    params[nparams].powern_del = atof(words[14]);
+    params[nparams].c1 = atof(words[15]);
+    params[nparams].c2 = atof(words[16]);
+    params[nparams].c3 = atof(words[17]);
+    params[nparams].c4 = atof(words[18]);
+    params[nparams].c5 = atof(words[19]);
+    params[nparams].c0 = atof(words[20]);
+
+    // currently only allow m exponent of 1
+
+    params[nparams].powermint = int(params[nparams].powerm);
+
+    if (
+	params[nparams].lam3 < 0.0 || params[nparams].powern < 0.0 ||
+	params[nparams].beta < 0.0 || params[nparams].lam2 < 0.0 ||
+	params[nparams].bigb < 0.0 || params[nparams].bigr < 0.0 ||
+	params[nparams].bigd < 0.0 ||
+                               params[nparams].bigd > params[nparams].bigr ||
+	params[nparams].lam3 < 0.0 || params[nparams].biga < 0.0 ||
+	params[nparams].powerm - params[nparams].powermint != 0.0 ||
+    (params[nparams].powermint != 3 && params[nparams].powermint != 1))
+      error->all(FLERR,"Illegal Tersoff parameter");
+
+    nparams++;
+  }
+
+  delete [] words;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void PairTersoffMODC::repulsive(Param *param, double rsq, double &fforce,
+                                int eflag, double &eng)
+{
+  double r,tmp_fc,tmp_fc_d,tmp_exp;
+
+  r = sqrt(rsq);
+  tmp_fc = ters_fc(r,param);
+  tmp_fc_d = ters_fc_d(r,param);
+  tmp_exp = exp(-param->lam1 * r);
+  fforce = -param->biga * tmp_exp * (tmp_fc_d - tmp_fc * param->lam1) / r - param->c0 * tmp_fc_d / r;
+  if (eflag) eng = tmp_fc * param->biga * tmp_exp + param->c0 * tmp_fc;
+}
+

src/MANYBODY/pair_tersoff_mod_c.h

@@ -0,0 +1,66 @@
+/* -*- c++ -*- ----------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   http://lammps.sandia.gov, Sandia National Laboratories
+   Steve Plimpton, sjplimp@sandia.gov
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#ifdef PAIR_CLASS
+
+PairStyle(tersoff/mod/c,PairTersoffMODC)
+
+#else
+
+#ifndef LMP_PAIR_TERSOFF_MOD_C_H
+#define LMP_PAIR_TERSOFF_MOD_C_H
+
+#include "pair_tersoff_mod.h"
+
+namespace LAMMPS_NS {
+
+class PairTersoffMODC : public PairTersoffMOD {
+ public:
+  PairTersoffMODC(class LAMMPS *lmp) : PairTersoffMOD(lmp) {};
+  ~PairTersoffMODC() {}
+
+ protected:
+  void read_file(char *);
+  void repulsive(Param *, double, double &, int, double &);
+};
+
+}
+
+#endif
+#endif
+
+/* ERROR/WARNING messages:
+
+E: Cannot open Tersoff potential file %s
+
+The specified potential file cannot be opened.  Check that the path
+and name are correct.
+
+E: Incorrect format in Tersoff potential file
+
+Incorrect number of words per line in the potential file.
+
+E: Illegal Tersoff parameter
+
+One or more of the coefficients defined in the potential file is
+invalid.
+
+E: Potential file has duplicate entry
+
+The potential file has more than one entry for the same element.
+
+E: Potential file is missing an entry
+
+The potential file does not have a needed entry.
+
+*/

src/USER-OMP/pair_tersoff_mod_c_omp.cpp

@@ -0,0 +1,253 @@
+/* ----------------------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   http://lammps.sandia.gov, Sandia National Laboratories
+   Steve Plimpton, sjplimp@sandia.gov
+
+   This software is distributed under the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+/* ----------------------------------------------------------------------
+   Contributing author: Axel Kohlmeyer (Temple U)
+------------------------------------------------------------------------- */
+
+#include <math.h>
+#include "pair_tersoff_mod_c_omp.h"
+#include "atom.h"
+#include "comm.h"
+#include "force.h"
+#include "neighbor.h"
+#include "neigh_list.h"
+
+#include "suffix.h"
+using namespace LAMMPS_NS;
+
+/* ---------------------------------------------------------------------- */
+
+PairTersoffMODCOMP::PairTersoffMODCOMP(LAMMPS *lmp) :
+  PairTersoffMODC(lmp), ThrOMP(lmp, THR_PAIR)
+{
+  suffix_flag |= Suffix::OMP;
+  respa_enable = 0;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void PairTersoffMODCOMP::compute(int eflag, int vflag)
+{
+  if (eflag || vflag) {
+    ev_setup(eflag,vflag);
+  } else evflag = vflag_fdotr = vflag_atom = 0;
+
+  const int nall = atom->nlocal + atom->nghost;
+  const int nthreads = comm->nthreads;
+  const int inum = list->inum;
+
+#if defined(_OPENMP)
+#pragma omp parallel default(none) shared(eflag,vflag)
+#endif
+  {
+    int ifrom, ito, tid;
+
+    loop_setup_thr(ifrom, ito, tid, inum, nthreads);
+    ThrData *thr = fix->get_thr(tid);
+    thr->timer(Timer::START);
+    ev_setup_thr(eflag, vflag, nall, eatom, vatom, thr);
+
+    if (evflag) {
+      if (eflag) {
+        if (vflag_atom) eval<1,1,1>(ifrom, ito, thr);
+        else eval<1,1,0>(ifrom, ito, thr);
+      } else {
+        if (vflag_atom) eval<1,0,1>(ifrom, ito, thr);
+        else eval<1,0,0>(ifrom, ito, thr);
+      }
+    } else eval<0,0,0>(ifrom, ito, thr);
+
+    thr->timer(Timer::PAIR);
+    reduce_thr(this, eflag, vflag, thr);
+  } // end of omp parallel region
+}
+
+template <int EVFLAG, int EFLAG, int VFLAG_ATOM>
+void PairTersoffMODCOMP::eval(int iifrom, int iito, ThrData * const thr)
+{
+  int i,j,k,ii,jj,kk,jnum;
+  tagint itag,jtag;
+  int itype,jtype,ktype,iparam_ij,iparam_ijk;
+  double xtmp,ytmp,ztmp,delx,dely,delz,evdwl,fpair;
+  double rsq,rsq1,rsq2;
+  double delr1[3],delr2[3],fi[3],fj[3],fk[3];
+  double zeta_ij,prefactor;
+  int *ilist,*jlist,*numneigh,**firstneigh;
+
+  evdwl = 0.0;
+
+  const dbl3_t * _noalias const x = (dbl3_t *) atom->x[0];
+  dbl3_t * _noalias const f = (dbl3_t *) thr->get_f()[0];
+  const tagint * _noalias const tag = atom->tag;
+  const int * _noalias const type = atom->type;
+  const int nlocal = atom->nlocal;
+
+  ilist = list->ilist;
+  numneigh = list->numneigh;
+  firstneigh = list->firstneigh;
+
+  double fxtmp,fytmp,fztmp;
+
+  // loop over full neighbor list of my atoms
+
+  for (ii = iifrom; ii < iito; ++ii) {
+
+    i = ilist[ii];
+    itag = tag[i];
+    itype = map[type[i]];
+    xtmp = x[i].x;
+    ytmp = x[i].y;
+    ztmp = x[i].z;
+    fxtmp = fytmp = fztmp = 0.0;
+
+    // two-body interactions, skip half of them
+
+    jlist = firstneigh[i];
+    jnum = numneigh[i];
+
+    for (jj = 0; jj < jnum; jj++) {
+      j = jlist[jj];
+      j &= NEIGHMASK;
+      jtag = tag[j];
+
+      if (itag > jtag) {
+        if ((itag+jtag) % 2 == 0) continue;
+      } else if (itag < jtag) {
+        if ((itag+jtag) % 2 == 1) continue;
+      } else {
+        if (x[j].z < ztmp) continue;
+        if (x[j].z == ztmp && x[j].y < ytmp) continue;
+        if (x[j].z == ztmp && x[j].y == ytmp && x[j].x < xtmp) continue;
+      }
+
+      jtype = map[type[j]];
+
+      delx = xtmp - x[j].x;
+      dely = ytmp - x[j].y;
+      delz = ztmp - x[j].z;
+      rsq = delx*delx + dely*dely + delz*delz;
+
+      iparam_ij = elem2param[itype][jtype][jtype];
+      if (rsq > params[iparam_ij].cutsq) continue;
+
+      repulsive(&params[iparam_ij],rsq,fpair,EFLAG,evdwl);
+
+      fxtmp += delx*fpair;
+      fytmp += dely*fpair;
+      fztmp += delz*fpair;
+      f[j].x -= delx*fpair;
+      f[j].y -= dely*fpair;
+      f[j].z -= delz*fpair;
+
+      if (EVFLAG) ev_tally_thr(this,i,j,nlocal,/* newton_pair */ 1,
+                               evdwl,0.0,fpair,delx,dely,delz,thr);
+    }
+
+    // three-body interactions
+    // skip immediately if I-J is not within cutoff
+    double fjxtmp,fjytmp,fjztmp;
+
+    for (jj = 0; jj < jnum; jj++) {
+      j = jlist[jj];
+      j &= NEIGHMASK;
+      jtype = map[type[j]];
+      iparam_ij = elem2param[itype][jtype][jtype];
+
+      delr1[0] = x[j].x - xtmp;
+      delr1[1] = x[j].y - ytmp;
+      delr1[2] = x[j].z - ztmp;
+      rsq1 = delr1[0]*delr1[0] + delr1[1]*delr1[1] + delr1[2]*delr1[2];
+      if (rsq1 > params[iparam_ij].cutsq) continue;
+
+      // accumulate bondorder zeta for each i-j interaction via loop over k
+
+      fjxtmp = fjytmp = fjztmp = 0.0;
+      zeta_ij = 0.0;
+
+      for (kk = 0; kk < jnum; kk++) {
+        if (jj == kk) continue;
+        k = jlist[kk];
+        k &= NEIGHMASK;
+        ktype = map[type[k]];
+        iparam_ijk = elem2param[itype][jtype][ktype];
+
+        delr2[0] = x[k].x - xtmp;
+        delr2[1] = x[k].y - ytmp;
+        delr2[2] = x[k].z - ztmp;
+        rsq2 = delr2[0]*delr2[0] + delr2[1]*delr2[1] + delr2[2]*delr2[2];
+        if (rsq2 > params[iparam_ijk].cutsq) continue;
+
+        zeta_ij += zeta(&params[iparam_ijk],rsq1,rsq2,delr1,delr2);
+      }
+
+      // pairwise force due to zeta
+
+      force_zeta(&params[iparam_ij],rsq1,zeta_ij,fpair,prefactor,EFLAG,evdwl);
+
+      fxtmp += delr1[0]*fpair;
+      fytmp += delr1[1]*fpair;
+      fztmp += delr1[2]*fpair;
+      fjxtmp -= delr1[0]*fpair;
+      fjytmp -= delr1[1]*fpair;
+      fjztmp -= delr1[2]*fpair;
+
+      if (EVFLAG) ev_tally_thr(this,i,j,nlocal,/* newton_pair */ 1,evdwl,0.0,
+                               -fpair,-delr1[0],-delr1[1],-delr1[2],thr);
+
+      // attractive term via loop over k
+
+      for (kk = 0; kk < jnum; kk++) {
+        if (jj == kk) continue;
+        k = jlist[kk];
+        k &= NEIGHMASK;
+        ktype = map[type[k]];
+        iparam_ijk = elem2param[itype][jtype][ktype];
+
+        delr2[0] = x[k].x - xtmp;
+        delr2[1] = x[k].y - ytmp;
+        delr2[2] = x[k].z - ztmp;
+        rsq2 = delr2[0]*delr2[0] + delr2[1]*delr2[1] + delr2[2]*delr2[2];
+        if (rsq2 > params[iparam_ijk].cutsq) continue;
+
+        attractive(&params[iparam_ijk],prefactor,
+                   rsq1,rsq2,delr1,delr2,fi,fj,fk);
+
+        fxtmp += fi[0];
+        fytmp += fi[1];
+        fztmp += fi[2];
+        fjxtmp += fj[0];
+        fjytmp += fj[1];
+        fjztmp += fj[2];
+        f[k].x += fk[0];
+        f[k].y += fk[1];
+        f[k].z += fk[2];
+
+        if (VFLAG_ATOM) v_tally3_thr(i,j,k,fj,fk,delr1,delr2,thr);
+      }
+      f[j].x += fjxtmp;
+      f[j].y += fjytmp;
+      f[j].z += fjztmp;
+    }
+    f[i].x += fxtmp;
+    f[i].y += fytmp;
+    f[i].z += fztmp;
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
+double PairTersoffMODCOMP::memory_usage()
+{
+  double bytes = memory_usage_thr();
+  bytes += PairTersoffMOD::memory_usage();
+
+  return bytes;
+}

src/USER-OMP/pair_tersoff_mod_c_omp.h

@@ -0,0 +1,43 @@
+/* -*- c++ -*- ----------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   http://lammps.sandia.gov, Sandia National Laboratories
+   Steve Plimpton, sjplimp@sandia.gov
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+/* ----------------------------------------------------------------------
+   Contributing author: Axel Kohlmeyer (Temple U)
+------------------------------------------------------------------------- */
+
+#ifdef PAIR_CLASS
+
+PairStyle(tersoff/mod/c/omp,PairTersoffMODCOMP)
+
+#else
+
+#ifndef LMP_PAIR_TERSOFF_MOD_C_OMP_H
+#define LMP_PAIR_TERSOFF_MOD_C_OMP_H
+
+#include "pair_tersoff_mod_c.h"
+#include "thr_omp.h"
+
+namespace LAMMPS_NS {
+
+class PairTersoffMODCOMP : public PairTersoffMODC, public ThrOMP {
+
+ public:
+  PairTersoffMODCOMP(class LAMMPS *);
+
+  virtual void compute(int, int);
+  virtual double memory_usage();
+
+ private:
+  template <int EVFLAG, int EFLAG, int VFLAG_ATOM>
+  void eval(int ifrom, int ito, ThrData * const thr);
+};
+
+}
+
+#endif
+#endif