diff --git a/examples/USER/flow_gauss/README b/examples/USER/flow_gauss/README
deleted file mode 100644
index ef7cc82d961fc056b16fd127b8dda5dfb5b38220..0000000000000000000000000000000000000000
--- a/examples/USER/flow_gauss/README
+++ /dev/null
@@ -1,45 +0,0 @@
-The input script in.GD is an example simulation using Gaussian dynamics (GD).
-The simulation is of a simple 2d Lennard-Jones fluid flowing through a pipe.
-For details see online LAMMPS documentation and
-Strong and Eaves, J. Phys. Chem. Lett. 7(10) 2016, p. 1907.
-
-Note that the run times and box size are chosen to allow a fast example run.
-They are not adequate for a real simulation.
-
-The script has the following parts:
-1) initialize variables
-    These can be modified to customize the simulation. Note that if the
-    pipe dimensions L or d are changed, the geometry should be checked
-    by visualizing the coordinates in all.init.lammpstrj.
-
-2) create box
-
-3) set up potential
-
-4) create atoms
-
-5) set up profile-unbiased thermostat (PUT)
-    see Evans and Morriss, Phys. Rev. Lett. 56(20) 1986, p. 2172
-    By default, this uses boxes which contain on average 8 molecules.
-
-6) equilibrate without GD
-
-7) initialize the center-of-mass velocity and run to achieve steady-state
-    The system is initialized with a uniform velocity profile, which
-    relaxes over the course of the simulation.
-
-8) collect data
-    The data is output in several files:
-    GD.out contains the force that GD applies, and the flux in the x- and
-        y- directions. The output Jx should be equal to the value of
-        J set in section 1, which is 0.1 by default.
-    x_profiles contains the velocity, density, and pressure profiles in
-        the x-direction. The pressure profile is given by
-        (-1/2V)*(c_spa[1] + c_spa[2]), where V is the volume of a
-        slice. The pressure profile is computed with IK1, see
-        Todd, Evans, and Davis, Phys. Rev. E 52(2) 1995, p. 1627.
-        Note that to compare with the pump method, or to
-        compute a pressure drop, you must correct this pressure
-        profile as described in Strong 2016 above.
-    Vy_profile is the velocity profile inside the pipe along the
-        y-direction, u_x(y).
diff --git a/examples/USER/flow_gauss/in.GD b/examples/USER/flow_gauss/in.GD
deleted file mode 100644
index bcff4d4c57bba23b014b7c3dae942884be8c84a2..0000000000000000000000000000000000000000
--- a/examples/USER/flow_gauss/in.GD
+++ /dev/null
@@ -1,262 +0,0 @@
-#LAMMPS input script
-#in.GD
-#see README for details
-
-###############################################################################
-#initialize variables
-clear
-
-#frequency for outputting info (timesteps)
-variable        dump_rate       equal 50
-variable        thermo_rate     equal 10
-
-#equilibration time (timesteps)
-variable        equil           equal 1000
-
-#stabilization time (timesteps to reach steady-state)
-variable        stabil          equal 1000
-
-#data collection time (timesteps)
-variable        run             equal 2000
-
-#length of pipe
-variable        L               equal 30
-
-#width of pipe
-variable        d               equal 20
-
-#flux (mass/sigma*tau)
-variable        J               equal 0.1
-
-#simulation box dimensions
-variable        Lx              equal 100
-variable        Ly              equal 40
-
-#bulk fluid density
-variable        dens            equal 0.8
-
-#lattice spacing for wall atoms
-variable        aWall           equal 1.0 #1.7472
-
-#timestep
-variable        ts              equal 0.001
-
-#temperature
-variable        T               equal 2.0
-
-#thermostat damping constant
-variable        tdamp           equal ${ts}*100
-
-units           lj
-dimension       2
-atom_style      atomic
-
-
-###############################################################################
-#create box
-
-#create lattice with the spacing aWall
-variable        rhoWall         equal ${aWall}^(-2)
-lattice         sq ${rhoWall}
-
-#modify input dimensions to be multiples of aWall
-variable        L1 equal round($L/${aWall})*${aWall}
-variable        d1 equal round($d/${aWall})*${aWall}
-variable        Ly1 equal round(${Ly}/${aWall})*${aWall}
-variable        Lx1 equal round(${Lx}/${aWall})*${aWall}
-
-#create simulation box
-variable        lx2 equal ${Lx1}/2
-variable        ly2 equal ${Ly1}/2
-region          simbox block -${lx2} ${lx2} -${ly2} ${ly2} 0 0.1 units box
-create_box      2 simbox
-
-#####################################################################
-#set up potential
-
-mass            1 1.0           #fluid atoms
-mass            2 1.0           #wall atoms
-
-pair_style      lj/cut 2.5
-pair_modify     shift yes
-pair_coeff      1 1 1.0 1.0 2.5
-pair_coeff      1 2 1.0 1.0 1.12246
-pair_coeff      2 2 0.0 0.0
-
-neigh_modify  exclude type 2 2
-
-timestep        ${ts}
-
-#####################################################################
-#create atoms
-
-#create wall atoms everywhere
-create_atoms    2 box
-
-#define region which is "walled off"
-variable        dhalf equal ${d1}/2
-variable        Lhalf equal ${L1}/2
-region          walltop block -${Lhalf} ${Lhalf} ${dhalf} EDGE -0.1 0.1 &
-                units box
-region          wallbot block -${Lhalf} ${Lhalf} EDGE -${dhalf} -0.1 0.1 &
-                units box
-region          outsidewall union 2 walltop wallbot side out
-
-#remove wall atoms outside wall region
-group           outside region outsidewall
-delete_atoms    group outside
-
-#remove wall atoms that aren't on edge of wall region
-variable        x1 equal ${Lhalf}-${aWall}
-variable        y1 equal ${dhalf}+${aWall}
-region          insideTop block -${x1} ${x1} ${y1} EDGE -0.1 0.1 units box
-region          insideBot block -${x1} ${x1} EDGE -${y1} -0.1 0.1 units box
-region          insideWall union 2 insideTop insideBot
-group           insideWall region insideWall
-delete_atoms    group insideWall
-
-#define new lattice, to give correct fluid density
-#y lattice const must be a multiple of aWall
-variable        atrue equal ${dens}^(-1/2)
-variable        ay equal round(${atrue}/${aWall})*${aWall}
-
-#choose x lattice const to give correct density
-variable        ax equal (${ay}*${dens})^(-1)
-
-#change Lx to be multiple of ax
-variable        Lx1 equal round(${Lx}/${ax})*${ax}
-variable        lx2 equal ${Lx1}/2
-change_box      all x final -${lx2} ${lx2} units box
-
-#define new lattice
-lattice         custom ${dens} &
-                a1 ${ax} 0.0 0.0 a2 0.0 ${ay} 0.0 a3 0.0 0.0 1.0 &
-                basis 0.0 0.0 0.0
-
-#fill in rest of box with bulk particles
-variable        delta equal 0.001
-variable        Ldelt equal ${Lhalf}+${delta}
-variable        dDelt equal ${dhalf}-${delta}
-region          left block EDGE -${Ldelt} EDGE EDGE -0.1 0.1 units box
-region          right block ${Ldelt} EDGE EDGE EDGE -0.1 0.1 units box
-region          pipe block -${Ldelt} ${Ldelt} -${dDelt} ${dDelt} -0.1 0.1 &
-                units box
-
-region          bulk union 3 left pipe right
-create_atoms    1 region bulk
-
-group           bulk type 1
-group           wall type 2
-
-#remove atoms that are too close to wall
-delete_atoms    overlap 0.9 bulk wall
-
-neighbor        0.3 bin
-neigh_modify    delay 0 every 1 check yes
-neigh_modify    exclude group wall wall
-
-velocity        bulk create $T 78915 dist gaussian rot yes mom yes loop geom
-
-#####################################################################
-#set up PUT
-#see Evans and Morriss, Phys. Rev. Lett. 56(20) 1986, p. 2172
-
-#average number of particles per box, Evans and Morriss used 2.0
-variable        NperBox equal 8.0
-
-#calculate box sizes
-variable        boxSide equal sqrt(${NperBox}/${dens})
-variable        nX equal round(lx/${boxSide})
-variable        nY equal round(ly/${boxSide})
-variable        dX equal lx/${nX}
-variable        dY equal ly/${nY}
-
-#temperature of fluid (excluding wall)
-compute         myT bulk temp
-
-#profile-unbiased temperature of fluid
-compute         myTp bulk temp/profile 1 1 0 xy ${nX} ${nY}
-
-#thermo setup
-thermo          ${thermo_rate}
-thermo_style    custom step c_myT c_myTp etotal press
-
-#dump initial configuration
-# dump            55 all custom 1 all.init.lammpstrj id type x y z vx vy vz
-# dump            56 wall custom 1 wall.init.lammpstrj id type x y z
-# dump_modify     55 sort id
-# dump_modify     56 sort id
-run             0
-# undump          55
-# undump          56
-
-#####################################################################
-#equilibrate without GD
-
-fix             nvt bulk nvt temp $T $T ${tdamp}
-fix_modify      nvt temp myTp
-fix             2 bulk enforce2d
-
-run             ${equil}
-
-#####################################################################
-#initialize the COM velocity and run to achieve steady-state
-
-#calculate velocity to add: V=J/rho_total
-variable        Vadd            equal $J*lx*ly/count(bulk)
-
-#first remove any COM velocity, then add back the streaming velocity
-velocity        bulk zero linear
-velocity        bulk set ${Vadd} 0.0 0.0 units box sum yes mom no
-
-fix             GD bulk flow/gauss 1 0 0 #energy yes
-#fix_modify     GD energy yes
-
-run             ${stabil}
-
-#####################################################################
-#collect data
-
-#print the applied force and total flux to ensure conservation of Jx
-variable        Fapp equal f_GD[1]
-compute         vxBulk bulk reduce sum vx
-compute         vyBulk bulk reduce sum vy
-variable        invVol equal 1.0/(lx*ly)
-variable        jx equal c_vxBulk*${invVol}
-variable        jy equal c_vyBulk*${invVol}
-variable        curr_step equal step
-variable        p_Fapp format Fapp %.3f
-variable        p_jx format jx %.5g
-variable        p_jy format jy %.5g
-fix             print_vCOM all print ${dump_rate} &
-                "${curr_step} ${p_Fapp} ${p_jx} ${p_jy}" file GD.out screen no &
-                title "timestep Fapp Jx Jy"
-
-#compute IK1 pressure profile
-#see Todd, Evans, and Davis, Phys. Rev. E 52(2) 1995, p. 1627
-#use profile-unbiased temperature to remove the streaming velocity
-#from the kinetic part of the pressure
-compute         spa bulk stress/atom myTp
-
-#for the pressure profile, use the same grid as the PUT
-compute         chunkX bulk chunk/atom bin/1d x lower ${dX} units box
-
-#output pressure profile and other profiles
-#the pressure profile is (-1/2V)*(c_spa[1] + c_spa[2]), where
-#V is the volume of a slice
-fix             profiles bulk ave/chunk 1 1 ${dump_rate} chunkX &
-                vx density/mass  c_spa[1] c_spa[2] &
-                file x_profiles ave running overwrite
-
-#compute velocity profile across the pipe with a finer grid
-variable        dYnew equal ${dY}/10
-compute         chunkY bulk chunk/atom bin/1d y center ${dYnew} units box &
-                region pipe
-fix             velYprof bulk ave/chunk 1 1 ${dump_rate} chunkY &
-                vx file Vy_profile ave running overwrite
-
-#full trajectory
-# dump          7 bulk custom ${dump_rate} bulk.lammpstrj id type x y z
-# dump_modify     7 sort id
-
-run             ${run}
diff --git a/examples/USER/misc/flow_gauss/README b/examples/USER/misc/flow_gauss/README
index 4966cd2dc919eee2bd33f3970da1cb1a0d075eba..ef7cc82d961fc056b16fd127b8dda5dfb5b38220 100644
--- a/examples/USER/misc/flow_gauss/README
+++ b/examples/USER/misc/flow_gauss/README
@@ -1,45 +1,45 @@
 The input script in.GD is an example simulation using Gaussian dynamics (GD).
 The simulation is of a simple 2d Lennard-Jones fluid flowing through a pipe.
