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

[Fangyong Yan]

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/20160711/b38f49c4/attachment.htm>