Re: OutputAccumulatedWork...where's the output?

From: Giacomo Fiorin (giacomo.fiorin_at_gmail.com)
Date: Mon Dec 14 2015 - 13:32:17 CST

Ok, that is a separate issue of equilibrium work vs. non-equilibrium work.
It is also the same issue you can find with both implementations of TMD, of
course.

I would examine the atomic trajectory (i.e. the DCD file), to see whether
the native structure is recovered (you start from a 8 Å RMSD). If you
ended up kinetically trapped in a non-native state, the energy can do
nothing but increase.

Note also that even if you do reach the RMSD = 0 state, there will always
be significant non-equilibrium work due to finite time (e.g. the 10 ns you
used).

Lastly, RMSD = 0 has only one microstate, and its entropy is extremely low,
i.e. you have a divergently high free energy. Rather than zero, I would
set as target value the RMSD that is observed in a MD simulation that is
initialized from the native structure.

Giacomo

On Mon, Dec 14, 2015 at 1:41 PM, Bryan Roessler <roessler_at_uab.edu> wrote:

> Giacomo,
>
> It is the same value, so it appears that the W-<variable name> is indeed
> the accumulating work.
>
> The reason I am confused is that the accumulated work seems to increase
> harmonically, whereas in other published experiments the accumulating work
> is dependent on the reaction coordinate (transition tube). As the protein
> progresses down its favorable transition state, the accumulating work
> should slow, and as it enters unfavorable metastable states, it should grow
> faster (to pull the protein out of those states towards the targeted
> structure). I'm seeing a perfectly harmonic curve. Could this be caused by
> a force constant that is too high?
>
> I'm working on a fairly large membrane transporter and am using a K of 200
> kcal/mol as is recommended by the NAMD manual as a good starting point for
> C-alpha TMD (except that I am using the colvar module and not the built-in
> method).
>
> Thanks,
>
> Bryan
>
> *Bryan Roessler | Graduate Research Assistant*
> UAB | The University of Alabama at Birmingham
> *uab.edu/cmdb <http://uab.edu/cmdb>*
> Knowledge that will change your world
>
> On Mon, Dec 14, 2015 at 10:55 AM, Giacomo Fiorin <giacomo.fiorin_at_gmail.com
> > wrote:
>
>> Hello Bryan, the output is in the column labeled as "W_<variable name>".
>>
>> Does the last number in that column correspond to what you have in the
>> state file?
>>
>> Giacomo
>>
>> On Sun, Dec 13, 2015 at 8:36 PM, Bryan Roessler <roessler_at_uab.edu> wrote:
>>
>>> Hello,
>>>
>>> I am performing an RMSD TMD using the colvar module. I have enabled
>>> outputaccumulatedwork, defined colvarsTrajFrequency in my colvars.conf
>>> file, and set a targetCenters as 0.
>>>
>>> #colvars.conf
>>> colvarsTrajFrequency 100
>>>
>>> colvar {
>>> name RMSD
>>> outputAppliedForce on # keep track of bias force on this variable
>>> rmsd {
>>> atoms {
>>> atomsfile target.pdb # Select biased atoms from this file
>>> atomsCol O # based on occupancy (all nonzero values)
>>> }
>>> refPositionsFile target.pdb # ref. positions are selected based on
>>> atom number
>>> }
>>> }
>>>
>>> harmonic { # Define a moving harmonic restraint
>>> colvars RMSD
>>> centers 8.108955 # go from x Angstrom RMSD...
>>> targetCenters 0.0 # ... to 0 Angstrom
>>> targetNumSteps 5000000 # in 5 million MD steps
>>> forceConstant 200. # in kcal/mol/A^2
>>> outputAccumulatedWork yes # keep track of how much energy will be
>>> added
>>> }
>>>
>>> In my colvars.traj file, I get the following result:
>>>
>>> # step RMSD fa_RMSD
>>> W_harmonic1
>>> 0 8.10893758151549e+00 3.48369690286177e-03
>>> -5.64982828378909e-09
>>> 100 8.11052619012773e+00 -3.46673845550072e-01
>>> 5.65938731303684e-05
>>> 200 8.11502663013762e+00 -1.27919766752989e+00
>>> 1.31113744891078e-04
>>> 300 8.11885738603493e+00 -2.07778466699651e+00
>>> 3.83688077215023e-04
>>> 400 8.12532180221536e+00 -3.40310372308537e+00
>>> 8.18258935522487e-04
>>> 500 8.12519728600043e+00 -3.41063630010119e+00
>>> 1.44251410460275e-03
>>> 600 8.10888978266427e+00 -1.81571452874252e-01
>>> 1.78822263431076e-03
>>> 700 8.10725589897312e+00 1.12769465353679e-01
>>> 1.71416299364519e-03
>>> 800 8.11735001417091e+00 -1.93848939420818e+00
>>> 2.02673879603268e-03
>>> 900 8.11555019832383e+00 -1.61096204479421e+00
>>> 2.34953782273829e-03
>>> 1000 8.11527862312204e+00 -1.58908282444123e+00
>>> 2.65996804746196e-03
>>> 1100 8.11598268331998e+00 -1.76233068403207e+00
>>> 2.98824694512575e-03
>>> 1200 8.11313380101662e+00 -1.22499004336234e+00
>>> 3.29052949104976e-03
>>> 1300 8.10521083280471e+00 3.27167779015269e-01
>>> 3.28998207179514e-03
>>> 1400 8.12048965688006e+00 -2.76103285605771e+00
>>> 3.37644752668391e-03
>>> 1500 8.11309309302196e+00 -1.31415590443957e+00
>>> 3.70217004390935e-03
>>> 1600 8.11441076220149e+00 -1.61012556035054e+00
>>> 3.90344676028449e-03
>>> 1700 8.11531243386224e+00 -1.82289571250216e+00
>>> 4.13469721710158e-03
>>> 1800 8.10407381960276e+00 3.92391319390129e-01
>>> 4.30633172564267e-03
>>> 1900 8.10543201783556e+00 8.83158528271366e-02
>>> 4.25614003062294e-03
>>> etc.
>>>
>>> and in traj.state:
>>>
>>> configuration {
>>> step 5000000
>>> dt 2.000000e+00
>>> }
>>>
>>> colvar {
>>> name RMSD
>>> x 2.25972252769051e+00
>>> }
>>>
>>> restraint {
>>> configuration {
>>> name harmonic1
>>> centers -5.85060571139554e-10
>>> centers_raw -5.85060571139554e-10
>>> accumulatedWork 1.21825523336397e+03
>>> }
>>> }
>>>
>>>
>>> So obviously the accumulated work is being calculated, but only output
>>> when the .state file is written (final step), unless I'm missing it
>>> elsewhere.
>>>
>>> Thanks,
>>> Bryan
>>>
>>> *Bryan Roessler | Graduate Research Assistant*
>>> UAB | The University of Alabama at Birmingham
>>> *uab.edu/cmdb <http://uab.edu/cmdb>*
>>> Knowledge that will change your world
>>>
>>>
>>
>>
>> --
>> 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
>> http://goo.gl/Q3TBQU
>> http://giacomofiorin.github.io/
>> 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
http://goo.gl/Q3TBQU
http://giacomofiorin.github.io/
https://icms.cst.temple.edu/members.html

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