small modifications of potential #839

doc/src/pair_coul_shield.txt

@@ -38,7 +38,7 @@ 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)"_#Leven2, "(Leven2)"_#Leven3 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 +82,11 @@ LAMMPS"_Section_start.html#start_2_3 section for more info.
 
 :line
 
+:link(Leven2)
+[(Leven)] I. Leven, I. Azuri, L. Kronik and O. Hod, J. Chem. Phys. 140, 104106 (2014).
+
+:link(Leven3)
+[(Leven)] 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).

doc/src/pair_ilp_graphene_hbn.txt

@@ -31,7 +31,7 @@ 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
+potential (ILP) potential as described in "(Leven1)"_#Leven1, "(Leven2)"_#Leven and
 "(Maaravi)"_#Maaravi2. The normals are calculated in the way as described
 in "(Kolmogorov)"_#Kolmogorov2.
 
@@ -40,7 +40,7 @@ in "(Kolmogorov)"_#Kolmogorov2.
 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 definitons of each parameter in the above
-equation can be found in "(Leven)"_#Leven and "(Maaravi)"_#Maaravi2.
+equation can be found in "(Leven)"_#Leven1 and "(Maaravi)"_#Maaravi2.
 
 It is important to include all the pairs to build the neighbor list for
 calculating the normals.
@@ -62,9 +62,10 @@ 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)"_#Ouyang.
-These parameters provide a good description in both short- and long-range interaction regime. 
+These parameters provide a good description in both short- and long-range interaction regimes. 
-While the old ILP parameters published in "(Leven)"_#Leven and "(Maaravi)"_#Maaravi2 are only suitable for long-range interaction regime.
+While the old ILP parameters published in "(Leven)"_#Leven and "(Maaravi)"_#Maaravi2 are
-This feature is essential for simulations in high pressure regime (i.e., interlayer distance smaller than the equilibrium distance).
+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.
@@ -110,14 +111,17 @@ units, if your simulation does not use {metal} units.
 
 :line
 
+:link(Leven1)
+[(Leven)] I. Leven, I. Azuri, L. Kronik and O. Hod, J. Chem. Phys. 140, 104106 (2014).
+
 :link(Leven)
-[(Leven)] I. Leven et al, J. Chem.Theory Comput. 12, 2896-905 (2016)
+[(Leven)] 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(Ouyang)
-[(Ouyang)] W. Ouyang, D. Mandelli,  M. Urbakh, O. Hod, arXiv:1806.09555 (2018).
+[(Ouyang)] W. Ouyang, D. Mandelli, M. Urbakh and O. Hod, Nano Lett. 18, 6009-6016 (2018).

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_style hybrid/overlay rebo kolmogorov/crespi/full 16.0 1
-pair_coeff * * rebo                    CH.airebo  C C
+pair_coeff * * rebo                    CH.airebo    C H
-pair_coeff * * kolmogorov/crespi/full  CC.KC      C C :pre
+pair_coeff * * kolmogorov/crespi/full  CH_taper.KC  C H :pre
 
 [Description:]
 
@@ -38,24 +38,30 @@ 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)"_#Ouyang.
+NOTE: Two new sets of parameters of KC potential for hydrocarbons, CH.KC (without the taper function)
-The parameters in CH.KC provides a good description in both short- and long-range interaction regime.
+and CH_taper.KC (with the taper function) are presented in "(Ouyang)"_#Ouyang1.
-While the original parameters (CC.KC) published in "(Kolmogorov)"_#Kolmogorov1 are only suitable for long-range interaction regime. 
+The energy for the KC potential with the taper function goes continuously to zero at the cutoff. 
-This feature is essential for simulations in high pressure regime (i.e., interlayer distance smaller than the equilibrium distance).
+The parameters in both CH.KC and CH_taper.KC provide a good description in
-The benchmark tests and comparison of these parameters can be found in "(Ouyang)"_#Ouyang.
+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}.
@@ -81,7 +87,7 @@ 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
+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
@@ -102,5 +108,5 @@ units.
 :link(Kolmogorov1)
 [(Kolmogorov)] A. N. Kolmogorov, V. H. Crespi, Phys. Rev. B 71, 235415 (2005)
 
-:link(Ouyang)
+:link(Ouyang1)
-[(Ouyang)] W. Ouyang, D. Mandelli,  M. Urbakh, O. Hod, arXiv:1806.09555 (2018).
+[(Ouyang)] W. Ouyang, D. Mandelli, M. Urbakh and O. Hod, Nano Lett. 18, 6009-6016 (2018).

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

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

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
+

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][0] -= fkcx + fprod2[0]*Tap;
+        f[j][1] -= fkcy + fprod2[1]*Tap;
-          f[j][1] -= fkcy + fprod2[1]*Tap;
+        f[j][2] -= fkcz + fprod2[2]*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
-            // 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[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[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;
-            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][0] += (-prodnorm2*dprodnorm2[0]*fpair2)*Tap;
+          f[l][1] += (-prodnorm2*dprodnorm2[1]*fpair2)*Tap;
-            f[l][1] += (-prodnorm2*dprodnorm2[1]*fpair2)*Tap;
+          f[l][2] += (-prodnorm2*dprodnorm2[2]*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");

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][0] -= fkcx + fprod2[0]*Tap;
+        f[j][1] -= fkcy + fprod2[1]*Tap;
-          f[j][1] -= fkcy + fprod2[1]*Tap;
+        f[j][2] -= fkcz + fprod2[2]*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
-            // 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[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[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;
-            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][0] += (-prodnorm2*dprodnorm2[0]*fpair2)*Tap;
+	  f[l][1] += (-prodnorm2*dprodnorm2[1]*fpair2)*Tap;
-	    f[l][1] += (-prodnorm2*dprodnorm2[1]*fpair2)*Tap;
+	  f[l][2] += (-prodnorm2*dprodnorm2[2]*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");