How to choose a proper rotation matrix and apply it to rotate unit cell about Cartesian axis to align I-Cr-I bonds parallel to Cartesian axis

Question

I was surveying this paper Remarkably enhanced Curie temperature in monolayer CrI3 by hydrogen and oxygen adsorption: A first-principles calculations, where they have rotated the unit cell (without changing lattice parameters) using a rotation matrix

$$ \begin{bmatrix} -0.8 & 0 & 0.6 \\ 0.4 & -0.7 & 0.6 \\ 0.4 & 0.7 & 0.6 \end{bmatrix} $$

to align $I-Cr-I$ bonds parallel to Cartesian axis in octahedrons.

My question is: How to calculate this rotation matrix?

In this paper they have mentioned "For this we first find the proper rotation matrix using the VESTA program. Then we calculated the corresponding rotation angles and rotated the $\ce{CrI3}$ unit cell using Atomsk." But I didn't understand how to find proper rotation matrix from VESTA. So how to do this? Moreover, I have written a code regarding this by analyzing POSCAR file (here of $\ce{CrI3}$), but unable to extract this reported rotation matrix, please do comment where I'm making mistake?

The code is,

import numpy as np

def fractional_to_cartesian(frac_coords, lattice):
    return np.dot(frac_coords, lattice)

def rotation_matrix_from_vectors(vec1, vec2):
    """ Find the rotation matrix that aligns vec1 to vec2 """
    a, b = (vec1 / np.linalg.norm(vec1)).reshape(3), (vec2 / np.linalg.norm(vec2)).reshape(3)
    v = np.cross(a, b)
    c = np.dot(a, b)
    s = np.linalg.norm(v)
    kmat = np.array([[0, -v[2], v[1]], [v[2], 0, -v[0]], [-v[1], v[0], 0]])
    rotation_matrix = np.eye(3) + kmat + kmat.dot(kmat) * ((1 - c) / (s ** 2))
    return rotation_matrix


with open('POSCAR', 'r') as file:
    poscar_contents = file.readlines()


lattice_vectors = np.array([list(map(float, poscar_contents[i].split())) for i in range(2, 5)])


num_atoms_per_type = list(map(int, poscar_contents[6].split()))
num_atoms = sum(num_atoms_per_type)

atom_positions = np.array([list(map(float, line.split())) for line in poscar_contents[8:8 + num_atoms]])


cartesian_positions = np.array([fractional_to_cartesian(pos, lattice_vectors) for pos in atom_positions])

# Indices for the chosen octahedron around Cr(0) #CrI3 structure
cr_index = 0
i_indices = [14, 12, 8, 20, 18, 6] # I atoms in an octahedral

# Created vectors for the bonds we need to align
v_x1 = cartesian_positions[8] - cartesian_positions[cr_index]
v_x2 = cartesian_positions[cr_index] - cartesian_positions[20]
v_y1 = cartesian_positions[6] - cartesian_positions[cr_index]
v_y2 = cartesian_positions[cr_index] - cartesian_positions[18]
v_z1 = cartesian_positions[14] - cartesian_positions[cr_index]
v_z2 = cartesian_positions[cr_index] - cartesian_positions[12]

# Define the target vectors along Cartesian axes
target_x = np.array([1, 0, 0])
target_y = np.array([0, 1, 0])
target_z = np.array([0, 0, 1])

# Calculate the average vectors for x, y, z directions
avg_v_x = (v_x1 + v_x2) / 2
avg_v_y = (v_y1 + v_y2) / 2
avg_v_z = (v_z1 + v_z2) / 2

# Normalize the average vectors
avg_v_x /= np.linalg.norm(avg_v_x)
avg_v_y /= np.linalg.norm(avg_v_y)
avg_v_z /= np.linalg.norm(avg_v_z)


R_x = rotation_matrix_from_vectors(avg_v_x, target_x)
R_y = rotation_matrix_from_vectors(avg_v_y, target_y)
R_z = rotation_matrix_from_vectors(avg_v_z, target_z)

print(" R_x:")

print(R_x)
print("R_y:")
print(R_y)
print("R_z:")
print(R_z)

R = R_z @ R_y @ R_x

print("Rotation Matrix R:")
print(R)


rotated_cartesian_positions = np.dot(cartesian_positions, R)

inv_lattice_vectors = np.linalg.inv(lattice_vectors)
rotated_fractional_positions = np.dot(rotated_cartesian_positions, inv_lattice_vectors)


new_poscar_contents = poscar_contents[:2]
new_lattice_vectors = np.dot(lattice_vectors, R)
new_poscar_contents += ['{:22.16f} {:22.16f} {:22.16f}\n'.format(*vec) for vec in new_lattice_vectors]
new_poscar_contents += poscar_contents[5:8]
new_poscar_contents += ['{:22.16f} {:22.16f} {:22.16f}\n'.format(*pos) for pos in rotated_fractional_positions]


with open('POSCAR_rotated', 'w') as file:
    file.writelines(new_poscar_contents)

print("New POSCAR file written to POSCAR_rotated")