-For details see online LAMMPS documentation and 
+For details see online LAMMPS documentation and
 Strong and Eaves, J. Phys. Chem. Lett. 7(10) 2016, p. 1907.
 
-Note that the run times and box size are chosen to allow a fast example run. 
-They are not adequate for a real simulation. 
+Note that the run times and box size are chosen to allow a fast example run.
+They are not adequate for a real simulation.
 
 The script has the following parts:
 1) initialize variables
-	These can be modified to customize the simulation. Note that if the 
-	pipe dimensions L or d are changed, the geometry should be checked 
-	by visualizing the coordinates in all.init.lammpstrj. 
+    These can be modified to customize the simulation. Note that if the
+    pipe dimensions L or d are changed, the geometry should be checked
+    by visualizing the coordinates in all.init.lammpstrj.
 
 2) create box
-	
+
 3) set up potential
 
 4) create atoms
 
 5) set up profile-unbiased thermostat (PUT)
-	see Evans and Morriss, Phys. Rev. Lett. 56(20) 1986, p. 2172
-	By default, this uses boxes which contain on average 8 molecules.
+    see Evans and Morriss, Phys. Rev. Lett. 56(20) 1986, p. 2172
+    By default, this uses boxes which contain on average 8 molecules.
 
 6) equilibrate without GD
