diff --git a/doc/pair_comb.html b/doc/pair_comb.html
index a09ca8732..38e59ffd2 100644
--- a/doc/pair_comb.html
+++ b/doc/pair_comb.html
@@ -1,198 +1,199 @@
 <HTML>
 <CENTER><A HREF = "http://lammps.sandia.gov">LAMMPS WWW Site</A> - <A HREF = "Manual.html">LAMMPS Documentation</A> - <A HREF = "Section_commands.html#comm">LAMMPS Commands</A> 
 </CENTER>
 
 
 
 
 
 
 <HR>
 
 <H3>pair_style comb command 
 </H3>
 <H3>pair_style comb/omp command 
 </H3>
 <H3>pair_style comb3 command 
 </H3>
 <P><B>Syntax:</B>
 </P>
 <PRE>pair_style comb
 pair_style comb3 keyword 
 </PRE>
 <PRE>keyword = <I>polar</I>
   <I>polar</I> value = <I>polar_on</I> or <I>polar_off</I> = whether or not to include atomic polarization 
 </PRE>
 
 </UL>
 <P><B>Examples:</B>
 </P>
 <PRE>pair_style comb
 pair_coeff * * ../potentials/ffield.comb Si
 pair_coeff * * ../potentials/ffield.comb Hf Si O 
 </PRE>
 <PRE>pair_style comb3 polar_off
 pair_coeff * * ../potentials/ffield.comb3 O Cu N C O 
 </PRE>
 <P><B>Description:</B>
 </P>
 <P>Style <I>comb</I> computes the second-generation variable charge COMB
 (Charge-Optimized Many-Body) potential.  Style <I>comb3</I> computes the
 third-generation COMB potential.  These COMB potentials are described
 in <A HREF = "#COMB">(COMB)</A> and <A HREF = "#COMB3">(COMB3)</A>.  Briefly, the total energy
 <I>E<sub>T</sub></I> of a system of atoms is given by
 </P>
 <CENTER><IMG SRC = "Eqs/pair_comb1.jpg">
 </CENTER>
 <P>where <I>E<sub>i</sub><sup>self</sup></I> is the self-energy of atom <I>i</I> 
 (including atomic ionization energies and electron affinities),
 <I>E<sub>ij</sub><sup>short</sup></I> is the bond-order potential between 
 atoms <I>i</I> and <I>j</I>,
 <I>E<sub>ij</sub><sup>Coul</sup></I> is the Coulomb interactions,
 <I>E<sup>polar</sup></I> is the polarization term for organic systems 
 (style <I>comb3</I> only), 
 <I>E<sup>vdW</sup></I> is the van der Waals energy (style <I>comb3</I> only), 
 <I>E<sup>barr</sup></I> is a charge barrier function, and 
 <I>E<sup>corr</sup></I> are angular correction terms.
 </P>
 <P>The COMB potentials (styles <I>comb</I> and <I>comb3</I>) are variable charge
 potentials.  The equilibrium charge on each atom is calculated by the
 electronegativity equalization (QEq) method.  See <A HREF = "#Rick">Rick</A> for
 further details.  This is implemented by the <A HREF = "fix_qeq_comb.html">fix
 qeq/comb</A> command, which should normally be
 specified in the input script when running a model with the COMB
 potential.  The <A HREF = "fix_qeq_comb.html">fix qeq/comb</A> command has options
 that determine how often charge equilibration is performed, its
 convergence criterion, and which atoms are included in the
 calculation.
 </P>
 <P>Only a single pair_coeff command is used with the <I>comb</I> and <I>comb3</I>
 styles which specifies the COMB potential file with parameters for all
 needed elements.  These are mapped to LAMMPS atom types by specifying
 N additional arguments after the potential file in the pair_coeff
 command, where N is the number of LAMMPS atom types.
 </P>
 <P>For example, if your LAMMPS simulation of a Si/SiO<sub>2</sub>/
 HfO<sub>2</sub> interface has 4 atom types, and you want the 1st and
 last to be Si, the 2nd to be Hf, and the 3rd to be O, and you would
 use the following pair_coeff command:
 </P>
 <PRE>pair_coeff * * ../potentials/ffield.comb Si Hf O Si 
 </PRE>
 <P>The first two arguments must be * * so as to span all LAMMPS atom
 types.  The first and last Si arguments map LAMMPS atom types 1 and 4
 to the Si element in the <I>ffield.comb</I> file.  The second Hf argument
 maps LAMMPS atom type 2 to the Hf element, and the third O argument
 maps LAMMPS atom type 3 to the O element in the potential file.  If a
 mapping value is specified as NULL, the mapping is not performed.
 This can be used when a <I>comb</I> potential is used as part of the
 <I>hybrid</I> pair style.  The NULL values are placeholders for atom types
 that will be used with other potentials.
 </P>
 <P>For style <I>comb</I>, the provided potential file <I>ffield.comb</I> contains
 all currently-available 2nd generation COMB parameterizations: for Si,
 Cu, Hf, Ti, O, their oxides and Zr, Zn and U metals.  For style
 <I>comb3</I>, the potential file <I>ffield.comb3</I> contains all
 currently-available 3rd generation COMB paramterizations: O, Cu, N, C,
