informal procedure:

I took a 1 tmao molec .gro file and did
  $ genbox-lionx_dp -ci tmao.1.gro -cs spc216.gro -o tmao.3.3.3.gro -box 3 3 3 -nmol 17
or
  $ genbox-lionx_dp -ci tmao.1.gro -cs spc216.gro -o tmao.4.4.4.gro -box 4 4 4 -nmol 39
depending on whether I was building a 3 nm or 4 nm edge box. for the pure water systems, I just did
  $ genbox-lionx_dp -cs spc216.gro -o water.#.#.#.gro -box # # #
where # is either 3 or 4.

From there, I ran an energy minimization if necessary
 -- deemed necessary if the next step didn't work
Then, I annealed the system in the NVT ensemble by going from
time (ns): 0    9    11    20
Temp (K) : 298  1000 1000  298
After that, I took the resulting configuration and ran a 5 ns NPT equilibration
From this 5 ns equilibration, I found the average volume, and extracted a frame close to this average volume
With these frames, I added 40 nm on to the z-dimension of the box and ran them for 100ns in the NVT ensemble
-- 3x3x3 boxes actually ran for : 5 ns eq + 50 ns prod
-- 4x4x4 boxes actually ran for : 5 ns eq + 25 ns prod

more formal procedure:

For the surface tension calculations, systems were initially constructed in cubic boxes of edge length 3 or 
4 nm. Energy minimization was performed (if necessary) using the steepest descent algorithm. The resulting 
configuration was then heated in the NVT ensemble using the stochastic dynamics integrator from 298K to 1000K 
over a 9ns window, simulated at 1000K for 2ns, and then cooled down to 298K over a 9ns window. The resultant 
configuration was then simulated in the NPT ensemble using the stochastic dynamics integrator and the 
Parrinello-Rahman barostat for 5 ns at an external pressure of 1 bar. From this 5 ns simulation, the average
volume was calculated, and a configuration was extracted that had a volume as close to the average as possible.
The z-dimensions of these configurations were extended an additional 40nm, and the resulting configurations 
simulated for ??? ns in the NVT ensemble using the stochastic dynamics integrator. All simulations used a 2 fs 
timestep, held the hydrogen bonds constant using the LINCS algorithm, used a 2 ps coupling constant for the 
thermostat, cutoff the LJ and short ranged electrostatic interactions at 1.3 nm, and used the PME method for 
long ranged electrostatic interactions. The NPT simulations used a coupling constant of 0.5 ps and a 
compressibility of 4.5e-05 bar^-1, i.e. that of pure water.

17 TMAO :  804 water 3x3x3
39 TMAO : 1967 water 4x4x4

Molarity calculations:

17 TMAO * 1 mol / 6.022E23 / 27 nm^3 * 10^24 nm^3/L = 1.046 mol/L
39 TMAO * 1 mol / 6.022E23 / 64 nm^3 * 10^24 nm^3/L = 1.012 mol/L

Molality calculations:

17 TMAO * 1 mol / 6.022E23 / 804 waters * 1 water / 18.0154 amu * 6.022E23 amu / g * 1000g / kg = 1.174 molal
39 TMAO * 1 mol / 6.022E23 / 1967 waters * 1 water / 18.0154 amu * 6.022E23 amu / g * 1000g / kg = 1.101 molal 

Surface Tension Calculations:

If I understand the GROMACS 4.5.3 bug with the surface tension calculations correctly, then g_energy SHOULD
give the correct average surface tension for all of my simulations because I did not set a value for 
nstcalcenergy. the default is -1, which sets it equal to nstlist, which I set as 1. So, since nstcalcenergy
is = 1, it not greater than the lowest common denominator of nstcomm, nsttcouple, and/or nstpcouple.

Furthermore, the value given by g_energy is an average over all steps, not just every 500th step. For this reason,
it is more accurate than the average obtained by opening the .xvg file and using xmgrace to find the average. Also,
these two answers will obviously be different. SO. All of the surface tensions I report will be the values obtained
from 
  g_energy -f XXX.edr -o ST.xvg -b 5000 -e #### 
where #### depends on how many nanoseconds I want to average over, which will depend on the system in question and
how far the simulation made it before it got killed after 4 days.