-	
+
 7) initialize the center-of-mass velocity and run to achieve steady-state
-	The system is initialized with a uniform velocity profile, which 
-	relaxes over the course of the simulation.
+    The system is initialized with a uniform velocity profile, which
+    relaxes over the course of the simulation.
 
 8) collect data
-	The data is output in several files:
-	GD.out contains the force that GD applies, and the flux in the x- and 
-		y- directions. The output Jx should be equal to the value of 
-		J set in section 1, which is 0.1 by default.
-	x_profiles contains the velocity, density, and pressure profiles in 
-		the x-direction. The pressure profile is given by 
-		(-1/2V)*(c_spa[1] + c_spa[2]), where V is the volume of a 
-		slice. The pressure profile is computed with IK1, see 
-		Todd, Evans, and Davis, Phys. Rev. E 52(2) 1995, p. 1627.
-		Note that to compare with the pump method, or to 
-		compute a pressure drop, you must correct this pressure 
-		profile as described in Strong 2016 above.
-	Vy_profile is the velocity profile inside the pipe along the 
-		y-direction, u_x(y).
+    The data is output in several files:
+    GD.out contains the force that GD applies, and the flux in the x- and
+        y- directions. The output Jx should be equal to the value of
+        J set in section 1, which is 0.1 by default.
+    x_profiles contains the velocity, density, and pressure profiles in
+        the x-direction. The pressure profile is given by
+        (-1/2V)*(c_spa[1] + c_spa[2]), where V is the volume of a
+        slice. The pressure profile is computed with IK1, see
+        Todd, Evans, and Davis, Phys. Rev. E 52(2) 1995, p. 1627.
+        Note that to compare with the pump method, or to
+        compute a pressure drop, you must correct this pressure
+        profile as described in Strong 2016 above.
+    Vy_profile is the velocity profile inside the pipe along the
+        y-direction, u_x(y).
diff --git a/examples/USER/misc/flow_gauss/in.GD b/examples/USER/misc/flow_gauss/in.GD
old mode 100755
new mode 100644
index 8117715c12705752c4efacaa04d6b562aef5e716..bcff4d4c57bba23b014b7c3dae942884be8c84a2
--- a/examples/USER/misc/flow_gauss/in.GD
+++ b/examples/USER/misc/flow_gauss/in.GD
@@ -7,83 +7,85 @@
 clear
 
 #frequency for outputting info (timesteps)
