<div>Hi,</div><div>Did you figure out, how to obtain the free energy profile? <br></div><div>I obtained the average force along the reaction coordinate..but I was wondering if there is any script to perform the thermodynamic integration?<br></div><div class="gmail_quote"><div dir="auto" class="gmail_attr">On Sunday, May 21, 2017 at 7:11:46 AM UTC+5:30 ya...@gmail.com wrote:<br/></div><blockquote class="gmail_quote" style="margin: 0 0 0 0.8ex; border-left: 1px solid rgb(204, 204, 204); padding-left: 1ex;"><div dir="ltr"><span style="font-size:12.8px">I am not talking about using classical MD to do reactions, I am just saying that in classical MD, there are also free energy mapping, which includes local minima, maxima, and saddle points.</span><div><br></div></div><div class="gmail_extra"><br><div class="gmail_quote">On Sat, May 20, 2017 at 9:24 PM, Fangyong Yan <span dir="ltr"><<a href data-email-masked rel="nofollow">ya...@gmail.com</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div dir="ltr">I think in order to use constrained MD, we need to locate the reaction coordinate, which is not an easy task. I think lots of people have been working on such problems, how to map the free energy surface, not only in ab-initio MD calculations, but also in classical MD calculations, how to locate the reactants, saddle points (transition states), and products. In quantum mechanics calculation, we use potential energy surface, in MD, we use free energy surface. <div><br></div><div>I am still trying to learn all these, but I think in order to understand our problems, we need to have a better understanding of our energy surfaces, which is a very complicated problem. </div></div><div><div><div class="gmail_extra"><br><div class="gmail_quote">On Fri, Mar 3, 2017 at 12:22 PM, Fangyong Yan <span dir="ltr"><<a href data-email-masked rel="nofollow">ya...@gmail.com</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div dir="ltr">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. </div><div><div><div class="gmail_extra"><br><div class="gmail_quote">On Sun, Feb 26, 2017 at 7:48 PM, Fangyong Yan <span dir="ltr"><<a href data-email-masked rel="nofollow">ya...@gmail.com</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div dir="ltr">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. <div><span><br>On Wednesday, July 6, 2016 at 12:05:25 PM UTC-4, Fangyong Yan wrote:</span><div><div><blockquote class="gmail_quote" style="margin:0;margin-left:0.8ex;border-left:1px #ccc solid;padding-left:1ex"><div dir="ltr">Hi,<div><br></div><div>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, </div><div><br></div><div>sigma = (ri - rj) ** 2 - dij ** 2, where dij is the constrain distance,</div><div>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, </div><div>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)).</div><div><br></div><div>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, </div><div>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 </div><div><br></div><div>dW / d Zeta' = < Z^(-1/2) * [ -lamda + kTG] > / < Z^(-1/2)></div><div><br></div><div>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 </div><div>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, </div><div>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 </div><div><br></div><div>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, <br></div><div>dW / d Zeta' = <-lamda><br></div><div><br></div><div>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. </div><div><br></div><div>Thanks for your patience for reading this, and I hope someone who can help me with this issue!</div><div><br></div><div>Fangyong</div><div><br></div><div><br></div><div><br></div><div><br></div><div><br></div><div><br></div><div><br></div><div><br></div><div><br></div><div><br></div><div> </div><div><br></div><div><br></div></div></blockquote></div></div></div></div><div><div> <p></p> -- <br> You received this message because you are subscribed to a topic in the Google Groups "cp2k" group.<br> To unsubscribe from this topic, visit <a href="https://groups.google.com/d/topic/cp2k/yWqahb93_38/unsubscribe" target="_blank" rel="nofollow" data-saferedirecturl="https://www.google.com/url?hl=en&q=https://groups.google.com/d/topic/cp2k/yWqahb93_38/unsubscribe&source=gmail&ust=1621594889469000&usg=AFQjCNG088WEasbi_aBxHAEAMhT9QFqzTw">https://groups.google.com/d/topic/cp2k/yWqahb93_38/unsubscribe</a>.<br> To unsubscribe from this group and all its topics, send an email to <a href data-email-masked rel="nofollow">cp...@googlegroups.com</a>.<br> To post to this group, send email to <a href data-email-masked rel="nofollow">c...@googlegroups.com</a>.<br> Visit this group at <a href="https://groups.google.com/group/cp2k" target="_blank" rel="nofollow" data-saferedirecturl="https://www.google.com/url?hl=en&q=https://groups.google.com/group/cp2k&source=gmail&ust=1621594889470000&usg=AFQjCNGtk66KWC5OgLnNfGT4FTs9lWzUGQ">https://groups.google.com/group/cp2k</a>.<br> For more options, visit <a href="https://groups.google.com/d/optout" target="_blank" rel="nofollow" data-saferedirecturl="https://www.google.com/url?hl=en&q=https://groups.google.com/d/optout&source=gmail&ust=1621594889470000&usg=AFQjCNEkKPBBWYTFOnxC4Rph2TuHbJ1TQg">https://groups.google.com/d/optout</a>.<br> </div></div></blockquote></div><br></div> </div></div></blockquote></div><br></div> </div></div></blockquote></div><br></div> </blockquote></div>