Merge branch 'master' into lammps-icms

doc/Eqs/pair_mie.tex

@@ -0,0 +1,10 @@
+\documentstyle[12pt]{article}
+
+\begin{document}
+
+$$
+ E =  C \epsilon \left[ \left(\frac{\sigma}{r}\right)^{\gamma_{rep}} - \left(\frac{\sigma}{r}\right)^{\gamma_{att}} \right]
+                       \qquad r < r_c
+$$
+
+\end{document}

doc/Eqs/pair_mie2.tex

@@ -0,0 +1,9 @@
+\documentstyle[12pt]{article}
+
+\begin{document}
+
+$$
+ C = \left(\frac{\gamma_{rep}}{\gamma_{rep}-\gamma_{att}}\right) \left(\frac{\gamma_{rep}}{\gamma_{att}}\right)^{\left(\frac{\gamma_{att}}{\gamma_{rep}-\gamma_{att}}\right)}
+$$
+
+\end{document}

doc/Section_commands.html

@@ -433,10 +433,11 @@ potentials.  Click on the style itself for a full description:
 <TR ALIGN="center"><TD ><A HREF = "pair_lj_long.html">lj/long/coul/long</A></TD><TD ><A HREF = "pair_lj.html">lj/cut/tip4p/long</A></TD><TD ><A HREF = "pair_lj_long.html">lj/long/tip4p/long</A></TD><TD ><A HREF = "pair_lj_expand.html">lj/expand</A></TD></TR>
 <TR ALIGN="center"><TD ><A HREF = "pair_gromacs.html">lj/gromacs</A></TD><TD ><A HREF = "pair_gromacs.html">lj/gromacs/coul/gromacs</A></TD><TD ><A HREF = "pair_lj_smooth.html">lj/smooth</A></TD><TD ><A HREF = "pair_lj_smooth_linear.html">lj/smooth/linear</A></TD></TR>
 <TR ALIGN="center"><TD ><A HREF = "pair_lj96.html">lj96/cut</A></TD><TD ><A HREF = "pair_lubricate.html">lubricate</A></TD><TD ><A HREF = "pair_lubricate.html">lubricate/poly</A></TD><TD ><A HREF = "pair_lubricateU.html">lubricateU</A></TD></TR>
-<TR ALIGN="center"><TD ><A HREF = "pair_lubricateU.html">lubricateU/poly</A></TD><TD ><A HREF = "pair_meam.html">meam</A></TD><TD ><A HREF = "pair_morse.html">morse</A></TD><TD ><A HREF = "pair_peri.html">peri/lps</A></TD></TR>
+<TR ALIGN="center"><TD ><A HREF = "pair_lubricateU.html">lubricateU/poly</A></TD><TD ><A HREF = "pair_meam.html">meam</A></TD><TD ><A HREF = "pair_mie.html">mie/cut</A></TD><TD ><A HREF = "pair_morse.html">morse</A></TD></TR>
-<TR ALIGN="center"><TD ><A HREF = "pair_peri.html">peri/pmb</A></TD><TD ><A HREF = "pair_reax.html">reax</A></TD><TD ><A HREF = "pair_airebo.html">rebo</A></TD><TD ><A HREF = "pair_resquared.html">resquared</A></TD></TR>
+<TR ALIGN="center"><TD ><A HREF = "pair_peri.html">peri/lps</A></TD><TD ><A HREF = "pair_peri.html">peri/pmb</A></TD><TD ><A HREF = "pair_reax.html">reax</A></TD><TD ><A HREF = "pair_airebo.html">rebo</A></TD></TR>
-<TR ALIGN="center"><TD ><A HREF = "pair_soft.html">soft</A></TD><TD ><A HREF = "pair_sw.html">sw</A></TD><TD ><A HREF = "pair_table.html">table</A></TD><TD ><A HREF = "pair_tersoff.html">tersoff</A></TD></TR>
+<TR ALIGN="center"><TD ><A HREF = "pair_resquared.html">resquared</A></TD><TD ><A HREF = "pair_soft.html">soft</A></TD><TD ><A HREF = "pair_sw.html">sw</A></TD><TD ><A HREF = "pair_table.html">table</A></TD></TR>
-<TR ALIGN="center"><TD ><A HREF = "pair_tersoff_zbl.html">tersoff/zbl</A></TD><TD ><A HREF = "pair_tri_lj.html">tri/lj</A></TD><TD ><A HREF = "pair_yukawa.html">yukawa</A></TD><TD ><A HREF = "pair_yukawa_colloid.html">yukawa/colloid</A> 
+<TR ALIGN="center"><TD ><A HREF = "pair_tersoff.html">tersoff</A></TD><TD ><A HREF = "pair_tersoff_zbl.html">tersoff/zbl</A></TD><TD ><A HREF = "pair_tri_lj.html">tri/lj</A></TD><TD ><A HREF = "pair_yukawa.html">yukawa</A></TD></TR>
+<TR ALIGN="center"><TD ><A HREF = "pair_yukawa_colloid.html">yukawa/colloid</A> 
 </TD></TR></TABLE></DIV>
 
 <P>These are pair styles contributed by users, which can be used if

doc/Section_commands.txt

@@ -692,6 +692,7 @@ potentials.  Click on the style itself for a full description:
 "lubricateU"_pair_lubricateU.html,
 "lubricateU/poly"_pair_lubricateU.html,
 "meam"_pair_meam.html,
+"mie/cut"_pair_mie.html,
 "morse"_pair_morse.html,
 "peri/lps"_pair_peri.html,
 "peri/pmb"_pair_peri.html,

doc/dihedral_fourier.html

@@ -18,7 +18,7 @@
 <P><B>Examples:</B>
 </P>
 <PRE>dihedral_style fourier
-dihedral_coeff 3 -0.846200 3 0 7.578800 1 0 0.138000 2 -180 
+dihedral_coeff 3 -0.846200 3 0.0 7.578800 1 0 0.138000 2 -180.0 
 </PRE>
 <P><B>Description:</B>
 </P>
@@ -31,15 +31,15 @@ dihedral_coeff 3 -0.846200 3 0 7.578800 1 0 0.138000 2 -180
 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:
 </P>
-<P>m (integer >=1)
+<UL><LI>m (integer >=1)
-K1 (energy)
+<LI>K1 (energy)
-n1 (integer >= 0)
+<LI>n1 (integer >= 0)
-d1 (integer value of degrees)
+<LI>d1 (degrees)
-....
+<LI>....
-Km (energy)
+<LI>Km (energy)
-nm (integer >= 0)
+<LI>nm (integer >= 0)
-dm (integer value of degrees)
+<LI>dm (degrees) 
-</P>
+</UL>
 <HR>
 
 <P><B>Restrictions:</B>

doc/dihedral_fourier.txt

@@ -15,7 +15,7 @@ dihedral_style fourier :pre
 [Examples:]
 
 dihedral_style fourier
-dihedral_coeff 3 -0.846200 3 0 7.578800 1 0 0.138000 2 -180 :pre
+dihedral_coeff 3 -0.846200 3 0.0 7.578800 1 0 0.138000 2 -180.0 :pre
 
 [Description:]
 
@@ -31,11 +31,11 @@ or "read_restart"_read_restart.html commands:
 m (integer >=1)
 K1 (energy)
 n1 (integer >= 0)
-d1 (integer value of degrees)
+d1 (degrees)
 ....
 Km (energy)
 nm (integer >= 0)
-dm (integer value of degrees)
+dm (degrees) :ul
 
 :line
 

doc/dihedral_quadratic.html

@@ -33,9 +33,9 @@ dihedral_coeff 100.0 80.0
 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:
 </P>
-<P>K (energy)
+<UL><LI>K (energy)
-phi0 (degrees)
+<LI>phi0 (degrees) 
-</P>
+</UL>
 <HR>
 
 <P><B>Restrictions:</B>

doc/dihedral_quadratic.txt

@@ -31,7 +31,7 @@ the data file or restart files read by the "read_data"_read_data.html
 or "read_restart"_read_restart.html commands:
 
 K (energy)
-phi0 (degrees)
+phi0 (degrees) :ul
 
 :line
 

doc/dump.html

@@ -308,12 +308,14 @@ minimization converges.  Note that this means a dump will not be
 performed on the initial timestep after the dump command is invoked,
 if the current timestep is not a multiple of N.  This behavior can be
 changed via the <A HREF = "dump_modify.html">dump_modify first</A> command, which
-can be useful if the dump command is invoked after a minimization
+can also be useful if the dump command is invoked after a minimization
 ended on an arbitrary timestep.  N can be changed between runs by
 using the <A HREF = "dump_modify.html">dump_modify every</A> command (not allowed
 for <I>dcd</I> style).  The <A HREF = "dump_modify.html">dump_modify every</A> command
 also allows a variable to be used to determine the sequence of
-timesteps on which dump files are written.
+timesteps on which dump files are written.  In this mode a dump on the
+first timestep of a run will also not be written unless the
+<A HREF = "dump_modify.html">dump_modify first</A> command is used.
 </P>
 <P>The specified filename determines how the dump file(s) is written.
 The default is to write one large text file, which is opened when the

doc/dump.txt

@@ -296,12 +296,14 @@ minimization converges.  Note that this means a dump will not be
 performed on the initial timestep after the dump command is invoked,
 if the current timestep is not a multiple of N.  This behavior can be
 changed via the "dump_modify first"_dump_modify.html command, which
-can be useful if the dump command is invoked after a minimization
+can also be useful if the dump command is invoked after a minimization
 ended on an arbitrary timestep.  N can be changed between runs by
 using the "dump_modify every"_dump_modify.html command (not allowed
 for {dcd} style).  The "dump_modify every"_dump_modify.html command
 also allows a variable to be used to determine the sequence of
-timesteps on which dump files are written.
+timesteps on which dump files are written.  In this mode a dump on the
+first timestep of a run will also not be written unless the
+"dump_modify first"_dump_modify.html command is used.
 
 The specified filename determines how the dump file(s) is written.
 The default is to write one large text file, which is opened when the

doc/dump_modify.html

@@ -365,16 +365,22 @@ which should be specified as v_name, where name is the variable name.
 </P>
 <P>In this case, the variable is evaluated at the beginning of a run to
 determine the next timestep at which a dump snapshot will be written
-out.  On that timestep, the variable will be evaluated again to
+out.  On that timestep the variable will be evaluated again to
 determine the next timestep, etc.  Thus the variable should return
 timestep values.  See the stagger() and logfreq() and stride() math
 functions for <A HREF = "variable.html">equal-style variables</A>, as examples of
 useful functions to use in this context.  Other similar math functions
 could easily be added as options for <A HREF = "variable.html">equal-style
-variables</A>.  When using the variable option with the
+variables</A>.  Also see the next() function, which allows
-<I>every</I> keyword, you also need to use the <I>first</I> option if you want
+use of a file-style variable which reads successive values from a
-an initial snapshot written to the dump file.  The <I>every</I> keyword
+file, each time the variable is evaluated.  Used with the <I>every</I>
-cannot be used with the dump <I>dcd</I> style.
+keyword, if the file contains a list of ascending timesteps, you can
+output snapshots whenever you wish.
+</P>
+<P>Note that when using the variable option with the <I>every</I> keyword, you
+need to use the <I>first</I> option if you want an initial snapshot written
+to the dump file.  The <I>every</I> keyword cannot be used with the dump
+<I>dcd</I> style.
 </P>
 <P>For example, the following commands will
 write snapshots at timesteps 0,10,20,30,100,200,300,1000,2000,etc:
@@ -383,6 +389,25 @@ write snapshots at timesteps 0,10,20,30,100,200,300,1000,2000,etc:
 dump		1 all atom 100 tmp.dump
 dump_modify	1 every v_s first yes 
 </PRE>
+<P>The following commands would write snapshots at the timesteps listed
+in file tmp.times:
+</P>
+<PRE>variable        f file tmp.times
+variable	s equal next(f)
+dump		1 all atom 100 tmp.dump
+dump_modify	1 every v_s 
+</PRE>
+<P>IMPORTANT NOTE: When using a file-style variable with the <I>every</I>
+keyword, the file of timesteps must list a first timestep that is
+beyond the current timestep (e.g. it cannot be 0).  And it must list
+one or more timesteps beyond the length of the run you perform.  This
+is because the dump command will generate an error if the next
+timestep it reads from the file is not a value greater than the
+current timestep.  Thus if you wanted output on steps 0,15,100 of a
+100-timestep run, the file should contain the values 15,100,101 and
+you should also use the dump_modify first command.  Any final value >
+100 could be used in place of 101.
+</P>
 <HR>
 
 <P>The <I>first</I> keyword determines whether a dump snapshot is written on

doc/dump_modify.txt

@@ -358,16 +358,22 @@ which should be specified as v_name, where name is the variable name.
 
 In this case, the variable is evaluated at the beginning of a run to
 determine the next timestep at which a dump snapshot will be written
-out.  On that timestep, the variable will be evaluated again to
+out.  On that timestep the variable will be evaluated again to
 determine the next timestep, etc.  Thus the variable should return
 timestep values.  See the stagger() and logfreq() and stride() math
 functions for "equal-style variables"_variable.html, as examples of
 useful functions to use in this context.  Other similar math functions
 could easily be added as options for "equal-style
-variables"_variable.html.  When using the variable option with the
+variables"_variable.html.  Also see the next() function, which allows
-{every} keyword, you also need to use the {first} option if you want
+use of a file-style variable which reads successive values from a
-an initial snapshot written to the dump file.  The {every} keyword
+file, each time the variable is evaluated.  Used with the {every}
-cannot be used with the dump {dcd} style.
+keyword, if the file contains a list of ascending timesteps, you can
+output snapshots whenever you wish.
+
+Note that when using the variable option with the {every} keyword, you
+need to use the {first} option if you want an initial snapshot written
+to the dump file.  The {every} keyword cannot be used with the dump
+{dcd} style.
 
 For example, the following commands will
 write snapshots at timesteps 0,10,20,30,100,200,300,1000,2000,etc:
@@ -376,6 +382,25 @@ variable	s equal logfreq(10,3,10)
 dump		1 all atom 100 tmp.dump
 dump_modify	1 every v_s first yes :pre
 
+The following commands would write snapshots at the timesteps listed
+in file tmp.times:
+
+variable        f file tmp.times
+variable	s equal next(f)
+dump		1 all atom 100 tmp.dump
+dump_modify	1 every v_s :pre
+
+IMPORTANT NOTE: When using a file-style variable with the {every}
+keyword, the file of timesteps must list a first timestep that is
+beyond the current timestep (e.g. it cannot be 0).  And it must list
+one or more timesteps beyond the length of the run you perform.  This
+is because the dump command will generate an error if the next
+timestep it reads from the file is not a value greater than the
+current timestep.  Thus if you wanted output on steps 0,15,100 of a
+100-timestep run, the file should contain the values 15,100,101 and
+you should also use the dump_modify first command.  Any final value >
+100 could be used in place of 101.
+
 :line
 
 The {first} keyword determines whether a dump snapshot is written on

doc/fix_gcmc.html

@@ -25,7 +25,7 @@
 
 <LI>M = number of MC displacements to attempt every N steps 
 
-<LI>type = atom type of exchanged particles 
+<LI>type = atom type or molecule ID of exchanged gas 
 
 <LI>seed = random # seed (positive integer) 
 
@@ -42,49 +42,61 @@
 <PRE>  <I>molecule</I> value = <I>no</I> or <I>yes</I>
   <I>region</I> value = region-ID
     region-ID = ID of region to use as an exchange/move volume 
+  <I>maxangle</I> value = maximum molecular rotation angle (degrees) 
 </PRE>
 
 </UL>
 <P><B>Examples:</B>
 </P>
-<PRE>fix 2 all gcmc 10 1000 1000 2 29494 298.0 -0.5 0.01
+<PRE>fix 2 gas gcmc 10 1000 1000 2 29494 298.0 -0.5 0.01
-fix 3 all gcmc 10 100 100 1 3456543 3.0 -2.5 0.1 molecule yes
+fix 3 Kr gcmc 10 100 100 1 3456543 3.0 -2.5 0.1 molecule yes maxrot 180
-fix 4 all gcmc 1 10 10 1 123456543 300.0 -12.5 1.0 region disk 
+fix 4 my_gas gcmc 1 10 10 1 123456543 300.0 -12.5 1.0 region disk 
 </PRE>
 <P><B>Description:</B>
 </P>
 <P>This fix performs grand canonical Monte Carlo (GCMC) exchanges of
-particles of the given type with an imaginary ideal gas reservoir at
+atoms or molecules of the given type with an imaginary ideal gas reservoir at
-the specified T and chemical potential (mu    ) as discussed in
+the specified T and chemical potential (mu) as discussed in
 <A HREF = "#Frenkel">(Frenkel)</A>. If used with the <A HREF = "fix_nh.html">fix nvt</A> command,
 simulations in the grand canonical enemble (muVT, constant chemical
 potential, constant volume, and constant temperature) can be
 performed.  Specific uses include computing isotherms in microporous
 materials, or computing vapor-liquid coexistence curves.
 </P>
-<P>Perform up to X exchanges of particles of the given type between the
+<P>Perform up to X exchanges of gas atoms or molecules of the given type 
-simulation domain and the imaginary reservoir every N timesteps. Also
+between the simulation domain and the imaginary reservoir every N 
-perform M Monte Carlo displacements of particles of the given type
+timesteps. Also perform M Monte Carlo displacements or rotations 
-within the simulation domain. M should typically be chosen to be
+(for molecules) of gas of the given type within the simulation domain.
-approximately equal to the expected number of particles of the given
+M should typically be chosen to be approximately equal to the expected
-type within the domain, which will result in roughly one MC
+number of gas atoms or molecules of the given type within the domain, 
-translation per particle per MC cycle.
+which will result in roughly one MC translation per atom or molecule 
+per MC cycle.
 </P>
-<P>This fix cannot be used to perform MC displacements of particles other
+<P>For MC moves of molecular gasses, rotations and translations are each 
-than the exchanged type. All particles in the simulation domain can be
+attempted with 50% probability. For MC moves of atomic gasses, 
-moved using regular time integration displacements, e.g. via
+translations are attempted 100% of the time. For MC exchanges of either
-<A HREF = "fix_nvt.html">fix_nvt</A>, resulting in a hybrid GCMC+MD simulation.
+molecular or atomic gasses, deletions and insertions are each attempted
+with 50% probability.
+</P>
+<P>This fix cannot be used to perform MC insertions of gas atoms or
+molecules other than the exchanged type, but MC deletions, translations,
+and rotations can be performed on any atom/molecule in the fix group. 
+All atoms in the simulation domain can be moved using regular time 
+integration displacements, e.g. via <A HREF = "fix_nvt.html">fix_nvt</A>, resulting
+in a hybrid GCMC+MD simulation. A smaller-than-usual timestep size 
+may be needed when running such a hybrid simulation, especially if 
+the inserted molecules are not well equilibrated.
 </P>
 <P>This command may optionally use the <I>region</I> keyword to define an 
 exchange and move volume.  The specified region must have been 
 previously defined with a <A HREF = "region.html">region</A> command.  It must be 
 defined with side = <I>in</I>.  Insertion attempts occur only within the 
 specified region. Move and deletion attempt candidates are selected 
-from particles within the region. If no candidate can be found within
+from gas atoms or molecules within the region. If no candidate can be 
-the specified region after randomly selecting candidates 1000 times,
+found within the specified region after randomly selecting candidates 
-the move or deletion attempt is considered a failure. Moves must start
+1000 times, the move or deletion attempt is considered a failure. Moves 
-within the specified region, but may move the particle slightly
+must start within the specified region, but may move the atom or molecule
-outside of the region.
+slightly outside of the region.
 </P>
 <P>If used with <A HREF = "fix_nvt.html">fix_nvt</A>, the temperature of the imaginary
 reservoir, T, should be set to be equivalent to the target temperature
@@ -95,17 +107,16 @@ will not be in thermal equilibrium with the simulation domain.
 invoked, so you should not set N to be too small.  However, periodic
 rebuilds are necessary in order to avoid dangerous rebuilds and missed
 interactions. Specifically, avoid performing so many MC displacements
-per timestep that a particle can move beyond the neighbor list skin
+per timestep that atoms can move beyond the neighbor list skin
 distance. See the <A HREF = "neighbor.html">neighbor</A> command for details.
 </P>
-<P>When a particle is to be inserted, its coordinates are chosen as a
+<P>When an atom or molecule is to be inserted, its center-of-mass 
-random position within the current simulation domain, and its velocity
+coordinates are chosen as a random position within the current 
-is randomly chosen from the specified temperature distribution given
+simulation domain, and new atom velocities are randomly chosen 
-by T.
+from the specified temperature distribution given by T.  Relative 
-</P>
+coordinates for atoms in a molecule are taken from the template
-<P>Exchanged particles have the specified atom type and are assigned to
+molecule provided by the user. A random initial rotation is used in the
-two groups: the default group "all" and the group specified in the fix
+case of molecule insertions.
-gcmc command (which can also be "all").
 </P>
 <P>If the setting for the <I>molecule</I> keyword is <I>no</I>, then only single
 atoms are exchanged.  In this case, you should ensure you do not
@@ -116,7 +127,28 @@ non-zero molecule ID, but does not check for this at the time of
 deletion.
 </P>
 <P>If the setting for the <I>molecule</I> keyword is <I>yes</I>, entire molecules
-are exchanged. This feature is not yet supported.
+are exchanged. The user must supply a model molecule in the data
+file to use as a template for exchanges, and that molecule's number
+must be given in the fix GCMC command as the "type" of the exchanged
+gas. The model molecule can then be immediately deleted using the
+<A HREF = "delete_atoms.html">delete_atoms</A> command.
+</P>
+<P>Optionally, users may specify the maximum rotation angle for 
+molecular rotations using the <I>maxangle</I> keyword and specifying
+the angle in degrees. The specified angle will apply to all three 
+Euler angles used internally to define the rotation matrix for
+molecular rotations. The max angle can be set to zero, but rotations
+will be pointless. Note that the default is ten degrees for each 
+Euler angle.
+</P>
+<P>For atomic gasses, inserted atoms have the specified atom type, but
+deleted atoms are any atoms that have been inserted or that belong 
+to the user-specified fix group. For molecular gasses, exchanged 
+molecules use the same atom types as in the template molecule 
+supplied by the user.  In both cases, exchanged
+atoms/molecules are assigned to two groups: the default group "all"
+and the group specified in the fix gcmc command (which can also be 
+"all"). 
 </P>
 <P>Use of this fix typically will cause the number of atoms to fluctuate,
 therefore, you will want to use the
@@ -150,10 +182,12 @@ values are the following global cummulative quantities:
 </P>
 <UL><LI>1 = displacement attempts
 <LI>2 = displacement successes
-<LI>3 = deletion attempts
+<LI>3 = insertion attempts
-<LI>4 = deletion successes
+<LI>4 = insertion successes 
-<LI>5 = insertion attempts
+<LI>5 = deletion attempts
-<LI>6 = insertion successes 
+<LI>6 = deletion successes
+<LI>7 = rotation attempts
+<LI>8 = rotation successes 
 </UL>
 <P>The vector values calculated by this fix are "extensive".
 </P>
@@ -170,9 +204,6 @@ LAMMPS</A> section for more info.
 <P>Do not set "neigh_modify once yes" or else this fix will never be
 called.  Reneighboring is required.
 </P>
-<P>You cannot currently exchange charged particles or molecules with a
-net charge.
-</P>
 <P>Only pairwise interactions, as defined by the
 <A HREF = "pair_style.html">pair_style</A> command, are included in this
 calculation.  Long-range interactions due to a
@@ -183,14 +214,25 @@ this fix.  For example, 3-body potentials, such as
 be used.  <A HREF = "pair_eam.html">EAM</A> potentials for metals only include the
 pair potential portion of the EAM interaction, not the embedding term.
 </P>
+<P>Can be run in parallel, but aspects of the GCMC part will not scale 
+well in parallel. Only usable for 3D simulations with orthogonal 
+simulation cells.
+</P>
+<P>Note that very lengthy simulations involving insertions/deletions of
+billions of gas molecules may run out of atom or molecule IDs and
+trigger an error, so it is better to run multiple shorter-duration 
+simulations. Likewise, very large molecules have not been tested
+and may turn out to be problematic.
+</P>
 <P><B>Related commands:</B>
 </P>
 <P><A HREF = "fix_nvt.html">fix_nvt</A>, <A HREF = "neighbor.html">neighbor</A>, 
-<A HREF = "fix_deposit.html">fix_deposit</A>, <A HREF = "fix_evaporate.html">fix_evaporate</A>
+<A HREF = "fix_deposit.html">fix_deposit</A>, <A HREF = "fix_evaporate.html">fix_evaporate</A>,
+<A HREF = "delete_atoms.html">delete_atoms</A>
 </P>
 <P><B>Default:</B>
 </P>
-<P>The option defaults are molecule = no.
+<P>The option defaults are molecule = no, maxangle = 10.
 </P>
 <HR>
 

doc/fix_gcmc.txt

@@ -17,7 +17,7 @@ gcmc = style name of this fix command :l
 N = invoke this fix every N steps :l
 X = number of exchanges to attempt every N steps :l
 M = number of MC displacements to attempt every N steps :l
-type = atom type of exchanged particles :l
+type = atom type or molecule ID of exchanged gas :l
 seed = random # seed (positive integer) :l
 T = temperature of the ideal gas reservoir (temperature units) :l
 mu = chemical potential of the ideal gas reservoir (energy units) :l
@@ -26,49 +26,61 @@ zero or more keyword/value pairs may be appended to args :l
 keyword = {molecule} or {region} :l
   {molecule} value = {no} or {yes}
   {region} value = region-ID
-    region-ID = ID of region to use as an exchange/move volume :pre
+    region-ID = ID of region to use as an exchange/move volume 
+  {maxangle} value = maximum molecular rotation angle (degrees) :pre
 :ule
 
 [Examples:]
 
