code: def quaternion_rmsd(P, Q):
      
      rot = quaternion_rotate(P, Q)
      P = np.dot(P, rot)
      return rmsd(P, Q)
docs: Rotate matrix P unto Q and calculate the RMSD
      based on doi:10.1016/1049-9660(91)90036-O
  
      Parameters
      ----------
      P : array
          (N,D) matrix, where N is points and D is dimension.
      Q : array
          (N,D) matrix, where N is points and D is dimension.
  
      Returns
      -------
      rmsd : float