Rotate object's mesh so it matches another object's rotation
If your objects are in Euler rotation mode, then the quaternion wont be set. For any mode using

ob.matrix_world.to_euler()
ob.matrix_world.to_quaternion()


will give the global euler and quaternion rotation respectively.

import bpy
import mathutils
from mathutils import  Quaternion, Vector, Matrix
import bmesh
import math
from math import radians

context = bpy.context
data = bpy.data

obj1 = data.objects['Cube']
obj2 = data.objects['Cube.001']

from_quat_rot = obj1.matrix_world.to_quaternion()
to_quat_rot = obj2.matrix_world.to_quaternion()
from_euler_rot = obj1.rotation_euler
to_euler_rot = obj2.rotation_euler

difference_quat = from_quat_rot.rotation_difference(to_quat_rot)
difference_angle = difference_quat.angle

print(f"difference angle: {difference_angle}")

vector_zero = Vector((0, 0, 0))
bm = bmesh.new()
bm.from_mesh(obj1.data)
bmesh.ops.rotate(bm, 
        cent=vector_zero, 
        verts=bm.verts, 
        matrix=Matrix.Rotation(difference_angle, 4, difference_quat.axis))
bm.to_mesh(obj1.data)
bm.free()


This rotates the mesh as expected.  Running again will have an accumulative effect.

could also

M = difference_quat.to_matrix()
obj1.data.transform(M)
obj1.data.update()