From 6295fc6908509f5b231f63445d88602a8a95d2e8 Mon Sep 17 00:00:00 2001
From: sjplimp <sjplimp@f3b2605a-c512-4ea7-a41b-209d697bcdaa>
Date: Thu, 16 Jul 2015 22:30:07 +0000
Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@13616
 f3b2605a-c512-4ea7-a41b-209d697bcdaa

---
 doc/compute_saed.html                 | 186 +++++++++++++++++++++++
 doc/compute_saed.txt                  | 173 +++++++++++++++++++++
 doc/compute_xrd.html                  | 210 ++++++++++++++++++++++++++
 doc/compute_xrd.txt                   | 196 ++++++++++++++++++++++++
 doc/fix_ave_histo.html                |  18 ++-
 doc/fix_ave_histo.txt                 |  22 ++-
 doc/fix_saed_vtk.html                 | 191 +++++++++++++++++++++++
 doc/fix_saed_vtk.txt                  | 181 ++++++++++++++++++++++
 src/MPIIO/dump_atom_mpiio.cpp         |   4 +-
 src/SRD/fix_wall_srd.cpp              |   2 +-
 src/STUBS/mpi.c                       |   2 +-
 src/USER-ATC/fix_atc.cpp              |   4 +
 src/USER-AWPMD/fix_nve_awpmd.cpp      |   4 -
 src/USER-AWPMD/pair_awpmd_cut.cpp     |   5 +-
 src/USER-COLVARS/colvarproxy_lammps.h |   2 +-
 src/USER-MISC/pair_tersoff_table.cpp  |   6 +-
 src/compute_cna_atom.cpp              |   2 +-
 17 files changed, 1187 insertions(+), 21 deletions(-)
 create mode 100644 doc/compute_saed.html
 create mode 100644 doc/compute_saed.txt
 create mode 100644 doc/compute_xrd.html
 create mode 100644 doc/compute_xrd.txt
 create mode 100644 doc/fix_saed_vtk.html
 create mode 100644 doc/fix_saed_vtk.txt