-fix 2 all gcmc 10 1000 1000 2 29494 298.0 -0.5 0.01
+fix 2 gas gcmc 10 1000 1000 2 29494 298.0 -0.5 0.01
-fix 3 all gcmc 10 100 100 1 3456543 3.0 -2.5 0.1 molecule yes
+fix 3 Kr gcmc 10 100 100 1 3456543 3.0 -2.5 0.1 molecule yes maxrot 180
-fix 4 all gcmc 1 10 10 1 123456543 300.0 -12.5 1.0 region disk :pre
+fix 4 my_gas gcmc 1 10 10 1 123456543 300.0 -12.5 1.0 region disk :pre
 
 [Description:]
 
 This fix performs grand canonical Monte Carlo (GCMC) exchanges of
-particles of the given type with an imaginary ideal gas reservoir at
+atoms or molecules of the given type with an imaginary ideal gas reservoir at
-the specified T and chemical potential (mu    ) as discussed in
+the specified T and chemical potential (mu) as discussed in
 "(Frenkel)"_#Frenkel. If used with the "fix nvt"_fix_nh.html command,
 simulations in the grand canonical enemble (muVT, constant chemical
 potential, constant volume, and constant temperature) can be
 performed.  Specific uses include computing isotherms in microporous
 materials, or computing vapor-liquid coexistence curves.
 
-Perform up to X exchanges of particles of the given type between the
+Perform up to X exchanges of gas atoms or molecules of the given type 
-simulation domain and the imaginary reservoir every N timesteps. Also
+between the simulation domain and the imaginary reservoir every N 
-perform M Monte Carlo displacements of particles of the given type
+timesteps. Also perform M Monte Carlo displacements or rotations 
-within the simulation domain. M should typically be chosen to be
+(for molecules) of gas of the given type within the simulation domain.
-approximately equal to the expected number of particles of the given
+M should typically be chosen to be approximately equal to the expected
-type within the domain, which will result in roughly one MC
+number of gas atoms or molecules of the given type within the domain, 
-translation per particle per MC cycle.
+which will result in roughly one MC translation per atom or molecule 
+per MC cycle.
 
-This fix cannot be used to perform MC displacements of particles other
+For MC moves of molecular gasses, rotations and translations are each 
-than the exchanged type. All particles in the simulation domain can be
+attempted with 50% probability. For MC moves of atomic gasses, 
-moved using regular time integration displacements, e.g. via
+translations are attempted 100% of the time. For MC exchanges of either
-"fix_nvt"_fix_nvt.html, resulting in a hybrid GCMC+MD simulation.
+molecular or atomic gasses, deletions and insertions are each attempted
+with 50% probability.
+
+This fix cannot be used to perform MC insertions of gas atoms or
+molecules other than the exchanged type, but MC deletions, translations,
+and rotations can be performed on any atom/molecule in the fix group. 
+All atoms in the simulation domain can be moved using regular time 
+integration displacements, e.g. via "fix_nvt"_fix_nvt.html, resulting
+in a hybrid GCMC+MD simulation. A smaller-than-usual timestep size 
+may be needed when running such a hybrid simulation, especially if 
+the inserted molecules are not well equilibrated.
 
 This command may optionally use the {region} keyword to define an 
 exchange and move volume.  The specified region must have been 
 previously defined with a "region"_region.html command.  It must be 
 defined with side = {in}.  Insertion attempts occur only within the 
 specified region. Move and deletion attempt candidates are selected 
-from particles within the region. If no candidate can be found within
+from gas atoms or molecules within the region. If no candidate can be 
-the specified region after randomly selecting candidates 1000 times,
+found within the specified region after randomly selecting candidates 
-the move or deletion attempt is considered a failure. Moves must start
+1000 times, the move or deletion attempt is considered a failure. Moves 
-within the specified region, but may move the particle slightly
+must start within the specified region, but may move the atom or molecule
-outside of the region.
+slightly outside of the region.
 
 If used with "fix_nvt"_fix_nvt.html, the temperature of the imaginary
 reservoir, T, should be set to be equivalent to the target temperature
@@ -79,17 +91,16 @@ Note that neighbor lists are re-built every timestep that this fix is
 invoked, so you should not set N to be too small.  However, periodic
 rebuilds are necessary in order to avoid dangerous rebuilds and missed
 interactions. Specifically, avoid performing so many MC displacements
-per timestep that a particle can move beyond the neighbor list skin
+per timestep that atoms can move beyond the neighbor list skin
 distance. See the "neighbor"_neighbor.html command for details.
 
-When a particle is to be inserted, its coordinates are chosen as a
+When an atom or molecule is to be inserted, its center-of-mass 
-random position within the current simulation domain, and its velocity
+coordinates are chosen as a random position within the current 
-is randomly chosen from the specified temperature distribution given
+simulation domain, and new atom velocities are randomly chosen 
-by T.
+from the specified temperature distribution given by T.  Relative 
-
+coordinates for atoms in a molecule are taken from the template
-Exchanged particles have the specified atom type and are assigned to
+molecule provided by the user. A random initial rotation is used in the
-two groups: the default group "all" and the group specified in the fix
+case of molecule insertions.
-gcmc command (which can also be "all").
 
 If the setting for the {molecule} keyword is {no}, then only single
 atoms are exchanged.  In this case, you should ensure you do not
@@ -100,7 +111,28 @@ non-zero molecule ID, but does not check for this at the time of
 deletion.
 
 If the setting for the {molecule} keyword is {yes}, entire molecules
-are exchanged. This feature is not yet supported.
+are exchanged. The user must supply a model molecule in the data
+file to use as a template for exchanges, and that molecule's number
+must be given in the fix GCMC command as the "type" of the exchanged
+gas. The model molecule can then be immediately deleted using the
+"delete_atoms"_delete_atoms.html command.
+
+Optionally, users may specify the maximum rotation angle for 
+molecular rotations using the {maxangle} keyword and specifying
+the angle in degrees. The specified angle will apply to all three 
+Euler angles used internally to define the rotation matrix for
+molecular rotations. The max angle can be set to zero, but rotations
+will be pointless. Note that the default is ten degrees for each 
+Euler angle.
+
+For atomic gasses, inserted atoms have the specified atom type, but
+deleted atoms are any atoms that have been inserted or that belong 
+to the user-specified fix group. For molecular gasses, exchanged 
+molecules use the same atom types as in the template molecule 
+supplied by the user.  In both cases, exchanged
+atoms/molecules are assigned to two groups: the default group "all"
+and the group specified in the fix gcmc command (which can also be 
+"all"). 
 
 Use of this fix typically will cause the number of atoms to fluctuate,
 therefore, you will want to use the
@@ -134,10 +166,12 @@ values are the following global cummulative quantities:
 
 1 = displacement attempts
 2 = displacement successes
-3 = deletion attempts
+3 = insertion attempts
-4 = deletion successes
+4 = insertion successes 
-5 = insertion attempts
+5 = deletion attempts
-6 = insertion successes :ul
+6 = deletion successes
+7 = rotation attempts
+8 = rotation successes :ul
 
 The vector values calculated by this fix are "extensive".
 
@@ -154,9 +188,6 @@ LAMMPS"_Section_start.html#start_3 section for more info.
 Do not set "neigh_modify once yes" or else this fix will never be
 called.  Reneighboring is required.
 
-You cannot currently exchange charged particles or molecules with a
-net charge.
-
 Only pairwise interactions, as defined by the
 "pair_style"_pair_style.html command, are included in this
 calculation.  Long-range interactions due to a
@@ -167,14 +198,25 @@ this fix.  For example, 3-body potentials, such as
 be used.  "EAM"_pair_eam.html potentials for metals only include the
 pair potential portion of the EAM interaction, not the embedding term.
 
+Can be run in parallel, but aspects of the GCMC part will not scale 
+well in parallel. Only usable for 3D simulations with orthogonal 
+simulation cells.
+
+Note that very lengthy simulations involving insertions/deletions of
+billions of gas molecules may run out of atom or molecule IDs and
+trigger an error, so it is better to run multiple shorter-duration 
+simulations. Likewise, very large molecules have not been tested
+and may turn out to be problematic.
+
 [Related commands:]
 
 "fix_nvt"_fix_nvt.html, "neighbor"_neighbor.html, 
-"fix_deposit"_fix_deposit.html, "fix_evaporate"_fix_evaporate.html
+"fix_deposit"_fix_deposit.html, "fix_evaporate"_fix_evaporate.html,
+"delete_atoms"_delete_atoms.html
 
 [Default:]
 
-The option defaults are molecule = no.
+The option defaults are molecule = no, maxangle = 10.
 
 :line
 

doc/fix_thermal_conductivity.html

@@ -91,9 +91,17 @@ number.
 swaps is computed by the fix and can be output.  Dividing this
 quantity by time and the cross-sectional area of the simulation box
 yields a heat flux.  The ratio of heat flux to the slope of the
-temperature profile is the thermal conductivity of the fluid,
+temperature profile is proportional to the thermal conductivity of the
-in appopriate units.  See the <A HREF = "#Muller-Plathe">Muller-Plathe paper</A> for
+fluid, in appropriate units.  See the <A HREF = "#Muller-Plathe">Muller-Plathe
-details.
+paper</A> for details.
+</P>
+<P>IMPORTANT NOTE: If your system is periodic in the direction of the
+heat flux, then the flux is going in 2 directions.  This means the
+effective heat flux in one direction is reduced by a factor of 2.  You
+will see this in the equations for thermal conductivity (kappa) in the
+Muller-Plathe paper.  LAMMPS is simply tallying kinetic energy which
+does not account for whether or not your system is periodic; you must
+use the value appropriately to yield a kappa for your system.
 </P>
 <P>IMPORTANT NOTE: After equilibration, if the temperature gradient you
 observe is not linear, then you are likely swapping energy too

doc/fix_thermal_conductivity.txt

@@ -81,9 +81,17 @@ As described below, the total kinetic energy transferred by these
 swaps is computed by the fix and can be output.  Dividing this
 quantity by time and the cross-sectional area of the simulation box
 yields a heat flux.  The ratio of heat flux to the slope of the
-temperature profile is the thermal conductivity of the fluid,
+temperature profile is proportional to the thermal conductivity of the
-in appopriate units.  See the "Muller-Plathe paper"_#Muller-Plathe for
+fluid, in appropriate units.  See the "Muller-Plathe
-details.
+paper"_#Muller-Plathe for details.
+
+IMPORTANT NOTE: If your system is periodic in the direction of the
+heat flux, then the flux is going in 2 directions.  This means the
+effective heat flux in one direction is reduced by a factor of 2.  You
+will see this in the equations for thermal conductivity (kappa) in the
+Muller-Plathe paper.  LAMMPS is simply tallying kinetic energy which
+does not account for whether or not your system is periodic; you must
+use the value appropriately to yield a kappa for your system.
 
 IMPORTANT NOTE: After equilibration, if the temperature gradient you
 observe is not linear, then you are likely swapping energy too

doc/fix_viscosity.html

@@ -92,9 +92,17 @@ sense.  This is why <I>Nbin</I> is restricted to being an even number.
 swaps is computed by the fix and can be output.  Dividing this
 quantity by time and the cross-sectional area of the simulation box
 yields a momentum flux.  The ratio of momentum flux to the slope of
-the shear velocity profile is the viscosity of the fluid, in
+the shear velocity profile is proportional to the viscosity of the
-appopriate units.  See the <A HREF = "#Muller-Plathe">Muller-Plathe paper</A> for
+fluid, in appropriate units.  See the <A HREF = "#Muller-Plathe">Muller-Plathe
-details.
+paper</A> for details.
+</P>
+<P>IMPORTANT NOTE: If your system is periodic in the direction of the
+momentum flux, then the flux is going in 2 directions.  This means the
+effective momentum flux in one direction is reduced by a factor of 2.
+You will see this in the equations for viscosity in the Muller-Plathe
+paper.  LAMMPS is simply tallying momentum which does not account for
+whether or not your system is periodic; you must use the value
+appropriately to yield a viscosity for your system.
 </P>
 <P>IMPORTANT NOTE: After equilibration, if the velocity profile you
 observe is not linear, then you are likely swapping momentum too

doc/fix_viscosity.txt

@@ -81,9 +81,17 @@ As described below, the total momentum transferred by these velocity
 swaps is computed by the fix and can be output.  Dividing this
 quantity by time and the cross-sectional area of the simulation box
 yields a momentum flux.  The ratio of momentum flux to the slope of
-the shear velocity profile is the viscosity of the fluid, in
+the shear velocity profile is proportional to the viscosity of the
-appopriate units.  See the "Muller-Plathe paper"_#Muller-Plathe for
+fluid, in appropriate units.  See the "Muller-Plathe
-details.
+paper"_#Muller-Plathe for details.
+
+IMPORTANT NOTE: If your system is periodic in the direction of the
+momentum flux, then the flux is going in 2 directions.  This means the
+effective momentum flux in one direction is reduced by a factor of 2.
+You will see this in the equations for viscosity in the Muller-Plathe
+paper.  LAMMPS is simply tallying momentum which does not account for
+whether or not your system is periodic; you must use the value
+appropriately to yield a viscosity for your system.
 
 IMPORTANT NOTE: After equilibration, if the velocity profile you
 observe is not linear, then you are likely swapping momentum too

doc/improper_fourier.html

@@ -59,11 +59,4 @@ section for more info on packages.
 </P>
 <P><B>Default:</B> none
 </P>
-<HR>
-
-<A NAME = "Mayo"></A>
-
-<P><B>(Mayo)</B> Mayo, Olfason, Goddard III, J Phys Chem, 94, 8897-8909
-(1990),
-</P>
 </HTML>

doc/improper_fourier.txt

@@ -56,8 +56,3 @@ section for more info on packages.
 
 [Default:] none
 
-:line
-
-:link(Mayo)
-[(Mayo)] Mayo, Olfason, Goddard III, J Phys Chem, 94, 8897-8909
-(1990),

doc/kspace_style.html

@@ -122,21 +122,22 @@ It should be used with <A HREF = "pair_style.html">pair styles</A> with a
 </P>
 <HR>
 
-<P>The <I>pppm/disp</I> and <I>pppm/disp/tip4p</I> styles add a long-range
+<P>The <I>pppm/disp</I> and <I>pppm/disp/tip4p</I> styles add a mesh-based long-range
-dispersion sum option for 1/r^6 potentials, similar to the
+dispersion sum option for 1/r^6 potentials <A HREF = "#Isele-Holder">(Isele-Holder)</A>,
-<I>ewald/disp</I> style.  The 1/r^6 capability means that Lennard-Jones or
+similar to the <I>ewald/disp</I> style. The 1/r^6 capability means
-Buckingham potentials can be used without a cutoff, i.e. they become
+that Lennard-Jones or Buckingham potentials can be used without a cutoff,
-full long-range potentials.
+i.e. they become full long-range potentials.
 </P>
 <P>For these styles, it is currently recommended that you set the
 dispersion mesh size and other parameters explicitly via the
 <A HREF = "kspace_modify.html">kspace_modify</A> command, rather than let LAMMPS set
 them automatically.  For example, a set of parameters that works well
-for TIP4P water models is a LJ cutoff of 10 Angstrom, interpolation
+for surface systems when using real units is a LJ cutoff of 10 Angstrom,
-order = 5 (the default), grid spacing = 4.17 Angstroms, and Ewald
+interpolation order = 5 (the default), grid spacing = 4.17 Angstroms,
-parameter = 0.28.  These parameters work well for the <I>ik</I>
+and Ewald parameter = 0.28.  These parameters work well for the <I>ik</I>
-differentiation.  For the <I>ad</I> setting, a smaller grid is needed,
+differentiation.  For the <I>ad</I> setting, a smaller grid spacing is needed,
-e.g. 3 Angstroms.
+e.g. 3 Angstroms. Further information on the influence of the parameters
+and how to choose them is described in <A HREF = "#Isele-Holder">(Isele-Holder)</A>.
 </P>
 <HR>
 
@@ -282,6 +283,10 @@ Adam Hilger, NY (1989).
 
 <P><B>(Veld)</B> In 't Veld, Ismail, Grest, J Chem Phys, in press (2007).
 </P>
+<A NAME = "Isele-Holder"></A>
+
+<P><B>(Isele-Holder)</B> Isele-Holder, Mitchell, Ismail, J Chem Phys, 137, 174107 (2012). 
+</P>
 <A NAME = "Hardy"></A>
 
 <P><B>(Hardy)</B> David Hardy thesis: Multilevel Summation for the Fast

doc/kspace_style.txt

@@ -115,21 +115,22 @@ It should be used with "pair styles"_pair_style.html with a
 
 :line
 
-The {pppm/disp} and {pppm/disp/tip4p} styles add a long-range
+The {pppm/disp} and {pppm/disp/tip4p} styles add a mesh-based long-range
-dispersion sum option for 1/r^6 potentials, similar to the
+dispersion sum option for 1/r^6 potentials "(Isele-Holder)"_#Isele-Holder,
-{ewald/disp} style.  The 1/r^6 capability means that Lennard-Jones or
+similar to the {ewald/disp} style. The 1/r^6 capability means
-Buckingham potentials can be used without a cutoff, i.e. they become
+that Lennard-Jones or Buckingham potentials can be used without a cutoff,
-full long-range potentials.
+i.e. they become full long-range potentials.
 
 For these styles, it is currently recommended that you set the
 dispersion mesh size and other parameters explicitly via the
 "kspace_modify"_kspace_modify.html command, rather than let LAMMPS set
 them automatically.  For example, a set of parameters that works well
-for TIP4P water models is a LJ cutoff of 10 Angstrom, interpolation
+for surface systems when using real units is a LJ cutoff of 10 Angstrom,
-order = 5 (the default), grid spacing = 4.17 Angstroms, and Ewald
+interpolation order = 5 (the default), grid spacing = 4.17 Angstroms,
-parameter = 0.28.  These parameters work well for the {ik}
+and Ewald parameter = 0.28.  These parameters work well for the {ik}
-differentiation.  For the {ad} setting, a smaller grid is needed,
+differentiation.  For the {ad} setting, a smaller grid spacing is needed,
-e.g. 3 Angstroms.
+e.g. 3 Angstroms. Further information on the influence of the parameters
+and how to choose them is described in "(Isele-Holder)"_#Isele-Holder.
 
 :line
 
@@ -268,6 +269,9 @@ Adam Hilger, NY (1989).
 :link(Veld)
 [(Veld)] In 't Veld, Ismail, Grest, J Chem Phys, in press (2007).
 
+:link(Isele-Holder)
+[(Isele-Holder)] Isele-Holder, Mitchell, Ismail, J Chem Phys, 137, 174107 (2012). 
+
 :link(Hardy)
 [(Hardy)] David Hardy thesis: Multilevel Summation for the Fast
 Evaluation of Forces for the Simulation of Biomolecules, University of

doc/next.html

@@ -38,14 +38,15 @@ be used in an input script command as $a or $z.  If it is multiple
 letters, it can be used as ${myTemp}.
 </P>
 <P>If multiple variables are used as arguments to the <I>next</I> command,
-then all must be of the same variable style: <I>index</I>, <I>loop</I>,
+then all must be of the same variable style: <I>index</I>, <I>loop</I>, <I>file</I>,
 <I>universe</I>, or <I>uloop</I>.  An exception is that <I>universe</I>- and
 <I>uloop</I>-style variables can be mixed in the same <I>next</I> command.
 </P>
 <P>All the variables specified with the next command are incremented by
-one value from their respective list or values.  <I>String-</I> or <I>atom</I>-
+one value from their respective list of values.  A <I>file</I>-style variable
-or <I>equal</I>- or <I>world</I>-style variables cannot be used with the the
+reads the next line from its associated file.  <I>String-</I> or <I>atom</I>- or
-next command, since they only store a single value.
+<I>equal</I>- or <I>world</I>-style variables cannot be used with the the next
+command, since they only store a single value.
 </P>
 <P>When any of the variables in the next command has no more values, a
 flag is set that causes the input script to skip the next
@@ -53,18 +54,19 @@ flag is set that causes the input script to skip the next
 a next command to exit.  As explained in the <A HREF = "variable.html">variable</A>
 command, the variable that has exhausted its values is also deleted.
 This allows it to be used and re-defined later in the input script.
+<I>File</I>-style variables are exhausted when the end-of-file is reached.
 </P>
-<P>When the next command is used with <I>index</I>- or <I>loop</I>-style variables,
+<P>When the next command is used with <I>index</I>- or <I>loop</I>-style or
-the next value is assigned to the variable for all processors.  When
+<I>file</I>-style variables, the next value is assigned to the variable for
-the next command is used with <I>universe</I>- or <I>uloop</I>-style variables,
+all processors.  When the next command is used with <I>universe</I>- or
-the next value is assigned to whichever processor partition executes
+<I>uloop</I>-style variables, the next value is assigned to whichever
-the command first.  All processors in the partition are assigned the
+processor partition executes the command first.  All processors in the
-same value.  Running LAMMPS on multiple partitions of processors via
+partition are assigned the same value.  Running LAMMPS on multiple
-the "-partition" command-line switch is described in <A HREF = "Section_start.html#start_7">this
+partitions of processors via the "-partition" command-line switch is
-section</A> of the manual.  <I>Universe</I>- and
+described in <A HREF = "Section_start.html#start_7">this section</A> of the manual.
-<I>uloop</I>-style variables are incremented using the files
+<I>Universe</I>- and <I>uloop</I>-style variables are incremented using the
-"tmp.lammps.variable" and "tmp.lammps.variable.lock" which you will
+files "tmp.lammps.variable" and "tmp.lammps.variable.lock" which you
-see in your directory during such a LAMMPS run.
+will see in your directory during such a LAMMPS run.
 </P>
 <P>Here is an example of running a series of simulations using the next
 command with an <I>index</I>-style variable.  If this input script is named

doc/next.txt

@@ -35,14 +35,15 @@ be used in an input script command as $a or $z.  If it is multiple
 letters, it can be used as $\{myTemp\}.
 
 If multiple variables are used as arguments to the {next} command,
-then all must be of the same variable style: {index}, {loop},
+then all must be of the same variable style: {index}, {loop}, {file},
 {universe}, or {uloop}.  An exception is that {universe}- and
 {uloop}-style variables can be mixed in the same {next} command.
 
 All the variables specified with the next command are incremented by
-one value from their respective list or values.  {String-} or {atom}-
+one value from their respective list of values.  A {file}-style variable
-or {equal}- or {world}-style variables cannot be used with the the
+reads the next line from its associated file.  {String-} or {atom}- or
-next command, since they only store a single value.
+{equal}- or {world}-style variables cannot be used with the the next
+command, since they only store a single value.
 
 When any of the variables in the next command has no more values, a
 flag is set that causes the input script to skip the next
@@ -50,18 +51,19 @@ flag is set that causes the input script to skip the next
 a next command to exit.  As explained in the "variable"_variable.html
 command, the variable that has exhausted its values is also deleted.
 This allows it to be used and re-defined later in the input script.
+{File}-style variables are exhausted when the end-of-file is reached.
 
-When the next command is used with {index}- or {loop}-style variables,
+When the next command is used with {index}- or {loop}-style or
-the next value is assigned to the variable for all processors.  When
+{file}-style variables, the next value is assigned to the variable for
-the next command is used with {universe}- or {uloop}-style variables,
+all processors.  When the next command is used with {universe}- or
-the next value is assigned to whichever processor partition executes
+{uloop}-style variables, the next value is assigned to whichever
-the command first.  All processors in the partition are assigned the
+processor partition executes the command first.  All processors in the
-same value.  Running LAMMPS on multiple partitions of processors via
+partition are assigned the same value.  Running LAMMPS on multiple
-the "-partition" command-line switch is described in "this
+partitions of processors via the "-partition" command-line switch is
-section"_Section_start.html#start_7 of the manual.  {Universe}- and
+described in "this section"_Section_start.html#start_7 of the manual.
-{uloop}-style variables are incremented using the files
+{Universe}- and {uloop}-style variables are incremented using the
-"tmp.lammps.variable" and "tmp.lammps.variable.lock" which you will
+files "tmp.lammps.variable" and "tmp.lammps.variable.lock" which you
-see in your directory during such a LAMMPS run.
+will see in your directory during such a LAMMPS run.
 
 Here is an example of running a series of simulations using the next
 command with an {index}-style variable.  If this input script is named

doc/pair_coeff.html

@@ -153,6 +153,7 @@ section of <A HREF = "Section_commands.html#cmd_5">this page</A>.
 <LI><A HREF = "pair_lubricateU.html">pair_style lubricateU</A> - hydrodynamic lubrication forces for Fast Lubrication Dynamics
 <LI><A HREF = "pair_lubricateU.html">pair_style lubricateU/poly</A> - hydrodynamic lubrication forces for Fast Lubrication Dynamics with polydispersity
 <LI><A HREF = "pair_meam.html">pair_style meam</A> - modified embedded atom method (MEAM)
+<LI><A HREF = "pair_mie.html">pair_style mie/cut</A> - Mie potential
 <LI><A HREF = "pair_morse.html">pair_style morse</A> - Morse potential
 <LI><A HREF = "pair_peri.html">pair_style peri/lps</A> - peridynamic LPS potential
 <LI><A HREF = "pair_peri.html">pair_style peri/pmb</A> - peridynamic PMB potential

doc/pair_coeff.txt

@@ -150,6 +150,7 @@ section of "this page"_Section_commands.html#cmd_5.
 "pair_style lubricateU"_pair_lubricateU.html - hydrodynamic lubrication forces for Fast Lubrication Dynamics
 "pair_style lubricateU/poly"_pair_lubricateU.html - hydrodynamic lubrication forces for Fast Lubrication Dynamics with polydispersity
 "pair_style meam"_pair_meam.html - modified embedded atom method (MEAM)
+"pair_style mie/cut"_pair_mie.html - Mie potential
 "pair_style morse"_pair_morse.html - Morse potential
 "pair_style peri/lps"_pair_peri.html - peridynamic LPS potential
 "pair_style peri/pmb"_pair_peri.html - peridynamic PMB potential

doc/pair_mie.html

