Finding the 8 outer corner vertices of an object
KDTree and the objects bounding box.
Each object has an axis aligned bounding box giving the 8 coordinates of the cuboid the object fits into. If a cuboid is perfectly axis aligned the bounding box corner coordinates exactly match those of the corner verts.
A KDTree can be used to quickly find the closest element to a coordinate
Putting it together in a script, run with object of interest in edit mode (vertex select mode)
Note, this method assumes the mesh is axis aligned

Result of running on default UV sphere
import bmesh
from mathutils.kdtree import KDTree
from bpy import context

ob = context.edit_object
me = ob.data

bm = bmesh.from_edit_mesh(me)

kd = KDTree(len(bm.verts))

for v in bm.verts:
    kd.insert(v.co, v.index)

kd.balance()

# get closest to 8 bbox corners

closest = [kd.find(tuple(c))[1] for c in ob.bound_box]
print(closest)
#select 

for v in bm.verts:
    v.select = v.index in closest

bm.select_flush_mode()    
bmesh.update_edit_mesh(me)

Ordering the result
The bounding box is &quot;wound&quot; in a consistent order. Looking at the default cube
&gt;&gt;&gt; bbox = [tuple(b) for b in C.object.bound_box]
&gt;&gt;&gt; for i, v in enumerate(bbox):
...     i, v
...     
(0, (-1.0, -1.0, -1.0))
(1, (-1.0, -1.0, 1.0))
(2, (-1.0, 1.0, 1.0))
(3, (-1.0, 1.0, -1.0))
(4, (1.0, -1.0, -1.0))
(5, (1.0, -1.0, 1.0))
(6, (1.0, 1.0, 1.0))
(7, (1.0, 1.0, -1.0))

Hence the front y=-1, top (z=1) , left (x=1) vert is index 5. The 4 verts that make up the top corners are at indices 1, 2, 5, 6
The closest list has the same order.
If a definite 8 corner cuboid
In a cuboid mesh the 8 corner verts are only linked to 3 edges, whereas all others are linked to 4.  Can use this.
import bmesh
from bpy import context

ob = context.edit_object
me = ob.data

bm = bmesh.from_edit_mesh(me)

for v in bm.verts:
    v.select = (
            len(v.link_edges) == 3
            )

bm.select_flush_mode()    
bmesh.update_edit_mesh(me)

If this is also axis aligned then the &quot;order&quot; could be determined by looking at the components of the vertex normals.  For example, the normals of the corner verts of default cube would be equivalent to their normalized coordinates