[gmx-users] Fwd: RMSF => still confused ?

Thu Oct 28 10:56:31 CEST 2010

Hi Lin,

Fair questions.

Worst case: imagine a straight path from A to B, chopped into n=10 stretches d:

A----------B

Something goes from A to B linearly and you make a hundred
observations of the position, ten at each stretch. For each stretch i,
the average position is i*d-d/2 and the fluctuation is calculated
about that position, but only includes deviations between -d/2 and
+d/2. The total fluctuation is around (A+B)/2, and includes deviations
with sizes up to +/- (B-A)/2. No surprise that the total RMSF in this
case is far larger than the per bin RMSF. In this example the RMSF
values would end up similar for each bin. If there are some obstacles
along the way, the motion won't be uniform and there will be
differences between bins, but as long as there's a trend going from A
to B (non-equilibrium situation), the total RMSF will be larger than
the RMSF per bin.

Hope this answers your questions.

Cheers,

Tsjerk

---------- Forwarded message ----------
From: Chih-Ying Lin <chihying2008 at gmail.com>
Date: Thu, Oct 28, 2010 at 8:28 AM
Subject: RMSF => still confused ?
To: Tsjerk Wassenaar <tsjerkw at gmail.com>

Hi Tsjerk:
Well.... i am still confused about the RMSF.

1. Consider a particle, it has been observed appearing at 100 positions.
    We have ten observations for one time zone. It took the duration
of 10 time zones to get 100 positions

2. The reference position is at the origin. (0,0,0)

3. For each time zones (10 positions), we calculate its standard deviation.
    => Then, I can get 10 different standard deviations for 10 time
zones, SD1, SD2, SD3.......  SD10

4. Further, consider the 10 time zone as a big-whole time zone.
    I calculate the standard deviation for the 100 positions, say SD-whole.

5. Then I found
    SD-whole > SD1,
    SD-whole > SD2
    SD-whole > SD3
    SD-whole > SD4
    SD-whole > SD5
    SD-whole > SD6
    SD-whole > SD7
    SD-whole > SD8
    SD-whole > SD9
    SD-whole > SD10

6. I suppose that
     min(SD1~SD10) <= SD-whole <= max(SD1~SD10)    ?
    As you said,
    "the mean squared deviation of the path from A to B around (A+B)/2,"

7.
Finally, I test the md-results downloading from your website.

the full set of results
http://nmr.chem.uu.nl/~tsjerk/course/md-tutorial/

THere are 4000 ps. For every 400 ps, I calculate RMSF.
Also, I calculate RMSF(from 0 to 4000 ps) as whole,

I found the same thing that
    For each alpha-carbon,
    SD-whole > SD1,
    SD-whole > SD2
    SD-whole > SD3
    SD-whole > SD4
    SD-whole > SD5
    SD-whole > SD6
    SD-whole > SD7
    SD-whole > SD8
    SD-whole > SD9
    SD-whole > SD10

WHY ???

THank you
Lin