@@ -0,0 +1,105 @@
+<HTML>
+<CENTER><A HREF = "http://lammps.sandia.gov">LAMMPS WWW Site</A> - <A HREF = "Manual.html">LAMMPS Documentation</A> - <A HREF = "Section_commands.html#comm">LAMMPS Commands</A> 
+</CENTER>
+
+
+
+
+
+
+<HR>
+
+<H3>pair_style mie/cut command 
+</H3>
+<P><B>Syntax:</B>
+</P>
+<PRE>pair_style mie/cut cutoff 
+</PRE>
+<UL><LI>cutoff = global cutoff for mie/cut interactions (distance units) 
+</UL>
+<P><B>Examples:</B>
+</P>
+<PRE>pair_style mie/cut 10.0
+pair_coeff 1 1 0.72 3.40 23.00 6.66
+pair_coeff 2 2 0.30 3.55 12.65 6.00
+pair_coeff 1 2 0.46 3.32 16.90 6.31 
+</PRE>
+<P><B>Description:</B>
+</P>
+<P>The <I>mie/cut</I> style computes the Mie potential, given by
+</P>
+<CENTER><IMG SRC = "Eqs/pair_mie.jpg">
+</CENTER>
+<P>Rc is the cutoff and C is a function that depends on the repulsive and 
+attractive exponents, given by:
+</P>
+<CENTER><IMG SRC = "Eqs/pair_mie2.jpg">
+</CENTER>
+<P>Note that for 12/6 exponents, C is equal to 4 and the formula is the
+same as the standard Lennard-Jones potential.
+</P>
+<P>The following coefficients must be defined for each pair of atoms
+types via the <A HREF = "pair_coeff.html">pair_coeff</A> command as in the examples
+above, 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>epsilon (energy units)
+<LI>sigma (distance units)
+<LI>gammaR
+<LI>gammaA
+<LI>cutoff (distance units) 
+</UL>
+<P>The last coefficient is optional.  If not specified, the global 
+cutoff specified in the pair_style command is used.
+</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 epsilon and sigma coefficients
+and cutoff distance for all of the mie/cut pair styles can be mixed.
+If not explicity defined, both the repulsive and attractive gamma
+exponents for different atoms will be calculated following the same
+mixing rule defined for distances.  The default mix value is
+<I>geometric</I>. See the "pair_modify" command for details.
+</P>
+<P>This pair style supports the <A HREF = "pair_modify.html">pair_modify</A> shift
+option for the energy of the pair interaction.
+</P>
+<P>This pair style supports the <A HREF = "pair_modify.html">pair_modify</A> tail
+option for adding a long-range tail correction to the energy and
+pressure of the pair interaction.
+</P>
+<P>This pair style writes its information to <A HREF = "restart.html">binary restart
+files</A>, so pair_style and pair_coeff commands do not need
+to be specified in an input script that reads a restart file.
+</P>
+<P>This pair style supports the use of the <I>inner</I>, <I>middle</I>, and <I>outer</I>
+keywords of the <A HREF = "run_style.html">run_style respa</A> command, meaning the
+pairwise forces can be partitioned by distance at different levels of
+the rRESPA hierarchy.  See the <A HREF = "run_style.html">run_style</A> command for
+details.
+</P>
+<HR>
+
+<P><B>Restrictions:</B> none
+</P>
+<P><B>Related commands:</B>
+</P>
+<P><A HREF = "pair_coeff.html">pair_coeff</A>
+</P>
+<P><B>Default:</B> none
+</P>
+<HR>
+
+<A NAME = "Mie"></A>
+
+<P><B>(Mie)</B> G. Mie, Ann Phys, 316, 657 (1903).
+</P>
+<A NAME = "Avendano"></A>
+
+<P><B>(Avendano)</B> C. Avendano, T. Lafitte, A. Galindo, C. S. Adjiman,
+G. Jackson, E. Muller, J Phys Chem B, 115, 11154 (2011).
+</P>
+</HTML>

doc/pair_mie.txt

