a question about constrained MD, shake algorithm and the free energy calculation using Lagrange multiplier

[Fangyong Yan]

Just as these papers refer to, a constrained MD simulation is to simulate a 
rare event. 

Suppose for a reaction, A + B -> C, and it takes hours to happen in 
reality, and the reason is, A irregularly hit B, and only 1/10000 times, 
that is, 10000 times, A hit B once. 

Now, the constrained MD does this, it constrained the distance A-B, say, 
2.5 Angstrom, and instead of 10000 times hit once, since it is constrained, 
in 10000 times, it hit 10000 times, and you get 10000 times more 
opportunity to make the reaction to happen, by the way, a reaction doesnot 
mean a "real reaction", it can mean a inter-molecular complex. 

Fangyong


On Monday, July 11, 2016 at 12:22:13 PM UTC-4, Fangyong Yan wrote:
>
> Finally with the help from my advisor, I figured out the issue, 
> shake uses zeta** 2 = (ri - rj) **2, thus, zeta = |ri - rj|, which is 
> exactly in Ciccotti's paper, now the force on atom i is,
>
> F_i_total = f_i - Lagrange multiplier * the first derivative of zeta 
> versus ri,
> let the lagrange multiplier equal to lamda, and 
> the first derivative of zeta versus ri is equal to, 
> d|ri - rj| / d(r_i), 
> since it is the first derivative of the absolute value, it takes more 
> care, let y = |ri - rj| = (u**2)**1/2, and u = (ri - rj)
> then d|ri - rj| / d(ri) = (dy / du) * (du / dri) = [1/2* 2* u * u' / 
> (u**2)**1/2] * du / d(ri) = [u * u' / (u**2)**1/2] * u' = (ri-rj) / |ri - 
> rj|,
> and d|ri - rj| / d(rj) = (rj - ri) / |ri - rj|,
>
> now the force F_i_total = f_i - lambda * (ri - rj) / |ri - rj|,
>                        F_j_total = f_j - lambda * (rj - ri) / |ri - rj|,
>
> and the total force on i and j, or on the constrained bond is, 
>
> F_j_total - F_i_total = (f_j - f_i) - lambda * 2 * (rj - ri) / |ri - rj|, 
> now projection the force on ij, we get this forces projection on the 
> vector ij, we get 
> |f_j - f_i| * cos_theta - lambda * 2, theta is the angle between force 
> vector f_j-f_i and position vector r_j-r_i, and cos_theta' for lambda is 1 
> since lambda lies at the same line with r_j-r_i. And it becomes total force 
> on the bond ij = |f_j - f_i| * cos_theta - 2 * lambda. 
>
> I am not good with math at all, so please correct me if you find some 
> mistake in these eqs. Thanks!
>
> Fangyong
>
> 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
>>
>>
>>
>>
>>
>>
>>
>>
>>
>>
>>  
>>
>>
>>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://lists.cp2k.org/archives/cp2k-user/attachments/20160817/b0dfd7da/attachment.htm>