-H, Ti and Zn.  The status of the optimization of the compounds, for
+H, Ti, Zn and Zr.  The status of the optimization of the compounds, for
 example Cu<sub>2</sub>O, TiN and hydrocarbons, are given in the
 following table:
 </P>
 <CENTER><IMG SRC = "Eqs/pair_comb2.jpg">
 </CENTER>
 <P>For style <I>comb3</I>, in addition to ffield.comb3, a special parameter
 file, <I>lib.comb3</I>, that is exclusively used for C/O/H systems, will be
 automatically loaded if carbon atom is detected in LAMMPS input
 structure.  Keyword <I>polar</I> indicates whether the force field includes
 the atomic polarization. Since the equilibration of the polarization
 has not yet been implemented, it can only set polar_off at present.
 </P>
 <P>IMPORTANT NOTE: You can not use potential file <I>ffield.comb</I> with
 style <I>comb3</I>, nor file <I>ffield.comb3</I> with style <I>comb</I>.
 </P>
 <HR>
 
 <P>Styles with a <I>cuda</I>, <I>gpu</I>, <I>omp</I>, or <I>opt</I> suffix are functionally
 the same as the corresponding style without the suffix.  They have
 been optimized to run faster, depending on your available hardware, as
 discussed in <A HREF = "Section_accelerate.html">Section_accelerate</A> of the
 manual.  The accelerated styles take the same arguments and should
 produce the same results, except for round-off and precision issues.
 </P>
 <P>These accelerated styles are part of the USER-CUDA, GPU, USER-OMP and OPT
 packages, respectively.  They are only enabled if LAMMPS was built with
 those packages.  See the <A HREF = "Section_start.html#start_3">Making LAMMPS</A>
 section for more info.
 </P>
 <P>You can specify the accelerated styles explicitly in your input script
 by including their suffix, or you can use the <A HREF = "Section_start.html#start_7">-suffix command-line
 switch</A> when you invoke LAMMPS, or you can
 use the <A HREF = "suffix.html">suffix</A> command in your input script.
 </P>
 <P>See <A HREF = "Section_accelerate.html">Section_accelerate</A> of the manual for
 more instructions on how to use the accelerated styles effectively.
 </P>
 <HR>
 
 <P><B>Mixing, shift, table, tail correction, restart, rRESPA info</B>:
 </P>
 <P>For atom type pairs I,J and I != J, where types I and J correspond to
 two different element types, mixing is performed by LAMMPS as
 described above from values in the potential file.
 </P>
 <P>These pair styles does not support the <A HREF = "pair_modify.html">pair_modify</A>
 shift, table, and tail options.
 </P>
 <P>These pair styles do not write its information to <A HREF = "restart.html">binary restart
 files</A>, since it is stored in potential files.  Thus, you
 need to re-specify the pair_style, pair_coeff, and <A HREF = "fix_qeq_comb.html">fix
 qeq/comb</A> commands in an input script that reads a
 restart file.
 </P>
 <P>These pair styles can only be used via the <I>pair</I> keyword of the
 <A HREF = "run_style.html">run_style respa</A> command.  It does not support the
 <I>inner</I>, <I>middle</I>, <I>outer</I> keywords.
 </P>
 <HR>
 
 <P><B>Restrictions:</B>
 </P>
 <P>These pair styles are part of the MANYBODY package.  It is only enabled
 if LAMMPS was built with that package (which it is by default).  See
 the <A HREF = "Section_start.html#start_3">Making LAMMPS</A> section for more info.
 </P>
 <P>These pair styles requires the <A HREF = "newton.html">newton</A> setting to be "on"
 for pair interactions.
 </P>
 <P>The COMB potentials in the <I>ffield.comb</I> and <I>ffield.comb3</I> files provided 
 with LAMMPS (see the potentials directory) are parameterized for metal
 <A HREF = "units.html">units</A>.  You can use the COMB potential with any LAMMPS
 units, but you would need to create your own COMB potential file with
 coefficients listed in the appropriate units if your simulation
 doesn't use "metal" units.
 </P>
 <P><B>Related commands:</B>
 </P>
 <P><A HREF = "pair_style.html">pair_style</A>, <A HREF = "pair_coeff.html">pair_coeff</A>, 
 <A HREF = "fix_qeq_comb.html">fix_qeq/comb</A>
 </P>
 <P><B>Default:</B> none
 </P>
 <HR>
 
 <A NAME = "COMB"></A>
 
