I didn't manage to solve the problem (edit: in the meantime an answer was posted which says a precise formula using radicals cannot be found), but I can post a proof that the line joining the points where the minimal/maximal distance is achieved is perpendicular to the tangent line at the circles in those contact points. (inspired by the comment of Gerhard Paseman)
To do this, choose $\vec{a}$ and $\vec{d}$ the position vectors of the centers and $\vec{b},\vec{c}$, respectively $\vec{e},\vec{f}$ be two pairs of orthogonal unit vectors which span the planes of the first and respectively the second circle. Denote by $r,s$ the radii of the two circles. Consider the circles parametrized as (in fact, the argument works for any parametrization)
$$ p(\theta)=\vec{a}+r\cos\theta\ \vec{b}+r\sin\theta\ \vec{c}, \ \theta \in [0,2\pi] $$
$$ q(\tau)=\vec{d}+s\cos\tau\ \vec{e}+s\sin\tau\ \vec{f}, \tau \in [0,2\pi]$$
and denote $F(\theta,\tau)=|p(\theta)-q(\tau)|^2$. Then the pair of points which realize the minimal/maximal distance must satisfy
$$ \frac{\partial F}{\partial \theta}=\frac{\partial F}{\partial \tau}=0. $$
We have
$$ \frac{\partial F}{\partial \theta}=2\sum_{i=1}^3 [p_i(\theta)-q_i(\tau)]p_i'(\theta)=2 (p(\theta)-q(\tau))\cdot p'(\theta) $$
$$ \frac{\partial F}{\partial \tau}=-2\sum_{i=1}^3 [ p_i(\theta)-q_i(\tau) ] q_i'(\tau)=-2 (p(\theta)-q(\tau))\cdot q'(\tau) $$
where "$\cdot$" is the usual dot product. Therefore when $\theta,\tau$ correspond to the minimum/maximum value, the partial derivatives vanish and $p'(\theta)\perp (p(\theta)-q(\tau))$ and $q(\tau)'\perp (p(\theta)-q(\tau))$ where $p'(\theta),q'(\tau)$ are the tangent vectors in the contact points and $p(\theta)-q(\tau)$ is the vector connecting the points where minimal/maximal distance is achieved.