@@ -0,0 +1,100 @@
+"LAMMPS WWW Site"_lws - "LAMMPS Documentation"_ld - "LAMMPS Commands"_lc :c
+
+:link(lws,http://lammps.sandia.gov)
+:link(ld,Manual.html)
+:link(lc,Section_commands.html#comm)
+
+:line
+
+pair_style mie/cut command :h3
+
+[Syntax:]
+
+pair_style mie/cut cutoff :pre
+
+cutoff = global cutoff for mie/cut interactions (distance units) :ul
+
+[Examples:]
+
+pair_style mie/cut 10.0
+pair_coeff 1 1 0.72 3.40 23.00 6.66
+pair_coeff 2 2 0.30 3.55 12.65 6.00
+pair_coeff 1 2 0.46 3.32 16.90 6.31 :pre
+
+[Description:]
+
+The {mie/cut} style computes the Mie potential, given by
+
+:c,image(Eqs/pair_mie.jpg)
+
+Rc is the cutoff and C is a function that depends on the repulsive and 
+attractive exponents, given by:
+
+:c,image(Eqs/pair_mie2.jpg)
+
+Note that for 12/6 exponents, C is equal to 4 and the formula is the
+same as the standard Lennard-Jones potential.
+
+The following coefficients must be defined for each pair of atoms
+types via the "pair_coeff"_pair_coeff.html command as in the examples
+above, 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:
+
+epsilon (energy units)
+sigma (distance units)
+gammaR
+gammaA
+cutoff (distance units) :ul
+
+The last coefficient is optional.  If not specified, the global 
+cutoff specified in the pair_style command is used.
+
+:line
+
+[Mixing, shift, table, tail correction, restart, rRESPA info]:
+
+For atom type pairs I,J and I != J, the epsilon and sigma coefficients
+and cutoff distance for all of the mie/cut pair styles can be mixed.
+If not explicity defined, both the repulsive and attractive gamma
+exponents for different atoms will be calculated following the same
+mixing rule defined for distances.  The default mix value is
+{geometric}. See the "pair_modify" command for details.
+
+This pair style supports the "pair_modify"_pair_modify.html shift
+option for the energy of the pair interaction.
+
+This pair style supports the "pair_modify"_pair_modify.html tail
+option for adding a long-range tail correction to the energy and
+pressure of the pair interaction.
+
+This pair style writes its information to "binary restart
+files"_restart.html, so pair_style and pair_coeff commands do not need
+to be specified in an input script that reads a restart file.
+
+This pair style supports the use of the {inner}, {middle}, and {outer}
+keywords of the "run_style respa"_run_style.html command, meaning the
+pairwise forces can be partitioned by distance at different levels of
+the rRESPA hierarchy.  See the "run_style"_run_style.html command for
+details.
+
+:line
+
+[Restrictions:] none
+
+[Related commands:]
+
+"pair_coeff"_pair_coeff.html
+
+[Default:] none
+
+:line
+
+:link(Mie)
+[(Mie)] G. Mie, Ann Phys, 316, 657 (1903).
+
+:link(Avendano)
+[(Avendano)] C. Avendano, T. Lafitte, A. Galindo, C. S. Adjiman,
+G. Jackson, E. Muller, J Phys Chem B, 115, 11154 (2011).
+
+

doc/pair_reax_c.html

@@ -302,7 +302,8 @@ appropriate units if your simulation doesn't use "real" units.
 </P>
 <P><B>Default:</B>
 </P>
-<P>The keyword defaults are checkqeq = yes, lgvdw = no, safezone = 1.2, mincap = 50.
+<P>The keyword defaults are checkqeq = yes, lgvdw = no, safezone = 1.2,
+mincap = 50.
 </P>
 <HR>
 
@@ -313,8 +314,8 @@ Journal of Physical Chemistry A, 112, 1040-1053 (2008).
 </P>
 <A NAME = "Aktulga"></A>
 
-<P><B>(Aktulga)</B> Aktulga, Fogarty, Pandit, Grama, Parallel Computing, to
+<P><B>(Aktulga)</B> Aktulga, Fogarty, Pandit, Grama, Parallel Computing, 38,
-appear (2011).
+245-259 (2012).
 </P>
 <A NAME = "Liu_2011"></A>
 

doc/pair_reax_c.txt

@@ -296,7 +296,8 @@ appropriate units if your simulation doesn't use "real" units.
 
 [Default:]
 
-The keyword defaults are checkqeq = yes, lgvdw = no, safezone = 1.2, mincap = 50.
+The keyword defaults are checkqeq = yes, lgvdw = no, safezone = 1.2,
+mincap = 50.
 
 :line
 
@@ -305,8 +306,8 @@ The keyword defaults are checkqeq = yes, lgvdw = no, safezone = 1.2, mincap = 50
 Journal of Physical Chemistry A, 112, 1040-1053 (2008).
 
 :link(Aktulga)
-[(Aktulga)] Aktulga, Fogarty, Pandit, Grama, Parallel Computing, to
+[(Aktulga)] Aktulga, Fogarty, Pandit, Grama, Parallel Computing, 38,
-appear (2011).
+245-259 (2012).
 
 :link(Liu_2011)
 [(Liu)] L. Liu, Y. Liu, S. V. Zybin, H. Sun and W. A. Goddard, Journal

doc/pair_style.html

@@ -155,6 +155,7 @@ section of <A HREF = "Section_commands.html#cmd_5">this page</A>.
 <LI><A HREF = "pair_lubricateU.html">pair_style lubricateU</A> - hydrodynamic lubrication forces for Fast Lubrication Dynamics
 <LI><A HREF = "pair_lubricateU.html">pair_style lubricateU/poly</A> - hydrodynamic lubrication forces for Fast Lubrication with polydispersity
 <LI><A HREF = "pair_meam.html">pair_style meam</A> - modified embedded atom method (MEAM)
+<LI><A HREF = "pair_mie.html">pair_style mie/cut</A> - Mie potential
 <LI><A HREF = "pair_morse.html">pair_style morse</A> - Morse potential
 <LI><A HREF = "pair_peri.html">pair_style peri/lps</A> - peridynamic LPS potential
 <LI><A HREF = "pair_peri.html">pair_style peri/pmb</A> - peridynamic PMB potential

doc/pair_style.txt

@@ -152,6 +152,7 @@ section of "this page"_Section_commands.html#cmd_5.
 "pair_style lubricateU"_pair_lubricateU.html - hydrodynamic lubrication forces for Fast Lubrication Dynamics
 "pair_style lubricateU/poly"_pair_lubricateU.html - hydrodynamic lubrication forces for Fast Lubrication with polydispersity
 "pair_style meam"_pair_meam.html - modified embedded atom method (MEAM)
+"pair_style mie/cut"_pair_mie.html - Mie potential
 "pair_style morse"_pair_morse.html - Morse potential
 "pair_style peri/lps"_pair_peri.html - peridynamic LPS potential
 "pair_style peri/pmb"_pair_peri.html - peridynamic PMB potential

doc/variable.html

@@ -17,7 +17,7 @@
 </PRE>
 <UL><LI>name = name of variable to define 
 
-<LI>style = <I>delete</I> or <I>index</I> or <I>loop</I> or <I>world</I> or <I>universe</I> or <I>uloop</I> or <I>string</I> or <I>equal</I> or <I>atom</I> 
+<LI>style = <I>delete</I> or <I>index</I> or <I>loop</I> or <I>world</I> or <I>universe</I> or <I>uloop</I> or <I>string</I> or <I>file</I> or <I>equal</I> or <I>atom</I> 
 
 <PRE>  <I>delete</I> = no args
   <I>index</I> args = one or more strings
@@ -39,6 +39,7 @@
     N = integer size of loop
     pad = all values will be same length, e.g. 001, 002, ..., 100
   <I>string</I> arg = one string
+  <I>file</I> arg = filename
   <I>equal</I> or <I>atom</I> args = one formula containing numbers, thermo keywords, math operations, group functions, atom values and vectors, compute/fix/variable references
     numbers = 0.0, 100, -5.4, 2.8e-4, etc
     constants = PI
@@ -59,7 +60,7 @@
 		      bound(group,xmin,region), gyration(group,region), ke(group,reigon),
 		      angmom(group,dim,region), torque(group,dim,region), 
 		      inertia(group,dimdim,region), omega(group,dim,region)
-    special functions = sum(x), min(x), max(x), ave(x), trap(x), gmask(x), rmask(x), grmask(x,y)
+    special functions = sum(x), min(x), max(x), ave(x), trap(x), gmask(x), rmask(x), grmask(x,y), next(x)
     atom value = mass[i], type[i], x[i], y[i], z[i], vx[i], vy[i], vz[i], fx[i], fy[i], fz[i]
     atom vector = mass, type, x, y, z, vx, vy, vz, fx, fy, fz
     compute references = c_ID, c_ID[i], c_ID[i][j]
@@ -79,6 +80,7 @@ variable b equal xcm(mol1,x)/2.0
 variable b equal c_myTemp
 variable b atom x*y/vol
 variable foo string myfile
+variable f file values.txt
 variable temp world 300.0 310.0 320.0 ${Tfinal}
 variable x universe 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
 variable x uloop 15 pad
@@ -154,13 +156,13 @@ the variable's string.  The variable name can be referenced as $x if
 the name "x" is a single character, or as ${LoopVar} if the name
 "LoopVar" is one or more characters.
 </P>
-<P>As described below, for variable styles <I>index</I>, <I>loop</I>, <I>universe</I>,
+<P>As described below, for variable styles <I>index</I>, <I>loop</I>, <I>file</I>,
-and <I>uloop</I>, which string is assigned to a variable can be incremented
+<I>universe</I>, and <I>uloop</I>, which string is assigned to a variable can be
-via the <A HREF = "next.html">next</A> command.  When there are no more strings to
+incremented via the <A HREF = "next.html">next</A> command.  When there are no more
-assign, the variable is exhausted and a flag is set that causes the
+strings to assign, the variable is exhausted and a flag is set that
-next <A HREF = "jump.html">jump</A> command encountered in the input script to be
+causes the next <A HREF = "jump.html">jump</A> command encountered in the input
-skipped.  This enables the construction of simple loops in the input
+script to be skipped.  This enables the construction of simple loops
-script that are iterated over and then exited from.
+in the input script that are iterated over and then exited from.
 </P>
 <P>As explained above, an exhausted variable can be re-used in an input
 script.  The <I>delete</I> style also removes the variable, the same as if
@@ -235,6 +237,32 @@ strings are the integers from 1 to N.  This allows generation of long
 list of runs (e.g. 1000) without having to list N strings in the input
 script.
 </P>
+<P>For the <I>string</I> style, a single string is assigned to the variable.
+The only difference between this and using the <I>index</I> style with a
+single string is that a variable with <I>string</I> style can be redefined.
+E.g. by another command later in the input script, or if the script is
+read again in a loop.
+</P>
+<P>For the <I>file</I> style, a filename is provided which contains a list of
+strings to assign to the variable, one per line.  The strings can be
+numeric values if desired; see the discussion of the next() function
+below for equal-style variables, which will convert the string of a
+file-style variable into a numeric value in a formula.
+</P>
+<P>When a file-style variable is defined, the file is opened and the
+string on the first line is read and stored with the variable.  This
+means the variable can then be evaluated as many times as desired and
+will return that string.  There are two ways to cause the next string
+from the file to be read: use the <A HREF = "next.html">next</A> command or the
+next() function in an equal- or atom-style variable, as discussed
+below.
+</P>
+<P>The rules for formatting the file are as follows.  A comment character
+"#" can be used anywhere on a line; text starting with the comment
+character is stripped.  Blank lines are skipped.  The first "word" of
+a non-blank line, delimited by white space, is the "string" assigned
+to the variable.
+</P>
 <HR>
 
 <P>For the <I>equal</I> and <I>atom</I> styles, a single string is specified which
@@ -279,7 +307,7 @@ references to other variables.
 <TR><TD >Math functions</TD><TD > sqrt(x), exp(x), ln(x), log(x), sin(x), cos(x), tan(x), asin(x), acos(x), atan(x), atan2(y,x), random(x,y,z), normal(x,y,z), ceil(x), floor(x), round(x), ramp(x,y), stagger(x,y), logfreq(x,y,z), stride(x,y,z), vdisplace(x,y), swiggle(x,y,z), cwiggle(x,y,z)</TD></TR>
 <TR><TD >Group functions</TD><TD > count(ID), mass(ID), charge(ID), xcm(ID,dim),                  vcm(ID,dim), fcm(ID,dim), bound(ID,dir), 		 gyration(ID), ke(ID), angmom(ID,dim), torque(ID,dim), 		 inertia(ID,dimdim), omega(ID,dim)</TD></TR>
 <TR><TD >Region functions</TD><TD > count(ID,IDR), mass(ID,IDR), charge(ID,IDR), 		  xcm(ID,dim,IDR), vcm(ID,dim,IDR), fcm(ID,dim,IDR), 		  bound(ID,dir,IDR), gyration(ID,IDR), ke(ID,IDR), 		  angmom(ID,dim,IDR), torque(ID,dim,IDR),                   inertia(ID,dimdim,IDR), omega(ID,dim,IDR)</TD></TR>
-<TR><TD >Special functions</TD><TD > sum(x), min(x), max(x), ave(x), trap(x), gmask(x), rmask(x), grmask(x,y)</TD></TR>
+<TR><TD >Special functions</TD><TD > sum(x), min(x), max(x), ave(x), trap(x), gmask(x), rmask(x), grmask(x,y), next(x)</TD></TR>
 <TR><TD >Atom values</TD><TD > mass[i], type[i], x[i], y[i], z[i],              vx[i], vy[i], vz[i], fx[i], fy[i], fz[i]</TD></TR>
 <TR><TD >Atom vectors</TD><TD > mass, type, x, y, z, vx, vy, vz, fx, fy, fz</TD></TR>
 <TR><TD >Compute references</TD><TD > c_ID, c_ID[i], c_ID[i][j]</TD></TR>
@@ -559,6 +587,20 @@ and the second is a region ID.  It can only be used in atom-style
 variables.  It returns a 1 for atoms that are in both the group and
 region, and a 0 for atoms that are not in both.
 </P>
+<P>The next(x) function takes 1 argument which is a variable ID (not
+"v_foo", just "foo").  It must be for a file-style variable.  Each
+time the next() function is invoked (i.e. each time the equal-style or
+atom-style variable is evaluated), the current string value stored by
+the file-style variable is converted to a numeric value returned by
+the function, and the next value from the file is read and stored.
+Note that if the line previously read from the file was not a numeric
+string, then it will typically evaluate to 0.0, which is likely not
+what you want.  Since the file-style variable reads and stores the
+first line of the file when it is defined in the input script, this is
+the value that will be returned the first time the next() function is
+invoked.  If next() is invoked more times than there are lines in the
+file, the value for the last line is repeatedly returned.
+</P>
 <HR>
 
 <H4>Atom Values and Vectors 
@@ -649,20 +691,26 @@ the doc pages for individual fix commands for details.
 
 <H4>Variable References 
 </H4>
-<P>Variable references access quantities calulated by other variables,
+<P>Variable references access quantities stored or calculated by other
-which will cause those variables to be evaluated.  The name in the
+variables, which will cause those variables to be evaluated.  The name
-reference should be replaced by the name of a variable defined
+in the reference should be replaced by the name of a variable defined
-elsewhere in the input script.  As discussed on this doc page,
+elsewhere in the input script. 
-atom-style variables generate a per-atom vector of values; all other
+</P>
-variable styles generate a global scalar value.  An equal-style
+<P>As discussed on this doc page, equal-style variables generate a global
-variable can only use scalar values, which means another equal-style
+scalar numeric value; atom-style variables generate a per-atom vector
-variable or an element of an atom-style variable.  Atom-style
+of numeric values; all other variables store a string.  The formula
-variables can use the same scalar values.  They can also use other
+for an equal-style variable can use any style of variable except an
-atom-style variables.
+atom-style (unless only a single value from the atom-style variable is
+accessed via a subscript).  If a string-storing variable is used,
+the string is converted to a numeric value.  Note that this
+will typically produce a 0.0 if the string is not a numeric string,
+which is likely not what you want.  The formula for an atom-style
+variable can use any style of variable, including other atom-style
+variables.
 </P>
 <P>Examples of different kinds of variable references are as follows.
 There is no ambiguity as to what a reference means, since variables
-produce only a global scalar or a per-atom vectors, never both.
+produce only a global scalar or a per-atom vector, never both.
 </P>
 <DIV ALIGN=center><TABLE  BORDER=1 >
 <TR><TD >v_name</TD><TD > scalar, or per-atom vector</TD></TR>

doc/variable.txt

@@ -13,7 +13,7 @@ variable command :h3
 variable name style args ... :pre
 
 name = name of variable to define :ulb,l
-style = {delete} or {index} or {loop} or {world} or {universe} or {uloop} or {string} or {equal} or {atom} :l
+style = {delete} or {index} or {loop} or {world} or {universe} or {uloop} or {string} or {file} or {equal} or {atom} :l
   {delete} = no args
   {index} args = one or more strings
   {loop} args = N
@@ -34,6 +34,7 @@ style = {delete} or {index} or {loop} or {world} or {universe} or {uloop} or {st
     N = integer size of loop
     pad = all values will be same length, e.g. 001, 002, ..., 100
   {string} arg = one string
+  {file} arg = filename
   {equal} or {atom} args = one formula containing numbers, thermo keywords, math operations, group functions, atom values and vectors, compute/fix/variable references
     numbers = 0.0, 100, -5.4, 2.8e-4, etc
     constants = PI
@@ -54,7 +55,7 @@ style = {delete} or {index} or {loop} or {world} or {universe} or {uloop} or {st
 		      bound(group,xmin,region), gyration(group,region), ke(group,reigon),
 		      angmom(group,dim,region), torque(group,dim,region), 
 		      inertia(group,dimdim,region), omega(group,dim,region)
-    special functions = sum(x), min(x), max(x), ave(x), trap(x), gmask(x), rmask(x), grmask(x,y)
+    special functions = sum(x), min(x), max(x), ave(x), trap(x), gmask(x), rmask(x), grmask(x,y), next(x)
     atom value = mass\[i\], type\[i\], x\[i\], y\[i\], z\[i\], vx\[i\], vy\[i\], vz\[i\], fx\[i\], fy\[i\], fz\[i\]
     atom vector = mass, type, x, y, z, vx, vy, vz, fx, fy, fz
     compute references = c_ID, c_ID\[i\], c_ID\[i\]\[j\]
@@ -73,6 +74,7 @@ variable b equal xcm(mol1,x)/2.0
 variable b equal c_myTemp
 variable b atom x*y/vol
 variable foo string myfile
+variable f file values.txt
 variable temp world 300.0 310.0 320.0 $\{Tfinal\}
 variable x universe 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
 variable x uloop 15 pad
@@ -148,13 +150,13 @@ the variable's string.  The variable name can be referenced as $x if
 the name "x" is a single character, or as $\{LoopVar\} if the name
 "LoopVar" is one or more characters.
 
-As described below, for variable styles {index}, {loop}, {universe},
+As described below, for variable styles {index}, {loop}, {file},
-and {uloop}, which string is assigned to a variable can be incremented
+{universe}, and {uloop}, which string is assigned to a variable can be
-via the "next"_next.html command.  When there are no more strings to
+incremented via the "next"_next.html command.  When there are no more
-assign, the variable is exhausted and a flag is set that causes the
+strings to assign, the variable is exhausted and a flag is set that
-next "jump"_jump.html command encountered in the input script to be
+causes the next "jump"_jump.html command encountered in the input
-skipped.  This enables the construction of simple loops in the input
+script to be skipped.  This enables the construction of simple loops
-script that are iterated over and then exited from.
+in the input script that are iterated over and then exited from.
 
 As explained above, an exhausted variable can be re-used in an input
 script.  The {delete} style also removes the variable, the same as if
@@ -229,6 +231,32 @@ strings are the integers from 1 to N.  This allows generation of long
 list of runs (e.g. 1000) without having to list N strings in the input
 script.
 
+For the {string} style, a single string is assigned to the variable.
+The only difference between this and using the {index} style with a
+single string is that a variable with {string} style can be redefined.
+E.g. by another command later in the input script, or if the script is
+read again in a loop.
+
+For the {file} style, a filename is provided which contains a list of
+strings to assign to the variable, one per line.  The strings can be
+numeric values if desired; see the discussion of the next() function
+below for equal-style variables, which will convert the string of a
+file-style variable into a numeric value in a formula.
+
+When a file-style variable is defined, the file is opened and the
+string on the first line is read and stored with the variable.  This
+means the variable can then be evaluated as many times as desired and
+will return that string.  There are two ways to cause the next string
+from the file to be read: use the "next"_next.html command or the
+next() function in an equal- or atom-style variable, as discussed
+below.
+
+The rules for formatting the file are as follows.  A comment character
+"#" can be used anywhere on a line; text starting with the comment
+character is stripped.  Blank lines are skipped.  The first "word" of
+a non-blank line, delimited by white space, is the "string" assigned
+to the variable.
+
 :line
 
 For the {equal} and {atom} styles, a single string is specified which
@@ -279,7 +307,7 @@ Region functions: count(ID,IDR), mass(ID,IDR), charge(ID,IDR), \
 		  bound(ID,dir,IDR), gyration(ID,IDR), ke(ID,IDR), \
 		  angmom(ID,dim,IDR), torque(ID,dim,IDR), \
                   inertia(ID,dimdim,IDR), omega(ID,dim,IDR)
-Special functions: sum(x), min(x), max(x), ave(x), trap(x), gmask(x), rmask(x), grmask(x,y)
+Special functions: sum(x), min(x), max(x), ave(x), trap(x), gmask(x), rmask(x), grmask(x,y), next(x)
 Atom values: mass\[i\], type\[i\], x\[i\], y\[i\], z\[i\], \
              vx\[i\], vy\[i\], vz\[i\], fx\[i\], fy\[i\], fz\[i\]
 Atom vectors: mass, type, x, y, z, vx, vy, vz, fx, fy, fz
@@ -559,6 +587,20 @@ and the second is a region ID.  It can only be used in atom-style
 variables.  It returns a 1 for atoms that are in both the group and
 region, and a 0 for atoms that are not in both.
 
+The next(x) function takes 1 argument which is a variable ID (not
+"v_foo", just "foo").  It must be for a file-style variable.  Each
+time the next() function is invoked (i.e. each time the equal-style or
+atom-style variable is evaluated), the current string value stored by
+the file-style variable is converted to a numeric value returned by
+the function, and the next value from the file is read and stored.
+Note that if the line previously read from the file was not a numeric
+string, then it will typically evaluate to 0.0, which is likely not
+what you want.  Since the file-style variable reads and stores the
+first line of the file when it is defined in the input script, this is
+the value that will be returned the first time the next() function is
+invoked.  If next() is invoked more times than there are lines in the
+file, the value for the last line is repeatedly returned.
+
 :line
 
 Atom Values and Vectors :h4
@@ -645,20 +687,26 @@ the doc pages for individual fix commands for details.
 
 Variable References :h4
 
-Variable references access quantities calulated by other variables,
+Variable references access quantities stored or calculated by other
-which will cause those variables to be evaluated.  The name in the
+variables, which will cause those variables to be evaluated.  The name
-reference should be replaced by the name of a variable defined
+in the reference should be replaced by the name of a variable defined
-elsewhere in the input script.  As discussed on this doc page,
+elsewhere in the input script. 
-atom-style variables generate a per-atom vector of values; all other
+
-variable styles generate a global scalar value.  An equal-style
+As discussed on this doc page, equal-style variables generate a global
-variable can only use scalar values, which means another equal-style
+scalar numeric value; atom-style variables generate a per-atom vector
-variable or an element of an atom-style variable.  Atom-style
+of numeric values; all other variables store a string.  The formula
-variables can use the same scalar values.  They can also use other
+for an equal-style variable can use any style of variable except an
-atom-style variables.
+atom-style (unless only a single value from the atom-style variable is
+accessed via a subscript).  If a string-storing variable is used,
+the string is converted to a numeric value.  Note that this
+will typically produce a 0.0 if the string is not a numeric string,
+which is likely not what you want.  The formula for an atom-style
+variable can use any style of variable, including other atom-style
+variables.
 
 Examples of different kinds of variable references are as follows.
 There is no ambiguity as to what a reference means, since variables
-produce only a global scalar or a per-atom vectors, never both.
+produce only a global scalar or a per-atom vector, never both.
 
 v_name: scalar, or per-atom vector
 v_name\[I\]: atom I's value in per-atom vector :tb(s=:)

python/lammps.py

@@ -153,7 +153,7 @@ class lammps:
   def gather_atoms(self,name,type,count):
     natoms = self.lib.lammps_get_natoms(self.lmp)
     if type == 0:
-      data = ((count*natoms)*c_double)()
+      data = ((count*natoms)*c_int)()
       self.lib.lammps_gather_atoms(self.lmp,name,type,count,data)
     elif type == 1:
       data = ((count*natoms)*c_double)()

src/CLASS2/pair_lj_class2_coul_long.cpp

@@ -42,6 +42,7 @@ using namespace MathConst;
 PairLJClass2CoulLong::PairLJClass2CoulLong(LAMMPS *lmp) : Pair(lmp)
 {
   ewaldflag = pppmflag = 1;
+  ftable = NULL;
 }
 
 /* ---------------------------------------------------------------------- */
@@ -68,8 +69,9 @@ PairLJClass2CoulLong::~PairLJClass2CoulLong()
 
 void PairLJClass2CoulLong::compute(int eflag, int vflag)
 {
-  int i,j,ii,jj,inum,jnum,itype,jtype;
+  int i,j,ii,jj,inum,jnum,itable,itype,jtype;
   double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,evdwl,ecoul,fpair;
+  double fraction,table;
   double rsq,r,rinv,r2inv,r3inv,r6inv,forcecoul,forcelj;
   double grij,expm2,prefactor,t,erfc;
   double factor_coul,factor_lj;
@@ -122,14 +124,29 @@ void PairLJClass2CoulLong::compute(int eflag, int vflag)
         r2inv = 1.0/rsq;
 
         if (rsq < cut_coulsq) {
-          r = sqrt(rsq);
+          if (!ncoultablebits || rsq <= tabinnersq) {
-          grij = g_ewald * r;
+            r = sqrt(rsq);
-          expm2 = exp(-grij*grij);
+            grij = g_ewald * r;
-          t = 1.0 / (1.0 + EWALD_P*grij);
+            expm2 = exp(-grij*grij);
-          erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2;
+            t = 1.0 / (1.0 + EWALD_P*grij);
-          prefactor = qqrd2e * qtmp*q[j]/r;
+            erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2;
-          forcecoul = prefactor * (erfc + EWALD_F*grij*expm2);
+            prefactor = qqrd2e * qtmp*q[j]/r;
-          if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor;
+            forcecoul = prefactor * (erfc + EWALD_F*grij*expm2);
+            if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor;
+          } else { 
+            union_int_float_t rsq_lookup; 
+            rsq_lookup.f = rsq; 
+            itable = rsq_lookup.i & ncoulmask; 
+            itable >>= ncoulshiftbits; 
+            fraction = (rsq_lookup.f - rtable[itable]) * drtable[itable]; 
+            table = ftable[itable] + fraction*dftable[itable]; 
+            forcecoul = qtmp*q[j] * table; 
+            if (factor_coul < 1.0) { 
+              table = ctable[itable] + fraction*dctable[itable]; 
+              prefactor = qtmp*q[j] * table; 
+              forcecoul -= (1.0-factor_coul)*prefactor; 
+            }
+          }
         } else forcecoul = 0.0;
 
         if (rsq < cut_ljsq[itype][jtype]) {
@@ -152,7 +169,12 @@ void PairLJClass2CoulLong::compute(int eflag, int vflag)
 
         if (eflag) {
           if (rsq < cut_coulsq) {
-            ecoul = prefactor*erfc;
+            if (!ncoultablebits || rsq <= tabinnersq) 
+              ecoul = prefactor*erfc; 
+            else { 
+              table = etable[itable] + fraction*detable[itable]; 
+              ecoul = qtmp*q[j] * table; 
+            }
             if (factor_coul < 1.0) ecoul -= (1.0-factor_coul)*prefactor;
           } else ecoul = 0.0;
           if (rsq < cut_ljsq[itype][jtype]) {
@@ -270,9 +292,12 @@ void PairLJClass2CoulLong::init_style()
 
   // insure use of KSpace long-range solver, set g_ewald
 
- if (force->kspace == NULL)
+  if (force->kspace == NULL)
     error->all(FLERR,"Pair style requires a KSpace style");
   g_ewald = force->kspace->g_ewald;
+
+  // setup force tables 
+  if (ncoultablebits) init_tables(cut_coul,NULL);
 }
 
 /* ----------------------------------------------------------------------
@@ -400,6 +425,8 @@ void PairLJClass2CoulLong::write_restart_settings(FILE *fp)
   fwrite(&offset_flag,sizeof(int),1,fp);
   fwrite(&mix_flag,sizeof(int),1,fp);
   fwrite(&tail_flag,sizeof(int),1,fp);
+  fwrite(&ncoultablebits,sizeof(int),1,fp);
+  fwrite(&tabinner,sizeof(double),1,fp);
 }
 
 /* ----------------------------------------------------------------------
@@ -414,12 +441,16 @@ void PairLJClass2CoulLong::read_restart_settings(FILE *fp)
     fread(&offset_flag,sizeof(int),1,fp);
     fread(&mix_flag,sizeof(int),1,fp);
     fread(&tail_flag,sizeof(int),1,fp);
+    fread(&ncoultablebits,sizeof(int),1,fp);
+    fread(&tabinner,sizeof(double),1,fp);
   }
   MPI_Bcast(&cut_lj_global,1,MPI_DOUBLE,0,world);
   MPI_Bcast(&cut_coul,1,MPI_DOUBLE,0,world);
   MPI_Bcast(&offset_flag,1,MPI_INT,0,world);
   MPI_Bcast(&mix_flag,1,MPI_INT,0,world);
   MPI_Bcast(&tail_flag,1,MPI_INT,0,world);
+  MPI_Bcast(&ncoultablebits,1,MPI_INT,0,world);
+  MPI_Bcast(&tabinner,1,MPI_DOUBLE,0,world);
 }
 
 /* ---------------------------------------------------------------------- */
@@ -430,18 +461,34 @@ double PairLJClass2CoulLong::single(int i, int j, int itype, int jtype,
                                     double &fforce)
 {
   double r2inv,r,rinv,r3inv,r6inv,grij,expm2,t,erfc,prefactor;
-  double forcecoul,forcelj,phicoul,philj;
+  double fraction,table,forcecoul,forcelj,phicoul,philj;
+  int itable;
 
   r2inv = 1.0/rsq;
   if (rsq < cut_coulsq) {
-    r = sqrt(rsq);
+    if (!ncoultablebits || rsq <= tabinnersq) {
-    grij = g_ewald * r;
+      r = sqrt(rsq);
-    expm2 = exp(-grij*grij);
+      grij = g_ewald * r;
-    t = 1.0 / (1.0 + EWALD_P*grij);
+      expm2 = exp(-grij*grij);
-    erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2;
+      t = 1.0 / (1.0 + EWALD_P*grij);
-    prefactor = force->qqrd2e * atom->q[i]*atom->q[j]/r;
+      erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2;
-    forcecoul = prefactor * (erfc + EWALD_F*grij*expm2);
+      prefactor = force->qqrd2e * atom->q[i]*atom->q[j]/r;
-    if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor;
+      forcecoul = prefactor * (erfc + EWALD_F*grij*expm2);
+      if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor;
+    } else {
+      union_int_float_t rsq_lookup;
+      rsq_lookup.f = rsq;
+      itable = rsq_lookup.i & ncoulmask;
+      itable >>= ncoulshiftbits;
+      fraction = (rsq_lookup.f - rtable[itable]) * drtable[itable];
+      table = ftable[itable] + fraction*dftable[itable];
+      forcecoul = atom->q[i]*atom->q[j] * table;
+      if (factor_coul < 1.0) {
+        table = ctable[itable] + fraction*dctable[itable];
+        prefactor = atom->q[i]*atom->q[j] * table;
+        forcecoul -= (1.0-factor_coul)*prefactor;
+      }
+    }
   } else forcecoul = 0.0;
   if (rsq < cut_ljsq[itype][jtype]) {
     rinv = sqrt(r2inv);
@@ -453,7 +500,12 @@ double PairLJClass2CoulLong::single(int i, int j, int itype, int jtype,
 
   double eng = 0.0;
   if (rsq < cut_coulsq) {
-    phicoul = prefactor*erfc;
+    if (!ncoultablebits || rsq <= tabinnersq)
+      phicoul = prefactor*erfc;
+    else {
+      table = etable[itable] + fraction*detable[itable];
+      phicoul = atom->q[i]*atom->q[j] * table;
+    }
     if (factor_coul < 1.0) phicoul -= (1.0-factor_coul)*prefactor;
     eng += phicoul;
   }

src/GPU/pair_dipole_sf_gpu.cpp

@@ -303,6 +303,7 @@ void PairDipoleSFGPU::cpu_compute(int start, int inum, int eflag, int vflag,
           if (mu[i][3] > 0.0 && q[j] != 0.0) { 
             r3inv = r2inv*rinv;
             r5inv = r3inv*r2inv;
+            pidotr = mu[i][0]*delx + mu[i][1]*dely + mu[i][2]*delz; 
             rcutcoul2inv=1.0/cut_coulsq[itype][jtype];
             pre1 = 3.0 * q[j] * r5inv * pidotr * (1-rsq*rcutcoul2inv);
             pqfac = 1.0 - 3.0*rsq*rcutcoul2inv + 
@@ -320,6 +321,7 @@ void PairDipoleSFGPU::cpu_compute(int start, int inum, int eflag, int vflag,
           if (mu[j][3] > 0.0 && qtmp != 0.0) { 
             r3inv = r2inv*rinv;
             r5inv = r3inv*r2inv;
+            pjdotr = mu[j][0]*delx + mu[j][1]*dely + mu[j][2]*delz;
             rcutcoul2inv=1.0/cut_coulsq[itype][jtype];
             pre1 = 3.0 * qtmp * r5inv * pjdotr * (1-rsq*rcutcoul2inv);
             qpfac = 1.0 - 3.0*rsq*rcutcoul2inv +

src/GPU/pppm_gpu.cpp

@@ -24,6 +24,7 @@
 #include "pppm_gpu.h"
 #include "atom.h"
 #include "comm.h"
+#include "commgrid.h"
 #include "neighbor.h"
 #include "force.h"
 #include "pair.h"
@@ -48,6 +49,9 @@ using namespace MathConst;
 #define LARGE 10000.0
 #define EPS_HOC 1.0e-7
 
+enum{REVERSE_RHO};
+enum{FORWARD_IK,FORWARD_AD,FORWARD_IK_PERATOM,FORWARD_AD_PERATOM};
+
 #ifdef FFT_SINGLE
 #define ZEROF 0.0f
 #define ONEF  1.0f
@@ -56,7 +60,7 @@ using namespace MathConst;
 #define ONEF  1.0
 #endif
 
-// External functions from cuda library for atom decomposition
+// external functions from cuda library for atom decomposition
 
 #ifdef FFT_SINGLE
 #define PPPM_GPU_API(api)  pppm_gpu_ ## api ## _f
@@ -111,10 +115,35 @@ PPPMGPU::~PPPMGPU()
 
 void PPPMGPU::init()
 {
+  // PPPM init manages all arrays except density_brick_gpu and vd_brick
+  //      thru its deallocate(), allocate()
+  // NOTE: could free density_brick and vdxyz_brick after PPPM allocates them,
+  //       before allocating db_gpu and vd_brick down below, if don't need,
+  //       if do this, make sure to set them to NULL
+  // NOTE: delete/realloc of cg necessary b/c packing 4 values per grid pt,
+  //       not 3 as PPPM does - probably a better way to account for this
+  //       in PPPM::init()
+
+  destroy_3d_offset(density_brick_gpu,nzlo_out,nylo_out);
+  destroy_3d_offset(vd_brick,nzlo_out,nylo_out);
+  density_brick_gpu = vd_brick = NULL;
+
   PPPM::init();
 
+  if (differentiation_flag == 0) {
+    delete cg;
+    int (*procneigh)[2] = comm->procneigh;
+    cg = new CommGrid(lmp,world,4,1,
+                      nxlo_in,nxhi_in,nylo_in,nyhi_in,nzlo_in,nzhi_in,
+                      nxlo_out,nxhi_out,nylo_out,nyhi_out,nzlo_out,nzhi_out,
+                      procneigh[0][0],procneigh[0][1],procneigh[1][0],
+                      procneigh[1][1],procneigh[2][0],procneigh[2][1]);
+  }
+
+  // unsupported option
+
   if (differentiation_flag == 1)
-    error->all(FLERR,"Cannot (yet) do analytic differentiation with pppm/gpu.");
+    error->all(FLERR,"Cannot (yet) do analytic differentiation with pppm/gpu");
 
   if (strcmp(update->integrate_style,"verlet/split") == 0) {
     kspace_split=true;
@@ -142,6 +171,8 @@ void PPPMGPU::init()
 
   GPU_EXTRA::check_flag(success,error,world);
 
+  // allocate density_brick_gpu and vd_brick
+
   density_brick_gpu =
     create_3d_offset(nzlo_out,nzhi_out,nylo_out,nyhi_out,
                      nxlo_out,nxhi_out,"pppm:density_brick_gpu",h_brick,1);
@@ -149,7 +180,7 @@ void PPPMGPU::init()
     create_3d_offset(nzlo_out,nzhi_out,nylo_out,nyhi_out,
                      nxlo_out,nxhi_out,"pppm:vd_brick",data,4);
 
-  poisson_time=0;
+  poisson_time = 0.0;
 }
 
 /* ----------------------------------------------------------------------
@@ -164,10 +195,8 @@ void PPPMGPU::compute(int eflag, int vflag)
     if (atom->nlocal > old_nlocal) {
       nago=0;
       old_nlocal = atom->nlocal;
-    } else
+    } else nago = 1;
-      nago=1;
+  } else nago = neighbor->ago;
-  } else
-    nago=neighbor->ago;
 
   // set energy/virial flags
   // invoke allocate_peratom() if needed for first time
@@ -211,12 +240,13 @@ void PPPMGPU::compute(int eflag, int vflag)
     domain->x2lamda(atom->nlocal);
   }
 
-  double t3=MPI_Wtime();
+  double t3 = MPI_Wtime();
 
   // all procs communicate density values from their ghost cells
   //   to fully sum contribution in their 3d bricks
   // remap from 3d decomposition to FFT decomposition
 
+  cg->reverse_comm(this,REVERSE_RHO);
   brick2fft();
 
   // compute potential gradient on my FFT grid and
@@ -225,23 +255,32 @@ void PPPMGPU::compute(int eflag, int vflag)
 
   poisson();
 
-  // all procs communicate E-field values to fill ghost cells
+  // all procs communicate E-field values
-  //   surrounding their 3d bricks
+  // to fill ghost cells surrounding their 3d bricks
 
-  fillbrick();
+  if (differentiation_flag == 1) cg->forward_comm(this,FORWARD_AD);
+  else cg->forward_comm(this,FORWARD_IK);
 
-  poisson_time+=MPI_Wtime()-t3;
+  // extra per-atom energy/virial communication
+
+  if (evflag_atom) {
+    if (differentiation_flag == 1 && vflag_atom) 
+      cg_peratom->forward_comm(this,FORWARD_AD_PERATOM);
+    else if (differentiation_flag == 0)
+      cg_peratom->forward_comm(this,FORWARD_IK_PERATOM);
+  }
+
+  poisson_time += MPI_Wtime()-t3;
 
   // calculate the force on my particles
 
   FFT_SCALAR qscale = force->qqrd2e * scale;
   PPPM_GPU_API(interp)(qscale);
 
-  // Compute per-atom energy/virial on host if requested
+  // per-atom energy/virial
+  // energy includes self-energy correction
 
   if (evflag_atom) {
-    fillbrick_peratom();
-    fieldforce_peratom();
     double *q = atom->q;
     int nlocal = atom->nlocal;
 
@@ -293,246 +332,13 @@ void PPPMGPU::compute(int eflag, int vflag)
 }
 
 /* ----------------------------------------------------------------------
-   allocate memory that depends on # of K-vectors and order
-------------------------------------------------------------------------- */
-
-void PPPMGPU::allocate()
-{
-  memory->create(density_fft,nfft_both,"pppm:density_fft");
-  memory->create(greensfn,nfft_both,"pppm:greensfn");
-  memory->create(work1,2*nfft_both,"pppm:work1");
-  memory->create(work2,2*nfft_both,"pppm:work2");
-  memory->create(vg,nfft_both,6,"pppm:vg");
-
-  memory->create1d_offset(fkx,nxlo_fft,nxhi_fft,"pppm:fkx");
-  memory->create1d_offset(fky,nylo_fft,nyhi_fft,"pppm:fky");
-  memory->create1d_offset(fkz,nzlo_fft,nzhi_fft,"pppm:fkz");
-
-  memory->create(buf1,nbuf,"pppm:buf1");
-  memory->create(buf2,nbuf,"pppm:buf2");
-
-  // summation coeffs
-
-  memory->create(gf_b,order,"pppm:gf_b");
-  memory->create2d_offset(rho1d,3,-order/2,order/2,"pppm:rho1d");
-  memory->create2d_offset(drho1d,3,-order/2,order/2,"pppm:drho1d");
-  memory->create2d_offset(rho_coeff,order,(1-order)/2,order/2,"pppm:rho_coeff");
-  memory->create2d_offset(drho_coeff,order,(1-order)/2,order/2,
-                          "pppm:drho_coeff");
-
-  // create 2 FFTs and a Remap
-  // 1st FFT keeps data in FFT decompostion
-  // 2nd FFT returns data in 3d brick decomposition
-  // remap takes data from 3d brick to FFT decomposition
-
-  int tmp;
-
-  fft1 = new FFT3d(lmp,world,nx_pppm,ny_pppm,nz_pppm,
-                   nxlo_fft,nxhi_fft,nylo_fft,nyhi_fft,nzlo_fft,nzhi_fft,
-                   nxlo_fft,nxhi_fft,nylo_fft,nyhi_fft,nzlo_fft,nzhi_fft,
-                   0,0,&tmp);
-
-  fft2 = new FFT3d(lmp,world,nx_pppm,ny_pppm,nz_pppm,
-                   nxlo_fft,nxhi_fft,nylo_fft,nyhi_fft,nzlo_fft,nzhi_fft,
-                   nxlo_in,nxhi_in,nylo_in,nyhi_in,nzlo_in,nzhi_in,
-                   0,0,&tmp);
-
-  remap = new Remap(lmp,world,
-                    nxlo_in,nxhi_in,nylo_in,nyhi_in,nzlo_in,nzhi_in,
-                    nxlo_fft,nxhi_fft,nylo_fft,nyhi_fft,nzlo_fft,nzhi_fft,
-                    1,0,0,FFT_PRECISION);
-}
-
-/* ----------------------------------------------------------------------
-   deallocate memory that depends on # of K-vectors and order
-------------------------------------------------------------------------- */
-
-void PPPMGPU::deallocate()
-{
-  destroy_3d_offset(density_brick_gpu,nzlo_out,nylo_out);
-  destroy_3d_offset(vd_brick,nzlo_out,nylo_out);
-
-  memory->destroy(density_fft);
-  memory->destroy(greensfn);
-  memory->destroy(work1);
-  memory->destroy(work2);
-  memory->destroy(vg);
-
-  memory->destroy1d_offset(fkx,nxlo_fft);
-  memory->destroy1d_offset(fky,nylo_fft);
-  memory->destroy1d_offset(fkz,nzlo_fft);
-
-  memory->destroy(buf1);
-  memory->destroy(buf2);
-
-  memory->destroy(gf_b);
-  memory->destroy2d_offset(rho1d,-order/2);
-  memory->destroy2d_offset(drho1d,-order/2);
-  memory->destroy2d_offset(rho_coeff,(1-order)/2);
-  memory->destroy2d_offset(drho_coeff,(1-order)/2);
-
-  delete fft1;
-  delete fft2;
-  delete remap;
-}
-
-
-/* ----------------------------------------------------------------------
-   ghost-swap to accumulate full density in brick decomposition
    remap density from 3d brick decomposition to FFT decomposition
 ------------------------------------------------------------------------- */
 
 void PPPMGPU::brick2fft()
 {
-  int i,n,ix,iy,iz;
+  int n,ix,iy,iz;
-  MPI_Request request;
-  MPI_Status status;
 
-  // pack my ghosts for +x processor
-  // pass data to self or +x processor
-  // unpack and sum recv data into my real cells
-
-  n = 0;
-  for (iz = nzlo_out; iz <= nzhi_out; iz++)
-    for (iy = nylo_out; iy <= nyhi_out; iy++)
-      for (ix = nxhi_in+1; ix <= nxhi_out; ix++)
-        buf1[n++] = density_brick_gpu[iz][iy][ix];
-
-  if (comm->procneigh[0][1] == me)
-    for (i = 0; i < n; i++) buf2[i] = buf1[i];
-  else {
-    MPI_Irecv(buf2,nbuf,MPI_FFT_SCALAR,comm->procneigh[0][0],0,world,&request);
-    MPI_Send(buf1,n,MPI_FFT_SCALAR,comm->procneigh[0][1],0,world);
-    MPI_Wait(&request,&status);
-  }
-
-  n = 0;
-  for (iz = nzlo_out; iz <= nzhi_out; iz++)
-    for (iy = nylo_out; iy <= nyhi_out; iy++)
-      for (ix = nxlo_in; ix < nxlo_in+nxlo_ghost; ix++)
-        density_brick_gpu[iz][iy][ix] += buf2[n++];
-
-  // pack my ghosts for -x processor
-  // pass data to self or -x processor
-  // unpack and sum recv data into my real cells
-
-  n = 0;
-  for (iz = nzlo_out; iz <= nzhi_out; iz++)
-    for (iy = nylo_out; iy <= nyhi_out; iy++)
-      for (ix = nxlo_out; ix < nxlo_in; ix++)
-        buf1[n++] = density_brick_gpu[iz][iy][ix];
-
-  if (comm->procneigh[0][0] == me)
-    for (i = 0; i < n; i++) buf2[i] = buf1[i];
-  else {
-    MPI_Irecv(buf2,nbuf,MPI_FFT_SCALAR,comm->procneigh[0][1],0,world,&request);
-    MPI_Send(buf1,n,MPI_FFT_SCALAR,comm->procneigh[0][0],0,world);
-    MPI_Wait(&request,&status);
-  }
-
-  n = 0;
-  for (iz = nzlo_out; iz <= nzhi_out; iz++)
-    for (iy = nylo_out; iy <= nyhi_out; iy++)
-      for (ix = nxhi_in-nxhi_ghost+1; ix <= nxhi_in; ix++)
-        density_brick_gpu[iz][iy][ix] += buf2[n++];
-
-  // pack my ghosts for +y processor
-  // pass data to self or +y processor
-  // unpack and sum recv data into my real cells
-
-  n = 0;
-  for (iz = nzlo_out; iz <= nzhi_out; iz++)
-    for (iy = nyhi_in+1; iy <= nyhi_out; iy++)
-      for (ix = nxlo_in; ix <= nxhi_in; ix++)
-        buf1[n++] = density_brick_gpu[iz][iy][ix];
-
-  if (comm->procneigh[1][1] == me)
-    for (i = 0; i < n; i++) buf2[i] = buf1[i];
-  else {
-    MPI_Irecv(buf2,nbuf,MPI_FFT_SCALAR,comm->procneigh[1][0],0,world,&request);
-    MPI_Send(buf1,n,MPI_FFT_SCALAR,comm->procneigh[1][1],0,world);
-    MPI_Wait(&request,&status);
-  }
-
-  n = 0;
-  for (iz = nzlo_out; iz <= nzhi_out; iz++)
-    for (iy = nylo_in; iy < nylo_in+nylo_ghost; iy++)
-      for (ix = nxlo_in; ix <= nxhi_in; ix++)
-        density_brick_gpu[iz][iy][ix] += buf2[n++];
-
-  // pack my ghosts for -y processor
-  // pass data to self or -y processor
-  // unpack and sum recv data into my real cells
-
-  n = 0;
-  for (iz = nzlo_out; iz <= nzhi_out; iz++)
-    for (iy = nylo_out; iy < nylo_in; iy++)
-      for (ix = nxlo_in; ix <= nxhi_in; ix++)
-        buf1[n++] = density_brick_gpu[iz][iy][ix];
-
-  if (comm->procneigh[1][0] == me)
-    for (i = 0; i < n; i++) buf2[i] = buf1[i];
-  else {
-    MPI_Irecv(buf2,nbuf,MPI_FFT_SCALAR,comm->procneigh[1][1],0,world,&request);
-    MPI_Send(buf1,n,MPI_FFT_SCALAR,comm->procneigh[1][0],0,world);
-    MPI_Wait(&request,&status);
-  }
-
-  n = 0;
-  for (iz = nzlo_out; iz <= nzhi_out; iz++)
-    for (iy = nyhi_in-nyhi_ghost+1; iy <= nyhi_in; iy++)
-      for (ix = nxlo_in; ix <= nxhi_in; ix++)
-        density_brick_gpu[iz][iy][ix] += buf2[n++];
-
-  // pack my ghosts for +z processor
-  // pass data to self or +z processor
-  // unpack and sum recv data into my real cells
-
-  n = 0;
-  for (iz = nzhi_in+1; iz <= nzhi_out; iz++)
-    for (iy = nylo_in; iy <= nyhi_in; iy++)
-      for (ix = nxlo_in; ix <= nxhi_in; ix++)
-        buf1[n++] = density_brick_gpu[iz][iy][ix];
-
-  if (comm->procneigh[2][1] == me)
-    for (i = 0; i < n; i++) buf2[i] = buf1[i];
-  else {
-    MPI_Irecv(buf2,nbuf,MPI_FFT_SCALAR,comm->procneigh[2][0],0,world,&request);
-    MPI_Send(buf1,n,MPI_FFT_SCALAR,comm->procneigh[2][1],0,world);
-    MPI_Wait(&request,&status);
-  }
-
-  n = 0;
-  for (iz = nzlo_in; iz < nzlo_in+nzlo_ghost; iz++)
-    for (iy = nylo_in; iy <= nyhi_in; iy++)
-      for (ix = nxlo_in; ix <= nxhi_in; ix++)
-        density_brick_gpu[iz][iy][ix] += buf2[n++];
-
-  // pack my ghosts for -z processor
-  // pass data to self or -z processor
-  // unpack and sum recv data into my real cells
-
-  n = 0;
-  for (iz = nzlo_out; iz < nzlo_in; iz++)
-    for (iy = nylo_in; iy <= nyhi_in; iy++)
-      for (ix = nxlo_in; ix <= nxhi_in; ix++)
-        buf1[n++] = density_brick_gpu[iz][iy][ix];
-
-  if (comm->procneigh[2][0] == me)
-    for (i = 0; i < n; i++) buf2[i] = buf1[i];
-  else {
-    MPI_Irecv(buf2,nbuf,MPI_FFT_SCALAR,comm->procneigh[2][1],0,world,&request);
-    MPI_Send(buf1,n,MPI_FFT_SCALAR,comm->procneigh[2][0],0,world);
-    MPI_Wait(&request,&status);
-  }
-
-  n = 0;
-  for (iz = nzhi_in-nzhi_ghost+1; iz <= nzhi_in; iz++)
-    for (iy = nylo_in; iy <= nyhi_in; iy++)
-      for (ix = nxlo_in; ix <= nxhi_in; ix++)
-        density_brick_gpu[iz][iy][ix] += buf2[n++];
-
-  // remap from 3d brick decomposition to FFT decomposition
   // copy grabs inner portion of density from 3d brick
   // remap could be done as pre-stage of FFT,
   //   but this works optimally on only double values, not complex values
@@ -546,216 +352,6 @@ void PPPMGPU::brick2fft()
   remap->perform(density_fft,density_fft,work1);
 }
 
-/* ----------------------------------------------------------------------
-   Same as base class - needed to call GPU version of fillbrick_.
-------------------------------------------------------------------------- */
-
-void PPPMGPU::fillbrick()
-{
-  if (differentiation_flag == 1) fillbrick_ad();
-  else fillbrick_ik();
-}
-
-/* ----------------------------------------------------------------------
-   ghost-swap to fill ghost cells of my brick with field values
-------------------------------------------------------------------------- */
-
-void PPPMGPU::fillbrick_ik()
-{
-  int i,n,ix,iy,iz;
-  MPI_Request request;
-  MPI_Status status;
-
-  // pack my real cells for +z processor
-  // pass data to self or +z processor
-  // unpack and sum recv data into my ghost cells
-
-  n = 0;
-  int x_lo = nxlo_in * 4;
-  int x_hi = nxhi_in * 4 + 1;
-  for (iz = nzhi_in-nzhi_ghost+1; iz <= nzhi_in; iz++)
-    for (iy = nylo_in; iy <= nyhi_in; iy++)
-      for (ix = x_lo; ix < x_hi; ix+=4) {
-        buf1[n++] = vd_brick[iz][iy][ix];
-        buf1[n++] = vd_brick[iz][iy][ix+1];
-        buf1[n++] = vd_brick[iz][iy][ix+2];
-      }
-
-  if (comm->procneigh[2][1] == me)
-    for (i = 0; i < n; i++) buf2[i] = buf1[i];
-  else {
-    MPI_Irecv(buf2,nbuf,MPI_FFT_SCALAR,comm->procneigh[2][0],0,world,&request);
-    MPI_Send(buf1,n,MPI_FFT_SCALAR,comm->procneigh[2][1],0,world);
-    MPI_Wait(&request,&status);
-  }
-
-  n = 0;
-  for (iz = nzlo_out; iz < nzlo_in; iz++)
-    for (iy = nylo_in; iy <= nyhi_in; iy++)
-      for (ix = x_lo; ix < x_hi; ix+=4) {
-        vd_brick[iz][iy][ix] = buf2[n++];
-        vd_brick[iz][iy][ix+1] = buf2[n++];
-        vd_brick[iz][iy][ix+2] = buf2[n++];
-      }
-
-  // pack my real cells for -z processor
-  // pass data to self or -z processor
-  // unpack and sum recv data into my ghost cells
-
-  n = 0;
-  for (iz = nzlo_in; iz < nzlo_in+nzlo_ghost; iz++)
-    for (iy = nylo_in; iy <= nyhi_in; iy++)
-      for (ix = x_lo; ix < x_hi; ix+=4) {
-        buf1[n++] = vd_brick[iz][iy][ix];
-        buf1[n++] = vd_brick[iz][iy][ix+1];
-        buf1[n++] = vd_brick[iz][iy][ix+2];
-      }
-
-  if (comm->procneigh[2][0] == me)
-    for (i = 0; i < n; i++) buf2[i] = buf1[i];
-  else {
-    MPI_Irecv(buf2,nbuf,MPI_FFT_SCALAR,comm->procneigh[2][1],0,world,&request);
-    MPI_Send(buf1,n,MPI_FFT_SCALAR,comm->procneigh[2][0],0,world);
-    MPI_Wait(&request,&status);
-  }
-
-  n = 0;
-  for (iz = nzhi_in+1; iz <= nzhi_out; iz++)
-    for (iy = nylo_in; iy <= nyhi_in; iy++)
-      for (ix = x_lo; ix < x_hi; ix+=4) {
-        vd_brick[iz][iy][ix] = buf2[n++];
-        vd_brick[iz][iy][ix+1] = buf2[n++];
-        vd_brick[iz][iy][ix+2] = buf2[n++];
-      }
-
-  // pack my real cells for +y processor
-  // pass data to self or +y processor
-  // unpack and sum recv data into my ghost cells
-
-  n = 0;
-  for (iz = nzlo_out; iz <= nzhi_out; iz++)
-    for (iy = nyhi_in-nyhi_ghost+1; iy <= nyhi_in; iy++)
-      for (ix = x_lo; ix < x_hi; ix+=4) {
-        buf1[n++] = vd_brick[iz][iy][ix];
-        buf1[n++] = vd_brick[iz][iy][ix+1];
-        buf1[n++] = vd_brick[iz][iy][ix+2];
-      }
-
-  if (comm->procneigh[1][1] == me)
-    for (i = 0; i < n; i++) buf2[i] = buf1[i];
-  else {
-    MPI_Irecv(buf2,nbuf,MPI_FFT_SCALAR,comm->procneigh[1][0],0,world,&request);
-    MPI_Send(buf1,n,MPI_FFT_SCALAR,comm->procneigh[1][1],0,world);
-    MPI_Wait(&request,&status);
-  }
-
-  n = 0;
-  for (iz = nzlo_out; iz <= nzhi_out; iz++)
-    for (iy = nylo_out; iy < nylo_in; iy++)
-      for (ix = x_lo; ix < x_hi; ix+=4) {
-        vd_brick[iz][iy][ix] = buf2[n++];
-        vd_brick[iz][iy][ix+1] = buf2[n++];
-        vd_brick[iz][iy][ix+2] = buf2[n++];
-      }
-
-  // pack my real cells for -y processor
-  // pass data to self or -y processor
-  // unpack and sum recv data into my ghost cells
-
-  n = 0;
-  for (iz = nzlo_out; iz <= nzhi_out; iz++)
-    for (iy = nylo_in; iy < nylo_in+nylo_ghost; iy++)
-      for (ix = x_lo; ix < x_hi; ix+=4) {
-        buf1[n++] = vd_brick[iz][iy][ix];
-        buf1[n++] = vd_brick[iz][iy][ix+1];
-        buf1[n++] = vd_brick[iz][iy][ix+2];
-      }
-
-  if (comm->procneigh[1][0] == me)
-    for (i = 0; i < n; i++) buf2[i] = buf1[i];
-  else {
-    MPI_Irecv(buf2,nbuf,MPI_FFT_SCALAR,comm->procneigh[1][1],0,world,&request);
-    MPI_Send(buf1,n,MPI_FFT_SCALAR,comm->procneigh[1][0],0,world);
-    MPI_Wait(&request,&status);
-  }
-
-  n = 0;
-  for (iz = nzlo_out; iz <= nzhi_out; iz++)
-    for (iy = nyhi_in+1; iy <= nyhi_out; iy++)
-      for (ix = x_lo; ix < x_hi; ix+=4) {
-        vd_brick[iz][iy][ix] = buf2[n++];
-        vd_brick[iz][iy][ix+1] = buf2[n++];
-        vd_brick[iz][iy][ix+2] = buf2[n++];
-      }
-
-  // pack my real cells for +x processor
-  // pass data to self or +x processor
-  // unpack and sum recv data into my ghost cells
-
-  n = 0;
-  x_lo = (nxhi_in-nxhi_ghost+1) * 4;
-  for (iz = nzlo_out; iz <= nzhi_out; iz++)
-    for (iy = nylo_out; iy <= nyhi_out; iy++)
-      for (ix = x_lo; ix < x_hi; ix+=4) {
-        buf1[n++] = vd_brick[iz][iy][ix];
-        buf1[n++] = vd_brick[iz][iy][ix+1];
-        buf1[n++] = vd_brick[iz][iy][ix+2];
-      }
-
-  if (comm->procneigh[0][1] == me)
-    for (i = 0; i < n; i++) buf2[i] = buf1[i];
-  else {
-    MPI_Irecv(buf2,nbuf,MPI_FFT_SCALAR,comm->procneigh[0][0],0,world,&request);
-    MPI_Send(buf1,n,MPI_FFT_SCALAR,comm->procneigh[0][1],0,world);
-    MPI_Wait(&request,&status);
-  }
-
-  n = 0;
-  x_lo = nxlo_out * 4;
-  x_hi = nxlo_in * 4;
-  for (iz = nzlo_out; iz <= nzhi_out; iz++)
-    for (iy = nylo_out; iy <= nyhi_out; iy++)
-      for (ix = x_lo; ix < x_hi; ix+=4) {
-        vd_brick[iz][iy][ix] = buf2[n++];
-        vd_brick[iz][iy][ix+1] = buf2[n++];
-        vd_brick[iz][iy][ix+2] = buf2[n++];
-      }
-
-  // pack my real cells for -x processor
-  // pass data to self or -x processor
-  // unpack and sum recv data into my ghost cells
-
-  n = 0;
-  x_lo = x_hi;
-  x_hi = (nxlo_in+nxlo_ghost) * 4;
-  for (iz = nzlo_out; iz <= nzhi_out; iz++)
-    for (iy = nylo_out; iy <= nyhi_out; iy++)
-      for (ix = x_lo; ix < x_hi; ix+=4) {
-        buf1[n++] = vd_brick[iz][iy][ix];
-        buf1[n++] = vd_brick[iz][iy][ix+1];
-        buf1[n++] = vd_brick[iz][iy][ix+2];
-      }
-
-  if (comm->procneigh[0][0] == me)
-    for (i = 0; i < n; i++) buf2[i] = buf1[i];
-  else {
-    MPI_Irecv(buf2,nbuf,MPI_FFT_SCALAR,comm->procneigh[0][1],0,world,&request);
-    MPI_Send(buf1,n,MPI_FFT_SCALAR,comm->procneigh[0][0],0,world);
-    MPI_Wait(&request,&status);
-  }
-
-  n = 0;
-  x_lo = (nxhi_in + 1) * 4;
-  x_hi = nxhi_out * 4 + 1;
-  for (iz = nzlo_out; iz <= nzhi_out; iz++)
-    for (iy = nylo_out; iy <= nyhi_out; iy++)
-      for (ix = x_lo; ix < x_hi; ix+=4) {
-        vd_brick[iz][iy][ix] = buf2[n++];
-        vd_brick[iz][iy][ix+1] = buf2[n++];
-        vd_brick[iz][iy][ix+2] = buf2[n++];
-      }
-}
-
 /* ----------------------------------------------------------------------
    Same code as base class - necessary to call GPU version of poisson_ik
 ------------------------------------------------------------------------- */
@@ -892,12 +488,154 @@ void PPPMGPU::poisson_ik()
 }
 
 /* ----------------------------------------------------------------------
-   Create array using offsets from pinned memory allocation
+   pack own values to buf to send to another proc
+------------------------------------------------------------------------- */
+
+void PPPMGPU::pack_forward(int flag, FFT_SCALAR *buf, int nlist, int *list)
+{
+  int n = 0;
+
+  if (flag == FORWARD_IK) {
+    int offset;
+    FFT_SCALAR *src = &vd_brick[nzlo_out][nylo_out][4*nxlo_out];
+    for (int i = 0; i < nlist; i++) {
+      offset = 4*list[i];
+      buf[n++] = src[offset++];
+      buf[n++] = src[offset++];
+      buf[n++] = src[offset];
+    }
+  } else if (flag == FORWARD_AD) {
+    FFT_SCALAR *src = &u_brick[nzlo_out][nylo_out][nxlo_out];
+    for (int i = 0; i < nlist; i++)
+      buf[i] = src[list[i]];
+  } else if (flag == FORWARD_IK_PERATOM) {
+    FFT_SCALAR *esrc = &u_brick[nzlo_out][nylo_out][nxlo_out];
+    FFT_SCALAR *v0src = &v0_brick[nzlo_out][nylo_out][nxlo_out];
+    FFT_SCALAR *v1src = &v1_brick[nzlo_out][nylo_out][nxlo_out];
+    FFT_SCALAR *v2src = &v2_brick[nzlo_out][nylo_out][nxlo_out];
+    FFT_SCALAR *v3src = &v3_brick[nzlo_out][nylo_out][nxlo_out];
+    FFT_SCALAR *v4src = &v4_brick[nzlo_out][nylo_out][nxlo_out];
+    FFT_SCALAR *v5src = &v5_brick[nzlo_out][nylo_out][nxlo_out];
+    for (int i = 0; i < nlist; i++) {
+      if (eflag_atom) buf[n++] = esrc[list[i]];
+      if (vflag_atom) {
+        buf[n++] = v0src[list[i]];
+        buf[n++] = v1src[list[i]];
+        buf[n++] = v2src[list[i]];
+        buf[n++] = v3src[list[i]];
+        buf[n++] = v4src[list[i]];
+        buf[n++] = v5src[list[i]];
+      }
+    }
+  } else if (flag == FORWARD_AD_PERATOM) {
+    FFT_SCALAR *v0src = &v0_brick[nzlo_out][nylo_out][nxlo_out];
+    FFT_SCALAR *v1src = &v1_brick[nzlo_out][nylo_out][nxlo_out];
+    FFT_SCALAR *v2src = &v2_brick[nzlo_out][nylo_out][nxlo_out];
+    FFT_SCALAR *v3src = &v3_brick[nzlo_out][nylo_out][nxlo_out];
+    FFT_SCALAR *v4src = &v4_brick[nzlo_out][nylo_out][nxlo_out];
+    FFT_SCALAR *v5src = &v5_brick[nzlo_out][nylo_out][nxlo_out];
+    for (int i = 0; i < nlist; i++) {
+      buf[n++] = v0src[list[i]];
+      buf[n++] = v1src[list[i]];
+      buf[n++] = v2src[list[i]];
+      buf[n++] = v3src[list[i]];
+      buf[n++] = v4src[list[i]];
+      buf[n++] = v5src[list[i]];
+    }
+  }
+}
+
+/* ----------------------------------------------------------------------
+   unpack another proc's own values from buf and set own ghost values
+------------------------------------------------------------------------- */
+
+void PPPMGPU::unpack_forward(int flag, FFT_SCALAR *buf, int nlist, int *list)
+{
+  int n = 0;
+
+  if (flag == FORWARD_IK) {
+    int offset;
+    FFT_SCALAR *dest = &vdx_brick[nzlo_out][nylo_out][4*nxlo_out];
+    for (int i = 0; i < nlist; i++) {
+      offset = 4*list[i];
+      dest[offset++] = buf[n++];
+      dest[offset++] = buf[n++];
+      dest[offset] = buf[n++];
+    }
+  } else if (flag == FORWARD_AD) {
+    FFT_SCALAR *dest = &u_brick[nzlo_out][nylo_out][nxlo_out];
+    for (int i = 0; i < nlist; i++)
+      dest[list[i]] = buf[n++];
+  } else if (flag == FORWARD_IK_PERATOM) {
+    FFT_SCALAR *esrc = &u_brick[nzlo_out][nylo_out][nxlo_out];
+    FFT_SCALAR *v0src = &v0_brick[nzlo_out][nylo_out][nxlo_out];
+    FFT_SCALAR *v1src = &v1_brick[nzlo_out][nylo_out][nxlo_out];
+    FFT_SCALAR *v2src = &v2_brick[nzlo_out][nylo_out][nxlo_out];
+    FFT_SCALAR *v3src = &v3_brick[nzlo_out][nylo_out][nxlo_out];
+    FFT_SCALAR *v4src = &v4_brick[nzlo_out][nylo_out][nxlo_out];
+    FFT_SCALAR *v5src = &v5_brick[nzlo_out][nylo_out][nxlo_out];
+    for (int i = 0; i < nlist; i++) {
+      if (eflag_atom) esrc[list[i]] = buf[n++];
+      if (vflag_atom) {
+        v0src[list[i]] = buf[n++];
+        v1src[list[i]] = buf[n++];
+        v2src[list[i]] = buf[n++];
+        v3src[list[i]] = buf[n++];
+        v4src[list[i]] = buf[n++];
+        v5src[list[i]] = buf[n++];
+      }
+    }
+  } else if (flag == FORWARD_AD_PERATOM) {
+    FFT_SCALAR *v0src = &v0_brick[nzlo_out][nylo_out][nxlo_out];
+    FFT_SCALAR *v1src = &v1_brick[nzlo_out][nylo_out][nxlo_out];
+    FFT_SCALAR *v2src = &v2_brick[nzlo_out][nylo_out][nxlo_out];
+    FFT_SCALAR *v3src = &v3_brick[nzlo_out][nylo_out][nxlo_out];
+    FFT_SCALAR *v4src = &v4_brick[nzlo_out][nylo_out][nxlo_out];
+    FFT_SCALAR *v5src = &v5_brick[nzlo_out][nylo_out][nxlo_out];
+    for (int i = 0; i < nlist; i++) {
+      v0src[list[i]] = buf[n++];
+      v1src[list[i]] = buf[n++];
+      v2src[list[i]] = buf[n++];
+      v3src[list[i]] = buf[n++];
+      v4src[list[i]] = buf[n++];
+      v5src[list[i]] = buf[n++];
+    }
+  }
+}
+
+/* ----------------------------------------------------------------------
+   pack ghost values into buf to send to another proc
+------------------------------------------------------------------------- */
+
+void PPPMGPU::pack_reverse(int flag, FFT_SCALAR *buf, int nlist, int *list)
+{
+  if (flag == REVERSE_RHO) {
+    FFT_SCALAR *src = &density_brick_gpu[nzlo_out][nylo_out][nxlo_out];
+    for (int i = 0; i < nlist; i++)
+      buf[i] = src[list[i]];
+  }
+}
+
+/* ----------------------------------------------------------------------
+   unpack another proc's ghost values from buf and add to own values
+------------------------------------------------------------------------- */
+
+void PPPMGPU::unpack_reverse(int flag, FFT_SCALAR *buf, int nlist, int *list)
+{
+  if (flag == REVERSE_RHO) {
+    FFT_SCALAR *dest = &density_brick_gpu[nzlo_out][nylo_out][nxlo_out];
+    for (int i = 0; i < nlist; i++)
+      dest[list[i]] += buf[i];
+  } 
+}
+
+/* ----------------------------------------------------------------------
+   create array using offsets from pinned memory allocation
 ------------------------------------------------------------------------- */
 
 FFT_SCALAR ***PPPMGPU::create_3d_offset(int n1lo, int n1hi, int n2lo, int n2hi,
-                                     int n3lo, int n3hi, const char *name,
+                                        int n3lo, int n3hi, const char *name,
-                                     FFT_SCALAR *data, int vec_length)
+                                        FFT_SCALAR *data, int vec_length)
 {
   int i,j;
   int n1 = n1hi - n1lo + 1;
@@ -942,19 +680,11 @@ void PPPMGPU::destroy_3d_offset(FFT_SCALAR ***array, int n1_offset,
 
 double PPPMGPU::memory_usage()
 {
-  double bytes = nmax*3 * sizeof(double);
+  double bytes = PPPM::memory_usage();
-  int nbrick = (nxhi_out-nxlo_out+1) * (nyhi_out-nylo_out+1) *
-    (nzhi_out-nzlo_out+1);
-  bytes += 4 * nbrick * sizeof(FFT_SCALAR);
-  bytes += 6 * nfft_both * sizeof(double);
-  bytes += nfft_both * sizeof(double);
-  bytes += nfft_both*5 * sizeof(FFT_SCALAR);
-  bytes += 2 * nbuf * sizeof(double);
 
-  if (peratom_allocate_flag) {
+  // NOTE: add tallying here for density_brick_gpu and vd_brick
-    bytes += 7 * nbrick * sizeof(FFT_SCALAR);
+  //       could subtract out density_brick and vdxyz_brick if freed them above
-    bytes += 2 * nbuf_peratom * sizeof(FFT_SCALAR);
+  //       it the net efffect is zero, do nothing
-  }
 
   return bytes + PPPM_GPU_API(bytes)();
 }

src/GPU/pppm_gpu.h

@@ -35,20 +35,19 @@ class PPPMGPU : public PPPM {
   double memory_usage();
 
  protected:
-
   FFT_SCALAR ***density_brick_gpu, ***vd_brick;
   bool kspace_split, im_real_space;
+  int old_nlocal;
+  double poisson_time;
 
-  void allocate();
-  void deallocate();
   void brick2fft();
-  void fillbrick();
-  void fillbrick_ik();
   void poisson();
   void poisson_ik();
 
-  int old_nlocal;
+  void pack_forward(int, FFT_SCALAR *, int, int *);
-  double poisson_time;
+  void unpack_forward(int, FFT_SCALAR *, int, int *);
+  void pack_reverse(int, FFT_SCALAR *, int, int *);
+  void unpack_reverse(int, FFT_SCALAR *, int, int *);
 
   FFT_SCALAR ***create_3d_offset(int, int, int, int, int, int, const char *,
                                  FFT_SCALAR *, int);

src/GRANULAR/pair_gran_hooke_history.h

@@ -41,7 +41,7 @@ class PairGranHookeHistory : public Pair {
   void write_restart_settings(FILE *);
   void read_restart_settings(FILE *);
   void reset_dt();
-  double single(int, int, int, int, double, double, double, double &);
+  virtual double single(int, int, int, int, double, double, double, double &);
   int pack_comm(int, int *, double *, int, int *);
   void unpack_comm(int, int, double *);
   void *extract(const char *, int &);

src/KSPACE/fft3d.h

@@ -314,14 +314,18 @@ struct fft_plan_3d {
 
 // function prototypes
 
-void fft_3d(FFT_DATA *, FFT_DATA *, int, struct fft_plan_3d *);
+extern "C" { 
-struct fft_plan_3d *fft_3d_create_plan(MPI_Comm, int, int, int,
+  void fft_3d(FFT_DATA *, FFT_DATA *, int, struct fft_plan_3d *);
-  int, int, int, int, int, int, int, int, int, int, int, int,
+  struct fft_plan_3d *fft_3d_create_plan(MPI_Comm, int, int, int,
-  int, int, int *);
+                                         int, int, int, int, int, 
-void fft_3d_destroy_plan(struct fft_plan_3d *);
+                                         int, int, int, int, int, int, int,
-void factor(int, int *, int *);
+                                         int, int, int *);
-void bifactor(int, int *, int *);
+  void fft_3d_destroy_plan(struct fft_plan_3d *);
-void fft_1d_only(FFT_DATA *, int, int, struct fft_plan_3d *);
+  void factor(int, int *, int *);
+  void bifactor(int, int *, int *);
+  void fft_1d_only(FFT_DATA *, int, int, struct fft_plan_3d *);
+}
+
 /* ERROR/WARNING messages:
 
 */

src/KSPACE/pair_born_coul_long.cpp

@@ -158,6 +158,7 @@ void PairBornCoulLong::compute(int eflag, int vflag)
         } else forcecoul = 0.0;
 
         if (rsq < cut_ljsq[itype][jtype]) {
+          r = sqrt(rsq);
           r6inv = r2inv*r2inv*r2inv;
           rexp = exp((sigma[itype][jtype]-r)*rhoinv[itype][jtype]);
           forceborn = born1[itype][jtype]*r*rexp - born2[itype][jtype]*r6inv
@@ -381,7 +382,8 @@ void PairBornCoulLong::init_style()
   neighbor->request(this);
   
   // setup force tables
-  if (ncoultablebits) init_tables(cut_coul,cut_respa);
+
+  if (ncoultablebits) init_tables(cut_coul,NULL);
 }
 
 /* ----------------------------------------------------------------------

src/KSPACE/pair_buck_coul_long.cpp

@@ -150,6 +150,7 @@ void PairBuckCoulLong::compute(int eflag, int vflag)
         } else forcecoul = 0.0;
 
         if (rsq < cut_ljsq[itype][jtype]) {
+          r = sqrt(rsq);
           r6inv = r2inv*r2inv*r2inv;
           rexp = exp(-r*rhoinv[itype][jtype]);
           forcebuck = buck1[itype][jtype]*r*rexp - buck2[itype][jtype]*r6inv;
@@ -358,7 +359,8 @@ void PairBuckCoulLong::init_style()
   neighbor->request(this);
   
   // setup force tables
-  if (ncoultablebits) init_tables(cut_coul,cut_respa);
+
+  if (ncoultablebits) init_tables(cut_coul,NULL);
 }
 
 /* ----------------------------------------------------------------------

src/KSPACE/pair_coul_long.cpp

@@ -240,15 +240,6 @@ void PairCoulLong::init_style()
 
   cut_coulsq = cut_coul * cut_coul;
 
-  // set & error check interior rRESPA cutoffs
-
-  if (strstr(update->integrate_style,"respa") &&
-      ((Respa *) update->integrate)->level_inner >= 0) {
-    cut_respa = ((Respa *) update->integrate)->cutoff;
-    if (cut_coul < cut_respa[3])
-      error->all(FLERR,"Pair cutoff < Respa interior cutoff");
-  } else cut_respa = NULL;
-
   // insure use of KSpace long-range solver, set g_ewald
 
  if (force->kspace == NULL)
@@ -257,7 +248,7 @@ void PairCoulLong::init_style()
 
   // setup force tables
 
-  if (ncoultablebits) init_tables(cut_coul,cut_respa);
+  if (ncoultablebits) init_tables(cut_coul,NULL);
 }
 
 /* ----------------------------------------------------------------------

src/KSPACE/pppm.h

@@ -150,10 +150,10 @@ class PPPM : public KSpace {
 
   // grid communication
 
-  void pack_forward(int, FFT_SCALAR *, int, int *);
+  virtual void pack_forward(int, FFT_SCALAR *, int, int *);
-  void unpack_forward(int, FFT_SCALAR *, int, int *);
+  virtual void unpack_forward(int, FFT_SCALAR *, int, int *);
-  void pack_reverse(int, FFT_SCALAR *, int, int *);
+  virtual void pack_reverse(int, FFT_SCALAR *, int, int *);
-  void unpack_reverse(int, FFT_SCALAR *, int, int *);
+  virtual void unpack_reverse(int, FFT_SCALAR *, int, int *);
 
   // group-group interactions
 

src/MC/fix_gcmc.h

@@ -32,49 +32,77 @@ class FixGCMC : public Fix {
   int setmask();
   void init();
   void pre_exchange();
-  void attempt_move();
+  void attempt_atomic_translation();
-  void attempt_deletion();
+  void attempt_atomic_deletion();
-  void attempt_insertion();
+  void attempt_atomic_insertion();
-  double energy(int, double *);
+  void attempt_molecule_translation();
+  void attempt_molecule_rotation();
+  void attempt_molecule_deletion();
+  void attempt_molecule_insertion();
+  double energy(int, int, int, double *);
+  int pick_random_gas_atom();
+  int pick_random_gas_atom_in_region();
+  int pick_random_gas_molecule();
+  int pick_random_gas_molecule_in_region();
+  double molecule_energy(int);
+  void get_rotation_matrix(double, double *);
+  void get_model_molecule();
+  void update_gas_atoms_list();
   double compute_vector(int);
   double memory_usage();
   void write_restart(FILE *);
   void restart(char *);
 
  private:
-  int ntype,nevery,seed;
+  int rotation_group,rotation_groupbit;
+  int rotation_inversegroupbit;
+  int ngcmc_type,nevery,seed;
   int ncycles,nexchanges,nmcmoves;
-  int ngas;           // # of gas molecules (or atoms) on all procs
+  int ngas;                 // # of gas atoms on all procs
-  int ngas_local;     // # of gas molecules (or atoms) on this proc
+  int ngas_local;           // # of gas atoms on this proc
-  int ngas_before;    // # of gas molecules (or atoms) on procs < this proc
+  int ngas_before;          // # of gas atoms on procs < this proc
-  int molflag;        // 0 = atomic, 1 = molecular system
+  int molflag;              // 0 = atomic, 1 = molecular system
-  int regionflag;     // 0 = anywhere, 1 = specific region
+  int regionflag;           // 0 = anywhere in box, 1 = specific region
-  int iregion;        // exchange/move region
+  int iregion;              // GCMC region
-  char *idregion;     // exchange/move region id
+  char *idregion;           // GCMC region id
 
-  double nmove_attempts;
+  int maxmol;               // largest molecule tag across all existing atoms
-  double nmove_successes;
+  int natoms_per_molecule;  // number of atoms in each gas molecule
-  double ndel_attempts;
-  double ndel_successes;
-  double ninsert_attempts;
-  double ninsert_successes;
 
-  int nmax;
+  double ntranslation_attempts;
+  double ntranslation_successes;
+  double nrotation_attempts;
+  double nrotation_successes;
+  double ndeletion_attempts;
+  double ndeletion_successes;
+  double ninsertion_attempts;
+  double ninsertion_successes;
+
+  int gcmc_nmax;
   int max_region_attempts;
+  double gas_mass;
   double reservoir_temperature;
   double chemical_potential;
   double displace;
+  double max_rotation_angle;
   double beta,zz,sigma,volume;
   double xlo,xhi,ylo,yhi,zlo,zhi;
   double region_xlo,region_xhi,region_ylo,region_yhi,region_zlo,region_zhi;
+  double region_volume;
   double *sublo,*subhi;
   int *local_gas_list;
   double **cutsq;
+  double **atom_coord;
+  double *model_atom_buf;
+  tagint imagetmp;
+
   class Pair *pair;
 
   class RanPark *random_equal;
   class RanPark *random_unequal;
   
+  class Atom *model_atom;
+
   void options(int, char **);
 };
 
@@ -100,26 +128,68 @@ E: Cannot do GCMC on atoms in atom_modify first group
 This is a restriction due to the way atoms are organized in a list to
 enable the atom_modify first command.
 
-W: Fix GCMC may delete atom with non-zero molecule ID
+E: Fix GCMC cannot exchange individual atoms belonging to a molecule
 
-This is probably an error, since you should not delete only one atom
+This is an error since you should not delete only one atom of a molecule.
-of a molecule. The GCMC molecule exchange feature does not yet work.
+The user has specified atomic (non-molecular) gas exchanges, but an atom
+belonging to a molecule could be deleted.
 
-E: Fix GCMC molecule command requires atom attribute molecule
+E: All mol IDs should be set for fix GCMC group atoms
+
+The molecule flag is on, yet not all molecule ids in the fix group have
+been set to non-zero positive values by the user. This is an error since
+all atoms in the fix GCMC group are eligible for deletion, rotation, and
+translation and therefore must have valid molecule ids.
+
+E: Cannot use fix GCMC in a 2d simulation
+
+Fix GCMC is set up to run in 3d only. No 2d simulations with fix GCMC
+are allowed.
+
+E: Cannot use fix GCMC with a triclinic box
+
+Fix GCMC is set up to run with othogonal boxes only. Simulations with
+triclinic boxes and fix GCMC are not allowed.
+
+E: Fix GCMC molecule command requires that atoms have molecule attributes
 
 Should not choose the GCMC molecule feature if no molecules are being
 simulated. The general molecule flag is off, but GCMC's molecule flag
 is on.
 
-E: Fix GCMC molecule feature does not yet work
+E: Fix GCMC could not find any atoms in the user-supplied template molecule
 
-Fix GCMC cannot (yet) be used to exchange molecules, only atoms.
+When using the molecule option with fix GCMC, the user must supply a 
+template molecule in the usual LAMMPS data file with its molecule id
+specified in the fix GCMC command as the "type" of the exchanged gas.
 
 E: Fix GCMC incompatible with given pair_style
 
 Some pair_styles do not provide single-atom energies, which are needed
 by fix GCMC.
 
+E: Fix GCMC incorrect number of atoms per molecule
+
+The number of atoms in each gas molecule was not computed correctly.
+
+E: Illegal fix GCMC gas mass <= 0
+
+The computed mass of the designated gas molecule or atom type was less 
+than or equal to zero.
+
+E: Fix GCMC ran out of available molecule IDs
+
+This is a code limitation where more than MAXSMALLINT (usually around
+two billion) molecules have been created. The code needs to be 
+modified to either allow molecule ID recycling or use bigger ints for
+molecule IDs. A work-around is to run shorter simulations.
+
+W: Fix GCMC fix group should be all
+
+Fix GCMC will ignore the fix group specified by the user. User should
+set the fix group to "all". Fix GCMC will overwrite the user-specified
+fix group with a group consisting of all GCMC gas atoms.
+
 E: Fix GCMC region does not support a bounding box 
  
 Not all regions represent bounded volumes.  You cannot use 
@@ -129,7 +199,7 @@ E: Fix GCMC region cannot be dynamic
  
 Only static regions can be used with fix GCMC. 
  
-E: Deposition region extends outside simulation box
+E: Fix GCMC region extends outside simulation box 
  
 Self-explanatory. 
  

src/USER-MISC/angle_fourier_simple.cpp

@@ -13,7 +13,7 @@
 
 /* ----------------------------------------------------------------------
    Contributing author: Loukas D. Peristeras (Scienomics SARL)
-   [ based on angle_cosine_quared.cpp Naveen Michaud-Agrawal (Johns Hopkins U)]
+   [ based on angle_cosine_squared.cpp Naveen Michaud-Agrawal (Johns Hopkins U)]
 ------------------------------------------------------------------------- */
 
 #include "math.h"
@@ -160,7 +160,7 @@ void AngleFourierSimple::compute(int eflag, int vflag)
     }
 
     if (evflag) ev_tally(i1,i2,i3,nlocal,newton_bond,eangle,f1,f3,
-			 delx1,dely1,delz1,delx2,dely2,delz2);
+                         delx1,dely1,delz1,delx2,dely2,delz2);
   }
 }
 

src/USER-MISC/dihedral_fourier.cpp

@@ -104,26 +104,22 @@ void DihedralFourier::compute(int eflag, int vflag)
     vb1x = x[i1][0] - x[i2][0];
     vb1y = x[i1][1] - x[i2][1];
     vb1z = x[i1][2] - x[i2][2];
-    domain->minimum_image(vb1x,vb1y,vb1z);
 
     // 2nd bond
 
     vb2x = x[i3][0] - x[i2][0];
     vb2y = x[i3][1] - x[i2][1];
     vb2z = x[i3][2] - x[i2][2];
-    domain->minimum_image(vb2x,vb2y,vb2z);
 
     vb2xm = -vb2x;
     vb2ym = -vb2y;
     vb2zm = -vb2z;
-    domain->minimum_image(vb2xm,vb2ym,vb2zm);
 
     // 3rd bond
 
     vb3x = x[i4][0] - x[i3][0];
     vb3y = x[i4][1] - x[i3][1];
     vb3z = x[i4][2] - x[i3][2];
-    domain->minimum_image(vb3x,vb3y,vb3z);
 
     ax = vb1y*vb2zm - vb1z*vb2ym;
     ay = vb1z*vb2xm - vb1x*vb2zm;
@@ -152,19 +148,19 @@ void DihedralFourier::compute(int eflag, int vflag)
       int me;
       MPI_Comm_rank(world,&me);
       if (screen) {
-	char str[128];
+        char str[128];
-	sprintf(str,"Dihedral problem: %d " BIGINT_FORMAT " %d %d %d %d",
+        sprintf(str,"Dihedral problem: %d " BIGINT_FORMAT " %d %d %d %d",
-		me,update->ntimestep,
+                me,update->ntimestep,
-		atom->tag[i1],atom->tag[i2],atom->tag[i3],atom->tag[i4]);
+                atom->tag[i1],atom->tag[i2],atom->tag[i3],atom->tag[i4]);
-	error->warning(FLERR,str,0);
+        error->warning(FLERR,str,0);
-	fprintf(screen,"  1st atom: %d %g %g %g\n",
+        fprintf(screen,"  1st atom: %d %g %g %g\n",
-		me,x[i1][0],x[i1][1],x[i1][2]);
+                me,x[i1][0],x[i1][1],x[i1][2]);
-	fprintf(screen,"  2nd atom: %d %g %g %g\n",
+        fprintf(screen,"  2nd atom: %d %g %g %g\n",
-		me,x[i2][0],x[i2][1],x[i2][2]);
+                me,x[i2][0],x[i2][1],x[i2][2]);
-	fprintf(screen,"  3rd atom: %d %g %g %g\n",
+        fprintf(screen,"  3rd atom: %d %g %g %g\n",
-		me,x[i3][0],x[i3][1],x[i3][2]);
+                me,x[i3][0],x[i3][1],x[i3][2]);
-	fprintf(screen,"  4th atom: %d %g %g %g\n",
+        fprintf(screen,"  4th atom: %d %g %g %g\n",
-		me,x[i4][0],x[i4][1],x[i4][2]);
+                me,x[i4][0],x[i4][1],x[i4][2]);
       }
     }
 
@@ -268,7 +264,7 @@ void DihedralFourier::compute(int eflag, int vflag)
 
     if (evflag)
       ev_tally(i1,i2,i3,i4,nlocal,newton_bond,edihedral,f1,f3,f4,
-	       vb1x,vb1y,vb1z,vb2x,vb2y,vb2z,vb3x,vb3y,vb3z);
+               vb1x,vb1y,vb1z,vb2x,vb2y,vb2z,vb3x,vb3y,vb3z);
 
   }
 }
@@ -283,12 +279,12 @@ void DihedralFourier::allocate()
   memory->create(nterms,n+1,"dihedral:nterms");
   k = new double * [n+1];
   multiplicity = new int * [n+1];
-  shift = new int * [n+1];
+  shift = new double * [n+1];
   cos_shift = new double * [n+1];
   sin_shift = new double * [n+1];
   for (int i = 1; i <= n; i++) {
-    k[i] = cos_shift[i] = sin_shift[i] = 0;
+    k[i] = shift[i] = cos_shift[i] = sin_shift[i] = 0;
-    multiplicity[i] = shift[i] = 0;
+    multiplicity[i] = 0;
   }
   
   memory->create(setflag,n+1,"dihedral:setflag");
@@ -313,7 +309,7 @@ void DihedralFourier::coeff(int narg, char **arg)
 
   double k_one;
   int multiplicity_one;
-  int shift_one;
+  double shift_one;
   int nterms_one = force->inumeric(arg[1]);
 
   if (nterms_one < 1)
@@ -327,14 +323,14 @@ void DihedralFourier::coeff(int narg, char **arg)
     nterms[i] = nterms_one;
     k[i] = new double [nterms_one];
     multiplicity[i] = new int [nterms_one];
-    shift[i] = new int [nterms_one];
+    shift[i] = new double [nterms_one];
     cos_shift[i] = new double [nterms_one]; 
     sin_shift[i] = new double [nterms_one]; 
     for (int j = 0; j<nterms_one; j++) {
       int offset = 1+3*j;
       k_one = force->numeric(arg[offset+1]);
       multiplicity_one = force->inumeric(arg[offset+2]);
-      shift_one = force->inumeric(arg[offset+3]);
+      shift_one = force->numeric(arg[offset+3]);
       k[i][j] = k_one;
       multiplicity[i][j] = multiplicity_one;
       shift[i][j] = shift_one;
@@ -359,7 +355,7 @@ void DihedralFourier::write_restart(FILE *fp)
   for(int i = 1; i <= atom->ndihedraltypes; i++) {
     fwrite(k[i],sizeof(double),nterms[i],fp);
     fwrite(multiplicity[i],sizeof(int),nterms[i],fp);
-    fwrite(shift[i],sizeof(int),nterms[i],fp);
+    fwrite(shift[i],sizeof(double),nterms[i],fp);
   }
 
 }
@@ -381,7 +377,7 @@ void DihedralFourier::read_restart(FILE *fp)
   for (int i=1; i<=atom->ndihedraltypes; i++) {
     k[i] = new double [nterms[i]];
     multiplicity[i] = new int [nterms[i]];
-    shift[i] = new int [nterms[i]];
+    shift[i] = new double [nterms[i]];
     cos_shift[i] = new double [nterms[i]];
     sin_shift[i] = new double [nterms[i]];
   }
@@ -390,14 +386,14 @@ void DihedralFourier::read_restart(FILE *fp)
     for (int i=1; i<=atom->ndihedraltypes; i++) {
       fread(k[i],sizeof(double),nterms[i],fp);
       fread(multiplicity[i],sizeof(int),nterms[i],fp);
-      fread(shift[i],sizeof(int),nterms[i],fp);
+      fread(shift[i],sizeof(double),nterms[i],fp);
     }
   }
 
   for (int i=1; i<=atom->ndihedraltypes; i++) {
     MPI_Bcast(k[i],nterms[i],MPI_DOUBLE,0,world);
     MPI_Bcast(multiplicity[i],nterms[i],MPI_INT,0,world);
-    MPI_Bcast(shift[i],nterms[i],MPI_INT,0,world);
+    MPI_Bcast(shift[i],nterms[i],MPI_DOUBLE,0,world);
   }
 
   for (int i=1; i <= atom->ndihedraltypes; i++) {

src/USER-MISC/dihedral_fourier.h

@@ -35,8 +35,8 @@ class DihedralFourier : public Dihedral {
   void read_restart(FILE *);
 
  protected:
-  double **k,**cos_shift,**sin_shift;
+  double **k,**cos_shift,**sin_shift,**shift;
-  int **multiplicity,**shift;
+  int **multiplicity;
   int *nterms;
   int implicit,weightflag;
 

src/USER-MISC/dihedral_nharmonic.cpp

@@ -155,19 +155,19 @@ void DihedralNHarmonic::compute(int eflag, int vflag)
       int me;
       MPI_Comm_rank(world,&me);
       if (screen) {
-	char str[128];
+        char str[128];
-	sprintf(str,"Dihedral problem: %d " BIGINT_FORMAT " %d %d %d %d",
+        sprintf(str,"Dihedral problem: %d " BIGINT_FORMAT " %d %d %d %d",
-		me,update->ntimestep,
+                me,update->ntimestep,
-		atom->tag[i1],atom->tag[i2],atom->tag[i3],atom->tag[i4]);
+                atom->tag[i1],atom->tag[i2],atom->tag[i3],atom->tag[i4]);
-	error->warning(FLERR,str);
+        error->warning(FLERR,str,0);
-	fprintf(screen,"  1st atom: %d %g %g %g\n",
+        fprintf(screen,"  1st atom: %d %g %g %g\n",
-		me,x[i1][0],x[i1][1],x[i1][2]);
+                me,x[i1][0],x[i1][1],x[i1][2]);
-	fprintf(screen,"  2nd atom: %d %g %g %g\n",
+        fprintf(screen,"  2nd atom: %d %g %g %g\n",
-		me,x[i2][0],x[i2][1],x[i2][2]);
+                me,x[i2][0],x[i2][1],x[i2][2]);
-	fprintf(screen,"  3rd atom: %d %g %g %g\n",
+        fprintf(screen,"  3rd atom: %d %g %g %g\n",
-		me,x[i3][0],x[i3][1],x[i3][2]);
+                me,x[i3][0],x[i3][1],x[i3][2]);
-	fprintf(screen,"  4th atom: %d %g %g %g\n",
+        fprintf(screen,"  4th atom: %d %g %g %g\n",
-		me,x[i4][0],x[i4][1],x[i4][2]);
+                me,x[i4][0],x[i4][1],x[i4][2]);
       }
     }
 
@@ -246,7 +246,7 @@ void DihedralNHarmonic::compute(int eflag, int vflag)
 
     if (evflag)
       ev_tally(i1,i2,i3,i4,nlocal,newton_bond,edihedral,f1,f3,f4,
-	       vb1x,vb1y,vb1z,vb2x,vb2y,vb2z,vb3x,vb3y,vb3z);
+               vb1x,vb1y,vb1z,vb2x,vb2y,vb2z,vb3x,vb3y,vb3z);
   }
 }
 

src/USER-MISC/improper_fourier.cpp

@@ -78,21 +78,18 @@ void ImproperFourier::compute(int eflag, int vflag)
     vb1x = x[i2][0] - x[i1][0];
     vb1y = x[i2][1] - x[i1][1];
     vb1z = x[i2][2] - x[i1][2];
-    domain->minimum_image(vb1x,vb1y,vb1z);
 
     // 2nd bond
 
     vb2x = x[i3][0] - x[i1][0];
     vb2y = x[i3][1] - x[i1][1];
     vb2z = x[i3][2] - x[i1][2];
-    domain->minimum_image(vb2x,vb2y,vb2z);
 
     // 3rd bond
 
     vb3x = x[i4][0] - x[i1][0];
     vb3y = x[i4][1] - x[i1][1];
     vb3z = x[i4][2] - x[i1][2];
-    domain->minimum_image(vb3x,vb3y,vb3z);
 
     addone(i1,i2,i3,i4, type,evflag,eflag,
            vb1x, vb1y, vb1z, 
@@ -165,10 +162,14 @@ void ImproperFourier::addone(const int &i1,const int &i2,const int &i3,const int
               me,update->ntimestep,
               atom->tag[i1],atom->tag[i2],atom->tag[i3],atom->tag[i4]);
       error->warning(FLERR,str,0);
-      fprintf(screen,"  1st atom: %d %g %g %g\n", me,x[i1][0],x[i1][1],x[i1][2]);
+      fprintf(screen,"  1st atom: %d %g %g %g\n",
-      fprintf(screen,"  2nd atom: %d %g %g %g\n", me,x[i2][0],x[i2][1],x[i2][2]);
+              me,x[i1][0],x[i1][1],x[i1][2]);
-      fprintf(screen,"  3rd atom: %d %g %g %g\n", me,x[i3][0],x[i3][1],x[i3][2]);
+      fprintf(screen,"  2nd atom: %d %g %g %g\n",
-      fprintf(screen,"  4th atom: %d %g %g %g\n", me,x[i4][0],x[i4][1],x[i4][2]);
+              me,x[i2][0],x[i2][1],x[i2][2]);
+      fprintf(screen,"  3rd atom: %d %g %g %g\n",
+              me,x[i3][0],x[i3][1],x[i3][2]);
+      fprintf(screen,"  4th atom: %d %g %g %g\n",
+              me,x[i4][0],x[i4][1],x[i4][2]);
     }
   }
 

src/compute_reduce.cpp

@@ -156,7 +156,7 @@ ComputeReduce::ComputeReduce(LAMMPS *lmp, int narg, char **arg) :
       int col2 = atoi(arg[iarg+2]) - 1;
       if (col1 < 0 || col1 >= nvalues || col2 < 0 || col2 >= nvalues)
         error->all(FLERR,"Illegal compute reduce command");
-      if (col1 == col2)        error->all(FLERR,"Illegal compute reduce command");
+      if (col1 == col2) error->all(FLERR,"Illegal compute reduce command");
       if (replace[col1] >= 0 || replace[col2] >= 0)
         error->all(FLERR,"Invalid replace values in compute reduce");
       replace[col1] = col2;

src/neigh_request.cpp

@@ -37,7 +37,6 @@ NeighRequest::NeighRequest(LAMMPS *lmp) : Pointers(lmp)
   gran = granhistory = 0;
   respainner = respamiddle = respaouter = 0;
   half_from_full = 0;
-  ghost = 0;
 
   // default is every reneighboring
   // default is use newton_pair setting in force

src/output.cpp

@@ -760,7 +760,7 @@ void Output::memory_usage()
   bytes += update->memory_usage();
   bytes += force->memory_usage();
   bytes += modify->memory_usage();
-  for (int i = 0; i < ndump; i++) dump[i]->memory_usage();
+  for (int i = 0; i < ndump; i++) bytes += dump[i]->memory_usage();
 
   double mbytes = bytes/1024.0/1024.0;
 

src/pair_mie_cut.cpp

@@ -0,0 +1,743 @@
+/* ----------------------------------------------------------------------
+   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.
+------------------------------------------------------------------------- */
+
+/* ----------------------------------------------------------------------
+   Contributing author: Cassiano Aimoli (aimoli@gmail.com)
+------------------------------------------------------------------------- */
+
+#include "math.h"
+#include "stdio.h"
+#include "stdlib.h"
+#include "string.h"
+#include "pair_mie_cut.h"
+#include "atom.h"
+#include "comm.h"
+#include "force.h"
+#include "neighbor.h"
+#include "neigh_list.h"
+#include "neigh_request.h"
+#include "update.h"
+#include "integrate.h"
+#include "respa.h"
+#include "math_const.h"
+#include "memory.h"
+#include "error.h"
+
+using namespace LAMMPS_NS;
+using namespace MathConst;
+
+/* ---------------------------------------------------------------------- */
+
+PairMIECut::PairMIECut(LAMMPS *lmp) : Pair(lmp)
+{
+  respa_enable = 1;
+}
+
+/* ---------------------------------------------------------------------- */
+
+PairMIECut::~PairMIECut()
+{
+  if (allocated) {
+    memory->destroy(setflag);
+    memory->destroy(cutsq);
+
+    memory->destroy(cut);
+    memory->destroy(epsilon);
+    memory->destroy(sigma);
+    memory->destroy(gamR);
+    memory->destroy(gamA);
+    memory->destroy(Cmie);
+    memory->destroy(mie1);
+    memory->destroy(mie2);
+    memory->destroy(mie3);
+    memory->destroy(mie4);
+    memory->destroy(offset);
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void PairMIECut::compute(int eflag, int vflag)
+{
+  int i,j,ii,jj,inum,jnum,itype,jtype;
+  double xtmp,ytmp,ztmp,delx,dely,delz,evdwl,fpair;
+  double rsq,r2inv,rgamR,rgamA,forcemie,factor_mie;
+  int *ilist,*jlist,*numneigh,**firstneigh;
+
+  evdwl = 0.0;
+  if (eflag || vflag) ev_setup(eflag,vflag);
+  else evflag = vflag_fdotr = 0;
+
+  double **x = atom->x;
+  double **f = atom->f;
+  int *type = atom->type;
+  int nlocal = atom->nlocal;
+  double *special_mie = force->special_lj;
+  int newton_pair = force->newton_pair;
+
+  inum = list->inum;
+  ilist = list->ilist;
+  numneigh = list->numneigh;
+  firstneigh = list->firstneigh;
+
+  // loop over neighbors of my atoms
+
+  for (ii = 0; ii < inum; ii++) {
+    i = ilist[ii];
+    xtmp = x[i][0];
+    ytmp = x[i][1];
+    ztmp = x[i][2];
+    itype = type[i];
+    jlist = firstneigh[i];
+    jnum = numneigh[i];
+
+    for (jj = 0; jj < jnum; jj++) {
+      j = jlist[jj];
+      factor_mie = special_mie[sbmask(j)];
+      j &= NEIGHMASK;
+
+      delx = xtmp - x[j][0];
+      dely = ytmp - x[j][1];
+      delz = ztmp - x[j][2];
+      rsq = delx*delx + dely*dely + delz*delz;
+      jtype = type[j];
+
+      if (rsq < cutsq[itype][jtype]) {
+        r2inv = 1.0/rsq;
+        rgamA = pow(r2inv,(gamA[itype][jtype]/2.0));
+        rgamR = pow(r2inv,(gamR[itype][jtype]/2.0));
+        forcemie =  (mie1[itype][jtype]*rgamR - mie2[itype][jtype]*rgamA);
+        fpair = factor_mie*forcemie*r2inv;
+
+        f[i][0] += delx*fpair;
+        f[i][1] += dely*fpair;
+        f[i][2] += delz*fpair;
+        if (newton_pair || j < nlocal) {
+          f[j][0] -= delx*fpair;
+          f[j][1] -= dely*fpair;
+          f[j][2] -= delz*fpair;
+        }
+
+        if (eflag) {
+          evdwl = (mie3[itype][jtype]*rgamR - mie4[itype][jtype]*rgamA) -
+            offset[itype][jtype];
+          evdwl *= factor_mie;
+        }
+
+        if (evflag) ev_tally(i,j,nlocal,newton_pair,
+                             evdwl,0.0,fpair,delx,dely,delz);
+      }
+    }
+  }
+
+  if (vflag_fdotr) virial_fdotr_compute();
+}
+
+/* ---------------------------------------------------------------------- */
+
+void PairMIECut::compute_inner()
+{
+  int i,j,ii,jj,inum,jnum,itype,jtype;
+  double xtmp,ytmp,ztmp,delx,dely,delz,fpair;
+  double rsq,r2inv,rgamR,rgamA,forcemie,factor_mie,rsw;
+  int *ilist,*jlist,*numneigh,**firstneigh;
+
+  double **x = atom->x;
+  double **f = atom->f;
+  int *type = atom->type;
+  int nlocal = atom->nlocal;
+  double *special_mie = force->special_lj;
+  int newton_pair = force->newton_pair;
+
+  inum = listinner->inum;
+  ilist = listinner->ilist;
+  numneigh = listinner->numneigh;
+  firstneigh = listinner->firstneigh;
+
+  double cut_out_on = cut_respa[0];
+  double cut_out_off = cut_respa[1];
+
+  double cut_out_diff = cut_out_off - cut_out_on;
+  double cut_out_on_sq = cut_out_on*cut_out_on;
+  double cut_out_off_sq = cut_out_off*cut_out_off;
+
+  // loop over neighbors of my atoms
+
+  for (ii = 0; ii < inum; ii++) {
+    i = ilist[ii];
+    xtmp = x[i][0];
+    ytmp = x[i][1];
+    ztmp = x[i][2];
+    itype = type[i];
+    jlist = firstneigh[i];
+    jnum = numneigh[i];
+
+    for (jj = 0; jj < jnum; jj++) {
+      j = jlist[jj];
+      factor_mie = special_mie[sbmask(j)];
+      j &= NEIGHMASK;
+
+      delx = xtmp - x[j][0];
+      dely = ytmp - x[j][1];
+      delz = ztmp - x[j][2];
+      rsq = delx*delx + dely*dely + delz*delz;
+
+      if (rsq < cut_out_off_sq) {
+        r2inv = 1.0/rsq;
+        rgamA = pow(r2inv,(gamA[itype][jtype]/2.0));
+        rgamR = pow(r2inv,(gamR[itype][jtype]/2.0));
+        jtype = type[j];
+        forcemie =  (mie1[itype][jtype]*rgamR - mie2[itype][jtype]*rgamA);
+        fpair = factor_mie*forcemie*r2inv;
+        if (rsq > cut_out_on_sq) {
+          rsw = (sqrt(rsq) - cut_out_on)/cut_out_diff;
+          fpair *= 1.0 - rsw*rsw*(3.0 - 2.0*rsw);
+        }
+
+        f[i][0] += delx*fpair;
+        f[i][1] += dely*fpair;
+        f[i][2] += delz*fpair;
+        if (newton_pair || j < nlocal) {
+          f[j][0] -= delx*fpair;
+          f[j][1] -= dely*fpair;
+          f[j][2] -= delz*fpair;
+        }
+      }
+    }
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void PairMIECut::compute_middle()
+{
+  int i,j,ii,jj,inum,jnum,itype,jtype;
+  double xtmp,ytmp,ztmp,delx,dely,delz,fpair;
+  double rsq,r2inv,rgamR,rgamA,forcemie,factor_mie,rsw;
+  int *ilist,*jlist,*numneigh,**firstneigh;
+
+  double **x = atom->x;
+  double **f = atom->f;
+  int *type = atom->type;
+  int nlocal = atom->nlocal;
+  double *special_mie = force->special_lj;
+  int newton_pair = force->newton_pair;
+
+  inum = listmiddle->inum;
+  ilist = listmiddle->ilist;
+  numneigh = listmiddle->numneigh;
+  firstneigh = listmiddle->firstneigh;
+
+  double cut_in_off = cut_respa[0];
+  double cut_in_on = cut_respa[1];
+  double cut_out_on = cut_respa[2];
+  double cut_out_off = cut_respa[3];
+
+  double cut_in_diff = cut_in_on - cut_in_off;
+  double cut_out_diff = cut_out_off - cut_out_on;
+  double cut_in_off_sq = cut_in_off*cut_in_off;
+  double cut_in_on_sq = cut_in_on*cut_in_on;
+  double cut_out_on_sq = cut_out_on*cut_out_on;
+  double cut_out_off_sq = cut_out_off*cut_out_off;
+
+  // loop over neighbors of my atoms
+
+  for (ii = 0; ii < inum; ii++) {
+    i = ilist[ii];
+    xtmp = x[i][0];
+    ytmp = x[i][1];
+    ztmp = x[i][2];
+    itype = type[i];
+    jlist = firstneigh[i];
+    jnum = numneigh[i];
+
+    for (jj = 0; jj < jnum; jj++) {
+      j = jlist[jj];
+      factor_mie = special_mie[sbmask(j)];
+      j &= NEIGHMASK;
+
+      delx = xtmp - x[j][0];
+      dely = ytmp - x[j][1];
+      delz = ztmp - x[j][2];
+      rsq = delx*delx + dely*dely + delz*delz;
+
+      if (rsq < cut_out_off_sq && rsq > cut_in_off_sq) {
+        r2inv = 1.0/rsq;
+        rgamA = pow(r2inv,(gamA[itype][jtype]/2.0));
+        rgamR = pow(r2inv,(gamR[itype][jtype]/2.0));
+        jtype = type[j];
+        forcemie =  (mie1[itype][jtype]*rgamR - mie2[itype][jtype]*rgamA);
+        fpair = factor_mie*forcemie*r2inv;
+        if (rsq < cut_in_on_sq) {
+          rsw = (sqrt(rsq) - cut_in_off)/cut_in_diff;
+          fpair *= rsw*rsw*(3.0 - 2.0*rsw);
+        }
+        if (rsq > cut_out_on_sq) {
+          rsw = (sqrt(rsq) - cut_out_on)/cut_out_diff;
+          fpair *= 1.0 + rsw*rsw*(2.0*rsw - 3.0);
+        }
+
+        f[i][0] += delx*fpair;
+        f[i][1] += dely*fpair;
+        f[i][2] += delz*fpair;
+        if (newton_pair || j < nlocal) {
+          f[j][0] -= delx*fpair;
+          f[j][1] -= dely*fpair;
+          f[j][2] -= delz*fpair;
+        }
+      }
+    }
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void PairMIECut::compute_outer(int eflag, int vflag)
+{
+  int i,j,ii,jj,inum,jnum,itype,jtype;
+  double xtmp,ytmp,ztmp,delx,dely,delz,evdwl,fpair;
+  double rsq,r2inv,rgamR,rgamA,forcemie,factor_mie,rsw;
+  int *ilist,*jlist,*numneigh,**firstneigh;
+
+  evdwl = 0.0;
+  if (eflag || vflag) ev_setup(eflag,vflag);
+  else evflag = 0;
+
+  double **x = atom->x;
+  double **f = atom->f;
+  int *type = atom->type;
+  int nlocal = atom->nlocal;
+  double *special_mie = force->special_lj;
+  int newton_pair = force->newton_pair;
+
+  inum = listouter->inum;
+  ilist = listouter->ilist;
+  numneigh = listouter->numneigh;
+  firstneigh = listouter->firstneigh;
+
+  double cut_in_off = cut_respa[2];
+  double cut_in_on = cut_respa[3];
+
+  double cut_in_diff = cut_in_on - cut_in_off;
+  double cut_in_off_sq = cut_in_off*cut_in_off;
+  double cut_in_on_sq = cut_in_on*cut_in_on;
+
+  // loop over neighbors of my atoms
+
+  for (ii = 0; ii < inum; ii++) {
+    i = ilist[ii];
+    xtmp = x[i][0];
+    ytmp = x[i][1];
+    ztmp = x[i][2];
+    itype = type[i];
+    jlist = firstneigh[i];
+    jnum = numneigh[i];
+
+    for (jj = 0; jj < jnum; jj++) {
+      j = jlist[jj];
+      factor_mie = special_mie[sbmask(j)];
+      j &= NEIGHMASK;
+
+      delx = xtmp - x[j][0];
+      dely = ytmp - x[j][1];
+      delz = ztmp - x[j][2];
+      rsq = delx*delx + dely*dely + delz*delz;
+      jtype = type[j];
+
+      if (rsq < cutsq[itype][jtype]) {
+        if (rsq > cut_in_off_sq) {
+        r2inv = 1.0/rsq;
+        rgamA = pow(r2inv,(gamA[itype][jtype]/2.0));
+        rgamR = pow(r2inv,(gamR[itype][jtype]/2.0));
+        forcemie =  (mie1[itype][jtype]*rgamR - mie2[itype][jtype]*rgamA);
+          fpair = factor_mie*forcemie*r2inv;
+          if (rsq < cut_in_on_sq) {
+            rsw = (sqrt(rsq) - cut_in_off)/cut_in_diff;
+            fpair *= rsw*rsw*(3.0 - 2.0*rsw);
+          }
+
+          f[i][0] += delx*fpair;
+          f[i][1] += dely*fpair;
+          f[i][2] += delz*fpair;
+          if (newton_pair || j < nlocal) {
+            f[j][0] -= delx*fpair;
+            f[j][1] -= dely*fpair;
+            f[j][2] -= delz*fpair;
+          }
+        }
+
+        if (eflag) {
+          r2inv = 1.0/rsq;
+        rgamA = pow(r2inv,(gamA[itype][jtype]/2.0));
+        rgamR = pow(r2inv,(gamR[itype][jtype]/2.0));
+          evdwl = (mie3[itype][jtype]*rgamR - mie4[itype][jtype]*rgamA) - 
+            offset[itype][jtype];
+          evdwl *= factor_mie;
+        }
+
+        if (vflag) {
+          if (rsq <= cut_in_off_sq) {
+        r2inv = 1.0/rsq;
+        rgamA = pow(r2inv,(gamA[itype][jtype]/2.0));
+        rgamR = pow(r2inv,(gamR[itype][jtype]/2.0));
+        forcemie =  (mie1[itype][jtype]*rgamR - mie2[itype][jtype]*rgamA);
+            fpair = factor_mie*forcemie*r2inv;
+          } else if (rsq < cut_in_on_sq)
+            fpair = factor_mie*forcemie*r2inv;
+        }
+
+        if (evflag) ev_tally(i,j,nlocal,newton_pair,
+                             evdwl,0.0,fpair,delx,dely,delz);
+      }
+    }
+  }
+}
+
+/* ----------------------------------------------------------------------
+   allocate all arrays
+------------------------------------------------------------------------- */
+
+void PairMIECut::allocate()
+{
+  allocated = 1;
+  int n = atom->ntypes;
+
+  memory->create(setflag,n+1,n+1,"pair:setflag");
+  for (int i = 1; i <= n; i++)
+    for (int j = i; j <= n; j++)
+      setflag[i][j] = 0;
+
+  memory->create(cutsq,n+1,n+1,"pair:cutsq");
+
+  memory->create(cut,n+1,n+1,"pair:cut");
+  memory->create(epsilon,n+1,n+1,"pair:epsilon");
+  memory->create(sigma,n+1,n+1,"pair:sigma");
+  memory->create(gamR,n+1,n+1,"pair:gamR");
+  memory->create(gamA,n+1,n+1,"pair:gamA");
+  memory->create(Cmie,n+1,n+1,"pair:Cmie");
+  memory->create(mie1,n+1,n+1,"pair:mie1");
+  memory->create(mie2,n+1,n+1,"pair:mie2");
+  memory->create(mie3,n+1,n+1,"pair:mie3");
+  memory->create(mie4,n+1,n+1,"pair:mie4");
+  memory->create(offset,n+1,n+1,"pair:offset");
+}
+
+/* ----------------------------------------------------------------------
+   global settings
+------------------------------------------------------------------------- */
+
+void PairMIECut::settings(int narg, char **arg)
+{
+  if (narg != 1) error->all(FLERR,"Illegal pair_style command");
+
+  cut_global = force->numeric(arg[0]);
+
+  // reset cutoffs that have been explicitly set
+
+  if (allocated) {
+    int i,j;
+    for (i = 1; i <= atom->ntypes; i++)
+      for (j = i+1; j <= atom->ntypes; j++)
+        if (setflag[i][j]) cut[i][j] = cut_global;
+  }
+}
+
+/* ----------------------------------------------------------------------
+   set coeffs for one or more type pairs
+------------------------------------------------------------------------- */
+
+void PairMIECut::coeff(int narg, char **arg)
+{
+  if (narg < 6 || narg > 7)
+    error->all(FLERR,"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 epsilon_one = force->numeric(arg[2]);
+  double sigma_one = force->numeric(arg[3]);
+  double gamR_one = force->numeric(arg[4]);
+  double gamA_one = force->numeric(arg[5]);
+
+  double cut_one = cut_global;
+  if (narg == 7) cut_one = force->numeric(arg[6]);
+
+  int count = 0;
+  for (int i = ilo; i <= ihi; i++) {
+    for (int j = MAX(jlo,i); j <= jhi; j++) {
+      epsilon[i][j] = epsilon_one;
+      sigma[i][j] = sigma_one;
+      gamR[i][j] = gamR_one;
+      gamA[i][j] = gamA_one;
+      cut[i][j] = cut_one;
+      setflag[i][j] = 1;
+      count++;
+    }
+  }
+
+  if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients");
+}
+
+/* ----------------------------------------------------------------------
+   init specific to this pair style
+------------------------------------------------------------------------- */
+
+void PairMIECut::init_style()
+{
+  // request regular or rRESPA neighbor lists
+
+  int irequest;
+
+  if (update->whichflag == 1 && strstr(update->integrate_style,"respa")) {
+    int respa = 0;
+    if (((Respa *) update->integrate)->level_inner >= 0) respa = 1;
+    if (((Respa *) update->integrate)->level_middle >= 0) respa = 2;
+
+    if (respa == 0) irequest = neighbor->request(this);
+    else if (respa == 1) {
+      irequest = neighbor->request(this);
+      neighbor->requests[irequest]->id = 1;
+      neighbor->requests[irequest]->half = 0;
+      neighbor->requests[irequest]->respainner = 1;
+      irequest = neighbor->request(this);
+      neighbor->requests[irequest]->id = 3;
+      neighbor->requests[irequest]->half = 0;
+      neighbor->requests[irequest]->respaouter = 1;
+    } else {
+      irequest = neighbor->request(this);
+      neighbor->requests[irequest]->id = 1;
+      neighbor->requests[irequest]->half = 0;
+      neighbor->requests[irequest]->respainner = 1;
+      irequest = neighbor->request(this);
+      neighbor->requests[irequest]->id = 2;
+      neighbor->requests[irequest]->half = 0;
+      neighbor->requests[irequest]->respamiddle = 1;
+      irequest = neighbor->request(this);
+      neighbor->requests[irequest]->id = 3;
+      neighbor->requests[irequest]->half = 0;
+      neighbor->requests[irequest]->respaouter = 1;
+    }
+
+  } else irequest = neighbor->request(this);
+
+  // set rRESPA cutoffs
+
+  if (strstr(update->integrate_style,"respa") &&
+      ((Respa *) update->integrate)->level_inner >= 0)
+    cut_respa = ((Respa *) update->integrate)->cutoff;
+  else cut_respa = NULL;
+}
+
+/* ----------------------------------------------------------------------
+   neighbor callback to inform pair style of neighbor list to use
+   regular or rRESPA
+------------------------------------------------------------------------- */
+
+void PairMIECut::init_list(int id, NeighList *ptr)
+{
+  if (id == 0) list = ptr;
+  else if (id == 1) listinner = ptr;
+  else if (id == 2) listmiddle = ptr;
+  else if (id == 3) listouter = ptr;
+}
+
+/* ----------------------------------------------------------------------
+   init for one type pair i,j and corresponding j,i
+------------------------------------------------------------------------- */
+
+double PairMIECut::init_one(int i, int j)
+{
+  if (setflag[i][j] == 0) {
+    epsilon[i][j] = mix_energy(epsilon[i][i],epsilon[j][j],
+                               sigma[i][i],sigma[j][j]);
+    sigma[i][j] = mix_distance(sigma[i][i],sigma[j][j]);
+    gamR[i][j] = mix_distance(gamR[i][i],gamR[j][j]);
+    gamA[i][j] = mix_distance(gamA[i][i],gamA[j][j]);
+    cut[i][j] = mix_distance(cut[i][i],cut[j][j]);
+  }
+  
+  gamA[j][i] = gamA[i][j];
+  gamR[j][i] = gamR[i][j];
+  Cmie[i][j] = (gamR[i][j]/(gamR[i][j]-gamA[i][j]) * 
+                pow((gamR[i][j]/gamA[i][j]),
+                    (gamA[i][j]/(gamR[i][j]-gamA[i][j]))));
+  mie1[i][j] = Cmie[i][j] * gamR[i][j]* epsilon[i][j] * 
+    pow(sigma[i][j],gamR[i][j]);
+  mie2[i][j] = Cmie[i][j] * gamA[i][j] * epsilon[i][j] * 
+    pow(sigma[i][j],gamA[i][j]);
+  mie3[i][j] = Cmie[i][j] * epsilon[i][j] * pow(sigma[i][j],gamR[i][j]);
+  mie4[i][j] = Cmie[i][j] * epsilon[i][j] * pow(sigma[i][j],gamA[i][j]);
+
+  if (offset_flag) {
+    double ratio = sigma[i][j] / cut[i][j];
+    offset[i][j] = Cmie[i][j] * epsilon[i][j] * 
+      (pow(ratio,gamR[i][j]) - pow(ratio,gamA[i][j]));
+  } else offset[i][j] = 0.0;
+
+  mie1[j][i] = mie1[i][j];
+  mie2[j][i] = mie2[i][j];
+  mie3[j][i] = mie3[i][j];
+  mie4[j][i] = mie4[i][j];
+  offset[j][i] = offset[i][j];
+
+  // check interior rRESPA cutoff
+
+  if (cut_respa && cut[i][j] < cut_respa[3])
+    error->all(FLERR,"Pair cutoff < Respa interior cutoff");
+
+  // compute I,J contribution to long-range tail correction
+  // count total # of atoms of type I and J via Allreduce
+
+  if (tail_flag) {
+    int *type = atom->type;
+    int nlocal = atom->nlocal;
+
+    double count[2],all[2];
+    count[0] = count[1] = 0.0;
+    for (int k = 0; k < nlocal; k++) {
+      if (type[k] == i) count[0] += 1.0;
+      if (type[k] == j) count[1] += 1.0;
+    }
+    MPI_Allreduce(count,all,2,MPI_DOUBLE,MPI_SUM,world);
+
+    double siggamA = pow(sigma[i][j],gamA[i][j]);
+    double siggamR = pow(sigma[i][j],gamR[i][j]);
+    double rcgamA = pow(cut[i][j],(gamA[i][j]-3.0));
+    double rcgamR = pow(cut[i][j],(gamR[i][j]-3.0));
+    etail_ij = Cmie[i][j]*2.0*MY_PI*all[0]*all[1]*epsilon[i][j]*
+      (siggamR/((gamR[i][j]-3.0)*rcgamR)-siggamA/((gamA[i][j]-3.0)*rcgamA));
+    ptail_ij = Cmie[i][j]*2.0*MY_PI*all[0]*all[1]*epsilon[i][j]/3.0*
+      ((gamR[i][j]/(gamR[i][j]-3.0))*siggamR/rcgamR-
+       (gamA[i][j]/(gamA[i][j]-3.0))*siggamA/rcgamA);
+  }
+
+  return cut[i][j];
+}
+
+/* ----------------------------------------------------------------------
+   proc 0 writes to restart file
+------------------------------------------------------------------------- */
+
+void PairMIECut::write_restart(FILE *fp)
+{
+  write_restart_settings(fp);
+
+  int i,j;
+  for (i = 1; i <= atom->ntypes; i++)
+    for (j = i; j <= atom->ntypes; j++) {
+      fwrite(&setflag[i][j],sizeof(int),1,fp);
+      if (setflag[i][j]) {
+        fwrite(&epsilon[i][j],sizeof(double),1,fp);
+        fwrite(&sigma[i][j],sizeof(double),1,fp);
+        fwrite(&gamR[i][j],sizeof(double),1,fp);
+        fwrite(&gamA[i][j],sizeof(double),1,fp);
+        fwrite(&cut[i][j],sizeof(double),1,fp);
+      }
+    }
+}
+
+/* ----------------------------------------------------------------------
+   proc 0 reads from restart file, bcasts
+------------------------------------------------------------------------- */
+
+void PairMIECut::read_restart(FILE *fp)
+{
+  read_restart_settings(fp);
+  allocate();
+
+  int i,j;
+  int me = comm->me;
+  for (i = 1; i <= atom->ntypes; i++)
+    for (j = i; j <= atom->ntypes; j++) {
+      if (me == 0) fread(&setflag[i][j],sizeof(int),1,fp);
+      MPI_Bcast(&setflag[i][j],1,MPI_INT,0,world);
+      if (setflag[i][j]) {
+        if (me == 0) {
+          fread(&epsilon[i][j],sizeof(double),1,fp);
+          fread(&sigma[i][j],sizeof(double),1,fp);
+          fread(&gamR[i][j],sizeof(double),1,fp);
+          fread(&gamA[i][j],sizeof(double),1,fp);
+          fread(&cut[i][j],sizeof(double),1,fp);
+        }
+        MPI_Bcast(&epsilon[i][j],1,MPI_DOUBLE,0,world);
+        MPI_Bcast(&sigma[i][j],1,MPI_DOUBLE,0,world);
+        MPI_Bcast(&gamR[i][j],1,MPI_DOUBLE,0,world);
+        MPI_Bcast(&gamA[i][j],1,MPI_DOUBLE,0,world);
+        MPI_Bcast(&cut[i][j],1,MPI_DOUBLE,0,world);
+      }
+    }
+}
+
+/* ----------------------------------------------------------------------
+   proc 0 writes to restart file
+------------------------------------------------------------------------- */
+
+void PairMIECut::write_restart_settings(FILE *fp)
+{
+  fwrite(&cut_global,sizeof(double),1,fp);
+  fwrite(&offset_flag,sizeof(int),1,fp);
+  fwrite(&mix_flag,sizeof(int),1,fp);
+  fwrite(&tail_flag,sizeof(int),1,fp);
+}
+
+/* ----------------------------------------------------------------------
+   proc 0 reads from restart file, bcasts
+------------------------------------------------------------------------- */
+
+void PairMIECut::read_restart_settings(FILE *fp)
+{
+  int me = comm->me;
+  if (me == 0) {
+    fread(&cut_global,sizeof(double),1,fp);
+    fread(&offset_flag,sizeof(int),1,fp);
+    fread(&mix_flag,sizeof(int),1,fp);
+    fread(&tail_flag,sizeof(int),1,fp);
+  }
+  MPI_Bcast(&cut_global,1,MPI_DOUBLE,0,world);
+  MPI_Bcast(&offset_flag,1,MPI_INT,0,world);
+  MPI_Bcast(&mix_flag,1,MPI_INT,0,world);
+  MPI_Bcast(&tail_flag,1,MPI_INT,0,world);
+}
+
+/* ---------------------------------------------------------------------- */
+
+double PairMIECut::single(int i, int j, int itype, int jtype, double rsq,
+                           double factor_coul, double factor_mie,
+                           double &fforce)
+{
+  double r2inv,rgamR,rgamA,forcemie,phimie;
+
+  r2inv = 1.0/rsq;
+  rgamA = pow(r2inv,(gamA[itype][jtype]/2.0));
+  rgamR = pow(r2inv,(gamR[itype][jtype]/2.0));
+  forcemie =  (mie1[itype][jtype]*rgamR - mie2[itype][jtype]*rgamA);
+  fforce = factor_mie*forcemie*r2inv;
+
+  phimie = (mie3[itype][jtype]*rgamR - mie4[itype][jtype]*rgamA) -
+    offset[itype][jtype];
+  return factor_mie*phimie;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void *PairMIECut::extract(const char *str, int &dim)
+{
+  dim = 2;
+  if (strcmp(str,"epsilon") == 0) return (void *) epsilon;
+  if (strcmp(str,"sigma") == 0) return (void *) sigma;
+  if (strcmp(str,"gamR") == 0) return (void *) gamR;
+  if (strcmp(str,"gamA") == 0) return (void *) gamA;
+  return NULL;
+}

src/pair_mie_cut.h

@@ -0,0 +1,81 @@
+/* ----------------------------------------------------------------------
+   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 PAIR_CLASS
+
+PairStyle(mie/cut,PairMIECut)
+
+#else
+
+#ifndef LMP_PAIR_MIE_CUT_H
+#define LMP_PAIR_MIE_CUT_H
+
+#include "pair.h"
+
+namespace LAMMPS_NS {
+
+class PairMIECut : public Pair {
+ public:
+  PairMIECut(class LAMMPS *);
+  virtual ~PairMIECut();
+  virtual void compute(int, int);
+  void settings(int, char **);
+  void coeff(int, char **);
+  void init_style();
+  void init_list(int, class NeighList *);
+  double init_one(int, int);
+  void write_restart(FILE *);
+  void read_restart(FILE *);
+  void write_restart_settings(FILE *);
+  void read_restart_settings(FILE *);
+  double single(int, int, int, int, double, double, double, double &);
+  void *extract(const char *, int &);
+
+  void compute_inner();
+  void compute_middle();
+  void compute_outer(int, int);
+
+ protected:
+  double cut_global;
+  double **cut;
+  double **epsilon,**sigma;
+  double **gamR,**gamA,**Cmie;
+  double **mie1,**mie2,**mie3,**mie4,**offset;
+  double *cut_respa;
+
+  void allocate();
+};
+
+}
+
+#endif
+#endif
+
+/* ERROR/WARNING messages:
+
+E: Illegal ... command
+
+Self-explanatory.  Check the input script syntax and compare to the
+documentation for the command.  You can use -echo screen as a
+command-line option when running LAMMPS to see the offending line.
+
+E: Incorrect args for pair coefficients
+
+Self-explanatory.  Check the input script or data file.
+
+E: Pair cutoff < Respa interior cutoff
+
+One or more pairwise cutoffs are too short to use with the specified
+rRESPA cutoffs.
+
+*/

src/read_data.cpp

@@ -552,6 +552,7 @@ void ReadData::atoms()
         if (eof == NULL) error->one(FLERR,"Unexpected end of data file");
         m += strlen(&buffer[m]);
       }
+      if (buffer[m-1] != '\n') strcpy(&buffer[m++],"\n");
       m++;
     }
     MPI_Bcast(&m,1,MPI_INT,0,world);
@@ -644,6 +645,7 @@ void ReadData::velocities()
         if (eof == NULL) error->one(FLERR,"Unexpected end of data file");
         m += strlen(&buffer[m]);
       }
+      if (buffer[m-1] != '\n') strcpy(&buffer[m++],"\n");
       m++;
     }
     MPI_Bcast(&m,1,MPI_INT,0,world);
@@ -695,6 +697,7 @@ void ReadData::bonus(bigint nbonus, AtomVec *ptr, const char *type)
         if (eof == NULL) error->one(FLERR,"Unexpected end of data file");
         m += strlen(&buffer[m]);
       }
+      if (buffer[m-1] != '\n') strcpy(&buffer[m++],"\n");
       m++;
     }
     MPI_Bcast(&m,1,MPI_INT,0,world);
@@ -734,6 +737,7 @@ void ReadData::bonds()
         if (eof == NULL) error->one(FLERR,"Unexpected end of data file");
         m += strlen(&buffer[m]);
       }
+      if (buffer[m-1] != '\n') strcpy(&buffer[m++],"\n");
       m++;
     }
     MPI_Bcast(&m,1,MPI_INT,0,world);
@@ -780,6 +784,7 @@ void ReadData::angles()
         if (eof == NULL) error->one(FLERR,"Unexpected end of data file");
         m += strlen(&buffer[m]);
       }
+      if (buffer[m-1] != '\n') strcpy(&buffer[m++],"\n");
       m++;
     }
     MPI_Bcast(&m,1,MPI_INT,0,world);
@@ -826,6 +831,7 @@ void ReadData::dihedrals()
         if (eof == NULL) error->one(FLERR,"Unexpected end of data file");
         m += strlen(&buffer[m]);
       }
+      if (buffer[m-1] != '\n') strcpy(&buffer[m++],"\n");
       m++;
     }
     MPI_Bcast(&m,1,MPI_INT,0,world);
@@ -872,6 +878,7 @@ void ReadData::impropers()
         if (eof == NULL) error->one(FLERR,"Unexpected end of data file");
         m += strlen(&buffer[m]);
       }
