Constrained Molecular Dynamics Method

The idea of using constrained MD to compute the PMF is simple: by using a Lagrange multiplier to fix the distance between two solutes, we are able to compute the required constraint force between two solutes. After a simple integration, modified by the volume entropy term, we obtain the PMF.

Starting with Lagrangian formulation of classical dynamics, we introduce a modified Lagrangian $L^{\prime}$

$$L^{\prime}=\sum_{i}\frac{1}{2}m_{i}\dot{\vec{R}_{i}}^{2}-U(\{\vec{R}_{i}\})-\lambda[(\vec{R}_{1}-\vec{R}_{2})^{2}-c^{2}]$$ (27)

where $\lambda$ is the Lagrange multiplier, $\vec{R}_{1}$ and $\vec{R}_{2}$ are the coordinates of the carbon atom in methanes. Let:

$$\sigma_{1}(\vec{R}_{1},\vec{R}_{2})=\vert\vec{R}_{1}-\vec{R}_{2}\vert^{2}-c^{2}=R_{12}^{2}-c^{2}.$$ (28)

From Lagrange's equation

$$\frac{d}{dt}\frac{\partial L}{\partial\dot{\vec{R}}_{i}}-\frac{\partial L}{\partial\vec{R}_{i}}=0,$$ (29)

we obtain the equation of motion

$$m_{i}\ddot{\vec{R}}_{i}=-\frac{\partial U}{\partial\vec{R}_{... ...artial\sigma_{1}}{\partial\vec{R}_{i}}=\vec{F}_{i}+\vec{G}_{i}$$ (30)

The constraint force $\vec{G}_{i}$ is given by

$$\vec{G}_{1}=-\lambda\cdot2(\vec{R}_{1}-\vec{R}_{2})$$ (31)

$$\vec{G}_{2}=+\lambda\cdot2(\vec{R}_{1}-\vec{R}_{2}).$$ (32)

By its construction, the following condition holds at any $t$:

$$\vert\vec{R}_{1}(t)-\vec{R}_{2}(t)\vert=c$$ (33)

Using Eq. (14), Verlet algorithm, we obtain

$$\vec{R}_{1}(t+\Delta t)=-\vec{R}_{1}(t-\Delta t)+2\vec{R}_{1... ...ac{\Delta t^{2}}{m_{1}}\lambda[(\vec{R}_{1}(t)-\vec{R}_{2}(t)]$$ (34)

$$\vec{R}_{2}(t+\Delta t)=-\vec{R}_{2}(t-\Delta t)+2\vec{R}_{2... ...c{\Delta t^{2}}{m_{2}}\lambda[(\vec{R}_{1}(t)-\vec{R}_{2}(t)].$$ (35)

Subtracting the last two equation, let

$$\vec{R}_{2}(t+\Delta t)-\vec{R}_{1}(t+\Delta t)=\overrightarrow{A}+\lambda\overrightarrow{B}$$ (36)

where

$$\overrightarrow{A}=2\left[\vec{R}_{2}(t)-\vec{R}_{1}(t)\righ... ...m_{2}}\vec{F}_{2}(t)-\frac{(\Delta t)^{2}}{m_{1}}\vec{F}_{1}(t)$$

$$\overrightarrow{B}=-2\frac{(\Delta t)^{2}}{m_{2}}\left[\vec{... ...elta t)^{2}}{m_{1}}\left[\vec{R}_{2}(t)-\vec{R}_{1}(t)\right].$$ (37)

Taking square on both sides, we get

$$c^{2}=A^{2}+2\lambda A\cdot B+\lambda^{2}B^{2}.$$ (38)

This tells us how to solve $\lambda$ from the known $\vec{A}(t,t-\Delta t)$ and $\vec{B}(t)$:

$$\lambda=\frac{-2\overrightarrow{A}\cdot\overrightarrow{B}\pm... ...w{A}\cdot\overrightarrow{B})^{2}-4B^{2}(A^{2}-c^{2})}}{2B^{2}}$$ (39)

There are two roots for Eq. (39), and we shall use the smaller one. The larger root corresponds to the swap of two atoms, not what we want in a dynamics. Once one solves $\lambda$, the constraint force on the first atom is given by Eq. (31).