From 9f1c507de30633c7d05971f79f2a767528174190 Mon Sep 17 00:00:00 2001
From: sjplimp <sjplimp@f3b2605a-c512-4ea7-a41b-209d697bcdaa>
Date: Fri, 18 Jun 2010 15:59:35 +0000
Subject: [PATCH 01/20] git-svn-id:
 svn://svn.icms.temple.edu/lammps-ro/trunk@4293
 f3b2605a-c512-4ea7-a41b-209d697bcdaa

---
 doc/Section_commands.html | 20 ++++++++--------
 doc/Section_commands.txt  |  1 +
 doc/fix.html              |  3 ++-
 doc/fix.txt               |  1 +
 doc/pair_soft.html        | 50 +++++++++++++++++++++++----------------
 doc/pair_soft.txt         | 50 +++++++++++++++++++++++----------------
 doc/run.html              |  3 +--
 doc/run.txt               |  3 +--
 8 files changed, 74 insertions(+), 57 deletions(-)

diff --git a/doc/Section_commands.html b/doc/Section_commands.html
index b91ad8fce7..62fa8d8b35 100644
--- a/doc/Section_commands.html
+++ b/doc/Section_commands.html
@@ -323,16 +323,16 @@ in the command's documentation.
 of each style or click on the style itself for a full description:
 </P>
 <DIV ALIGN=center><TABLE  BORDER=1 >
-<TR ALIGN="center"><TD ><A HREF = "fix_addforce.html">addforce</A></TD><TD ><A HREF = "fix_aveforce.html">aveforce</A></TD><TD ><A HREF = "fix_ave_atom.html">ave/atom</A></TD><TD ><A HREF = "fix_ave_histo.html">ave/histo</A></TD><TD ><A HREF = "fix_ave_spatial.html">ave/spatial</A></TD><TD ><A HREF = "fix_ave_time.html">ave/time</A></TD><TD ><A HREF = "fix_bond_break.html">bond/break</A></TD><TD ><A HREF = "fix_bond_create.html">bond/create</A></TD></TR>
-<TR ALIGN="center"><TD ><A HREF = "fix_bond_swap.html">bond/swap</A></TD><TD ><A HREF = "fix_box_relax.html">box/relax</A></TD><TD ><A HREF = "fix_deform.html">deform</A></TD><TD ><A HREF = "fix_deposit.html">deposit</A></TD><TD ><A HREF = "fix_drag.html">drag</A></TD><TD ><A HREF = "fix_dt_reset.html">dt/reset</A></TD><TD ><A HREF = "fix_efield.html">efield</A></TD><TD ><A HREF = "fix_enforce2d.html">enforce2d</A></TD></TR>
-<TR ALIGN="center"><TD ><A HREF = "fix_evaporate.html">evaporate</A></TD><TD ><A HREF = "fix_freeze.html">freeze</A></TD><TD ><A HREF = "fix_gravity.html">gravity</A></TD><TD ><A HREF = "fix_heat.html">heat</A></TD><TD ><A HREF = "fix_indent.html">indent</A></TD><TD ><A HREF = "fix_langevin.html">langevin</A></TD><TD ><A HREF = "fix_lineforce.html">lineforce</A></TD><TD ><A HREF = "fix_momentum.html">momentum</A></TD></TR>
-<TR ALIGN="center"><TD ><A HREF = "fix_move.html">move</A></TD><TD ><A HREF = "fix_msst.html">msst</A></TD><TD ><A HREF = "fix_nh.html">nph</A></TD><TD ><A HREF = "fix_nph_asphere.html">nph/asphere</A></TD><TD ><A HREF = "fix_nph_sphere.html">nph/sphere</A></TD><TD ><A HREF = "fix_nh.html">npt</A></TD><TD ><A HREF = "fix_npt_asphere.html">npt/asphere</A></TD><TD ><A HREF = "fix_npt_sphere.html">npt/sphere</A></TD></TR>
-<TR ALIGN="center"><TD ><A HREF = "fix_nve.html">nve</A></TD><TD ><A HREF = "fix_nve_asphere.html">nve/asphere</A></TD><TD ><A HREF = "fix_nve_limit.html">nve/limit</A></TD><TD ><A HREF = "fix_nve_noforce.html">nve/noforce</A></TD><TD ><A HREF = "fix_nve_sphere.html">nve/sphere</A></TD><TD ><A HREF = "fix_nh.html">nvt</A></TD><TD ><A HREF = "fix_nvt_asphere.html">nvt/asphere</A></TD><TD ><A HREF = "fix_nvt_sllod.html">nvt/sllod</A></TD></TR>
-<TR ALIGN="center"><TD ><A HREF = "fix_nvt_sphere.html">nvt/sphere</A></TD><TD ><A HREF = "fix_orient_fcc.html">orient/fcc</A></TD><TD ><A HREF = "fix_planeforce.html">planeforce</A></TD><TD ><A HREF = "fix_poems.html">poems</A></TD><TD ><A HREF = "fix_pour.html">pour</A></TD><TD ><A HREF = "fix_press_berendsen.html">press/berendsen</A></TD><TD ><A HREF = "fix_print.html">print</A></TD><TD ><A HREF = "fix_reax_bonds.html">reax/bonds</A></TD></TR>
-<TR ALIGN="center"><TD ><A HREF = "fix_recenter.html">recenter</A></TD><TD ><A HREF = "fix_rigid.html">rigid</A></TD><TD ><A HREF = "fix_rigid.html">rigid/nve</A></TD><TD ><A HREF = "fix_rigid.html">rigid/nvt</A></TD><TD ><A HREF = "fix_setforce.html">setforce</A></TD><TD ><A HREF = "fix_shake.html">shake</A></TD><TD ><A HREF = "fix_spring.html">spring</A></TD><TD ><A HREF = "fix_spring_rg.html">spring/rg</A></TD></TR>
-<TR ALIGN="center"><TD ><A HREF = "fix_spring_self.html">spring/self</A></TD><TD ><A HREF = "fix_store_force.html">store/force</A></TD><TD ><A HREF = "fix_store_state.html">store/state</A></TD><TD ><A HREF = "fix_temp_berendsen.html">temp/berendsen</A></TD><TD ><A HREF = "fix_temp_rescale.html">temp/rescale</A></TD><TD ><A HREF = "fix_thermal_conductivity.html">thermal/conductivity</A></TD><TD ><A HREF = "fix_tmd.html">tmd</A></TD><TD ><A HREF = "fix_ttm.html">ttm</A></TD></TR>
-<TR ALIGN="center"><TD ><A HREF = "fix_viscosity.html">viscosity</A></TD><TD ><A HREF = "fix_viscous.html">viscous</A></TD><TD ><A HREF = "fix_wall.html">wall/colloid</A></TD><TD ><A HREF = "fix_wall_gran.html">wall/gran</A></TD><TD ><A HREF = "fix_wall.html">wall/harmonic</A></TD><TD ><A HREF = "fix_wall.html">wall/lj126</A></TD><TD ><A HREF = "fix_wall.html">wall/lj93</A></TD><TD ><A HREF = "fix_wall_reflect.html">wall/reflect</A></TD></TR>
-<TR ALIGN="center"><TD ><A HREF = "fix_wall_region.html">wall/region</A> 
+<TR ALIGN="center"><TD ><A HREF = "fix_adapt.html">adapt</A></TD><TD ><A HREF = "fix_addforce.html">addforce</A></TD><TD ><A HREF = "fix_aveforce.html">aveforce</A></TD><TD ><A HREF = "fix_ave_atom.html">ave/atom</A></TD><TD ><A HREF = "fix_ave_histo.html">ave/histo</A></TD><TD ><A HREF = "fix_ave_spatial.html">ave/spatial</A></TD><TD ><A HREF = "fix_ave_time.html">ave/time</A></TD><TD ><A HREF = "fix_bond_break.html">bond/break</A></TD></TR>
+<TR ALIGN="center"><TD ><A HREF = "fix_bond_create.html">bond/create</A></TD><TD ><A HREF = "fix_bond_swap.html">bond/swap</A></TD><TD ><A HREF = "fix_box_relax.html">box/relax</A></TD><TD ><A HREF = "fix_deform.html">deform</A></TD><TD ><A HREF = "fix_deposit.html">deposit</A></TD><TD ><A HREF = "fix_drag.html">drag</A></TD><TD ><A HREF = "fix_dt_reset.html">dt/reset</A></TD><TD ><A HREF = "fix_efield.html">efield</A></TD></TR>
+<TR ALIGN="center"><TD ><A HREF = "fix_enforce2d.html">enforce2d</A></TD><TD ><A HREF = "fix_evaporate.html">evaporate</A></TD><TD ><A HREF = "fix_freeze.html">freeze</A></TD><TD ><A HREF = "fix_gravity.html">gravity</A></TD><TD ><A HREF = "fix_heat.html">heat</A></TD><TD ><A HREF = "fix_indent.html">indent</A></TD><TD ><A HREF = "fix_langevin.html">langevin</A></TD><TD ><A HREF = "fix_lineforce.html">lineforce</A></TD></TR>
+<TR ALIGN="center"><TD ><A HREF = "fix_momentum.html">momentum</A></TD><TD ><A HREF = "fix_move.html">move</A></TD><TD ><A HREF = "fix_msst.html">msst</A></TD><TD ><A HREF = "fix_nh.html">nph</A></TD><TD ><A HREF = "fix_nph_asphere.html">nph/asphere</A></TD><TD ><A HREF = "fix_nph_sphere.html">nph/sphere</A></TD><TD ><A HREF = "fix_nh.html">npt</A></TD><TD ><A HREF = "fix_npt_asphere.html">npt/asphere</A></TD></TR>
+<TR ALIGN="center"><TD ><A HREF = "fix_npt_sphere.html">npt/sphere</A></TD><TD ><A HREF = "fix_nve.html">nve</A></TD><TD ><A HREF = "fix_nve_asphere.html">nve/asphere</A></TD><TD ><A HREF = "fix_nve_limit.html">nve/limit</A></TD><TD ><A HREF = "fix_nve_noforce.html">nve/noforce</A></TD><TD ><A HREF = "fix_nve_sphere.html">nve/sphere</A></TD><TD ><A HREF = "fix_nh.html">nvt</A></TD><TD ><A HREF = "fix_nvt_asphere.html">nvt/asphere</A></TD></TR>
+<TR ALIGN="center"><TD ><A HREF = "fix_nvt_sllod.html">nvt/sllod</A></TD><TD ><A HREF = "fix_nvt_sphere.html">nvt/sphere</A></TD><TD ><A HREF = "fix_orient_fcc.html">orient/fcc</A></TD><TD ><A HREF = "fix_planeforce.html">planeforce</A></TD><TD ><A HREF = "fix_poems.html">poems</A></TD><TD ><A HREF = "fix_pour.html">pour</A></TD><TD ><A HREF = "fix_press_berendsen.html">press/berendsen</A></TD><TD ><A HREF = "fix_print.html">print</A></TD></TR>
+<TR ALIGN="center"><TD ><A HREF = "fix_reax_bonds.html">reax/bonds</A></TD><TD ><A HREF = "fix_recenter.html">recenter</A></TD><TD ><A HREF = "fix_rigid.html">rigid</A></TD><TD ><A HREF = "fix_rigid.html">rigid/nve</A></TD><TD ><A HREF = "fix_rigid.html">rigid/nvt</A></TD><TD ><A HREF = "fix_setforce.html">setforce</A></TD><TD ><A HREF = "fix_shake.html">shake</A></TD><TD ><A HREF = "fix_spring.html">spring</A></TD></TR>
+<TR ALIGN="center"><TD ><A HREF = "fix_spring_rg.html">spring/rg</A></TD><TD ><A HREF = "fix_spring_self.html">spring/self</A></TD><TD ><A HREF = "fix_store_force.html">store/force</A></TD><TD ><A HREF = "fix_store_state.html">store/state</A></TD><TD ><A HREF = "fix_temp_berendsen.html">temp/berendsen</A></TD><TD ><A HREF = "fix_temp_rescale.html">temp/rescale</A></TD><TD ><A HREF = "fix_thermal_conductivity.html">thermal/conductivity</A></TD><TD ><A HREF = "fix_tmd.html">tmd</A></TD></TR>
+<TR ALIGN="center"><TD ><A HREF = "fix_ttm.html">ttm</A></TD><TD ><A HREF = "fix_viscosity.html">viscosity</A></TD><TD ><A HREF = "fix_viscous.html">viscous</A></TD><TD ><A HREF = "fix_wall.html">wall/colloid</A></TD><TD ><A HREF = "fix_wall_gran.html">wall/gran</A></TD><TD ><A HREF = "fix_wall.html">wall/harmonic</A></TD><TD ><A HREF = "fix_wall.html">wall/lj126</A></TD><TD ><A HREF = "fix_wall.html">wall/lj93</A></TD></TR>
+<TR ALIGN="center"><TD ><A HREF = "fix_wall_reflect.html">wall/reflect</A></TD><TD ><A HREF = "fix_wall_region.html">wall/region</A> 
 </TD></TR></TABLE></DIV>
 
 <P>These are fix styles contributed by users, which can be used if
diff --git a/doc/Section_commands.txt b/doc/Section_commands.txt
index f31218c96e..c701c89cdb 100644
--- a/doc/Section_commands.txt
+++ b/doc/Section_commands.txt
@@ -382,6 +382,7 @@ Fix styles :h4
 See the "fix"_fix.html command for one-line descriptions
 of each style or click on the style itself for a full description:
 
+"adapt"_fix_adapt.html,
 "addforce"_fix_addforce.html,
 "aveforce"_fix_aveforce.html,
 "ave/atom"_fix_ave_atom.html,
diff --git a/doc/fix.html b/doc/fix.html
index 71c151b381..fc6e8828cc 100644
--- a/doc/fix.html
+++ b/doc/fix.html
@@ -160,7 +160,8 @@ for further info.
 arguments and what it does, as listed below.  Here is an alphabetic
 list of fix styles available in LAMMPS:
 </P>
-<UL><LI><A HREF = "fix_addforce.html">addforce</A> - add a force to each atom
+<UL><LI><A HREF = "fix_adapt.html">adapt</A> - change a simulation parameter over time
+<LI><A HREF = "fix_addforce.html">addforce</A> - add a force to each atom
 <LI><A HREF = "fix_aveforce.html">aveforce</A> - add an averaged force to each atom
 <LI><A HREF = "fix_ave_atom.html">ave/atom</A> - compute per-atom time-averaged quantities
 <LI><A HREF = "fix_ave_atom.html">ave/histo</A> - compute/output time-averaged histograms
diff --git a/doc/fix.txt b/doc/fix.txt
index a9bbf59241..0ec33264f2 100644
--- a/doc/fix.txt
+++ b/doc/fix.txt
@@ -155,6 +155,7 @@ Each fix style has its own documentation page which describes its
 arguments and what it does, as listed below.  Here is an alphabetic
 list of fix styles available in LAMMPS:
 
+"adapt"_fix_adapt.html - change a simulation parameter over time
 "addforce"_fix_addforce.html - add a force to each atom
 "aveforce"_fix_aveforce.html - add an averaged force to each atom
 "ave/atom"_fix_ave_atom.html - compute per-atom time-averaged quantities
diff --git a/doc/pair_soft.html b/doc/pair_soft.html
index bcacfea72f..d2aab51c72 100644
--- a/doc/pair_soft.html
+++ b/doc/pair_soft.html
@@ -20,8 +20,8 @@
 <P><B>Examples:</B>
 </P>
 <PRE>pair_style soft 2.5
-pair_coeff * * 0.0 60.0
-pair_coeff 1 1 0.0 60.0 3.0 
+pair_coeff * * 10.0
+pair_coeff 1 1 10.0 3.0 
 </PRE>
 <P><B>Description:</B>
 </P>
@@ -30,11 +30,12 @@ pair_coeff 1 1 0.0 60.0 3.0
 <CENTER><IMG SRC = "Eqs/pair_soft.jpg">
 </CENTER>
 <P>It is useful for pushing apart overlapping atoms, since it does not
-blow up as r goes to 0.  A is a pre-factor that varies in time from
-the start to the end of the run.  The <A HREF = "run.html">run</A> command documents
-how to make the ramping take place across multiple runs.  Rc is the
-cutoff.  See the <A HREF = "fix_nve_limit.html">fix nve/limit</A> command for
-another way to push apart overlapping atoms.
+blow up as r goes to 0.  A is a pre-factor that can be made to vary in
+time from the start to the end of the run (see discussion below),
+e.g. to start with a very soft potential and slowly harden the
+interactions over time.  Rc is the cutoff.  See the <A HREF = "fix_nve_limit.html">fix
+nve/limit</A> command for another way to push apart
+overlapping atoms.
 </P>
 <P>The following coefficients must be defined for each pair of atom types
 via the <A HREF = "pair_coeff.html">pair_coeff</A> command as in the examples above,
@@ -42,28 +43,34 @@ or in the data file or restart files read by the
 <A HREF = "read_data.html">read_data</A> or <A HREF = "read_restart.html">read_restart</A>
 commands, or by mixing as described below:
 </P>
-<UL><LI>Astart (energy units)
-<LI>Astop (energy units)
+<UL><LI>A (energy units)
 <LI>cutoff (distance units) 
 </UL>
-<P>Astart and Astop are the values of the prefactor at the start and end
-of the next run.  At intermediate times the value of A will be ramped
-between these 2 values.  Note that before performing a 2nd run, you
-will want to adjust the values of Astart and Astop for all type pairs,
-or switch to a new pair style.
-</P>
 <P>The last coefficient is optional.  If not specified, the global soft
 cutoff is used.
 </P>
