[gmx-users] Three different diffusivities by "gmx msd" /

[Johannes Hörmann]

Dear GROMACS community,

evaluating the diffusivity of as single dodecyl sulfate (DS) anion, I 
use a small trajectory chunk "md_00001.trr" 5 ps (500 frames) long for 
testing purposes. The .trr has been stripped of any atoms except the 42 
ones forming the DS chain (12H25SO4). What is more, I prepared an index 
file DS_selection.ndx with the content

     [ Surfactant ]
        1    2    3    4    5    6    7    8    9   10   11   12 13   
14   15
       16   17   18   19   20   21   22   23   24   25   26   27 28   
29   30
       31   32   33   34   35   36   37   38   39   40   41   42

in order to select the molecule of interest explicitly. As far as I 
understand the GMX 2016.4 manual and help, the commad

     gmx msd -f md_00001.trr -s md.tpr -o atom_wise_msd.xvg -n 
DS_selection.ndx

should yield the atom-wise evaluated MSD averaged over all 42 atoms for 
each timestep. "atom_wise_msd.xvg" contains the comment

     ....
     # MSD gathered over 4.99 ps with 1 restarts
     # Diffusion constants fitted from time 0.5 to 4.49 ps
     # D[Surfactant] = 1.8305 (+/- 1.3704) (1e-5 cm^2/s)
     ....

which I understand  as the atom-averaged diffusivity. On the other hand,

     gmx msd -f md_00001.trr -s md.tpr -o mol_wise_msd.xvg -mol 
mol_wise_D.xvg -n DS_selection.ndx

should yield MSD (and diffusivity D) for the molecule's COM, and 
consequently a different value D

     ....
     # MSD gathered over 4.99 ps with 1 restarts
     # Diffusion constants fitted from time 0.5 to 4.49 ps
     # D[Surfactant] = 1.5575 (+/- 1.3521) (1e-5 cm^2/s)
     ....

is found within "mol_wise_msd.xvg". So far, so good, my understanding is 
correct up until here?

Now, looking at the latter command's stdout, I have

     <D> = 1.1395 Std. Dev. = -nan Error = -nan
     Fitting from 0.5 to 4.49 ps

     D[Surfactant] 1.5575 (+/- 1.3521) 1e-5 cm^2/s

and in the file "mol_wise_D.xvg" of course just one entry for the single 
DS molecule

     ....
     @    title "Diffusion Coefficients / Molecule"
     @    xaxis  label "Molecule"
     @    yaxis  label "D"
     @TYPE xy
          0     1.13953

Here, I am confused. How does the third value for D arise? Naturally, I 
would have expected "mol_wise_D.xvg" to list the same value of D = 
1.5575 as found within the comment of "mol_wise_msd.xvg".

I am looking forward to a clarifying answer, as I could not find any 
notes on this discrepancy in the documentation. Is it maybe possible 
that the <D> = 1.1395 is evaluated over all 5ps, whereas D[Surfactant] = 
1.5575is fitted only on the 80% trajectory from 0.5 to 4.49 ps? If so, 
some clear hint towards this fact would be appreciated in the "gmx help 
msd" text.

A different topic: I could not find any option for "gmx msd" to evaluate 
a trajectory frame-wise without loading the whole file into memory. When 
applied onto a several GB .trr, even with relatively small "-trestart”， 
GMX would continuously eat up memory until the job would be canceled due 
to running out of the node's 128 GB RAM. Maybe any undocumented option 
hidden somewhere?

Here the output of "gmx --version"

     GROMACS:      gmx, version 2016.4
     Executable: 
/work/ws/nemo/fr_lp1029-IMTEK_SIMULATION-0/gromacs_2016.4_gnu/bin/gmx
     Data prefix: 
/work/ws/nemo/fr_lp1029-IMTEK_SIMULATION-0/gromacs_2016.4_gnu
     Working dir: 
/work/ws/nemo/fr_jh1130-201708-0/jobs/gmxlab/sds/201801/farafonov2017developing/1ds/rtp/GAFF-SPC
     Command line:
       gmx --version

     GROMACS version:    2016.4
     Precision:          single
     Memory model:       64 bit
     MPI library:        thread_mpi
     OpenMP support:     enabled (GMX_OPENMP_MAX_THREADS = 32)
     GPU support:        disabled
     SIMD instructions:  AVX2_256
     FFT library:        fftw-3.3.7-sse2-avx-avx2-avx2_128
     RDTSCP usage:       enabled
     TNG support:        enabled
     Hwloc support:      hwloc-1.11.0
     Tracing support:    disabled
     Built on:           Wed Dec 13 10:46:41 CET 2017
     Built by:           fr_jh1130 at int02.nemo.privat [CMAKE]
     Build OS/arch:      Linux 3.10.0-693.5.2.el7.x86_64 x86_64
     Build CPU vendor:   Intel
     Build CPU brand:    Intel(R) Xeon(R) CPU E5-2630 v4 @ 2.20GHz
     Build CPU family:   6   Model: 79   Stepping: 1
     Build CPU features: aes apic avx avx2 clfsh cmov cx8 cx16 f16c fma 
hle htt lahf mmx msr nonstop_tsc pcid    pclmuldq pdcm pdpe1gb popcnt 
pse rdrnd rdtscp rtm sse2 sse3 sse4.1 sse4.2 ssse3 tdt x2apic
     C compiler: /opt/bwhpc/common/compiler/gnu/5.2.0/bin/gcc GNU 5.2.0
     C compiler flags:    -march=core-avx2     -O3 -DNDEBUG 
-funroll-all-loops -fexcess-precision=fast
     C++ compiler: /opt/bwhpc/common/compiler/gnu/5.2.0/bin/g++ GNU 5.2.0
     C++ compiler flags:  -march=core-avx2    -std=c++0x   -O3 -DNDEBUG 
-funroll-all-loops -fexcess-precision=fast

Thanks for your help.

Best regards,

-- 
Johannes Hörmann

University of Freiburg
IMTEK - Department of Microsystems Engineering
Laboratory for Simulation
Georges-Köhler-Allee 103, Room 03 020
79110 Freiburg, Germany

Office: +49-761-203-67726
Email: johannes.hoermann at imtek.uni-freiburg.de