Hi,
Lin, please think your questions over thoroughly in stead of flushing
every thought right to the mailing list. It also helps to stick to a
certain subject (reply) to make sure everything ends up in the same
thread. Maybe it's not a bad idea to read over
http://www.catb.org/esr/faqs/smart-questions.html (again).
Anyway, this question's quite valid :), and Marks answer is a bit
off... The reference structure for computing the RMSF is only used for
fitting; you'd also have seen this behaviour if you'd taken an average
structure from an equilibrium simulation. The fluctuations are around
the mean structure. The larger "fluctuations" in the first part of the
simulation are caused by the relaxation from the starting structure.
<oversimplification>If you think of the starting structure and the
'equilibrium structure' as position A and B, then the first part shows
you the mean squared deviation of the path from A to B around (A+B)/2,
whereas the next parts give you the mean squared deviation around
B</oversimplification>. Note that this gives you a means to assess
whether you've reached equilibrium, albeit not a sufficient measure.
Hope it helps,
Tsjerk
On Tue, Oct 26, 2010 at 8:24 AM, Mark Abraham <mark.abraham at anu.edu.au> wrote:
> If the reference structure from which the fluctuations are measured is taken
> from the -s file (check the documentation or code), then a non-equilibrium
> trajectory could well have this property w.r.t. a crystal structure.
>
> Mark
>
> ----- Original Message -----
> From: Chih-Ying Lin <chihying2008 at gmail.com>
> Date: Tuesday, October 26, 2010 15:55
> Subject: [gmx-users] g_rmsf => average over # of time frames ???
> To: gmx-users at gromacs.org
>
>>
>>
>> Hi
>>
>> From source code => gmx_rmsf.c
>
>> "g_rmsf computes the root mean square fluctuation (RMSF, i.e. standard ",
>
>>     "deviation) of atomic positions ",
>
>>
>
>>     if (devfn) {
>
>>       /* Calculate RMS Deviation */
>
>>       for(i=0;(i<isize);i++) {
>
>> aid = index[i];
>
>> for(d=0;(d<DIM);d++) {
>
>>   rmsd_x[i][d] += sqr(x[aid][d]-xref[aid][d]);
>
>> }
>
>>       }
>
>>     }
>
>>     count += 1.0;
>
>>
>
>> rmsf[i] = (rmsd_x[i][XX]+rmsd_x[i][YY]+rmsd_x[i][ZZ])/count;
>
>>
>
>>
>
>> Therefore, g_rmsf is the average of structure deviation over the time
>> frames.
>
>>
>
>> However, I issued the commands  ( C-alpha is selected )
>
>> g_rmsf -f abc.xtc -b 0 -e 100 -s abc-crystal.tpr -o RMSF-abc-0th-100th.xvg
>
>> g_rmsf -f abc.xtc -b 100 -e 200 -s abc-crystal.tpr -o
>> RMSF-abc-100th-200th.xvg
>
>> g_rmsf -f abc.xtc -b 200 -e 300 -s abc-crystal.tpr -o
>> RMSF-abc-200th-300th.xvg
>
>> g_rmsf -f abc.xtc -b 300 -e 400 -s abc-crystal.tpr -o
>> RMSF-abc-300th-400th.xvg
>
>> g_rmsf -f abc.xtc -b 400 -e 500 -s abc-crystal.tpr -o
>> RMSF-abc-400th-500th.xvg
>
>> g_rmsf -f abc.xtc -b 500 -e 600 -s abc-crystal.tpr -o
>> RMSF-abc-500th-600th.xvg
>
>> g_rmsf -f abc.xtc -b 600 -e 700 -s abc-crystal.tpr -o
>> RMSF-abc-600th-700th.xvg
>
>> g_rmsf -f abc.xtc -b 700 -e 800 -s abc-crystal.tpr -o
>> RMSF-abc-700th-800th.xvg
>
>> g_rmsf -f abc.xtc -b 800 -e 900 -s abc-crystal.tpr -o
>> RMSF-abc-800th-900th.xvg
>
>> g_rmsf -f abc.xtc -b 900 -e 1000 -s abc-crystal.tpr -o
>> RMSF-abc-900th-1000th.xvg
>
>>
>
>>
>
>> Also, (  C-alpha is selected )
>
>> g_rmsf -f abc.xtc -b 0 -e 1000 -s abc-crystal.tpr -o
>> RMSF-abc-0th-1000th.xvg
>
>>
>
>>
>
>>
>
>> Then I ploted,
>
>> RMSF-abc-0th-100th.xvg
>
>> RMSF-abc-100th-200th.xvg
>
>> RMSF-abc-200th-300th.xvg
>
>> RMSF-abc-300th-400th.xvg
>
>> RMSF-abc-400th-500th.xvg
>
>> RMSF-abc-500th-600th.xvg
>
>> RMSF-abc-600th-700th.xvg
>
>> RMSF-abc-700th-800th.xvg
>
>> RMSF-abc-800th-900th.xvg
>
>> RMSF-abc-900th-1000th.xvg
>
>>
>
>>
>
>> Also, I ploted
>
>> RMSF-abc-0th-1000th.xvg
>
>>
>
>>
>
>>
>
>> The PLOT RMSF-abc-0th-1000th.xvg has all RMSF-Values much higher than
>> those from RMSF-abc-0th-100th.xvg / RMSF-abc-100th-200th.xvg
>> / RMSF-abc-200th-300th.xvg / ......   / RMSF-abc-900th-1000th.xvg.......
>
>>
>
>>
>
>> It does not make sense... I supposed.....
>
>> Did I misunderstand something ?
>
>>
>
>>
>
>>
>
>> THank you
>
>> Lin
>
>>
>
>>
>
>>
>
>>
>
>>
>
>>
>
>>
>
>> --
>> gmx-users mailing list    gmx-users at gromacs.org
>> http://lists.gromacs.org/mailman/listinfo/gmx-users
>> Please search the archive at
>> http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
>> Please don't post (un)subscribe requests to the list. Use the
>> www interface or send it to gmx-users-request at gromacs.org.
>> Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
> --
> gmx-users mailing list    gmx-users at gromacs.org
> http://lists.gromacs.org/mailman/listinfo/gmx-users
> Please search the archive at
> http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
> Please don't post (un)subscribe requests to the list. Use the
> www interface or send it to gmx-users-request at gromacs.org.
> Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>

--
Tsjerk A. Wassenaar, Ph.D.
post-doctoral researcher
Molecular Dynamics Group
* Groningen Institute for Biomolecular Research and Biotechnology
* Zernike Institute for Advanced Materials
University of Groningen
The Netherlands

-- 
Tsjerk A. Wassenaar, Ph.D.

post-doctoral researcher
Molecular Dynamics Group
* Groningen Institute for Biomolecular Research and Biotechnology
* Zernike Institute for Advanced Materials
University of Groningen
The Netherlands