Finally updated the github tutorial

doc/src/Eqs/pair_coul_shield.tex

+\documentclass[aps,pr,onecolumn,superscriptaddress,noshowpacs,a4paper,15pt]{revtex4}
+\pdfoutput=1
+\bibliographystyle{apsrev4}
+\usepackage{color}
+\usepackage{dcolumn} %Align table columns on decimal point
+\usepackage{amssymb}
+\usepackage{amsmath}
+\usepackage{amsthm}
+\usepackage{graphicx}
+\usepackage[pdftex]{hyperref}
+\hypersetup{colorlinks=true,citecolor=blue,linkcolor=red,urlcolor=blue}
+\usepackage[all]{hypcap}
+\newcommand{\red}{\color{red}}
+\newcommand{\blue}{\color{blue}}
+\definecolor{green}{rgb}{0,0.5,0}
+\newcommand{\green}{\color{green}}
+\newcommand{\white}{\color{white}}
+%\newcommand{\cite}[1]{\hspace{-1 ex} % \nocite{#1}\citenum{#1}}
+\thickmuskip=0.5\thickmuskip %shorter spaces in math
+
+\begin{document}
+\begingroup
+\Large
+\begin{eqnarray*}
+  E & = & \frac{1}{2} \sum_i \sum_{j \neq i} V_{ij} \\[15pt]
+  V_{ij} & = & {\rm Tap}(r_{ij})\frac{\kappa q_i q_j}{\sqrt[3]{r_{ij}^3+(1/\lambda_{ij})^3}}\\[15pt]
+  {\rm Tap}(r_{ij}) & = & 20\left ( \frac{r_{ij}}{R_{cut}} \right )^7 -
+                          70\left ( \frac{r_{ij}}{R_{cut}} \right )^6 +
+                          84\left ( \frac{r_{ij}}{R_{cut}} \right )^5 -
+                          35\left ( \frac{r_{ij}}{R_{cut}} \right )^4 + 1
+\end{eqnarray*}
+\endgroup
+\end{document}

doc/src/Eqs/pair_ilp_gr_hBN.tex

+\documentclass[aps,pr,onecolumn,superscriptaddress,noshowpacs,a4paper,15pt]{revtex4}
+\pdfoutput=1
+\bibliographystyle{apsrev4}
+\usepackage{color}
+\usepackage{dcolumn} %Align table columns on decimal point
+\usepackage{amssymb}
+\usepackage{amsmath}
+\usepackage{amsthm}
+\usepackage{graphicx}
+\usepackage[pdftex]{hyperref}
+\hypersetup{colorlinks=true,citecolor=blue,linkcolor=red,urlcolor=blue}
+\usepackage[all]{hypcap}
+\newcommand{\red}{\color{red}}
+\newcommand{\blue}{\color{blue}}
+\definecolor{green}{rgb}{0,0.5,0}
+\newcommand{\green}{\color{green}}
+\newcommand{\white}{\color{white}}
+%\newcommand{\cite}[1]{\hspace{-1 ex} % \nocite{#1}\citenum{#1}}
+\thickmuskip=0.5\thickmuskip %shorter spaces in math
+%
+\begin{document}
+%
+\begingroup
+\Large
+\begin{eqnarray*}
+  E & = & \frac{1}{2} \sum_i \sum_{j \neq i} V_{ij} \\[15pt]
+  V_{ij} & = & {\rm Tap}(r_{ij})\left \{ e^{-\alpha (r_{ij}/\beta -1)} 
+               \left [ \epsilon + f(\rho_{ij}) + f(\rho_{ji})\right ] - 
+                \frac{1}{1+e^{-d\left [ \left ( r_{ij}/\left (s_R \cdot r^{eff} \right ) \right )-1 \right ]}}
+                \cdot \frac{C_6}{r^6_{ij}} \right \}\\[15pt]
+  \rho_{ij}^2 & = & r_{ij}^2 - ({\bf r}_{ij} \cdot {\bf n}_i)^2 \\[15pt]
+  \rho_{ji}^2 & = & r_{ij}^2 - ({\bf r}_{ij} \cdot {\bf n}_j)^2 \\[15pt]
+  f(\rho) & = &  C e^{ -( \rho / \delta )^2 }\\[15pt]
+  {\rm Tap}(r_{ij}) & = & 20\left ( \frac{r_{ij}}{R_{cut}} \right )^7 -
+                          70\left ( \frac{r_{ij}}{R_{cut}} \right )^6 +
+                          84\left ( \frac{r_{ij}}{R_{cut}} \right )^5 -
+                          35\left ( \frac{r_{ij}}{R_{cut}} \right )^4 + 1
+\end{eqnarray*}
+\endgroup
+%
+\end{document}
+%

doc/src/Eqs/pair_kolmogorov_crespi_full.tex

+\documentclass[aps,pr,onecolumn,superscriptaddress,noshowpacs,a4paper,15pt]{revtex4}
+\pdfoutput=1
+\bibliographystyle{apsrev4}
+\usepackage{color}
+\usepackage{dcolumn} %Align table columns on decimal point
+\usepackage{amssymb}
+\usepackage{amsmath}
+\usepackage{amsthm}
+\usepackage{graphicx}
+\usepackage[pdftex]{hyperref}
+\hypersetup{colorlinks=true,citecolor=blue,linkcolor=red,urlcolor=blue}
+\usepackage[all]{hypcap}
+\newcommand{\red}{\color{red}}
+\newcommand{\blue}{\color{blue}}
+\definecolor{green}{rgb}{0,0.5,0}
+\newcommand{\green}{\color{green}}
+\newcommand{\white}{\color{white}}
+%\newcommand{\cite}[1]{\hspace{-1 ex} % \nocite{#1}\citenum{#1}}
+\thickmuskip=0.5\thickmuskip %shorter spaces in math
+
+\begin{document}
+
+\begingroup
+\Large
+\begin{eqnarray*}
+  E & = & \frac{1}{2} \sum_i \sum_{j \neq i} V_{ij} \\[15pt]
+  V_{ij} & = & e^{-\lambda (r_{ij} -z_0)} \left [ C + f(\rho_{ij}) + f(\rho_{ji}) - A \left ( \frac{r_{ij}}{z_0}\right )^{-6} \right ] \\
+  \rho_{ij}^2 & = & r_{ij}^2 - ({\bf r}_{ij}\cdot {\bf n}_{i})^2 \\[15pt]
+  \rho_{ji}^2 & = & r_{ij}^2 - ({\bf r}_{ij}\cdot  {\bf n}_{j})^2 \\[15pt]
+  f(\rho) & = &  e^{-(\rho/\delta)^2} \sum_{n=0}^2 C_{2n} { \rho/\delta }^{2n}
+\end{eqnarray*}
+\endgroup
+\end{document}

doc/src/pair_coul_shield.txt

+"LAMMPS WWW Site"_lws - "LAMMPS Documentation"_ld - "LAMMPS Commands"_lc :c
+
+:link(lws,http://lammps.sandia.gov)
+:link(ld,Manual.html)
+:link(lc,Section_commands.html#comm)
+
+:line
+
+pair_style coul/shield command :h3
+
+[Syntax:]
+
+pair_style coul/shield cutoff tap_flag :pre
+
+cutoff = global cutoff (distance units)
+tap_flag = 0/1 to turn off/on the Tapper function
+
+[Examples:]
+
+pair_style coul/shield 16.0 1
+pair_coeff 1 2 0.70 :pre
+
+
+[Description:]
+
+Style {coul/shield} computes a Coulomb interaction for boron and nitrigon
+atoms locate in different layers of hexagonal boron nitride. This potential must
+be used in combination with the potential "pair_style pair_ilp_gr_hBN"_pair_ilp_gr_hBN.html
+
+NOTE: This potential is intended for electrostatic interactions between two different
+layers of hexagonal boron nitride. Therefore, to avoid interaction within the same layers,
+each layer should have a separate molecule id and is recommended to use
+"full" atom style in the data file.
+
+:c,image(Eqs/pair_coul_shield.jpg)
+
+Where Tap(r_ij) is the tapper function which provides a continuous cutoff (up to third derivative)
+for interatomic separations larger than r_c "(Maaravi)"_#Maaravi. Here \lambda is the shielding
+parameter that eliminates the short-range sigularity of the classical monopolar electrostatic
+interaction expression "(Maaravi)"_#Maaravi.
+
+The shielding parameter \lambda (1/distance units) must be defined for each pair of atom
+types via the "pair_coeff"_pair_coeff.html command as in the example
+above, or in the data file or restart files read by the
+"read_data"_read_data.html or "read_restart"_read_restart.html
+commands:
+
+The global cutoff (r_c) specified in the pair_style command is used.
+
+:line
+
+[Mixing, shift, table, tail correction, restart, rRESPA info]:
+
+This pair style does not support parameter mixing. Coefficients must be given explicitly for each type of particle pairs.
+
+The "pair_modify"_pair_modify.html table option is not relevant
+for this pair style.
+
+This pair style does not support the "pair_modify"_pair_modify.html
+tail option for adding long-range tail corrections to energy and
+pressure.
+
+This pair style can only be used via the {pair} keyword of the
+"run_style respa"_run_style.html command.  It does not support the
+{inner}, {middle}, {outer} keywords.
+
+[Restrictions:]
+
+This style is part of the "user-misc" package.  It is only enabled
+if LAMMPS was built with that package.  See the "Making
+LAMMPS"_Section_start.html#start_2_3 section for more info.
+
+[Related commands:]
+
+"pair_coeff"_pair_coeff.html
+"pair_style pair_ilp_gr_hBN"_pair_ilp_gr_hBN.html
+
+[Default:] tap_flag = 1
+
+:line
+
+:link(Maaravi)
+[(Maaravi)] T. Maaravi et al, J. Phys. Chem. C 121, 22826-22835 (2017).

doc/src/pair_ilp_gr_hBN.txt

+"LAMMPS WWW Site"_lws - "LAMMPS Documentation"_ld - "LAMMPS Commands"_lc :c
+
+:link(lws,http://lammps.sandia.gov)
+:link(ld,Manual.html)
+:link(lc,Section_commands.html#comm)
+
+:line
+
+pair_style ILP/graphene/hBN command :h3
+
+[Syntax:]
+
+pair_style hybrid/overlay ILP/graphene/hBN cutoff tap_flag :pre
+
+cutoff = global cutoff (distance units)
+tap_flag = 0/1 to turn off/on the Tapper function 
+
+[Examples:]
+
+pair_style hybrid/overlay ILP/graphene/hBN 16.0 1
+pair_coeff * * ILP/graphene/hBN  BNCH.ILP B N C :pre
+
+pair_style  hybrid/overlay rebo tersoff ILP/graphene/hBN 16.0 coul/shield 16.0
+pair_coeff  * * rebo                 CH.airebo   NULL NULL C
+pair_coeff  * * tersoff              BNC.tersoff B    N    NULL
+pair_coeff  * * ILP/graphene/hBN     BNCH.ILP    B    N    C :pre
+pair_coeff  1 1 coul/shield 0.70
+pair_coeff  1 2 coul/shield 0.69498201415576216335
+pair_coeff  2 2 coul/shield 0.69
+
+[Description:]
+
+The {ILP/graphene/hBN} style computes the registry-dependent interlayer potential (RDILP)
+potential as described in "(Leven)"_#Leven and "(Maaravi)"_#Maaravi. The normals are calculated in 
+the way as described in "(Kolmogorov)"_#Kolmogorov.
+
+:c,image(Eqs/pair_ilp_gr_hBN.jpg)
+
+Where Tap(r_ij) is the tapper function which provides a continuous cutoff (up to third derivative)
+for interatomic separations larger than r_c "(Maaravi)"_#Maaravi. The definitons of each parameter
+in the above equation can be found in "(Leven)"_#Leven and "(Maaravi)"_#Maaravi.
+
+It is important to include all the pairs to build the neighbor list for
+calculating the normals.
+
+NOTE: This potential is intended for interactions between two different layers
+of graphene or hexagonal boron nitride. Therefore, to avoid interaction within
+the same layers, each layer should have a separate molecule id and is recommended
+to use "full" atom style in the data file.
+
+The parameter file (e.g. BNCH.ILP), is intended for use with metal
+"units"_units.html, with energies in meV. Two additional parameters, {S},
+and {rcut} are included in the parameter file. {S} is designed to 
+facilitate scaling of energies. {rcut} is designed to build the neighbor
+list for calculating the normals for each atom pair.
+
+NOTE: The parameters presented in the parameter file (e.g. BNCH.ILP), are
+fitted with tapper function by setting the cutoff equal to 16.0 Angstrom.
+Using different cutoff or tapper function should be careful.
+
+This potential must be used in combination with hybrid/overlay.
+Other interactions can be set to zero using pair_style {none}.
+
+:line
+
+[Mixing, shift, table, tail correction, restart, rRESPA info]:
+
+This pair style does not support the pair_modify mix, shift, table, and tail options.
+
+This pair style does not write their information to binary restart files, 
+since it is stored in potential files. Thus, you need to re-specify the 
+pair_style and pair_coeff commands in an input script that reads a restart file.
+
+[Restrictions:]
+
+This fix is part of the USER-MISC package.  It is only enabled if
+LAMMPS was built with that package.  See the "Making
+LAMMPS"_Section_start.html#start_3 section for more info.
+
+This pair potential requires the newton setting to be “on” for pair interactions.
+
+The BNCH.ILP potential file provided with LAMMPS (see the potentials directory) are parameterized for metal units. 
+You can use this potential with any LAMMPS units, but you would need to create your BNCH.ILP potential file with 
+coefficients listed in the appropriate units, if your simulation doesn’t use “metal” units.
+
+[Related commands:]
+
+"pair_coeff"_pair_coeff.html
+"pair_none"_pair_none.html
+"pair_style hybrid/overlay"_pair_hybrid.html
+"pair_style pair_kolmogorov_crespi_z"_pair_kolmogorov_crespi_z.html
+"pair_style pair_kolmogorov_crespi_full"_pair_kolmogorov_crespi_full.html
+"pair_style pair_coul_shield"_pair_coul_shield.html
+
+[Default:] tap_flag = 1
+
+:line
+
+:link(Leven)
+[(Leven)] I. Leven et al, J. Chem.Theory Comput. 12, 2896-905 (2016)
+
+:link(Maaravi)
+[(Maaravi)] T. Maaravi et al, J. Phys. Chem. C 121, 22826-22835 (2017).
+
+:link(Kolmogorov)
+[(Kolmogorov)] A. N. Kolmogorov, V. H. Crespi, Phys. Rev. B 71, 235415 (2005)

doc/src/pair_kolmogorov_crespi_full.txt

+"LAMMPS WWW Site"_lws - "LAMMPS Documentation"_ld - "LAMMPS Commands"_lc :c
+
+:link(lws,http://lammps.sandia.gov)
+:link(ld,Manual.html)
+:link(lc,Section_commands.html#comm)
+
+:line
+
+pair_style kolmogorov/crespi/full command :h3
+
+[Syntax:]
+
+pair_style hybrid/overlay kolmogorov/crespi/full cutoff tap_flag:pre
+
+cutoff = global cutoff (distance units)
+tap_flag = 0/1 to turn off/on the Tapper function 
+
+[Examples:]
+
+pair_style hybrid/overlay kolmogorov/crespi/full 20.0 0
+pair_coeff * * none
+pair_coeff * * kolmogorov/crespi/full  CC.KC   C C :pre
+
+pair_style hybrid/overlay rebo kolmogorov/crespi/full 16.0
+pair_coeff * * rebo                    CH.airebo  C C
+pair_coeff * * kolmogorov/crespi/full  CC.KC      C C :pre
+
+[Description:]
+
+The {kolmogorov/crespi/full} style computes the Kolmogorov-Crespi interaction
+potential as described in "(Kolmogorov)"_#Kolmogorov. No simplification is made,
+
+:c,image(Eqs/pair_kolmogorov_crespi_full.jpg)
+
+It is important to have a suffiently large cutoff to ensure smooth forces and
+to include all the pairs to build the neighbor list for calculating the normals.
+Energies are shifted so that they go continously to zero at the cutoff assuming
+that the exponential part of {Vij} (first term) decays sufficiently fast.
+This shift is achieved by the last term in the equation for {Vij} above.
+
+NOTE: This potential is intended for interactions between two different graphene
+layers. Therefore, to avoid interaction within the same layers, each layer
+should have a separate molecule id and is recommended to use
+"full" atom style in the data file.
+
+The parameter file (e.g. CC.KC), is intended for use with metal
+"units"_units.html, with energies in meV. Two additional parameters, {S},
+and {rcut} are included in the parameter file. {S} is designed to
+facilitate scaling of energies. {rcut} is designed to build the neighbor
+list for calculating the normals for each atom pair.
+
+This potential must be used in combination with hybrid/overlay.
+Other interactions can be set to zero using pair_style {none}.
+
+:line
+
+[Mixing, shift, table, tail correction, restart, rRESPA info]:
+
+This pair style does not support the pair_modify mix, shift, table, and tail options.
+
+This pair style does not write their information to binary restart files, 
+since it is stored in potential files. Thus, you need to re-specify the
+pair_style and pair_coeff commands in an input script that reads a restart file.
+
+[Restrictions:]
+
+This fix is part of the USER-MISC package.  It is only enabled if
+LAMMPS was built with that package.  See the "Making
+LAMMPS"_Section_start.html#start_3 section for more info.
+
+This pair potential requires the newton setting to be “on” for pair interactions.
+
+The CC.KC potential file provided with LAMMPS (see the potentials directory) are parameterized for metal units. 
+You can use this potential with any LAMMPS units, but you would need to create your CC.KC potential file with 
+coefficients listed in the appropriate units, if your simulation doesn’t use “metal” units.
+
+[Related commands:]
+
+"pair_coeff"_pair_coeff.html
+"pair_none"_pair_none.html
+"pair_style hybrid/overlay"_pair_hybrid.html
+"pair_style pair_kolmogorov_crespi_z"_pair_kolmogorov_crespi_z.html
+"pair_style pair_ilp_gr_hBN"_pair_ilp_gr_hBN.html
+
+[Default:] tap_flag = 0
+
+:line
+
+:link(Kolmogorov)
+[(Kolmogorov)] A. N. Kolmogorov, V. H. Crespi, Phys. Rev. B 71, 235415 (2005)

src/USER-MISC/pair_coul_shield.cpp

+/* ----------------------------------------------------------------------
+   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: Wengen Ouyang (Tel Aviv University)
+   e-mail: w.g.ouyang at gmail dot com
+  
+   This is a Coulomb potential described in
+   [Maaravi et al, J. Phys. Chem. C 121, 22826-22835 (2017)]
+------------------------------------------------------------------------- */
+
+#include <math.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include "pair_coul_shield.h"
+#include "atom.h"
+#include "comm.h"
+#include "force.h"
+#include "neighbor.h"
+#include "neigh_list.h"
+#include "memory.h"
+#include "error.h"
+
+using namespace LAMMPS_NS;
+
+/* ---------------------------------------------------------------------- */
+
+PairCoulshield::PairCoulshield(LAMMPS *lmp) : Pair(lmp) {
+Tapflag = 1;
+}
+
+/* ---------------------------------------------------------------------- */
+
+PairCoulshield::~PairCoulshield()
+{
+  if (allocated) {
+    memory->destroy(setflag);
+    memory->destroy(sigmae);
+    memory->destroy(offset);
+    memory->destroy(cutsq);
+    memory->destroy(cut);
+    allocated = 0;
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void PairCoulshield::compute(int eflag, int vflag)
+{
+  int i,j,ii,jj,inum,jnum,itype,jtype;
+  double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,ecoul,fpair,Tap,dTap;
+  double rsq,r,r3,rarg,th,depsdr,epsr,forcecoul,factor_coul,Vc,fvc;
+  int *ilist,*jlist,*numneigh,**firstneigh;
+
+  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;
+  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_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];
+
+      // only include the interation between different layers
+      if (rsq < cutsq[itype][jtype] && atom->molecule[i] != atom->molecule[j]) {
+        r = sqrt(rsq);
+	r3 = rsq*r;
+        rarg = 1.0/sigmae[itype][jtype];
+        th = r3 + pow(rarg,3.0);
+        epsr = 1.0/pow(th,0.333333333333333333333333);
+        depsdr = pow(epsr,4.0);
+	Vc = qqrd2e*qtmp*q[j]*epsr;
+
+	// turn on/off Tapper function
+        if (Tapflag) {
+          Tap = calc_Tap(r,sqrt(cutsq[itype][jtype]));
+          dTap = calc_dTap(r,sqrt(cutsq[itype][jtype]));
+        } else {Tap = 1.0; dTap = 0.0;}
+
+        forcecoul = qqrd2e*qtmp*q[j]*r*depsdr;
+        fvc = forcecoul*Tap - Vc*dTap/r;
+        fpair = factor_coul*fvc;
+
+        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 (Tapflag) ecoul = Vc*Tap;
+	  else ecoul = Vc - offset[itype][jtype];
+          ecoul *= factor_coul;
+        }
+
+        if (evflag) ev_tally(i,j,nlocal,newton_pair,0.0,
+                             ecoul,fpair,delx,dely,delz);
+      }
+    }
+  }
+
+  if (vflag_fdotr) virial_fdotr_compute();
+}
+
+/* ----------------------------------------------------------------------
+   allocate all arrays
+------------------------------------------------------------------------- */
+
+void PairCoulshield::allocate()
+{
+  allocated = 1;
+  int n = atom->ntypes;
+
+  memory->create(setflag,n+1,n+1,"pair:setflag");
+  for (int i = 1; i <= n; i++)
+    for (int j = i; j <= n; j++)
+      setflag[i][j] = 0;
+
+  memory->create(cutsq,n+1,n+1,"pair:cutsq");
+  memory->create(cut,n+1,n+1,"pair:cut");
+  memory->create(sigmae,n+1,n+1,"pair:sigmae");
+  //memory->create(rme,n+1,n+1,"pair:rme");
+  memory->create(offset,n+1,n+1,"pair:offset");
+}
+
+/* ----------------------------------------------------------------------
+   global settings
+------------------------------------------------------------------------- */
+
+void PairCoulshield::settings(int narg, char **arg)
+{
+  if (narg < 1 || narg > 2) error->all(FLERR,"Illegal pair_style command");
+
+  cut_global = force->numeric(FLERR,arg[0]);
+  if (narg == 2) Tapflag = force->numeric(FLERR,arg[1]);
+
+  // reset cutoffs that have been explicitly set
+
+  if (allocated) {
+    int i,j;
+    for (i = 1; i <= atom->ntypes; i++)
+      for (j = i+1; j <= atom->ntypes; j++)
+        if (setflag[i][j]) cut[i][j] = cut_global;
+  }
+}
+
+/* ----------------------------------------------------------------------
+   set coeffs for one or more type pairs
+------------------------------------------------------------------------- */
+
+void PairCoulshield::coeff(int narg, char **arg)
+{
+  if (narg < 3 || narg > 4) error->all(FLERR,"Incorrect args for pair coefficients");
+  if (!allocated) allocate();
+
+  int ilo,ihi,jlo,jhi;
+  force->bounds(FLERR,arg[0],atom->ntypes,ilo,ihi);
+  force->bounds(FLERR,arg[1],atom->ntypes,jlo,jhi);
+
+  double sigmae_one = force->numeric(FLERR,arg[2]);
+
+  double cut_one = cut_global;
+  if (narg == 4) cut_one = force->numeric(FLERR,arg[3]);
+
+  int count = 0;
+  for (int i = ilo; i <= ihi; i++) {
+    for (int j = MAX(jlo,i); j <= jhi; j++) {
+      sigmae[i][j] = sigmae_one;
+      cut[i][j] = cut_one;
+      setflag[i][j] = 1;
+      count++;
+    }
+  }
+
+  if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients");
+}
+
+
+/* ----------------------------------------------------------------------
+   init specific to this pair style
+------------------------------------------------------------------------- */
+
+void PairCoulshield::init_style()
+{
+  if (!atom->q_flag)
+    error->all(FLERR,"Pair style coul/shield requires atom attribute q");
+
+  neighbor->request(this,instance_me);
+}
+
+/* ----------------------------------------------------------------------
+   init for one type pair i,j and corresponding j,i
+------------------------------------------------------------------------- */
+
+double PairCoulshield::init_one(int i, int j)
+{
+  if (setflag[i][j] == 0) {
+    error->all(FLERR,"for pair style coul/shield, parameters need to be set explicitly for all pairs.");
+  }
+
+  double *q = atom->q;
+  double qqrd2e = force->qqrd2e;
+  double r,r3,rarg,th,epsr;
+
+  if (offset_flag) {
+     r = cut[i][j];
+     r3 = r*r*r;
+     rarg = 1.0/sigmae[i][j];
+     th = r3 + pow(rarg,3.0);
+     epsr = 1.0/pow(th,0.333333333333333333);
+     offset[i][j] = qqrd2e*q[i]*q[j]*epsr;
+  } else offset[i][j] = 0.0;
+
+
+  sigmae[j][i] = sigmae[i][j];
+  offset[j][i] = offset[i][j];
+  cut[j][i] = cut[i][j];
+
+  return cut[i][j];
+}
+
+/* ----------------------------------------------------------------------
+  proc 0 writes to restart file
+------------------------------------------------------------------------- */
+
+void PairCoulshield::write_restart(FILE *fp)
+{
+  write_restart_settings(fp);
+
+  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);
+      if (setflag[i][j]) {
+        fwrite(&sigmae[i][j],sizeof(double),1,fp);
+        fwrite(&cut[i][j],sizeof(double),1,fp);
+      }
+    }
+}
+
+/* ----------------------------------------------------------------------
+  proc 0 reads from restart file, bcasts
+------------------------------------------------------------------------- */
+
+void PairCoulshield::read_restart(FILE *fp)
+{
+  read_restart_settings(fp);
+  allocate();
+
+  int i,j;
+  int me = comm->me;
+  for (i = 1; i <= atom->ntypes; i++)
+    for (j = i; j <= atom->ntypes; j++) {
+      if (setflag[i][j]) {
+        if (me == 0) {
+          fread(&sigmae[i][j],sizeof(double),1,fp);
+          fread(&cut[i][j],sizeof(double),1,fp);
+        }
+        MPI_Bcast(&sigmae[i][j],1,MPI_DOUBLE,0,world);
+        MPI_Bcast(&cut[i][j],1,MPI_DOUBLE,0,world);
+      }
+    }
+}
+
+/* ----------------------------------------------------------------------
+  proc 0 writes to restart file
+------------------------------------------------------------------------- */
+
+void PairCoulshield::write_restart_settings(FILE *fp)
+{
+  fwrite(&cut_global,sizeof(double),1,fp);
+  fwrite(&offset_flag,sizeof(int),1,fp);
+  fwrite(&mix_flag,sizeof(int),1,fp);
+}
+
+/* ----------------------------------------------------------------------
+  proc 0 reads from restart file, bcasts
+------------------------------------------------------------------------- */
+
+void PairCoulshield::read_restart_settings(FILE *fp)
+{
+  if (comm->me == 0) {
+    fread(&cut_global,sizeof(double),1,fp);
+    fread(&offset_flag,sizeof(int),1,fp);
+    fread(&mix_flag,sizeof(int),1,fp);
+  }
+  MPI_Bcast(&cut_global,1,MPI_DOUBLE,0,world);
+  MPI_Bcast(&offset_flag,1,MPI_INT,0,world);
+  MPI_Bcast(&mix_flag,1,MPI_INT,0,world);
+}
+
+/* ---------------------------------------------------------------------- */
+
+double PairCoulshield::single(int i, int j, int itype, int jtype,
+                           double rsq, double factor_coul, double factor_lj,
+                           double &fforce)
+{
+  double r, rarg,Vc,fvc,forcecoul,phishieldec;
+  double r3,th,epsr,depsdr,Tap,dTap;
+  double *q = atom->q;
+  double qqrd2e = force->qqrd2e;
+
+   r = sqrt(rsq);
+   r3 = rsq*r;
+   rarg = 1.0/sigmae[itype][jtype];
+   th = r3 + pow(rarg,3.0);
+   epsr = 1.0/pow(th,0.333333333333333333);
+   depsdr = pow(epsr,4.0);
+   Vc = qqrd2e*q[i]*q[j]*epsr;
+
+   // turn on/off Tapper function
+   if (Tapflag) {
+     Tap = calc_Tap(r,sqrt(cutsq[itype][jtype]));
+     dTap = calc_dTap(r,sqrt(cutsq[itype][jtype]));
+   } else {Tap = 1.0; dTap = 0.0;}
+
+   forcecoul = qqrd2e*q[i]*q[j]*r*depsdr;
+   fvc = forcecoul*Tap - Vc*dTap/r;
+   fforce = factor_coul*fvc;
+
+  if (Tapflag) phishieldec = factor_coul*Vc*Tap;
+  else phishieldec = Vc - offset[itype][jtype];
+  return factor_coul*phishieldec;
+}

src/USER-MISC/pair_coul_shield.h

+/* -*- 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(coul/shield,PairCoulshield)
+
+#else
+
+#ifndef LMP_PAIR_COUL_SHIELD_H
+#define LMP_PAIR_COUL_SHIELD_H
+
+#include "pair.h"
+
+namespace LAMMPS_NS {
+
+class PairCoulshield : public Pair {
+ public:
+  PairCoulshield(class LAMMPS *);
+  virtual ~PairCoulshield();
+
+  virtual void compute(int, int);
+
+  virtual void settings(int, char **);
+  virtual void coeff(int, char **);
+  virtual void init_style();
+  virtual double init_one(int, int);
+  virtual void write_restart(FILE *);
+  virtual void read_restart(FILE *);
+  virtual void write_restart_settings(FILE *);
+  virtual void read_restart_settings(FILE *);
+
+  virtual double single(int, int, int, int, double, double, double, double &);
+
+ protected:
+  double cut_global;
+  double **cut;
+  double **sigmae, **offset;
+  //double a_eps, b_eps, eps_s;
+  int Tapflag;
+
+  void allocate();
+
+/* ----Calculate the long-range cutoff term */
+  inline double calc_Tap(double r_ij, double Rcut) {
+    double Tap,r;
+    //int Tap_coeff[8] = {1,0,0,0,-35,84,-70,20};
+    double Tap_coeff[8] = {1.0,0.0,0.0,0.0,-35.0,84.0,-70.0,20.0};
+
+    r = r_ij/Rcut;
+    Tap = 0.0;
+
+    Tap = Tap_coeff[7] * r + Tap_coeff[6];
+    Tap = Tap * r  + Tap_coeff[5];
+    Tap = Tap * r  + Tap_coeff[4];
+    Tap = Tap * r  + Tap_coeff[3];
+    Tap = Tap * r  + Tap_coeff[2];
+    Tap = Tap * r  + Tap_coeff[1];
+    Tap = Tap * r  + Tap_coeff[0];
+
+    return(Tap);
+  }
+
+ /* ----Calculate the derivatives of long-range cutoff term */
+  inline double calc_dTap(double r_ij, double Rcut) {
+    double dTap,r;
+    double Tap_coeff[8] = {1.0,0.0,0.0,0.0,-35.0,84.0,-70.0,20.0};
+
+    r = r_ij/Rcut;
+    dTap = 0.0;
+
+    dTap = 7.0*Tap_coeff[7] * r + 6.0*Tap_coeff[6];
+    dTap = dTap * r  + 5.0*Tap_coeff[5];
+    dTap = dTap * r  + 4.0*Tap_coeff[4];
+    dTap = dTap * r  + 3.0*Tap_coeff[3];
+    dTap = dTap * r  + 2.0*Tap_coeff[2];
+    dTap = dTap * r  + Tap_coeff[1];
+    dTap = dTap/Rcut;
+
+    return(dTap);
+  }
+};
+
+}
+
+#endif
+#endif

src/USER-MISC/pair_ilp_gr_hBN.h

+/* ----------------------------------------------------------------------
+   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(ILP/graphene/hBN,PairILP)
+
+#else
+
+#ifndef LMP_PAIR_ILP_H
+#define LMP_PAIR_ILP_H
+
+#include "pair.h"
+#include "my_page.h"
+#include <math.h>
+
+namespace LAMMPS_NS {
+
+class PairILP : public Pair {
+ public:
+  PairILP(class LAMMPS *);
+  virtual ~PairILP();
+
+  virtual void compute(int, int);
+  void settings(int, char **);
+  void coeff(int, char **);
+  double init_one(int, int);
+  void init_style();
+  void calc_normal();
+  int pack_forward_comm(int, int *, double *, int, int *);
+  void unpack_forward_comm(int, int, double *);
+  double single(int, int, int, int, double, double, double, double &);
+
+ protected:
+  int me;
+  int maxlocal;                    // size of numneigh, firstneigh arrays
+  int pgsize;                      // size of neighbor page
+  int oneatom;                     // max # of neighbors for one atom
+  MyPage<int> *ipage;              // neighbor list pages
+  int *ILP_numneigh;                // # of pair neighbors for each atom
+  int **ILP_firstneigh;             // ptr to 1st neighbor of each atom
+  int tap_flag;			   // flag to turn on/off Tapper function
+
+  struct Param {
+    double z0,alpha,epsilon,C,delta,d,sR,reff,C6,S;
+    double delta2inv,seff,lambda,rcut;
+    int ielement,jelement;
+  };
+  Param *params;       // parameter set for I-J interactions
+  char **elements;     // names of unique elements
+  int **elem2param;    // mapping from element pairs to parameters
+  int *map;            // mapping from atom types to elements
+  int nelements;       // # of unique elements
+  int nparams;         // # of stored parameter sets
+  int maxparam;        // max # of parameter sets
+  int nmax;            // max # of atoms
+
+  double cut_global;
+  double cut_normal;
+  double **cut;
+  double **cutILPsq;    // mapping the cutoff from element pairs to parameters
+  double **offset;
+  double **normal;
+  double ***dnormdri;
+  double ****dnormal;
+  
+  void read_file( char * );
+  void allocate();
+  void ILP_neigh();
+  
+  /* ----Calculate the long-range cutoff term */
+  inline double calc_Tap(double r_ij, double Rcut) {
+    double Tap,r;
+    double Tap_coeff[8] = {1.0,0.0,0.0,0.0,-35.0,84.0,-70.0,20.0};
+  
+    r = r_ij/Rcut; 
+    if(r >= 1.0) {Tap = 0.0;}
+    else {
+      Tap = Tap_coeff[7] * r + Tap_coeff[6];
+      Tap = Tap * r  + Tap_coeff[5];
+      Tap = Tap * r  + Tap_coeff[4];
+      Tap = Tap * r  + Tap_coeff[3];
+      Tap = Tap * r  + Tap_coeff[2];
+      Tap = Tap * r  + Tap_coeff[1];
+      Tap = Tap * r  + Tap_coeff[0];
+    }
+    return(Tap);
+}
+
+  /* ----Calculate the derivatives of long-range cutoff term */
+  inline double calc_dTap(double r_ij, double Rcut) {
+    double dTap,r;
+    double Tap_coeff[8] = {1.0,0.0,0.0,0.0,-35.0,84.0,-70.0,20.0};
+
+    r = r_ij/Rcut; 
+    if(r >= 1.0) {dTap = 0.0;}
+    else {
+      dTap = 7.0*Tap_coeff[7] * r + 6.0*Tap_coeff[6];
+      dTap = dTap * r  + 5.0*Tap_coeff[5];
+      dTap = dTap * r  + 4.0*Tap_coeff[4];
+      dTap = dTap * r  + 3.0*Tap_coeff[3];
+      dTap = dTap * r  + 2.0*Tap_coeff[2];
+      dTap = dTap * r  + Tap_coeff[1];
+      dTap = dTap/Rcut;
+    }
+    return(dTap);
+  }
+};
+
+}
+
+#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.
+
+*/
+

src/USER-MISC/pair_kolmogorov_crespi_full.h

+/* ----------------------------------------------------------------------
+   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(kolmogorov/crespi/full,PairKolmogorovCrespiFull)
+
+#else
+
+#ifndef LMP_PAIR_KolmogorovCrespi_FULL_H
+#define LMP_PAIR_KolmogorovCrespi_FULL_H
+
+#include "pair.h"
+#include "my_page.h"
+#include <math.h>
+
+namespace LAMMPS_NS {
+
+class PairKolmogorovCrespiFull : public Pair {
+ public:
+  PairKolmogorovCrespiFull(class LAMMPS *);
+  virtual ~PairKolmogorovCrespiFull();
+
+  virtual void compute(int, int);
+  void settings(int, char **);
+  void coeff(int, char **);
+  double init_one(int, int);
+  void init_style();
+  void calc_normal();
+  int pack_forward_comm(int, int *, double *, int, int *);
+  void unpack_forward_comm(int, int, double *);
+  double single(int, int, int, int, double, double, double, double &);
+
+ protected:
+  int me;
+  int maxlocal;                    // size of numneigh, firstneigh arrays
+  int pgsize;                      // size of neighbor page
+  int oneatom;                     // max # of neighbors for one atom
+  MyPage<int> *ipage;              // neighbor list pages
+  int *KC_numneigh;                // # of pair neighbors for each atom
+  int **KC_firstneigh;             // ptr to 1st neighbor of each atom
+  int tap_flag;			   // flag to turn on/off Tapper function
+
+
+  struct Param {
+    double z0,C0,C2,C4,C,delta,lambda,A,S;
+    double delta2inv,z06,rcut;
+    int ielement,jelement;
+  };
+  Param *params;       // parameter set for I-J interactions
+  char **elements;     // names of unique elements
+  int **elem2param;    // mapping from element pairs to parameters
+  int *map;            // mapping from atom types to elements
+  int nelements;       // # of unique elements
+  int nparams;         // # of stored parameter sets
+  int maxparam;        // max # of parameter sets
+  int nmax;            // max # of atoms
+
+  double cut_global;
+  double cut_normal;
+  double **cut;
+  double **cutKCsq;
+  double **offset;
+  double **normal;
+  double ***dnormdri;
+  double ****dnormal;
+  
+  void read_file( char * );
+  void allocate();
+  void KC_neigh();
+
+  
+  /* ----Calculate the long-range cutoff term */
+  inline double calc_Tap(double r_ij, double Rcut) {
+    double Tap,r;
+    double Tap_coeff[8] = {1.0,0.0,0.0,0.0,-35.0,84.0,-70.0,20.0};
+  
+    r = r_ij/Rcut; 
+    if(r >= 1.0) {Tap = 0.0;}
+    else{
+      Tap = Tap_coeff[7] * r + Tap_coeff[6];
+      Tap = Tap * r  + Tap_coeff[5];
+      Tap = Tap * r  + Tap_coeff[4];
+      Tap = Tap * r  + Tap_coeff[3];
+      Tap = Tap * r  + Tap_coeff[2];
+      Tap = Tap * r  + Tap_coeff[1];
+      Tap = Tap * r  + Tap_coeff[0];
+    }
+
+    return(Tap);
+  }
+
+  /* ----Calculate the derivatives of long-range cutoff term */
+  inline double calc_dTap(double r_ij, double Rcut) {
+    double dTap,r;
+    double Tap_coeff[8] = {1.0,0.0,0.0,0.0,-35.0,84.0,-70.0,20.0};
+
+    r = r_ij/Rcut; 
+    if(r >= 1.0) {dTap = 0.0;}
+    else {
+      dTap = 7.0*Tap_coeff[7] * r + 6.0*Tap_coeff[6];
+      dTap = dTap * r  + 5.0*Tap_coeff[5];
+      dTap = dTap * r  + 4.0*Tap_coeff[4];
+      dTap = dTap * r  + 3.0*Tap_coeff[3];
+      dTap = dTap * r  + 2.0*Tap_coeff[2];
+      dTap = dTap * r  + Tap_coeff[1];
+      dTap = dTap/Rcut;
+    }
+
+    return(dTap);
+  }
+};
+
+}
+
+#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.
+
+*/
+