-variable	dump_rate 	equal 50
-variable	thermo_rate	equal 10
+variable        dump_rate       equal 50
+variable        thermo_rate     equal 10
 
 #equilibration time (timesteps)
-variable	equil		equal 1000
+variable        equil           equal 1000
 
 #stabilization time (timesteps to reach steady-state)
-variable	stabil		equal 1000
+variable        stabil          equal 1000
 
 #data collection time (timesteps)
-variable	run		equal 2000
+variable        run             equal 2000
 
 #length of pipe
-variable	L		equal 30
+variable        L               equal 30
 
 #width of pipe
-variable	d		equal 20
+variable        d               equal 20
 
 #flux (mass/sigma*tau)
-variable	J		equal 0.1
+variable        J               equal 0.1
 
 #simulation box dimensions
-variable	Lx		equal 100
-variable	Ly		equal 40
+variable        Lx              equal 100
+variable        Ly              equal 40
 
 #bulk fluid density
-variable	dens		equal 0.8
+variable        dens            equal 0.8
 
 #lattice spacing for wall atoms
-variable	aWall		equal 1.0 #1.7472
+variable        aWall           equal 1.0 #1.7472
 
 #timestep
-variable	ts		equal 0.001
+variable        ts              equal 0.001
 
 #temperature
-variable	T		equal 2.0
+variable        T               equal 2.0
 
 #thermostat damping constant
-variable	tdamp		equal ${ts}*100
+variable        tdamp           equal ${ts}*100
 
-units 		lj
-dimension	2
-atom_style	atomic
+units           lj
+dimension       2
+atom_style      atomic
 
 
 ###############################################################################
 #create box
 
 #create lattice with the spacing aWall
-variable	rhoWall		equal ${aWall}^(-2)
-lattice		sq ${rhoWall}
+variable        rhoWall         equal ${aWall}^(-2)
+lattice         sq ${rhoWall}
 
 #modify input dimensions to be multiples of aWall
-variable	L1 equal round($L/${aWall})*${aWall}
-variable	d1 equal round($d/${aWall})*${aWall}
-variable	Ly1 equal round(${Ly}/${aWall})*${aWall}
-variable	Lx1 equal round(${Lx}/${aWall})*${aWall}
+variable        L1 equal round($L/${aWall})*${aWall}
+variable        d1 equal round($d/${aWall})*${aWall}
+variable        Ly1 equal round(${Ly}/${aWall})*${aWall}
+variable        Lx1 equal round(${Lx}/${aWall})*${aWall}
 
 #create simulation box
-variable	lx2 equal ${Lx1}/2
-variable	ly2 equal ${Ly1}/2
-region		simbox block -${lx2} ${lx2} -${ly2} ${ly2} 0 0.1 units box
-create_box	2 simbox
+variable        lx2 equal ${Lx1}/2
+variable        ly2 equal ${Ly1}/2
+region          simbox block -${lx2} ${lx2} -${ly2} ${ly2} 0 0.1 units box
+create_box      2 simbox
 
 #####################################################################
 #set up potential
 
-mass 		1 1.0  		#fluid atoms
-mass 		2 1.0		#wall atoms
+mass            1 1.0           #fluid atoms
+mass            2 1.0           #wall atoms
 
-pair_style	lj/cut 2.5
-pair_modify	shift yes 
-pair_coeff	1 1 1.0 1.0 2.5
-pair_coeff	1 2 1.0 1.0 1.12246 
-pair_coeff	2 2 0.0 0.0 0.0
+pair_style      lj/cut 2.5
+pair_modify     shift yes
+pair_coeff      1 1 1.0 1.0 2.5
+pair_coeff      1 2 1.0 1.0 1.12246
+pair_coeff      2 2 0.0 0.0
 
-timestep 	${ts}
+neigh_modify  exclude type 2 2
+
+timestep        ${ts}
 
 #####################################################################
 #create atoms
@@ -92,167 +94,169 @@ timestep 	${ts}
 create_atoms    2 box
 
 #define region which is "walled off"
-variable 	dhalf equal ${d1}/2
-variable 	Lhalf equal ${L1}/2
-region		walltop block -${Lhalf} ${Lhalf} ${dhalf} EDGE -0.1 0.1 &
-		units box
-region		wallbot block -${Lhalf} ${Lhalf} EDGE -${dhalf} -0.1 0.1 &
-		units box
-region 		outsidewall union 2 walltop wallbot side out
+variable        dhalf equal ${d1}/2
+variable        Lhalf equal ${L1}/2
+region          walltop block -${Lhalf} ${Lhalf} ${dhalf} EDGE -0.1 0.1 &
+                units box
+region          wallbot block -${Lhalf} ${Lhalf} EDGE -${dhalf} -0.1 0.1 &
+                units box
+region          outsidewall union 2 walltop wallbot side out
 
 #remove wall atoms outside wall region
