Cross posted on Math SE

I am trying to understand virtual sites in MD simulations, and I came across this configuration:

Virtual sites diagram

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} $$

  • $\begingroup$ I think it's already mentioned in gromacs manual and the subsequent paper. $\endgroup$ Oct 5, 2023 at 20:22
  • $\begingroup$ @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! $\endgroup$ Oct 5, 2023 at 21:38
  • $\begingroup$ 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? $\endgroup$ Oct 6, 2023 at 0:09
  • $\begingroup$ 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. $\endgroup$ Oct 6, 2023 at 6:10
  • $\begingroup$ @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. $\endgroup$
    – Vasista
    Oct 6, 2023 at 6:57

2 Answers 2


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:

Geometry of atoms k, j, i, and virtual site s, showing the vector is as the sum of basis vectors parallel and perpendicular to in.

  • $\begingroup$ Useful relations to help understand this answer: $\sin(\theta - \pi)=-\sin(\theta),\cos(\theta-\pi)=-\cos(\theta).$ $\endgroup$ 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.


You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .