Fwd: Regarding system forces in NAMD

From: sruthi c k (cksruthikvr_at_gmail.com)
Date: Fri Jan 08 2016 - 03:38:26 CST

Hi Jerome and Giacomo,

Thanks for your clarifications.

I think my question is not clear enough, so let me clarify it.

1) I did 200ns long UNBIASED NPT simulation.
2) We do not do ABF. We use COLVAR module and system forces to postprocess
the forces for another purpose. But before doing so, we ran into confusions
discussed below.
3) Since we do not use ABF, HideJacobian option was NOT activated and from
the manual I understand that its default value is "no"
4) Hence I assume the system forces output includes Jacobian forces also.
5) I then calculated the histogram. The range of Rg was divided into 30
small bins and the number of Rg appeared in each bin was counted to get the
histogram.
6) The forces corresponding to Rg values falling in each bin was averaged
to get the <F_system(Rg)>
7) Regarding convergence of the simulation: The average forces look
*converged* when I calculate it using data obtained from the 150ns
simulation and compare it with average forces obtained in 200ns simulation.

*After following this procedure, here are two questions I have:A1)*. PMF1
defined as -kbT log(histogram) and that the PMF2 defined by -(integral of
<F_system(Rg)>) are NOT MATCHING.
I am attaching the plot of histogram corresponding to PMF1 and the *<F_system>
*corresponding to PMF2 with this mail.

I have another question related to system force which I had NOT mentioned
in the previous mail.

*B1)* I tried to calculate the system forces by taking the numerical
derivative of colvar velocities appearing in the traj file.
*B2)* Later I found that I am missing a term "m_ksi" as defined in Eqn (7)
in the paper "Adaptive biasing force method for scalar and vector free
energy calculations" by Eric Darve et al.

m_ksi is defined as
                    m_ksi = SUM((1/m_k)* (derivative of ksi wrt cartesian
coordinates)^2)

*B3)* When including m_ksi in my calculations, I was not sure how to handle
the derivatives relative to all possible cartesian coordinates, as some of
them are related by the covalent bond lengths in the system (rigid bonds
option "ALL" was used in the namd input file)

*B4)* I have read in the paper titled "The Adaptive Biasing Force Method:
Everything You Always Wanted To Know but Were Afraid To Ask"
 by Jeffrey Comer et al. that "If only the heavy atom of the bonded pair is
involved in the transition coordinate, then the adaptive biasing force
algorithm will include the constraint force emanating from the hydrogen in
addition to the thermodynamic force, thus contaminating its estimate. A
straightforward solution is to include both the hydrogen and its parent
atom in the collective variable(s) defining the transition coordinate,
causing the constraint forces to cancel each other"
*B5)* Since I calculate Radius of gyration even including hydrogen atoms
can I be not worried about the constraint forces? This question also came
up because the PMF in Question *A1* above was not coming out as expected

Thank you very much,
Sruthi

On Thu, Jan 7, 2016 at 11:13 PM, sruthi c k <cksruthikvr_at_gmail.com> wrote:

