diff --git a/doc/src/Eqs/pair_entropy.jpg b/doc/src/Eqs/pair_entropy.jpg
new file mode 100644
index 0000000000000000000000000000000000000000..44161e6aa0e74a05544f8f29aab3ca942ef305a5
Binary files /dev/null and b/doc/src/Eqs/pair_entropy.jpg differ
diff --git a/doc/src/Eqs/pair_entropy.tex b/doc/src/Eqs/pair_entropy.tex
new file mode 100644
index 0000000000000000000000000000000000000000..304c9d6138500061a55df40a6d8fb9b7863b45d7
--- /dev/null
+++ b/doc/src/Eqs/pair_entropy.tex
@@ -0,0 +1,10 @@
+\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}
diff --git a/doc/src/Eqs/pair_entropy2.jpg b/doc/src/Eqs/pair_entropy2.jpg
new file mode 100644
index 0000000000000000000000000000000000000000..db3dcb0ca75f444a6640e704c61326d7c8c598a9
Binary files /dev/null and b/doc/src/Eqs/pair_entropy2.jpg differ
diff --git a/doc/src/Eqs/pair_entropy2.tex b/doc/src/Eqs/pair_entropy2.tex
new file mode 100644
index 0000000000000000000000000000000000000000..c91004840fee1a3998702a2ce61cd14c85048b71
--- /dev/null
+++ b/doc/src/Eqs/pair_entropy2.tex
@@ -0,0 +1,10 @@
+\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}
diff --git a/doc/src/Eqs/pair_entropy3.jpg b/doc/src/Eqs/pair_entropy3.jpg
new file mode 100644
index 0000000000000000000000000000000000000000..0d5e0025b7200f9ca4255c48bd2bbabb291553a2
Binary files /dev/null and b/doc/src/Eqs/pair_entropy3.jpg differ
diff --git a/doc/src/Eqs/pair_entropy3.tex b/doc/src/Eqs/pair_entropy3.tex
new file mode 100644
index 0000000000000000000000000000000000000000..b0b3c763585dc98de3a66cde754a2bcd90db91f3
--- /dev/null
+++ b/doc/src/Eqs/pair_entropy3.tex
@@ -0,0 +1,10 @@
+\documentclass[12pt]{article}
+
+\begin{document}
+\thispagestyle{empty}
+
+$$
+ \bar{s}_S^i = \frac{\sum_j s_S^j + s_S^i}{N + 1} ,
+$$
+
+\end{document}
diff --git a/doc/src/compute_pair_entropy_atom.txt b/doc/src/compute_pair_entropy_atom.txt
new file mode 100644
index 0000000000000000000000000000000000000000..bed36ed4b681e29b8f0a533ad1d5864ae5e840c4
--- /dev/null
+++ b/doc/src/compute_pair_entropy_atom.txt
@@ -0,0 +1,121 @@
+"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).
diff --git a/src/compute_pair_entropy_atom.cpp b/src/compute_pair_entropy_atom.cpp
index c4e411f761fef434f2afc4a0143f2dd418eb14b6..97461b4eec837b404c312efa33286c114df092d8 100644
--- a/src/compute_pair_entropy_atom.cpp
+++ b/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) ;
}