+<P>The <A HREF = "fix_adapt.html">fix adapt</A> command can be used to vary A over the
+course of a simulation.  For example these commands will vary the
+prefactor A for all pairwise interactions from 0.0 at the beginning to
+30.0 at the end of a 10,000 step run:
+</P>
+<PRE>variable prefactor equal 30.0*elapsed/10000
+fix 1 all adapt 1 pair soft a * * prefactor 
+</PRE>
+<P>Note that a formula defined by an <A HREF = "variable.html">equal-style variable</A>
+can use the current timestep, elapsed time in the current run, elapsed
+time since the beginning of a series of runs, as well as access other
+variables.
+</P>
 <HR>
 
 <P><B>Mixing, shift, table, tail correction, restart, rRESPA info</B>:
 </P>
-<P>For atom type pairs I,J and I != J, the Astart, Astop coefficients and
-cutoff distance for this pair style can be mixed.  Astart and Atop are
-always mixed via a <I>geometric</I> rule.  The cutoff is mixed according to
-the pair_modify mix value.  The default mix value is <I>geometric</I>.  See
-the "pair_modify" command for details.
+<P>For atom type pairs I,J and I != J, the A coefficient and cutoff
+distance for this pair style can be mixed.  A is always mixed via a
+<I>geometric</I> rule.  The cutoff is mixed according to the pair_modify
+mix value.  The default mix value is <I>geometric</I>.  See the
+"pair_modify" command for details.
 </P>
 <P>This pair style does not support the <A HREF = "pair_modify.html">pair_modify</A>
 shift option, since the pair interaction goes to 0.0 at the cutoff.
@@ -85,7 +92,8 @@ to be specified in an input script that reads a restart file.
 </P>
 <P><B>Related commands:</B>
 </P>
-<P><A HREF = "pair_coeff.html">pair_coeff</A>, <A HREF = "fix_nve_limit.html">fix nve/limit</A>
+<P><A HREF = "pair_coeff.html">pair_coeff</A>, <A HREF = "fix_nve_limit.html">fix nve/limit</A>, <A HREF = "fix_adapt.html">fix
+adapt</A>
 </P>
 <P><B>Default:</B> none
 </P>
diff --git a/doc/pair_soft.txt b/doc/pair_soft.txt
index a5be96d086..86425b648f 100644
--- a/doc/pair_soft.txt
+++ b/doc/pair_soft.txt
@@ -17,8 +17,8 @@ cutoff = global cutoff for soft interactions (distance units) :ul
 [Examples:]
 
 pair_style soft 2.5
-pair_coeff * * 0.0 60.0
-pair_coeff 1 1 0.0 60.0 3.0 :pre
+pair_coeff * * 10.0
+pair_coeff 1 1 10.0 3.0 :pre
 
 [Description:]
 
@@ -27,11 +27,12 @@ Style {soft} computes pairwise interactions with the formula
 :c,image(Eqs/pair_soft.jpg)
 
 It is useful for pushing apart overlapping atoms, since it does not
-blow up as r goes to 0.  A is a pre-factor that varies in time from
-the start to the end of the run.  The "run"_run.html command documents
-how to make the ramping take place across multiple runs.  Rc is the
-cutoff.  See the "fix nve/limit"_fix_nve_limit.html command for
-another way to push apart overlapping atoms.
+blow up as r goes to 0.  A is a pre-factor that can be made to vary in
+time from the start to the end of the run (see discussion below),
+e.g. to start with a very soft potential and slowly harden the
+interactions over time.  Rc is the cutoff.  See the "fix
+nve/limit"_fix_nve_limit.html command for another way to push apart
+overlapping atoms.
 
 The following coefficients must be defined for each pair of atom types
 via the "pair_coeff"_pair_coeff.html command as in the examples above,
@@ -39,28 +40,34 @@ or in the data file or restart files read by the
 "read_data"_read_data.html or "read_restart"_read_restart.html
 commands, or by mixing as described below:
 
-Astart (energy units)
-Astop (energy units)
+A (energy units)
 cutoff (distance units) :ul
 
-Astart and Astop are the values of the prefactor at the start and end
-of the next run.  At intermediate times the value of A will be ramped
-between these 2 values.  Note that before performing a 2nd run, you
-will want to adjust the values of Astart and Astop for all type pairs,
-or switch to a new pair style.
-
 The last coefficient is optional.  If not specified, the global soft
 cutoff is used.
 
+The "fix adapt"_fix_adapt.html command can be used to vary A over the
+course of a simulation.  For example these commands will vary the
+prefactor A for all pairwise interactions from 0.0 at the beginning to
+30.0 at the end of a 10,000 step run:
+
+variable prefactor equal 30.0*elapsed/10000
+fix 1 all adapt 1 pair soft a * * prefactor :pre
+
+Note that a formula defined by an "equal-style variable"_variable.html
+can use the current timestep, elapsed time in the current run, elapsed
+time since the beginning of a series of runs, as well as access other
+variables.
+
 :line
 
 [Mixing, shift, table, tail correction, restart, rRESPA info]:
 
-For atom type pairs I,J and I != J, the Astart, Astop coefficients and
-cutoff distance for this pair style can be mixed.  Astart and Atop are
-always mixed via a {geometric} rule.  The cutoff is mixed according to
-the pair_modify mix value.  The default mix value is {geometric}.  See
-the "pair_modify" command for details.
+For atom type pairs I,J and I != J, the A coefficient and cutoff
+distance for this pair style can be mixed.  A is always mixed via a
+{geometric} rule.  The cutoff is mixed according to the pair_modify
+mix value.  The default mix value is {geometric}.  See the
+"pair_modify" command for details.
 
 This pair style does not support the "pair_modify"_pair_modify.html
 shift option, since the pair interaction goes to 0.0 at the cutoff.
@@ -82,6 +89,7 @@ This pair style can only be used via the {pair} keyword of the
 
 [Related commands:]
 
-"pair_coeff"_pair_coeff.html, "fix nve/limit"_fix_nve_limit.html
+"pair_coeff"_pair_coeff.html, "fix nve/limit"_fix_nve_limit.html, "fix
+adapt"_fix_adapt.html
 
 [Default:] none
diff --git a/doc/run.html b/doc/run.html
index 05524ee43b..a93d102276 100644
--- a/doc/run.html
+++ b/doc/run.html
@@ -69,8 +69,7 @@ performed and you want a <A HREF = "fix.html">fix</A> command that changes some
 value over time (e.g. temperature) to make the change across the
 entire set of runs and not just a single run.  See the doc page for
 individual fixes to see which ones can be used with the <I>start/stop</I>
-keywords.  The <A HREF = "pair_style.html">pair_style soft</A> potential also
-changes its pair potential coefficients in this manner.
+keywords.
 </P>
 <P>For example, consider this fix followed by 10 run commands:
 </P>
diff --git a/doc/run.txt b/doc/run.txt
index c80f308dc6..56f47be800 100644
--- a/doc/run.txt
+++ b/doc/run.txt
@@ -62,8 +62,7 @@ performed and you want a "fix"_fix.html command that changes some
 value over time (e.g. temperature) to make the change across the
 entire set of runs and not just a single run.  See the doc page for
 individual fixes to see which ones can be used with the {start/stop}
-keywords.  The "pair_style soft"_pair_style.html potential also
-changes its pair potential coefficients in this manner.
+keywords.
 
 For example, consider this fix followed by 10 run commands:
 

From 672b6728e8df9be90a0dcfe26f995c44a6ac0b5f Mon Sep 17 00:00:00 2001
From: sjplimp <sjplimp@f3b2605a-c512-4ea7-a41b-209d697bcdaa>
Date: Fri, 18 Jun 2010 15:59:41 +0000
Subject: [PATCH 02/20] git-svn-id:
 svn://svn.icms.temple.edu/lammps-ro/trunk@4294
 f3b2605a-c512-4ea7-a41b-209d697bcdaa

---
 doc/pair_lubricate.html | 4 ++++
 doc/pair_lubricate.txt  | 4 ++++
 2 files changed, 8 insertions(+)

diff --git a/doc/pair_lubricate.html b/doc/pair_lubricate.html
index 6f47436ec4..fbd8aa7fee 100644
--- a/doc/pair_lubricate.html
+++ b/doc/pair_lubricate.html
@@ -43,6 +43,10 @@ their relative velocities in the presence of a background solvent with
 viscosity mu.  Note that this is dynamic viscosity which has units of
 mass/distance/time, not kinematic viscosity.
 </P>
+<P>The viscosity mu can be varied in a time-dependent manner over the
+course of a simluation.  See the <A HREF = "fix_adapt.html">fix adapt</A> command
+for details.
+</P>
 <P>Rc is the outer cutoff specified in the pair_style command, the
 translational velocities of the 2 particles are v1 and v2, the angular
 velocities are w1 and w2, and n is the unit vector in the direction
diff --git a/doc/pair_lubricate.txt b/doc/pair_lubricate.txt
index aa1caae94b..b793a7f546 100644
--- a/doc/pair_lubricate.txt
+++ b/doc/pair_lubricate.txt
@@ -40,6 +40,10 @@ their relative velocities in the presence of a background solvent with
 viscosity mu.  Note that this is dynamic viscosity which has units of
 mass/distance/time, not kinematic viscosity.
 
+The viscosity mu can be varied in a time-dependent manner over the
+course of a simluation.  See the "fix adapt"_fix_adapt.html command
+for details.
+
 Rc is the outer cutoff specified in the pair_style command, the
 translational velocities of the 2 particles are v1 and v2, the angular
 velocities are w1 and w2, and n is the unit vector in the direction

From f3404df0f437d9f73fa6f35b657275fbbd736f46 Mon Sep 17 00:00:00 2001
From: sjplimp <sjplimp@f3b2605a-c512-4ea7-a41b-209d697bcdaa>
Date: Fri, 18 Jun 2010 16:00:04 +0000
Subject: [PATCH 03/20] git-svn-id:
 svn://svn.icms.temple.edu/lammps-ro/trunk@4295
 f3b2605a-c512-4ea7-a41b-209d697bcdaa

---
 src/pair.h | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/src/pair.h b/src/pair.h
