I am trying to understand virtual sites in MD simulations, and I came across this configuration:
Here, coordinate $\mathbf{s}$ represents the virtual site, which is formed by three other atoms $\mathbf{i}$, $\mathbf{j}$, and $\mathbf{k}$. The distance between atom $\mathbf{i}$ and the virtual site $\mathbf{s}$ is $|\mathbf{d}|$. The position of atoms $\mathbf{i}$, $\mathbf{j}$, and $\mathbf{k}$ are $\mathbf{r}_i$, $\mathbf{r}_j$, and $\mathbf{r}_k$, respectively.
In this case, the virtual site ($\mathbf{r}_s$) is in the plane of the other three particles at a distance of $|\mathbf{d}|$ from $\mathbf{i}$ at an angle of $\theta$ from $\mathbf{r}_{ij}$. Atom $\mathbf{k}$ defines the plane and direction of the angle.
How should I get the position of $\mathbf{r}_s$ using the other three atom positions and an angle?
I know the equation for $\mathbf{r}_s$, but I couldn't understand how it is derived. Any help is appreciated.
$$
\mathbf{r}_s = \mathbf{r}_i
+ d \cos\theta\, \frac{\mathbf{r}_{ij}}{|\mathbf{r}_{ij}|}
+ d \sin\theta\, \frac{\mathbf{r}_{\perp}}{|\mathbf{r}_{\perp}|}
$$
where
$$
\mathbf{r}_{\perp} = \mathbf{r}_{jk}
- \frac{\mathbf{r}_{ij} \cdot \mathbf{r}_{jk}}{\mathbf{r}_{ij} \cdot \mathbf{r}_{ij}}\, \mathbf{r}_{ij}
$$
Suppose that $d=|\mathbf{d}|$ and $\mathbf{r}_{ij}\ne \mathbf{0}$. Then in the plane of the other three particles with the coordinates origin at $\mathbf{r}_i$ and the basis of consisting of the unit vector $\frac{\mathbf{r}_{ij}}{|\mathbf{r}_{ij}|}$ and some orthogonal to it unit vector $\mathbf{r'}$, the virtual site has coordinates $(d\cos\theta,d\sin\theta)$. Thus $$\mathbf{r}_s = \mathbf{r}_i + d \cos\theta\, \frac{\mathbf{r}_{ij}}{|\mathbf{r}_{ij}|} + d \sin\theta\, \mathbf{r'}.$$
The vector $\mathbf{r'}$ belongs to the plane spanned by $\mathbf{r}_{ij}$ and $\mathbf{r}_{jk}$, so $\mathbf{r'}=\lambda_i\mathbf{r}_{ij}+\lambda_k\mathbf{r}_{jk}$ for some real $\lambda_i$ and $\lambda_k$. Moreover, $0=\mathbf{r'}\cdot \mathbf{r}_{ij}=\lambda_i\mathbf{r}_{ij}\cdot \mathbf{r}_{ij}+\lambda_k\mathbf{r}_{jk}\cdot \mathbf{r}_{ij}$. So $\lambda_i=-\lambda_k\frac{\mathbf{r}_{jk}\cdot \mathbf{r}_{ij}}{\mathbf{r}_{ij}\cdot \mathbf{r}_{ij}}$. Thus $$\mathbf{r'}=\lambda_k\left(-\frac{\mathbf{r}_{jk}\cdot \mathbf{r}_{ij}}{\mathbf{r}_{ij}\cdot \mathbf{r}_{ij}}\mathbf{r}_{ij} +\mathbf{r}_{jk}\right)=-\lambda_k\mathbf{r}_{\perp}.$$
Since $|\mathbf{r'}|=1$, we have $|\lambda_k|=\frac{1}{|\mathbf{r}_{\perp}|}$. It remains to determine the sign of $\lambda_k$. It corresponds to the direction of the angle $\theta$ between $\mathbf{r}_{ij}$ and $\mathbf{r}_{s}-\mathbf{r}_{i}$. In particular, when $\theta=\frac{\pi}{2}$, we have $\mathbf{r'}=\mathbf{r}_s-\mathbf{r}_i$. But this direction is not very clear to me, maybe, because there is no coordinate system in the supplementary picture.