diff --git a/doc/compute_saed.html b/doc/compute_saed.html
new file mode 100644
index 0000000000..178de056b4
--- /dev/null
+++ b/doc/compute_saed.html
@@ -0,0 +1,186 @@
+<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>compute saed command 
+</H3>
+<P><B>Syntax:</B>
+</P>
+<PRE>compute ID group-ID saed lambda type1 type2 ... typeN keyword value ... 
+</PRE>
+<UL><LI>ID, group-ID are documented in <A HREF = "compute.html">compute</A> command 
+
+<LI>saed = style name of this compute command 
+
+<LI>lambda = wavelength of incident radiation (length units) 
+
+<LI>type1 type2 ... typeN = chemical symbol of each atom type (see valid 
+options below) 
+
+<LI>zero or more keyword/value pairs may be appended 
+
+<LI>keyword = <I>Kmax</I> or <I>Zone</I> or <I>dR_Ewald</I> or <I>c</I> or <I>manual</I> or <I>echo</I> 
+
+<PRE>  <I>Kmax</I> value = Maximum distance explored from reciprocal space origin 
+                 (inverse length units)
+  <I>Zone</I> values = z1 z2 z3
+    z1,z2,z3 = Zone axis of incident radiation. If z1=z2=z3=0 all 
+               reciprocal space will be meshed up to <I>Kmax</I>
+  <I>dR_Ewald</I> value = Thickness of Ewald sphere slice intercepting 
+                     reciprocal space (inverse length units)
+  <I>c</I> values = c1 c2 c3
+    c1,c2,c3 = parameters to adjust the spacing of the reciprocal 
+               lattice nodes in the h, k, and l directions respectively
+  <I>manual</I> = flag to use manual spacing of reciprocal lattice points   
+             based on the values of the <I>c</I> parameters 
+  <I>echo</I> = flag to provide extra output for debugging purposes 
+</PRE>
+
+</UL>
+<P><B>Examples:</B>
+</P>
+<P>compute 1 all saed 0.0251 Al O Kmax 1.70 Zone 0 0 1 dR_Ewald 0.01 c 0.5 0.5 0.5
+compute 2 all saed 0.0251 Ni Kmax 1.70 Zone 0 0 0 c 0.05 0.05 0.05 manual echo
+</P>
+<PRE>fix saed/vtk 1 1 1 c_1 file Al2O3_001.saed
+fix saed/vtk 1 1 1 c_2 file Ni_000.saed 
+</PRE>
+<P><B>Description:</B>
+</P>
+<P>Define a computation that calculates electron diffraction intensity as 
+described in <A HREF = "#Coleman">(Coleman)</A> on a mesh of reciprocal lattice nodes 
+defined by the entire simulation domain (or manually) using simulated 
+radiation of wavelength lambda. 
+</P>
+<P>The electron diffraction intensity I at each reciprocal lattice point 
+is computed from the structure factor F using the equations:
+</P>
+<CENTER><IMG SRC = "Eqs/compute_saed1.jpg">
+</CENTER>
+<CENTER><IMG SRC = "Eqs/compute_saed2.jpg">
+</CENTER>
+<P>Here, K is the location of the reciprocal lattice node, rj is the 
+position of each atom, fj are atomic scattering factors.
+</P>
+<P>Diffraction intensities are calculated on a three-dimensional mesh of 
+reciprocal lattice nodes. The mesh spacing is defined either (I)  by 
+the entire simulation domain or (II) manually using selected values.
+</P>
+<P>For a mesh defined by the simulation domain, a rectilinear grid is
+constructed with spacing <I>c</I>*inv(A) along each reciprocal lattice
+axis. Where A are the vectors corresponding to the edges of the
+simulation cell. If one or two directions has non-periodic boundary
+conditions, then the spacing in these directions is defined from the
+average of the (inversed) box lengths with periodic boundary conditions.
+Meshes defined by the simulation domain must contain at least one periodic
+boundary.
+</P>
+<P>If the <I>manual</I> flag is included, the mesh of reciprocal lattice nodes 
+will defined using the <I>c</I> values for the spacing along each reciprocal 
+lattice axis. Note that manual mapping of the reciprocal space mesh is 
+good for comparing diffraction results from  multiple simulations; however 
+it can reduce the likelihood that Bragg reflections will be satisfied 
+unless small spacing parameters <B><0.05 Angstrom^(-1)</B> are implemented. 
+Meshes with manual spacing do not require a periodic boundary.
+</P>
+<P>The limits of the reciprocal lattice mesh are determined by the use of 
+the <I>Kmax</I>, <I>Zone</I>, and <I>dR_Ewald</I> parameters.  The rectilinear mesh 
+created about the origin of reciprocal space is terminated at the 
+boundary of a sphere of radius <I>Kmax</I> centered at the origin.  If 
+<I>Zone</I> parameters z1=z2=z3=0 are used, diffraction intensities are 
+computed throughout the entire spherical volume - note this can greatly 
+increase the cost of computation.  Otherwise, <I>Zone</I> parameters will 
+denote the z1=h, z2=k, and z3=l (in a global since) zone axis of an 
+intersecting Ewald sphere.  Diffraction intensities will only be 
+computed at the intersection of the reciprocal lattice mesh and a 
+<I>dR_Ewald</I> thick surface of the Ewald sphere.
+</P>
+<P>The atomic scattering factors, fj, accounts for the reduction in 
+diffraction intensity due to Compton scattering.  Compute saed uses 
+analytical approximations of the atomic scattering factors that vary 
+for each atom type (type1 type2 ... typeN) and angle of diffraction.  
+The analytic approximation is computed using the formula
+<A HREF = "#Brown">(Brown)</A>:
+</P>
+<CENTER><IMG SRC = "Eqs/compute_saed3.jpg">
+</CENTER>
+<P>Coefficients parameterized by <A HREF = "#Fox">(Fox)</A> are assigned for each 
+atom type designating the chemical symbol and charge of each atom 
+type. Valid chemical symbols for compute saed are:
+</P>
+<P>         H:      He:      Li:      Be:       B:
+         C:       N:       O:       F:      Ne:
+        Na:      Mg:      Al:      Si:       P:
+         S:      Cl:      Ar:       K:      Ca:
+        Sc:      Ti:       V:      Cr:      Mn:
+        Fe:      Co:      Ni:      Cu:      Zn:
+        Ga:      Ge:      As:      Se:      Br:
+        Kr:      Rb:      Sr:       Y:      Zr:
+        Nb:      Mo:      Tc:      Ru:      Rh:
+        Pd:      Ag:      Cd:      In:      Sn:
+        Sb:      Te:       I:      Xe:      Cs:
+        Ba:      La:      Ce:      Pr:      Nd:
+        Pm:      Sm:      Eu:      Gd:      Tb:
+        Dy:      Ho:      Er:      Tm:      Yb:
+        Lu:      Hf:      Ta:       W:      Re:
+        Os:      Ir:      Pt:      Au:      Hg:
+        Tl:      Pb:      Bi:      Po:      At:
+        Rn:      Fr:      Ra:      Ac:      Th:
+        Pa:       U:      Np:      Pu:      Am:
+        Cm:      Bk:      Cf:
+</P>
+<P>If the <I>echo</I> keyword is specified, compute saed will provide extra 
+reporting information to the screen.  
+</P>
+<P><B>Output info:</B>
+</P>
+<P>This compute calculates a global vector.  The length of the vector is 
+the number of reciprocal lattice nodes that are explored by the mesh.  
+The entries of the global vector are the computed diffraction 
+intensities as described above.
+</P>
+<P>The vector can be accessed by any command that uses global values
+from a compute as input.  See <A HREF = "Section_howto.html#howto_15">this
+section</A> for an overview of LAMMPS output
+options.
+</P>
+<P>All array values calculated by this compute are "intensive".  
+</P>
+<P><B>Restrictions:</B> 
+</P>
+<P>The compute_saed command does not work for triclinic cells. 
+</P>
+<P><B>Related commands:</B> 
+</P>
+<P><A HREF = "fix_saed_vtk.html">fix saed_vtk</A>
+<A HREF = "compute_xrd.html">compute xrd</A>
+</P>
+<P><B>Default:</B> 
+</P>
+<P>The option defaults are Kmax = 1.70, Zone 1 0 0, c 1 1 1, dR_Ewald = 0.01
+</P>
+<HR>
+
+<A NAME = "Coleman"></A>
+
+<P><B>(Coleman)</B> Coleman, Spearot, Capolungo, MSMSE, 21, 055020
+(2013).
+</P>
+<A NAME = "Brown"></A>
+
+<P><B>(Brown)</B> Brown et al. International Tables for Crystallography 
+Volume C: Mathematical and Chemical Tables, 554-95 (2004).
+</P>
+<A NAME = "Fox"></A>
+
+<P><B>(Fox)</B> Fox, O'Keefe, Tabbernor, Acta Crystallogr. A, 45, 786-93
+(1989).
+</P>
+</HTML>
diff --git a/doc/compute_saed.txt b/doc/compute_saed.txt
new file mode 100644
index 0000000000..d2aaae2056
--- /dev/null
+++ b/doc/compute_saed.txt
@@ -0,0 +1,173 @@
+"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
+
+compute saed command :h3
+
+[Syntax:]
+
+compute ID group-ID saed lambda type1 type2 ... typeN keyword value ... :pre
+
+ID, group-ID are documented in "compute"_compute.html command :ulb,l
+saed = style name of this compute command :l
+lambda = wavelength of incident radiation (length units) :l
+type1 type2 ... typeN = chemical symbol of each atom type (see valid 
+options below) :l
+
+zero or more keyword/value pairs may be appended :l
+keyword = {Kmax} or {Zone} or {dR_Ewald} or {c} or {manual} or {echo} :l
+  {Kmax} value = Maximum distance explored from reciprocal space origin 
+                 (inverse length units)
+  {Zone} values = z1 z2 z3
+    z1,z2,z3 = Zone axis of incident radiation. If z1=z2=z3=0 all 
+               reciprocal space will be meshed up to {Kmax}
+  {dR_Ewald} value = Thickness of Ewald sphere slice intercepting 
+                     reciprocal space (inverse length units)
+  {c} values = c1 c2 c3
+    c1,c2,c3 = parameters to adjust the spacing of the reciprocal 
+               lattice nodes in the h, k, and l directions respectively
+  {manual} = flag to use manual spacing of reciprocal lattice points   
+             based on the values of the {c} parameters 
+  {echo} = flag to provide extra output for debugging purposes :pre
+:ule
+
+[Examples:]
+
+compute 1 all saed 0.0251 Al O Kmax 1.70 Zone 0 0 1 dR_Ewald 0.01 c 0.5 0.5 0.5
+compute 2 all saed 0.0251 Ni Kmax 1.70 Zone 0 0 0 c 0.05 0.05 0.05 manual echo
+
+fix saed/vtk 1 1 1 c_1 file Al2O3_001.saed
+fix saed/vtk 1 1 1 c_2 file Ni_000.saed :pre
+
+[Description:]
+
+Define a computation that calculates electron diffraction intensity as 
+described in "(Coleman)"_#Coleman on a mesh of reciprocal lattice nodes 
+defined by the entire simulation domain (or manually) using simulated 
+radiation of wavelength lambda. 
+
+The electron diffraction intensity I at each reciprocal lattice point 
+is computed from the structure factor F using the equations:
+
+:c,image(Eqs/compute_saed1.jpg) 
+:c,image(Eqs/compute_saed2.jpg)
+
+Here, K is the location of the reciprocal lattice node, rj is the 
+position of each atom, fj are atomic scattering factors.
+
+Diffraction intensities are calculated on a three-dimensional mesh of 
+reciprocal lattice nodes. The mesh spacing is defined either (I)  by 
+the entire simulation domain or (II) manually using selected values.
+
+For a mesh defined by the simulation domain, a rectilinear grid is
+constructed with spacing {c}*inv(A) along each reciprocal lattice
+axis. Where A are the vectors corresponding to the edges of the
+simulation cell. If one or two directions has non-periodic boundary
+conditions, then the spacing in these directions is defined from the
+average of the (inversed) box lengths with periodic boundary conditions.
+Meshes defined by the simulation domain must contain at least one periodic
+boundary.
+
+If the {manual} flag is included, the mesh of reciprocal lattice nodes 
+will defined using the {c} values for the spacing along each reciprocal 
+lattice axis. Note that manual mapping of the reciprocal space mesh is 
+good for comparing diffraction results from  multiple simulations; however 
+it can reduce the likelihood that Bragg reflections will be satisfied 
+unless small spacing parameters [<0.05 Angstrom^(-1)] are implemented. 
+Meshes with manual spacing do not require a periodic boundary.
+
+The limits of the reciprocal lattice mesh are determined by the use of 
+the {Kmax}, {Zone}, and {dR_Ewald} parameters.  The rectilinear mesh 
+created about the origin of reciprocal space is terminated at the 
+boundary of a sphere of radius {Kmax} centered at the origin.  If 
+{Zone} parameters z1=z2=z3=0 are used, diffraction intensities are 
+computed throughout the entire spherical volume - note this can greatly 
+increase the cost of computation.  Otherwise, {Zone} parameters will 
+denote the z1=h, z2=k, and z3=l (in a global since) zone axis of an 
+intersecting Ewald sphere.  Diffraction intensities will only be 
+computed at the intersection of the reciprocal lattice mesh and a 
+{dR_Ewald} thick surface of the Ewald sphere.
+
+The atomic scattering factors, fj, accounts for the reduction in 
+diffraction intensity due to Compton scattering.  Compute saed uses 
+analytical approximations of the atomic scattering factors that vary 
+for each atom type (type1 type2 ... typeN) and angle of diffraction.  
+The analytic approximation is computed using the formula
+"(Brown)"_#Brown:
+
+:c,image(Eqs/compute_saed3.jpg)
+
+Coefficients parameterized by "(Fox)"_#Fox are assigned for each 
+atom type designating the chemical symbol and charge of each atom 
+type. Valid chemical symbols for compute saed are:
+
+         H:      He:      Li:      Be:       B:
+         C:       N:       O:       F:      Ne:
+        Na:      Mg:      Al:      Si:       P:
+         S:      Cl:      Ar:       K:      Ca:
+        Sc:      Ti:       V:      Cr:      Mn:
+        Fe:      Co:      Ni:      Cu:      Zn:
+        Ga:      Ge:      As:      Se:      Br:
+        Kr:      Rb:      Sr:       Y:      Zr:
+        Nb:      Mo:      Tc:      Ru:      Rh:
+        Pd:      Ag:      Cd:      In:      Sn:
+        Sb:      Te:       I:      Xe:      Cs:
+        Ba:      La:      Ce:      Pr:      Nd:
+        Pm:      Sm:      Eu:      Gd:      Tb:
+        Dy:      Ho:      Er:      Tm:      Yb:
+        Lu:      Hf:      Ta:       W:      Re:
+        Os:      Ir:      Pt:      Au:      Hg:
+        Tl:      Pb:      Bi:      Po:      At:
+        Rn:      Fr:      Ra:      Ac:      Th:
+        Pa:       U:      Np:      Pu:      Am:
+        Cm:      Bk:      Cf:
+
+
+If the {echo} keyword is specified, compute saed will provide extra 
+reporting information to the screen.  
+
+[Output info:]
+
+This compute calculates a global vector.  The length of the vector is 
+the number of reciprocal lattice nodes that are explored by the mesh.  
+The entries of the global vector are the computed diffraction 
+intensities as described above.
+
+The vector can be accessed by any command that uses global values
+from a compute as input.  See "this
+section"_Section_howto.html#howto_15 for an overview of LAMMPS output
+options.
+
+All array values calculated by this compute are "intensive".  
+
+[Restrictions:] 
+
+The compute_saed command does not work for triclinic cells. 
+
+[Related commands:] 
+
+"fix saed_vtk"_fix_saed_vtk.html
+"compute xrd"_compute_xrd.html
+
+[Default:] 
+
+The option defaults are Kmax = 1.70, Zone 1 0 0, c 1 1 1, dR_Ewald = 0.01
+
+:line
+
+:link(Coleman)
+[(Coleman)] Coleman, Spearot, Capolungo, MSMSE, 21, 055020
+(2013).
+
+:link(Brown)
+[(Brown)] Brown et al. International Tables for Crystallography 
+Volume C: Mathematical and Chemical Tables, 554-95 (2004).
+
+:link(Fox)
+[(Fox)] Fox, O'Keefe, Tabbernor, Acta Crystallogr. A, 45, 786-93
+(1989).
+
diff --git a/doc/compute_xrd.html b/doc/compute_xrd.html
new file mode 100644
index 0000000000..bbcd1ddfab
--- /dev/null
+++ b/doc/compute_xrd.html
@@ -0,0 +1,210 @@
+<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>compute xrd command 
+</H3>
+<P><B>Syntax:</B>
+</P>
+<PRE>compute ID group-ID xrd lambda type1 type2 ... typeN keyword value ... 
+</PRE>
+<UL><LI>ID, group-ID are documented in <A HREF = "compute.html">compute</A> command 
+
+<LI>xrd = style name of this compute command 
+
+<LI>lambda = wavelength of incident radiation (length units) 
+
+<LI>type1 type2 ... typeN = chemical symbol of each atom type (see valid 
+options below) 
+
+<LI>zero or more keyword/value pairs may be appended 
+
+<LI>keyword = <I>2Theta</I> or <I>c</I> or <I>LP</I> or <I>manual</I> or <I>echo</I> 
+
+<PRE>  <I>2Theta</I> values = Min2Theta Max2Theta
+    Min2Theta,Max2Theta = minimum and maximum 2 theta range to explore 
+    (radians or degrees)
+  <I>c</I> values = c1 c2 c3
+    c1,c2,c3 = parameters to adjust the spacing of the reciprocal 
+               lattice nodes in the h, k, and l directions respectively
+  <I>LP</I> value = switch to apply Lorentz-polarization factor
+    0/1 = off/on
+  <I>manual</I> = flag to use manual spacing of reciprocal lattice points 
+             based on the values of the <I>c</I> parameters 
+  <I>echo</I> = flag to provide extra output for debugging purposes 
+</PRE>
+
+</UL>
+<P><B>Examples:</B>
+</P>
+<P>compute 1 all xrd 1.541838 Al O 2Theta 0.087 0.87 c 1 1 1 LP 1 echo 
+compute 2 all xrd 1.541838 Al O 2Theta 10 100 c 0.05 0.05 0.05 LP 1 manual 
+</P>
+<P>fix 1 all ave/histo 1 1 1 0.087 0.87 250 c_1<B>1</B> mode vector weights c_1<B>2</B> file Rad2Theta.xrd
+fix 2 all ave/histo 1 1 1 10 100 250 c_2<B>1</B> mode vector weights c_2<B>2</B> file Deg2Theta.xrd
+</P>
+<PRE>
+</PRE>
+<P><B>Description:</B>
+</P>
+<P>Define a computation that calculates x-ray diffraction intensity as described
+in <A HREF = "#Coleman">(Coleman)</A> on a mesh of reciprocal lattice nodes defined 
+by the entire simulation domain (or manually) using simulated radiation
+of wavelength lambda. 
+</P>
+<P>The x-ray diffraction intensity I at each reciprocal lattice point 
+is computed from the structure factor F using the equations:
+</P>
+<CENTER><IMG SRC = "Eqs/compute_xrd1.jpg">
+</CENTER>
+<CENTER><IMG SRC = "Eqs/compute_xrd2.jpg">
+</CENTER>
+<CENTER><IMG SRC = "Eqs/compute_xrd3.jpg">
+</CENTER>
+<CENTER><IMG SRC = "Eqs/compute_xrd4.jpg">
+</CENTER>
+<P>Here, K is the location of the reciprocal lattice node, rj is the 
+position of each atom, fj are atomic scattering factors, LP is the 
+Lorentz-polarization factor, and theta is the scattering angle of 
+diffraction.  The Lorentz-polarization factor can be turned off using 
+the <I>LP</I> switch.
+</P>
+<P>Diffraction intensities are calculated on a three-dimensional mesh of 
+reciprocal lattice nodes. The mesh spacing is defined either (I) 
+by the entire simulation domain or (II) manually using selected values.
+</P>
+<P>For a mesh defined by the simulation domain, a rectilinear grid is
+constructed with spacing <I>c</I>*inv(A) along each reciprocal lattice
+axis. Where A are the vectors corresponding to the edges of the
+simulation cell. If one or two directions has non-periodic boundary
+conditions, then the spacing in these directions is defined from the
+average of the (inversed) box lengths with periodic boundary conditions.
+Meshes defined by the simulation domain must contain at least one periodic
+boundary.
+</P>
+<P>If the <I>manual</I> flag is included, the mesh of reciprocal lattice nodes 
+will defined using the <I>c</I> values for the spacing along each reciprocal 
+lattice axis. Note that manual mapping of the reciprocal space mesh is 
+good for comparing diffraction results from  multiple simulations; however 
+it can reduce the likelihood that Bragg reflections will be satisfied 
+unless small spacing parameters <B><0.05 Angstrom^(-1)</B> are implemented. 
+Meshes with manual spacing do not require a periodic boundary.
+</P>
+<P>The limits of the reciprocal lattice mesh are determined by range of 
+scattering angles explored.  The <I>2Theta</I> parameters allows the user to 
+reduce the scattering angle range to only the region of interest which 
+reduces the cost of the computation.
+</P>
+<P>The atomic scattering factors, fj, accounts for the reduction in 
+diffraction intensity due to Compton scattering.  Compute xrd uses 
+analytical approximations of the atomic scattering factors that vary 
+for each atom type (type1 type2 ... typeN) and angle of diffraction.  
+The analytic approximation is computed using the formula
+<A HREF = "#Colliex">(Colliex)</A>:
+</P>
+<CENTER><IMG SRC = "Eqs/compute_xrd5.jpg">
+</CENTER>
+<P>Coefficients parameterized by <A HREF = "#Peng">(Peng)</A> are assigned for each 
+atom type designating the chemical symbol and charge of each atom 
+type. Valid chemical symbols for compute xrd are:
+</P>
+<P>         H:    He1-:      He:      Li:    Li1+:
+        Be:    Be2+:       B:       C:    Cval:
+         N:       O:     O1-:       F:     F1-:
+        Ne:      Na:    Na1+:      Mg:    Mg2+:
+        Al:    Al3+:      Si:    Sival:   Si4+:
+         P:       S:      Cl:    Cl1-:      Ar:
+         K:      Ca:    Ca2+:      Sc:    Sc3+:
+        Ti:    Ti2+:    Ti3+:    Ti4+:       V:
+       V2+:     V3+:     V5+:      Cr:    Cr2+:
+      Cr3+:      Mn:    Mn2+:    Mn3+:    Mn4+:
+        Fe:    Fe2+:    Fe3+:      Co:    Co2+:
+        Co:      Ni:    Ni2+:    Ni3+:      Cu:
+      Cu1+:    Cu2+:      Zn:    Zn2+:      Ga:
+      Ga3+:      Ge:    Ge4+:      As:      Se:
+        Br:    Br1-:      Kr:      Rb:    Rb1+:
+        Sr:    Sr2+:       Y:     Y3+:      Zr:
+      Zr4+:      Nb:    Nb3+:    Nb5+:      Mo:
+      Mo3+:    Mo5+:    Mo6+:      Tc:      Ru:
+      Ru3+:    Ru4+:      Rh:    Rh3+:    Rh4+:
+        Pd:    Pd2+:    Pd4+:      Ag:    Ag1+:
+      Ag2+:      Cd:    Cd2+:      In:    In3+:
+        Sn:    Sn2+:    Sn4+:      Sb:    Sb3+:
+      Sb5+:      Te:       I:     I1-:      Xe:
+        Cs:    Cs1+:      Ba:    Ba2+:      La:
+      La3+:      Ce:    Ce3+:    Ce4+:      Pr:
+      Pr3+:    Pr4+:      Nd:    Nd3+:      Pm:
+      Pm3+:      Sm:    Sm3+:      Eu:    Eu2+:
+      Eu3+:      Gd:    Gd3+:      Tb:    Tb3+:
+        Dy:    Dy3+:      Ho:    Ho3+:      Er:
+      Er3+:      Tm:    Tm3+:      Yb:    Yb2+:
+      Yb3+:      Lu:    Lu3+:      Hf:    Hf4+:
+        Ta:    Ta5+:       W:     W6+:      Re:
+        Os:    Os4+:      Ir:    Ir3+:    Ir4+:
+        Pt:    Pt2+:    Pt4+:      Au:    Au1+:
+      Au3+:      Hg:    Hg1+:    Hg2+:      Tl:
+      Tl1+:    Tl3+:      Pb:    Pb2+:    Pb4+:
+        Bi:    Bi3+:    Bi5+:      Po:      At:
+        Rn:      Fr:      Ra:    Ra2+:      Ac:
+      Ac3+:      Th:    Th4+:      Pa:       U:
+       U3+:     U4+:     U6+:      Np:    Np3+:
+      Np4+:    Np6+:      Pu:    Pu3+:    Pu4+:
+      Pu6+:      Am:      Cm:      Bk:      Cf:
+</P>
+<P>If the <I>echo</I> keyword is specified, compute xrd will provide extra 
+reporting information to the screen.  
+</P>
+<P><B>Output info:</B>
+</P>
+<P>This compute calculates a global array.  The number of rows in the 
+array is the number of reciprocal lattice nodes that are explored 
+which by the mesh.  The global array has 2 columns.  
+</P>
+<P>The first column contains the diffraction angle in the units (radians
+or degrees) provided with the <I>2Theta</I> values. The second column contains 
+the computed diffraction intensities as described above.
+</P>
+<P>The array can be accessed by any command that uses global values
+from a compute as input.  See <A HREF = "Section_howto.html#howto_15">this section</A> for an overview of LAMMPS output
+options.
+</P>
+<P>All array values calculated by this compute are "intensive".  
+</P>
+<P><B>Restrictions:</B> 
+</P>
+<P>The compute_xrd command does not work for triclinic cells. 
+</P>
+<P><B>Related commands:</B> 
+</P>
+<P><A HREF = "compute_ave_histo.html">fix ave/histo</A>
+<A HREF = "compute_saed.html">compute saed</A>
+</P>
+<P><B>Default:</B> 
+</P>
+<P>The option defaults are 2Theta = 1 179 (degrees), c = 1 1 1, LP = 1, no manual flag, no echo flag
+</P>
+<HR>
+
+<A NAME = "Coleman"></A>
+
+<P><B>(Coleman)</B> Coleman, Spearot, Capolungo, MSMSE, 21, 055020
+(2013).
+</P>
+<A NAME = "Colliex"></A>
+
+<P><B>(Colliex)</B> Colliex et al. International Tables for Crystallography 
+Volume C: Mathematical and Chemical Tables, 249-429 (2004).
+</P>
+<A NAME = "Peng"></A>
+
+<P><B>(Peng)</B> Peng, Ren, Dudarev, Whelan, Acta Crystallogr. A, 52, 257-76
+(1996).
+</P>
+</HTML>
diff --git a/doc/compute_xrd.txt b/doc/compute_xrd.txt
new file mode 100644
index 0000000000..4e03ef5709
--- /dev/null
+++ b/doc/compute_xrd.txt
@@ -0,0 +1,196 @@
+"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
+
+compute xrd command :h3
+
+[Syntax:]
+
+compute ID group-ID xrd lambda type1 type2 ... typeN keyword value ... :pre
+
+ID, group-ID are documented in "compute"_compute.html command :ulb,l
+xrd = style name of this compute command :l
+lambda = wavelength of incident radiation (length units) :l
+type1 type2 ... typeN = chemical symbol of each atom type (see valid 
+options below) :l
+
+zero or more keyword/value pairs may be appended :l
+keyword = {2Theta} or {c} or {LP} or {manual} or {echo} :l
+  {2Theta} values = Min2Theta Max2Theta
+    Min2Theta,Max2Theta = minimum and maximum 2 theta range to explore 
+    (radians or degrees)
+  {c} values = c1 c2 c3
+    c1,c2,c3 = parameters to adjust the spacing of the reciprocal 
+               lattice nodes in the h, k, and l directions respectively
+  {LP} value = switch to apply Lorentz-polarization factor
+    0/1 = off/on
+  {manual} = flag to use manual spacing of reciprocal lattice points 
+             based on the values of the {c} parameters 
+  {echo} = flag to provide extra output for debugging purposes :pre
+:ule
+
+[Examples:]
+
+compute 1 all xrd 1.541838 Al O 2Theta 0.087 0.87 c 1 1 1 LP 1 echo 
+compute 2 all xrd 1.541838 Al O 2Theta 10 100 c 0.05 0.05 0.05 LP 1 manual 
+
+fix 1 all ave/histo 1 1 1 0.087 0.87 250 c_1[1] mode vector weights c_1[2] file Rad2Theta.xrd
+fix 2 all ave/histo 1 1 1 10 100 250 c_2[1] mode vector weights c_2[2] file Deg2Theta.xrd
+
+:pre
+
+[Description:]
+
+Define a computation that calculates x-ray diffraction intensity as described
+in "(Coleman)"_#Coleman on a mesh of reciprocal lattice nodes defined 
+by the entire simulation domain (or manually) using simulated radiation
+of wavelength lambda. 
+
+The x-ray diffraction intensity I at each reciprocal lattice point 
+is computed from the structure factor F using the equations:
+
+:c,image(Eqs/compute_xrd1.jpg)
+:c,image(Eqs/compute_xrd2.jpg)
+:c,image(Eqs/compute_xrd3.jpg)
+:c,image(Eqs/compute_xrd4.jpg)
+
+Here, K is the location of the reciprocal lattice node, rj is the 
+position of each atom, fj are atomic scattering factors, LP is the 
+Lorentz-polarization factor, and theta is the scattering angle of 
+diffraction.  The Lorentz-polarization factor can be turned off using 
+the {LP} switch.
+
+Diffraction intensities are calculated on a three-dimensional mesh of 
+reciprocal lattice nodes. The mesh spacing is defined either (I) 
+by the entire simulation domain or (II) manually using selected values.
+
+For a mesh defined by the simulation domain, a rectilinear grid is
+constructed with spacing {c}*inv(A) along each reciprocal lattice
+axis. Where A are the vectors corresponding to the edges of the
+simulation cell. If one or two directions has non-periodic boundary
+conditions, then the spacing in these directions is defined from the
+average of the (inversed) box lengths with periodic boundary conditions.
+Meshes defined by the simulation domain must contain at least one periodic
+boundary.
+
+If the {manual} flag is included, the mesh of reciprocal lattice nodes 
+will defined using the {c} values for the spacing along each reciprocal 
+lattice axis. Note that manual mapping of the reciprocal space mesh is 
+good for comparing diffraction results from  multiple simulations; however 
+it can reduce the likelihood that Bragg reflections will be satisfied 
+unless small spacing parameters [<0.05 Angstrom^(-1)] are implemented. 
+Meshes with manual spacing do not require a periodic boundary.
+
+The limits of the reciprocal lattice mesh are determined by range of 
+scattering angles explored.  The {2Theta} parameters allows the user to 
+reduce the scattering angle range to only the region of interest which 
+reduces the cost of the computation.
+
+The atomic scattering factors, fj, accounts for the reduction in 
+diffraction intensity due to Compton scattering.  Compute xrd uses 
+analytical approximations of the atomic scattering factors that vary 
+for each atom type (type1 type2 ... typeN) and angle of diffraction.  
+The analytic approximation is computed using the formula
+"(Colliex)"_#Colliex:
+
+:c,image(Eqs/compute_xrd5.jpg)
+
+Coefficients parameterized by "(Peng)"_#Peng are assigned for each 
+atom type designating the chemical symbol and charge of each atom 
+type. Valid chemical symbols for compute xrd are:
+
+         H:    He1-:      He:      Li:    Li1+:
+        Be:    Be2+:       B:       C:    Cval:
+         N:       O:     O1-:       F:     F1-:
+        Ne:      Na:    Na1+:      Mg:    Mg2+:
+        Al:    Al3+:      Si:    Sival:   Si4+:
+         P:       S:      Cl:    Cl1-:      Ar:
+         K:      Ca:    Ca2+:      Sc:    Sc3+:
+        Ti:    Ti2+:    Ti3+:    Ti4+:       V:
+       V2+:     V3+:     V5+:      Cr:    Cr2+:
+      Cr3+:      Mn:    Mn2+:    Mn3+:    Mn4+:
+        Fe:    Fe2+:    Fe3+:      Co:    Co2+:
+        Co:      Ni:    Ni2+:    Ni3+:      Cu:
+      Cu1+:    Cu2+:      Zn:    Zn2+:      Ga:
+      Ga3+:      Ge:    Ge4+:      As:      Se:
+        Br:    Br1-:      Kr:      Rb:    Rb1+:
+        Sr:    Sr2+:       Y:     Y3+:      Zr:
+      Zr4+:      Nb:    Nb3+:    Nb5+:      Mo:
+      Mo3+:    Mo5+:    Mo6+:      Tc:      Ru:
+      Ru3+:    Ru4+:      Rh:    Rh3+:    Rh4+:
+        Pd:    Pd2+:    Pd4+:      Ag:    Ag1+:
+      Ag2+:      Cd:    Cd2+:      In:    In3+:
+        Sn:    Sn2+:    Sn4+:      Sb:    Sb3+:
+      Sb5+:      Te:       I:     I1-:      Xe:
+        Cs:    Cs1+:      Ba:    Ba2+:      La:
+      La3+:      Ce:    Ce3+:    Ce4+:      Pr:
+      Pr3+:    Pr4+:      Nd:    Nd3+:      Pm:
+      Pm3+:      Sm:    Sm3+:      Eu:    Eu2+:
+      Eu3+:      Gd:    Gd3+:      Tb:    Tb3+:
+        Dy:    Dy3+:      Ho:    Ho3+:      Er:
+      Er3+:      Tm:    Tm3+:      Yb:    Yb2+:
+      Yb3+:      Lu:    Lu3+:      Hf:    Hf4+:
+        Ta:    Ta5+:       W:     W6+:      Re:
+        Os:    Os4+:      Ir:    Ir3+:    Ir4+:
+        Pt:    Pt2+:    Pt4+:      Au:    Au1+:
+      Au3+:      Hg:    Hg1+:    Hg2+:      Tl:
+      Tl1+:    Tl3+:      Pb:    Pb2+:    Pb4+:
+        Bi:    Bi3+:    Bi5+:      Po:      At:
+        Rn:      Fr:      Ra:    Ra2+:      Ac:
+      Ac3+:      Th:    Th4+:      Pa:       U:
+       U3+:     U4+:     U6+:      Np:    Np3+:
+      Np4+:    Np6+:      Pu:    Pu3+:    Pu4+:
+      Pu6+:      Am:      Cm:      Bk:      Cf:
+
+
+If the {echo} keyword is specified, compute xrd will provide extra 
+reporting information to the screen.  
+
+[Output info:]
+
+This compute calculates a global array.  The number of rows in the 
+array is the number of reciprocal lattice nodes that are explored 
+which by the mesh.  The global array has 2 columns.  
+
+The first column contains the diffraction angle in the units (radians
+or degrees) provided with the {2Theta} values. The second column contains 
+the computed diffraction intensities as described above.
+
+The array can be accessed by any command that uses global values
+from a compute as input.  See "this section"_Section_howto.html#howto_15 for an overview of LAMMPS output
+options.
+
+All array values calculated by this compute are "intensive".  
+
+[Restrictions:] 
+
+The compute_xrd command does not work for triclinic cells. 
+
+[Related commands:] 
+
+"fix ave/histo"_compute_ave_histo.html
+"compute saed"_compute_saed.html
+
+[Default:] 
+
+The option defaults are 2Theta = 1 179 (degrees), c = 1 1 1, LP = 1, no manual flag, no echo flag
+
+:line
+
+:link(Coleman)
+[(Coleman)] Coleman, Spearot, Capolungo, MSMSE, 21, 055020
+(2013).
+
+:link(Colliex)
+[(Colliex)] Colliex et al. International Tables for Crystallography 
+Volume C: Mathematical and Chemical Tables, 249-429 (2004).
+
+:link(Peng)
+[(Peng)] Peng, Ren, Dudarev, Whelan, Acta Crystallogr. A, 52, 257-76
+(1996).
+
+
diff --git a/doc/fix_ave_histo.html b/doc/fix_ave_histo.html
index b4d3879253..7533cc70a5 100644
--- a/doc/fix_ave_histo.html
+++ b/doc/fix_ave_histo.html
@@ -42,7 +42,7 @@
 </PRE>
 <LI>zero or more keyword/arg pairs may be appended 
 