+      if (buffer[m-1] != '\n') strcpy(&buffer[m++],"\n");
       m++;
     }
     MPI_Bcast(&m,1,MPI_INT,0,world);
@@ -914,6 +921,7 @@ void ReadData::mass()
       eof = fgets(&buf[m],MAXLINE,fp);
       if (eof == NULL) error->one(FLERR,"Unexpected end of data file");
       m += strlen(&buf[m]);
+      if (buf[m-1] != '\n') strcpy(&buf[m++],"\n");
       buf[m-1] = '\0';
     }
   }
@@ -943,6 +951,7 @@ void ReadData::paircoeffs()
       eof = fgets(&buf[m],MAXLINE,fp);
       if (eof == NULL) error->one(FLERR,"Unexpected end of data file");
       m += strlen(&buf[m]);
+      if (buf[m-1] != '\n') strcpy(&buf[m++],"\n");
       buf[m-1] = '\0';
     }
   }
@@ -974,6 +983,7 @@ void ReadData::bondcoeffs()
       eof = fgets(&buf[m],MAXLINE,fp);
       if (eof == NULL) error->one(FLERR,"Unexpected end of data file");
       m += strlen(&buf[m]);
+      if (buf[m-1] != '\n') strcpy(&buf[m++],"\n");
       buf[m-1] = '\0';
     }
   }
