Define $u := \frac{1}{d} (x_2 - x_1, y_2 - y_1, z_2 - z_1)^T$, in which $d$ is the distance between your two points. This means, $u$ is the unit vector in the direction of the line segment. You are looking for a rotation that maps $e_x = (1, 0, 0)$ to $u$.
The solution is not unique, as there are (infinitely) many rotations that achieve this, but you can single out the one that goes the most "direct" route - which is the one rotating about the axis perpendicular to both $u$ and $e_x$.
If you define $v := u \times e_x$ and $w := v \times u$ (using the cross product), the rotation matrix $R = (u, v, w)$ achieves exactly this (here, the vectors are the columns of $R$, which is a 3x3 matrix).
This matrix allows you to perform the desired rotation, and if all you need to do is rotate stuff, you should stick with it. However, if you really need to extract the rotation angles about the individual axes, it gets a bit ugly, bu see for example here for a way to do it.
As a first step, I'd move $P_1$ to the origin, so that the points become $P_1=(0,0,0)$, $P_2=(x_2-x_1, y_2-y_1, z_2-z_1)$, $P_3=(x_3-x_1, y_3-y_1, z_3-z_1)$. From there it's just a question of applying a rotation matrix.
A rotation matrix $\mathbf{R}$ which preserves all distances and shapes etc is
given by
$$
\mathbf{R}=\left(
\begin{array}
[c]{ccc}%
e_{1}^{x} & e_{2}^{x} & e_{3}^{x}\\
e_{1}^{y} & e_{2}^{y} & e_{3}^{y}\\
e_{1}^{z} & e_{2}^{z} & e_{3}^{z}%
\end{array}
\right)
$$
where the vectors $\left( e_{1}^{\alpha},e_{2}^{\alpha},e_{3}^{\alpha
}\right) $ for $\alpha\in\left\{ x,y,z\right\} $ are orthogonal and of unit
length, i.e.
$$
\sum_{i=1}^{3}e_{i}^{\alpha}e_{i}^{\beta}=\left\{
\begin{array}
[c]{ccc}%
1 & & \alpha=\beta\\
0 & & \alpha\neq\beta
\end{array}
\right. .%
$$
Applying $\mathbf{R}$, your rotated points $P_{i}^{\prime}=\left(
x_{i}^{\prime},y_{i}^{\prime},z_{i}^{\prime}\right) $ are given by
$$
\left(
\begin{array}
[c]{c}%
x_{i}^{\prime}\\
y_{i}^{\prime}\\
z_{i}^{\prime}%
\end{array}
\right) =\mathbf{R}\left(
\begin{array}
[c]{c}%
x_{i}\\
y_{i}\\
z_{i}%
\end{array}
\right)
$$
I suggest making the top row of $\mathbf{R}$ the unit
vector in the direction from $P_{1}$ to $P_{2}$, i.e.
$$
\left(
\begin{array}
[c]{c}%
e_{1}^{x}\\
e_{2}^{x}\\
e_{3}^{x}%
\end{array}
\right) =\frac{1}{\sqrt{\left( x_{1}-x_{2}\right) ^{2}+\left( y_{1}%
-y_{2}\right) ^{2}+\left( z_{1}-z_{2}\right) ^{2}}}\left(
\begin{array}
[c]{c}%
x_{1}-x_{2}\\
y_{1}-y_{2}\\
z_{1}-z_{2}%
\end{array}
\right)
$$
which would mean that the line from $P_{1}$ to $P_{2}$ gets mapped into the
$x$-axis.
Working out what you could use for the second and third rows of $\mathbf{R}$ is now just a question of solving some simple linear equations, to make sure that the line from $P_{1}$ to $P_{3}$ has no z component.
To solve for $\left( e_{1}^{z},e_{2}^{z},e_{3}^{z}\right)$, the vector
$e_{i}^{z}$ must have a zero dot-product with both $e_{i}^{x}$ and $\left(
x_{3}-x_{1},y_{3}-y_{1},z_{3}-z_{1}\right) $, so the equations are
$$
\begin{align*}
e_{1}^{x}e_{1}^{z}+e_{2}^{x}e_{2}^{z}+e_{3}^{x}e_{3}^{z} & =0\\
\left( x_{3}-x_{1}\right) e_{1}^{z}+\left( y_{3}-y_{1}\right) e_{2}%
^{z}+\left( z_{3}-z_{1}\right) e_{3}^{z} & =0
\end{align*}
$$
Hence
$$
\left(
\begin{array}
[c]{c}%
e_{1}^{z}\\
e_{2}^{z}\\
e_{3}^{z}%
\end{array}
\right) =\lambda_{z}\left(
\begin{array}
[c]{c}%
\left( z_{3}-z_{1}\right) e_{2}^{x}-\left( y_{3}-y_{1}\right) e_{3}^{x}\\
\left( x_{3}-x_{1}\right) e_{3}^{x}-\left( z_{3}-z_{1}\right) e_{1}^{x}\\
\left( y_{3}-y_{1}\right) e_{1}^{x}-\left( x_{3}-x_{1}\right) e_{2}^{x}%
\end{array}
\right)
$$
for some $\lambda_z$, which should be determined so that $e_{i}^{z}$ is a vector
of unit length. Finally $e_{i}^{y}$ can be determined as the vector which is
perpendicular to both $e_{i}^{x}$ and $e_{i}^{z}$, i.e.
$$
\left(
\begin{array}
[c]{c}%
e_{1}^{y}\\
e_{2}^{y}\\
e_{3}^{y}%
\end{array}
\right) =\lambda_{y}\left(
\begin{array}
[c]{c}%
e_{3}^{z}e_{2}^{x}-e_{2}^{z}e_{3}^{x}\\
e_{1}^{z}e_{3}^{x}-e_{3}^{z}e_{1}^{x}\\
e_{2}^{z}e_{1}^{x}-e_{1}^{z}e_{2}^{x}%
\end{array}
\right)
$$
where again $\lambda_{y}$ is determined so that $e_{i}^{y}$ is a vector of
unit length.
Best Answer
Given the point $A$ you shall imagine to draw a parallelepiped with a vertex on it, the opposite vertex in $O$ and with the faces parallel to the coordinate planes, as sketched.
Then a rotation around the axis $z$ will move any point in the space keeping fixed the plane orthogonal to z in which they lay , i.e. the $z$ coordinate, along a circle lying in that plane and centered on the $z$ axis, for an angle measured in that plane.
So $A=(x_1, y_1,z_1)$ will move to $A'=(x_2, y_2,z_1)$ (note z_1 is unchanged) , along a circle centered at $P$, with radius $r=\sqrt{{x_1}^2+{y_1}^2}$, and angle $\theta$ measured in that plane $z=z_1$.
$S$ will move in that same plane, with the same angle, but on a circle of radius $x_1$.
$Q$ will do the same, on a circle with the same radius, but on the plane $z=0$ along a circle centered at $O$, with an angle $\theta$ in that plane.
That is to say that the whole parallelepiped $OPSQA$ will rotate keeping fixed the edge $OP$ and thus rotate to $OPS'Q'A'$, with $$S'=(x_1 \cos \theta, x_1 \sin \theta, z_1), \; Q'=(x_1 \cos \theta, x_1 \sin \theta, 0)$$ while the expression for $A'$ is a little more complicated.