updated documentation and ev_tally

doc/src/Eqs/pair_kolmogorov_crespi_z.tex

@@ -5,9 +5,9 @@
 
 \begin{eqnarray*}
   E & = & \frac{1}{2} \sum_i \sum_{j \neq i} V_{ij} \\
-  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 & = & x_{ij}^2 + y_{ij}^2 \\
-  f(\rho) & = &  e^{-(\rho/\delta)^2} \sum_{n=0}^2 C_{2n} { \rho/\delta }^{2n}
+  V_{ij} & = & e^{-\lambda(r_{ij} -z_0}) \left[ C + f(\rho_{ij}) + f(\rho_{ji}) \right] - A \left( \frac{r_{ij}}{z_0}\right)^{-6} + A \left( \frac{\textrm{cutoff}}{z_0}\right)^{-6} \\
+  \rho_{ij}^2 = \rho_{ji}^2 & = &  x_{ij}^2 + y_{ij}^2 ~\hspace{2cm} (\mathbf{n_i}\equiv\hat \mathbf{z})\\
+  f(\rho) & = &  e^{-(\rho/\delta)^2} \sum_{n=0}^2 C_{2n} \left( \rho/\delta \right) ^{2n}
 \end{eqnarray*}
 
 \end{document}

doc/src/pair_kolmogorov_crespi_z.txt

@@ -10,7 +10,7 @@ pair_style kolmogorov/crespi/z command :h3
 
 [Syntax:]
 
-pair_style kolmogorov/crespi/z cutoff :pre
+pair_style hybrid/overlay kolmogorov/crespi/z cutoff :pre
 
 [Examples:]
 
@@ -24,14 +24,26 @@ pair_coeff 1 2 kolmogorov/crespi/z  CC.KC      C C :pre
 
 [Description:]
 
-The {kolmogorov/crespi/z} style computes the Kolmogorov-Crespi interaction potential as 
-described in "(KC)"_#KC. An important simplification is made, which is to 
-take all normals along the z-axis. 
+The {kolmogorov/crespi/z} style computes the Kolmogorov-Crespi interaction 
+potential as described in "(KC05)"_#KC05. An important simplification is made, 
+which is to take all normals along the z-axis. 
 
 :c,image(Eqs/pair_kolmogorov_crespi_z.jpg)
 
+It is important to have a suffiently large cutoff to ensure smooth forces. 
+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.
+
+This potential is intended for interactions between two layers of graphene. 
+Therefore, to avoid interaction between layers in multi-layered materials, 
+each layer should have a separate atom type and interactions should only 
+be computed between atom types of neighbouring layers.
+
 The parameter file (e.g. CC.KC), is intended for use with metal 
-"units"_units.html, with energies in meV.
+"units"_units.html, with energies in meV. An additional parameter, {S}, 
+is available to facilitate scaling of energies in accordance with 
+"(vanWijk)"_#vanWijk.
 
 This potential must be used in combination with hybrid/overlay.
 Other interactions can be set to zero using pair_style {none}.
@@ -46,10 +58,15 @@ LAMMPS"_Section_start.html#start_3 section for more info.
 
 "pair_coeff"_pair_coeff.html
 "pair_none"_pair_none.html
+"pair_style hybrid/overlay"_pair_hybrid.html
 
 [Default:] none
 
 :line
 
-:link(KC)
-[(KC)] Kolmogorov, Crespi, Phys. Rev. B 71, 235415 (2005)
+:link(KC05) 
+[(KC05)] A. N. Kolmogorov, V. H. Crespi, Phys. Rev. B 71, 235415 (2005)
+
+:link(vanWijk) 
+[(vanWijk)] M. M. van Wijk, A. Schuring, M. I. Katsnelson, and A. Fasolino, 
+Physical Review Letters, 113, 135504 (2014)

potentials/CC.KC

@@ -2,7 +2,7 @@
 #
 # Cite as A.N. Kolmogorov & V. H. Crespi, 
 # Registry-dependent interlayer potential for graphitic systems
-# Physical Review B 71, 2354150 (2005)
+# Physical Review B 71, 235415 (2005)
 # 
 #      z0     C0     C2      C4       C    delta  lambda      A     S
 C  C  3.34  15.71  12.29   4.933   3.030   0.578   3.629   10.238  1.0

src/USER-MISC/pair_kolmogorov_crespi_z.cpp

@@ -158,11 +158,12 @@ void PairKolmogorovCrespiZ::compute(int eflag, int vflag)
         }
 
         if (eflag) {
-          evdwl = -p.A*p.z06/r6+ exp1*sumCff+offset[itype][jtype];
+          evdwl = -p.A*p.z06/r6+ exp1*sumCff - offset[itype][jtype];
         }
 
         if (evflag){
-          ev_tally(i,j,nlocal,newton_pair,evdwl,0.0,fpair,delx,dely,delz);
+          ev_tally_xyz(i,j,nlocal,newton_pair,evdwl,0,
+                       fsum,fsum,fpair,delx,dely,delz);
         }
       }
     }
@@ -286,7 +287,8 @@ double PairKolmogorovCrespiZ::init_one(int i, int j)
 
   if (offset_flag) {
     int iparam_ij = elem2param[map[i]][map[j]];
-    offset[i][j] = params[iparam_ij].A*pow(params[iparam_ij].z0/cut[i][j],6);
+    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];