You mean how the Python code works, what the different statements mean?
# Calculate the mean
v_z = numpy.apply_over_axes(numpy.mean,potential[:,:,:],[0,1]).flatten()
v_z *= potential[:,:,:].unit()
# Print out the result
dX, dY, dZ = potential.volumeElement()
dz = dZ.norm()
shape = potential.shape()
for i in range(shape[2]):
print dz*i, v_z[i]
Line 2 does almost all the work. potential[:,:,:] could also be written potential.toArray() which is easier to read (but the method was not implemented when the example was written
). Line 3 can also be written slightly simpler (see below).
Then the numpy "mean" operator is applied to the X and Y directions (axis 0 and 1).
The result is then "flattened" to a 1D array, viz. the averaged potential in the Z direction. The flattening is needed because the whole operation returns the vector inside another vector (twice, since we average over 2 directions).
The rest of the code is just to put the unit back on and compute the positions in the Z direction. It could be made a bit simpler, since the averaging really only makes sense for the case when XY is perpendicular to Z. Here's another version of the code that does the same thing (in that case, at least), which also makes the units a bit more obvious (the original example prints in Bohr and Hartree):
# Calculate the mean
v_z = numpy.apply_over_axes(numpy.mean,potential.toArray(),[0,1]).flatten()
v_z *= potential.unit()
# Print out the result
dz = potential.volumeElement()[2][2].inUnitsOf(Angstrom)
print "Z / Ang potential / eV"
for i,v in enumerate(v_z.inUnitsOf(eV)):
print dz*i, v
Which of the construction for the final "print" part you prefer is more a matter of taste