-group 		outside region outsidewall
-delete_atoms 	group outside
+group           outside region outsidewall
+delete_atoms    group outside
 
 #remove wall atoms that aren't on edge of wall region
-variable	x1 equal ${Lhalf}-${aWall}
-variable	y1 equal ${dhalf}+${aWall}
-region 		insideTop block -${x1} ${x1} ${y1} EDGE -0.1 0.1 units box
-region 		insideBot block -${x1} ${x1} EDGE -${y1} -0.1 0.1 units box
-region		insideWall union 2 insideTop insideBot
-group 		insideWall region insideWall
-delete_atoms	group insideWall
+variable        x1 equal ${Lhalf}-${aWall}
+variable        y1 equal ${dhalf}+${aWall}
+region          insideTop block -${x1} ${x1} ${y1} EDGE -0.1 0.1 units box
+region          insideBot block -${x1} ${x1} EDGE -${y1} -0.1 0.1 units box
+region          insideWall union 2 insideTop insideBot
+group           insideWall region insideWall
+delete_atoms    group insideWall
 
 #define new lattice, to give correct fluid density
 #y lattice const must be a multiple of aWall
-variable	atrue equal ${dens}^(-1/2)
-variable	ay equal round(${atrue}/${aWall})*${aWall}
+variable        atrue equal ${dens}^(-1/2)
+variable        ay equal round(${atrue}/${aWall})*${aWall}
 
 #choose x lattice const to give correct density
-variable	ax equal (${ay}*${dens})^(-1)
+variable        ax equal (${ay}*${dens})^(-1)
 
 #change Lx to be multiple of ax
-variable	Lx1 equal round(${Lx}/${ax})*${ax}
-variable	lx2 equal ${Lx1}/2
-change_box 	all x final -${lx2} ${lx2} units box
+variable        Lx1 equal round(${Lx}/${ax})*${ax}
+variable        lx2 equal ${Lx1}/2
+change_box      all x final -${lx2} ${lx2} units box
 
 #define new lattice
-lattice 	custom ${dens} &
-		a1 ${ax} 0.0 0.0 a2 0.0 ${ay} 0.0 a3 0.0 0.0 1.0 &
-		basis 0.0 0.0 0.0
+lattice         custom ${dens} &
+                a1 ${ax} 0.0 0.0 a2 0.0 ${ay} 0.0 a3 0.0 0.0 1.0 &
+                basis 0.0 0.0 0.0
 
 #fill in rest of box with bulk particles
-variable	delta equal 0.001
-variable	Ldelt equal ${Lhalf}+${delta}
-variable	dDelt equal ${dhalf}-${delta}
-region		left block EDGE -${Ldelt} EDGE EDGE -0.1 0.1 units box
-region 		right block ${Ldelt} EDGE EDGE EDGE -0.1 0.1 units box
-region		pipe block -${Ldelt} ${Ldelt} -${dDelt} ${dDelt} -0.1 0.1 &
-		units box
+variable        delta equal 0.001
+variable        Ldelt equal ${Lhalf}+${delta}
+variable        dDelt equal ${dhalf}-${delta}
+region          left block EDGE -${Ldelt} EDGE EDGE -0.1 0.1 units box
+region          right block ${Ldelt} EDGE EDGE EDGE -0.1 0.1 units box
+region          pipe block -${Ldelt} ${Ldelt} -${dDelt} ${dDelt} -0.1 0.1 &
+                units box
 
-region 		bulk union 3 left pipe right 
-create_atoms	1 region bulk
+region          bulk union 3 left pipe right
+create_atoms    1 region bulk
 
-group 		bulk type 1
-group		wall type 2
+group           bulk type 1
+group           wall type 2
 
 #remove atoms that are too close to wall
 delete_atoms    overlap 0.9 bulk wall
 
-neighbor	0.3 bin
-neigh_modify	delay 0 every 1 check yes
+neighbor        0.3 bin
+neigh_modify    delay 0 every 1 check yes
 neigh_modify    exclude group wall wall
 
