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

Fangyong Yan yanfa... at gmail.com
Wed Jul 6 16:05:24 UTC 2016

```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>
```

More information about the CP2K-user mailing list