How can I multiply a quaternion and a list of vectors that represent xyz vertices of an object
Here's the answer that @batFINGER helped me figure out.

import bpy
import mathutils

ob = bpy.context.object

old_verts = [v.co for v in ob.data.vertices]
print("old_verts=%s" % old_verts)

#Do calculation 
quat1 = mathutils.Quaternion((1, 2, 3, 4))
new_verts = [quat1 @ v.co for v in ob.data.vertices]


print("--------------")
print("new_verts=%s" % new_verts)


Output:

old_verts=[Vector((1.0, 1.0, 1.0)), Vector((1.0, 1.0, -1.0)), Vector((1.0, -1.0, 1.0)), Vector((1.0, -1.0, -1.0)), Vector((-1.0, 1.0, 1.0)), Vector((-1.0, 1.0, -1.0)), Vector((-1.0, -1.0, 1.0)), Vector((-1.0, -1.0, -1.0))]
--------------
new_verts=[Vector((6.0, 30.0, 42.0)), Vector((-38.0, -10.0, 34.0)), Vector((-2.0, 50.0, -14.0)), Vector((-46.0, 10.0, -22.0)), Vector((46.0, -10.0, 22.0)), Vector((2.0, -50.0, 14.0)), Vector((38.0, 10.0, -34.0)), Vector((-6.0, -30.0, -42.0))]