diff --git a/doc/src/pair_python.rst b/doc/src/pair_python.rst
index 65e1cd1611..2c511a16c9 100644
--- a/doc/src/pair_python.rst
+++ b/doc/src/pair_python.rst
@@ -20,6 +20,9 @@ Examples
    pair_style python 2.5
    pair_coeff * * py_pot.LJCutMelt lj
 
+   pair_style python 10.0
+   pair_coeff * * py_pot.HarmonicCut A B
+
    pair_style hybrid/overlay coul/long 12.0 python 12.0
    pair_coeff * * coul/long
    pair_coeff * * python py_pot.LJCutSPCE OW NULL
@@ -159,9 +162,66 @@ Following the *LJCutMelt* example, here are the two functions:
    to be  multiplied by delta x, delta y, and delta z to conveniently obtain
    the three components of the force vector between these two atoms.
 
+Below is a more complex example using *real* units and defines an interaction
+equivalent to:
+
+.. code-block:: LAMMPS
+
+   units real
+   pair_style harmonic/cut
+   pair_coeff 1 1 0.2 9.0
+   pair_coeff 2 2 0.4 9.0
+
+This uses the default geometric mixing.  The equivalent input with pair
+style *python* is:
+
+.. code-block:: LAMMPS
+
+   units real
+   pair_style python 10.0
+   pair_coeff * * py_pot.Harmonic A B
+
+Note that while for pair style *harmonic/cut* the cutoff is implicitly
+set to the minimum of the harmonic potential, for pair style python a
+global cutoff must be set and it must be equal or larger to the implicit
+cutoff of the potential in python, which has to explicitly return zero
+force and energy beyond the cutoff.  Also, the mixed parameters have to
+be explicitly provided.  The corresponding python code is:
+
+.. code-block:: python
+
+   class Harmonic(LAMMPSPairPotential):
+       def __init__(self):
+           super(Harmonic,self).__init__()
+           self.units = 'real'
+           # set coeffs: K, r0
+           self.coeff = {'A'  : {'A'  : (0.2,9.0),
+                                 'B'  : (math.sqrt(0.2*0.4),9.0)},
+                         'B'  : {'A'  : (math.sqrt(0.2*0.4),9.0),
+                                 'B'  : (0.4,9.0)}}
+
+       def compute_force(self,rsq,itype,jtype):
+           coeff = self.coeff[self.pmap[itype]][self.pmap[jtype]]
+           r = math.sqrt(rsq)
+           delta = coeff[1]-r
+           if (r < coeff[1]):
+               return 2.0*delta*coeff[0]/r
+           else:
+               return 0.0
+
+       def compute_energy(self,rsq,itype,jtype):
+           coeff = self.coeff[self.pmap[itype]][self.pmap[jtype]]
+           r = math.sqrt(rsq)
+           delta = coeff[1]-r
+           if (r < coeff[1]):
+               return delta*delta*coeff[0]
+           else:
+               return 0.0
+
 ----------
 
-.. note::
+.. admonition:: Performance Impact
+   :class: note
 
    The evaluation of scripted python code will slow down the
    computation pairwise interactions quite significantly. However, this