How to do a relative rotation of a cuboid such that the longest edge is upright?

I am using pybullet.

I have a cuboid with half extents [4,2,1] corresponding to x,y,z dimensions. I rotate this cuboid around the z world axis theta and also translate it by [0, 0, .65].

I take a pointcloud of the cuboid, giving me an $N times 3$ vector where $N$ is the number of points, and 3 is the xyz dimensionality. Let’s call it $P$.

How do I apply a second rotation to this cuboid such that the longest edge becomes aligned with the z-world axis?

I have been doing:
$P’ = P @ begin{bmatrix}0 & 0 & 1 & 0 \ 0 & 1 & 0 & 0 \ 1 & 0 & 0 & 0 \ 0 & 0& 0 & 1end{bmatrix}$

(where @ is the matrix multiply operator in torch/numpy)

But this is giving me totally wrong results: https://i.snipboard.io/4qbup0.jpg

(Red is the pointcloud after the second rotation is applied, black is the pointcloud after the first rotation is applied. The longest edge of the red one should be aligned with the z-axis.)