diff --git a/doc/src/Eqs/pair_coul_shield.jpg b/doc/src/Eqs/pair_coul_shield.jpg
new file mode 100644
index 0000000000000000000000000000000000000000..83eb059ac800f604eeb2f8628735fa468914c946
Binary files /dev/null and b/doc/src/Eqs/pair_coul_shield.jpg differ
diff --git a/doc/src/Eqs/pair_coul_shield.tex b/doc/src/Eqs/pair_coul_shield.tex
new file mode 100644
index 0000000000000000000000000000000000000000..bd92fb680165fcc59eaf88ce8faead7171e79f91
--- /dev/null
+++ b/doc/src/Eqs/pair_coul_shield.tex
@@ -0,0 +1,33 @@
+\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}
diff --git a/doc/src/Eqs/pair_ilp_gr_hBN.jpg b/doc/src/Eqs/pair_ilp_gr_hBN.jpg
new file mode 100644
index 0000000000000000000000000000000000000000..b28c0c6dc92c3fb7316ac941eba7bb1622495334
Binary files /dev/null and b/doc/src/Eqs/pair_ilp_gr_hBN.jpg differ
diff --git a/doc/src/Eqs/pair_ilp_gr_hBN.tex b/doc/src/Eqs/pair_ilp_gr_hBN.tex
new file mode 100644
index 0000000000000000000000000000000000000000..47ce30ac7c9f130a4c4c777346690c0041b5a9fc
--- /dev/null
+++ b/doc/src/Eqs/pair_ilp_gr_hBN.tex
@@ -0,0 +1,42 @@
+\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}
+%
diff --git a/doc/src/Eqs/pair_kolmogorov_crespi_full.jpg b/doc/src/Eqs/pair_kolmogorov_crespi_full.jpg
new file mode 100644
index 0000000000000000000000000000000000000000..e9d306d5050faa7957b9a827b2546d03e0efd10b
Binary files /dev/null and b/doc/src/Eqs/pair_kolmogorov_crespi_full.jpg differ
diff --git a/doc/src/Eqs/pair_kolmogorov_crespi_full.tex b/doc/src/Eqs/pair_kolmogorov_crespi_full.tex
new file mode 100644
index 0000000000000000000000000000000000000000..606aa8b5fe1142464ce290f2f38268961839ee1f
--- /dev/null
+++ b/doc/src/Eqs/pair_kolmogorov_crespi_full.tex
@@ -0,0 +1,33 @@
+\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}
diff --git a/doc/src/pair_coul_shield.txt b/doc/src/pair_coul_shield.txt
new file mode 100644
index 0000000000000000000000000000000000000000..80a79cd383265400fe76517fea0a1fa2b1168461
--- /dev/null
+++ b/doc/src/pair_coul_shield.txt
@@ -0,0 +1,83 @@
+"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).
diff --git a/doc/src/pair_ilp_gr_hBN.txt b/doc/src/pair_ilp_gr_hBN.txt
new file mode 100644
index 0000000000000000000000000000000000000000..9f0f0c5d17b62370095897cd324262f529feddef
--- /dev/null
+++ b/doc/src/pair_ilp_gr_hBN.txt
@@ -0,0 +1,106 @@
+"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)
diff --git a/doc/src/pair_kolmogorov_crespi_full.txt b/doc/src/pair_kolmogorov_crespi_full.txt
new file mode 100644
index 0000000000000000000000000000000000000000..f921fc043a6fe8ef8fe21a8f83c24da67ea2c8f9
--- /dev/null
+++ b/doc/src/pair_kolmogorov_crespi_full.txt
@@ -0,0 +1,90 @@
+"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)
diff --git a/src/USER-MISC/pair_coul_shield.cpp b/src/USER-MISC/pair_coul_shield.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..cf75b6c6c828e69ac5d1ac3b4b01abf79ae3ddd6
--- /dev/null
+++ b/src/USER-MISC/pair_coul_shield.cpp
@@ -0,0 +1,370 @@
+/* ----------------------------------------------------------------------
+ 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;
+}
diff --git a/src/USER-MISC/pair_coul_shield.h b/src/USER-MISC/pair_coul_shield.h
new file mode 100644
index 0000000000000000000000000000000000000000..1d2b7b526a937ca0a0cbddbddbe50ea61f81b464
--- /dev/null
+++ b/src/USER-MISC/pair_coul_shield.h
@@ -0,0 +1,97 @@
+/* -*- 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
diff --git a/src/USER-MISC/pair_ilp_gr_hBN.cpp b/src/USER-MISC/pair_ilp_gr_hBN.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..c60262d1d5ceb656ad0e5e0fc9e23d52d4786a8c
--- /dev/null
+++ b/src/USER-MISC/pair_ilp_gr_hBN.cpp
@@ -0,0 +1,1117 @@
+/* ----------------------------------------------------------------------
+ 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 full version of the potential described in
+ [Maaravi et al, J. Phys. Chem. C 121, 22826-22835 (2017)]
+ The definition of normals are the same as that in
+ [Kolmogorov & Crespi, Phys. Rev. B 71, 235415 (2005)]
+------------------------------------------------------------------------- */
+
+#include <math.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <mpi.h>
+#include "pair_ilp_gr_hBN.h"
+#include "atom.h"
+#include "comm.h"
+#include "force.h"
+#include "neighbor.h"
+#include "neigh_list.h"
+#include "neigh_request.h"
+#include "my_page.h"
+#include "memory.h"
+#include "error.h"
+
+using namespace LAMMPS_NS;
+
+#define MAXLINE 1024
+#define DELTA 4
+#define PGDELTA 1
+
+/* ---------------------------------------------------------------------- */
+
+PairILP::PairILP(LAMMPS *lmp) : Pair(lmp)
+{
+ writedata = 1;
+
+ // initialize element to parameter maps
+ nelements = 0;
+ elements = NULL;
+ nparams = maxparam = 0;
+ params = NULL;
+ elem2param = NULL;
+ cutILPsq = NULL;
+ map = NULL;
+
+ nmax = 0;
+ maxlocal = 0;
+ ILP_numneigh = NULL;
+ ILP_firstneigh = NULL;
+ ipage = NULL;
+ pgsize = oneatom = 0;
+
+ normal = NULL;
+ dnormal = NULL;
+ dnormdri = NULL;
+
+ // always compute energy offset
+ offset_flag = 1;
+
+ // set comm size needed by this Pair
+ comm_forward = 39;
+ tap_flag = 1;
+}
+
+/* ---------------------------------------------------------------------- */
+
+PairILP::~PairILP()
+{
+ memory->destroy(ILP_numneigh);
+ memory->sfree(ILP_firstneigh);
+ delete [] ipage;
+ memory->destroy(normal);
+ memory->destroy(dnormal);
+ memory->destroy(dnormdri);
+
+ if (allocated) {
+ memory->destroy(setflag);
+ memory->destroy(cutsq);
+ memory->destroy(cut);
+ memory->destroy(offset);
+ }
+
+ if (elements)
+ for (int i = 0; i < nelements; i++) delete [] elements[i];
+ delete [] elements;
+ memory->destroy(params);
+ memory->destroy(elem2param);
+ memory->destroy(cutILPsq);
+ if (allocated) delete [] map;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void PairILP::compute(int eflag, int vflag)
+{
+ int i,j,ii,jj,inum,jnum,itype,jtype,id,k,l,kk,ll;
+ tagint itag,jtag;
+ double prodnorm1,prodnorm2,fkcx,fkcy,fkcz;
+ double xtmp,ytmp,ztmp,delx,dely,delz,evdwl,fpair,fpair1,fpair2;
+ double rsq,r,Rcut,rhosq1,rhosq2,exp0,exp1,exp2,r2inv,r6inv,r8inv,Tap,dTap,Vkc;
+ double frho1,frho2,TSvdw,TSvdw2inv,Erep,fsum,rdsq1,rdsq2;
+ int *ilist,*jlist,*numneigh,**firstneigh;
+ int *ILP_neighs_i,*ILP_neighs_j;
+
+ evdwl = 0.0;
+ if (eflag || vflag) ev_setup(eflag,vflag);
+ else evflag = vflag_fdotr = 0;
+
+ double **x = atom->x;
+ double **f = atom->f;
+ int *type = atom->type;
+ tagint *tag = atom->tag;
+ int nlocal = atom->nlocal;
+ int newton_pair = force->newton_pair;
+ double dprodnorm1[3] = {0.0, 0.0, 0.0};
+ double dprodnorm2[3] = {0.0, 0.0, 0.0};
+ double fp1[3] = {0.0, 0.0, 0.0};
+ double fp2[3] = {0.0, 0.0, 0.0};
+ double fprod1[3] = {0.0, 0.0, 0.0};
+ double fprod2[3] = {0.0, 0.0, 0.0};
+
+ inum = list->inum;
+ ilist = list->ilist;
+ numneigh = list->numneigh;
+ firstneigh = list->firstneigh;
+ // Build full neighbor list
+ ILP_neigh();
+ // Calculate the normals
+ calc_normal();
+
+ // communicate the normal vector and its derivatives
+ comm->forward_comm_pair(this);
+
+ // loop over neighbors of my atoms
+ for (ii = 0; ii < inum; ii++) {
+ i = ilist[ii];
+ itag = tag[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];
+ j &= NEIGHMASK;
+ jtype = type[j];
+ jtag = tag[j];
+
+ // two-body interactions from full neighbor list, skip half of them
+ if (itag > jtag) {
+ if ((itag+jtag) % 2 == 0) continue;
+ } else if (itag < jtag) {
+ if ((itag+jtag) % 2 == 1) continue;
+ } else {
+ if (x[j][2] < ztmp) continue;
+ if (x[j][2] == ztmp && x[j][1] < ytmp) continue;
+ if (x[j][2] == ztmp && x[j][1] == ytmp && x[j][0] < xtmp) continue;
+ }
+
+ delx = xtmp - x[j][0];
+ dely = ytmp - x[j][1];
+ delz = ztmp - x[j][2];
+ rsq = delx*delx + dely*dely + delz*delz;
+
+ // only include the interation between different layers
+ if (rsq < cutsq[itype][jtype] && atom->molecule[i] != atom->molecule[j]) {
+
+ int iparam_ij = elem2param[map[itype]][map[jtype]];
+ Param& p = params[iparam_ij];
+
+ r = sqrt(rsq);
+ r2inv = 1.0/rsq;
+ r6inv = r2inv*r2inv*r2inv;
+ r8inv = r6inv*r2inv;
+ // turn on/off Tapper function
+ if (tap_flag) {
+ Rcut = sqrt(cutsq[itype][jtype]);
+ Tap = calc_Tap(r,Rcut);
+ dTap = calc_dTap(r,Rcut);
+ } else {Tap = 1.0; dTap = 0.0;}
+
+ // Calculate the transverse distance
+ // note that rho_ij does not equal to rho_ji except when normals are all along z
+ prodnorm1 = normal[i][0]*delx + normal[i][1]*dely + normal[i][2]*delz;
+ prodnorm2 = normal[j][0]*delx + normal[j][1]*dely + normal[j][2]*delz;
+ rhosq1 = rsq - prodnorm1*prodnorm1; // rho_ij
+ rhosq2 = rsq - prodnorm2*prodnorm2; // rho_ji
+ rdsq1 = rhosq1*p.delta2inv; // (rho_ij/delta)^2
+ rdsq2 = rhosq2*p.delta2inv; // (rho_ji/delta)^2
+
+ // store exponents
+ exp0 = exp(-p.lambda*(r-p.z0));
+ exp1 = exp(-rdsq1);
+ exp2 = exp(-rdsq2);
+
+ TSvdw = 1.0 + exp(-p.d*(r/p.seff - 1.0));
+ TSvdw2inv = 1.0/pow(TSvdw,2.0);
+ frho1 = exp1*p.C;
+ frho2 = exp2*p.C;
+ Erep = p.epsilon + frho1 + frho2;
+ Vkc = -p.C6*r6inv/TSvdw + exp0*Erep;
+
+ // derivatives
+ fpair = -6.0*p.C6*r8inv/TSvdw + p.d/p.seff*p.C6*(TSvdw-1.0)*TSvdw2inv*r8inv*r + p.lambda*exp0/r*Erep;
+ fpair1 = 2.0*exp0*frho1*p.delta2inv;
+ fpair2 = 2.0*exp0*frho2*p.delta2inv;
+ fsum = fpair + fpair1 + fpair2;
+ // derivatives of the product of rij and ni, the result is a vector
+ dprodnorm1[0] = dnormdri[0][0][i]*delx + dnormdri[1][0][i]*dely + dnormdri[2][0][i]*delz;
+ dprodnorm1[1] = dnormdri[0][1][i]*delx + dnormdri[1][1][i]*dely + dnormdri[2][1][i]*delz;
+ dprodnorm1[2] = dnormdri[0][2][i]*delx + dnormdri[1][2][i]*dely + dnormdri[2][2][i]*delz;
+ // derivatives of the product of rji and nj, the result is a vector
+ dprodnorm2[0] = dnormdri[0][0][j]*delx + dnormdri[1][0][j]*dely + dnormdri[2][0][j]*delz;
+ dprodnorm2[1] = dnormdri[0][1][j]*delx + dnormdri[1][1][j]*dely + dnormdri[2][1][j]*delz;
+ dprodnorm2[2] = dnormdri[0][2][j]*delx + dnormdri[1][2][j]*dely + dnormdri[2][2][j]*delz;
+ fp1[0] = prodnorm1*normal[i][0]*fpair1;
+ fp1[1] = prodnorm1*normal[i][1]*fpair1;
+ fp1[2] = prodnorm1*normal[i][2]*fpair1;
+ fp2[0] = prodnorm2*normal[j][0]*fpair2;
+ fp2[1] = prodnorm2*normal[j][1]*fpair2;
+ fp2[2] = prodnorm2*normal[j][2]*fpair2;
+ fprod1[0] = prodnorm1*dprodnorm1[0]*fpair1;
+ fprod1[1] = prodnorm1*dprodnorm1[1]*fpair1;
+ fprod1[2] = prodnorm1*dprodnorm1[2]*fpair1;
+ fprod2[0] = prodnorm2*dprodnorm2[0]*fpair2;
+ fprod2[1] = prodnorm2*dprodnorm2[1]*fpair2;
+ fprod2[2] = prodnorm2*dprodnorm2[2]*fpair2;
+ fkcx = (delx*fsum - fp1[0] - fp2[0])*Tap - Vkc*dTap*delx/r;
+ fkcy = (dely*fsum - fp1[1] - fp2[1])*Tap - Vkc*dTap*dely/r;
+ fkcz = (delz*fsum - fp1[2] - fp2[2])*Tap - Vkc*dTap*delz/r;
+
+ f[i][0] += fkcx - fprod1[0]*Tap;
+ f[i][1] += fkcy - fprod1[1]*Tap;
+ f[i][2] += fkcz - fprod1[2]*Tap;
+ if (newton_pair || j < nlocal) {
+ f[j][0] -= fkcx + fprod2[0]*Tap;
+ f[j][1] -= fkcy + fprod2[1]*Tap;
+ f[j][2] -= fkcz + fprod2[2]*Tap;
+ }
+
+ // calculate the forces acted on the neighbors of atom i from atom j
+ ILP_neighs_i = ILP_firstneigh[i];
+ for (kk = 0; kk < ILP_numneigh[i]; kk++) {
+ k = ILP_neighs_i[kk];
+ if (k == i) continue;
+ // derivatives of the product of rij and ni respect to rk, k=0,1,2, where atom k is the neighbors of atom i
+ dprodnorm1[0] = dnormal[0][0][kk][i]*delx + dnormal[1][0][kk][i]*dely + dnormal[2][0][kk][i]*delz;
+ dprodnorm1[1] = dnormal[0][1][kk][i]*delx + dnormal[1][1][kk][i]*dely + dnormal[2][1][kk][i]*delz;
+ dprodnorm1[2] = dnormal[0][2][kk][i]*delx + dnormal[1][2][kk][i]*dely + dnormal[2][2][kk][i]*delz;
+ f[k][0] += (-prodnorm1*dprodnorm1[0]*fpair1)*Tap;
+ f[k][1] += (-prodnorm1*dprodnorm1[1]*fpair1)*Tap;
+ f[k][2] += (-prodnorm1*dprodnorm1[2]*fpair1)*Tap;
+ }
+
+ // calculate the forces acted on the neighbors of atom j from atom i
+ ILP_neighs_j = ILP_firstneigh[j];
+ for (ll = 0; ll < ILP_numneigh[j]; ll++) {
+ l = ILP_neighs_j[ll];
+ if (l == j) continue;
+ if (newton_pair || l < nlocal) {
+ // derivatives of the product of rji and nj respect to rl, l=0,1,2, where atom l is the neighbors of atom j
+ dprodnorm2[0] = dnormal[0][0][ll][j]*delx + dnormal[1][0][ll][j]*dely + dnormal[2][0][ll][j]*delz;
+ dprodnorm2[1] = dnormal[0][1][ll][j]*delx + dnormal[1][1][ll][j]*dely + dnormal[2][1][ll][j]*delz;
+ dprodnorm2[2] = dnormal[0][2][ll][j]*delx + dnormal[1][2][ll][j]*dely + dnormal[2][2][ll][j]*delz;
+ f[l][0] += (-prodnorm2*dprodnorm2[0]*fpair2)*Tap;
+ f[l][1] += (-prodnorm2*dprodnorm2[1]*fpair2)*Tap;
+ f[l][2] += (-prodnorm2*dprodnorm2[2]*fpair2)*Tap;
+ }
+ }
+
+ if (eflag) {
+ if (tap_flag) evdwl = Tap*Vkc;
+ else evdwl = Vkc - offset[itype][jtype];
+ }
+
+ if (evflag){
+ ev_tally_xyz(i,j,nlocal,newton_pair,evdwl,0,
+ fkcx,fkcy,fkcz,delx,dely,delz);
+ }
+ }
+ }
+ }
+ if (vflag_fdotr) virial_fdotr_compute();
+}
+
+/* ----------------------------------------------------------------------
+ Calculate the normals for each atom
+------------------------------------------------------------------------- */
+void PairILP::calc_normal()
+{
+ int i,j,ii,jj,inum,jnum,itype,jtype;
+ int cont,k,kk,cont1,id,ip,m;
+ tagint itag,jtag;
+ double dist,nn,xtp,ytp,ztp,delx,dely,delz,nn2;
+ int *ilist,*jlist;
+ double pv12[3],pv31[3],pv23[3],n1[3],dni[3],dnn[3][3],vet[3][3],dpvdri[3][3];
+ double dn1[3][3][3],dpv12[3][3][3],dpv23[3][3][3],dpv31[3][3][3];
+
+ double **x = atom->x;
+ int *type = atom->type;
+ int ntypes = atom->ntypes;
+ int nlocal = atom->nlocal;
+ tagint *tag = atom->tag;
+
+ // grow normal array if necessary
+
+ if (atom->nmax > nmax) {
+ memory->destroy(normal);
+ memory->destroy(dnormal);
+ memory->destroy(dnormdri);
+ nmax = atom->nmax;
+ memory->create(normal,nmax,3,"ILP:normal");
+ memory->create(dnormdri,3,3,nmax,"ILP:dnormdri");
+ memory->create(dnormal,3,3,3,nmax,"ILP:dnormal");
+ }
+
+ inum = list->inum;
+ ilist = list->ilist;
+ //Calculate normals
+ for (ii = 0; ii < inum; ii++) {
+ i = ilist[ii];
+ xtp = x[i][0];
+ ytp = x[i][1];
+ ztp = x[i][2];
+ itype = type[i];
+
+ // Initialize the arrays
+ for (id = 0; id < 3; id++){
+ pv12[id] = 0.0;
+ pv31[id] = 0.0;
+ pv23[id] = 0.0;
+ n1[id] = 0.0;
+ dni[id] = 0.0;
+ normal[i][id] = 0.0;
+ for (ip = 0; ip < 3; ip++){
+ vet[ip][id] = 0.0;
+ dnn[ip][id] = 0.0;
+ dpvdri[ip][id] = 0.0;
+ dnormdri[ip][id][i] = 0.0;
+ for (m = 0; m < 3; m++){
+ dpv12[ip][id][m] = 0.0;
+ dpv31[ip][id][m] = 0.0;
+ dpv23[ip][id][m] = 0.0;
+ dn1[ip][id][m] = 0.0;
+ dnormal[ip][id][m][i] = 0.0;
+ }
+ }
+ }
+
+ cont = 0;
+ jlist = ILP_firstneigh[i];
+ jnum = ILP_numneigh[i];
+ for (jj = 0; jj < jnum; jj++) {
+ j = jlist[jj];
+ j &= NEIGHMASK;
+
+ delx = x[j][0] - xtp;
+ dely = x[j][1] - ytp;
+ delz = x[j][2] - ztp;
+ vet[cont][0] = delx;
+ vet[cont][1] = dely;
+ vet[cont][2] = delz;
+ cont++;
+ }
+
+ if (cont <= 1) {
+ normal[i][0] = 0.0;
+ normal[i][1] = 0.0;
+ normal[i][2] = 1.0;
+ for (id = 0; id < 3; id++){
+ for (ip = 0; ip < 3; ip++){
+ dnormdri[id][ip][i] = 0.0;
+ for (m = 0; m < 3; m++){
+ dnormal[id][ip][m][i] = 0.0;
+ }
+ }
+ }
+ }
+ else if (cont == 2) {
+ pv12[0] = vet[0][1]*vet[1][2] - vet[1][1]*vet[0][2];
+ pv12[1] = vet[0][2]*vet[1][0] - vet[1][2]*vet[0][0];
+ pv12[2] = vet[0][0]*vet[1][1] - vet[1][0]*vet[0][1];
+ // derivatives of pv12[0] to ri
+ dpvdri[0][0] = 0.0;
+ dpvdri[0][1] = vet[0][2]-vet[1][2];
+ dpvdri[0][2] = vet[1][1]-vet[0][1];
+ // derivatives of pv12[1] to ri
+ dpvdri[1][0] = vet[1][2]-vet[0][2];
+ dpvdri[1][1] = 0.0;
+ dpvdri[1][2] = vet[0][0]-vet[1][0];
+ // derivatives of pv12[2] to ri
+ dpvdri[2][0] = vet[0][1]-vet[1][1];
+ dpvdri[2][1] = vet[1][0]-vet[0][0];
+ dpvdri[2][2] = 0.0;
+
+ dpv12[0][0][0] = 0.0;
+ dpv12[0][1][0] = vet[1][2];
+ dpv12[0][2][0] = -vet[1][1];
+ dpv12[1][0][0] = -vet[1][2];
+ dpv12[1][1][0] = 0.0;
+ dpv12[1][2][0] = vet[1][0];
+ dpv12[2][0][0] = vet[1][1];
+ dpv12[2][1][0] = -vet[1][0];
+ dpv12[2][2][0] = 0.0;
+
+ // derivatives respect to the second neighbor, atom l
+ dpv12[0][0][1] = 0.0;
+ dpv12[0][1][1] = -vet[0][2];
+ dpv12[0][2][1] = vet[0][1];
+ dpv12[1][0][1] = vet[0][2];
+ dpv12[1][1][1] = 0.0;
+ dpv12[1][2][1] = -vet[0][0];
+ dpv12[2][0][1] = -vet[0][1];
+ dpv12[2][1][1] = vet[0][0];
+ dpv12[2][2][1] = 0.0;
+
+ // derivatives respect to the third neighbor, atom n
+ // derivatives of pv12 to rn is zero
+ for (id = 0; id < 3; id++){
+ for (ip = 0; ip < 3; ip++){
+ dpv12[id][ip][2] = 0.0;
+ }
+ }
+
+ n1[0] = pv12[0];
+ n1[1] = pv12[1];
+ n1[2] = pv12[2];
+ // the magnitude of the normal vector
+ nn2 = n1[0]*n1[0] + n1[1]*n1[1] + n1[2]*n1[2];
+ nn = sqrt(nn2);
+ if (nn == 0) error->all(FLERR,"The magnitude of the normal vector is zero");
+ // the unit normal vector
+ normal[i][0] = n1[0]/nn;
+ normal[i][1] = n1[1]/nn;
+ normal[i][2] = n1[2]/nn;
+ // derivatives of nn, dnn:3x1 vector
+ dni[0] = (n1[0]*dpvdri[0][0] + n1[1]*dpvdri[1][0] + n1[2]*dpvdri[2][0])/nn;
+ dni[1] = (n1[0]*dpvdri[0][1] + n1[1]*dpvdri[1][1] + n1[2]*dpvdri[2][1])/nn;
+ dni[2] = (n1[0]*dpvdri[0][2] + n1[1]*dpvdri[1][2] + n1[2]*dpvdri[2][2])/nn;
+ // derivatives of unit vector ni respect to ri, the result is 3x3 matrix
+ for (id = 0; id < 3; id++){
+ for (ip = 0; ip < 3; ip++){
+ dnormdri[id][ip][i] = dpvdri[id][ip]/nn - n1[id]*dni[ip]/nn2;
+ }
+ }
+ // derivatives of non-normalized normal vector, dn1:3x3x3 array
+ for (id = 0; id < 3; id++){
+ for (ip = 0; ip < 3; ip++){
+ for (m = 0; m < 3; m++){
+ dn1[id][ip][m] = dpv12[id][ip][m];
+ }
+ }
+ }
+ // derivatives of nn, dnn:3x3 vector
+ // dnn[id][m]: the derivative of nn respect to r[id][m], id,m=0,1,2
+ // r[id][m]: the id's component of atom m
+ for (m = 0; m < 3; m++){
+ for (id = 0; id < 3; id++){
+ dnn[id][m] = (n1[0]*dn1[0][id][m] + n1[1]*dn1[1][id][m] + n1[2]*dn1[2][id][m])/nn;
+ }
+ }
+ // dnormal[id][ip][m][i]: the derivative of normal[id] respect to r[ip][m], id,ip=0,1,2
+ // for atom m, which is a neighbor atom of atom i, m=0,jnum-1
+ for (m = 0; m < 3; m++){
+ for (id = 0; id < 3; id++){
+ for (ip = 0; ip < 3; ip++){
+ dnormal[id][ip][m][i] = dn1[id][ip][m]/nn - n1[id]*dnn[ip][m]/nn2;
+ }
+ }
+ }
+ }
+//##############################################################################################
+
+ else if(cont == 3) {
+ pv12[0] = vet[0][1]*vet[1][2] - vet[1][1]*vet[0][2];
+ pv12[1] = vet[0][2]*vet[1][0] - vet[1][2]*vet[0][0];
+ pv12[2] = vet[0][0]*vet[1][1] - vet[1][0]*vet[0][1];
+ // derivatives respect to the first neighbor, atom k
+ dpv12[0][0][0] = 0.0;
+ dpv12[0][1][0] = vet[1][2];
+ dpv12[0][2][0] = -vet[1][1];
+ dpv12[1][0][0] = -vet[1][2];
+ dpv12[1][1][0] = 0.0;
+ dpv12[1][2][0] = vet[1][0];
+ dpv12[2][0][0] = vet[1][1];
+ dpv12[2][1][0] = -vet[1][0];
+ dpv12[2][2][0] = 0.0;
+ // derivatives respect to the second neighbor, atom l
+ dpv12[0][0][1] = 0.0;
+ dpv12[0][1][1] = -vet[0][2];
+ dpv12[0][2][1] = vet[0][1];
+ dpv12[1][0][1] = vet[0][2];
+ dpv12[1][1][1] = 0.0;
+ dpv12[1][2][1] = -vet[0][0];
+ dpv12[2][0][1] = -vet[0][1];
+ dpv12[2][1][1] = vet[0][0];
+ dpv12[2][2][1] = 0.0;
+
+ // derivatives respect to the third neighbor, atom n
+ for (id = 0; id < 3; id++){
+ for (ip = 0; ip < 3; ip++){
+ dpv12[id][ip][2] = 0.0;
+ }
+ }
+
+ pv31[0] = vet[2][1]*vet[0][2] - vet[0][1]*vet[2][2];
+ pv31[1] = vet[2][2]*vet[0][0] - vet[0][2]*vet[2][0];
+ pv31[2] = vet[2][0]*vet[0][1] - vet[0][0]*vet[2][1];
+ // derivatives respect to the first neighbor, atom k
+ dpv31[0][0][0] = 0.0;
+ dpv31[0][1][0] = -vet[2][2];
+ dpv31[0][2][0] = vet[2][1];
+ dpv31[1][0][0] = vet[2][2];
+ dpv31[1][1][0] = 0.0;
+ dpv31[1][2][0] = -vet[2][0];
+ dpv31[2][0][0] = -vet[2][1];
+ dpv31[2][1][0] = vet[2][0];
+ dpv31[2][2][0] = 0.0;
+ // derivatives respect to the third neighbor, atom n
+ dpv31[0][0][2] = 0.0;
+ dpv31[0][1][2] = vet[0][2];
+ dpv31[0][2][2] = -vet[0][1];
+ dpv31[1][0][2] = -vet[0][2];
+ dpv31[1][1][2] = 0.0;
+ dpv31[1][2][2] = vet[0][0];
+ dpv31[2][0][2] = vet[0][1];
+ dpv31[2][1][2] = -vet[0][0];
+ dpv31[2][2][2] = 0.0;
+ // derivatives respect to the second neighbor, atom l
+ for (id = 0; id < 3; id++){
+ for (ip = 0; ip < 3; ip++){
+ dpv31[id][ip][1] = 0.0;
+ }
+ }
+
+ pv23[0] = vet[1][1]*vet[2][2] - vet[2][1]*vet[1][2];
+ pv23[1] = vet[1][2]*vet[2][0] - vet[2][2]*vet[1][0];
+ pv23[2] = vet[1][0]*vet[2][1] - vet[2][0]*vet[1][1];
+ // derivatives respect to the second neighbor, atom k
+ for (id = 0; id < 3; id++){
+ for (ip = 0; ip < 3; ip++){
+ dpv23[id][ip][0] = 0.0;
+ }
+ }
+ // derivatives respect to the second neighbor, atom l
+ dpv23[0][0][1] = 0.0;
+ dpv23[0][1][1] = vet[2][2];
+ dpv23[0][2][1] = -vet[2][1];
+ dpv23[1][0][1] = -vet[2][2];
+ dpv23[1][1][1] = 0.0;
+ dpv23[1][2][1] = vet[2][0];
+ dpv23[2][0][1] = vet[2][1];
+ dpv23[2][1][1] = -vet[2][0];
+ dpv23[2][2][1] = 0.0;
+ // derivatives respect to the third neighbor, atom n
+ dpv23[0][0][2] = 0.0;
+ dpv23[0][1][2] = -vet[1][2];
+ dpv23[0][2][2] = vet[1][1];
+ dpv23[1][0][2] = vet[1][2];
+ dpv23[1][1][2] = 0.0;
+ dpv23[1][2][2] = -vet[1][0];
+ dpv23[2][0][2] = -vet[1][1];
+ dpv23[2][1][2] = vet[1][0];
+ dpv23[2][2][2] = 0.0;
+
+//############################################################################################
+ // average the normal vectors by using the 3 neighboring planes
+ n1[0] = (pv12[0] + pv31[0] + pv23[0])/cont;
+ n1[1] = (pv12[1] + pv31[1] + pv23[1])/cont;
+ n1[2] = (pv12[2] + pv31[2] + pv23[2])/cont;
+ // the magnitude of the normal vector
+ nn2 = n1[0]*n1[0] + n1[1]*n1[1] + n1[2]*n1[2];
+ nn = sqrt(nn2);
+ if (nn == 0) error->all(FLERR,"The magnitude of the normal vector is zero");
+ // the unit normal vector
+ normal[i][0] = n1[0]/nn;
+ normal[i][1] = n1[1]/nn;
+ normal[i][2] = n1[2]/nn;
+
+ // for the central atoms, dnormdri is always zero
+ for (id = 0; id < 3; id++){
+ for (ip = 0; ip < 3; ip++){
+ dnormdri[id][ip][i] = 0.0;
+ }
+ }
+
+ // derivatives of non-normalized normal vector, dn1:3x3x3 array
+ for (id = 0; id < 3; id++){
+ for (ip = 0; ip < 3; ip++){
+ for (m = 0; m < 3; m++){
+ dn1[id][ip][m] = (dpv12[id][ip][m] + dpv23[id][ip][m] + dpv31[id][ip][m])/cont;
+ }
+ }
+ }
+ // derivatives of nn, dnn:3x3 vector
+ // dnn[id][m]: the derivative of nn respect to r[id][m], id,m=0,1,2
+ // r[id][m]: the id's component of atom m
+ for (m = 0; m < 3; m++){
+ for (id = 0; id < 3; id++){
+ dnn[id][m] = (n1[0]*dn1[0][id][m] + n1[1]*dn1[1][id][m] + n1[2]*dn1[2][id][m])/nn;
+ }
+ }
+ // dnormal[id][ip][m][i]: the derivative of normal[id] respect to r[ip][m], id,ip=0,1,2
+ // for atom m, which is a neighbor atom of atom i, m=0,jnum-1
+ for (m = 0; m < 3; m++){
+ for (id = 0; id < 3; id++){
+ for (ip = 0; ip < 3; ip++){
+ dnormal[id][ip][m][i] = dn1[id][ip][m]/nn - n1[id]*dnn[ip][m]/nn2;
+ }
+ }
+ }
+ }
+ else {
+ error->all(FLERR,"There are too many neighbors for calculating normals");
+ }
+
+//##############################################################################################
+ }
+}
+
+/* ----------------------------------------------------------------------
+ init specific to this pair style
+------------------------------------------------------------------------- */
+
+void PairILP::init_style()
+{
+ if (force->newton_pair == 0)
+ error->all(FLERR,"Pair style ILP requires newton pair on");
+
+ // need a full neighbor list, including neighbors of ghosts
+
+ int irequest = neighbor->request(this,instance_me);
+ neighbor->requests[irequest]->half = 0;
+ neighbor->requests[irequest]->full = 1;
+ neighbor->requests[irequest]->ghost = 1;
+
+ // local ILP neighbor list
+ // create pages if first time or if neighbor pgsize/oneatom has changed
+
+ int create = 0;
+ if (ipage == NULL) create = 1;
+ if (pgsize != neighbor->pgsize) create = 1;
+ if (oneatom != neighbor->oneatom) create = 1;
+
+ if (create) {
+ delete [] ipage;
+ pgsize = neighbor->pgsize;
+ oneatom = neighbor->oneatom;
+
+ int nmypage= comm->nthreads;
+ ipage = new MyPage<int>[nmypage];
+ for (int i = 0; i < nmypage; i++)
+ ipage[i].init(oneatom,pgsize,PGDELTA);
+ }
+}
+
+
+/* ----------------------------------------------------------------------
+ create ILP neighbor list from main neighbor list
+ ILP neighbor list stores neighbors of ghost atoms
+ ILP_numneigh for calcualting normals and
+ ILP_pair_numneigh for calculating force
+------------------------------------------------------------------------- */
+
+void PairILP::ILP_neigh()
+{
+ int i,j,ii,jj,n,np,allnum,inum,jnum,itype,jtype;
+ double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
+ int *ilist,*jlist,*numneigh,**firstneigh;
+ int *neighptr;
+
+ double **x = atom->x;
+ int *type = atom->type;
+
+ if (atom->nmax > maxlocal) {
+ maxlocal = atom->nmax;
+ memory->destroy(ILP_numneigh);
+ memory->sfree(ILP_firstneigh);
+ memory->create(ILP_numneigh,maxlocal,"ILP:numneigh");
+ ILP_firstneigh = (int **) memory->smalloc(maxlocal*sizeof(int *),"ILP:firstneigh");
+ }
+
+ inum = list->inum;
+ allnum = list->inum + list->gnum;
+ ilist = list->ilist;
+ numneigh = list->numneigh;
+ firstneigh = list->firstneigh;
+
+ // store all ILP neighs of owned and ghost atoms
+ // scan full neighbor list of I
+
+ ipage->reset();
+
+ for (ii = 0; ii < allnum; ii++) {
+ i = ilist[ii];
+
+ n = 0;
+ neighptr = ipage->vget();
+
+ xtmp = x[i][0];
+ ytmp = x[i][1];
+ ztmp = x[i][2];
+ itype = map[type[i]];
+ jlist = firstneigh[i];
+ jnum = numneigh[i];
+
+ for (jj = 0; jj < jnum; jj++) {
+ j = jlist[jj];
+ j &= NEIGHMASK;
+ jtype = map[type[j]];
+ delx = xtmp - x[j][0];
+ dely = ytmp - x[j][1];
+ delz = ztmp - x[j][2];
+ rsq = delx*delx + dely*dely + delz*delz;
+
+ if (rsq != 0 && rsq < cutILPsq[itype][jtype] && atom->molecule[i] == atom->molecule[j]) {
+ neighptr[n++] = j;
+ }
+ }
+
+ ILP_firstneigh[i] = neighptr;
+ ILP_numneigh[i] = n;
+ ipage->vgot(n);
+ if (ipage->status())
+ error->one(FLERR,"Neighbor list overflow, boost neigh_modify one");
+ }
+}
+
+
+/* ----------------------------------------------------------------------
+ allocate all arrays
+------------------------------------------------------------------------- */
+
+void PairILP::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(offset,n+1,n+1,"pair:offset");
+ map = new int[atom->ntypes+1];
+}
+
+/* ----------------------------------------------------------------------
+ global settings
+------------------------------------------------------------------------- */
+
+void PairILP::settings(int narg, char **arg)
+{
+ if (narg < 1 || narg > 2) error->all(FLERR,"Illegal pair_style command");
+ if (strcmp(force->pair_style,"hybrid/overlay")!=0)
+ error->all(FLERR,"ERROR: requires hybrid/overlay pair_style");
+
+ cut_global = force->numeric(FLERR,arg[0]);
+ if (narg == 2) tap_flag = 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; j <= atom->ntypes; j++)
+ if (setflag[i][j]) cut[i][j] = cut_global;
+ }
+}
+
+/* ----------------------------------------------------------------------
+ set coeffs for one or more type pairs
+------------------------------------------------------------------------- */
+
+void PairILP::coeff(int narg, char **arg)
+{
+ int i,j,n;
+
+ if (narg != 3 + atom->ntypes)
+ 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);
+
+ // read args that map atom types to elements in potential file
+ // map[i] = which element the Ith atom type is, -1 if NULL
+ // nelements = # of unique elements
+ // elements = list of element names
+
+ if (elements) {
+ for (i = 0; i < nelements; i++) delete [] elements[i];
+ delete [] elements;
+ }
+ elements = new char*[atom->ntypes];
+ for (i = 0; i < atom->ntypes; i++) elements[i] = NULL;
+
+ nelements = 0;
+ for (i = 3; i < narg; i++) {
+ if (strcmp(arg[i],"NULL") == 0) {
+ map[i-2] = -1;
+ continue;
+ }
+ for (j = 0; j < nelements; j++)
+ if (strcmp(arg[i],elements[j]) == 0) break;
+ map[i-2] = j;
+ if (j == nelements) {
+ n = strlen(arg[i]) + 1;
+ elements[j] = new char[n];
+ strcpy(elements[j],arg[i]);
+ nelements++;
+ }
+ }
+
+
+ read_file(arg[2]);
+
+ double cut_one = cut_global;
+
+ int count = 0;
+ for (int i = ilo; i <= ihi; i++) {
+ for (int j = MAX(jlo,i); j <= jhi; j++) {
+ cut[i][j] = cut_one;
+ setflag[i][j] = 1;
+ count++;
+ }
+ }
+
+ if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients");
+}
+
+
+/* ----------------------------------------------------------------------
+ init for one type pair i,j and corresponding j,i
+------------------------------------------------------------------------- */
+
+double PairILP::init_one(int i, int j)
+{
+ if (setflag[i][j] == 0) error->all(FLERR,"All pair coeffs are not set");
+
+ if (offset_flag && (cut[i][j] > 0.0)) {
+ int iparam_ij = elem2param[map[i]][map[j]];
+ Param& p = params[iparam_ij];
+ offset[i][j] = -p.C6*pow(1.0/cut[i][j],6)/(1.0 + exp(-p.d*(cut[i][j]/p.seff - 1.0)));
+ } else offset[i][j] = 0.0;
+ offset[j][i] = offset[i][j];
+
+ return cut[i][j];
+}
+
+/* ----------------------------------------------------------------------
+ read Interlayer potential file
+------------------------------------------------------------------------- */
+
+void PairILP::read_file(char *filename)
+{
+ int params_per_line = 13;
+ 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(filename);
+ if (fp == NULL) {
+ char str[128];
+ sprintf(str,"Cannot open ILP potential file %s",filename);
+ 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 i,j,n,m,nwords,ielement,jelement;
+ 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,"Insufficient format in ILP 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 = 1st args
+ // if these 2 args are in element list, then parse this line
+ // else skip to next line (continue)
+
+ 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;
+
+ // 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].z0 = atof(words[2]);
+ params[nparams].alpha = atof(words[3]);
+ params[nparams].delta = atof(words[4]);
+ params[nparams].epsilon = atof(words[5]);
+ params[nparams].C = atof(words[6]);
+ params[nparams].d = atof(words[7]);
+ params[nparams].sR = atof(words[8]);
+ params[nparams].reff = atof(words[9]);
+ params[nparams].C6 = atof(words[10]);
+ // S provides a convenient scaling of all energies
+ params[nparams].S = atof(words[11]);
+ params[nparams].rcut = atof(words[12]);
+
+ // energies in meV further scaled by S
+ // S = 43.3634 meV = 1 kcal/mol
+ double meV = 1e-3*params[nparams].S;
+ params[nparams].C *= meV;
+ params[nparams].C6 *= meV;
+ params[nparams].epsilon *= meV;
+
+ // precompute some quantities
+ params[nparams].delta2inv = pow(params[nparams].delta,-2.0);
+ params[nparams].lambda = params[nparams].alpha/params[nparams].z0;
+ params[nparams].seff = params[nparams].sR * params[nparams].reff;
+
+ nparams++;
+ }
+ memory->destroy(elem2param);
+ memory->destroy(cutILPsq);
+ memory->create(elem2param,nelements,nelements,"pair:elem2param");
+ memory->create(cutILPsq,nelements,nelements,"pair:cutILPsq");
+ for (i = 0; i < nelements; i++) {
+ for (j = 0; j < nelements; j++) {
+ n = -1;
+ for (m = 0; m < nparams; m++) {
+ if (i == params[m].ielement && j == params[m].jelement) {
+ if (n >= 0) error->all(FLERR,"Potential file has duplicate entry");
+ n = m;
+ }
+ }
+ if (n < 0) error->all(FLERR,"Potential file is missing an entry");
+ elem2param[i][j] = n;
+ cutILPsq[i][j] = params[n].rcut*params[n].rcut;
+ }
+ }
+ delete [] words;
+}
+
+/* ---------------------------------------------------------------------- */
+
+double PairILP::single(int i, int j, int itype, int jtype, double rsq,
+ double factor_coul, double factor_lj,
+ double &fforce)
+{
+ double r,r2inv,r6inv,r8inv,forcelj,philj,fpair;
+ double Tap,dTap,Vkc,TSvdw,TSvdw2inv;
+
+ int iparam_ij = elem2param[map[itype]][map[jtype]];
+ Param& p = params[iparam_ij];
+
+ r = sqrt(rsq);
+ // turn on/off Tapper function
+ if (tap_flag) {
+ Tap = calc_Tap(r,sqrt(cutsq[itype][jtype]));
+ dTap = calc_dTap(r,sqrt(cutsq[itype][jtype]));
+ } else {Tap = 1.0; dTap = 0.0;}
+
+ r2inv = 1.0/rsq;
+ r6inv = r2inv*r2inv*r2inv;
+ r8inv = r2inv*r6inv;
+
+ TSvdw = 1.0 + exp(-p.d*(r/p.seff - 1.0));
+ TSvdw2inv = pow(TSvdw,-2.0);
+ Vkc = -p.C6*r6inv/TSvdw;
+ // derivatives
+ fpair = -6.0*p.C6*r8inv/TSvdw + p.d/p.seff*p.C6*(TSvdw - 1.0)*r6inv*TSvdw2inv/r;
+ forcelj = fpair;
+ fforce = factor_lj*(forcelj*Tap - Vkc*dTap/r);
+
+ philj = Vkc*Tap;
+ return factor_lj*philj;
+}
+
+
+/* ---------------------------------------------------------------------- */
+
+int PairILP::pack_forward_comm(int n, int *list, double *buf,
+ int pbc_flag, int *pbc)
+{
+ int i,j,m,id,ip,l;
+
+ m = 0;
+ for (i = 0; i < n; i++) {
+ j = list[i];
+ buf[m++] = normal[j][0];
+ buf[m++] = normal[j][1];
+ buf[m++] = normal[j][2];
+ buf[m++] = dnormdri[0][0][j];
+ buf[m++] = dnormdri[0][1][j];
+ buf[m++] = dnormdri[0][2][j];
+ buf[m++] = dnormdri[1][0][j];
+ buf[m++] = dnormdri[1][1][j];
+ buf[m++] = dnormdri[1][2][j];
+ buf[m++] = dnormdri[2][0][j];
+ buf[m++] = dnormdri[2][1][j];
+ buf[m++] = dnormdri[2][2][j];
+ for (l = 0; l < 3; l++){
+ for (id = 0; id < 3; id++){
+ for (ip = 0; ip < 3; ip++){
+ buf[m++] = dnormal[id][ip][l][j];
+ }
+ }
+ }
+ }
+ return m;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void PairILP::unpack_forward_comm(int n, int first, double *buf)
+{
+ int i,m,last,id,ip,l;
+
+ m = 0;
+ last = first + n;
+ for (i = first; i < last; i++) {
+ normal[i][0] = buf[m++];
+ normal[i][1] = buf[m++];
+ normal[i][2] = buf[m++];
+ dnormdri[0][0][i] = buf[m++];
+ dnormdri[0][1][i] = buf[m++];
+ dnormdri[0][2][i] = buf[m++];
+ dnormdri[1][0][i] = buf[m++];
+ dnormdri[1][1][i] = buf[m++];
+ dnormdri[1][2][i] = buf[m++];
+ dnormdri[2][0][i] = buf[m++];
+ dnormdri[2][1][i] = buf[m++];
+ dnormdri[2][2][i] = buf[m++];
+ for (l = 0; l < 3; l++){
+ for (id = 0; id < 3; id++){
+ for (ip = 0; ip < 3; ip++){
+ dnormal[id][ip][l][i] = buf[m++];
+ }
+ }
+ }
+ }
+}
+
+/* ---------------------------------------------------------------------- */
diff --git a/src/USER-MISC/pair_ilp_gr_hBN.h b/src/USER-MISC/pair_ilp_gr_hBN.h
new file mode 100644
index 0000000000000000000000000000000000000000..e25910d929eed85e3245852986f2e9d12758fd30
--- /dev/null
+++ b/src/USER-MISC/pair_ilp_gr_hBN.h
@@ -0,0 +1,143 @@
+/* ----------------------------------------------------------------------
+ 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.
+
+*/
+
diff --git a/src/USER-MISC/pair_kolmogorov_crespi_full.cpp b/src/USER-MISC/pair_kolmogorov_crespi_full.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..0068144e4a7056d7566e51988044eaa3521b8413
--- /dev/null
+++ b/src/USER-MISC/pair_kolmogorov_crespi_full.cpp
@@ -0,0 +1,1121 @@
+/* ----------------------------------------------------------------------
+ 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
+ based on previous versions by Jaap Kroes
+
+ This is a complete version of the potential described in
+ [Kolmogorov & Crespi, Phys. Rev. B 71, 235415 (2005)]
+------------------------------------------------------------------------- */
+
+#include <math.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <mpi.h>
+#include "pair_kolmogorov_crespi_full.h"
+#include "atom.h"
+#include "comm.h"
+#include "force.h"
+#include "neighbor.h"
+#include "neigh_list.h"
+#include "neigh_request.h"
+#include "my_page.h"
+#include "memory.h"
+#include "error.h"
+
+using namespace LAMMPS_NS;
+
+#define MAXLINE 1024
+#define DELTA 4
+#define PGDELTA 1
+
+/* ---------------------------------------------------------------------- */
+
+PairKolmogorovCrespiFull::PairKolmogorovCrespiFull(LAMMPS *lmp) : Pair(lmp)
+{
+ writedata = 1;
+
+ // initialize element to parameter maps
+ nelements = 0;
+ elements = NULL;
+ nparams = maxparam = 0;
+ params = NULL;
+ elem2param = NULL;
+ cutKCsq = NULL;
+ map = NULL;
+
+ nmax = 0;
+ maxlocal = 0;
+ KC_numneigh = NULL;
+ KC_firstneigh = NULL;
+ ipage = NULL;
+ pgsize = oneatom = 0;
+
+ normal = NULL;
+ dnormal = NULL;
+ dnormdri = NULL;
+
+ // always compute energy offset
+ offset_flag = 1;
+
+ // set comm size needed by this Pair
+ comm_forward = 39;
+ tap_flag = 0;
+}
+
+/* ---------------------------------------------------------------------- */
+
+PairKolmogorovCrespiFull::~PairKolmogorovCrespiFull()
+{
+ memory->destroy(KC_numneigh);
+ memory->sfree(KC_firstneigh);
+ delete [] ipage;
+ memory->destroy(normal);
+ memory->destroy(dnormal);
+ memory->destroy(dnormdri);
+
+ if (allocated) {
+ memory->destroy(setflag);
+ memory->destroy(cutsq);
+ memory->destroy(cut);
+ memory->destroy(offset);
+ }
+
+ if (elements)
+ for (int i = 0; i < nelements; i++) delete [] elements[i];
+ delete [] elements;
+ memory->destroy(params);
+ memory->destroy(elem2param);
+ memory->destroy(cutKCsq);
+ if (allocated) delete [] map;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void PairKolmogorovCrespiFull::compute(int eflag, int vflag)
+{
+ int i,j,ii,jj,inum,jnum,itype,jtype,id,k,l,kk,ll;
+ tagint itag,jtag;
+ double prodnorm1,prodnorm2,fkcx,fkcy,fkcz;
+ double xtmp,ytmp,ztmp,delx,dely,delz,evdwl,fpair,fpair1,fpair2;
+ double rsq,r,rhosq1,rhosq2,exp0,exp1,exp2,r2inv,r6inv,r8inv,Tap,dTap,Vkc;
+ double frho1,frho2,sumC1,sumC2,sumC11,sumC22,sumCff,fsum,rdsq1,rdsq2;
+ int *ilist,*jlist,*numneigh,**firstneigh;
+ int *KC_neighs_i,*KC_neighs_j;
+
+ evdwl = 0.0;
+ if (eflag || vflag) ev_setup(eflag,vflag);
+ else evflag = vflag_fdotr = 0;
+
+ double **x = atom->x;
+ double **f = atom->f;
+ int *type = atom->type;
+ tagint *tag = atom->tag;
+ int nlocal = atom->nlocal;
+ int newton_pair = force->newton_pair;
+ double dprodnorm1[3] = {0.0, 0.0, 0.0};
+ double dprodnorm2[3] = {0.0, 0.0, 0.0};
+ double fp1[3] = {0.0, 0.0, 0.0};
+ double fp2[3] = {0.0, 0.0, 0.0};
+ double fprod1[3] = {0.0, 0.0, 0.0};
+ double fprod2[3] = {0.0, 0.0, 0.0};
+
+ inum = list->inum;
+ ilist = list->ilist;
+ numneigh = list->numneigh;
+ firstneigh = list->firstneigh;
+ // Build full neighbor list
+ KC_neigh();
+ // Calculate the normals
+ calc_normal();
+
+ // communicate the normal vector and its derivatives
+ comm->forward_comm_pair(this);
+
+ // loop over neighbors of my atoms
+ for (ii = 0; ii < inum; ii++) {
+ i = ilist[ii];
+ itag = tag[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];
+ j &= NEIGHMASK;
+ jtype = type[j];
+ jtag = tag[j];
+
+ // two-body interactions from full neighbor list, skip half of them
+ if (itag > jtag) {
+ if ((itag+jtag) % 2 == 0) continue;
+ } else if (itag < jtag) {
+ if ((itag+jtag) % 2 == 1) continue;
+ } else {
+ if (x[j][2] < ztmp) continue;
+ if (x[j][2] == ztmp && x[j][1] < ytmp) continue;
+ if (x[j][2] == ztmp && x[j][1] == ytmp && x[j][0] < xtmp) continue;
+ }
+
+ delx = xtmp - x[j][0];
+ dely = ytmp - x[j][1];
+ delz = ztmp - x[j][2];
+ rsq = delx*delx + dely*dely + delz*delz;
+
+ // only include the interation between different layers
+ if (rsq < cutsq[itype][jtype] && atom->molecule[i] != atom->molecule[j]) {
+
+ int iparam_ij = elem2param[map[itype]][map[jtype]];
+ Param& p = params[iparam_ij];
+
+ r = sqrt(rsq);
+ r2inv = 1.0/rsq;
+ r6inv = r2inv*r2inv*r2inv;
+ r8inv = r2inv*r6inv;
+ // turn on/off Tapper function
+ if (tap_flag) {
+ Tap = calc_Tap(r,sqrt(cutsq[itype][jtype]));
+ dTap = calc_dTap(r,sqrt(cutsq[itype][jtype]));
+ } else {Tap = 1.0; dTap = 0.0;}
+
+ // Calculate the transverse distance
+ // note that rho_ij does not equal to rho_ji except when normals are all along z
+ prodnorm1 = normal[i][0]*delx + normal[i][1]*dely + normal[i][2]*delz;
+ prodnorm2 = normal[j][0]*delx + normal[j][1]*dely + normal[j][2]*delz;
+ rhosq1 = rsq - prodnorm1*prodnorm1; // rho_ij
+ rhosq2 = rsq - prodnorm2*prodnorm2; // rho_ji
+ rdsq1 = rhosq1*p.delta2inv; // (rho_ij/delta)^2
+ rdsq2 = rhosq2*p.delta2inv; // (rho_ji/delta)^2
+
+ // store exponents
+ exp0 = exp(-p.lambda*(r-p.z0));
+ exp1 = exp(-rdsq1);
+ exp2 = exp(-rdsq2);
+
+ sumC1 = p.C0 + p.C2*rdsq1 + p.C4*rdsq1*rdsq1;
+ sumC2 = p.C0 + p.C2*rdsq2 + p.C4*rdsq2*rdsq2;
+ sumC11 = (p.C2 + 2.0*p.C4*rdsq1)*p.delta2inv;
+ sumC22 = (p.C2 + 2.0*p.C4*rdsq2)*p.delta2inv;
+ frho1 = exp1*sumC1;
+ frho2 = exp2*sumC2;
+ sumCff = p.C + frho1 + frho2;
+ Vkc = -p.A*p.z06*r6inv + exp0*sumCff;
+
+ // derivatives
+ fpair = -6.0*p.A*p.z06*r8inv + p.lambda*exp0/r*sumCff;
+ fpair1 = 2.0*exp0*exp1*(p.delta2inv*sumC1 - sumC11);
+ fpair2 = 2.0*exp0*exp2*(p.delta2inv*sumC2 - sumC22);
+ fsum = fpair + fpair1 + fpair2;
+ // derivatives of the product of rij and ni, the result is a vector
+ dprodnorm1[0] = dnormdri[0][0][i]*delx + dnormdri[1][0][i]*dely + dnormdri[2][0][i]*delz;
+ dprodnorm1[1] = dnormdri[0][1][i]*delx + dnormdri[1][1][i]*dely + dnormdri[2][1][i]*delz;
+ dprodnorm1[2] = dnormdri[0][2][i]*delx + dnormdri[1][2][i]*dely + dnormdri[2][2][i]*delz;
+ // derivatives of the product of rji and nj, the result is a vector
+ dprodnorm2[0] = dnormdri[0][0][j]*delx + dnormdri[1][0][j]*dely + dnormdri[2][0][j]*delz;
+ dprodnorm2[1] = dnormdri[0][1][j]*delx + dnormdri[1][1][j]*dely + dnormdri[2][1][j]*delz;
+ dprodnorm2[2] = dnormdri[0][2][j]*delx + dnormdri[1][2][j]*dely + dnormdri[2][2][j]*delz;
+ fp1[0] = prodnorm1*normal[i][0]*fpair1;
+ fp1[1] = prodnorm1*normal[i][1]*fpair1;
+ fp1[2] = prodnorm1*normal[i][2]*fpair1;
+ fp2[0] = prodnorm2*normal[j][0]*fpair2;
+ fp2[1] = prodnorm2*normal[j][1]*fpair2;
+ fp2[2] = prodnorm2*normal[j][2]*fpair2;
+ fprod1[0] = prodnorm1*dprodnorm1[0]*fpair1;
+ fprod1[1] = prodnorm1*dprodnorm1[1]*fpair1;
+ fprod1[2] = prodnorm1*dprodnorm1[2]*fpair1;
+ fprod2[0] = prodnorm2*dprodnorm2[0]*fpair2;
+ fprod2[1] = prodnorm2*dprodnorm2[1]*fpair2;
+ fprod2[2] = prodnorm2*dprodnorm2[2]*fpair2;
+ fkcx = (delx*fsum - fp1[0] - fp2[0])*Tap - Vkc*dTap*delx/r;
+ fkcy = (dely*fsum - fp1[1] - fp2[1])*Tap - Vkc*dTap*dely/r;
+ fkcz = (delz*fsum - fp1[2] - fp2[2])*Tap - Vkc*dTap*delz/r;
+
+ f[i][0] += fkcx - fprod1[0]*Tap;
+ f[i][1] += fkcy - fprod1[1]*Tap;
+ f[i][2] += fkcz - fprod1[2]*Tap;
+ if (newton_pair || j < nlocal) {
+ f[j][0] -= fkcx + fprod2[0]*Tap;
+ f[j][1] -= fkcy + fprod2[1]*Tap;
+ f[j][2] -= fkcz + fprod2[2]*Tap;
+ }
+
+ // calculate the forces acted on the neighbors of atom i from atom j
+ KC_neighs_i = KC_firstneigh[i];
+ for (kk = 0; kk < KC_numneigh[i]; kk++) {
+ k = KC_neighs_i[kk];
+ if (k == i) continue;
+ // derivatives of the product of rij and ni respect to rk, k=0,1,2, where atom k is the neighbors of atom i
+ dprodnorm1[0] = dnormal[0][0][kk][i]*delx + dnormal[1][0][kk][i]*dely + dnormal[2][0][kk][i]*delz;
+ dprodnorm1[1] = dnormal[0][1][kk][i]*delx + dnormal[1][1][kk][i]*dely + dnormal[2][1][kk][i]*delz;
+ dprodnorm1[2] = dnormal[0][2][kk][i]*delx + dnormal[1][2][kk][i]*dely + dnormal[2][2][kk][i]*delz;
+ f[k][0] += (-prodnorm1*dprodnorm1[0]*fpair1)*Tap;
+ f[k][1] += (-prodnorm1*dprodnorm1[1]*fpair1)*Tap;
+ f[k][2] += (-prodnorm1*dprodnorm1[2]*fpair1)*Tap;
+ }
+
+ // calculate the forces acted on the neighbors of atom j from atom i
+ KC_neighs_j = KC_firstneigh[j];
+ for (ll = 0; ll < KC_numneigh[j]; ll++) {
+ l = KC_neighs_j[ll];
+ if (l == j) continue;
+ if (newton_pair || l < nlocal) {
+ // derivatives of the product of rji and nj respect to rl, l=0,1,2, where atom l is the neighbors of atom j
+ dprodnorm2[0] = dnormal[0][0][ll][j]*delx + dnormal[1][0][ll][j]*dely + dnormal[2][0][ll][j]*delz;
+ dprodnorm2[1] = dnormal[0][1][ll][j]*delx + dnormal[1][1][ll][j]*dely + dnormal[2][1][ll][j]*delz;
+ dprodnorm2[2] = dnormal[0][2][ll][j]*delx + dnormal[1][2][ll][j]*dely + dnormal[2][2][ll][j]*delz;
+ f[l][0] += (-prodnorm2*dprodnorm2[0]*fpair2)*Tap;
+ f[l][1] += (-prodnorm2*dprodnorm2[1]*fpair2)*Tap;
+ f[l][2] += (-prodnorm2*dprodnorm2[2]*fpair2)*Tap;
+ }
+ }
+
+ if (eflag) {
+ if (tap_flag) evdwl = Tap*Vkc;
+ else evdwl = Vkc - offset[itype][jtype];
+ }
+
+ if (evflag){
+ ev_tally_xyz(i,j,nlocal,newton_pair,evdwl,0,
+ fkcx,fkcy,fkcz,delx,dely,delz);
+ }
+ }
+ }
+ }
+ if (vflag_fdotr) virial_fdotr_compute();
+}
+
+/* ----------------------------------------------------------------------
+ Calculate the normals for each atom
+------------------------------------------------------------------------- */
+void PairKolmogorovCrespiFull::calc_normal()
+{
+ int i,j,ii,jj,inum,jnum,itype,jtype;
+ int cont,k,kk,cont1,id,ip,m;
+ tagint itag,jtag;
+ double dist,nn,xtp,ytp,ztp,delx,dely,delz,nn2;
+ int *ilist,*jlist;
+ double pv12[3],pv31[3],pv23[3],n1[3],dni[3],dnn[3][3],vet[3][3],dpvdri[3][3];
+ double dn1[3][3][3],dpv12[3][3][3],dpv23[3][3][3],dpv31[3][3][3];
+
+ double **x = atom->x;
+ int *type = atom->type;
+ int ntypes = atom->ntypes;
+ int nlocal = atom->nlocal;
+ tagint *tag = atom->tag;
+
+ // grow normal array if necessary
+
+ if (atom->nmax > nmax) {
+ memory->destroy(normal);
+ memory->destroy(dnormal);
+ memory->destroy(dnormdri);
+ nmax = atom->nmax;
+ memory->create(normal,nmax,3,"KolmogorovCrespiFull:normal");
+ memory->create(dnormdri,3,3,nmax,"KolmogorovCrespiFull:dnormdri");
+ memory->create(dnormal,3,3,3,nmax,"KolmogorovCrespiFull:dnormal");
+ }
+
+ inum = list->inum;
+ ilist = list->ilist;
+ //Calculate normals
+ for (ii = 0; ii < inum; ii++) {
+ i = ilist[ii];
+ xtp = x[i][0];
+ ytp = x[i][1];
+ ztp = x[i][2];
+ itype = type[i];
+
+ // Initialize the arrays
+ for (id = 0; id < 3; id++){
+ pv12[id] = 0.0;
+ pv31[id] = 0.0;
+ pv23[id] = 0.0;
+ n1[id] = 0.0;
+ dni[id] = 0.0;
+ normal[i][id] = 0.0;
+ for (ip = 0; ip < 3; ip++){
+ vet[ip][id] = 0.0;
+ dnn[ip][id] = 0.0;
+ dpvdri[ip][id] = 0.0;
+ dnormdri[ip][id][i] = 0.0;
+ for (m = 0; m < 3; m++){
+ dpv12[ip][id][m] = 0.0;
+ dpv31[ip][id][m] = 0.0;
+ dpv23[ip][id][m] = 0.0;
+ dn1[ip][id][m] = 0.0;
+ dnormal[ip][id][m][i] = 0.0;
+ }
+ }
+ }
+
+ cont = 0;
+ jlist = KC_firstneigh[i];
+ jnum = KC_numneigh[i];
+ for (jj = 0; jj < jnum; jj++) {
+ j = jlist[jj];
+ j &= NEIGHMASK;
+
+ delx = x[j][0] - xtp;
+ dely = x[j][1] - ytp;
+ delz = x[j][2] - ztp;
+ vet[cont][0] = delx;
+ vet[cont][1] = dely;
+ vet[cont][2] = delz;
+ cont++;
+ }
+
+ if (cont <= 1) {
+ normal[i][0] = 0.0;
+ normal[i][1] = 0.0;
+ normal[i][2] = 1.0;
+ // derivatives of normal vector is zero
+ for (id = 0; id < 3; id++){
+ for (ip = 0; ip < 3; ip++){
+ dnormdri[id][ip][i] = 0.0;
+ for (m = 0; m < 3; m++){
+ dnormal[id][ip][m][i] = 0.0;
+ }
+ }
+ }
+ }
+ else if (cont == 2) {
+ // for the atoms at the edge who has only two neighbor atoms
+ pv12[0] = vet[0][1]*vet[1][2] - vet[1][1]*vet[0][2];
+ pv12[1] = vet[0][2]*vet[1][0] - vet[1][2]*vet[0][0];
+ pv12[2] = vet[0][0]*vet[1][1] - vet[1][0]*vet[0][1];
+ dpvdri[0][0] = 0.0;
+ dpvdri[0][1] = vet[0][2]-vet[1][2];
+ dpvdri[0][2] = vet[1][1]-vet[0][1];
+ dpvdri[1][0] = vet[1][2]-vet[0][2];
+ dpvdri[1][1] = 0.0;
+ dpvdri[1][2] = vet[0][0]-vet[1][0];
+ dpvdri[2][0] = vet[0][1]-vet[1][1];
+ dpvdri[2][1] = vet[1][0]-vet[0][0];
+ dpvdri[2][2] = 0.0;
+
+ // derivatives respect to the first neighbor, atom k
+ dpv12[0][0][0] = 0.0;
+ dpv12[0][1][0] = vet[1][2];
+ dpv12[0][2][0] = -vet[1][1];
+ dpv12[1][0][0] = -vet[1][2];
+ dpv12[1][1][0] = 0.0;
+ dpv12[1][2][0] = vet[1][0];
+ dpv12[2][0][0] = vet[1][1];
+ dpv12[2][1][0] = -vet[1][0];
+ dpv12[2][2][0] = 0.0;
+
+ // derivatives respect to the second neighbor, atom l
+ dpv12[0][0][1] = 0.0;
+ dpv12[0][1][1] = -vet[0][2];
+ dpv12[0][2][1] = vet[0][1];
+ dpv12[1][0][1] = vet[0][2];
+ dpv12[1][1][1] = 0.0;
+ dpv12[1][2][1] = -vet[0][0];
+ dpv12[2][0][1] = -vet[0][1];
+ dpv12[2][1][1] = vet[0][0];
+ dpv12[2][2][1] = 0.0;
+
+ // derivatives respect to the third neighbor, atom n
+ for (id = 0; id < 3; id++){
+ for (ip = 0; ip < 3; ip++){
+ dpv12[id][ip][2] = 0.0;
+ }
+ }
+
+ n1[0] = pv12[0];
+ n1[1] = pv12[1];
+ n1[2] = pv12[2];
+ // the magnitude of the normal vector
+ nn2 = n1[0]*n1[0] + n1[1]*n1[1] + n1[2]*n1[2];
+ nn = sqrt(nn2);
+ if (nn == 0) error->all(FLERR,"The magnitude of the normal vector is zero");
+ // the unit normal vector
+ normal[i][0] = n1[0]/nn;
+ normal[i][1] = n1[1]/nn;
+ normal[i][2] = n1[2]/nn;
+ // derivatives of nn, dnn:3x1 vector
+ dni[0] = (n1[0]*dpvdri[0][0] + n1[1]*dpvdri[1][0] + n1[2]*dpvdri[2][0])/nn;
+ dni[1] = (n1[0]*dpvdri[0][1] + n1[1]*dpvdri[1][1] + n1[2]*dpvdri[2][1])/nn;
+ dni[2] = (n1[0]*dpvdri[0][2] + n1[1]*dpvdri[1][2] + n1[2]*dpvdri[2][2])/nn;
+ // derivatives of unit vector ni respect to ri, the result is 3x3 matrix
+ for (id = 0; id < 3; id++){
+ for (ip = 0; ip < 3; ip++){
+ dnormdri[id][ip][i] = dpvdri[id][ip]/nn - n1[id]*dni[ip]/nn2;
+ }
+ }
+
+ // derivatives of non-normalized normal vector, dn1:3x3x3 array
+ for (id = 0; id < 3; id++){
+ for (ip = 0; ip < 3; ip++){
+ for (m = 0; m < 3; m++){
+ dn1[id][ip][m] = dpv12[id][ip][m];
+ }
+ }
+ }
+ // derivatives of nn, dnn:3x3 vector
+ // dnn[id][m]: the derivative of nn respect to r[id][m], id,m=0,1,2
+ // r[id][m]: the id's component of atom m
+ for (m = 0; m < 3; m++){
+ for (id = 0; id < 3; id++){
+ dnn[id][m] = (n1[0]*dn1[0][id][m] + n1[1]*dn1[1][id][m] + n1[2]*dn1[2][id][m])/nn;
+ }
+ }
+ // dnormal[id][ip][m][i]: the derivative of normal[id] respect to r[ip][m], id,ip=0,1,2
+ // for atom m, which is a neighbor atom of atom i, m=0,jnum-1
+ for (m = 0; m < 3; m++){
+ for (id = 0; id < 3; id++){
+ for (ip = 0; ip < 3; ip++){
+ dnormal[id][ip][m][i] = dn1[id][ip][m]/nn - n1[id]*dnn[ip][m]/nn2;
+ }
+ }
+ }
+ }
+//##############################################################################################
+
+ else if(cont == 3) {
+ // for the atoms at the edge who has only two neighbor atoms
+ pv12[0] = vet[0][1]*vet[1][2] - vet[1][1]*vet[0][2];
+ pv12[1] = vet[0][2]*vet[1][0] - vet[1][2]*vet[0][0];
+ pv12[2] = vet[0][0]*vet[1][1] - vet[1][0]*vet[0][1];
+ // derivatives respect to the first neighbor, atom k
+ dpv12[0][0][0] = 0.0;
+ dpv12[0][1][0] = vet[1][2];
+ dpv12[0][2][0] = -vet[1][1];
+ dpv12[1][0][0] = -vet[1][2];
+ dpv12[1][1][0] = 0.0;
+ dpv12[1][2][0] = vet[1][0];
+ dpv12[2][0][0] = vet[1][1];
+ dpv12[2][1][0] = -vet[1][0];
+ dpv12[2][2][0] = 0.0;
+ // derivatives respect to the second neighbor, atom l
+ dpv12[0][0][1] = 0.0;
+ dpv12[0][1][1] = -vet[0][2];
+ dpv12[0][2][1] = vet[0][1];
+ dpv12[1][0][1] = vet[0][2];
+ dpv12[1][1][1] = 0.0;
+ dpv12[1][2][1] = -vet[0][0];
+ dpv12[2][0][1] = -vet[0][1];
+ dpv12[2][1][1] = vet[0][0];
+ dpv12[2][2][1] = 0.0;
+
+ // derivatives respect to the third neighbor, atom n
+ for (id = 0; id < 3; id++){
+ for (ip = 0; ip < 3; ip++){
+ dpv12[id][ip][2] = 0.0;
+ }
+ }
+
+ pv31[0] = vet[2][1]*vet[0][2] - vet[0][1]*vet[2][2];
+ pv31[1] = vet[2][2]*vet[0][0] - vet[0][2]*vet[2][0];
+ pv31[2] = vet[2][0]*vet[0][1] - vet[0][0]*vet[2][1];
+ // derivatives respect to the first neighbor, atom k
+ dpv31[0][0][0] = 0.0;
+ dpv31[0][1][0] = -vet[2][2];
+ dpv31[0][2][0] = vet[2][1];
+ dpv31[1][0][0] = vet[2][2];
+ dpv31[1][1][0] = 0.0;
+ dpv31[1][2][0] = -vet[2][0];
+ dpv31[2][0][0] = -vet[2][1];
+ dpv31[2][1][0] = vet[2][0];
+ dpv31[2][2][0] = 0.0;
+ // derivatives respect to the third neighbor, atom n
+ dpv31[0][0][2] = 0.0;
+ dpv31[0][1][2] = vet[0][2];
+ dpv31[0][2][2] = -vet[0][1];
+ // derivatives of pv13[1] to rn
+ dpv31[1][0][2] = -vet[0][2];
+ dpv31[1][1][2] = 0.0;
+ dpv31[1][2][2] = vet[0][0];
+ // derivatives of pv13[2] to rn
+ dpv31[2][0][2] = vet[0][1];
+ dpv31[2][1][2] = -vet[0][0];
+ dpv31[2][2][2] = 0.0;
+
+ // derivatives respect to the second neighbor, atom l
+ for (id = 0; id < 3; id++){
+ for (ip = 0; ip < 3; ip++){
+ dpv31[id][ip][1] = 0.0;
+ }
+ }
+
+ pv23[0] = vet[1][1]*vet[2][2] - vet[2][1]*vet[1][2];
+ pv23[1] = vet[1][2]*vet[2][0] - vet[2][2]*vet[1][0];
+ pv23[2] = vet[1][0]*vet[2][1] - vet[2][0]*vet[1][1];
+ // derivatives respect to the second neighbor, atom k
+ for (id = 0; id < 3; id++){
+ for (ip = 0; ip < 3; ip++){
+ dpv23[id][ip][0] = 0.0;
+ }
+ }
+ // derivatives respect to the second neighbor, atom l
+ dpv23[0][0][1] = 0.0;
+ dpv23[0][1][1] = vet[2][2];
+ dpv23[0][2][1] = -vet[2][1];
+ dpv23[1][0][1] = -vet[2][2];
+ dpv23[1][1][1] = 0.0;
+ dpv23[1][2][1] = vet[2][0];
+ dpv23[2][0][1] = vet[2][1];
+ dpv23[2][1][1] = -vet[2][0];
+ dpv23[2][2][1] = 0.0;
+ // derivatives respect to the third neighbor, atom n
+ dpv23[0][0][2] = 0.0;
+ dpv23[0][1][2] = -vet[1][2];
+ dpv23[0][2][2] = vet[1][1];
+ dpv23[1][0][2] = vet[1][2];
+ dpv23[1][1][2] = 0.0;
+ dpv23[1][2][2] = -vet[1][0];
+ dpv23[2][0][2] = -vet[1][1];
+ dpv23[2][1][2] = vet[1][0];
+ dpv23[2][2][2] = 0.0;
+
+//############################################################################################
+ // average the normal vectors by using the 3 neighboring planes
+ n1[0] = (pv12[0] + pv31[0] + pv23[0])/cont;
+ n1[1] = (pv12[1] + pv31[1] + pv23[1])/cont;
+ n1[2] = (pv12[2] + pv31[2] + pv23[2])/cont;
+ // the magnitude of the normal vector
+ nn2 = n1[0]*n1[0] + n1[1]*n1[1] + n1[2]*n1[2];
+ nn = sqrt(nn2);
+ if (nn == 0) error->all(FLERR,"The magnitude of the normal vector is zero");
+ // the unit normal vector
+ normal[i][0] = n1[0]/nn;
+ normal[i][1] = n1[1]/nn;
+ normal[i][2] = n1[2]/nn;
+
+ // for the central atoms, dnormdri is always zero
+ for (id = 0; id < 3; id++){
+ for (ip = 0; ip < 3; ip++){
+ dnormdri[id][ip][i] = 0.0;
+ }
+ } // end of derivatives of normals respect to atom i
+
+ // derivatives of non-normalized normal vector, dn1:3x3x3 array
+ for (id = 0; id < 3; id++){
+ for (ip = 0; ip < 3; ip++){
+ for (m = 0; m < 3; m++){
+ dn1[id][ip][m] = (dpv12[id][ip][m] + dpv23[id][ip][m] + dpv31[id][ip][m])/cont;
+ }
+ }
+ }
+ // derivatives of nn, dnn:3x3 vector
+ // dnn[id][m]: the derivative of nn respect to r[id][m], id,m=0,1,2
+ // r[id][m]: the id's component of atom m
+ for (m = 0; m < 3; m++){
+ for (id = 0; id < 3; id++){
+ dnn[id][m] = (n1[0]*dn1[0][id][m] + n1[1]*dn1[1][id][m] + n1[2]*dn1[2][id][m])/nn;
+ }
+ }
+ // dnormal[id][ip][m][i]: the derivative of normal[id] respect to r[ip][m], id,ip=0,1,2
+ // for atom m, which is a neighbor atom of atom i, m=0,jnum-1
+ for (m = 0; m < 3; m++){
+ for (id = 0; id < 3; id++){
+ for (ip = 0; ip < 3; ip++){
+ dnormal[id][ip][m][i] = dn1[id][ip][m]/nn - n1[id]*dnn[ip][m]/nn2;
+ }
+ }
+ }
+ }
+ else {
+ error->all(FLERR,"There are too many neighbors for calculating normals");
+ }
+
+//##############################################################################################
+ }
+}
+
+/* ----------------------------------------------------------------------
+ init specific to this pair style
+------------------------------------------------------------------------- */
+
+void PairKolmogorovCrespiFull::init_style()
+{
+ if (force->newton_pair == 0)
+ error->all(FLERR,"Pair style KC requires newton pair on");
+
+ // need a full neighbor list, including neighbors of ghosts
+
+ int irequest = neighbor->request(this,instance_me);
+ neighbor->requests[irequest]->half = 0;
+ neighbor->requests[irequest]->full = 1;
+ neighbor->requests[irequest]->ghost = 1;
+
+ // local KC neighbor list
+ // create pages if first time or if neighbor pgsize/oneatom has changed
+
+ int create = 0;
+ if (ipage == NULL) create = 1;
+ if (pgsize != neighbor->pgsize) create = 1;
+ if (oneatom != neighbor->oneatom) create = 1;
+
+ if (create) {
+ delete [] ipage;
+ pgsize = neighbor->pgsize;
+ oneatom = neighbor->oneatom;
+
+ int nmypage= comm->nthreads;
+ ipage = new MyPage<int>[nmypage];
+ for (int i = 0; i < nmypage; i++)
+ ipage[i].init(oneatom,pgsize,PGDELTA);
+ }
+}
+
+
+/* ----------------------------------------------------------------------
+ create KC neighbor list from main neighbor list
+ KC neighbor list stores neighbors of ghost atoms
+ KC_numneigh for calcualting normals and
+ KC_pair_numneigh for calculating force
+------------------------------------------------------------------------- */
+
+void PairKolmogorovCrespiFull::KC_neigh()
+{
+ int i,j,ii,jj,n,np,allnum,inum,jnum,itype,jtype;
+ double xtmp,ytmp,ztmp,delx,dely,delz,rsq,dS;
+ int *ilist,*jlist,*numneigh,**firstneigh;
+ int *neighptr;
+
+ double **x = atom->x;
+ int *type = atom->type;
+
+ if (atom->nmax > maxlocal) {
+ maxlocal = atom->nmax;
+ memory->destroy(KC_numneigh);
+ memory->sfree(KC_firstneigh);
+ memory->create(KC_numneigh,maxlocal,"KolmogorovCrespiFull:numneigh");
+ KC_firstneigh = (int **) memory->smalloc(maxlocal*sizeof(int *),
+ "KolmogorovCrespiFull:firstneigh");
+ }
+
+ inum = list->inum;
+ allnum = list->inum + list->gnum;
+ ilist = list->ilist;
+ numneigh = list->numneigh;
+ firstneigh = list->firstneigh;
+
+ // store all KC neighs of owned and ghost atoms
+ // scan full neighbor list of I
+
+ ipage->reset();
+
+ for (ii = 0; ii < allnum; ii++) {
+ i = ilist[ii];
+
+ n = 0;
+ neighptr = ipage->vget();
+
+ xtmp = x[i][0];
+ ytmp = x[i][1];
+ ztmp = x[i][2];
+ itype = map[type[i]];
+ jlist = firstneigh[i];
+ jnum = numneigh[i];
+
+ for (jj = 0; jj < jnum; jj++) {
+ j = jlist[jj];
+ j &= NEIGHMASK;
+ jtype = map[type[j]];
+ delx = xtmp - x[j][0];
+ dely = ytmp - x[j][1];
+ delz = ztmp - x[j][2];
+ rsq = delx*delx + dely*dely + delz*delz;
+
+ if (rsq != 0 && rsq < cutKCsq[itype][jtype] && atom->molecule[i] == atom->molecule[j]) {
+ neighptr[n++] = j;
+ }
+ }
+
+ KC_firstneigh[i] = neighptr;
+ KC_numneigh[i] = n;
+ ipage->vgot(n);
+ if (ipage->status())
+ error->one(FLERR,"Neighbor list overflow, boost neigh_modify one");
+ }
+}
+
+
+/* ----------------------------------------------------------------------
+ allocate all arrays
+------------------------------------------------------------------------- */
+
+void PairKolmogorovCrespiFull::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(offset,n+1,n+1,"pair:offset");
+ map = new int[atom->ntypes+1];
+}
+
+/* ----------------------------------------------------------------------
+ global settings
+------------------------------------------------------------------------- */
+
+void PairKolmogorovCrespiFull::settings(int narg, char **arg)
+{
+ if (narg < 1 || narg > 2) error->all(FLERR,"Illegal pair_style command");
+ if (strcmp(force->pair_style,"hybrid/overlay")!=0)
+ error->all(FLERR,"ERROR: requires hybrid/overlay pair_style");
+
+ cut_global = force->numeric(FLERR,arg[0]);
+ if (narg == 2) tap_flag = 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; j <= atom->ntypes; j++)
+ if (setflag[i][j]) cut[i][j] = cut_global;
+ }
+}
+
+/* ----------------------------------------------------------------------
+ set coeffs for one or more type pairs
+------------------------------------------------------------------------- */
+
+void PairKolmogorovCrespiFull::coeff(int narg, char **arg)
+{
+ int i,j,n;
+
+ if (narg != 3 + atom->ntypes)
+ 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);
+
+ // read args that map atom types to elements in potential file
+ // map[i] = which element the Ith atom type is, -1 if NULL
+ // nelements = # of unique elements
+ // elements = list of element names
+
+ if (elements) {
+ for (i = 0; i < nelements; i++) delete [] elements[i];
+ delete [] elements;
+ }
+ elements = new char*[atom->ntypes];
+ for (i = 0; i < atom->ntypes; i++) elements[i] = NULL;
+
+ nelements = 0;
+ for (i = 3; i < narg; i++) {
+ if (strcmp(arg[i],"NULL") == 0) {
+ map[i-2] = -1;
+ continue;
+ }
+ for (j = 0; j < nelements; j++)
+ if (strcmp(arg[i],elements[j]) == 0) break;
+ map[i-2] = j;
+ if (j == nelements) {
+ n = strlen(arg[i]) + 1;
+ elements[j] = new char[n];
+ strcpy(elements[j],arg[i]);
+ nelements++;
+ }
+ }
+
+
+ read_file(arg[2]);
+
+ double cut_one = cut_global;
+
+ int count = 0;
+ for (int i = ilo; i <= ihi; i++) {
+ for (int j = MAX(jlo,i); j <= jhi; j++) {
+ cut[i][j] = cut_one;
+ setflag[i][j] = 1;
+ count++;
+ }
+ }
+
+ if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients");
+}
+
+
+/* ----------------------------------------------------------------------
+ init for one type pair i,j and corresponding j,i
+------------------------------------------------------------------------- */
+
+double PairKolmogorovCrespiFull::init_one(int i, int j)
+{
+ if (setflag[i][j] == 0) error->all(FLERR,"All pair coeffs are not set");
+
+ if (offset_flag && (cut[i][j] > 0.0)) {
+ int iparam_ij = elem2param[map[i]][map[j]];
+ Param& p = params[iparam_ij];
+ offset[i][j] = -p.A*pow(p.z0/cut[i][j],6);
+ } else offset[i][j] = 0.0;
+ offset[j][i] = offset[i][j];
+
+ return cut[i][j];
+}
+
+/* ----------------------------------------------------------------------
+ read Kolmogorov-Crespi potential file
+------------------------------------------------------------------------- */
+
+void PairKolmogorovCrespiFull::read_file(char *filename)
+{
+ int params_per_line = 12;
+ 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(filename);
+ if (fp == NULL) {
+ char str[128];
+ sprintf(str,"Cannot open KC potential file %s",filename);
+ 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 i,j,n,m,nwords,ielement,jelement;
+ 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,"Insufficient format in KC 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 = 1st args
+ // if these 2 args are in element list, then parse this line
+ // else skip to next line (continue)
+
+ 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;
+
+ // 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].z0 = atof(words[2]);
+ params[nparams].C0 = atof(words[3]);
+ params[nparams].C2 = atof(words[4]);
+ params[nparams].C4 = atof(words[5]);
+ params[nparams].C = atof(words[6]);
+ params[nparams].delta = atof(words[7]);
+ params[nparams].lambda = atof(words[8]);
+ params[nparams].A = atof(words[9]);
+ // S provides a convenient scaling of all energies
+ params[nparams].S = atof(words[10]);
+ params[nparams].rcut = atof(words[11]);
+
+ // energies in meV further scaled by S
+ double meV = 1.0e-3*params[nparams].S;
+ params[nparams].C *= meV;
+ params[nparams].A *= meV;
+ params[nparams].C0 *= meV;
+ params[nparams].C2 *= meV;
+ params[nparams].C4 *= meV;
+
+ // precompute some quantities
+ params[nparams].delta2inv = pow(params[nparams].delta,-2);
+ params[nparams].z06 = pow(params[nparams].z0,6);
+
+ nparams++;
+ //if(nparams >= pow(atom->ntypes,3)) break;
+ }
+ memory->destroy(elem2param);
+ memory->destroy(cutKCsq);
+ memory->create(elem2param,nelements,nelements,"pair:elem2param");
+ memory->create(cutKCsq,nelements,nelements,"pair:cutKCsq");
+ for (i = 0; i < nelements; i++) {
+ for (j = 0; j < nelements; j++) {
+ n = -1;
+ for (m = 0; m < nparams; m++) {
+ if (i == params[m].ielement && j == params[m].jelement) {
+ if (n >= 0) error->all(FLERR,"Potential file has duplicate entry");
+ n = m;
+ }
+ }
+ if (n < 0) error->all(FLERR,"Potential file is missing an entry");
+ elem2param[i][j] = n;
+ cutKCsq[i][j] = params[n].rcut*params[n].rcut;
+ }
+ }
+ delete [] words;
+}
+
+/* ---------------------------------------------------------------------- */
+
+double PairKolmogorovCrespiFull::single(int i, int j, int itype, int jtype, double rsq,
+ double factor_coul, double factor_lj,
+ double &fforce)
+{
+ double r,r2inv,r6inv,r8inv,forcelj,philj;
+ double Tap,dTap,Vkc,fpair;
+
+ int iparam_ij = elem2param[map[itype]][map[jtype]];
+ Param& p = params[iparam_ij];
+
+ r = sqrt(rsq);
+ // turn on/off Tapper function
+ if (tap_flag) {
+ Tap = calc_Tap(r,sqrt(cutsq[itype][jtype]));
+ dTap = calc_dTap(r,sqrt(cutsq[itype][jtype]));
+ } else {Tap = 1.0; dTap = 0.0;}
+
+ r2inv = 1.0/rsq;
+ r6inv = r2inv*r2inv*r2inv;
+ r8inv = r2inv*r6inv;
+
+ Vkc = -p.A*p.z06*r6inv;
+ // derivatives
+ fpair = -6.0*p.A*p.z06*r8inv;
+ forcelj = fpair;
+ fforce = factor_lj*(forcelj*Tap - Vkc*dTap/r);
+
+ if (tap_flag) philj = Vkc*Tap;
+ else philj = Vkc - offset[itype][jtype];
+ return factor_lj*philj;
+}
+
+/* ---------------------------------------------------------------------- */
+
+int PairKolmogorovCrespiFull::pack_forward_comm(int n, int *list, double *buf,
+ int pbc_flag, int *pbc)
+{
+ int i,j,m,l,ip,id;
+
+ m = 0;
+ for (i = 0; i < n; i++) {
+ j = list[i];
+ buf[m++] = normal[j][0];
+ buf[m++] = normal[j][1];
+ buf[m++] = normal[j][2];
+ buf[m++] = dnormdri[0][0][j];
+ buf[m++] = dnormdri[0][1][j];
+ buf[m++] = dnormdri[0][2][j];
+ buf[m++] = dnormdri[1][0][j];
+ buf[m++] = dnormdri[1][1][j];
+ buf[m++] = dnormdri[1][2][j];
+ buf[m++] = dnormdri[2][0][j];
+ buf[m++] = dnormdri[2][1][j];
+ buf[m++] = dnormdri[2][2][j];
+ for (l = 0; l < 3; l++){
+ for (id = 0; id < 3; id++){
+ for (ip = 0; ip < 3; ip++){
+ buf[m++] = dnormal[id][ip][l][j];
+ }
+ }
+ }
+ }
+ return m;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void PairKolmogorovCrespiFull::unpack_forward_comm(int n, int first, double *buf)
+{
+ int i,m,last,l,ip,id;
+
+ m = 0;
+ last = first + n;
+ for (i = first; i < last; i++) {
+ normal[i][0] = buf[m++];
+ normal[i][1] = buf[m++];
+ normal[i][2] = buf[m++];
+ dnormdri[0][0][i] = buf[m++];
+ dnormdri[0][1][i] = buf[m++];
+ dnormdri[0][2][i] = buf[m++];
+ dnormdri[1][0][i] = buf[m++];
+ dnormdri[1][1][i] = buf[m++];
+ dnormdri[1][2][i] = buf[m++];
+ dnormdri[2][0][i] = buf[m++];
+ dnormdri[2][1][i] = buf[m++];
+ dnormdri[2][2][i] = buf[m++];
+ for (l = 0; l < 3; l++){
+ for (id = 0; id < 3; id++){
+ for (ip = 0; ip < 3; ip++){
+ dnormal[id][ip][l][i] = buf[m++];
+ }
+ }
+ }
+ }
+}
+
+/* ---------------------------------------------------------------------- */
diff --git a/src/USER-MISC/pair_kolmogorov_crespi_full.h b/src/USER-MISC/pair_kolmogorov_crespi_full.h
new file mode 100644
index 0000000000000000000000000000000000000000..74ba3c4bc8336aa11a2cc2358020f0cdf13af1e5
--- /dev/null
+++ b/src/USER-MISC/pair_kolmogorov_crespi_full.h
@@ -0,0 +1,147 @@
+/* ----------------------------------------------------------------------
+ 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.
+
+*/
+