diff --git a/doc/src/compute_bond_local.rst b/doc/src/compute_bond_local.rst
index 0368f33d29..24b0943484 100644
--- a/doc/src/compute_bond_local.rst
+++ b/doc/src/compute_bond_local.rst
@@ -65,7 +65,7 @@ of some larger molecule.
 
 The value *dist* is the current length of the bond.
 The values *dx*, *dy*, and *dz* are the xyz components of the
-*distance* between the pair of atoms. This value is always the 
+*distance* between the pair of atoms. This value is always the
 distance from the atom of lower to the one with the higher id.
 
 The value *engpot* is the potential energy for the bond,
diff --git a/doc/src/compute_pair_local.rst b/doc/src/compute_pair_local.rst
index 5cff338758..38953d203c 100644
--- a/doc/src/compute_pair_local.rst
+++ b/doc/src/compute_pair_local.rst
@@ -58,7 +58,7 @@ commands.
 
 The value *dist* is the distance between the pair of atoms.
 The values *dx*, *dy*, and *dz* are the xyz components of the
-*distance* between the pair of atoms. This value is always the 
+*distance* between the pair of atoms. This value is always the
 distance from the atom of lower to the one with the higher id.
 
 The value *eng* is the interaction energy for the pair of atoms.
@@ -93,9 +93,9 @@ from the second of the two sub-styles.  If the referenced *pN*
 is not computed for the specific pairwise interaction (based on
 atom types), then the output will be 0.0.
 
-The value *dist*, *dx*, *dy* and *dz* will be in distance :doc:`units <units>`.  
-The value *eng* will be in energy :doc:`units <units>`.  
-The values *force*, *fx*, *fy*, and *fz* will be in force :doc:`units <units>`.  
+The value *dist*, *dx*, *dy* and *dz* will be in distance :doc:`units <units>`.
+The value *eng* will be in energy :doc:`units <units>`.
+The values *force*, *fx*, *fy*, and *fz* will be in force :doc:`units <units>`.
 The values *pN* will be in whatever units the pair style defines.
 
 The optional *cutoff* keyword determines how the force cutoff distance
diff --git a/doc/src/pair_hybrid.rst b/doc/src/pair_hybrid.rst
index 541cdc1911..1460927add 100644
--- a/doc/src/pair_hybrid.rst
+++ b/doc/src/pair_hybrid.rst
@@ -74,14 +74,17 @@ atoms interact with each other via an *eam* potential, the surface atoms
 interact with each other via a *lj/cut* potential, and the metal/surface
 interaction is also computed via a *lj/cut* potential.  The
 *hybrid/overlay* style could be used as in the second example above,
-where multiple potentials are superposed in an additive fashion to
+where multiple potentials are superimposed in an additive fashion to
 compute the interaction between atoms.  In this example, using *lj/cut*
 and *coul/long* together gives the same result as if the
 *lj/cut/coul/long* potential were used by itself.  In this case, it
 would be more efficient to use the single combined potential, but in
 general any combination of pair potentials can be used together in to
 produce an interaction that is not encoded in any single pair_style
-file, e.g. adding Coulombic forces between granular particles.
+file, e.g. adding Coulombic forces between granular particles.  Another
+limitation of using the *hybrid/overlay* variant, that it does not generate
+*lj/cut* parameters for mixed atom types from a mixing rule due to
+restrictions discussed below.
 
 If the *hybrid/scaled* style is used instead of *hybrid/overlay*,
 contributions from sub-styles are weighted by their scale factors, which
@@ -150,10 +153,14 @@ with Tersoff, and the cross-interactions with Lennard-Jones:
    pair_coeff * * tersoff 2 C.tersoff NULL C
    pair_coeff 1 2 lj/cut 1.0 1.5
 
-If pair coefficients are specified in the data file read via the
-:doc:`read_data <read_data>` command, then the same rule applies.
-E.g. "eam/alloy" or "lj/cut" must be added after the atom type, for
-each line in the "Pair Coeffs" section, e.g.
+
+It is not recommended to read pair coefficients for a hybrid style from a "Pair Coeffs"
+or "PairIJ Coeffs" section of a data file via the :doc:`read_data <read_data>` command,
+since those sections expect a fixed number of lines, either one line per atom type or
+one line pair pair of atom types, respectively.  When reading from a data file, the
+lines of the "Pair Coeffs" and "PairIJ Coeffs" are changed in the same way as the *pair_coeff*
+command, i.e. the name of the pair style to which the parameters apply must follow the
+atom type (or atom types), e.g.
 
 .. parsed-literal::
 
@@ -162,6 +169,11 @@ each line in the "Pair Coeffs" section, e.g.
    1 lj/cut/coul/cut 1.0 1.0
    ...
 
