One way to go about it, using answer from https://blender.stackexchange.com/a/53976/15543.
A cube has six sides, with orthogonal set of normals, all with equal area.
Below is code to split all faces into "sides" with same face normal and calculates the sides area and volume. The volume should be zero, if not the faces aren't coplanar and don't make one side of a cube.
Run the test code in edit mode.
If end up with 6 sides, equal area, all coplanar with orthogonal (either dot product is zero or sum length is zero ( v == -v) )
from math import radians
from mathutils import Vector
context = bpy.context
def coplanar(faces, TOL=0.00001):
bm2 = bmesh.new()
verts = [v.co for f in faces for v in f.verts]
o = sum(verts, Vector((0,0,0))) / len(verts)
for v in verts:
# origin to center (matters for convex_hull)
bm2.verts.new(v - o)
#calculate a convex hull
# if the convex hull has no volume its on same plane.
volume = bm2.calc_volume()
# create a bmesh
TOL = 0.01
bm = bmesh.from_edit_mesh(context.edit_object.data)
# tried (unsuccesfully) to calculate orthonormal set from face normal
f = bm.faces
normal = f.normal
v = Vector((1, 0, 0))
if normal.dot(v) < 0.000001: # parallel
v = Vector((0, 1, 0))
#get the axis
x = normal
y = x.cross(v).normalized()
z = x.cross(y).normalized()
print("-" * 72)
faces = bm.faces[:]
sides = 
f = faces
axis = f.normal
# has to be at least one to match f.normal
side_faces = [f for f in faces if f.normal.angle(axis) < TOL]
area = sum(f.calc_area() for f in side_faces)
volume = coplanar(side_faces)
for f in side_faces:
sides.append((axis, area, volume))
print("Axis:", axis, "Area %5.2f Volume %5.2f" % (area, volume))
Won't be a cube if
len(sides) != 6 and so on.
Better if an orthogonal axis set was calculated from the first face normal, but couldn't work that out.