[CP2K-user] [CP2K:8746] Re: a question about constrained MD, shake algorithm and the free energy calculation using Lagrange multiplier

[Annesha Debroy]

Hi,
Did you figure out, how to obtain the free energy profile? 
I obtained the average force along the reaction coordinate..but I was 
wondering if there is any script to perform the thermodynamic integration?
On Sunday, May 21, 2017 at 7:11:46 AM UTC+5:30 ya... at gmail.com wrote:

> 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.
>
>
> On Sat, May 20, 2017 at 9:24 PM, Fangyong Yan <ya... at gmail.com> wrote:
>
>> 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. 
>>
>> 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. 
>>
>> On Fri, Mar 3, 2017 at 12:22 PM, Fangyong Yan <ya... at gmail.com> wrote:
>>
>>> 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. 
>>>
>>> On Sun, Feb 26, 2017 at 7:48 PM, Fangyong Yan <ya... at gmail.com> 
>>> wrote:
>>>
>>>> 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.  
>>>>
>>>> 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
>>>>>
>>>>>
>>>>>
>>>>>
>>>>>
>>>>>
>>>>>
>>>>>
>>>>>
>>>>>
>>>>>  
>>>>>
>>>>>
>>>>> -- 
>>>> You received this message because you are subscribed to a topic in the 
>>>> Google Groups "cp2k" group.
>>>> To unsubscribe from this topic, visit 
>>>> https://groups.google.com/d/topic/cp2k/yWqahb93_38/unsubscribe.
>>>> To unsubscribe from this group and all its topics, send an email to 
>>>> cp... at googlegroups.com.
>>>> To post to this group, send email to c... at googlegroups.com.
>>>> Visit this group at https://groups.google.com/group/cp2k.
>>>> For more options, visit https://groups.google.com/d/optout.
>>>>
>>>
>>>
>>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://lists.cp2k.org/archives/cp2k-user/attachments/20210520/b28c55f8/attachment.htm>