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

[Fangyong Yan]

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/20160706/5f468f3a/attachment.htm>