+   PairIJ Coeffs
+
+   1 1 lj/cut/coul/cut 1.0 1.0
+   ...
+
 Note that the pair_coeff command for some potentials such as
 :doc:`pair_style eam/alloy <pair_eam>` includes a mapping specification
 of elements to all atom types, which in the hybrid case, can include
@@ -208,12 +220,22 @@ examples above, or in the data file read by the :doc:`read_data
 <read_data>`, or by mixing as described below.  Also all sub-styles
 must be used at least once in a :doc:`pair_coeff <pair_coeff>` command.
 
-.. note::
+.. warning::
 
-   LAMMPS never performs mixing of parameters from different sub-styles,
-   **even** if they use the same type of coefficients, e.g. contain
-   a Lennard-Jones potential variant.  Those parameters must be provided
-   explicitly.
+   With hybrid pair styles the use of mixing to generate pair
+   coefficients is significantly limited compared to the individual pair
+   styles.  LAMMPS **never** performs mixing of parameters from
+   different sub-styles, **even** if they use the same type of
+   coefficients, e.g. contain a Lennard-Jones potential variant.  Those
+   parameters must be provided explicitly.  Also for *hybrid/overlay*
+   and *hybrid/scaled* mixing is **only** performed for pairs of atom
+   types for which only a single pair style is assigned.
+
+   Thus it is strongly recommended to provide all mixed terms
+   explicitly.  For non-hybrid styles those could be generated and
+   written out using the :doc:`write_coeff command <write_coeff>` and
+   then edited as needed to comply with the requirements for hybrid
+   styles as explained above.
 
 If you want there to be no interactions between a particular pair of
 atom types, you have 3 choices.  You can assign the pair of atom types
diff --git a/doc/src/pair_modify.rst b/doc/src/pair_modify.rst
index 1a62a4c1a0..4941693fbd 100644
--- a/doc/src/pair_modify.rst
+++ b/doc/src/pair_modify.rst
@@ -71,21 +71,23 @@ The *mix* keyword affects pair coefficients for interactions between
 atoms of type I and J, when I != J and the coefficients are not
 explicitly set in the input script.  Note that coefficients for I = J
 must be set explicitly, either in the input script via the
-:doc:`pair_coeff <pair_coeff>` command or in the "Pair Coeffs" section of the
-:doc:`data file <read_data>`.  For some pair styles it is not
+:doc:`pair_coeff <pair_coeff>` command or in the "Pair Coeffs" or "PairIJ Coeffs"
+sections of the :doc:`data file <read_data>`.  For some pair styles it is not
 necessary to specify coefficients when I != J, since a "mixing" rule
 will create them from the I,I and J,J settings.  The pair_modify
 *mix* value determines what formulas are used to compute the mixed
 coefficients.  In each case, the cutoff distance is mixed the same way
 as sigma.
 
-Note that not all pair styles support mixing and some mix options
-are not available for certain pair styles. Also, there are additional
-restrictions when using :doc:`pair style hybrid or hybrid/overlay <pair_hybrid>`.
-See the page for individual pair styles for those restrictions.  Note also that the
-:doc:`pair_coeff <pair_coeff>` command also can be used to directly set
-coefficients for a specific I != J pairing, in which case no mixing is
-performed.
+Note that not all pair styles support mixing and some mix options are
+not available for certain pair styles. Also, there are additional
+restrictions when using :doc:`pair style hybrid or hybrid/overlay
+<pair_hybrid>`.  See the page for individual pair styles for those
+restrictions.  Note also that the :doc:`pair_coeff <pair_coeff>` command
+also can be used to directly set coefficients for a specific I != J
+pairing, in which case no mixing is performed.  If possible, LAMMPS will
+print an informational message about how many of the mixed pair
+coefficients were generated and which mixing rule was applied.
 
 - mix *geometric*
 
diff --git a/src/pair.cpp b/src/pair.cpp
index 1039875718..0bc698f27e 100644
--- a/src/pair.cpp
+++ b/src/pair.cpp
@@ -41,6 +41,7 @@ using namespace LAMMPS_NS;
 using namespace MathConst;
 
 enum{NONE,RLINEAR,RSQ,BMP};
+static const std::string mixing_rule_names[Pair::SIXTHPOWER+1] = {"geometric", "arithmetic", "sixthpower" };
 
 // allocate space for static class instance variable and initialize it
 
@@ -217,11 +218,9 @@ void Pair::init()
   if (tail_flag && domain->nonperiodic && comm->me == 0)
     error->warning(FLERR,"Using pair tail corrections with non-periodic system");
   if (!compute_flag && tail_flag && comm->me == 0)