> Hi Giacomo,
>
> Thanks for your reply.
>
> As the hide Jacobian option is NOT activated, I think the system forces
> output includes the Jacobian forces also.If I understand correctly hence
> the free energy surface obtained as -kbT log(histogram) and that by
> integrating <F_system> should match.
>
> Regarding convergence: The average forces look converged when I calculate
> it using data obtained from the 150ns simulation and compare it with
> average forces obtained in 200ns simulation.
>
> Thanks
> Sruthi
>
>
> On Thu, Jan 7, 2016 at 8:51 PM, sruthi c k <cksruthikvr_at_gmail.com> wrote:
>
>> Dear Sir,
>> I asked about the system forces in the namd mailing list and I am
>> forwarding that mail to you.
>> Thanks
>> Sruthi
>>
>> ---------- Forwarded message ----------
>> From: Giacomo Fiorin <giacomo.fiorin_at_gmail.com>
>> Date: Thu, Jan 7, 2016 at 8:44 PM
>> Subject: Re: namd-l: Regarding system forces in NAMD
>> To: NAMD list <namd-l_at_ks.uiuc.edu>, Jérôme Hénin <jerome.henin_at_ibpc.fr>
>> Cc: sruthi c k <cksruthikvr_at_gmail.com>
>>
>>
>> Sorry, one clarification: I meant to say that the Jacobian force is not
>> calculated *explicitly* within the histogram, but it is implicitly
>> included in it and the histogram carries full thermodynamic information (at
>> least for that variable).
>>
>> Giacomo
>>
>>
>> On Thu, Jan 7, 2016 at 10:10 AM, Giacomo Fiorin <giacomo.fiorin_at_gmail.com
>> > wrote:
>>
>>> Hi, a small addition about the comparison with the histogram. The
>>> Jacobian force (which is very strong near zero for many distance-related
>>> variables as Jérôme said) is only calculated explicitly for system forces,
>>> but is not implemented for histograms.
>>>
>>> You can either:
>>> 1) use the hideJacobian option of ABF (and you must use ABF with
>>> "applyBias off" to calculate the system-force profile); or
>>> 2) add instead the Jacobian force to the derivative of the histogram,
>>> using the formula for that specific colvar, as listed in the tables of our
>>> reference paper.
>>>
>>> Lastly, note that until full convergence is achieved, it is perfectly
>>> normal to have different PMFs from the system force vs the histogram.
>>>
>>> Giacomo
>>>
>>> On Thu, Jan 7, 2016 at 9:37 AM, Jérôme Hénin <jerome.henin_at_ibpc.fr>
>>> wrote:
>>>
>>>> Hi Sruthi,
>>>>
>>>> If you are looking at a radius of gyration, keep in mind that that
>>>> coordinate experiences a strong Jacobian force especially at small values.
>>>> Essentially, entropy tends to pull the particles apart, while cohesive
>>>> forces from the force field keep it together - those would be almost always
>>>> negative. Did you activate the hideJacobian option?
>>>>
>>>> Jerome
>>>>
>>>> On 7 January 2016 at 13:49, sruthi c k <cksruthikvr_at_gmail.com> wrote:
>>>>
>>>>> Dear NAMD developers,
>>>>>
>>>>> I have some queries about the system forces that are calculated in the
>>>>> COLVAR module implemented in NAMD.
>>>>>
>>>>> I was trying to generate the free-energy surface by integrating the
>>>>> <F_system> on the colvar obtained from an unconstrained and unbiased MD
>>>>> simulation.
>>>>>
>>>>> I did 200 ns equilibrium molecular dynamics (NPT) simulation (NO BIAS
>>>>> APPLIED) and collected the system force acting on the colvar. The
>>>>> <F_system> for a given bin of the colvar was then calculated by taking the
>>>>> average of the system forces corresponding to the rg values falling in that
>>>>> bin.
>>>>>
>>>>> What I observe is that <F_system> is a straight-line but it is all
>>>>> negative numbers. I was expecting that <F_system> will be equal to zero
>>>>> where the histogram of colvar has a peak and +ve and -ve on either side of
>>>>> it.
>>>>>
>>>>> Could you please help me in understanding this?
>>>>> Thanks in advance
>>>>>
>>>>> Regards
>>>>> Sruthi
>>>>>
>>>>>
>>>>>
>>>>>
>>>>>
>>>>>
>>>>>
>>>>>
>>>>>
>>>>>
>>>>>
>>>>>
>>>>
>>>
>>>
>>> --
>>> Giacomo Fiorin
>>> Assistant Professor of Research
>>> Institute for Computational Molecular Science (ICMS)
>>> College of Science and Technology, Temple University
>>> 1925 North 12th Street (035-07), Room 704D
>>> Philadelphia, PA 19122-1801
>>> Phone: +1-215-204-4213
>>>
>>> Scholar: http://goo.gl/Q3TBQU
>>> Personal: http://giacomofiorin.github.io/
>>> Lab page: https://icms.cst.temple.edu/members.html
>>>
>>>
>>>
>>
>>
>> --
>> Giacomo Fiorin
>> Assistant Professor of Research
>> Institute for Computational Molecular Science (ICMS)
>> College of Science and Technology, Temple University
>> 1925 North 12th Street (035-07), Room 704D
>> Philadelphia, PA 19122-1801
>> Phone: +1-215-204-4213
>>
>> Scholar: http://goo.gl/Q3TBQU
>> Personal: http://giacomofiorin.github.io/
>> Lab page: https://icms.cst.temple.edu/members.html
>>
>>
>>
>>
>


This archive was generated by hypermail 2.1.6 : Tue Dec 27 2016 - 23:21:43 CST