-<P><B>(COMB)</B> S. R. Phillpot and S. B. Sinnott, Science 325, 1634-1635 (2009)
+<P><B>(COMB)</B>  T.-R. Shan, B. D. Devine, T. W. Kemper, S. B. Sinnott, and 
+S. R. Phillpot, Phys. Rev. B 81, 125328 (2010)
 </P>
 <A NAME = "COMB3"></A>
 
 <P><B>(COMB3)</B> T. Liang, T.-R. Shan, Y.-T. Cheng, B. D. Devine, M. Noordhoek, 
-Y. Li, Z. Lu, S. R. Phillpot, and S. B. Sinnott, Mat. Sci. & Eng: R, 
-in press (DOI: 10.1016/j.mser.2013.07.001)
+Y. Li, Z. Lu, S. R. Phillpot, and S. B. Sinnott, Mat. Sci. & Eng: R 74, 
+255-279 (2013).
 </P>
 <A NAME = "Rick"></A>
 
 <P><B>(Rick)</B> S. W. Rick, S. J. Stuart, B. J. Berne, J Chem Phys 101, 6141
 (1994).
 </P>
 </HTML>
diff --git a/doc/pair_comb.txt b/doc/pair_comb.txt
index bc92b12db..9b87404d5 100644
--- a/doc/pair_comb.txt
+++ b/doc/pair_comb.txt
@@ -1,186 +1,187 @@
 "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
 
 pair_style comb command :h3
 pair_style comb/omp command :h3
 pair_style comb3 command :h3
 
 [Syntax:]
 
 pair_style comb
 pair_style comb3 keyword :pre
 
 keyword = {polar}
   {polar} value = {polar_on} or {polar_off} = whether or not to include atomic polarization :pre
 :ule
 
 [Examples:]
 
 pair_style comb
 pair_coeff * * ../potentials/ffield.comb Si
 pair_coeff * * ../potentials/ffield.comb Hf Si O :pre
 pair_style comb3 polar_off
 pair_coeff * * ../potentials/ffield.comb3 O Cu N C O :pre
 
 [Description:]
 
 Style {comb} computes the second-generation variable charge COMB
 (Charge-Optimized Many-Body) potential.  Style {comb3} computes the
 third-generation COMB potential.  These COMB potentials are described
 in "(COMB)"_#COMB and "(COMB3)"_#COMB3.  Briefly, the total energy
 {E<sub>T</sub>} of a system of atoms is given by
 
 :c,image(Eqs/pair_comb1.jpg)
 
 where {E<sub>i</sub><sup>self</sup>} is the self-energy of atom {i} 
 (including atomic ionization energies and electron affinities),
 {E<sub>ij</sub><sup>short</sup>} is the bond-order potential between 
 atoms {i} and {j},
 {E<sub>ij</sub><sup>Coul</sup>} is the Coulomb interactions,
 {E<sup>polar</sup>} is the polarization term for organic systems 
 (style {comb3} only), 
 {E<sup>vdW</sup>} is the van der Waals energy (style {comb3} only), 
 {E<sup>barr</sup>} is a charge barrier function, and 
 {E<sup>corr</sup>} are angular correction terms.
 
 The COMB potentials (styles {comb} and {comb3}) are variable charge
 potentials.  The equilibrium charge on each atom is calculated by the
 electronegativity equalization (QEq) method.  See "Rick"_#Rick for
 further details.  This is implemented by the "fix
 qeq/comb"_fix_qeq_comb.html command, which should normally be
 specified in the input script when running a model with the COMB
 potential.  The "fix qeq/comb"_fix_qeq_comb.html command has options
 that determine how often charge equilibration is performed, its
 convergence criterion, and which atoms are included in the
 calculation.
 
 Only a single pair_coeff command is used with the {comb} and {comb3}
 styles which specifies the COMB potential file with parameters for all
 needed elements.  These are mapped to LAMMPS atom types by specifying
 N additional arguments after the potential file in the pair_coeff
 command, where N is the number of LAMMPS atom types.
 
 For example, if your LAMMPS simulation of a Si/SiO<sub>2</sub>/
 HfO<sub>2</sub> interface has 4 atom types, and you want the 1st and
 last to be Si, the 2nd to be Hf, and the 3rd to be O, and you would
 use the following pair_coeff command:
 
 pair_coeff * * ../potentials/ffield.comb Si Hf O Si :pre
 
 The first two arguments must be * * so as to span all LAMMPS atom
 types.  The first and last Si arguments map LAMMPS atom types 1 and 4
 to the Si element in the {ffield.comb} file.  The second Hf argument
 maps LAMMPS atom type 2 to the Hf element, and the third O argument
 maps LAMMPS atom type 3 to the O element in the potential file.  If a
 mapping value is specified as NULL, the mapping is not performed.
 This can be used when a {comb} potential is used as part of the
 {hybrid} pair style.  The NULL values are placeholders for atom types
 that will be used with other potentials.
 
 For style {comb}, the provided potential file {ffield.comb} contains
 all currently-available 2nd generation COMB parameterizations: for Si,
 Cu, Hf, Ti, O, their oxides and Zr, Zn and U metals.  For style
 {comb3}, the potential file {ffield.comb3} contains all
 currently-available 3rd generation COMB paramterizations: O, Cu, N, C,
