The surface tension is commonly calculated from MD simulations as
gamma = L_z/2 * <P_zz - 1/2 * (P_xx + P_zz) >
You can easily get that quantity from the md_trajectory via
# Get the pressures tensors of all snapshots.
# Possibly discard the first snapshots for equilibration.
p_tensors = md_trajectory.pressureTensors()
# Strip the units to be able to use numpy.mean().
p_tensors = p_tensors.inUnitsOf(eV/Ang**3)
# Get the height L_z
cell_height = md_trajectory.cells()[0, 2, 2]
# Average the pressure components over the snapshots.
surface_tension = numpy.mean(p_tensors[:, 2, 2] - 0.5*(p_tensors[:, 0, 0] + p_tensors[:, 1, 1]))
# Apply the units again and normalize with respect to the two surfaces and the height.
surface_tension = surface_tension*eV/Ang**3*cell_height/2.0
print surface_tension