diff --git a/doc/src/pair_coul_shield.txt b/doc/src/pair_coul_shield.txt
index 62ac157e76cf4b78797dfc2a073d9007c7eee270..df04e76de9d296a28b10723e15ee7b6c005c05ad 100644
--- a/doc/src/pair_coul_shield.txt
+++ b/doc/src/pair_coul_shield.txt
@@ -38,7 +38,8 @@ charge and molecule ID information is included.
Where Tap(r_ij) is the taper function which provides a continuous cutoff
(up to third derivative) for inter-atomic separations larger than r_c
-"(Maaravi)"_#Maaravi1. Here {lambda} is the shielding parameter that
+"(Leven1)"_#Leven3, "(Leven2)"_#Leven4 and "(Maaravi)"_#Maaravi1.
+Here {lambda} is the shielding parameter that
eliminates the short-range singularity of the classical mono-polar
electrostatic interaction expression "(Maaravi)"_#Maaravi1.
@@ -82,5 +83,11 @@ package"_Build_package.html doc page for more info.
:line
+:link(Leven3)
+[(Leven1)] I. Leven, I. Azuri, L. Kronik and O. Hod, J. Chem. Phys. 140, 104106 (2014).
+
+:link(Leven4)
+[(Leven2)] I. Leven et al, J. Chem.Theory Comput. 12, 2896-905 (2016).
+
:link(Maaravi1)
[(Maaravi)] T. Maaravi et al, J. Phys. Chem. C 121, 22826-22835 (2017).
diff --git a/doc/src/pair_ilp_graphene_hbn.txt b/doc/src/pair_ilp_graphene_hbn.txt
index c74028faa91fe9ff80ef42560c13d1bf74b3c685..f048b16ccf0dad6828c9f52330e1f790ee694c28 100644
--- a/doc/src/pair_ilp_graphene_hbn.txt
+++ b/doc/src/pair_ilp_graphene_hbn.txt
@@ -25,22 +25,24 @@ 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
pair_coeff 1 1 coul/shield 0.70
-pair_coeff 1 2 coul/shield 0.69498201415576216335
+pair_coeff 1 2 coul/shield 0.695
pair_coeff 2 2 coul/shield 0.69 :pre
[Description:]
The {ilp/graphene/hbn} style computes the registry-dependent interlayer
-potential (ILP) potential as described in "(Leven)"_#Leven and
-"(Maaravi)"_#Maaravi2. The normals are calculated in the way as described
+potential (ILP) potential as described in "(Leven1)"_#Leven1,
+"(Leven2)"_#Leven2 and "(Maaravi)"_#Maaravi2.
+The normals are calculated in the way as described
in "(Kolmogorov)"_#Kolmogorov2.
:c,image(Eqs/pair_ilp_graphene_hbn.jpg)
Where Tap(r_ij) is the taper function which provides a continuous
cutoff (up to third derivative) for interatomic separations larger than
-r_c "(Maaravi)"_#Maaravi2. The definitions of each parameter in the above
-equation can be found in "(Leven)"_#Leven and "(Maaravi)"_#Maaravi2.
+r_c "(Maaravi)"_#Maaravi2. The definitons of each parameter in the above
+equation can be found in "(Leven1)"_#Leven1 and "(Maaravi)"_#Maaravi2.
+
It is important to include all the pairs to build the neighbor list for
calculating the normals.
@@ -61,13 +63,15 @@ NOTE: The parameters presented in the parameter file (e.g. BNCH.ILP),
are fitted with taper function by setting the cutoff equal to 16.0
Angstrom. Using different cutoff or taper function should be careful.
-NOTE: Two new sets of parameters of ILP for two-dimensional hexagonal Materials
-are presented in "(Ouyang)"_#Ouyang1. These parameters provide a good description
-in both short- and long-range interaction regime, while the old ILP parameters
-published in "(Leven)"_#Leven and "(Maaravi)"_#Maaravi2 are only suitable for
-long-range interaction regime. This feature is essential for simulations in
-high-pressure regime (i.e., the interlayer distance smaller than the equilibrium distance).
-The benchmark tests and comparison of these parameters can be found in "(Ouyang)"_#Ouyang1.
+NOTE: Two new sets of parameters of ILP for two-dimensional hexagonal
+Materials are presented in "(Ouyang)"_#Ouyang. These parameters provide
+a good description in both short- and long-range interaction regimes.
+While the old ILP parameters published in "(Leven2)"_#Leven2 and
+"(Maaravi)"_#Maaravi2 are only suitable for long-range interaction
+regime. This feature is essential for simulations in high pressure
+regime (i.e., the interlayer distance is smaller than the equilibrium
+distance). The benchmark tests and comparison of these parameters can
+be found in "(Ouyang)"_#Ouyang.
This potential must be used in combination with hybrid/overlay.
Other interactions can be set to zero using pair_style {none}.
@@ -112,14 +116,17 @@ units, if your simulation does not use {metal} units.
:line
-:link(Leven)
-[(Leven)] I. Leven et al, J. Chem.Theory Comput. 12, 2896-905 (2016)
+:link(Leven1)
+[(Leven1)] I. Leven, I. Azuri, L. Kronik and O. Hod, J. Chem. Phys. 140, 104106 (2014).
+
+:link(Leven2)
+[(Leven2)] I. Leven et al, J. Chem.Theory Comput. 12, 2896-905 (2016).
:link(Maaravi2)
[(Maaravi)] T. Maaravi et al, J. Phys. Chem. C 121, 22826-22835 (2017).
:link(Kolmogorov2)
-[(Kolmogorov)] A. N. Kolmogorov, V. H. Crespi, Phys. Rev. B 71, 235415 (2005)
+[(Kolmogorov)] A. N. Kolmogorov, V. H. Crespi, Phys. Rev. B 71, 235415 (2005).
-:link(Ouyang1)
-[(Ouyang)] W. Ouyang, D. Mandelli, M. Urbakh, O. Hod, arXiv:1806.09555 (2018).
+:link(Ouyang)
+[(Ouyang)] W. Ouyang, D. Mandelli, M. Urbakh and O. Hod, Nano Lett. 18, 6009-6016 (2018).
diff --git a/doc/src/pair_kolmogorov_crespi_full.txt b/doc/src/pair_kolmogorov_crespi_full.txt
index 358d2168df277c2ce74e00992dfb87fe403817e5..df9a9696be4108d7ddaccddd5884405a3ce4a338 100644
--- a/doc/src/pair_kolmogorov_crespi_full.txt
+++ b/doc/src/pair_kolmogorov_crespi_full.txt
@@ -19,11 +19,11 @@ tap_flag = 0/1 to turn off/on the taper function
pair_style hybrid/overlay kolmogorov/crespi/full 20.0 0
pair_coeff * * none
-pair_coeff * * kolmogorov/crespi/full CC.KC C C :pre
+pair_coeff * * kolmogorov/crespi/full CH.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
+pair_style hybrid/overlay rebo kolmogorov/crespi/full 16.0 1
+pair_coeff * * rebo CH.airebo C H
+pair_coeff * * kolmogorov/crespi/full CH_taper.KC C H :pre
[Description:]
@@ -38,27 +38,32 @@ forces and to include all the pairs to build the neighbor list for
calculating the normals. Energies are shifted so that they go
continuously 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.
+the last term in the equation for {Vij} above. This is essential only
+when the tapper function is turned off. The formula of taper function
+can be found in pair style "ilp/graphene/hbn"_pair_ilp_graphene_hbn.html.
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}
+The parameter file (e.g. CH.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.
-NOTE: A new set of parameters of KC potential for hydrocarbons (CH.KC)
-is presented in "(Ouyang)"_#Ouyang2. The parameters in CH.KC provides
-a good description in both short- and long-range interaction regime,
-while the original parameters (CC.KC) published in "(Kolmogorov)"_#Kolmogorov1
-are only suitable for long-range interaction regime.
-This feature is essential for simulations in high-pressure regime
-(i.e., the interlayer distance smaller than the equilibrium distance).
-The benchmark tests and comparison of these parameters can be found in "(Ouyang)"_#Ouyang2.
+NOTE: Two new sets of parameters of KC potential for hydrocarbons, CH.KC
+(without the taper function) and CH_taper.KC (with the taper function)
+are presented in "(Ouyang)"_#Ouyang1. The energy for the KC potential
+with the taper function goes continuously to zero at the cutoff. The
+parameters in both CH.KC and CH_taper.KC provide a good description in
+both short- and long-range interaction regimes. While the original
+parameters (CC.KC) published in "(Kolmogorov)"_#Kolmogorov1 are only
+suitable for long-range interaction regime. This feature is essential
+for simulations in high pressure regime (i.e., the interlayer distance
+is smaller than the equilibrium distance). The benchmark tests and
+comparison of these parameters can be found in "(Ouyang)"_#Ouyang1.
This potential must be used in combination with hybrid/overlay.
Other interactions can be set to zero using pair_style {none}.
@@ -84,7 +89,7 @@ package"_Build_package.html doc page 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
+The CH.KC potential file provided with LAMMPS (see the potentials
folder) are parameterized for metal units. You can use this potential
with any LAMMPS units, but you would need to create your own custom
CC.KC potential file with all coefficients converted to the appropriate
@@ -105,5 +110,5 @@ units.
:link(Kolmogorov1)
[(Kolmogorov)] A. N. Kolmogorov, V. H. Crespi, Phys. Rev. B 71, 235415 (2005)
-:link(Ouyang2)
-[(Ouyang)] W. Ouyang, D. Mandelli, M. Urbakh, O. Hod, arXiv:1806.09555 (2018).
+:link(Ouyang1)
+[(Ouyang)] W. Ouyang, D. Mandelli, M. Urbakh and O. Hod, Nano Lett. 18, 6009-6016 (2018).
diff --git a/potentials/BNCH.ILP b/potentials/BNCH.ILP
index 3f38b5e35aea826c1732f61e68830cbedc874db0..6ca56ca18620a049ab2d5227f33a61039e7bced9 100644
--- a/potentials/BNCH.ILP
+++ b/potentials/BNCH.ILP
@@ -1,6 +1,6 @@
# Interlayer Potential (ILP) for graphene/graphene, graphene/hBN and hBN/hBN junctions
#
-# Cite as Wengen Ouyang,Davide Mandelli, Michael Urbakh, Oded Hod, arXiv:1806.09555 (2018).
+# Cite as W. Ouyang, D. Mandelli, M. Urbakh and O. Hod, Nano Letters 18, 6009-6016 (2018).
#
# ----------------- Repulsion Potential ------------------++++++++++++++ Vdw Potential ++++++++++++++++************
# beta(A) alpha delta(A) epsilon(meV) C(meV) d sR reff(A) C6(meV*A^6) S rcut
diff --git a/potentials/CH.KC b/potentials/CH.KC
index d2fe543cf8f11cd21581767edcb0ad0a1e93121f..029f682f3494725ce1b5f16bc2ea7f8fa72c0325 100644
--- a/potentials/CH.KC
+++ b/potentials/CH.KC
@@ -1,6 +1,6 @@
-# Refined parameters for Kolmogorov-Crespi Potential
+# Refined parameters for Kolmogorov-Crespi Potential without taper function
#
-# Cite as Wengen Ouyang,Davide Mandelli, Michael Urbakh, Oded Hod, arXiv:1806.09555 (2018).
+# Cite as W. Ouyang, D. Mandelli, M. Urbakh and O. Hod, Nano Letters 18, 6009-6016 (2018).
#
# z0 C0 C2 C4 C delta lambda A S rcut
C C 3.328819 21.847167 12.060173 4.711099 6.678908e-4 0.7718101 3.143921 12.660270 1.0 2.0
diff --git a/potentials/CH_taper.KC b/potentials/CH_taper.KC
new file mode 100644
index 0000000000000000000000000000000000000000..0f08b2e6dd3d7ad069f3b70aeb50dde7e7ae2c4d
--- /dev/null
+++ b/potentials/CH_taper.KC
@@ -0,0 +1,10 @@
+# Refined parameters for Kolmogorov-Crespi Potential with taper function
+#
+# Cite as W. Ouyang, D. Mandelli, M. Urbakh and O. Hod, Nano Letters 18, 6009-6016 (2018).
+#
+# z0 C0 C2 C4 C delta lambda A S rcut
+C C 3.416084 20.021583 10.9055107 4.2756354 1.0010836E-2 0.8447122 2.9360584 14.3132588 1.0 2.0
+C H 2.849054 72.557245 1.0164169E-2 65.923312 8.7962504E-5 0.3349237 3.0402632 14.7533201 1.0 1.5
+H H 2.187478 3.915802E-5 5.0896431E-5 3.6657827 1.5373722446 0.9633581 0.4249989 1.570737E-4 1.0 1.2
+H C 2.849054 72.557245 1.0164169E-2 65.923312 8.7962504E-5 0.3349237 3.0402632 14.7533201 1.0 2.2
+
diff --git a/src/USER-MISC/pair_ilp_graphene_hbn.cpp b/src/USER-MISC/pair_ilp_graphene_hbn.cpp
index 6b925954726f094bcc10e8d517e4c1c5e214dbec..f613f0aff058f70888dd12dd649f907f46202fec 100644
--- a/src/USER-MISC/pair_ilp_graphene_hbn.cpp
+++ b/src/USER-MISC/pair_ilp_graphene_hbn.cpp
@@ -249,11 +249,9 @@ void PairILPGrapheneHBN::compute(int eflag, int vflag)
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;
- }
+ 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];
@@ -274,15 +272,13 @@ void PairILPGrapheneHBN::compute(int eflag, int vflag)
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;
- }
+ // 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) {
@@ -729,7 +725,8 @@ void PairILPGrapheneHBN::ILP_neigh()
ILP_firstneigh[i] = neighptr;
ILP_numneigh[i] = n;
- if (n > 3) error->all(FLERR,"There are too many neighbors for some atoms, please reduce the cutoff for normals");
+ if (n == 0) error->all(FLERR,"Could not build neighbor list to calculate normals, please check your configuration");
+ if (n > 3) error->all(FLERR,"There are too many neighbors for some atoms, please check your configuration");
ipage->vgot(n);
if (ipage->status())
error->one(FLERR,"Neighbor list overflow, boost neigh_modify one");
diff --git a/src/USER-MISC/pair_kolmogorov_crespi_full.cpp b/src/USER-MISC/pair_kolmogorov_crespi_full.cpp
index 18870ef69a71605daba3fa6894c9ce4d42ef67f4..2f868ae4183a144f6846febd36149f16c461f9af 100644
--- a/src/USER-MISC/pair_kolmogorov_crespi_full.cpp
+++ b/src/USER-MISC/pair_kolmogorov_crespi_full.cpp
@@ -249,11 +249,9 @@ void PairKolmogorovCrespiFull::compute(int eflag, int vflag)
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;
- }
+ 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];
@@ -274,15 +272,13 @@ void PairKolmogorovCrespiFull::compute(int eflag, int vflag)
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;
- }
+ // 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) {
@@ -734,7 +730,8 @@ void PairKolmogorovCrespiFull::KC_neigh()
KC_firstneigh[i] = neighptr;
KC_numneigh[i] = n;
- if (n > 3) error->all(FLERR,"There are too many neighbors for some atoms, please reduce the cutoff for normals");
+ if (n == 0) error->all(FLERR,"Could not build neighbor list to calculate normals, please check your configuration");
+ if (n > 3) error->all(FLERR,"There are too many neighbors for some atoms, please check your configuration");
ipage->vgot(n);
if (ipage->status())
error->one(FLERR,"Neighbor list overflow, boost neigh_modify one");