@@ -1005,6 +1015,7 @@ void ReadData::anglecoeffs(int which)
       eof = fgets(&buf[m],MAXLINE,fp);
       if (eof == NULL) error->one(FLERR,"Unexpected end of data file");
       m += strlen(&buf[m]);
+      if (buf[m-1] != '\n') strcpy(&buf[m++],"\n");
       buf[m-1] = '\0';
     }
   }
@@ -1038,6 +1049,7 @@ void ReadData::dihedralcoeffs(int which)
       eof = fgets(&buf[m],MAXLINE,fp);
       if (eof == NULL) error->one(FLERR,"Unexpected end of data file");
       m += strlen(&buf[m]);
+      if (buf[m-1] != '\n') strcpy(&buf[m++],"\n");
       buf[m-1] = '\0';
     }
   }
@@ -1074,6 +1086,7 @@ void ReadData::impropercoeffs(int which)
       eof = fgets(&buf[m],MAXLINE,fp);
       if (eof == NULL) error->one(FLERR,"Unexpected end of data file");
       m += strlen(&buf[m]);
+      if (buf[m-1] != '\n') strcpy(&buf[m++],"\n");
       buf[m-1] = '\0';
     }
   }
@@ -1113,6 +1126,7 @@ void ReadData::fix(int ifix, char *line, bigint nlines)
         if (eof == NULL) error->one(FLERR,"Unexpected end of data file");
         m += strlen(&buffer[m]);
       }