-    error->warning(FLERR,"Using pair tail corrections with "
-                   "pair_modify compute no");
+    error->warning(FLERR,"Using pair tail corrections with pair_modify compute no");
   if (!compute_flag && offset_flag && comm->me == 0)
-    error->warning(FLERR,"Using pair potential shift with "
-                   "pair_modify compute no");
+    error->warning(FLERR,"Using pair potential shift with pair_modify compute no");
 
   // for manybody potentials
   // check if bonded exclusions could invalidate the neighbor list
@@ -259,13 +258,18 @@ void Pair::init()
   etail = ptail = 0.0;
   mixed_flag = 1;
   double cut;
+  int mixed_count = 0;
 
   for (i = 1; i <= atom->ntypes; i++)
     for (j = i; j <= atom->ntypes; j++) {
-      if ((i != j) && setflag[i][j]) mixed_flag = 0;
+      did_mix = false;
       cut = init_one(i,j);
       cutsq[i][j] = cutsq[j][i] = cut*cut;
       cutforce = MAX(cutforce,cut);
+      if (i != j) {
+        if (setflag[i][j]) mixed_flag = 0;
+        if (did_mix) ++mixed_count;
+      }
       if (tail_flag) {
         etail += etail_ij;
         ptail += ptail_ij;
@@ -275,6 +279,12 @@ void Pair::init()
         }
       }
     }
+
+  if (!manybody_flag && (comm->me == 0)) {
+    const int num_mixed_pairs = atom->ntypes * (atom->ntypes - 1) / 2;
+    utils::logmesg(lmp,"  generated {} of {} mixed pair_coeff terms from {} mixing rule\n",
+                   mixed_count, num_mixed_pairs, mixing_rule_names[mix_flag]);
+  }
 }
 
 /* ----------------------------------------------------------------------
@@ -681,6 +691,7 @@ void Pair::free_disp_tables()
 
 double Pair::mix_energy(double eps1, double eps2, double sig1, double sig2)
 {
+  did_mix = true;
   if (mix_flag == GEOMETRIC)
     return sqrt(eps1*eps2);
   else if (mix_flag == ARITHMETIC)
@@ -688,7 +699,8 @@ double Pair::mix_energy(double eps1, double eps2, double sig1, double sig2)
   else if (mix_flag == SIXTHPOWER)
     return (2.0 * sqrt(eps1*eps2) *
       pow(sig1,3.0) * pow(sig2,3.0) / (pow(sig1,6.0) + pow(sig2,6.0)));
-  else return 0.0;
+  else did_mix = false;
+  return 0.0;
 }
 
 /* ----------------------------------------------------------------------
diff --git a/src/pair.h b/src/pair.h
index f37c0732ed..00e6734773 100644
--- a/src/pair.h
+++ b/src/pair.h
@@ -110,6 +110,7 @@ class Pair : protected Pointers {
                        //       public so external driver can check
   int compute_flag;    // 0 if skip compute()
   int mixed_flag;      // 1 if all itype != jtype coeffs are from mixing
+  bool did_mix;        // set to true by mix_energy() to indicate that mixing was performed
 
   enum { GEOMETRIC, ARITHMETIC, SIXTHPOWER };    // mixing options
 
diff --git a/src/pair_hybrid.cpp b/src/pair_hybrid.cpp
index b5daa111da..e962e02c9e 100644
--- a/src/pair_hybrid.cpp
+++ b/src/pair_hybrid.cpp
@@ -706,11 +706,10 @@ double PairHybrid::init_one(int i, int j)
   for (int k = 0; k < nmap[i][j]; k++) {
     map[j][i][k] = map[i][j][k];
     double cut = styles[map[i][j][k]]->init_one(i,j);
-    styles[map[i][j][k]]->cutsq[i][j] =
-      styles[map[i][j][k]]->cutsq[j][i] = cut*cut;
+    if (styles[map[i][j][k]]->did_mix) did_mix = true;
+    styles[map[i][j][k]]->cutsq[i][j] = styles[map[i][j][k]]->cutsq[j][i] = cut*cut;
     if (styles[map[i][j][k]]->ghostneigh)
-      cutghost[i][j] = cutghost[j][i] =
-        MAX(cutghost[i][j],styles[map[i][j][k]]->cutghost[i][j]);
+      cutghost[i][j] = cutghost[j][i] = MAX(cutghost[i][j],styles[map[i][j][k]]->cutghost[i][j]);
     if (tail_flag) {
       etail_ij += styles[map[i][j][k]]->etail_ij;
       ptail_ij += styles[map[i][j][k]]->ptail_ij;