[CP2K:8746] Re: a question about constrained MD, shake algorithm and the free energy calculation using Lagrange multiplier

Fangyong Yan yanfa... at gmail.com
Fri Mar 3 17:22:47 UTC 2017

```however, at this moment, I dont think I can find a better way to simulate
high free energy barrier for chemical reactions. Ab-initio constrained MD
seems to be the only choice, with enough sampling and good guess for the
reaction coordinate, the trajectory made up ab-initio constrained MD will
be reasonable.

On Sun, Feb 26, 2017 at 7:48 PM, Fangyong Yan <yanfa... at gmail.com> wrote:

> personally thinking, the constrained md is kind of biasing to the reaction
> coordinates, so I personally prefer running a 1000 ps unconstrained MD, to
> running the same length of constrained md. For some reactions with high
> free energy barrier, unconstrained MD cannot simulate the reaction. In this
> case, maybe I can try other methods.
>
> On Wednesday, July 6, 2016 at 12:05:25 PM UTC-4, Fangyong Yan wrote:
>
>> Hi,
>>
>> 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!
>>
>> Fangyong
>>
>>
>>
>>
>>
>>
>>
>>
>>
>>
>>
>>
>>
>> --
> You received this message because you are subscribed to a topic in the
> To unsubscribe from this topic, visit https://groups.google.com/d/
> topic/cp2k/yWqahb93_38/unsubscribe.
> To unsubscribe from this group and all its topics, send an email to