# How is the equation for the position of a virtual site derived?

Cross posted on Math SE

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}|}\tag{1},$$ in which $$\mathbf{r}_{\perp} = \mathbf{r}_{jk} - \frac{\mathbf{r}_{ij} \cdot \mathbf{r}_{jk}}{\mathbf{r}_{ij} \cdot \mathbf{r}_{ij}}\, \mathbf{r}_{ij}. \tag{2}$$

• I think it's already mentioned in gromacs manual and the subsequent paper. Oct 5, 2023 at 20:22
• @RoshanShrestha why not show where in the manual and where in the paper it's answered? Some people might not even know which paper to seek! Oct 5, 2023 at 21:38
• Wouldn't it just be $\mathbf{r_i} + \frac{d \left(\mathbf{r_j}- \mathbf{r_k}\right)}{|\mathbf{r_j}- \mathbf{r_k}|}$? If you already have $\mathbf{r_i}$, then why would an angle $\theta$ that has nothing to do with the length $d$ matter? Oct 6, 2023 at 0:09
• My bad, thanks @NikeDattani So, in the book Molecular Liquids Dynamics and Interactions, there is a chapter written by H.J.C. Berendsen and W.F. van Gunsteren from Page 475-500 on Molecular dynamics simulations: techniques and approaches. Oct 6, 2023 at 6:10
• @RoshanShrestha, Yes, you are correct. I went through that book thoroughly. I did manage to derive for simple cases (as given in the book), but for this configuration, I couldn't. Oct 6, 2023 at 6:57

The formula (1) in the question:

$$\mathbf{r}_s = \underset{\text{term A}}{\mathbf{r}_i} + \underset{\text{term B}}{d \cos\theta\, \frac{\mathbf{r}_{ij}}{|\mathbf{r}_{ij}|}} + \underset{\text{term C}}{d \sin\theta\, \frac{\mathbf{r}_{\perp}}{|\mathbf{r}_{\perp}|}}\tag{1}$$

can be understood as going to the position of $$s$$ by starting at the position of $$i$$ (term A) and adding the two legs (terms B and C) of the right-angled triangle with $$\mathbf{r}_{is}$$ as its hypotenuse (equivalently, stating $$\mathbf{r}_{is}$$ in the orthonormal basis vectors of terms B and C).

Term B is the unit vector $$\hat{\mathbf{r}}_{ij}$$, and term C is the unit vector coplanar and perpendicular to $$\hat{\mathbf{r}}_{ij}$$. The (not yet unit) vector $$\mathbf{r}_\perp$$ in term C is constructed to lie in the desired plane but be perpendicular to $$\mathbf{r}_{ij}$$ by subtracting from $$\mathbf{r}_{jk}$$ its projection into $$\mathbf{r}_{ij}$$ (a Gram-Schmidt orthogonalization). This is the formula (2) in the question.

The below diagram may help visually:

• Useful relations to help understand this answer: $\sin(\theta - \pi)=-\sin(\theta),\cos(\theta-\pi)=-\cos(\theta).$ Oct 13, 2023 at 15:34

I have posted the same question in Mathematics StackExchange and here is the reply.

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.