a question about constrained MD, shake algorithm and the free energy calculation using Lagrange multiplier
yanfa... at gmail.com
Sat Jan 7 22:28:16 CET 2017
constrained MD simulations belong to potential of mean force (free energy)
calculation, there are many very good reviews for free energy calculation,
"FREE ENERGY VIA MOLECULAR SIMULATION: Applications to Chemical and
Biomolecular Systems D. L. Beveridge and F. M. DiCapua", Annu. Rev.
Biophys. Biophys. Chem. 1989. 18:431-92, and many others.
And Dann Frenkel and Berend Smit's book is also helpful for understand
details of MD/MC simulations.
I have lots of fun in reading and understanding, even though I am a very
On Wednesday, July 6, 2016 at 12:05:25 PM UTC-4, Fangyong Yan wrote:
> I have a question about the free energy calculation using the constrained
> MD. For the simplest case, such as constraining a inter-molecular distance
> between two atoms, i and j, In the constrain MD in the NVT ensemble, CP2K
> uses shake algorithm to update the position and velocity, where the
> constrain follows the holonomic constrain,
> sigma = (ri - rj) ** 2 - dij ** 2, where dij is the constrain distance,
> and the total force is equal to F_i + G_i, G_i is the constrained force
> and is equal to, lamda * the first derivative of simga versus r_i, thus,
> G_i = -2 * lamda * r_i, (where these eq. borrows from the original shake
> paper, JEAN-PAUL RYCKAERT, GIOVANNI CICCOTTI, AND HERMAN J. C. BERENDSEN,
> JOURNAL OF COMPUTATIONAL. PHYSICS 23, 321-341 (1977)).
> In the free energy calculation, I think CP2K uses the eq. derived by
> Michiel Sprik and GIOVANNI CICCOTTI, Free energy from constrained molecular
> dynamics, J. Chem. Phys., Vol. 109, No. 18, 8 November 1998, where in this
> paper, the free energy uses a different constrain,
> where constrain is equal to |ri - rj| - dij = 0, "| |" represents the
> absolute value, and in this case, the constrained force G_i = - lamda, (see
> eq. 13 in the paper). The free energy is equal to
> dW / d Zeta' = < Z^(-1/2) * [ -lamda + kTG] > / < Z^(-1/2)>
> W is the free energy, Zeta is the constrained eq., in this case is equal
> to |ri - rj| - dij = 0, Zeta' represent different Zeta's; < > is the
> ensemble average, Z is a factor arises from the requirement that when Zeta
> is equal to zero for all times, the first derivative of Zeta (the velocity
> of this constrain) is also equal to zero for all times. (from E.A. CARTER,
> Giovanni CICCOTTI, James T. HYNES, Raymond KAPRAL, Chem. Phys. Lett. 156,
> 472 ~1989.); G is equal to
> G = (1 / Z^2) * (1/m_i * 1/m_j) * the first derivative of Zeta versus r_i
> * the second derivative of Zeta versus r_i and r_j * the first derivative
> of Zeta versus r_j,
> when Zeta = |r_i - r_j| - dij, the first derivative of Zeta versus r_i =
> the first derivative of Zeta versus r_j = 1, the second derivative of Zeta
> versus r_i and r_j = 0, thus, the free energy is equal to
> dW / d Zeta' = < Z^(-1/2) * [ -lamda + kTG] > / < Z^(-1/2)> = < Z^(-1/2) *
> [ -lamda] > / < Z^(-1/2)>, and Z is a constant in this simple case, thus,
> dW / d Zeta' = <-lamda>
> Now my question is, since shake uses Zeta = (r_i - r_J) ** 2 - dij**2 = 0,
> in this case, G wont disappear, and the constrained force G_i = - 2 * lamda
> * r_i. Since CP2K does use SHAKE algorithm, how does CP2K do the free
> energy calculation, do CP2K uses Zeta = (r_i - r_j) ** 2 - dij**2 =0, or
> Zeta = |r_i - r_j| - dij = 0, since these two cases the lagrange multiplier
> is different.
> Thanks for your patience for reading this, and I hope someone who can help
> me with this issue!
-------------- next part --------------
An HTML attachment was scrubbed...
More information about the CP2K-user