-<LI>keyword = <I>mode</I> or <I>file</I> or <I>ave</I> or <I>start</I> or <I>beyond</I> or <I>overwrite</I> or <I>title1</I> or <I>title2</I> or <I>title3</I> 
+<LI>keyword = <I>mode</I> or <I>file</I> or <I>ave</I> or <I>start</I> or <I>beyond</I> or <I>overwrite</I> or <I>title1</I> or <I>title2</I> or <I>title3</I> or <I>weights</I> 
 
 <PRE>  <I>mode</I> arg = <I>scalar</I> or <I>vector</I>
     scalar = all input values are scalars
@@ -66,6 +66,12 @@
     string = text to print as 2nd line of output file
   <I>title3</I> arg = string
     string = text to print as 3rd line of output file, only for vector mode 
+  <I>weights</I> arg = c_ID, c_ID[N], f_ID, f_ID[N], v_name 
+    c_ID = scalar or vector calculated by a compute with ID
+    c_ID[I] = Ith component of vector or Ith column of array calculated by a compute with ID
+    f_ID = scalar or vector calculated by a fix with ID
+    f_ID[I] = Ith component of vector or Ith column of array calculated by a fix with ID
+    v_name = value(s) calculated by an equal-style or atom-style variable with name 
 </PRE>
 
 </UL>
@@ -74,6 +80,7 @@
 <PRE>fix 1 all ave/histo 100 5 1000 0.5 1.5 50 c_myTemp file temp.histo ave running
 fix 1 all ave/histo 100 5 1000 -5 5 100 c_thermo_press[2] c_thermo_press[3] title1 "My output values"
 fix 1 all ave/histo 1 100 1000 -2.0 2.0 18 vx vy vz mode vector ave running beyond extra 
+fix 1 all ave/histo 1 1 1 10 100 2000 c_XRD<B>1</B> mode vector weights c_XRD<B>2</B> 
 </PRE>
 <P><B>Description:</B>
 </P>
@@ -287,6 +294,15 @@ describes the six values that are printed at the first of each section
 of output.  The third describes the 4 values printed for each bin in
 the histogram.
 </P>
+<P>If the <I>weights</I> keyword is specified, the fix will compute a weighted
+histogram using per-bin weights specified by the <I>weights</I> argument.  As
+normal, the bin locations will be will be generated based off value1.
+However, instead of each binned value contributing 1 to the bin
+location, the contributing ammount is assigned the weights
+argument. Only a single value1 and weights argument pair can be can be
+used for each fix ave/histo.  The length of value1 must match the
+lenght of the weights arguemnt.
+</P>
 <HR>
 
 <P><B>Restart, fix_modify, output, run start/stop, minimize info:</B>