index 8bb36b13c8..3b3e3f9f36 100644
--- a/src/pair.h
+++ b/src/pair.h
@@ -97,6 +97,8 @@ class Pair : protected Pointers {
   virtual void swap_eam(double *, double **) {}
   virtual void reset_dt() {}
   virtual void min_pointers(double **, double **) {}
+  virtual int pre_adapt(char *, int, int, int, int) {return -1;}
+  virtual void adapt(int, int, int, int, int, double) {}
 
  protected:
   int allocated;                       // 0/1 = whether arrays are allocated

From c3414bf5f612dbe3e280b529dfc363cffc52cb04 Mon Sep 17 00:00:00 2001
From: sjplimp <sjplimp@f3b2605a-c512-4ea7-a41b-209d697bcdaa>
Date: Fri, 18 Jun 2010 16:00:10 +0000
Subject: [PATCH 04/20] git-svn-id:
 svn://svn.icms.temple.edu/lammps-ro/trunk@4296
 f3b2605a-c512-4ea7-a41b-209d697bcdaa

---
 src/fix_adapt.cpp | 196 ++++++++++++++++++++++++++++++++++++++++++++++
 src/fix_adapt.h   |  50 ++++++++++++
 2 files changed, 246 insertions(+)
 create mode 100644 src/fix_adapt.cpp
 create mode 100644 src/fix_adapt.h

diff --git a/src/fix_adapt.cpp b/src/fix_adapt.cpp
new file mode 100644
index 0000000000..65397b1894
--- /dev/null
+++ b/src/fix_adapt.cpp
@@ -0,0 +1,196 @@
+/* ----------------------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   http://lammps.sandia.gov, Sandia National Laboratories
+   Steve Plimpton, sjplimp@sandia.gov
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under 
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#include "string.h"
+#include "stdlib.h"
+#include "fix_adapt.h"
+#include "atom.h"
+#include "force.h"
+#include "pair.h"
+#include "input.h"
+#include "variable.h"
+#include "error.h"
+
+using namespace LAMMPS_NS;
+
+enum{PAIR,ATOM};
+enum{DIAMETER};
+
+/* ---------------------------------------------------------------------- */
+
+FixAdapt::FixAdapt(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg)
+{
+  if (narg < 4) error->all("Illegal fix adapt command");
+  nevery = atoi(arg[3]);
+  if (nevery <= 0) error->all("Illegal fix adapt command");
+
+  // count # of adaptations
+
+  nadapt = 0;
+
+  int iarg = 4;
+  while (iarg < narg) {
+    if (strcmp(arg[iarg],"pair") == 0) {
+      if (iarg+6 > narg) error->all("Illegal fix adapt command");
+      nadapt++;
+      iarg += 6;
+    } else if (strcmp(arg[iarg],"atom") == 0) {
+      if (iarg+3 > narg) error->all("Illegal fix adapt command");
+      nadapt++;
+      iarg += 3;
+    } else error->all("Illegal fix adapt command");
+  }
+
+  // allocate per-adapt vectors
+
+  which = new int[nadapt];
+  pair = new char*[nadapt];
+  param = new char*[nadapt];
+  ilo = new int[nadapt];
+  ihi = new int[nadapt];
+  jlo = new int[nadapt];
+  jhi = new int[nadapt];
+  var = new char*[nadapt];
+  ivar = new int[nadapt];
+  pairptr = new Pair*[nadapt];
+  pairindex = new int[nadapt];
+  awhich = new int[nadapt];
+
+  // parse keywords
+
+  nadapt = 0;
+
+  iarg = 4;
+  while (iarg < narg) {
+    if (strcmp(arg[iarg],"pair") == 0) {
+      if (iarg+6 > narg) error->all("Illegal fix adapt command");
+      which[nadapt] = PAIR;
+      int n = strlen(arg[iarg+1]) + 1;
+      pair[nadapt] = new char[n];
+      strcpy(pair[nadapt],arg[iarg+1]);
+      n = strlen(arg[iarg+2]) + 1;
+      param[nadapt] = new char[n];
+      strcpy(param[nadapt],arg[iarg+2]);
+      force->bounds(arg[iarg+3],atom->ntypes,ilo[nadapt],ihi[nadapt]);
+      force->bounds(arg[iarg+4],atom->ntypes,jlo[nadapt],jhi[nadapt]);
+      n = strlen(arg[iarg+5]) + 1;
+      var[nadapt] = new char[n];
+      strcpy(var[nadapt],arg[iarg+5]);
+      nadapt++;
+      iarg += 6;
+    } else if (strcmp(arg[iarg],"atom") == 0) {
+      if (iarg+3 > narg) error->all("Illegal fix adapt command");
+      which[nadapt] = ATOM;
+      int n = strlen(arg[iarg+1]) + 1;
+      param[nadapt] = new char[n];
+      strcpy(param[nadapt],arg[iarg+1]);
+      n = strlen(arg[iarg+2]) + 1;
+      var[nadapt] = new char[n];
+      strcpy(var[nadapt],arg[iarg+2]);
+      nadapt++;
+      iarg += 3;
+    } else error->all("Illegal fix adapt command");
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
+FixAdapt::~FixAdapt()
+{
+  for (int i = 0; i < nadapt; i++) {
+    if (which[i] == PAIR) delete [] pair[i];
+    delete [] param[i];
+    delete [] var[i];
+  }
+  delete [] which;
+  delete [] pair;
+  delete [] param;
+  delete [] ilo;
+  delete [] ihi;
+  delete [] jlo;
+  delete [] jhi;
+  delete [] var;
+  delete [] ivar;
+  delete [] pairptr;
+  delete [] pairindex;
+  delete [] awhich;
+}
+
+/* ---------------------------------------------------------------------- */
+
+int FixAdapt::setmask()
+{
+  int mask = 0;
+  mask |= PRE_FORCE;
+  return mask;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixAdapt::init()
+{
+  // error checks
+
+  for (int m = 0; m < nadapt; m++) {
+    if (which[m] == PAIR) {
+      pairptr[m] = force->pair_match(pair[m],1);
+      if (pairptr[m] == NULL) 
+	error->all("Fix adapt pair style does not exist");
+      pairindex[m] = 
+	pairptr[m]->pre_adapt(param[m],ilo[m],ihi[m],jlo[m],jhi[m]);
+      if (pairindex[m] < 0)
+	error->all("Fix adapt pair style is not compatible");
+    } else if (which[m] == ATOM) {
+      if (strcmp(param[m],"diameter") == 0) {
+	awhich[m] = DIAMETER;
+	if (!atom->radius_flag)
+	  error->all("Fix adapt requires atom attribute radius");
+      } else error->all("Invalid fix adapt atom attribute");
+    }
+
+    ivar[m] = input->variable->find(var[m]);
+    if (ivar[m] < 0) error->all("Variable name for fix adapt does not exist");
+    if (!input->variable->equalstyle(ivar[m]))
+      error->all("Variable for fix adapt is not equal style");
+  }
+
+  // set params to values for initial force calculation
+  // needs to happen here in init() instead of setup()
+  // because modify->setup() is called after pre-Verlet forces are computed
+
+  pre_force(0);
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixAdapt::pre_force(int vflag)
+{
+  for (int m = 0; m < nadapt; m++) {
+    double value = input->variable->compute_equal(ivar[m]);
+
+    if (which[m] == PAIR)
+      pairptr[m]->adapt(pairindex[m],ilo[m],ihi[m],jlo[m],jhi[m],value);
+
+    else if (which[m] == ATOM) {
+      double *radius = atom->radius;
+      int *mask = atom->mask;
+      int nlocal = atom->nlocal;
+
+      if (awhich[m] == DIAMETER) {
+	for (int i = 0; i < nlocal; i++)
+	  if (mask[i] & groupbit)
+	    radius[i] = 0.5*value;
+      }
+    }
+  }
+}
diff --git a/src/fix_adapt.h b/src/fix_adapt.h
new file mode 100644
index 0000000000..1411708e71
--- /dev/null
+++ b/src/fix_adapt.h
@@ -0,0 +1,50 @@
+/* ----------------------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   http://lammps.sandia.gov, Sandia National Laboratories
+   Steve Plimpton, sjplimp@sandia.gov
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under 
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#ifdef FIX_CLASS
+
+FixStyle(adapt,FixAdapt)
+
+#else
+
+#ifndef LMP_FIX_ADAPT_H
+#define LMP_FIX_ADAPT_H
+
+#include "fix.h"
+
+namespace LAMMPS_NS {
+
+class FixAdapt : public Fix {
+ public:
+  FixAdapt(class LAMMPS *, int, char **);
+  ~FixAdapt();
+  int setmask();
+  void init();
+  void pre_force(int);
+
+ private:
+  int nadapt;
+  int *which;
+  char **pair,**param,**var;
+  int *ilo,*ihi,*jlo,*jhi;
+
+  int *ivar;
+  class Pair **pairptr;
+  int *pairindex;
+  int *awhich;
+};
+
+}
+
+#endif
+#endif

From 7e402f387132befa26016d73a437d8585aaca2d7 Mon Sep 17 00:00:00 2001
From: sjplimp <sjplimp@f3b2605a-c512-4ea7-a41b-209d697bcdaa>
Date: Fri, 18 Jun 2010 16:33:00 +0000
Subject: [PATCH 05/20] git-svn-id:
 svn://svn.icms.temple.edu/lammps-ro/trunk@4297
 f3b2605a-c512-4ea7-a41b-209d697bcdaa

---
 src/fix_adapt.cpp         | 30 ++++++++++----
 src/pair.cpp              |  8 ----
 src/pair_hybrid.h         |  2 +-
 src/pair_hybrid_overlay.h |  1 +
 src/pair_soft.cpp         | 87 +++++++++++++++++++--------------------
 src/pair_soft.h           |  3 +-
 6 files changed, 69 insertions(+), 62 deletions(-)

diff --git a/src/fix_adapt.cpp b/src/fix_adapt.cpp
index 65397b1894..f9a5ba5899 100644
--- a/src/fix_adapt.cpp
+++ b/src/fix_adapt.cpp
@@ -11,6 +11,7 @@
    See the README file in the top-level LAMMPS directory.
 ------------------------------------------------------------------------- */
 
+#include "math.h"
 #include "string.h"
 #include "stdlib.h"
 #include "fix_adapt.h"
@@ -148,14 +149,17 @@ void FixAdapt::init()
 	error->all("Fix adapt pair style does not exist");
       pairindex[m] = 
 	pairptr[m]->pre_adapt(param[m],ilo[m],ihi[m],jlo[m],jhi[m]);
-      if (pairindex[m] < 0)
-	error->all("Fix adapt pair style is not compatible");
+      if (pairindex[m] == -1)
+	error->all("Fix adapt pair parameter is not recognized");
+      if (pairindex[m] == -2)
+	error->all("Fix adapt pair types are not valid");
+
     } else if (which[m] == ATOM) {
       if (strcmp(param[m],"diameter") == 0) {
 	awhich[m] = DIAMETER;
 	if (!atom->radius_flag)
 	  error->all("Fix adapt requires atom attribute radius");
-      } else error->all("Invalid fix adapt atom attribute");
+      } else error->all("Fix adapt atom attribute is not recognized");
     }
 
     ivar[m] = input->variable->find(var[m]);
@@ -182,14 +186,26 @@ void FixAdapt::pre_force(int vflag)
       pairptr[m]->adapt(pairindex[m],ilo[m],ihi[m],jlo[m],jhi[m],value);
 
     else if (which[m] == ATOM) {
-      double *radius = atom->radius;
-      int *mask = atom->mask;
-      int nlocal = atom->nlocal;
+
+      // set radius from diameter
+      // set rmass if both rmass and density are defined
 
       if (awhich[m] == DIAMETER) {
+	int mflag = 0;
+	if (atom->rmass_flag && atom->density_flag) mflag = 1;
+	double PI = 4.0*atan(1.0);
+
+	double *radius = atom->radius;
+	double *rmass = atom->rmass;
+	double *density = atom->density;
+	int *mask = atom->mask;
+	int nlocal = atom->nlocal;
+
 	for (int i = 0; i < nlocal; i++)
-	  if (mask[i] & groupbit)
+	  if (mask[i] & groupbit) {
 	    radius[i] = 0.5*value;
+	    rmass[i] = 4.0*PI/3.0 * radius[i]*radius[i]*radius[i] * density[i];
+	  }
       }
     }
   }
diff --git a/src/pair.cpp b/src/pair.cpp
index d72566b4e2..e108c712a1 100644
--- a/src/pair.cpp
+++ b/src/pair.cpp
@@ -23,7 +23,6 @@
 #include "stdlib.h"
 #include "string.h"
 #include "pair.h"
-#include "pair_soft.h"
 #include "atom.h"
 #include "neighbor.h"
 #include "neigh_list.h"
@@ -895,13 +894,6 @@ void Pair::write_file(int narg, char **arg)
 
   force->init();
 
-  // if pair style = soft, set prefactor to final value
-
-  Pair *spair = force->pair_match("soft",1);
-  if (spair)
-    ((PairSoft *) spair)->prefactor[itype][jtype] =
-      ((PairSoft *) spair)->prestop[itype][jtype];
-
   // if pair style = any of EAM, swap in dummy fp vector
 
   double eamfp[2];
diff --git a/src/pair_hybrid.h b/src/pair_hybrid.h
index 6b66ccd5ec..9554d0595d 100644
--- a/src/pair_hybrid.h
+++ b/src/pair_hybrid.h
@@ -32,7 +32,7 @@ class PairHybrid : public Pair {
   char **keywords;              // style name of each Pair style
 
   PairHybrid(class LAMMPS *);
-  ~PairHybrid();
+  virtual ~PairHybrid();
   void compute(int, int);
   void settings(int, char **);
   virtual void coeff(int, char **);
diff --git a/src/pair_hybrid_overlay.h b/src/pair_hybrid_overlay.h
index 770e0e066c..dea84bae4c 100644
--- a/src/pair_hybrid_overlay.h
+++ b/src/pair_hybrid_overlay.h
@@ -27,6 +27,7 @@ namespace LAMMPS_NS {
 class PairHybridOverlay : public PairHybrid {
  public:
   PairHybridOverlay(class LAMMPS *);
+  ~PairHybridOverlay() {}
   void coeff(int, char **);
 
  private:
diff --git a/src/pair_soft.cpp b/src/pair_soft.cpp
index 92fab52cf6..8add30e678 100644
--- a/src/pair_soft.cpp
+++ b/src/pair_soft.cpp
@@ -14,6 +14,7 @@
 #include "math.h"
 #include "stdio.h"
 #include "stdlib.h"
+#include "string.h"
 #include "pair_soft.h"
 #include "atom.h"
 #include "comm.h"
@@ -43,8 +44,6 @@ PairSoft::~PairSoft()
     memory->destroy_2d_int_array(setflag);
     memory->destroy_2d_double_array(cutsq);
 
-    memory->destroy_2d_double_array(prestart);
-    memory->destroy_2d_double_array(prestop);
     memory->destroy_2d_double_array(prefactor);
     memory->destroy_2d_double_array(cut);
   }
@@ -63,21 +62,6 @@ void PairSoft::compute(int eflag, int vflag)
   if (eflag || vflag) ev_setup(eflag,vflag);
   else evflag = vflag_fdotr = 0;
 
-  // set current prefactor
-  // for minimization, set to prestop
-  // for dynamics, ramp from prestart to prestop
-  // for 0-step dynamics, set to prestart
-
-  double delta = update->ntimestep - update->beginstep;
-  if (update->whichflag == 2) delta = 1.0;
-  else if (update->nsteps) delta /= update->endstep - update->beginstep;
-  else delta = 0.0;
-  int ntypes = atom->ntypes;
-  for (i = 1; i <= ntypes; i++)
-    for (j = 1; j <= ntypes; j++)
-      prefactor[i][j] = prestart[i][j] + 
-	delta * (prestop[i][j] - prestart[i][j]);
-
   double **x = atom->x;
   double **f = atom->f;
   int *type = atom->type;
@@ -161,20 +145,8 @@ void PairSoft::allocate()
 
   cutsq = memory->create_2d_double_array(n+1,n+1,"pair:cutsq");
 
-  prestart = memory->create_2d_double_array(n+1,n+1,"pair:prestart");
-  prestop = memory->create_2d_double_array(n+1,n+1,"pair:prestop");
   prefactor = memory->create_2d_double_array(n+1,n+1,"pair:prefactor");
   cut = memory->create_2d_double_array(n+1,n+1,"pair:cut");
-
-  // init prestart and prestop to 0.0
-  // since pair_hybrid can use all types even if pair_soft sub-class
-  //   never sets them
-
-  for (int i = 1; i <= n; i++)
-    for (int j = i; j <= n; j++) {
-      prestart[i][j] = 0.0;
-      prestop[i][j] = 0.0;
-    }
 }
 
 /* ----------------------------------------------------------------------
@@ -203,24 +175,22 @@ void PairSoft::settings(int narg, char **arg)
 
 void PairSoft::coeff(int narg, char **arg)
 {
-  if (narg < 4 || narg > 5) error->all("Incorrect args for pair coefficients");
+  if (narg < 3 || narg > 4) error->all("Incorrect args for pair coefficients");
   if (!allocated) allocate();
 
   int ilo,ihi,jlo,jhi;
   force->bounds(arg[0],atom->ntypes,ilo,ihi);
   force->bounds(arg[1],atom->ntypes,jlo,jhi);
 
-  double prestart_one = force->numeric(arg[2]);
-  double prestop_one = force->numeric(arg[3]);
+  double prefactor_one = force->numeric(arg[2]);
 
   double cut_one = cut_global;
-  if (narg == 5) cut_one = force->numeric(arg[4]);
+  if (narg == 4) cut_one = force->numeric(arg[3]);
 
   int count = 0;
   for (int i = ilo; i <= ihi; i++) {
     for (int j = MAX(jlo,i); j <= jhi; j++) {
-      prestart[i][j] = prestart_one;
-      prestop[i][j] = prestop_one;
+      prefactor[i][j] = prefactor_one;
       cut[i][j] = cut_one;
       setflag[i][j] = 1;
       count++;
@@ -239,13 +209,11 @@ double PairSoft::init_one(int i, int j)
   // always mix prefactors geometrically
 
   if (setflag[i][j] == 0) {
-    prestart[i][j] = sqrt(prestart[i][i]*prestart[j][j]);
-    prestop[i][j] = sqrt(prestop[i][i]*prestop[j][j]);
+    prefactor[i][j] = sqrt(prefactor[i][i]*prefactor[j][j]);
     cut[i][j] = mix_distance(cut[i][i],cut[j][j]);
   }
 
-  prestart[j][i] = prestart[i][j];
-  prestop[j][i] = prestop[i][j];
+  prefactor[j][i] = prefactor[i][j];
   cut[j][i] = cut[i][j];
 
   return cut[i][j];
@@ -264,8 +232,7 @@ void PairSoft::write_restart(FILE *fp)
     for (j = i; j <= atom->ntypes; j++) {
       fwrite(&setflag[i][j],sizeof(int),1,fp);
       if (setflag[i][j]) {
-	fwrite(&prestart[i][j],sizeof(double),1,fp);
-	fwrite(&prestop[i][j],sizeof(double),1,fp);
+	fwrite(&prefactor[i][j],sizeof(double),1,fp);
 	fwrite(&cut[i][j],sizeof(double),1,fp);
       }
     }
@@ -289,12 +256,10 @@ void PairSoft::read_restart(FILE *fp)
       MPI_Bcast(&setflag[i][j],1,MPI_INT,0,world);
       if (setflag[i][j]) {
 	if (me == 0) {
-	  fread(&prestart[i][j],sizeof(double),1,fp);
-	  fread(&prestop[i][j],sizeof(double),1,fp);
+	  fread(&prefactor[i][j],sizeof(double),1,fp);
 	  fread(&cut[i][j],sizeof(double),1,fp);
 	}
-	MPI_Bcast(&prestart[i][j],1,MPI_DOUBLE,0,world);
-	MPI_Bcast(&prestop[i][j],1,MPI_DOUBLE,0,world);
+	MPI_Bcast(&prefactor[i][j],1,MPI_DOUBLE,0,world);
 	MPI_Bcast(&cut[i][j],1,MPI_DOUBLE,0,world);
       }
     }
@@ -324,6 +289,38 @@ void PairSoft::read_restart_settings(FILE *fp)
   MPI_Bcast(&mix_flag,1,MPI_INT,0,world);
 }
 
+/* ----------------------------------------------------------------------
+   check if name is recognized, return integer index for that name
+   if name not recognized, return -1
+   if type pair setting, return -2 if no type pairs are set
+------------------------------------------------------------------------- */
+
+int PairSoft::pre_adapt(char *name, int ilo, int ihi, int jlo, int jhi)
+{
+  int count = 0;
+  for (int i = ilo; i <= ihi; i++)
+    for (int j = MAX(jlo,i); j <= jhi; j++)
+      count++;
+  if (count == 0) return -2;
+
+  if (strcmp(name,"a") == 0) return 0;
+  return -1;
+}
+
+/* ----------------------------------------------------------------------
+   adapt parameter indexed by which
+   change all pair variables affected by the reset parameter
+   if type pair setting, set I-J and J-I coeffs
+------------------------------------------------------------------------- */
+
+void PairSoft::adapt(int which, int ilo, int ihi, int jlo, int jhi,
+		     double value)
+{
+  for (int i = ilo; i <= ihi; i++)
+    for (int j = MAX(jlo,i); j <= jhi; j++)
+      prefactor[i][j] = prefactor[j][i] = value;
+}
+
 /* ---------------------------------------------------------------------- */
 
 double PairSoft::single(int i, int j, int itype, int jtype, double rsq,
diff --git a/src/pair_soft.h b/src/pair_soft.h
index 026c351f64..eceb98ec6f 100644
--- a/src/pair_soft.h
+++ b/src/pair_soft.h
@@ -38,12 +38,13 @@ class PairSoft : public Pair {
   void read_restart(FILE *);
   void write_restart_settings(FILE *);
   void read_restart_settings(FILE *);
+  int pre_adapt(char *, int, int, int, int);
+  void adapt(int, int, int, int, int, double);
   double single(int, int, int, int, double, double, double, double &);
 
  private:
   double PI;
   double cut_global;
-  double **prestart,**prestop;
   double **prefactor;
   double **cut;
 

From d8b96e90270335b9bcc9771a54d6fdeb2a101bb7 Mon Sep 17 00:00:00 2001
From: sjplimp <sjplimp@f3b2605a-c512-4ea7-a41b-209d697bcdaa>
Date: Fri, 18 Jun 2010 16:33:24 +0000
Subject: [PATCH 06/20] git-svn-id:
 svn://svn.icms.temple.edu/lammps-ro/trunk@4298
 f3b2605a-c512-4ea7-a41b-209d697bcdaa

---
 src/COLLOID/pair_lubricate.cpp | 24 ++++++++++++++++++++++++
 src/COLLOID/pair_lubricate.h   |  2 ++
 2 files changed, 26 insertions(+)

diff --git a/src/COLLOID/pair_lubricate.cpp b/src/COLLOID/pair_lubricate.cpp
index b17be877ec..69a207430b 100644
--- a/src/COLLOID/pair_lubricate.cpp
+++ b/src/COLLOID/pair_lubricate.cpp
@@ -533,3 +533,27 @@ void PairLubricate::read_restart_settings(FILE *fp)
   delete random;
   random = new RanMars(lmp,seed + comm->me);
 }
+
+/* ----------------------------------------------------------------------
+   check if name is recognized, return integer index for that name
+   if name not recognized, return -1
+   if type pair setting, return -2 if no type pairs are set
+------------------------------------------------------------------------- */
+
+int PairLubricate::pre_adapt(char *name, int ilo, int ihi, int jlo, int jhi)
+{
+  if (strcmp(name,"mu") == 0) return 0;
+  return -1;
+}
+
+/* ----------------------------------------------------------------------
+   adapt parameter indexed by which
+   change all pair variables affected by the reset parameter
+   if type pair setting, set I-J and J-I coeffs
+------------------------------------------------------------------------- */
+
+void PairLubricate::adapt(int which, int ilo, int ihi, int jlo, int jhi,
+			  double value)
+{
+  mu = value;
+}
diff --git a/src/COLLOID/pair_lubricate.h b/src/COLLOID/pair_lubricate.h
index b4e76c2b6a..5206c6d5fc 100644
--- a/src/COLLOID/pair_lubricate.h
+++ b/src/COLLOID/pair_lubricate.h
@@ -37,6 +37,8 @@ class PairLubricate : public Pair {
   void read_restart(FILE *);
   void write_restart_settings(FILE *);
   void read_restart_settings(FILE *);
+  int pre_adapt(char *, int, int, int, int);
+  void adapt(int, int, int, int, int, double);
 
  protected:
   double cut_inner_global,cut_global;

From aadd0209bf88bf7c1207b48f363b5b13e13daf78 Mon Sep 17 00:00:00 2001
From: sjplimp <sjplimp@f3b2605a-c512-4ea7-a41b-209d697bcdaa>
Date: Fri, 18 Jun 2010 16:57:46 +0000
Subject: [PATCH 07/20] git-svn-id:
 svn://svn.icms.temple.edu/lammps-ro/trunk@4299
 f3b2605a-c512-4ea7-a41b-209d697bcdaa

---
 examples/micelle/in.micelle | 9 +++++++--
 1 file changed, 7 insertions(+), 2 deletions(-)

diff --git a/examples/micelle/in.micelle b/examples/micelle/in.micelle
index 128e4539ef..841e55ab2c 100644
--- a/examples/micelle/in.micelle
+++ b/examples/micelle/in.micelle
@@ -13,20 +13,25 @@ read_data	data.micelle
 special_bonds	fene
 
 pair_style	soft 1.12246
-pair_coeff	* * 1.0 20.0 1.12246
+pair_coeff	* * 0.0 1.12246
 
 bond_style 	harmonic
 bond_coeff	1 50.0 0.75
 
 velocity	all create 0.45 2349852
 
+variable	prefactor equal 1.0+elapsed*(20.0-1.0)/1000
+
 fix		1 all nve
 fix		2 all temp/rescale 100 0.45 0.45 0.02 1.0
-fix		3 all enforce2d
+fix		3 all adapt 1 pair soft a * * prefactor
+fix		4 all enforce2d
 
 thermo		50
 run		1000
 
+unfix		3
+
 # Main run
 
 pair_style	lj/cut 2.5

From 2b95d21d58d2acf38d24df5ab147bcc756c37b84 Mon Sep 17 00:00:00 2001
From: sjplimp <sjplimp@f3b2605a-c512-4ea7-a41b-209d697bcdaa>
Date: Fri, 18 Jun 2010 17:09:50 +0000
Subject: [PATCH 08/20] git-svn-id:
 svn://svn.icms.temple.edu/lammps-ro/trunk@4300
 f3b2605a-c512-4ea7-a41b-209d697bcdaa

---
 src/version.h | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/version.h b/src/version.h
index 9aa8998568..064954de83 100644
--- a/src/version.h
+++ b/src/version.h
@@ -1 +1 @@
-#define LAMMPS_VERSION "19 Jun 2010"
+#define LAMMPS_VERSION "20 Jun 2010"

From 2110a29427480ec58d03e3077077767fb8ea007b Mon Sep 17 00:00:00 2001
From: sjplimp <sjplimp@f3b2605a-c512-4ea7-a41b-209d697bcdaa>
Date: Fri, 18 Jun 2010 17:17:47 +0000
Subject: [PATCH 09/20] git-svn-id:
 svn://svn.icms.temple.edu/lammps-ro/trunk@4302
 f3b2605a-c512-4ea7-a41b-209d697bcdaa

---
 src/version.h | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/version.h b/src/version.h
index 064954de83..0621eb2550 100644
--- a/src/version.h
+++ b/src/version.h
@@ -1 +1 @@
-#define LAMMPS_VERSION "20 Jun 2010"
+#define LAMMPS_VERSION "21 Jun 2010"

From b3039e9e2aee5788a18cd8d7eb53bd313fd1272c Mon Sep 17 00:00:00 2001
From: sjplimp <sjplimp@f3b2605a-c512-4ea7-a41b-209d697bcdaa>
Date: Fri, 18 Jun 2010 17:22:51 +0000
Subject: [PATCH 10/20] git-svn-id:
 svn://svn.icms.temple.edu/lammps-ro/trunk@4304
 f3b2605a-c512-4ea7-a41b-209d697bcdaa

---
 src/Makefile             | 2 +-
 src/fix_wall_reflect.cpp | 8 ++++----
 2 files changed, 5 insertions(+), 5 deletions(-)

diff --git a/src/Makefile b/src/Makefile
index d172e4dc14..3d4843a1bd 100755
--- a/src/Makefile
+++ b/src/Makefile
@@ -14,7 +14,7 @@ OBJ = 	$(SRC:.cpp=.o)
 # Package variables
 
 PACKAGE = asphere class2 colloid dipole dsmc gpu granular \
-	  kspace manybody meam molecule opt peri poems prd reax xtc
+	  kspace manybody meam molecule opt peri poems prd reax shock xtc
 
 PACKUSER = user-ackland user-atc user-cd-eam user-cg-cmm user-ewaldn \
 	   user-imd user-smd
diff --git a/src/fix_wall_reflect.cpp b/src/fix_wall_reflect.cpp
index 7e93c2da0e..26d35980e6 100644
--- a/src/fix_wall_reflect.cpp
+++ b/src/fix_wall_reflect.cpp
@@ -74,7 +74,7 @@ FixWallReflect::FixWallReflect(LAMMPS *lmp, int narg, char **arg) :
       if (iarg+2 > narg) error->all("Illegal fix wall/reflect command");
       if (strcmp(arg[iarg+1],"EDGE") == 0) {
 	yloflag = EDGE;
-	ylo = domain->boxlo[0];
+	ylo = domain->boxlo[1];
       } else if (strstr(arg[iarg+1],"v_") == arg[iarg+1]) {
 	yloflag = VARIABLE;
 	int n = strlen(&arg[iarg+1][2]) + 1;
@@ -89,7 +89,7 @@ FixWallReflect::FixWallReflect(LAMMPS *lmp, int narg, char **arg) :
       if (iarg+2 > narg) error->all("Illegal fix wall/reflect command");
       if (strcmp(arg[iarg+1],"EDGE") == 0) {
 	yhiflag = EDGE;
-	yhi = domain->boxhi[0];
+	yhi = domain->boxhi[1];
       } else if (strstr(arg[iarg+1],"v_") == arg[iarg+1]) {
 	yhiflag = VARIABLE;
 	int n = strlen(&arg[iarg+1][2]) + 1;
@@ -104,7 +104,7 @@ FixWallReflect::FixWallReflect(LAMMPS *lmp, int narg, char **arg) :
       if (iarg+2 > narg) error->all("Illegal fix wall/reflect command");
       if (strcmp(arg[iarg+1],"EDGE") == 0) {
 	zloflag = EDGE;
-	zlo = domain->boxlo[0];
+	zlo = domain->boxlo[2];
       } else if (strstr(arg[iarg+1],"v_") == arg[iarg+1]) {
 	zloflag = VARIABLE;
 	int n = strlen(&arg[iarg+1][2]) + 1;
@@ -119,7 +119,7 @@ FixWallReflect::FixWallReflect(LAMMPS *lmp, int narg, char **arg) :
       if (iarg+2 > narg) error->all("Illegal fix wall/reflect command");
       if (strcmp(arg[iarg+1],"EDGE") == 0) {
 	zhiflag = EDGE;
-	zhi = domain->boxhi[0];
+	zhi = domain->boxhi[2];
       } else if (strstr(arg[iarg+1],"v_") == arg[iarg+1]) {
 	zhiflag = VARIABLE;
 	int n = strlen(&arg[iarg+1][2]) + 1;

From 0038e68b2a84b61b2d2d89dfda496fb3ad49ccdd Mon Sep 17 00:00:00 2001
From: sjplimp <sjplimp@f3b2605a-c512-4ea7-a41b-209d697bcdaa>
Date: Fri, 18 Jun 2010 17:23:19 +0000
Subject: [PATCH 11/20] git-svn-id:
 svn://svn.icms.temple.edu/lammps-ro/trunk@4305
 f3b2605a-c512-4ea7-a41b-209d697bcdaa

---
 src/version.h | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/version.h b/src/version.h
index 0621eb2550..ae07ab14fe 100644
--- a/src/version.h
+++ b/src/version.h
@@ -1 +1 @@
-#define LAMMPS_VERSION "21 Jun 2010"
+#define LAMMPS_VERSION "22 Jun 2010"

From f8c2900a1d13bc7fe1bee60ca2a9e4ee92584f8d Mon Sep 17 00:00:00 2001
From: sjplimp <sjplimp@f3b2605a-c512-4ea7-a41b-209d697bcdaa>
Date: Fri, 18 Jun 2010 18:58:45 +0000
Subject: [PATCH 12/20] git-svn-id:
 svn://svn.icms.temple.edu/lammps-ro/trunk@4316
 f3b2605a-c512-4ea7-a41b-209d697bcdaa

---
 src/compute_property_local.cpp | 43 +++++++++++++++++++++++-----------
 src/compute_property_local.h   |  2 +-
 2 files changed, 30 insertions(+), 15 deletions(-)

diff --git a/src/compute_property_local.cpp b/src/compute_property_local.cpp
index 94e4142b8d..392bfef96e 100644
--- a/src/compute_property_local.cpp
+++ b/src/compute_property_local.cpp
@@ -26,7 +26,7 @@
 
 using namespace LAMMPS_NS;
 
-enum{NONE,PAIR,BOND,ANGLE,DIHEDRAL,IMPROPER};
+enum{NONE,NEIGH,PAIR,BOND,ANGLE,DIHEDRAL,IMPROPER};
 
 #define DELTA 10000
 
@@ -50,7 +50,18 @@ ComputePropertyLocal::ComputePropertyLocal(LAMMPS *lmp, int narg, char **arg) :
   for (int iarg = 3; iarg < narg; iarg++) {
     i = iarg-3;
 
-    if (strcmp(arg[iarg],"patom1") == 0) {
+    if (strcmp(arg[iarg],"natom1") == 0) {
+      pack_choice[i] = &ComputePropertyLocal::pack_patom1;
+      if (kindflag != NONE && kindflag != NEIGH) 
+	error->all("Compute property/local cannot use these inputs together");
+      kindflag = NEIGH;
+    } else if (strcmp(arg[iarg],"natom2") == 0) {
+      pack_choice[i] = &ComputePropertyLocal::pack_patom2;
+      if (kindflag != NONE && kindflag != NEIGH) 
+	error->all("Compute property/local cannot use these inputs together");
+      kindflag = NEIGH;
+
+    } else if (strcmp(arg[iarg],"patom1") == 0) {
       pack_choice[i] = &ComputePropertyLocal::pack_patom1;
       if (kindflag != NONE && kindflag != PAIR) 
 	error->all("Compute property/local cannot use these inputs together");
@@ -184,16 +195,16 @@ ComputePropertyLocal::~ComputePropertyLocal()
 
 void ComputePropertyLocal::init()
 {
-  if (kindflag == PAIR) {
+  if (kindflag == NEIGH || kindflag == PAIR) {
     if (force->pair == NULL) 
       error->all("No pair style is defined for compute property/local");
     if (force->pair->single_enable == 0)
       error->all("Pair style does not support compute property/local");
   }
 
-  // for PAIR< need an occasional half neighbor list
+  // for NEIGH/PAIR need an occasional half neighbor list
 
-  if (kindflag == PAIR) {
+  if (kindflag == NEIGH || kindflag == PAIR) {
     int irequest = neighbor->request((void *) this);
     neighbor->requests[irequest]->pair = 0;
     neighbor->requests[irequest]->compute = 1;
@@ -201,9 +212,10 @@ void ComputePropertyLocal::init()
   }
 
   // do initial memory allocation so that memory_usage() is correct
-  // cannot be done yet for PAIR, since neigh list does not exist
+  // cannot be done yet for NEIGH/PAIR, since neigh list does not exist
 
-  if (kindflag == PAIR) ncount = 0;
+  if (kindflag == NEIGH) ncount = 0;
+  else if (kindflag == PAIR) ncount = 0;
   else if (kindflag == BOND) ncount = count_bonds(0);
   else if (kindflag == ANGLE) ncount = count_angles(0);
   else if (kindflag == DIHEDRAL) ncount = count_dihedrals(0);
@@ -228,7 +240,8 @@ void ComputePropertyLocal::compute_local()
 
   // count local entries and generate list of indices
 
-  if (kindflag == PAIR) ncount = count_pairs(0);
+  if (kindflag == NEIGH) ncount = count_pairs(0,0);
+  else if (kindflag == PAIR) ncount = count_pairs(0,1);
   else if (kindflag == BOND) ncount = count_bonds(0);
   else if (kindflag == ANGLE) ncount = count_angles(0);
   else if (kindflag == DIHEDRAL) ncount = count_dihedrals(0);
@@ -237,7 +250,8 @@ void ComputePropertyLocal::compute_local()
   if (ncount > nmax) reallocate(ncount);
   size_local_rows = ncount;
 
-  if (kindflag == PAIR) ncount = count_pairs(1);
+  if (kindflag == NEIGH) ncount = count_pairs(1,0);
+  else if (kindflag == PAIR) ncount = count_pairs(1,1);
   else if (kindflag == BOND) ncount = count_bonds(1);
   else if (kindflag == ANGLE) ncount = count_angles(1);
   else if (kindflag == DIHEDRAL) ncount = count_dihedrals(1);
@@ -259,10 +273,11 @@ void ComputePropertyLocal::compute_local()
    count pairs and compute pair info on this proc
    only count pair once if newton_pair is off
    both atom I,J must be in group
-   if flag is set, compute requested info about pair
+   if allflag is set, compute requested info about pair
+   if forceflag = 1, pair must be within force cutoff, else neighbor cutoff
 ------------------------------------------------------------------------- */
 
-int ComputePropertyLocal::count_pairs(int flag)
+int ComputePropertyLocal::count_pairs(int allflag, int forceflag)
 {
   int i,j,m,n,ii,jj,inum,jnum,itype,jtype;
   double xtmp,ytmp,ztmp,delx,dely,delz;
@@ -281,7 +296,7 @@ int ComputePropertyLocal::count_pairs(int flag)
 
   // invoke half neighbor list (will copy or build if necessary)
 
-  if (flag == 0) neighbor->build_one(list->index);
+  if (allflag == 0) neighbor->build_one(list->index);
 
   inum = list->inum;
   ilist = list->ilist;
@@ -324,9 +339,9 @@ int ComputePropertyLocal::count_pairs(int flag)
       delz = ztmp - x[j][2];
       rsq = delx*delx + dely*dely + delz*delz;
       jtype = type[j];
-      if (rsq >= cutsq[itype][jtype]) continue;
+      if (forceflag && rsq >= cutsq[itype][jtype]) continue;
 
-      if (flag) {
+      if (allflag) {
 	indices[m][0] = tag[i];
 	indices[m][1] = tag[j];
       }
diff --git a/src/compute_property_local.h b/src/compute_property_local.h
index 59935acfe1..7517c674bc 100644
--- a/src/compute_property_local.h
+++ b/src/compute_property_local.h
@@ -46,7 +46,7 @@ class ComputePropertyLocal : public Compute {
   int ncount;
   int **indices;
 
-  int count_pairs(int);
+  int count_pairs(int, int);
   int count_bonds(int);
   int count_angles(int);
   int count_dihedrals(int);

From 8ffc7c464a9efec7ad9c31f13dfebfd90011237d Mon Sep 17 00:00:00 2001
From: sjplimp <sjplimp@f3b2605a-c512-4ea7-a41b-209d697bcdaa>
Date: Fri, 18 Jun 2010 18:58:59 +0000
Subject: [PATCH 13/20] git-svn-id:
 svn://svn.icms.temple.edu/lammps-ro/trunk@4317
 f3b2605a-c512-4ea7-a41b-209d697bcdaa

---
 doc/compute_property_local.html | 23 ++++++++++++++---------
 doc/compute_property_local.txt  | 23 ++++++++++++++---------
 2 files changed, 28 insertions(+), 18 deletions(-)

diff --git a/doc/compute_property_local.html b/doc/compute_property_local.html
index 0e8c0fae83..9583f4749d 100644
--- a/doc/compute_property_local.html
+++ b/doc/compute_property_local.html
@@ -21,13 +21,15 @@
 
 <LI>input = one or more attributes 
 
-<PRE>  possible attributes = patom1 patom2
+<PRE>  possible attributes = natom1 natom2
+		        patom1 patom2
                         batom1 batom2 btype
                         aatom1 aatom2 aatom3 atype
                         datom1 datom2 datom3 dtype
                         iatom1 iatom2 iatom3 itype 
 </PRE>
-<PRE>     patom1, patom2 = IDs of 2 atoms in each pair
+<PRE>     natom1, natom2 = IDs of 2 atoms in each pair (within neighbor cutoff)
+     patom1, patom2 = IDs of 2 atoms in each pair (within force cutoff)
      batom1, batom2 = IDs of 2 atoms in each bond
      btype = bond type of each bond
      aatom1, aatom2, aatom3 = IDs of 3 atoms in each angle
@@ -63,10 +65,13 @@ property/local command.
 <P>If the inputs are pair attributes, the local data is generated by
 looping over the pairwise neighbor list.  Info about an individual
 pairwise interaction will only be included if both atoms in the pair
-are in the specified compute group, and if the current pairwise
-distance is less than the force cutoff distance for that interaction,
-as defined by the <A HREF = "pair_style.html">pair_style</A> and
-<A HREF = "pair_coeff.html">pair_coeff</A> commands.
+are in the specified compute group.  For <I>natom1</I> and <I>natom2</I>, all
+atom pairs in the neighbor list are considered (out to the neighbor
+cutoff = force cutoff + <A HREF = "neighbor.html">neighbor skin</A>).  For <I>patom1</I>
+and <I>patom2</I>, the distance between the atoms must be less than the
+force cutoff distance for that pair to be included, as defined by the
+<A HREF = "pair_style.html">pair_style</A> and <A HREF = "pair_coeff.html">pair_coeff</A>
+commands.
 </P>
 <P>If the inputs are bond, angle, etc attributes, the local data is
 generated by looping over all the atoms owned on a processor and
@@ -84,9 +89,9 @@ bond/local</A> command can be combined with bond
 atom indices from this command and output by the <A HREF = "dump.html">dump
 local</A> command in a consistent way.
 </P>
-<P>The <I>patom1</I> and <I>patom2</I> attributes refer to the atom IDs of the 2
-atoms in each pairwise interaction computed by the
-<A HREF = "pair_style.html">pair_style</A> command. 
+<P>The <I>natom1</I> and <I>natom2</I>, or <I>patom1</I> and <I>patom2</I> attributes refer
+to the atom IDs of the 2 atoms in each pairwise interaction computed
+by the <A HREF = "pair_style.html">pair_style</A> command.
 </P>
 <P>IMPORTANT NOTE: For pairs, if two atoms I,J are involved in 1-2, 1-3,
 1-4 interactions within the molecular topology, their pairwise
diff --git a/doc/compute_property_local.txt b/doc/compute_property_local.txt
index d0409b31fd..1baaedd069 100644
--- a/doc/compute_property_local.txt
+++ b/doc/compute_property_local.txt
@@ -15,13 +15,15 @@ compute ID group-ID property/local input1 input2 ... :pre
 ID, group-ID are documented in "compute"_compute.html command :ulb,l
 property/local = style name of this compute command :l
 input = one or more attributes :l
-  possible attributes = patom1 patom2
+  possible attributes = natom1 natom2
+		        patom1 patom2
                         batom1 batom2 btype
                         aatom1 aatom2 aatom3 atype
                         datom1 datom2 datom3 dtype
                         iatom1 iatom2 iatom3 itype :pre
 
-     patom1, patom2 = IDs of 2 atoms in each pair
+     natom1, natom2 = IDs of 2 atoms in each pair (within neighbor cutoff)
+     patom1, patom2 = IDs of 2 atoms in each pair (within force cutoff)
      batom1, batom2 = IDs of 2 atoms in each bond
      btype = bond type of each bond
      aatom1, aatom2, aatom3 = IDs of 3 atoms in each angle
@@ -56,10 +58,13 @@ property/local command.
 If the inputs are pair attributes, the local data is generated by
 looping over the pairwise neighbor list.  Info about an individual
 pairwise interaction will only be included if both atoms in the pair
-are in the specified compute group, and if the current pairwise
-distance is less than the force cutoff distance for that interaction,
-as defined by the "pair_style"_pair_style.html and
-"pair_coeff"_pair_coeff.html commands.
+are in the specified compute group.  For {natom1} and {natom2}, all
+atom pairs in the neighbor list are considered (out to the neighbor
+cutoff = force cutoff + "neighbor skin"_neighbor.html).  For {patom1}
+and {patom2}, the distance between the atoms must be less than the
+force cutoff distance for that pair to be included, as defined by the
+"pair_style"_pair_style.html and "pair_coeff"_pair_coeff.html
+commands.
 
 If the inputs are bond, angle, etc attributes, the local data is
 generated by looping over all the atoms owned on a processor and
@@ -77,9 +82,9 @@ bond/local"_compute_bond_local.html command can be combined with bond
 atom indices from this command and output by the "dump
 local"_dump.html command in a consistent way.
 
-The {patom1} and {patom2} attributes refer to the atom IDs of the 2
-atoms in each pairwise interaction computed by the
-"pair_style"_pair_style.html command. 
+The {natom1} and {natom2}, or {patom1} and {patom2} attributes refer
+to the atom IDs of the 2 atoms in each pairwise interaction computed
+by the "pair_style"_pair_style.html command.
 
 IMPORTANT NOTE: For pairs, if two atoms I,J are involved in 1-2, 1-3,
 1-4 interactions within the molecular topology, their pairwise

From 8190e79a1c067c405eafcf4d5f692067f3c78e05 Mon Sep 17 00:00:00 2001
From: sjplimp <sjplimp@f3b2605a-c512-4ea7-a41b-209d697bcdaa>
Date: Fri, 18 Jun 2010 19:00:23 +0000
Subject: [PATCH 14/20] git-svn-id:
 svn://svn.icms.temple.edu/lammps-ro/trunk@4318
 f3b2605a-c512-4ea7-a41b-209d697bcdaa

---
 src/version.h | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/version.h b/src/version.h
index ae07ab14fe..c88bed11d5 100644
--- a/src/version.h
+++ b/src/version.h
@@ -1 +1 @@
-#define LAMMPS_VERSION "22 Jun 2010"
+#define LAMMPS_VERSION "23 Jun 2010"

From 9600268de94736ee33f86da18cbb615935b1a91b Mon Sep 17 00:00:00 2001
From: sjplimp <sjplimp@f3b2605a-c512-4ea7-a41b-209d697bcdaa>
Date: Fri, 18 Jun 2010 19:49:12 +0000
Subject: [PATCH 15/20] git-svn-id:
 svn://svn.icms.temple.edu/lammps-ro/trunk@4320
 f3b2605a-c512-4ea7-a41b-209d697bcdaa

---
 src/atom.cpp                   |  2 +-
 src/change_box.cpp             |  1 +
 src/compute_erotate_sphere.cpp |  6 +++--
 src/compute_improper_local.cpp |  1 -
 src/compute_property_local.cpp |  5 +---
 src/compute_temp_profile.cpp   |  4 ---
 src/read_restart.cpp           |  4 ---
 src/region.cpp                 |  2 --
 src/region_block.cpp           |  1 -
 src/region_cone.cpp            | 46 +++++++++++++++++-----------------
 src/region_cylinder.cpp        | 22 ++++++++--------
 src/region_prism.cpp           |  2 +-
 src/special.cpp                |  1 -
 src/temper.cpp                 |  2 +-
 14 files changed, 43 insertions(+), 56 deletions(-)

diff --git a/src/atom.cpp b/src/atom.cpp
index ecebf7a52b..00334ca3b5 100644
--- a/src/atom.cpp
+++ b/src/atom.cpp
@@ -1336,7 +1336,7 @@ void Atom::first_reorder()
 
 void Atom::sort()
 {
-  int i,m,n,ix,iy,iz,ibin,empty,ndone;
+  int i,m,n,ix,iy,iz,ibin,empty;
 
   // re-setup sort bins if needed
 
diff --git a/src/change_box.cpp b/src/change_box.cpp
index 9b9d477c2a..2bb5ad45fe 100644
--- a/src/change_box.cpp
+++ b/src/change_box.cpp
@@ -38,6 +38,7 @@ void ChangeBox::command(int narg, char **arg)
   int style;
   if (strcmp(arg[0],"ortho") == 0) style = ORTHO;
   else if (strcmp(arg[0],"triclinic") == 0) style = TRICLINIC;
+  else error->all("Illegal change_box command");
 
   if (style == ORTHO && domain->triclinic == 0)
     error->all("Change_box operation is invalid");
diff --git a/src/compute_erotate_sphere.cpp b/src/compute_erotate_sphere.cpp
index 7e6b46528d..1558196cce 100644
--- a/src/compute_erotate_sphere.cpp
+++ b/src/compute_erotate_sphere.cpp
@@ -102,12 +102,14 @@ double ComputeERotateSphere::compute_scalar()
 		      omega[i][2]*omega[i][2]) * radius[i]*radius[i]*rmass[i];
     } else {
       for (i = 0; i < nlocal; i++)
-	if (mask[i] & groupbit)
+	if (mask[i] & groupbit) {
+	  itype = type[i];
 	  erotate += (omega[i][0]*omega[i][0] + omega[i][1]*omega[i][1] + 
 		      omega[i][2]*omega[i][2]) * 
 	    radius[i]*radius[i]*mass[itype];
+	}
     }
-
+    
   } else {
     if (rmass) {
       for (i = 0; i < nlocal; i++) 
diff --git a/src/compute_improper_local.cpp b/src/compute_improper_local.cpp
index 1ac5d84c93..00153b66ac 100644
--- a/src/compute_improper_local.cpp
+++ b/src/compute_improper_local.cpp
@@ -117,7 +117,6 @@ int ComputeImproperLocal::compute_impropers(int flag)
   int **improper_atom2 = atom->improper_atom2;
   int **improper_atom3 = atom->improper_atom3;
   int **improper_atom4 = atom->improper_atom4;
-  int **improper_type = atom->improper_type;
   int *tag = atom->tag;
   int *mask = atom->mask;
   int nlocal = atom->nlocal;
diff --git a/src/compute_property_local.cpp b/src/compute_property_local.cpp
index 392bfef96e..d7c814eafb 100644
--- a/src/compute_property_local.cpp
+++ b/src/compute_property_local.cpp
@@ -281,7 +281,7 @@ int ComputePropertyLocal::count_pairs(int allflag, int forceflag)
 {
   int i,j,m,n,ii,jj,inum,jnum,itype,jtype;
   double xtmp,ytmp,ztmp,delx,dely,delz;
-  double rsq,eng,fpair,factor_coul,factor_lj;
+  double rsq,factor_coul,factor_lj;
   int *ilist,*jlist,*numneigh,**firstneigh;
 
   double **x = atom->x;
@@ -306,7 +306,6 @@ int ComputePropertyLocal::count_pairs(int allflag, int forceflag)
   // loop over neighbors of my atoms
   // skip if I or J are not in group
 
-  Pair *pair = force->pair;
   double **cutsq = force->pair->cutsq;
 
   m = n = 0;
@@ -453,7 +452,6 @@ int ComputePropertyLocal::count_dihedrals(int flag)
   int **dihedral_atom2 = atom->dihedral_atom2;
   int **dihedral_atom3 = atom->dihedral_atom3;
   int **dihedral_atom4 = atom->dihedral_atom4;
-  int **dihedral_type = atom->dihedral_type;
   int *tag = atom->tag;
   int *mask = atom->mask;
   int nlocal = atom->nlocal;
@@ -497,7 +495,6 @@ int ComputePropertyLocal::count_impropers(int flag)
   int **improper_atom2 = atom->improper_atom2;
   int **improper_atom3 = atom->improper_atom3;
   int **improper_atom4 = atom->improper_atom4;
-  int **improper_type = atom->improper_type;
   int *tag = atom->tag;
   int *mask = atom->mask;
   int nlocal = atom->nlocal;
diff --git a/src/compute_temp_profile.cpp b/src/compute_temp_profile.cpp
index 7e98a8145f..6d3bb452a4 100644
--- a/src/compute_temp_profile.cpp
+++ b/src/compute_temp_profile.cpp
@@ -171,7 +171,6 @@ double ComputeTempProfile::compute_scalar()
 
   bin_average();
 
-  double **x = atom->x;
   double **v = atom->v;
   double *mass = atom->mass;
   double *rmass = atom->rmass;
@@ -215,7 +214,6 @@ void ComputeTempProfile::compute_vector()
 
   bin_average();
 
-  double **x = atom->x;
   double **v = atom->v;
   double *mass = atom->mass;
   double *rmass = atom->rmass;
@@ -277,7 +275,6 @@ void ComputeTempProfile::remove_bias(int i, double *v)
 
 void ComputeTempProfile::remove_bias_all()
 {
-  double **x = atom->x;
   double **v = atom->v;
   int *mask = atom->mask;
   int nlocal = atom->nlocal;
@@ -367,7 +364,6 @@ void ComputeTempProfile::bin_average()
 
   // sum each particle's velocity to appropriate bin
 
-  double **x = atom->x;
   double **v = atom->v;
   int *mask = atom->mask;
   int nlocal = atom->nlocal;
diff --git a/src/read_restart.cpp b/src/read_restart.cpp
index 33a3e01d3e..bf908caaba 100644
--- a/src/read_restart.cpp
+++ b/src/read_restart.cpp
@@ -481,10 +481,6 @@ void ReadRestart::header()
 	}
       }
 
-      int flag = 0;
-      if ((xperiodic != domain->xperiodic || yperiodic != domain->yperiodic ||
-	   zperiodic != domain->zperiodic)) flag = 1;
-
       domain->boundary[0][0] = boundary[0][0];
       domain->boundary[0][1] = boundary[0][1];
       domain->boundary[1][0] = boundary[1][0];
diff --git a/src/region.cpp b/src/region.cpp
index bb41c50f04..e67a2b71c5 100644
--- a/src/region.cpp
+++ b/src/region.cpp
@@ -181,8 +181,6 @@ int Region::dynamic_check()
 
 int Region::match(double x, double y, double z)
 {
-  double a[3],b[3],c[3],d[3];
-
   if (dynamic) {
     double delta = (update->ntimestep - time_origin) * dt;
     if (dynamic == VELOCITY) {
diff --git a/src/region_block.cpp b/src/region_block.cpp
index 21321d25c5..dada041049 100644
--- a/src/region_block.cpp
+++ b/src/region_block.cpp
@@ -197,7 +197,6 @@ int RegBlock::surface_interior(double *x, double cutoff)
 int RegBlock::surface_exterior(double *x, double cutoff)
 {
   double xp,yp,zp;
-  double delta;
 
   // x is far enough from block that there is no contact
   // x is interior to block
diff --git a/src/region_cone.cpp b/src/region_cone.cpp
index 0c45de6adb..c6b94c295c 100644
--- a/src/region_cone.cpp
+++ b/src/region_cone.cpp
@@ -262,18 +262,18 @@ int RegCone::surface_interior(double *x, double cutoff)
     }
 
   } else if (axis == 'y') {
-    delx = x[0] - c1;
-    delz = x[2] - c2;
-    r = sqrt(delx*delx + delz*delz);
+    del1 = x[0] - c1;
+    del2 = x[2] - c2;
+    r = sqrt(del1*del1 + del2*del2);
     currentradius = radiuslo + (x[1]-lo)*(radiushi-radiuslo)/(hi-lo);
 
-    // x is exterior to cone
+    // y is exterior to cone
     
     if (r > currentradius || x[1] < lo || x[1] > hi) return 0;
 
-    // x is interior to cone or on its surface
-    // surflo = pt on outer circle of bottom end plane, same dir as x vs axis
-    // surfhi = pt on outer circle of top end plane, same dir as x vs axis
+    // y is interior to cone or on its surface
+    // surflo = pt on outer circle of bottom end plane, same dir as y vs axis
+    // surfhi = pt on outer circle of top end plane, same dir as y vs axis
 
     if (r > 0.0) {
       surflo[0] = c1 + del1*radiuslo/r;
@@ -312,18 +312,18 @@ int RegCone::surface_interior(double *x, double cutoff)
     }
 
   } else {
-    delx = x[0] - c1;
-    dely = x[1] - c2;
-    r = sqrt(delx*delx + dely*dely);
+    del1 = x[0] - c1;
+    del2 = x[1] - c2;
+    r = sqrt(del1*del1 + del2*del2);
     currentradius = radiuslo + (x[2]-lo)*(radiushi-radiuslo)/(hi-lo);
 
     // x is exterior to cone
     
     if (r > currentradius || x[2] < lo || x[2] > hi) return 0;
 
-    // x is interior to cone or on its surface
-    // surflo = pt on outer circle of bottom end plane, same dir as x vs axis
-    // surfhi = pt on outer circle of top end plane, same dir as x vs axis
+    // z is interior to cone or on its surface
+    // surflo = pt on outer circle of bottom end plane, same dir as z vs axis
+    // surfhi = pt on outer circle of top end plane, same dir as z vs axis
 
     if (r > 0.0) {
       surflo[0] = c1 + del1*radiuslo/r;
@@ -427,17 +427,17 @@ int RegCone::surface_exterior(double *x, double cutoff)
     r = sqrt(del1*del1 + del2*del2);
     currentradius = radiuslo + (x[1]-lo)*(radiushi-radiuslo)/(hi-lo);
 
-    // x is far enough from cone that there is no contact
-    // x is interior to cone
+    // y is far enough from cone that there is no contact
+    // y is interior to cone
 
     if (r >= maxradius+cutoff || 
 	x[1] <= lo-cutoff || x[1] >= hi+cutoff) return 0;
     if (r < currentradius && x[1] > lo && x[1] < hi) return 0;
 
-    // x is exterior to cone or on its surface
-    // corner1234 = 4 corner pts of half trapezoid = cone surf in plane of x
+    // y is exterior to cone or on its surface
+    // corner1234 = 4 corner pts of half trapezoid = cone surf in plane of y
     // project x to 3 line segments in half trapezoid (4th is axis of cone)
-    // nearest = point on surface of cone that x is closest to
+    // nearest = point on surface of cone that y is closest to
     //           could be edge of cone
     // do not add contact point if r >= cutoff
 
@@ -472,17 +472,17 @@ int RegCone::surface_exterior(double *x, double cutoff)
     r = sqrt(del1*del1 + del2*del2);
     currentradius = radiuslo + (x[2]-lo)*(radiushi-radiuslo)/(hi-lo);
 
-    // x is far enough from cone that there is no contact
-    // x is interior to cone
+    // z is far enough from cone that there is no contact
+    // z is interior to cone
 
     if (r >= maxradius+cutoff || 
 	x[2] <= lo-cutoff || x[2] >= hi+cutoff) return 0;
     if (r < currentradius && x[2] > lo && x[2] < hi) return 0;
 
-    // x is exterior to cone or on its surface
-    // corner1234 = 4 corner pts of half trapezoid = cone surf in plane of x
+    // z is exterior to cone or on its surface
+    // corner1234 = 4 corner pts of half trapezoid = cone surf in plane of z
     // project x to 3 line segments in half trapezoid (4th is axis of cone)
-    // nearest = point on surface of cone that x is closest to
+    // nearest = point on surface of cone that z is closest to
     //           could be edge of cone
     // do not add contact point if r >= cutoff
 
diff --git a/src/region_cylinder.cpp b/src/region_cylinder.cpp
index 1c342052e1..58bb51f212 100644
--- a/src/region_cylinder.cpp
+++ b/src/region_cylinder.cpp
@@ -231,11 +231,11 @@ int RegCylinder::surface_interior(double *x, double cutoff)
     del2 = x[2] - c2;
     r = sqrt(del1*del1 + del2*del2);
 
-    // x is exterior to cylinder
+    // y is exterior to cylinder
     
     if (r > radius || x[1] < lo || x[1] > hi) return 0;
 
-    // x is interior to cylinder or on its surface
+    // y is interior to cylinder or on its surface
 
     delta = radius - r;
     if (delta < cutoff && r > 0.0) {
@@ -265,11 +265,11 @@ int RegCylinder::surface_interior(double *x, double cutoff)
     del2 = x[1] - c2;
     r = sqrt(del1*del1 + del2*del2);
 
-    // x is exterior to cylinder
+    // z is exterior to cylinder
     
     if (r > radius || x[2] < lo || x[2] > hi) return 0;
 
-    // x is interior to cylinder or on its surface
+    // z is interior to cylinder or on its surface
 
     delta = radius - r;
     if (delta < cutoff && r > 0.0) {
@@ -306,7 +306,7 @@ int RegCylinder::surface_interior(double *x, double cutoff)
 
 int RegCylinder::surface_exterior(double *x, double cutoff)
 {
-  double del1,del2,r,delta;
+  double del1,del2,r;
   double xp,yp,zp;
 
   if (axis == 'x') {
@@ -345,13 +345,13 @@ int RegCylinder::surface_exterior(double *x, double cutoff)
     del2 = x[2] - c2;
     r = sqrt(del1*del1 + del2*del2);
 
-    // x is far enough from cylinder that there is no contact
-    // x is interior to cylinder
+    // y is far enough from cylinder that there is no contact
+    // y is interior to cylinder
 
     if (r >= radius+cutoff || x[1] <= lo-cutoff || x[1] >= hi+cutoff) return 0;
     if (r < radius && x[1] > lo && x[1] < hi) return 0;
 
-    // x is exterior to cylinder or on its surface
+    // y is exterior to cylinder or on its surface
     // xp,yp,zp = point on surface of cylinder that x is closest to
     //            could be edge of cylinder
     // do not add contact point if r >= cutoff
@@ -376,13 +376,13 @@ int RegCylinder::surface_exterior(double *x, double cutoff)
     del2 = x[1] - c2;
     r = sqrt(del1*del1 + del2*del2);
 
-    // x is far enough from cylinder that there is no contact
-    // x is interior to cylinder
+    // z is far enough from cylinder that there is no contact
+    // z is interior to cylinder
 
     if (r >= radius+cutoff || x[2] <= lo-cutoff || x[2] >= hi+cutoff) return 0;
     if (r < radius && x[2] > lo && x[2] < hi) return 0;
 
-    // x is exterior to cylinder or on its surface
+    // z is exterior to cylinder or on its surface
     // xp,yp,zp = point on surface of cylinder that x is closest to
     //            could be edge of cylinder
     // do not add contact point if r >= cutoff
diff --git a/src/region_prism.cpp b/src/region_prism.cpp
index 0be3a83b94..26d65ab4f7 100644
--- a/src/region_prism.cpp
+++ b/src/region_prism.cpp
@@ -304,7 +304,7 @@ int RegPrism::surface_interior(double *x, double cutoff)
 int RegPrism::surface_exterior(double *x, double cutoff)
 {
   int i;
-  double dot,delta;
+  double dot;
   double *corner;
   double xp,yp,zp;
 
diff --git a/src/special.cpp b/src/special.cpp
index f786257efa..e23afc8005 100644
--- a/src/special.cpp
+++ b/src/special.cpp
@@ -705,7 +705,6 @@ void Special::dihedral_trim()
   MPI_Request request;
   MPI_Status status;
 
-  int *tag = atom->tag;
   int *num_dihedral = atom->num_dihedral;
   int **dihedral_atom1 = atom->dihedral_atom1;
   int **dihedral_atom4 = atom->dihedral_atom4;
diff --git a/src/temper.cpp b/src/temper.cpp
index ded3111764..ef4ba79d79 100644
--- a/src/temper.cpp
+++ b/src/temper.cpp
@@ -176,7 +176,7 @@ void Temper::command(int narg, char **arg)
   // setup tempering runs
 
   int i,which,partner,swap,partner_set_temp,partner_world;
-  double pe,pe_partner,ke,boltz_factor,new_temp;
+  double pe,pe_partner,boltz_factor,new_temp;
   MPI_Status status;
 
   if (me_universe == 0 && universe->uscreen) 

From 6d6a62213e1235e3320aa5f4fb84fcdca324b307 Mon Sep 17 00:00:00 2001
From: sjplimp <sjplimp@f3b2605a-c512-4ea7-a41b-209d697bcdaa>
Date: Fri, 18 Jun 2010 19:49:25 +0000
Subject: [PATCH 16/20] git-svn-id:
 svn://svn.icms.temple.edu/lammps-ro/trunk@4321
 f3b2605a-c512-4ea7-a41b-209d697bcdaa

---
 src/ASPHERE/compute_temp_asphere.cpp | 1 -
 src/COLLOID/pair_colloid.cpp         | 1 +
 src/MOLECULE/angle_table.cpp         | 8 ++++----
 src/MOLECULE/bond_table.cpp          | 4 ++--
 src/PRD/prd.cpp                      | 2 +-
 5 files changed, 8 insertions(+), 8 deletions(-)

diff --git a/src/ASPHERE/compute_temp_asphere.cpp b/src/ASPHERE/compute_temp_asphere.cpp
index c79fa3cf84..8f4af98aba 100755
--- a/src/ASPHERE/compute_temp_asphere.cpp
+++ b/src/ASPHERE/compute_temp_asphere.cpp
@@ -166,7 +166,6 @@ double ComputeTempAsphere::compute_scalar()
   double **quat = atom->quat;
   double **angmom = atom->angmom;
   double *mass = atom->mass;
-  double **shape = atom->shape;
   int *type = atom->type;
   int *mask = atom->mask;
   int nlocal = atom->nlocal;
diff --git a/src/COLLOID/pair_colloid.cpp b/src/COLLOID/pair_colloid.cpp
index 6d83e3bd07..ef9e8aa83b 100644
--- a/src/COLLOID/pair_colloid.cpp
+++ b/src/COLLOID/pair_colloid.cpp
@@ -471,6 +471,7 @@ double PairColloid::single(int i, int j, int itype, int jtype, double rsq,
     c2 = a2[itype][jtype];
     K[1] = c2*c2;
     K[2] = rsq;
+    K[0] = K[1] - rsq;
     K[4] = rsq*rsq;
     K[3] = K[1] - K[2];
     K[3] *= K[3]*K[3];
diff --git a/src/MOLECULE/angle_table.cpp b/src/MOLECULE/angle_table.cpp
index aab49790c8..1e6a5fcdfe 100644
--- a/src/MOLECULE/angle_table.cpp
+++ b/src/MOLECULE/angle_table.cpp
@@ -64,10 +64,10 @@ AngleTable::~AngleTable()
 
 void AngleTable::compute(int eflag, int vflag)
 {
-  int i1,i2,i3,n,type,factor;
+  int i1,i2,i3,n,type;
   double eangle,f1[3],f3[3];
   double delx1,dely1,delz1,delx2,dely2,delz2;
-  double rsq1,rsq2,r1,r2,c,s,a,a11,a12,a22,vx1,vx2,vy1,vy2,vz1,vz2;
+  double rsq1,rsq2,r1,r2,c,s,a,a11,a12,a22;
   double theta,u,mdu; //mdu: minus du, -du/dx=f
 
   eangle = 0.0;
@@ -620,7 +620,7 @@ double AngleTable::splint(double *xa, double *ya, double *y2a, int n, double x)
 void AngleTable::uf_lookup(int type, double x, double &u, double &f)
 {
   int itable;
-  double fraction,value,a,b;
+  double fraction,a,b;
 
   Table *tb = &tables[tabindex[type]];
   
@@ -651,7 +651,7 @@ void AngleTable::uf_lookup(int type, double x, double &u, double &f)
 void AngleTable::u_lookup(int type, double x, double &u)
 {
   int itable;
-  double fraction,value,a,b;
+  double fraction,a,b;
 
   Table *tb = &tables[tabindex[type]];
   
diff --git a/src/MOLECULE/bond_table.cpp b/src/MOLECULE/bond_table.cpp
index 8d0a7a11dd..7ffce66f1f 100644
--- a/src/MOLECULE/bond_table.cpp
+++ b/src/MOLECULE/bond_table.cpp
@@ -544,7 +544,7 @@ double BondTable::splint(double *xa, double *ya, double *y2a, int n, double x)
 void BondTable::uf_lookup(int type, double x, double &u, double &f)
 {
   int itable;
-  double fraction,value,a,b;
+  double fraction,a,b;
 
   Table *tb = &tables[tabindex[type]];
   x = MAX(x,tb->lo);
@@ -578,7 +578,7 @@ void BondTable::uf_lookup(int type, double x, double &u, double &f)
 void BondTable::u_lookup(int type, double x, double &u)
 {
   int itable;
-  double fraction,value,a,b;
+  double fraction,a,b;
 
   Table *tb = &tables[tabindex[type]];
   x = MAX(x,tb->lo);
diff --git a/src/PRD/prd.cpp b/src/PRD/prd.cpp
index 29cec75100..d9e19571fb 100644
--- a/src/PRD/prd.cpp
+++ b/src/PRD/prd.cpp
@@ -59,7 +59,7 @@ PRD::PRD(LAMMPS *lmp) : Pointers(lmp) {}
 
 void PRD::command(int narg, char **arg)
 {
-  int i,flag,allflag,ireplica;
+  int flag,ireplica;
 
   // error checks
 

From 79d83814e3d77818e414fffc6012ab33827250f0 Mon Sep 17 00:00:00 2001
From: sjplimp <sjplimp@f3b2605a-c512-4ea7-a41b-209d697bcdaa>
Date: Fri, 18 Jun 2010 20:09:44 +0000
Subject: [PATCH 17/20] git-svn-id:
 svn://svn.icms.temple.edu/lammps-ro/trunk@4322
 f3b2605a-c512-4ea7-a41b-209d697bcdaa

---
 src/compute_reduce.cpp | 113 ++++++++++++++++++++---------------------
 1 file changed, 54 insertions(+), 59 deletions(-)

diff --git a/src/compute_reduce.cpp b/src/compute_reduce.cpp
index 29386bb299..32ff36a97a 100644
--- a/src/compute_reduce.cpp
+++ b/src/compute_reduce.cpp
@@ -488,8 +488,8 @@ double ComputeReduce::compute_one(int m, int flag)
   // only include atoms in group for atom properties and per-atom quantities
 
   index = -1;
-  int n = value2index[m];
-  int j = argindex[m];
+  int vidx = value2index[m];
+  int aidx = argindex[m];
 
   int *mask = atom->mask;
   int nlocal = atom->nlocal;
@@ -504,50 +504,50 @@ double ComputeReduce::compute_one(int m, int flag)
     double **x = atom->x;
     if (flag < 0) {
       for (i = 0; i < nlocal; i++)
-	if (mask[i] & groupbit) combine(one,x[i][j],i);
-    } else one = x[flag][j];
+	if (mask[i] & groupbit) combine(one,x[i][aidx],i);
+    } else one = x[flag][aidx];
   } else if (which[m] == V) {
     double **v = atom->v;
     if (flag < 0) {
       for (i = 0; i < nlocal; i++)
-	if (mask[i] & groupbit) combine(one,v[i][j],i);
-    } else one = v[flag][j];
+	if (mask[i] & groupbit) combine(one,v[i][aidx],i);
+    } else one = v[flag][aidx];
   } else if (which[m] == F) {
     double **f = atom->f;
     if (flag < 0) {
       for (i = 0; i < nlocal; i++)
-	if (mask[i] & groupbit) combine(one,f[i][j],i);
-    } else one = f[flag][j];
+	if (mask[i] & groupbit) combine(one,f[i][aidx],i);
+    } else one = f[flag][aidx];
     
   // invoke compute if not previously invoked
 
   } else if (which[m] == COMPUTE) {
-    Compute *compute = modify->compute[n];
+    Compute *compute = modify->compute[vidx];
     
     if (flavor[m] == GLOBAL) {
-      if (j == 0) {
+      if (aidx == 0) {
 	if (!(compute->invoked_flag & INVOKED_VECTOR)) {
 	  compute->compute_vector();
 	  compute->invoked_flag |= INVOKED_VECTOR;
 	}
-	double *compute_vector = compute->vector;
+	double *comp_vec = compute->vector;
 	int n = compute->size_vector;
 	if (flag < 0) 
 	  for (i = 0; i < n; i++)
-	    combine(one,compute_vector[i],i);
-	else one = compute_vector[flag];
+	    combine(one,comp_vec[i],i);
+	else one = comp_vec[flag];
       } else {
 	if (!(compute->invoked_flag & INVOKED_ARRAY)) {
 	  compute->compute_array();
 	  compute->invoked_flag |= INVOKED_ARRAY;
 	}
-	double **compute_array = compute->array;
+	double **carray = compute->array;
 	int n = compute->size_array_rows;
-	int jm1 = j - 1;
+	int aidx_1 = aidx - 1;
 	if (flag < 0)
 	  for (i = 0; i < n; i++)
-	    combine(one,compute_array[i][jm1],i);
-	else one = compute_array[flag][jm1];
+	    combine(one,carray[i][aidx_1],i);
+	else one = carray[flag][aidx_1];
       }
 
     } else if (flavor[m] == PERATOM) {
@@ -556,21 +556,21 @@ double ComputeReduce::compute_one(int m, int flag)
 	compute->invoked_flag |= INVOKED_PERATOM;
       }
 
-      if (j == 0) {
-	double *compute_vector = compute->vector_atom;
+      if (aidx == 0) {
+	double *comp_vec = compute->vector_atom;
 	int n = nlocal;
 	if (flag < 0) {
 	  for (i = 0; i < n; i++)
-	    if (mask[i] & groupbit) combine(one,compute_vector[i],i);
-	} else one = compute_vector[flag];
+	    if (mask[i] & groupbit) combine(one,comp_vec[i],i);
+	} else one = comp_vec[flag];
       } else {
-	double **compute_array = compute->array_atom;
+	double **carray_atom = compute->array_atom;
 	int n = nlocal;
-	int jm1 = j - 1;
+	int aidxm1 = aidx - 1;
 	if (flag < 0) {
 	  for (i = 0; i < n; i++)
-	    if (mask[i] & groupbit) combine(one,compute_array[i][jm1],i);
-	} else one = compute_array[flag][jm1];
+	    if (mask[i] & groupbit) combine(one,carray_atom[i][aidxm1],i);
+	} else one = carray_atom[flag][aidxm1];
       }
 
     } else if (flavor[m] == LOCAL) {
@@ -579,49 +579,48 @@ double ComputeReduce::compute_one(int m, int flag)
 	compute->invoked_flag |= INVOKED_LOCAL;
       }
 
-      if (j == 0) {
-	double *compute_vector = compute->vector_local;
+      if (aidx == 0) {
+	double *comp_vec = compute->vector_local;
 	int n = compute->size_local_rows;
 	if (flag < 0)
 	  for (i = 0; i < n; i++)
-	    combine(one,compute_vector[i],i);
-	else one = compute_vector[flag];
+	    combine(one,comp_vec[i],i);
+	else one = comp_vec[flag];
       } else {
-	double **compute_array = compute->array_local;
+	double **carray_local = compute->array_local;
 	int n = compute->size_local_rows;
-	int jm1 = j - 1;
+	int aidxm1 = aidx - 1;
 	if (flag < 0)
 	  for (i = 0; i < n; i++)
-	    combine(one,compute_array[i][jm1],i);
-	else one = compute_array[flag][jm1];
+	    combine(one,carray_local[i][aidxm1],i);
+	else one = carray_local[flag][aidxm1];
       }
     }
 
   // access fix fields, check if fix frequency is a match
 
   } else if (which[m] == FIX) {
-    if (update->ntimestep % modify->fix[n]->peratom_freq)
+    if (update->ntimestep % modify->fix[vidx]->peratom_freq)
       error->all("Fix used in compute reduce not computed at compatible time");
-    Fix *fix = modify->fix[n];
+    Fix *fix = modify->fix[vidx];
 
     if (flavor[m] == GLOBAL) {
-      if (j == 0) {
+      if (aidx == 0) {
 	int n = fix->size_vector;
 	if (flag < 0)
 	  for (i = 0; i < n; i++)
 	    combine(one,fix->compute_vector(i),i);
 	else one = fix->compute_vector(flag);
       } else {
-	int n = fix->size_array_rows;
-	int jm1 = j - 1;
+	int aidxm1 = aidx - 1;
 	if (flag < 0)
 	  for (i = 0; i < nlocal; i++)
-	    combine(one,fix->compute_array(i,jm1),i);
-	else one = fix->compute_array(flag,jm1);
+	    combine(one,fix->compute_array(i,aidxm1),i);
+	else one = fix->compute_array(flag,aidxm1);
       }
 
     } else if (flavor[m] == PERATOM) {
-      if (j == 0) {
+      if (aidx == 0) {
 	double *fix_vector = fix->vector_atom;
 	int n = nlocal;
 	if (flag < 0) {
@@ -630,16 +629,15 @@ double ComputeReduce::compute_one(int m, int flag)
 	} else one = fix_vector[flag];
       } else {
 	double **fix_array = fix->array_atom;
-	int n = nlocal;
-	int jm1 = j - 1;
+	int aidxm1 = aidx - 1;
 	if (flag < 0) {
 	  for (i = 0; i < nlocal; i++)
-	    if (mask[i] & groupbit) combine(one,fix_array[i][jm1],i);
-	} else one = fix_array[flag][jm1];
+	    if (mask[i] & groupbit) combine(one,fix_array[i][aidxm1],i);
+	} else one = fix_array[flag][aidxm1];
       }
 
     } else if (flavor[m] == LOCAL) {
-      if (j == 0) {
+      if (aidx == 0) {
 	double *fix_vector = fix->vector_local;
 	int n = fix->size_local_rows;
 	if (flag < 0)
@@ -649,11 +647,11 @@ double ComputeReduce::compute_one(int m, int flag)
       } else {
 	double **fix_array = fix->array_local;
 	int n = fix->size_local_rows;
-	int jm1 = j - 1;
+	int aidxm1 = aidx - 1;
 	if (flag < 0)
 	  for (i = 0; i < n; i++)
-	    combine(one,fix_array[i][jm1],i);
-	else one = fix_array[flag][jm1];
+	    combine(one,fix_array[i][aidxm1],i);
+	else one = fix_array[flag][aidxm1];
       }
     }
     
@@ -667,7 +665,7 @@ double ComputeReduce::compute_one(int m, int flag)
 	memory->smalloc(maxatom*sizeof(double),"reduce:varatom");
     }
 
-    input->variable->compute_atom(n,igroup,varatom,1,0);
+    input->variable->compute_atom(vidx,igroup,varatom,1,0);
     if (flag < 0) {
       for (i = 0; i < nlocal; i++)
 	if (mask[i] & groupbit) combine(one,varatom[i],i);
@@ -681,18 +679,15 @@ double ComputeReduce::compute_one(int m, int flag)
 
 double ComputeReduce::count(int m)
 {
-  int n = value2index[m];
-  int j = argindex[m];
-
-  int *mask = atom->mask;
-  int nlocal = atom->nlocal;
+  int vidx = value2index[m];
+  int aidx = argindex[m];
 
   if (which[m] == X || which[m] == V || which[m] == F)
     return group->count(igroup);
   else if (which[m] == COMPUTE) {
-    Compute *compute = modify->compute[n];
+    Compute *compute = modify->compute[vidx];
     if (flavor[m] == GLOBAL) {
-      if (j == 0) return(1.0*compute->size_vector);
+      if (aidx == 0) return(1.0*compute->size_vector);
       else return(1.0*compute->size_array_rows);
     } else if (flavor[m] == PERATOM) {
       return group->count(igroup);
@@ -703,9 +698,9 @@ double ComputeReduce::count(int m)
       return ncountall;
     }
   } else if (which[m] == FIX) {
-    Fix *fix = modify->fix[n];
+    Fix *fix = modify->fix[vidx];
     if (flavor[m] == GLOBAL) {
-      if (j == 0) return(1.0*fix->size_vector);
+      if (aidx == 0) return(1.0*fix->size_vector);
       else return(1.0*fix->size_array_rows);
     } else if (flavor[m] == PERATOM) {
       return group->count(igroup);

From 4eda3e3688867050bd75f1c8e4b0f2155646d838 Mon Sep 17 00:00:00 2001
From: sjplimp <sjplimp@f3b2605a-c512-4ea7-a41b-209d697bcdaa>
Date: Fri, 18 Jun 2010 20:39:29 +0000
Subject: [PATCH 18/20] git-svn-id:
 svn://svn.icms.temple.edu/lammps-ro/trunk@4323
 f3b2605a-c512-4ea7-a41b-209d697bcdaa

---
 src/fix_aveforce.cpp | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/src/fix_aveforce.cpp b/src/fix_aveforce.cpp
index 913b82a4d2..7abece0d37 100644
--- a/src/fix_aveforce.cpp
+++ b/src/fix_aveforce.cpp
@@ -40,6 +40,8 @@ FixAveForce::FixAveForce(LAMMPS *lmp, int narg, char **arg) :
   global_freq = 1;
   extvector = 1;
 
+  xstr = ystr = zstr = NULL;
+
   if (strstr(arg[3],"v_") == arg[3]) {
     int n = strlen(&arg[3][2]) + 1;
     xstr = new char[n];

From c71f0cf9f8a1d27b0bc1ec5b830e3257950d1c26 Mon Sep 17 00:00:00 2001
From: sjplimp <sjplimp@f3b2605a-c512-4ea7-a41b-209d697bcdaa>
Date: Fri, 18 Jun 2010 20:40:43 +0000
Subject: [PATCH 19/20] git-svn-id:
 svn://svn.icms.temple.edu/lammps-ro/trunk@4324
 f3b2605a-c512-4ea7-a41b-209d697bcdaa

---
 src/version.h | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/version.h b/src/version.h
index c88bed11d5..62e2d576b7 100644
--- a/src/version.h
+++ b/src/version.h
@@ -1 +1 @@
-#define LAMMPS_VERSION "23 Jun 2010"
+#define LAMMPS_VERSION "24 Jun 2010"

From b3de0db9db300e7e3fc0e9b284e9b825603825af Mon Sep 17 00:00:00 2001
From: sjplimp <sjplimp@f3b2605a-c512-4ea7-a41b-209d697bcdaa>
Date: Fri, 18 Jun 2010 21:54:42 +0000
Subject: [PATCH 20/20] git-svn-id:
 svn://svn.icms.temple.edu/lammps-ro/trunk@4326
 f3b2605a-c512-4ea7-a41b-209d697bcdaa

---
 src/MOLECULE/angle_table.cpp   | 38 +++++++--------
 src/MOLECULE/angle_table.h     |  2 +-
 src/MOLECULE/bond_table.cpp    | 38 +++++++--------
 src/MOLECULE/bond_table.h      |  2 +-
 src/compute_dihedral_local.cpp | 69 ++++++++++-----------------
 src/pair_table.cpp             | 85 +++++++++++++++++-----------------
 src/pair_table.h               |  2 +-
 7 files changed, 109 insertions(+), 127 deletions(-)

diff --git a/src/MOLECULE/angle_table.cpp b/src/MOLECULE/angle_table.cpp
index 1e6a5fcdfe..034db79e7e 100644
--- a/src/MOLECULE/angle_table.cpp
+++ b/src/MOLECULE/angle_table.cpp
@@ -189,8 +189,8 @@ void AngleTable::settings(int narg, char **arg)
   else if (strcmp(arg[0],"spline") == 0) tabstyle = SPLINE;
   else error->all("Unknown table style in angle style table");
 
-  n = atoi(arg[1]);
-  nm1 = n - 1;
+  tablength = force->inumeric(arg[1]);
+  if (tablength < 2) error->all("Illegal number of angle table entries");
 
   // delete old tables, since cannot just change settings
 
@@ -282,7 +282,7 @@ double AngleTable::equilibrium_angle(int i)
 void AngleTable::write_restart(FILE *fp)
 {
   fwrite(&tabstyle,sizeof(int),1,fp);
-  fwrite(&n,sizeof(int),1,fp);
+  fwrite(&tablength,sizeof(int),1,fp);
 }
 
 /* ----------------------------------------------------------------------
@@ -293,11 +293,10 @@ void AngleTable::read_restart(FILE *fp)
 {
   if (comm->me == 0) {
     fread(&tabstyle,sizeof(int),1,fp);
-    fread(&n,sizeof(int),1,fp);
+    fread(&tablength,sizeof(int),1,fp);
   }
   MPI_Bcast(&tabstyle,1,MPI_DOUBLE,0,world);
-  MPI_Bcast(&n,1,MPI_INT,0,world);
-  nm1 = n - 1;
+  MPI_Bcast(&tablength,1,MPI_INT,0,world);
 
   allocate();
 }
@@ -451,7 +450,8 @@ void AngleTable::compute_table(Table *tb)
 {
   // delta = table spacing in angle for N-1 bins
 
-  tb->delta = PI/ nm1;
+  int tlm1 = tablength-1;
+  tb->delta = PI/ tlm1;
   tb->invdelta = 1.0/tb->delta;
   tb->deltasq6 = tb->delta*tb->delta / 6.0;
   
@@ -460,31 +460,31 @@ void AngleTable::compute_table(Table *tb)
   // de,df values = delta values of e,f
   // ang,e,f are N in length so de,df arrays can compute difference
 
-  tb->ang = (double *) memory->smalloc(n*sizeof(double),"angle:ang");
-  tb->e = (double *) memory->smalloc(n*sizeof(double),"angle:e");
-  tb->de = (double *) memory->smalloc(nm1*sizeof(double),"angle:de");
-  tb->f = (double *) memory->smalloc(n*sizeof(double),"angle:f");
-  tb->df = (double *) memory->smalloc(nm1*sizeof(double),"angle:df");
-  tb->e2 = (double *) memory->smalloc(n*sizeof(double),"angle:e2");
-  tb->f2 = (double *) memory->smalloc(n*sizeof(double),"angle:f2");
+  tb->ang = (double *) memory->smalloc(tablength*sizeof(double),"angle:ang");
+  tb->e = (double *) memory->smalloc(tablength*sizeof(double),"angle:e");
+  tb->de = (double *) memory->smalloc(tlm1*sizeof(double),"angle:de");
+  tb->f = (double *) memory->smalloc(tablength*sizeof(double),"angle:f");
+  tb->df = (double *) memory->smalloc(tlm1*sizeof(double),"angle:df");
+  tb->e2 = (double *) memory->smalloc(tablength*sizeof(double),"angle:e2");
+  tb->f2 = (double *) memory->smalloc(tablength*sizeof(double),"angle:f2");
 
   double a;
-  for (int i = 0; i < n; i++) {
+  for (int i = 0; i < tablength; i++) {
     a = i*tb->delta;
     tb->ang[i] = a;
 	  tb->e[i] = splint(tb->afile,tb->efile,tb->e2file,tb->ninput,a);
 	  tb->f[i] = splint(tb->afile,tb->ffile,tb->f2file,tb->ninput,a);
   }
         
-  for (int i = 0; i < nm1; i++) {
+  for (int i = 0; i < tlm1; i++) {
     tb->de[i] = tb->e[i+1] - tb->e[i];
     tb->df[i] = tb->f[i+1] - tb->f[i];
   }
      
   double ep0 = - tb->f[0];
-  double epn = - tb->f[nm1];
-  spline(tb->ang,tb->e,n,ep0,epn,tb->e2);  
-  spline(tb->ang,tb->f,n,tb->fplo,tb->fphi,tb->f2);
+  double epn = - tb->f[tlm1];
+  spline(tb->ang,tb->e,tablength,ep0,epn,tb->e2);  
+  spline(tb->ang,tb->f,tablength,tb->fplo,tb->fphi,tb->f2);
 }
 
 /* ----------------------------------------------------------------------
diff --git a/src/MOLECULE/angle_table.h b/src/MOLECULE/angle_table.h
index a4895196d8..452ae13e13 100644
--- a/src/MOLECULE/angle_table.h
+++ b/src/MOLECULE/angle_table.h
@@ -38,7 +38,7 @@ class AngleTable : public Angle {
   double single(int, int, int, int);
 
  private:
-  int tabstyle,n,nm1;
+  int tabstyle,tablength;
   double *theta0;
 
   struct Table {
diff --git a/src/MOLECULE/bond_table.cpp b/src/MOLECULE/bond_table.cpp
index 7ffce66f1f..28e28f17f2 100644
--- a/src/MOLECULE/bond_table.cpp
+++ b/src/MOLECULE/bond_table.cpp
@@ -141,8 +141,8 @@ void BondTable::settings(int narg, char **arg)
   else if (strcmp(arg[0],"spline") == 0) tabstyle = SPLINE;
   else error->all("Unknown table style in bond style table");
 
-  n = atoi(arg[1]);
-  nm1 = n - 1;
+  tablength = force->inumeric(arg[1]);
+  if (tablength < 2) error->all("Illegal number of bond table entries");
 
   // delete old tables, since cannot just change settings
 
@@ -224,7 +224,7 @@ double BondTable::equilibrium_distance(int i)
 void BondTable::write_restart(FILE *fp)
 {
   fwrite(&tabstyle,sizeof(int),1,fp);
-  fwrite(&n,sizeof(int),1,fp);
+  fwrite(&tablength,sizeof(int),1,fp);
 }
 
 /* ----------------------------------------------------------------------
@@ -235,11 +235,10 @@ void BondTable::read_restart(FILE *fp)
 {
   if (comm->me == 0) {
     fread(&tabstyle,sizeof(int),1,fp);
-    fread(&n,sizeof(int),1,fp);
+    fread(&tablength,sizeof(int),1,fp);
   }
   MPI_Bcast(&tabstyle,1,MPI_DOUBLE,0,world);
-  MPI_Bcast(&n,1,MPI_INT,0,world);
-  nm1 = n - 1;
+  MPI_Bcast(&tablength,1,MPI_INT,0,world);
 
   allocate();
 }
@@ -374,8 +373,9 @@ void BondTable::spline_table(Table *tb)
 void BondTable::compute_table(Table *tb)
 {
   // delta = table spacing for N-1 bins
+  int tlm1 = tablength-1;
 
-  tb->delta = (tb->hi - tb->lo)/ nm1;
+  tb->delta = (tb->hi - tb->lo)/ tlm1;
   tb->invdelta = 1.0/tb->delta;
   tb->deltasq6 = tb->delta*tb->delta / 6.0;
   
@@ -384,31 +384,31 @@ void BondTable::compute_table(Table *tb)
   // de,df values = delta values of e,f
   // r,e,f are N in length so de,df arrays can compute difference
 
-  tb->r = (double *) memory->smalloc(n*sizeof(double),"bond:r");
-  tb->e = (double *) memory->smalloc(n*sizeof(double),"bond:e");
-  tb->de = (double *) memory->smalloc(nm1*sizeof(double),"bond:de");
-  tb->f = (double *) memory->smalloc(n*sizeof(double),"bond:f");
-  tb->df = (double *) memory->smalloc(nm1*sizeof(double),"bond:df");
-  tb->e2 = (double *) memory->smalloc(n*sizeof(double),"bond:e2");
-  tb->f2 = (double *) memory->smalloc(n*sizeof(double),"bond:f2");
+  tb->r = (double *) memory->smalloc(tablength*sizeof(double),"bond:r");
+  tb->e = (double *) memory->smalloc(tablength*sizeof(double),"bond:e");
+  tb->de = (double *) memory->smalloc(tlm1*sizeof(double),"bond:de");
+  tb->f = (double *) memory->smalloc(tablength*sizeof(double),"bond:f");
+  tb->df = (double *) memory->smalloc(tlm1*sizeof(double),"bond:df");
+  tb->e2 = (double *) memory->smalloc(tablength*sizeof(double),"bond:e2");
+  tb->f2 = (double *) memory->smalloc(tablength*sizeof(double),"bond:f2");
 
   double a;
-  for (int i = 0; i < n; i++) {
+  for (int i = 0; i < tablength; i++) {
     a = tb->lo + i*tb->delta;
     tb->r[i] = a;
     tb->e[i] = splint(tb->rfile,tb->efile,tb->e2file,tb->ninput,a);
     tb->f[i] = splint(tb->rfile,tb->ffile,tb->f2file,tb->ninput,a);
   }
         
-  for (int i = 0; i < nm1; i++) {
+  for (int i = 0; i < tlm1; i++) {
     tb->de[i] = tb->e[i+1] - tb->e[i];
     tb->df[i] = tb->f[i+1] - tb->f[i];
   }
      
   double ep0 = - tb->f[0];
-  double epn = - tb->f[nm1];
-  spline(tb->r,tb->e,n,ep0,epn,tb->e2);  
-  spline(tb->r,tb->f,n,tb->fplo,tb->fphi,tb->f2);
+  double epn = - tb->f[tlm1];
+  spline(tb->r,tb->e,tablength,ep0,epn,tb->e2);  
+  spline(tb->r,tb->f,tablength,tb->fplo,tb->fphi,tb->f2);
 }
 
 /* ----------------------------------------------------------------------
diff --git a/src/MOLECULE/bond_table.h b/src/MOLECULE/bond_table.h
index dda3290d04..f0d289f055 100644
--- a/src/MOLECULE/bond_table.h
+++ b/src/MOLECULE/bond_table.h
@@ -38,7 +38,7 @@ class BondTable : public Bond {
   double single(int, double, int, int);
 
  private:
-  int tabstyle,n,nm1;
+  int tabstyle,tablength;
   double *r0;
 
   struct Table {
diff --git a/src/compute_dihedral_local.cpp b/src/compute_dihedral_local.cpp
index 1aa41e9a01..a8f2eddc50 100644
--- a/src/compute_dihedral_local.cpp
+++ b/src/compute_dihedral_local.cpp
@@ -109,9 +109,8 @@ int ComputeDihedralLocal::compute_dihedrals(int flag)
 {
   int i,m,n,atom1,atom2,atom3,atom4;
   double vb1x,vb1y,vb1z,vb2x,vb2y,vb2z,vb3x,vb3y,vb3z,vb2xm,vb2ym,vb2zm;
-  double sb1,sb2,sb3,rb1,rb3,c0,b1mag2,b1mag,b2mag2;
-  double b2mag,b3mag2,b3mag,ctmp,r12c1,c1mag,r12c2;
-  double c2mag,sin2,sc1,sc2,s1,s2,s12,c;
+  double ax,ay,az,bx,by,bz,rasq,rbsq,rgsq,rg,rginv,ra2inv,rb2inv,rabinv;
+  double s,c;
   double *pbuf;
 
   double **x = atom->x;
@@ -148,71 +147,53 @@ int ComputeDihedralLocal::compute_dihedrals(int flag)
 
       if (flag) {
 
-	// phi calculation from dihedral style OPLS
+	// phi calculation from dihedral style harmonic
 
 	if (pflag >= 0) {
 	  vb1x = x[atom1][0] - x[atom2][0];
 	  vb1y = x[atom1][1] - x[atom2][1];
 	  vb1z = x[atom1][2] - x[atom2][2];
 	  domain->minimum_image(vb1x,vb1y,vb1z);
-
+	  
 	  vb2x = x[atom3][0] - x[atom2][0];
 	  vb2y = x[atom3][1] - x[atom2][1];
 	  vb2z = x[atom3][2] - x[atom2][2];
 	  domain->minimum_image(vb2x,vb2y,vb2z);
-
+	  
 	  vb2xm = -vb2x;
 	  vb2ym = -vb2y;
 	  vb2zm = -vb2z;
 	  domain->minimum_image(vb2xm,vb2ym,vb2zm);
-
+	  
 	  vb3x = x[atom4][0] - x[atom3][0];
 	  vb3y = x[atom4][1] - x[atom3][1];
 	  vb3z = x[atom4][2] - x[atom3][2];
 	  domain->minimum_image(vb3x,vb3y,vb3z);
-
-	  sb1 = 1.0 / (vb1x*vb1x + vb1y*vb1y + vb1z*vb1z);
-	  sb2 = 1.0 / (vb2x*vb2x + vb2y*vb2y + vb2z*vb2z);
-	  sb3 = 1.0 / (vb3x*vb3x + vb3y*vb3y + vb3z*vb3z);
-
-	  rb1 = sqrt(sb1);
-	  rb3 = sqrt(sb3);
-
-	  c0 = (vb1x*vb3x + vb1y*vb3y + vb1z*vb3z) * rb1*rb3;
-
-	  b1mag2 = vb1x*vb1x + vb1y*vb1y + vb1z*vb1z;
-	  b1mag = sqrt(b1mag2);
-	  b2mag2 = vb2x*vb2x + vb2y*vb2y + vb2z*vb2z;
-	  b2mag = sqrt(b2mag2);
-	  b3mag2 = vb3x*vb3x + vb3y*vb3y + vb3z*vb3z;
-	  b3mag = sqrt(b3mag2);
 	  
-	  ctmp = vb1x*vb2x + vb1y*vb2y + vb1z*vb2z;
-	  r12c1 = 1.0 / (b1mag*b2mag);
-	  c1mag = ctmp * r12c1;
+	  ax = vb1y*vb2zm - vb1z*vb2ym;
+	  ay = vb1z*vb2xm - vb1x*vb2zm;
+	  az = vb1x*vb2ym - vb1y*vb2xm;
+	  bx = vb3y*vb2zm - vb3z*vb2ym;
+	  by = vb3z*vb2xm - vb3x*vb2zm;
+	  bz = vb3x*vb2ym - vb3y*vb2xm;
 	  
-	  ctmp = vb2xm*vb3x + vb2ym*vb3y + vb2zm*vb3z;
-	  r12c2 = 1.0 / (b2mag*b3mag);
-	  c2mag = ctmp * r12c2;
+	  rasq = ax*ax + ay*ay + az*az;
+	  rbsq = bx*bx + by*by + bz*bz;
+	  rgsq = vb2xm*vb2xm + vb2ym*vb2ym + vb2zm*vb2zm;
+	  rg = sqrt(rgsq);
+	  
+	  rginv = ra2inv = rb2inv = 0.0;
+	  if (rg > 0) rginv = 1.0/rg;
+	  if (rasq > 0) ra2inv = 1.0/rasq;
+	  if (rbsq > 0) rb2inv = 1.0/rbsq;
+	  rabinv = sqrt(ra2inv*rb2inv);
 
-	  sin2 = MAX(1.0 - c1mag*c1mag,0.0);
-	  sc1 = sqrt(sin2);
-	  if (sc1 < SMALL) sc1 = SMALL;
-	  sc1 = 1.0/sc1;
-
-	  sin2 = MAX(1.0 - c2mag*c2mag,0.0);
-	  sc2 = sqrt(sin2);
-	  if (sc2 < SMALL) sc2 = SMALL;
-	  sc2 = 1.0/sc2;
-
-	  s1 = sc1 * sc1;
-	  s2 = sc2 * sc2;
-	  s12 = sc1 * sc2;
-	  c = (c0 + c1mag*c2mag) * s12;
+	  c = (ax*bx + ay*by + az*bz)*rabinv;
+	  s = rg*rabinv*(ax*vb3x + ay*vb3y + az*vb3z);
 
 	  if (c > 1.0) c = 1.0;
 	  if (c < -1.0) c = -1.0;
-	  pbuf[n] = 180.0*acos(c)/PI;
+	  pbuf[n] = 180.0*atan2(s,c)/PI;
 	}
 	n += nvalues;
       }
diff --git a/src/pair_table.cpp b/src/pair_table.cpp
index c652a4148d..a116857108 100644
--- a/src/pair_table.cpp
+++ b/src/pair_table.cpp
@@ -76,6 +76,7 @@ void PairTable::compute(int eflag, int vflag)
   Table *tb;
   
   union_int_float_t rsq_lookup;
+  int tlm1 = tablength - 1;
 
   evdwl = 0.0;
   if (eflag || vflag) ev_setup(eflag,vflag);
@@ -127,19 +128,19 @@ void PairTable::compute(int eflag, int vflag)
  
 	if (tabstyle == LOOKUP) {
 	  itable = static_cast<int> ((rsq - tb->innersq) * tb->invdelta);
-	  if (itable >= nm1)
+	  if (itable >= tlm1)
 	    error->one("Pair distance > table outer cutoff");
 	  fpair = factor_lj * tb->f[itable];
 	} else if (tabstyle == LINEAR) {
 	  itable = static_cast<int> ((rsq - tb->innersq) * tb->invdelta);
-	  if (itable >= nm1)
+	  if (itable >= tlm1)
 	    error->one("Pair distance > table outer cutoff");
 	  fraction = (rsq - tb->rsq[itable]) * tb->invdelta;
 	  value = tb->f[itable] + fraction*tb->df[itable];
 	  fpair = factor_lj * value;
 	} else if (tabstyle == SPLINE) {
 	  itable = static_cast<int> ((rsq - tb->innersq) * tb->invdelta);
-	  if (itable >= nm1)
+	  if (itable >= tlm1)
 	    error->one("Pair distance > table outer cutoff");
 	  b = (rsq - tb->rsq[itable]) * tb->invdelta;
 	  a = 1.0 - b;
@@ -220,8 +221,8 @@ void PairTable::settings(int narg, char **arg)
   else if (strcmp(arg[0],"bitmap") == 0) tabstyle = BITMAP;
   else error->all("Unknown table style in pair_style command");
 
-  n = force->inumeric(arg[1]);
-  nm1 = n - 1;
+  tablength = force->inumeric(arg[1]);
+  if (tablength < 2) error->all("Illegal number of pair table entries");
 
   // delete old tables, since cannot just change settings
 
@@ -286,15 +287,13 @@ void PairTable::coeff(int narg, char **arg)
   // match = 1 if don't need to spline read-in tables
   // this is only the case if r values needed by final tables
   //   exactly match r values read from file
+  // for tabstyle SPLINE, always need to build spline tables
 
   tb->match = 0;
-  if (tabstyle == LINEAR && tb->ninput == n && 
+  if (tabstyle == LINEAR && tb->ninput == tablength && 
       tb->rflag == RSQ && tb->rhi == tb->cut) tb->match = 1;
-  if (tabstyle == SPLINE && tb->ninput == n && 
-      tb->rflag == RSQ && tb->rhi == tb->cut) tb->match = 1;
-  if (tabstyle == BITMAP && tb->ninput == 1 << n && 
+  if (tabstyle == BITMAP && tb->ninput == 1 << tablength && 
       tb->rflag == BMP && tb->rhi == tb->cut) tb->match = 1;
-
   if (tb->rflag == BMP && tb->match == 0)
     error->all("Bitmapped table in file does not match requested table");
 
@@ -537,6 +536,8 @@ void PairTable::param_extract(Table *tb, char *line)
 
 void PairTable::compute_table(Table *tb)
 {
+  int tlm1 = tablength-1;
+
   // inner = inner table bound
   // cut = outer table bound
   // delta = table spacing in rsq for N-1 bins
@@ -545,7 +546,7 @@ void PairTable::compute_table(Table *tb)
   if (tb->rflag) inner = tb->rlo;
   else inner = tb->rfile[0];
   tb->innersq = inner*inner;
-  tb->delta = (tb->cut*tb->cut - tb->innersq) / nm1;
+  tb->delta = (tb->cut*tb->cut - tb->innersq) / tlm1;
   tb->invdelta = 1.0/tb->delta;
 
   // direct lookup tables
@@ -556,11 +557,11 @@ void PairTable::compute_table(Table *tb)
   // e,f are never a match to read-in values, always computed via spline interp
 
   if (tabstyle == LOOKUP) {
-    tb->e = (double *) memory->smalloc(nm1*sizeof(double),"pair:e");
-    tb->f = (double *) memory->smalloc(nm1*sizeof(double),"pair:f");
+    tb->e = (double *) memory->smalloc(tlm1*sizeof(double),"pair:e");
+    tb->f = (double *) memory->smalloc(tlm1*sizeof(double),"pair:f");
 
     double r,rsq;
-    for (int i = 0; i < nm1; i++) {
+    for (int i = 0; i < tlm1; i++) {
       rsq = tb->innersq + (i+0.5)*tb->delta;
       r = sqrt(rsq);
       tb->e[i] = splint(tb->rfile,tb->efile,tb->e2file,tb->ninput,r);
@@ -577,14 +578,14 @@ void PairTable::compute_table(Table *tb)
   // e,f can match read-in values, else compute via spline interp
 
   if (tabstyle == LINEAR) {
-    tb->rsq = (double *) memory->smalloc(n*sizeof(double),"pair:rsq");
-    tb->e = (double *) memory->smalloc(n*sizeof(double),"pair:e");
-    tb->f = (double *) memory->smalloc(n*sizeof(double),"pair:f");
-    tb->de = (double *) memory->smalloc(nm1*sizeof(double),"pair:de");
-    tb->df = (double *) memory->smalloc(nm1*sizeof(double),"pair:df");
+    tb->rsq = (double *) memory->smalloc(tablength*sizeof(double),"pair:rsq");
+    tb->e = (double *) memory->smalloc(tablength*sizeof(double),"pair:e");
+    tb->f = (double *) memory->smalloc(tablength*sizeof(double),"pair:f");
+    tb->de = (double *) memory->smalloc(tlm1*sizeof(double),"pair:de");
+    tb->df = (double *) memory->smalloc(tlm1*sizeof(double),"pair:df");
 
     double r,rsq;
-    for (int i = 0; i < n; i++) {
+    for (int i = 0; i < tablength; i++) {
       rsq = tb->innersq + i*tb->delta;
       r = sqrt(rsq);
       tb->rsq[i] = rsq;
@@ -597,7 +598,7 @@ void PairTable::compute_table(Table *tb)
       }
     }
     
-    for (int i = 0; i < nm1; i++) {
+    for (int i = 0; i < tlm1; i++) {
       tb->de[i] = tb->e[i+1] - tb->e[i];
       tb->df[i] = tb->f[i+1] - tb->f[i];
     }
@@ -612,16 +613,16 @@ void PairTable::compute_table(Table *tb)
   // e,f can match read-in values, else compute via spline interp
 
   if (tabstyle == SPLINE) {
-    tb->rsq = (double *) memory->smalloc(n*sizeof(double),"pair:rsq");
-    tb->e = (double *) memory->smalloc(n*sizeof(double),"pair:e");
-    tb->f = (double *) memory->smalloc(n*sizeof(double),"pair:f");
-    tb->e2 = (double *) memory->smalloc(n*sizeof(double),"pair:e2");
-    tb->f2 = (double *) memory->smalloc(n*sizeof(double),"pair:f2");
+    tb->rsq = (double *) memory->smalloc(tablength*sizeof(double),"pair:rsq");
+    tb->e = (double *) memory->smalloc(tablength*sizeof(double),"pair:e");
+    tb->f = (double *) memory->smalloc(tablength*sizeof(double),"pair:f");
+    tb->e2 = (double *) memory->smalloc(tablength*sizeof(double),"pair:e2");
+    tb->f2 = (double *) memory->smalloc(tablength*sizeof(double),"pair:f2");
 
     tb->deltasq6 = tb->delta*tb->delta / 6.0;
 
     double r,rsq;
-    for (int i = 0; i < n; i++) {
+    for (int i = 0; i < tablength; i++) {
       rsq = tb->innersq + i*tb->delta;
       r = sqrt(rsq);
       tb->rsq[i] = rsq;
@@ -637,8 +638,8 @@ void PairTable::compute_table(Table *tb)
     // ep0,epn = dE/dr at inner and at cut
 
     double ep0 = - tb->f[0];
-    double epn = - tb->f[nm1];
-    spline(tb->rsq,tb->e,n,ep0,epn,tb->e2);
+    double epn = - tb->f[tlm1];
+    spline(tb->rsq,tb->e,tablength,ep0,epn,tb->e2);
 
     // fp0,fpn = dh/dg at inner and at cut
     // h(r) = f(r)/r and g(r) = r^2
@@ -657,17 +658,17 @@ void PairTable::compute_table(Table *tb)
     }
 
     if (tb->fpflag && tb->cut == tb->rfile[tb->ninput-1]) fpn =
-      (tb->fphi/tb->cut - tb->f[nm1]/(tb->cut*tb->cut)) / (2.0 * tb->cut);
+      (tb->fphi/tb->cut - tb->f[tlm1]/(tb->cut*tb->cut)) / (2.0 * tb->cut);
     else {
       double rsq2 = tb->cut * tb->cut;
       double rsq1 = rsq2 - secant_factor*tb->delta;
-      fpn = (tb->f[nm1] / sqrt(rsq2) - 
+      fpn = (tb->f[tlm1] / sqrt(rsq2) - 
 	     splint(tb->rfile,tb->ffile,tb->f2file,tb->ninput,sqrt(rsq1)) /
 	     sqrt(rsq1)) / (secant_factor*tb->delta);
     }
 
-    for (int i = 0; i < n; i++) tb->f[i] /= sqrt(tb->rsq[i]);
-    spline(tb->rsq,tb->f,n,fp0,fpn,tb->f2);
+    for (int i = 0; i < tablength; i++) tb->f[i] /= sqrt(tb->rsq[i]);
+    spline(tb->rsq,tb->f,tablength,fp0,fpn,tb->f2);
   }
 
   // bitmapped linear tables
@@ -683,8 +684,8 @@ void PairTable::compute_table(Table *tb)
     // linear lookup tables of length ntable = 2^n
     // stored value = value at lower edge of bin
 	
-    init_bitmap(inner,tb->cut,n,masklo,maskhi,tb->nmask,tb->nshiftbits);
-    int ntable = 1 << n;
+    init_bitmap(inner,tb->cut,tablength,masklo,maskhi,tb->nmask,tb->nshiftbits);
+    int ntable = 1 << tablength;
     int ntablem1 = ntable - 1;
 
     tb->rsq = (double *) memory->smalloc(ntable*sizeof(double),"pair:rsq");
@@ -885,7 +886,7 @@ void PairTable::read_restart(FILE *fp)
 void PairTable::write_restart_settings(FILE *fp)
 {
   fwrite(&tabstyle,sizeof(int),1,fp);
-  fwrite(&n,sizeof(int),1,fp);
+  fwrite(&tablength,sizeof(int),1,fp);
 }
 
 /* ----------------------------------------------------------------------
@@ -896,11 +897,10 @@ void PairTable::read_restart_settings(FILE *fp)
 {
   if (comm->me == 0) {
     fread(&tabstyle,sizeof(int),1,fp);
-    fread(&n,sizeof(int),1,fp);
+    fread(&tablength,sizeof(int),1,fp);
   }
   MPI_Bcast(&tabstyle,1,MPI_DOUBLE,0,world);
-  MPI_Bcast(&n,1,MPI_INT,0,world);
-  nm1 = n - 1;
+  MPI_Bcast(&tablength,1,MPI_INT,0,world);
 }
 
 /* ---------------------------------------------------------------------- */
@@ -911,23 +911,24 @@ double PairTable::single(int i, int j, int itype, int jtype, double rsq,
 {
   int itable;
   double fraction,value,a,b,phi;
+  int tlm1 = tablength - 1;
 
   Table *tb = &tables[tabindex[itype][jtype]];
   if (rsq < tb->innersq) error->one("Pair distance < table inner cutoff");
 
   if (tabstyle == LOOKUP) {
     itable = static_cast<int> ((rsq-tb->innersq) * tb->invdelta);
-    if (itable >= nm1) error->one("Pair distance > table outer cutoff");
+    if (itable >= tlm1) error->one("Pair distance > table outer cutoff");
     fforce = factor_lj * tb->f[itable];
   } else if (tabstyle == LINEAR) {
     itable = static_cast<int> ((rsq-tb->innersq) * tb->invdelta);
-    if (itable >= nm1) error->one("Pair distance > table outer cutoff");
+    if (itable >= tlm1) error->one("Pair distance > table outer cutoff");
     fraction = (rsq - tb->rsq[itable]) * tb->invdelta;
     value = tb->f[itable] + fraction*tb->df[itable];
     fforce = factor_lj * value;
   } else if (tabstyle == SPLINE) {
     itable = static_cast<int> ((rsq-tb->innersq) * tb->invdelta);
-    if (itable >= nm1) error->one("Pair distance > table outer cutoff");
+    if (itable >= tlm1) error->one("Pair distance > table outer cutoff");
     b = (rsq - tb->rsq[itable]) * tb->invdelta;
     a = 1.0 - b;
     value = a * tb->f[itable] + b * tb->f[itable+1] + 
diff --git a/src/pair_table.h b/src/pair_table.h
index d56c1b8d40..da38a4a62d 100644
--- a/src/pair_table.h
+++ b/src/pair_table.h
@@ -40,7 +40,7 @@ class PairTable : public Pair {
   void *extract(char *);
 
  private:
-  int tabstyle,n,nm1;
+  int tabstyle,tablength;
   struct Table {
     int ninput,rflag,fpflag,match,ntablebits;
     int nshiftbits,nmask;