-velocity 	bulk create $T 78915 dist gaussian rot yes mom yes loop geom
+velocity        bulk create $T 78915 dist gaussian rot yes mom yes loop geom
 
 #####################################################################
 #set up PUT
 #see Evans and Morriss, Phys. Rev. Lett. 56(20) 1986, p. 2172
 
 #average number of particles per box, Evans and Morriss used 2.0
-variable	NperBox equal 8.0
+variable        NperBox equal 8.0
 
 #calculate box sizes
-variable	boxSide equal sqrt(${NperBox}/${dens})
-variable	nX equal round(lx/${boxSide})
-variable	nY equal round(ly/${boxSide})
-variable	dX equal lx/${nX}
-variable	dY equal ly/${nY}
+variable        boxSide equal sqrt(${NperBox}/${dens})
+variable        nX equal round(lx/${boxSide})
+variable        nY equal round(ly/${boxSide})
+variable        dX equal lx/${nX}
+variable        dY equal ly/${nY}
 
 #temperature of fluid (excluding wall)
-compute 	myT bulk temp
+compute         myT bulk temp
 
 #profile-unbiased temperature of fluid
-compute 	myTp bulk temp/profile 1 1 0 xy ${nX} ${nY} 
+compute         myTp bulk temp/profile 1 1 0 xy ${nX} ${nY}
 
 #thermo setup
-thermo		${thermo_rate}	
-thermo_style 	custom step c_myT c_myTp etotal press 
+thermo          ${thermo_rate}
+thermo_style    custom step c_myT c_myTp etotal press
 
 #dump initial configuration
-dump 		55 all custom 1 all.init.lammpstrj id type x y z vx vy vz
-dump 		56 wall custom 1 wall.init.lammpstrj id type x y z
-dump_modify   	55 sort id
-dump_modify   	56 sort id
-run 		0
-undump 		55
-undump 		56
+# dump            55 all custom 1 all.init.lammpstrj id type x y z vx vy vz
+# dump            56 wall custom 1 wall.init.lammpstrj id type x y z
+# dump_modify     55 sort id
+# dump_modify     56 sort id
+run             0
+# undump          55
+# undump          56
 
 #####################################################################
 #equilibrate without GD
 
-fix 		nvt bulk nvt temp $T $T ${tdamp}
-fix_modify	nvt temp myTp
-fix		2 bulk enforce2d
+fix             nvt bulk nvt temp $T $T ${tdamp}
+fix_modify      nvt temp myTp
+fix             2 bulk enforce2d
 
-run		${equil}
+run             ${equil}
 
 #####################################################################
 #initialize the COM velocity and run to achieve steady-state
 
 #calculate velocity to add: V=J/rho_total
-variable	Vadd		equal $J*lx*ly/count(bulk)
+variable        Vadd            equal $J*lx*ly/count(bulk)
 
 #first remove any COM velocity, then add back the streaming velocity
 velocity        bulk zero linear
-velocity 	bulk set ${Vadd} 0.0 0.0 units box sum yes mom no
+velocity        bulk set ${Vadd} 0.0 0.0 units box sum yes mom no
 
-fix		GD bulk flow/gauss 1 0 0 #energy yes
-#fix_modify	GD energy yes
+fix             GD bulk flow/gauss 1 0 0 #energy yes
+#fix_modify     GD energy yes
 
-run		${stabil}
+run             ${stabil}
 
 #####################################################################
 #collect data
 
 #print the applied force and total flux to ensure conservation of Jx
-variable	Fapp equal f_GD[1]
-compute 	vxBulk bulk reduce sum vx
-compute 	vyBulk bulk reduce sum vy
+variable        Fapp equal f_GD[1]
+compute         vxBulk bulk reduce sum vx
+compute         vyBulk bulk reduce sum vy
 variable        invVol equal 1.0/(lx*ly)
