Denote the two points as
$$
P_1:(x_1,y_1,z_1), \quad P_2: (x_2,y_2,z_2)
$$
and the center as $P_0: (x_0,y_0,z_0)$. First of all, note that the rotation matrix corresponding to the rotation about $P_0$ will take the vector $v_1 = P_1 - P_0$ to the vector $v_2 = P_2 - P_0$.
Note that a suitable rotation will only exist if $\|v_1\| = \|v_2\|$, i.e. both points are equidistant from the center of the rotation. If that is the case, then it suffices to find rotation matrix $R$ satisfying $Rv_1 = v_2$. Once this is done, the desired affine rotation is the same rotation about $P_0$, that is,
$$
T(v) = R(v - P_0) + P_0 = Rv + (I-R)P_0,
$$
corresponding to the homogeneous coordinates transformation matrix
$$
\pmatrix{R & (I-R)P_0\\0 & 1}.
$$
There are infinitely many choices for a rotation (about the center) from $v_1$ to $v_2$. However, I will suppose that we want the rotation by the smallest possible angle that achieves this result. I will also assume that $v_1$ and $v_2$ are linearly independent; the edge cases of $v_1 = \pm v_2$ can be handled separately.
I claim (without proof, for now) that the rotation that achieves this is the smaller of the two rotations along the axis orthogonal to both $v_1$ and $v_2$. In particular, it is the counterclockwise rotation by angle
$$
\theta = \arccos\left(\frac{v_1 \cdot v_2}{\|v_1\|\,\|v_2\|}\right)
$$
about the axis $v_1 \times v_2$. Let $u = (x_u,y_u,z_u)$ denote the unit vector parallel to $v_1 \times v_2$, $u = \frac{v_1\times v_2}{\|v_1 \times v_2\|}$. With Rodrigues' rotation formula, we find that this rotation has the matrix $R = I + \sin\theta K + (1-\cos \theta) K^2$, where $K$ denotes the "cross-product matrix" for the vector $u$,
$$
K = \pmatrix{
0 & -z_u & y_u\\
z_u&0&x_u\\
-y_u&x_u&0}.
$$
From there, the $z$-$x$-$Z$ extrinsic Euler angles for the rotation can be obtained using the formulas given here:
\begin{align}
\phi &= \operatorname{arctan2}\left(R_{31}, R_{32}\right),\\
\theta &= \arccos\left(R_{33}\right),\\
\psi &= -\operatorname{arctan2}\left(R_{13}, R_{23}\right).
\end{align}
As it turns out, we can neatly rewrite the formula for $R$ using the fact that $u = \frac 1{\sin\theta}v_1 \times v_2$. I claim (again, without proof for now) that if $v_1,v_2$ are normalized to have length $1$, then
$$
R = (v_1^\top v_2) I + (v_2v_1^\top - v_1v_2^\top) + (1 + v_1^\top v_2)(v_1 \times v_2)(v_1 \times v_2)^\top.
$$
Another nice reformulation:
$$
M = \pmatrix{v_1 & v_2 & v_1\times v_2}, \\
R = (v_1^\top v_2) I + M\pmatrix{0&-1&0\\1&0&0\\0&0&1 + v_1^\top v_2}M^\top.
$$