+      if (buffer[m-1] != '\n') strcpy(&buffer[m++],"\n");
       m++;
     }
     MPI_Bcast(&m,1,MPI_INT,0,world);

src/thermo.cpp

@@ -1490,7 +1490,7 @@ void Thermo::compute_spcpu()
 
 void Thermo::compute_atoms()
 {
-  bivalue = natoms;
+  bivalue = atom->natoms;
 }
 
 /* ---------------------------------------------------------------------- */

src/variable.cpp

@@ -39,10 +39,11 @@ using namespace MathConst;
 
 #define VARDELTA 4
 #define MAXLEVEL 4
+#define MAXLINE 256
 
 #define MYROUND(a) (( a-floor(a) ) >= .5) ? ceil(a) : floor(a)
 
-enum{INDEX,LOOP,WORLD,UNIVERSE,ULOOP,STRING,EQUAL,ATOM};
+enum{INDEX,LOOP,WORLD,UNIVERSE,ULOOP,STRING,FILEVAR,EQUAL,ATOM};
 enum{ARG,OP};
 
 // customize by adding a function
@@ -56,7 +57,7 @@ enum{DONE,ADD,SUBTRACT,MULTIPLY,DIVIDE,CARAT,UNARY,
 
 // customize by adding a special function
 
-enum{SUM,XMIN,XMAX,AVE,TRAP};
+enum{SUM,XMIN,XMAX,AVE,TRAP,NEXT};
 
 #define INVOKED_SCALAR 1
 #define INVOKED_VECTOR 2
@@ -77,6 +78,7 @@ Variable::Variable(LAMMPS *lmp) : Pointers(lmp)
   num = NULL;
   which = NULL;
   pad = NULL;
+  reader = NULL;
   data = NULL;
 
   eval_in_progress = NULL;
@@ -101,6 +103,7 @@ Variable::~Variable()
 {
   for (int i = 0; i < nvar; i++) {
     delete [] names[i];
+    delete reader[i];
     if (style[i] == LOOP || style[i] == ULOOP) delete [] data[i][0];
     else for (int j = 0; j < num[i]; j++) delete [] data[i][j];
     delete [] data[i];
@@ -110,6 +113,7 @@ Variable::~Variable()
   memory->destroy(num);
   memory->destroy(which);
   memory->destroy(pad);
+  memory->sfree(reader);
   memory->sfree(data);
 
   memory->destroy(eval_in_progress);
@@ -245,7 +249,8 @@ void Variable::set(int narg, char **arg)
     for (int jvar = 0; jvar < nvar; jvar++)
       if (num[jvar] && (style[jvar] == UNIVERSE || style[jvar] == ULOOP) &&
           num[nvar] != num[jvar])
-        error->all(FLERR,"All universe/uloop variables must have same # of values");
+        error->all(FLERR,
+                   "All universe/uloop variables must have same # of values");
 
   // STRING
   // remove pre-existing var if also style STRING (allows it to be reset)
@@ -267,6 +272,23 @@ void Variable::set(int narg, char **arg)
     data[nvar] = new char*[num[nvar]];
     copy(1,&arg[2],data[nvar]);
 
+  // FILEVAR for strings or numbers
+  // which = 1st value
+
+  } else if (strcmp(arg[1],"file") == 0) {
+    if (narg != 3) error->all(FLERR,"Illegal variable command");
+    if (find(arg[0]) >= 0) return;
+    if (nvar == maxvar) grow();
+    style[nvar] = FILEVAR;
+    num[nvar] = 1;
+    which[nvar] = 0;
+    pad[nvar] = 0;
+    data[nvar] = new char*[num[nvar]];
+    data[nvar][0] = new char[MAXLINE];
+    reader[nvar] = new VarReader(lmp,arg[2]);
+    int flag = reader[nvar]->read(data[nvar][0]);
+    if (flag) error->all(FLERR,"File variable could not read value");
+
   // EQUAL
   // remove pre-existing var if also style EQUAL (allows it to be reset)
   // num = 2, which = 1st value
@@ -385,6 +407,17 @@ int Variable::next(int narg, char **arg)
       }
     }
 
+  } else if (istyle == FILEVAR) {
+
+    for (int iarg = 0; iarg < narg; iarg++) {
+      ivar = find(arg[iarg]);
+      int done = reader[ivar]->read(data[ivar][0]);
+      if (done) {
+        flag = 1;
+        remove(ivar);
+      }
+    }
+
   } else if (istyle == UNIVERSE || istyle == ULOOP) {
 
     // wait until lock file can be created and owned by proc 0 of this world
@@ -445,7 +478,8 @@ char *Variable::retrieve(char *name)
 
   char *str;
   if (style[ivar] == INDEX || style[ivar] == WORLD ||
-      style[ivar] == UNIVERSE || style[ivar] == STRING) {
+      style[ivar] == UNIVERSE || style[ivar] == STRING || 
+      style[ivar] == FILEVAR) {
     str = data[ivar][which[ivar]];
   } else if (style[ivar] == LOOP || style[ivar] == ULOOP) {
     char result[16];
@@ -575,6 +609,7 @@ void Variable::remove(int n)
     num[i-1] = num[i];
     which[i-1] = which[i];
     pad[i-1] = pad[i];
+    reader[i-1] = reader[i];
     data[i-1] = data[i];
   }
   nvar--;
@@ -586,14 +621,20 @@ void Variable::remove(int n)
 
 void Variable::grow()
 {
+  int old = maxvar;
   maxvar += VARDELTA;
-  names = (char **)
+  names = (char **) memory->srealloc(names,maxvar*sizeof(char *),"var:names");
-  memory->srealloc(names,maxvar*sizeof(char *),"var:names");
   memory->grow(style,maxvar,"var:style");
   memory->grow(num,maxvar,"var:num");
   memory->grow(which,maxvar,"var:which");
   memory->grow(pad,maxvar,"var:pad");
+
+  reader = (VarReader **) 
+    memory->srealloc(reader,maxvar*sizeof(VarReader *),"var:reader");
+  for (int i = old; i < maxvar; i++) reader[i] = NULL;
+
   data = (char ***) memory->srealloc(data,maxvar*sizeof(char **),"var:data");
+
   memory->grow(eval_in_progress,maxvar,"var:eval_in_progress");
   for (int i = 0; i < maxvar; i++) eval_in_progress[i] = 0;
 }
@@ -743,14 +784,16 @@ double Variable::evaluate(char *str, Tree **tree)
 
       if (strncmp(word,"c_",2) == 0) {
         if (domain->box_exist == 0)
-          error->all(FLERR,"Variable evaluation before simulation box is defined");
+          error->all(FLERR,
+                     "Variable evaluation before simulation box is defined");
 
         n = strlen(word) - 2 + 1;
         char *id = new char[n];
         strcpy(id,&word[2]);
 
         int icompute = modify->find_compute(id);
-        if (icompute < 0) error->all(FLERR,"Invalid compute ID in variable formula");
+        if (icompute < 0) 
+          error->all(FLERR,"Invalid compute ID in variable formula");
         Compute *compute = modify->compute[icompute];
         delete [] id;
 
@@ -2826,7 +2869,7 @@ int Variable::special_function(char *word, char *contents, Tree **tree,
 
   if (strcmp(word,"sum") && strcmp(word,"min") && strcmp(word,"max") &&
       strcmp(word,"ave") && strcmp(word,"trap") && strcmp(word,"gmask") &&
-      strcmp(word,"rmask") && strcmp(word,"grmask"))
+      strcmp(word,"rmask") && strcmp(word,"grmask") && strcmp(word,"next"))
     return 0;
 
   // parse contents for arg1,arg2,arg3 separated by commas
@@ -2872,7 +2915,8 @@ int Variable::special_function(char *word, char *contents, Tree **tree,
     else if (strcmp(word,"ave") == 0) method = AVE;
     else if (strcmp(word,"trap") == 0) method = TRAP;
 
-    if (narg != 1) error->all(FLERR,"Invalid special function in variable formula");
+    if (narg != 1) 
+      error->all(FLERR,"Invalid special function in variable formula");
 
     Compute *compute = NULL;
     Fix *fix = NULL;
@@ -2887,12 +2931,14 @@ int Variable::special_function(char *word, char *contents, Tree **tree,
       } else index = 0;
 
       int icompute = modify->find_compute(&arg1[2]);
-      if (icompute < 0) error->all(FLERR,"Invalid compute ID in variable formula");
+      if (icompute < 0) 
+        error->all(FLERR,"Invalid compute ID in variable formula");
       compute = modify->compute[icompute];
       if (index == 0 && compute->vector_flag) {
         if (update->whichflag == 0) {
           if (compute->invoked_vector != update->ntimestep)
-            error->all(FLERR,"Compute used in variable between runs is not current");
+            error->all(FLERR,
+                       "Compute used in variable between runs is not current");
         } else if (!(compute->invoked_flag & INVOKED_VECTOR)) {
           compute->compute_vector();
           compute->invoked_flag |= INVOKED_VECTOR;
@@ -2961,7 +3007,7 @@ int Variable::special_function(char *word, char *contents, Tree **tree,
         else if (method == XMAX) value = MAX(value,vec[j]);
         else if (method == AVE) value += vec[j];
         else if (method == TRAP) {
-          if (i > 0 && i < nvec-1) value += vec[j];
+          if (i > 1 && i < nvec-1) value += vec[j];
           else value += 0.5*vec[j];
         }
         j += nstride;
@@ -2978,7 +3024,7 @@ int Variable::special_function(char *word, char *contents, Tree **tree,
         else if (method == XMAX) value = MAX(value,one);
         else if (method == AVE) value += one;
         else if (method == TRAP) {
-          if (i > 1 && i < nvec) value += one;
+          if (i > 1 && i < nvec-1) value += one;
           else value += 0.5*one;
         }
       }
@@ -2986,10 +3032,6 @@ int Variable::special_function(char *word, char *contents, Tree **tree,
 
     if (method == AVE) value /= nvec;
 
-    delete [] arg1;
-    delete [] arg2;
-    delete [] arg3;
-
     // save value in tree or on argstack
 
     if (tree) {
@@ -3005,7 +3047,8 @@ int Variable::special_function(char *word, char *contents, Tree **tree,
   } else if (strcmp(word,"gmask") == 0) {
     if (tree == NULL)
       error->all(FLERR,"Gmask function in equal-style variable formula");
-    if (narg != 1) error->all(FLERR,"Invalid special function in variable formula");
+    if (narg != 1) 
+      error->all(FLERR,"Invalid special function in variable formula");
 
     int igroup = group->find(arg1);
     if (igroup == -1)
@@ -3020,7 +3063,8 @@ int Variable::special_function(char *word, char *contents, Tree **tree,
   } else if (strcmp(word,"rmask") == 0) {
     if (tree == NULL)
       error->all(FLERR,"Rmask function in equal-style variable formula");
-    if (narg != 1) error->all(FLERR,"Invalid special function in variable formula");
+    if (narg != 1) 
+      error->all(FLERR,"Invalid special function in variable formula");
 
     int iregion = region_function(arg1);
 
@@ -3033,7 +3077,8 @@ int Variable::special_function(char *word, char *contents, Tree **tree,
   } else if (strcmp(word,"grmask") == 0) {
     if (tree == NULL)
       error->all(FLERR,"Grmask function in equal-style variable formula");
-    if (narg != 2) error->all(FLERR,"Invalid special function in variable formula");
+    if (narg != 2) 
+      error->all(FLERR,"Invalid special function in variable formula");
 
     int igroup = group->find(arg1);
     if (igroup == -1)
@@ -3046,8 +3091,37 @@ int Variable::special_function(char *word, char *contents, Tree **tree,
     newtree->ivalue2 = iregion;
     newtree->left = newtree->middle = newtree->right = NULL;
     treestack[ntreestack++] = newtree;
+
+  // file variable special function
+
+  } else if (strcmp(word,"next") == 0) {
+    if (narg != 1) 
+      error->all(FLERR,"Invalid special function in variable formula");
+
+    int ivar = find(arg1);
+    if (ivar == -1)
+      error->all(FLERR,"Variable ID in variable formula does not exist");
+    if (style[ivar] != FILEVAR)
+      error->all(FLERR,"Invalid variable in special function next");
+
+    double value = atof(data[ivar][0]);
+    reader[ivar]->read(data[ivar][0]);
+
+    // save value in tree or on argstack
+
+    if (tree) {
+      Tree *newtree = new Tree();
+      newtree->type = VALUE;
+      newtree->value = value;
+      newtree->left = newtree->middle = newtree->right = NULL;
+      treestack[ntreestack++] = newtree;
+    } else argstack[nargstack++] = value;
   }
 
+  delete [] arg1;
+  delete [] arg2;
+  delete [] arg3;
+
   return 1;
 }
 
@@ -3067,7 +3141,8 @@ void Variable::peratom2global(int flag, char *word,
                               double *argstack, int &nargstack)
 {
   if (atom->map_style == 0)
-    error->all(FLERR,"Indexed per-atom vector in variable formula without atom map");
+    error->all(FLERR,
+               "Indexed per-atom vector in variable formula without atom map");
 
   int index = atom->map(id);
 
@@ -3210,7 +3285,8 @@ double Variable::numeric(char *str)
     if (isdigit(str[i])) continue;
     if (str[i] == '-' || str[i] == '+' || str[i] == '.') continue;
     if (str[i] == 'e' || str[i] == 'E') continue;
-    error->all(FLERR,"Expected floating point parameter in variable definition");
+    error->all(FLERR,
+               "Expected floating point parameter in variable definition");
   }
 
   return atof(str);
@@ -3296,7 +3372,8 @@ double Variable::evaluate_boolean(char *str)
     // ----------------
 
     else if (onechar == '(') {
-      if (expect == OP) error->all(FLERR,"Invalid Boolean syntax in if command");
+      if (expect == OP) 
+        error->all(FLERR,"Invalid Boolean syntax in if command");
       expect = OP;
 
       char *contents;
@@ -3314,7 +3391,8 @@ double Variable::evaluate_boolean(char *str)
     // ----------------
 
     } else if (isdigit(onechar) || onechar == '.' || onechar == '-') {
-      if (expect == OP) error->all(FLERR,"Invalid Boolean syntax in if command");
+      if (expect == OP) 
+        error->all(FLERR,"Invalid Boolean syntax in if command");
       expect = OP;
 
       // istop = end of number, including scientific notation
@@ -3383,7 +3461,8 @@ double Variable::evaluate_boolean(char *str)
         continue;
       }
 
-      if (expect == ARG) error->all(FLERR,"Invalid Boolean syntax in if command");
+      if (expect == ARG) 
+        error->all(FLERR,"Invalid Boolean syntax in if command");
       expect = ARG;
 
       // evaluate stack as deep as possible while respecting precedence
@@ -3516,3 +3595,58 @@ unsigned int Variable::data_mask(char *str)
 
   return datamask;
 }
+
+/* ----------------------------------------------------------------------
+   class to read variable values from a file, line by line
+------------------------------------------------------------------------- */
+
+VarReader::VarReader(LAMMPS *lmp, char *file) : Pointers(lmp)
+{
+  MPI_Comm_rank(world,&me);
+
+  if (me == 0) {
+    fp = fopen(file,"r");
+    if (fp == NULL) {
+      char str[128];
+      sprintf(str,"Cannot open file variable file %s",file);
+      error->one(FLERR,str);
+    }
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
+VarReader::~VarReader()
+{
+  if (me == 0) fclose(fp);
+}
+
+/* ----------------------------------------------------------------------
+   read next value from file into str
+   strip comments, skip blank lines
+   return 0 if successful, 1 if end-of-file
+------------------------------------------------------------------------- */
+
+int VarReader::read(char *str)
+{
+  int n;
+  char *ptr;
+
+  if (me == 0) {
+    while (1) {
+      if (fgets(str,MAXLINE,fp) == NULL) n = 0;
+      else n = strlen(str);
+      if (n == 0) break;                                 // end of file
+      str[n-1] = '\0';                                   // strip newline
+      if (ptr = strchr(str,'#')) *ptr = '\0';            // strip comment
+      if (strtok(str," \t\n\r\f") == NULL) continue;     // skip if blank
+      n = strlen(str) + 1;
+      break;
+    }
+  }
+
+  MPI_Bcast(&n,1,MPI_INT,0,world);
+  if (n == 0) return 1;
+  MPI_Bcast(str,n,MPI_CHAR,0,world);
+  return 0;
+}

src/variable.h

@@ -14,6 +14,7 @@
 #ifndef LMP_VARIABLE_H
 #define LMP_VARIABLE_H
 
+#include "stdlib.h"
 #include "pointers.h"
 
 namespace LAMMPS_NS {
@@ -45,7 +46,9 @@ class Variable : protected Pointers {
   int *num;                // # of values for each variable
   int *which;              // next available value for each variable
   int *pad;                // 1 = pad loop/uloop variables with 0s, 0 = no pad
+  class VarReader **reader;   // variable that reads lines from file
   char ***data;            // str value of each variable's values
+
   int *eval_in_progress;   // flag if evaluation of variable is in progress
 
   class RanMars *randomequal;   // random number generator for equal-style vars
@@ -90,6 +93,17 @@ class Variable : protected Pointers {
   void print_tree(Tree *, int);
 };
 
+class VarReader : protected Pointers {
+ public:
+  VarReader(class LAMMPS *, char *);
+  ~VarReader();
+  int read(char *);
+
+ private:
+  int me;
+  FILE *fp;
+};
+
 }
 
 #endif

src/verlet.cpp

@@ -331,7 +331,7 @@ void Verlet::cleanup()
 
 /* ----------------------------------------------------------------------
    clear force on own & ghost atoms
-   setup and clear other arrays as needed
+   clear other arrays as needed
 ------------------------------------------------------------------------- */
 
 void Verlet::force_clear()

src/version.h

@@ -1 +1 @@
-#define LAMMPS_VERSION "28 Oct 2012"
+#define LAMMPS_VERSION "21 Nov 2012"