Added documentation for compute_pair_entropy_atom

doc/src/Eqs/pair_entropy.tex

+\documentclass[12pt]{article}
+
+\begin{document}
+\thispagestyle{empty}
+
+$$
+   s_S^i=-2\pi\rho k_B \int\limits_0^{r_m} \left [ g(r) \ln g(r) - g(r) + 1 \right ] r^2 dr ,
+$$
+
+\end{document}

doc/src/Eqs/pair_entropy2.tex

+\documentclass[12pt]{article}
+
+\begin{document}
+\thispagestyle{empty}
+
+$$
+   g_m^i(r) = \frac{1}{4 \pi \rho r^2} \sum\limits_{j} \frac{1}{\sqrt{2 \pi \sigma^2}} e^{-(r-r_{ij})^2/(2\sigma^2)} ,
+$$
+
+\end{document}

doc/src/Eqs/pair_entropy3.tex

+\documentclass[12pt]{article}
+
+\begin{document}
+\thispagestyle{empty}
+
+$$
+  \bar{s}_S^i  = \frac{\sum_j s_S^j + s_S^i}{N + 1} ,
+$$
+
+\end{document}

doc/src/compute_pair_entropy_atom.txt

+"LAMMPS WWW Site"_lws - "LAMMPS Documentation"_ld - "LAMMPS Commands"_lc :c
+
+:link(lws,http://lammps.sandia.gov)
+:link(ld,Manual.html)
+:link(lc,Section_commands.html#comm)
+
+:line
+
+compute pentropy/atom command :h3
+
+[Syntax:]
+
+compute ID group-ID pentropy/atom sigma cutoff  :l
+compute ID group-ID pentropy/atom sigma cutoff avg yes/no cutoff2 ... :pre
+
+ID, group-ID are documented in "compute"_compute.html command :l
+pentropy/atom = style name of this compute command :l
+sigma = width of gaussians used in the g(r) smoothening :l
+cutoff = cutoff for the g(r) calculation :l
+one or more keyword/value pairs may be appended :l
+{avg} yes/no cutoff2
+  avg = {no} or {yes}
+    {no} = do not average the pair entropy over neighbors
+    {yes} = average the pair entropy over neighbors
+  cutoff2 = cutoff for the averaging over neighbors :pre
+:ule
+
+[Examples:]
+
+compute 1 all pentropy/atom 0.25 5. :pre
+compute 1 all pentropy/atom 0.25 5. avg yes 5. :pre
+
+[Description:]
+
+Define a computation that calculates the pair entropy fingerprint for
+each atom in the group. The fingerprint is useful to distinguish between
+ordered and disordered environments, for instance liquid and solid-like 
+environments, or glassy and crystalline-like environments. Some 
+ applications could be the identification of grain boundaries, a 
+ melt-solid interface, or a solid cluster emerging from the melt. 
+The advantage of this parameter over others is that no a priori 
+ information about the solid structure is required.
+
+This parameter for atom i is computed using the following formula from
+"(Piaggi)"_#Piaggi
+
+:c,image(Eqs/pair_entropy.jpg)
+
+where r is a distance, g(r) is the radial distribution function of atom 
+ i and rho is the density of the system. The g(r) computed for each 
+ atom i can be noisy and therefore it is smoothened using:
+
+:c,image(Eqs/pair_entropy2.jpg)
+
+where the sum in j goes through the neighbors of atom i, and sigma is a
+parameter to control the smoothening.
+
+The input parameters are {sigma} the smoothening parameter, and the
+  {cutoff} for the calculation of g(r). 
+
+If the keyword {avg} has the setting {yes}, then this compute also
+ averages the parameter over the neighbors  of atom i according to:
+
+:c,image(Eqs/pair_entropy3.jpg)
+
+where the sum j goes over the neighbors of atom i and N is the number
+ of neighbors. This procedure provides a sharper distinction between
+order and disorder environments. In this case the input parameter 
+ {cutoff2} is the cutoff for the averaging over the neighbors and 
+ must also be specified.
+
+If the {avg yes} option is used, the effective cutoff of the neighbor
+ list should be {cutoff}+{cutoff2} and therefore it might be necessary 
+ to increase the skin of the neighbor list with:
+
+neighbor skin bin :pre
+
+See "neighbor"_neighbor.html for details.
+
+The neighbor list needed to compute this quantity is constructed each
+time the calculation is performed (e.g. each time a snapshot of atoms
+is dumped).  Thus it can be inefficient to compute/dump this quantity
+too frequently or to have multiple compute/dump commands, each with a
+{centro/atom} style.
+
+[Output info:]
+
+By default, this compute calculates the pair entropy value for each
+atom as a per-atom vector, which can be accessed by any command that
+uses per-atom values from a compute as input.  See "Section
+6.15"_Section_howto.html#howto_15 for an overview of LAMMPS output
+options.
+
+The pair entropy values have units of the Boltzmann constant. They are 
+ always negative, and lower values (lower entropy) correspond to more
+ ordered environments.
+
+Here are typical input parameters for fcc aluminum (lattice 
+ constant 4.05 Angstroms),
+
+compute 1 all pentropy/atom 0.25 5.7 avg yes 3.7 :pre
+
+and for bcc sodium (lattice constant 4.23 Angstroms),
+
+compute 1 all pentropy/atom 0.25 7.3 avg yes 5.1 :pre
+
+[Restrictions:] none
+
+[Related commands:]
+
+"compute cna/atom"_compute_cna_atom.html
+"compute centro/atom"_compute_centro_atom.html
+
+[Default:]
+
+The default value for the optional keyword is avg = no.
+
+:line
+
+:link(Piaggi)
+[(Piaggi)] Piaggi and Parrinello, J Chem Phys, 147, 114112 (2017).

src/compute_pair_entropy_atom.cpp

@@ -47,6 +47,13 @@ ComputePairEntropyAtom::ComputePairEntropyAtom(LAMMPS *lmp, int narg, char **arg
   if (narg < 5 || narg > 8)
     error->all(FLERR,"Illegal compute pentropy/atom command");
 
+  // Arguments are: sigma cutoff avg yes/no cutoff2
+  //   sigma is the gaussian width 
+  //   cutoff is the cutoff for the calculation of g(r) 
+  //   avg is optional and it means averaginf the pair entropy over the neighbors
+  //   the next argument should be yes or no 
+  //   cutoff2 is the cutoff for the averaging
+
   sigma = force->numeric(FLERR,arg[3]);
   cutoff = force->numeric(FLERR,arg[4]);
   if (cutoff < 0.0) error->all(FLERR,"Illegal compute pentropy/atom command; negative cutoff");
@@ -74,7 +81,7 @@ ComputePairEntropyAtom::ComputePairEntropyAtom(LAMMPS *lmp, int narg, char **arg
   nbin = static_cast<int>(cutoff / sigma) + 1;
   nmax = 0;
   maxneigh = 0;
-  deltabin = 2;
+  deltabin = 2; // 2 seems a good compromise between speed and good mollification
   deltar = sigma; 
   peratom_flag = 1;
   size_peratom_cols = 0;
@@ -234,7 +241,7 @@ void ComputePairEntropyAtom::compute_peratom()
           maxbin=bin +  deltabin;
           if (maxbin > (nbin-1)) maxbin=nbin-1;
           for(int k=minbin;k<maxbin+1;k++) {
-            double invNormKernel=invNormConstantBase/rbinsq[bin];
+            double invNormKernel=invNormConstantBase/rbinsq[k];
             double distance = r - rbin[k];
             gofr[k] += invNormKernel*exp(-distance*distance/sigmasq2) ; 
           }