flag -mol in g_msd
The tool g_msd gives inconsistent results with the flag -mol. If the diffusion coefficient is calculated for a single molecule using a corresponding index file, the result differs from an analysis for the same molecule, but using the
#1 Updated by Simon Butler almost 9 years ago
I've just posted this to gmx-users but I guess it should probably go here too. My tests suggest these patches resolve both problems with the -mol flag.
Regarding the excessive memory usage of g_msd when using the mol flag, as reported recently, I believe I've identified the source of the problem: the array of per-molecule MSD data in curr->lsq. In particular, it appears that the call to gmx_stats_add_point on line 476 is the key offender.
I haven't the time at the moment to dig any further than this, unfortunately, but the attached patch (gmx_msd.mem.patch) may be of temporary use to users of g_msd. It merely comments out the call to gmx_stats_init and the printmol routine. If anybody wishes to retain the calculation of D for each individual molecule, they will need to recode the data accumulation in the style of that for the overall calculation I expect.
There is also the problem of the large difference in the calculated results when using the -mol flag, as reported by Florian. I think the cause of this is an error in the logic within the corr_loop routine. The original sequence is:
1. If first iteration, copy current frame to previous
2. If -mol, make molecules whole
3. Remove PBC jumps
4. If -mol, calculate molecule COMs and copy to xa array
I believe the correct order for this sequence should be 2, 4, 1, 3. The original order results in xa[prev] containing zero for every position on step 1 and also in the PBC step operating on an array that hasn't been repopulated yet, with the result that jumps are not removed correctly.
I've tested the new sequence for a single PF6 molecule and for the corresponding P atom (which is almost exactly at the COM) and obtained very nearly identical results for each (which wasn't the case previously). I've attached a separate patch to reorder these function calls (gmx_msd.mol.patch). Both patches are for the v4.5.4 copy of gmx_msd.c, by the way.