diff --git a/doc/fix_ave_histo.txt b/doc/fix_ave_histo.txt
index 62e7c837b4..a954c29325 100644
--- a/doc/fix_ave_histo.txt
+++ b/doc/fix_ave_histo.txt
@@ -29,7 +29,7 @@ value = x, y, z, vx, vy, vz, fx, fy, fz, c_ID, c_ID\[N\], f_ID, f_ID\[N\], v_nam
   v_name = value(s) calculated by an equal-style or atom-style variable with name :pre
 
 zero or more keyword/arg pairs may be appended :l
-keyword = {mode} or {file} or {ave} or {start} or {beyond} or {overwrite} or {title1} or {title2} or {title3} :l
+keyword = {mode} or {file} or {ave} or {start} or {beyond} or {overwrite} or {title1} or {title2} or {title3} or {weights} :l
   {mode} arg = {scalar} or {vector}
     scalar = all input values are scalars
     vector = all input values are vectors
@@ -51,14 +51,21 @@ keyword = {mode} or {file} or {ave} or {start} or {beyond} or {overwrite} or {ti
   {title2} arg = string
     string = text to print as 2nd line of output file
   {title3} arg = string
-    string = text to print as 3rd line of output file, only for vector mode :pre
+    string = text to print as 3rd line of output file, only for vector mode 
+  {weights} arg = c_ID, c_ID\[N\], f_ID, f_ID\[N\], v_name 
+    c_ID = scalar or vector calculated by a compute with ID
+    c_ID\[I\] = Ith component of vector or Ith column of array calculated by a compute with ID
+    f_ID = scalar or vector calculated by a fix with ID
+    f_ID\[I\] = Ith component of vector or Ith column of array calculated by a fix with ID
+    v_name = value(s) calculated by an equal-style or atom-style variable with name :pre    
 :ule
 
 [Examples:]
 
 fix 1 all ave/histo 100 5 1000 0.5 1.5 50 c_myTemp file temp.histo ave running
 fix 1 all ave/histo 100 5 1000 -5 5 100 c_thermo_press\[2\] c_thermo_press\[3\] title1 "My output values"
-fix 1 all ave/histo 1 100 1000 -2.0 2.0 18 vx vy vz mode vector ave running beyond extra :pre
+fix 1 all ave/histo 1 100 1000 -2.0 2.0 18 vx vy vz mode vector ave running beyond extra 
+fix 1 all ave/histo 1 1 1 10 100 2000 c_XRD[1] mode vector weights c_XRD[2] :pre
 
 [Description:]
 
@@ -272,6 +279,15 @@ describes the six values that are printed at the first of each section
 of output.  The third describes the 4 values printed for each bin in
 the histogram.
 
+If the {weights} keyword is specified, the fix will compute a weighted
+histogram using per-bin weights specified by the {weights} argument.  As
+normal, the bin locations will be will be generated based off value1.
+However, instead of each binned value contributing 1 to the bin
+location, the contributing ammount is assigned the weights
+argument. Only a single value1 and weights argument pair can be can be
+used for each fix ave/histo.  The length of value1 must match the
+lenght of the weights arguemnt.
+
 :line
 
 [Restart, fix_modify, output, run start/stop, minimize info:]
diff --git a/doc/fix_saed_vtk.html b/doc/fix_saed_vtk.html
new file mode 100644
index 0000000000..55c15e8893
--- /dev/null
+++ b/doc/fix_saed_vtk.html
@@ -0,0 +1,191 @@
+<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>fix saed/vtk command 
+</H3>
+<P><B>Syntax:</B>
+</P>
+<PRE>fix ID group-ID saed/vtk Nevery Nrepeat Nfreak c_ID attribute args ... keyword args ... 
+</PRE>
+<UL><LI>ID, group-ID are documented in <A HREF = "fix.html">fix</A> command 
+
+<LI>ave/time/saed = style name of this fix command 
+
+<LI>Nevery = use input values every this many timesteps 
+
+<LI>Nrepeat = # of times to use input values for calculating averages 
+
+<LI>Nfreq = calculate averages every this many timesteps 
+
+<LI>c_ID = saed compute ID 
+
+<PRE>keyword = <I>file</I> or <I>ave</I> or <I>start</I> or <I>file</I> or <I>overwrite</I>:l
+  <I>ave</I> args = <I>one</I> or <I>running</I> or <I>window M</I>
+    one = output a new average value every Nfreq steps
+    running = output cumulative average of all previous Nfreq steps
+    window M = output average of M most recent Nfreq steps
+  <I>start</I> args = Nstart
+    Nstart = start averaging on this timestep
+  <I>file</I> arg = filename
+    filename = name of file to output time averages to 
+  <I>overwrite</I> arg = none = overwrite output file with only latest output 
+</PRE>
+
+</UL>
+<P><B>Examples:</B>
+</P>
+<P>compute 1 all saed 0.0251 Al O Kmax 1.70 Zone 0 0 1 dR_Ewald 0.01 c 0.5 0.5 0.5
+compute 2 all saed 0.0251 Ni Kmax 1.70 Zone 0 0 0 c 0.05 0.05 0.05 manual echo
+</P>
+<PRE>fix saed/vtk 1 1 1 c_1 file Al2O3_001.saed
+fix saed/vtk 1 1 1 c_2 file Ni_000.saed 
+</PRE>
+<P><B>Description:</B>
+</P>
+<P>Time average computed intensities from <A HREF = "compute_saed.txt">compute_saed</A> and 
+write output to a file in the 3rd generation vtk image data format for 
+visualization directly in parallelized visualization software packages 
+like ParaView and VisIt. Note that if no time averaging is done, this 
+command can be used as a convenient way to simply output diffraction 
+intensities at a single snapshot.
+</P>
+<P>To produce output in the image data vtk format ghost data is added
+outside the <I>Kmax</I> range assigned in the compute saed. The ghost data is 
+assigned a value of -1 and can be removed setting a minimum isovolume 
+of 0 within the vizualiziton software. SAED images can be created by 
+visualizing a spherical slice of the data that is centered at 
+R_Ewald*<B>h k l</B>/norm(<B>h k l</B>), where R_Ewald=1/lambda.
+</P>
+<P>The group specified within this command is ignored. However, note that 
+specified values may represent calculations performed by saed computes
+which store their own "group" definitions.
+</P>
+<P>Fix saed/vtk is designed to work only with <A HREF = "compute_saed.txt">compute_saed</A> 
+values.
+</P>
+<P>compute 3 top saed 0.0251 Al O 
+fix saed/vtk 1 1 1 c_3 file Al2O3_001.saed
+</P>
+<HR>
+
+<P>The <I>Nevery</I>, <I>Nrepeat</I>, and <I>Nfreq</I> arguments specify on what
+timesteps the input values will be used in order to contribute to the
+average.  The final averaged quantities are generated on timesteps
+that are a multiple of <I>Nfreq</I>.  The average is over <I>Nrepeat</I>
+quantities, computed in the preceding portion of the simulation every
+<I>Nevery</I> timesteps.  <I>Nfreq</I> must be a multiple of <I>Nevery</I> and
+<I>Nevery</I> must be non-zero even if <I>Nrepeat</I> is 1.  Also, the timesteps
+contributing to the average value cannot overlap, i.e. Nfreq >
+(Nrepeat-1)*Nevery is required.
+</P>
+<P>For example, if Nevery=2, Nrepeat=6, and Nfreq=100, then values on
+timesteps 90,92,94,96,98,100 will be used to compute the final average
+on timestep 100.  Similarly for timesteps 190,192,194,196,198,200 on
+timestep 200, etc.  If Nrepeat=1 and Nfreq = 100, then no time
+averaging is done; values are simply generated on timesteps
+100,200,etc.
+</P>
+<HR>
+
+<P>The output for fix ave/time/saed is a file writen with the 3rd generation
+vtk image data formatting.  The filename assigned by the <I>file</I> keyword is 
+appended with _N.vtk where N is an index (0,1,2...) to account for multiple 
+diffraction intensity outputs.
+</P>
+<P>By default the header contains the following information (with example data):
+</P>
+<P># vtk DataFile Version 3.0 c_SAED
+Image data set
+ASCII
+DATASET STRUCTURED_POINTS
+DIMENSIONS 337 219 209
+ASPECT_RATIO 0.00507953 0.00785161 0.00821458
+ORIGIN -0.853361 -0.855826 -0.854316 \n POINT_DATA 15424827
+SCALARS intensity float
+LOOKUP_TABLE default
+...data
+</P>
+<HR>
+
+<P>Additional optional keywords also affect the operation of this fix.
+</P>
+<P>The <I>ave</I> keyword determines how the values produced every <I>Nfreq</I>
+steps are averaged with values produced on previous steps that were
+multiples of <I>Nfreq</I>, before they are accessed by another output
+command or written to a file.
+</P>
+<P>If the <I>ave</I> setting is <I>one</I>, then the values produced on timesteps
+that are multiples of <I>Nfreq</I> are independent of each other; they are
+output as-is without further averaging.
+</P>
+<P>If the <I>ave</I> setting is <I>running</I>, then the values produced on
+timesteps that are multiples of <I>Nfreq</I> are summed and averaged in a
+cumulative sense before being output.  Each output value is thus the
+average of the value produced on that timestep with all preceding
+values.  This running average begins when the fix is defined; it can
+only be restarted by deleting the fix via the <A HREF = "unfix.html">unfix</A>
+command, or by re-defining the fix by re-specifying it.
+</P>
+<P>If the <I>ave</I> setting is <I>window</I>, then the values produced on
+timesteps that are multiples of <I>Nfreq</I> are summed and averaged within
+a moving "window" of time, so that the last M values are used to
+produce the output.  E.g. if M = 3 and Nfreq = 1000, then the output
+on step 10000 will be the average of the individual values on steps
+8000,9000,10000.  Outputs on early steps will average over less than M
+values if they are not available.
+</P>
+<P>The <I>start</I> keyword specifies what timestep averaging will begin on.
+The default is step 0.  Often input values can be 0.0 at time 0, so
+setting <I>start</I> to a larger value can avoid including a 0.0 in a
+running or windowed average.
+</P>
+<P>The <I>file</I> keyword allows a filename to be specified.  Every <I>Nfreq</I>
+steps, the vector of saed intensity data is written to a new file using 
+the 3rd generation vtk format.  The base of each file is assigned by 
+the <I>file</I> keyword and this string is appended with _N.vtk where N is 
+an index (0,1,2...) to account for situations with multiple diffraction
+intensity outputs.
+</P>
+<P>The <I>overwrite</I> keyword will continuously overwrite the output file
+with the latest output, so that it only contains one timestep worth of
+output.  This option can only be used with the <I>ave running</I> setting.
+</P>
+<P><B>Restart, fix_modify, output, run start/stop, minimize info:</B>
+</P>
+<P>No information about this fix is written to <A HREF = "restart.html">binary restart
+files</A>.  None of the <A HREF = "fix_modify.html">fix_modify</A> options
+are relevant to this fix.
+</P>
+<P>No parameter of this fix can be used with the <I>start/stop</I> keywords of
+the <A HREF = "run.html">run</A> command.  This fix is not invoked during <A HREF = "minimize.html">energy
+minimization</A>.
+</P>
+<P><B>Restrictions:</B> 
+</P>
+<P>The attributes for fix_saed_vtk must match the values assigned in the
+associated <A HREF = "compute_saed.txt">compute_saed</A> command.
+</P>
+<P><B>Related commands:</B> 
+</P>
+<P><A HREF = "compute_saed.html">compute_saed</A>
+</P>
+<P><B>Default:</B> 
+</P>
+<P>The option defaults are ave = one, start = 0, no file output.
+</P>
+<HR>
+
+<A NAME = "Coleman"></A>
+
+<P><B>(Coleman)</B> Coleman, Spearot, Capolungo, MSMSE, 21, 055020
+(2013).
+</P>
+</HTML>
diff --git a/doc/fix_saed_vtk.txt b/doc/fix_saed_vtk.txt
new file mode 100644
index 0000000000..67c940b192
--- /dev/null
+++ b/doc/fix_saed_vtk.txt
@@ -0,0 +1,181 @@
+"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
+
+fix saed/vtk command :h3
+
+[Syntax:]
+
+fix ID group-ID saed/vtk Nevery Nrepeat Nfreak c_ID attribute args ... keyword args ... :pre
+
+ID, group-ID are documented in "fix"_fix.html command :ulb,l
+ave/time/saed = style name of this fix command :l
+Nevery = use input values every this many timesteps :l
+Nrepeat = # of times to use input values for calculating averages :l
+Nfreq = calculate averages every this many timesteps :l
+c_ID = saed compute ID :l
+
+keyword = {file} or {ave} or {start} or {file} or {overwrite}:l
+  {ave} args = {one} or {running} or {window M}
+    one = output a new average value every Nfreq steps
+    running = output cumulative average of all previous Nfreq steps
+    window M = output average of M most recent Nfreq steps
+  {start} args = Nstart
+    Nstart = start averaging on this timestep
+  {file} arg = filename
+    filename = name of file to output time averages to 
+  {overwrite} arg = none = overwrite output file with only latest output :pre
+
+:ule
+
+[Examples:]
+
+compute 1 all saed 0.0251 Al O Kmax 1.70 Zone 0 0 1 dR_Ewald 0.01 c 0.5 0.5 0.5
+compute 2 all saed 0.0251 Ni Kmax 1.70 Zone 0 0 0 c 0.05 0.05 0.05 manual echo
+
+fix saed/vtk 1 1 1 c_1 file Al2O3_001.saed
+fix saed/vtk 1 1 1 c_2 file Ni_000.saed :pre
+
+[Description:]
+
+Time average computed intensities from "compute_saed"_compute_saed.txt and 
+write output to a file in the 3rd generation vtk image data format for 
+visualization directly in parallelized visualization software packages 
+like ParaView and VisIt. Note that if no time averaging is done, this 
+command can be used as a convenient way to simply output diffraction 
+intensities at a single snapshot.
+
+To produce output in the image data vtk format ghost data is added
+outside the {Kmax} range assigned in the compute saed. The ghost data is 
+assigned a value of -1 and can be removed setting a minimum isovolume 
+of 0 within the vizualiziton software. SAED images can be created by 
+visualizing a spherical slice of the data that is centered at 
+R_Ewald*[h k l]/norm([h k l]), where R_Ewald=1/lambda.
+
+The group specified within this command is ignored. However, note that 
+specified values may represent calculations performed by saed computes
+which store their own "group" definitions.
+
+Fix saed/vtk is designed to work only with "compute_saed"_compute_saed.txt 
+values.
+
+compute 3 top saed 0.0251 Al O 
+fix saed/vtk 1 1 1 c_3 file Al2O3_001.saed
+
+:line
+
+The {Nevery}, {Nrepeat}, and {Nfreq} arguments specify on what
+timesteps the input values will be used in order to contribute to the
+average.  The final averaged quantities are generated on timesteps
+that are a multiple of {Nfreq}.  The average is over {Nrepeat}
+quantities, computed in the preceding portion of the simulation every
+{Nevery} timesteps.  {Nfreq} must be a multiple of {Nevery} and
+{Nevery} must be non-zero even if {Nrepeat} is 1.  Also, the timesteps
+contributing to the average value cannot overlap, i.e. Nfreq >
+(Nrepeat-1)*Nevery is required.
+
+For example, if Nevery=2, Nrepeat=6, and Nfreq=100, then values on
+timesteps 90,92,94,96,98,100 will be used to compute the final average
+on timestep 100.  Similarly for timesteps 190,192,194,196,198,200 on
+timestep 200, etc.  If Nrepeat=1 and Nfreq = 100, then no time
+averaging is done; values are simply generated on timesteps
+100,200,etc.
+
+:line
+
+The output for fix ave/time/saed is a file writen with the 3rd generation
+vtk image data formatting.  The filename assigned by the {file} keyword is 
+appended with _N.vtk where N is an index (0,1,2...) to account for multiple 
+diffraction intensity outputs.
+
+By default the header contains the following information (with example data):
+
+# vtk DataFile Version 3.0 c_SAED
+Image data set
+ASCII
+DATASET STRUCTURED_POINTS
+DIMENSIONS 337 219 209
+ASPECT_RATIO 0.00507953 0.00785161 0.00821458
+ORIGIN -0.853361 -0.855826 -0.854316 \n POINT_DATA 15424827
+SCALARS intensity float
+LOOKUP_TABLE default
+...data
+
+:line
+
+Additional optional keywords also affect the operation of this fix.
+
+The {ave} keyword determines how the values produced every {Nfreq}
+steps are averaged with values produced on previous steps that were
+multiples of {Nfreq}, before they are accessed by another output
+command or written to a file.
+
+If the {ave} setting is {one}, then the values produced on timesteps
+that are multiples of {Nfreq} are independent of each other; they are
+output as-is without further averaging.
+
+If the {ave} setting is {running}, then the values produced on
+timesteps that are multiples of {Nfreq} are summed and averaged in a
+cumulative sense before being output.  Each output value is thus the
+average of the value produced on that timestep with all preceding
+values.  This running average begins when the fix is defined; it can
+only be restarted by deleting the fix via the "unfix"_unfix.html
+command, or by re-defining the fix by re-specifying it.
+
+If the {ave} setting is {window}, then the values produced on
+timesteps that are multiples of {Nfreq} are summed and averaged within
+a moving "window" of time, so that the last M values are used to
+produce the output.  E.g. if M = 3 and Nfreq = 1000, then the output
+on step 10000 will be the average of the individual values on steps
+8000,9000,10000.  Outputs on early steps will average over less than M
+values if they are not available.
+
+The {start} keyword specifies what timestep averaging will begin on.
+The default is step 0.  Often input values can be 0.0 at time 0, so
+setting {start} to a larger value can avoid including a 0.0 in a
+running or windowed average.
+
+The {file} keyword allows a filename to be specified.  Every {Nfreq}
+steps, the vector of saed intensity data is written to a new file using 
+the 3rd generation vtk format.  The base of each file is assigned by 
+the {file} keyword and this string is appended with _N.vtk where N is 
+an index (0,1,2...) to account for situations with multiple diffraction
+intensity outputs.
+
+The {overwrite} keyword will continuously overwrite the output file
+with the latest output, so that it only contains one timestep worth of
+output.  This option can only be used with the {ave running} setting.
+
+[Restart, fix_modify, output, run start/stop, minimize info:]
+
+No information about this fix is written to "binary restart
+files"_restart.html.  None of the "fix_modify"_fix_modify.html options
+are relevant to this fix.
+
+No parameter of this fix can be used with the {start/stop} keywords of
+the "run"_run.html command.  This fix is not invoked during "energy
+minimization"_minimize.html.
+
+[Restrictions:] 
+
+The attributes for fix_saed_vtk must match the values assigned in the
+associated "compute_saed"_compute_saed.txt command.
+
+[Related commands:] 
+
+"compute_saed"_compute_saed.html
+
+[Default:] 
+
+The option defaults are ave = one, start = 0, no file output.
+
+:line
+
+:link(Coleman)
+[(Coleman)] Coleman, Spearot, Capolungo, MSMSE, 21, 055020
+(2013).
+
diff --git a/src/MPIIO/dump_atom_mpiio.cpp b/src/MPIIO/dump_atom_mpiio.cpp
index ca8c4f542b..8ac03829ee 100644
--- a/src/MPIIO/dump_atom_mpiio.cpp
+++ b/src/MPIIO/dump_atom_mpiio.cpp
@@ -211,8 +211,8 @@ void DumpAtomMPIIO::init_style()
     strcat(format,"\n");
   } else {
     char *str;
-    if (image_flag == 0) str = (char *) "%d %d %g %g %g";
-    else str = (char *) "%d %d %g %g %g %d %d %d";
+    if (image_flag == 0) str = (char *) TAGINT_FORMAT " %d %g %g %g";
+    else str = (char *) TAGINT_FORMAT " %d %g %g %g %d %d %d";
     int n = strlen(str) + 2;
     format = new char[n];
     strcpy(format,str);
diff --git a/src/SRD/fix_wall_srd.cpp b/src/SRD/fix_wall_srd.cpp
index 4c9655a278..d192b55c4a 100644
--- a/src/SRD/fix_wall_srd.cpp
+++ b/src/SRD/fix_wall_srd.cpp
@@ -60,7 +60,7 @@ FixWallSRD::FixWallSRD(LAMMPS *lmp, int narg, char **arg) :
       else if (strcmp(arg[iarg],"zlo") == 0) newwall = ZLO;
       else if (strcmp(arg[iarg],"zhi") == 0) newwall = ZHI;
 
-      for (int m = 0; m < nwall; m++)
+      for (int m = 0; (m < nwall) && (m < 6); m++)
         if (newwall == wallwhich[m])
           error->all(FLERR,"Wall defined twice in fix wall/srd command");
 
diff --git a/src/STUBS/mpi.c b/src/STUBS/mpi.c
index 987479bb7b..fc2d60429a 100644
--- a/src/STUBS/mpi.c
+++ b/src/STUBS/mpi.c
@@ -361,7 +361,7 @@ int MPI_Cart_rank(MPI_Comm comm, int *coords, int *rank)
 int MPI_Type_contiguous(int count, MPI_Datatype oldtype, 
                         MPI_Datatype *newtype)
 {
-  if (nextra_datatype = MAXEXTRA_DATATYPE) return -1;
+  if (nextra_datatype == MAXEXTRA_DATATYPE) return -1;
   ptr_datatype[nextra_datatype] = newtype;
   index_datatype[nextra_datatype] = -(nextra_datatype + 1);
   size_datatype[nextra_datatype] = count * stubtypesize(oldtype);
diff --git a/src/USER-ATC/fix_atc.cpp b/src/USER-ATC/fix_atc.cpp
index ebf2b96088..b13719add3 100644
--- a/src/USER-ATC/fix_atc.cpp
+++ b/src/USER-ATC/fix_atc.cpp
@@ -43,6 +43,10 @@ using namespace LAMMPS_NS;
 using namespace FixConst;
 using std::string;
 
+#ifdef LAMMPS_BIGBIG
+#error "The USER-ATC package is not compatible with -DLAMMPS_BIGBIG"
+#endif
+
 // main page of doxygen documentation
 /*! \mainpage AtC : Atom-to-Continuum methods
     fix commands:
diff --git a/src/USER-AWPMD/fix_nve_awpmd.cpp b/src/USER-AWPMD/fix_nve_awpmd.cpp
index 377d42b37e..c4a3820e6b 100644
--- a/src/USER-AWPMD/fix_nve_awpmd.cpp
+++ b/src/USER-AWPMD/fix_nve_awpmd.cpp
@@ -88,12 +88,8 @@ void FixNVEAwpmd::initial_integrate(int vflag)
   double *erforce = atom->erforce;
   double *vforce=atom->vforce;
   double *ervelforce=atom->ervelforce;
-  double *cs=atom->cs;
-  double *csforce=atom->csforce;
-
 
   double *mass = atom->mass;
-  int *spin = atom->spin;
   int *type = atom->type;
   int *mask = atom->mask;
   int nlocal = atom->nlocal;
diff --git a/src/USER-AWPMD/pair_awpmd_cut.cpp b/src/USER-AWPMD/pair_awpmd_cut.cpp
index 9813c15d59..e9ca0352f2 100644
--- a/src/USER-AWPMD/pair_awpmd_cut.cpp
+++ b/src/USER-AWPMD/pair_awpmd_cut.cpp
@@ -74,8 +74,8 @@ PairAWPMDCut::~PairAWPMDCut()
 
 
 struct cmp_x{
-  double tol;
   double **xx;
+  double tol;
   cmp_x(double **xx_=NULL, double tol_=1e-12):xx(xx_),tol(tol_){}
   bool operator()(const pair<int,int> &left, const pair<int,int> &right) const {
     if(left.first==right.first){
@@ -116,8 +116,6 @@ void PairAWPMDCut::compute(int eflag, int vflag)
   double **x = atom->x;
   double **f = atom->f;
   double *q = atom->q;
-  double *erforce = atom->erforce;
-  double *eradius = atom->eradius;
   int *spin = atom->spin;
   int *type = atom->type;
   int *etag = atom->etag;
@@ -128,7 +126,6 @@ void PairAWPMDCut::compute(int eflag, int vflag)
   int ntot=nlocal+nghost;
 
   int newton_pair = force->newton_pair;
-  double qqrd2e = force->qqrd2e;
 
   int inum = list->inum;
   int *ilist = list->ilist;
diff --git a/src/USER-COLVARS/colvarproxy_lammps.h b/src/USER-COLVARS/colvarproxy_lammps.h
index 74f7193784..316d953f42 100644
--- a/src/USER-COLVARS/colvarproxy_lammps.h
+++ b/src/USER-COLVARS/colvarproxy_lammps.h
@@ -29,7 +29,7 @@ inline std::ostream & operator<< (std::ostream &out, const commdata &cd)
   out << " (" << cd.tag << "/" << cd.type << ": "
       << cd.x << ", " << cd.y << ", " << cd.z << ") ";
   return out;
-};
+}
 
 /// \brief Communication between colvars and LAMMPS
 /// (implementation of \link colvarproxy \endlink)
diff --git a/src/USER-MISC/pair_tersoff_table.cpp b/src/USER-MISC/pair_tersoff_table.cpp
index 69c5057cea..77002a5ba3 100644
--- a/src/USER-MISC/pair_tersoff_table.cpp
+++ b/src/USER-MISC/pair_tersoff_table.cpp
@@ -557,7 +557,7 @@ void PairTersoffTable::deallocateGrids()
 
 void PairTersoffTable::allocateGrids(void)
 {
-  int   i, j, l;
+  int   i, j, k, l;
 
   int     numGridPointsExponential, numGridPointsGtetaFunction, numGridPointsOneCutoffFunction;
   int     numGridPointsNotOneCutoffFunction, numGridPointsCutoffFunction, numGridPointsBetaZetaPower;
@@ -634,9 +634,9 @@ void PairTersoffTable::allocateGrids(void)
     zeta_max = MAX(zeta_max,numGridPointsBetaZetaPower);
 
     for (j=0; j<nelements; j++) {
-      for (j=0; j<nelements; j++) {
+      for (k=0; k<nelements; k++) {
 
-        int ijparam = elem2param[i][j][j];
+        int ijparam = elem2param[i][j][k];
         double cutoffR = params[ijparam].cutoffR;
         double cutoffS = params[ijparam].cutoffS;
 
diff --git a/src/compute_cna_atom.cpp b/src/compute_cna_atom.cpp
index 2bf33fec3d..69425457e8 100644
--- a/src/compute_cna_atom.cpp
+++ b/src/compute_cna_atom.cpp
@@ -262,7 +262,7 @@ void ComputeCNAAtom::compute_peratom()
         firstflag = 1;
         ncommon = 0;
         for (inear = 0; inear < nnearest[i]; inear++)
-          for (jnear = 0; jnear < n; jnear++)
+          for (jnear = 0; (jnear < n) && (n < MAXNEAR); jnear++)
             if (nearest[i][inear] == onenearest[jnear]) {
               if (ncommon < MAXCOMMON) common[ncommon++] = nearest[i][inear];
               else if (firstflag) {
-- 
GitLab