This method may not work for some cases since your question is not detailed enough (just in my opinion :) )
I am selecting 2 vertices (1 from the Cube and 1 from the Plane) that I assumed to be alligned by looking from your sight (as you see in the blender scene ) and getting their world coorinates. After that I am taking the difference of them to slide Cube. (Also I am assuming our objects are not rotated after created.)

import bpy, bmesh
from mathutils import Vector
bpy.ops.object.select_all(action='DESELECT')
# set active object to Cube
bpy.context.scene.objects.active = bpy.data.objects["Cube"]
# go to EDIT mode
bpy.ops.object.editmode_toggle()
# get Cube's updated vertex data
obj = bpy.data.objects["Cube"]
bm = bmesh.new()
bm = bmesh.from_edit_mesh(obj.data)
bm.verts.ensure_lookup_table()
# get world coordinate matrix for that object
mat = mat = obj.matrix_world
# get the right-bottom-closer vertex's world coordinates and save it to v1. By your sight , this vertex must have minimum x, y and z coordinates among the others.
# Additional info : if it is that vertex, its x coordinate is lower than upper adjacent , its y coordinate is lower than left adjacent ,its z coordinate is lower than further adjacent so I am searching these in the if and elif statements
for vert in bm.verts:
minx = False
miny = False
minz = False
for edge in vert.link_edges:
vert2 = edge.other_vert(vert)
if (mat*vert.co)[0]<(mat*vert2.co)[0]:
minx = True
elif (mat*vert.co)[1]<(mat*vert2.co)[1]:
miny = True
elif (mat*vert.co)[2]<(mat*vert2.co)[2]:
minz = True
if minx == True and miny == True and minz == True:
v1 = mat*vert.co
bm.free()
# switch object mode and set active object to the plane to get its vertex data
bpy.ops.object.editmode_toggle()
bpy.ops.object.select_all(action='DESELECT')
bpy.context.scene.objects.active = bpy.data.objects["Plane"]
bpy.ops.object.editmode_toggle()
obj = bpy.data.objects["Plane"]
bm = bmesh.new()
bm = bmesh.from_edit_mesh(obj.data)
bm.verts.ensure_lookup_table()
mat = mat = obj.matrix_world
# get the left-closer vertex's world coordinates and save it to v2. By your sight , this vertex must have minimum x and maximum y coordinates among the others.
# Additional info : if it is that vertex, its x coordinate is lower than upper adjacent , its y coordinate is higher than left adjacent ,its x coordinate is lower than further adjacent so I am searching these in the if and elif statements
for vert in bm.verts:
minx = False
maxy = False
for edge in vert.link_edges:
vert2 = edge.other_vert(vert)
if (mat*vert.co)[0]<(mat*vert2.co)[0]:
minx = True
elif (mat*vert.co)[1]>(mat*vert2.co)[1]:
maxy = True
if minx == True and maxy == True :
v2 = mat*vert.co
bm.free()
# switch to object mode to move all cube at once , since we got every information we need
bpy.ops.object.editmode_toggle()
# move cube by their difference vector, so it will slide just by their difference value
bpy.data.objects["Cube"].location+= (v2-v1)
EDIT:
I changed my code according to specifications now it gets 4 indexes (2 from cube(or glass), 2 from plane(or face) is enough to allign them)
import bpy, bmesh
from mathutils import Vector
import numpy as np , math
index_1_from_cube = 4
index_2_from_cube = 6
index_1_from_plane = 2
index_2_from_plane = 3
bpy.ops.object.select_all(action='DESELECT')
# set active object to Cube
bpy.context.scene.objects.active = bpy.data.objects["Cube"]
# go to EDIT mode
bpy.ops.object.editmode_toggle()
# get Cube's updated vertex data
obj = bpy.data.objects["Cube"]
bm = bmesh.new()
bm = bmesh.from_edit_mesh(obj.data)
bm.verts.ensure_lookup_table()
# get world coordinate matrix for that object
mat = mat = obj.matrix_world
# get world coordinates of 2 vertices from cube that is planned to allign
v1 = mat*bm.verts[index_1_from_cube].co
v1_for_rotate = mat*bm.verts[index_2_from_cube].co
bm.free()
# switch object mode and set active object to the plane to get its vertex data
bpy.ops.object.editmode_toggle()
bpy.ops.object.select_all(action='DESELECT')
bpy.context.scene.objects.active = bpy.data.objects["Plane"]
bpy.ops.object.editmode_toggle()
obj = bpy.data.objects["Plane"]
bm = bmesh.new()
bm = bmesh.from_edit_mesh(obj.data)
bm.verts.ensure_lookup_table()
mat = mat = obj.matrix_world
# get world coordinates of 2 vertices from plane that is planned to allign
v2 = mat*bm.verts[index_1_from_plane].co
v2_for_rotate = mat*bm.verts[index_2_from_plane].co
bm.free()
# switch to object mode to move all cube at once , since we got every information we need
bpy.ops.object.editmode_toggle()
# move cube by their difference vector, so it will slide just by their difference value
# but after this operation only the 2 vertices will be alligned,
# so we need to rotate it according to second set of points.
bpy.data.objects["Cube"].location+= (v2-v1)
#select the cube again
bpy.context.scene.objects.active = bpy.data.objects["Cube"]
bpy.data.objects["Cube"].select = True
# we will rotate the cube around the first point since the first points are already alligned (we don't need to change this coordinate of cube)
bpy.context.scene.cursor_location = v2
bpy.ops.object.origin_set(type='ORIGIN_CURSOR')
# get unit vectors in order to find a angle and rotation axis
unit_vec1 = (v1_for_rotate - v1 ) / np.linalg.norm(v1_for_rotate - v1 )
unit_vec2 = (v2_for_rotate - v2 ) / np.linalg.norm(v2_for_rotate - v2)
rotation_axis = (np.cross(unit_vec1, unit_vec2))
rotation_degree = np.arccos(np.clip(np.dot(unit_vec1, unit_vec2), -1.0, 1.0))
# rotate the object around first alligned points
bpy.ops.transform.rotate(value=rotation_degree, axis=rotation_axis, constraint_axis=(False, False, False))