-variable	jx equal c_vxBulk*${invVol}
-variable	jy equal c_vyBulk*${invVol}
-variable	curr_step equal step
-fix		print_vCOM all print ${dump_rate} &
-		"${curr_step} ${Fapp} ${jx} ${jy}" file GD.out screen no &
-		title "timestep Fapp Jx Jy"
-
-#compute IK1 pressure profile 
+variable        jx equal c_vxBulk*${invVol}
+variable        jy equal c_vyBulk*${invVol}
+variable        curr_step equal step
+variable        p_Fapp format Fapp %.3f
+variable        p_jx format jx %.5g
+variable        p_jy format jy %.5g
+fix             print_vCOM all print ${dump_rate} &
+                "${curr_step} ${p_Fapp} ${p_jx} ${p_jy}" file GD.out screen no &
+                title "timestep Fapp Jx Jy"
+
+#compute IK1 pressure profile
 #see Todd, Evans, and Davis, Phys. Rev. E 52(2) 1995, p. 1627
 #use profile-unbiased temperature to remove the streaming velocity
 #from the kinetic part of the pressure
-compute		spa bulk stress/atom myTp
+compute         spa bulk stress/atom myTp
 
 #for the pressure profile, use the same grid as the PUT
-compute		chunkX bulk chunk/atom bin/1d x lower ${dX} units box
+compute         chunkX bulk chunk/atom bin/1d x lower ${dX} units box
 
 #output pressure profile and other profiles
 #the pressure profile is (-1/2V)*(c_spa[1] + c_spa[2]), where
 #V is the volume of a slice
-fix		profiles bulk ave/chunk 1 1 ${dump_rate} chunkX &
-		vx density/mass  c_spa[1] c_spa[2] &
-		file x_profiles ave running overwrite
+fix             profiles bulk ave/chunk 1 1 ${dump_rate} chunkX &
+                vx density/mass  c_spa[1] c_spa[2] &
+                file x_profiles ave running overwrite
 
 #compute velocity profile across the pipe with a finer grid
-variable	dYnew equal ${dY}/10
-compute		chunkY bulk chunk/atom bin/1d y center ${dYnew} units box & 
-		region pipe
-fix		velYprof bulk ave/chunk 1 1 ${dump_rate} chunkY &
-		vx file Vy_profile ave running overwrite
+variable        dYnew equal ${dY}/10
+compute         chunkY bulk chunk/atom bin/1d y center ${dYnew} units box &
+                region pipe
+fix             velYprof bulk ave/chunk 1 1 ${dump_rate} chunkY &
+                vx file Vy_profile ave running overwrite
 
 #full trajectory
-dump 		7 bulk custom ${dump_rate} bulk.lammpstrj &
-		id type x y z 
-dump_modify	7 sort id
+# dump          7 bulk custom ${dump_rate} bulk.lammpstrj id type x y z
+# dump_modify     7 sort id
 
-run 		${run}
+run             ${run}
diff --git a/examples/USER/flow_gauss/log.6Jul17.GD.g++.1 b/examples/USER/misc/flow_gauss/log.6Jul17.GD.g++.1
similarity index 100%
rename from examples/USER/flow_gauss/log.6Jul17.GD.g++.1
rename to examples/USER/misc/flow_gauss/log.6Jul17.GD.g++.1
diff --git a/examples/USER/flow_gauss/log.6Jul17.GD.g++.4 b/examples/USER/misc/flow_gauss/log.6Jul17.GD.g++.4
similarity index 100%
rename from examples/USER/flow_gauss/log.6Jul17.GD.g++.4
rename to examples/USER/misc/flow_gauss/log.6Jul17.GD.g++.4
diff --git a/examples/USER/flow_gauss/output-files/GD.out b/examples/USER/misc/flow_gauss/output-files/GD.out
similarity index 100%
rename from examples/USER/flow_gauss/output-files/GD.out
rename to examples/USER/misc/flow_gauss/output-files/GD.out
diff --git a/examples/USER/flow_gauss/output-files/Vy_profile b/examples/USER/misc/flow_gauss/output-files/Vy_profile
similarity index 100%
rename from examples/USER/flow_gauss/output-files/Vy_profile
rename to examples/USER/misc/flow_gauss/output-files/Vy_profile
diff --git a/examples/USER/flow_gauss/output-files/x_profiles b/examples/USER/misc/flow_gauss/output-files/x_profiles
similarity index 100%
rename from examples/USER/flow_gauss/output-files/x_profiles
rename to examples/USER/misc/flow_gauss/output-files/x_profiles