-H, Ti and Zn.  The status of the optimization of the compounds, for
+H, Ti, Zn and Zr.  The status of the optimization of the compounds, for
 example Cu<sub>2</sub>O, TiN and hydrocarbons, are given in the
 following table:
 
 :c,image(Eqs/pair_comb2.jpg)
 
 For style {comb3}, in addition to ffield.comb3, a special parameter
 file, {lib.comb3}, that is exclusively used for C/O/H systems, will be
 automatically loaded if carbon atom is detected in LAMMPS input
 structure.  Keyword {polar} indicates whether the force field includes
 the atomic polarization. Since the equilibration of the polarization
 has not yet been implemented, it can only set polar_off at present.
 
 IMPORTANT NOTE: You can not use potential file {ffield.comb} with
 style {comb3}, nor file {ffield.comb3} with style {comb}.
 
 :line
 
 Styles with a {cuda}, {gpu}, {omp}, or {opt} suffix are functionally
 the same as the corresponding style without the suffix.  They have
 been optimized to run faster, depending on your available hardware, as
 discussed in "Section_accelerate"_Section_accelerate.html of the
 manual.  The accelerated styles take the same arguments and should
 produce the same results, except for round-off and precision issues.
 
 These accelerated styles are part of the USER-CUDA, GPU, USER-OMP and OPT
 packages, respectively.  They are only enabled if LAMMPS was built with
 those packages.  See the "Making LAMMPS"_Section_start.html#start_3
 section for more info.
 
 You can specify the accelerated styles explicitly in your input script
 by including their suffix, or you can use the "-suffix command-line
 switch"_Section_start.html#start_7 when you invoke LAMMPS, or you can
 use the "suffix"_suffix.html command in your input script.
 
 See "Section_accelerate"_Section_accelerate.html of the manual for
 more instructions on how to use the accelerated styles effectively.
 
 :line
 
 [Mixing, shift, table, tail correction, restart, rRESPA info]:
 
 For atom type pairs I,J and I != J, where types I and J correspond to
 two different element types, mixing is performed by LAMMPS as
 described above from values in the potential file.
 
 These pair styles does not support the "pair_modify"_pair_modify.html
 shift, table, and tail options.
 
 These pair styles do not write its information to "binary restart
 files"_restart.html, since it is stored in potential files.  Thus, you
 need to re-specify the pair_style, pair_coeff, and "fix
 qeq/comb"_fix_qeq_comb.html commands in an input script that reads a
 restart file.
 
 These pair styles can only be used via the {pair} keyword of the
 "run_style respa"_run_style.html command.  It does not support the
 {inner}, {middle}, {outer} keywords.
 
 :line
 
 [Restrictions:]
 
 These pair styles are part of the MANYBODY package.  It is only enabled
 if LAMMPS was built with that package (which it is by default).  See
 the "Making LAMMPS"_Section_start.html#start_3 section for more info.
 
 These pair styles requires the "newton"_newton.html setting to be "on"
 for pair interactions.
 
 The COMB potentials in the {ffield.comb} and {ffield.comb3} files provided 
 with LAMMPS (see the potentials directory) are parameterized for metal
 "units"_units.html.  You can use the COMB potential with any LAMMPS
 units, but you would need to create your own COMB potential file with
 coefficients listed in the appropriate units if your simulation
 doesn't use "metal" units.
 
 [Related commands:]
 
 "pair_style"_pair_style.html, "pair_coeff"_pair_coeff.html, 
 "fix_qeq/comb"_fix_qeq_comb.html
 
 [Default:] none
 
 :line
 
 :link(COMB)
-[(COMB)] S. R. Phillpot and S. B. Sinnott, Science 325, 1634-1635 (2009)
+[(COMB)]  T.-R. Shan, B. D. Devine, T. W. Kemper, S. B. Sinnott, and 
+S. R. Phillpot, Phys. Rev. B 81, 125328 (2010)
 
 :link(COMB3)
 [(COMB3)] T. Liang, T.-R. Shan, Y.-T. Cheng, B. D. Devine, M. Noordhoek, 
-Y. Li, Z. Lu, S. R. Phillpot, and S. B. Sinnott, Mat. Sci. & Eng: R, 
-in press (DOI: 10.1016/j.mser.2013.07.001)
+Y. Li, Z. Lu, S. R. Phillpot, and S. B. Sinnott, Mat. Sci. & Eng: R 74, 
+255-279 (2013).
 
 :link(Rick)
 [(Rick)] S. W. Rick, S. J. Stuart, B. J. Berne, J Chem Phys 101, 6141
 (1994).