2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5 * Copyright (c) 2001-2004, The GROMACS development team.
6 * Copyright (c) 2013,2014,2015,2016,2017 by the GROMACS development team.
7 * Copyright (c) 2018,2019,2020, by the GROMACS development team, led by
8 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
9 * and including many others, as listed in the AUTHORS file in the
10 * top-level source directory and at http://www.gromacs.org.
12 * GROMACS is free software; you can redistribute it and/or
13 * modify it under the terms of the GNU Lesser General Public License
14 * as published by the Free Software Foundation; either version 2.1
15 * of the License, or (at your option) any later version.
17 * GROMACS is distributed in the hope that it will be useful,
18 * but WITHOUT ANY WARRANTY; without even the implied warranty of
19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 * Lesser General Public License for more details.
22 * You should have received a copy of the GNU Lesser General Public
23 * License along with GROMACS; if not, see
24 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
25 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
27 * If you want to redistribute modifications to GROMACS, please
28 * consider that scientific software is very special. Version
29 * control is crucial - bugs must be traceable. We will be happy to
30 * consider code for inclusion in the official distribution, but
31 * derived work must not be called official GROMACS. Details are found
32 * in the README & COPYING files - if they are missing, get the
33 * official version at http://www.gromacs.org.
35 * To help us fund GROMACS development, we humbly ask that you cite
36 * the research papers on the package. Check out http://www.gromacs.org.
44 #include "gromacs/commandline/pargs.h"
45 #include "gromacs/commandline/viewit.h"
46 #include "gromacs/fileio/confio.h"
47 #include "gromacs/fileio/matio.h"
48 #include "gromacs/fileio/trxio.h"
49 #include "gromacs/fileio/xvgr.h"
50 #include "gromacs/gmxana/gmx_ana.h"
51 #include "gromacs/gmxana/gstat.h"
52 #include "gromacs/math/functions.h"
53 #include "gromacs/math/invertmatrix.h"
54 #include "gromacs/math/utilities.h"
55 #include "gromacs/math/vec.h"
56 #include "gromacs/pbcutil/pbc.h"
57 #include "gromacs/topology/index.h"
58 #include "gromacs/topology/topology.h"
59 #include "gromacs/utility/arraysize.h"
60 #include "gromacs/utility/cstringutil.h"
61 #include "gromacs/utility/futil.h"
62 #include "gromacs/utility/gmxassert.h"
63 #include "gromacs/utility/smalloc.h"
66 int gmx_vanhove(int argc
, char* argv
[])
68 const char* desc
[] = {
69 "[THISMODULE] computes the Van Hove correlation function.",
70 "The Van Hove G(r,t) is the probability that a particle that is at r[SUB]0[sub]",
71 "at time zero can be found at position r[SUB]0[sub]+r at time t.",
72 "[THISMODULE] determines G not for a vector r, but for the length of r.",
73 "Thus it gives the probability that a particle moves a distance of r",
75 "Jumps across the periodic boundaries are removed.",
76 "Corrections are made for scaling due to isotropic",
77 "or anisotropic pressure coupling.",
79 "With option [TT]-om[tt] the whole matrix can be written as a function",
80 "of t and r or as a function of [SQRT]t[sqrt] and r (option [TT]-sqrt[tt]).",
82 "With option [TT]-or[tt] the Van Hove function is plotted for one",
83 "or more values of t. Option [TT]-nr[tt] sets the number of times,",
84 "option [TT]-fr[tt] the number spacing between the times.",
85 "The binwidth is set with option [TT]-rbin[tt]. The number of bins",
86 "is determined automatically.",
88 "With option [TT]-ot[tt] the integral up to a certain distance",
89 "(option [TT]-rt[tt]) is plotted as a function of time.",
91 "For all frames that are read the coordinates of the selected particles",
92 "are stored in memory. Therefore the program may use a lot of memory.",
93 "For options [TT]-om[tt] and [TT]-ot[tt] the program may be slow.",
94 "This is because the calculation scales as the number of frames times",
95 "[TT]-fm[tt] or [TT]-ft[tt].",
96 "Note that with the [TT]-dt[tt] option the memory usage and calculation",
97 "time can be reduced."
99 static int fmmax
= 0, ftmax
= 0, nlev
= 81, nr
= 1, fshift
= 0;
100 static real sbin
= 0, rmax
= 2, rbin
= 0.01, mmax
= 0, rint
= 0;
106 "Use [SQRT]t[sqrt] on the matrix axis which binspacing # in [SQRT]ps[sqrt]" },
107 { "-fm", FALSE
, etINT
, { &fmmax
}, "Number of frames in the matrix, 0 is plot all" },
108 { "-rmax", FALSE
, etREAL
, { &rmax
}, "Maximum r in the matrix (nm)" },
109 { "-rbin", FALSE
, etREAL
, { &rbin
}, "Binwidth in the matrix and for [TT]-or[tt] (nm)" },
114 "Maximum density in the matrix, 0 is calculate (1/nm)" },
115 { "-nlevels", FALSE
, etINT
, { &nlev
}, "Number of levels in the matrix" },
116 { "-nr", FALSE
, etINT
, { &nr
}, "Number of curves for the [TT]-or[tt] output" },
117 { "-fr", FALSE
, etINT
, { &fshift
}, "Frame spacing for the [TT]-or[tt] output" },
118 { "-rt", FALSE
, etREAL
, { &rint
}, "Integration limit for the [TT]-ot[tt] output (nm)" },
123 "Number of frames in the [TT]-ot[tt] output, 0 is plot all" }
125 #define NPA asize(pa)
128 { efTRX
, nullptr, nullptr, ffREAD
}, { efTPS
, nullptr, nullptr, ffREAD
},
129 { efNDX
, nullptr, nullptr, ffOPTRD
}, { efXPM
, "-om", "vanhove", ffOPTWR
},
130 { efXVG
, "-or", "vanhove_r", ffOPTWR
}, { efXVG
, "-ot", "vanhove_t", ffOPTWR
}
132 #define NFILE asize(fnm)
134 gmx_output_env_t
* oenv
;
135 const char * matfile
, *otfile
, *orfile
;
138 matrix boxtop
, box
, *sbox
, avbox
, corr
;
139 rvec
* xtop
, *x
, **sx
;
140 int isize
, nalloc
, nallocn
;
144 int nfr
, f
, ff
, i
, m
, mat_nx
= 0, nbin
= 0, bin
, mbin
, fbin
;
145 real
* time
, t
, invbin
= 0, rmax2
= 0, rint2
= 0, d2
;
146 real invsbin
= 0, matmax
, normfac
, dt
, *tickx
, *ticky
;
147 char buf
[STRLEN
], **legend
;
148 real
** mat
= nullptr;
149 int * pt
= nullptr, **pr
= nullptr, *mcount
= nullptr, *tcount
= nullptr, *rcount
= nullptr;
151 t_rgb rlo
= { 1, 1, 1 }, rhi
= { 0, 0, 0 };
153 if (!parse_common_args(&argc
, argv
, PCA_CAN_VIEW
| PCA_CAN_TIME
, NFILE
, fnm
, asize(pa
), pa
,
154 asize(desc
), desc
, 0, nullptr, &oenv
))
159 matfile
= opt2fn_null("-om", NFILE
, fnm
);
160 if (opt2parg_bSet("-fr", NPA
, pa
))
162 orfile
= opt2fn("-or", NFILE
, fnm
);
166 orfile
= opt2fn_null("-or", NFILE
, fnm
);
168 if (opt2parg_bSet("-rt", NPA
, pa
))
170 otfile
= opt2fn("-ot", NFILE
, fnm
);
174 otfile
= opt2fn_null("-ot", NFILE
, fnm
);
177 if (!matfile
&& !otfile
&& !orfile
)
179 fprintf(stderr
, "For output set one (or more) of the output file options\n");
183 read_tps_conf(ftp2fn(efTPS
, NFILE
, fnm
), &top
, &ePBC
, &xtop
, nullptr, boxtop
, FALSE
);
184 get_index(&top
.atoms
, ftp2fn_null(efNDX
, NFILE
, fnm
), 1, &isize
, &index
, &grpname
);
192 read_first_x(oenv
, &status
, ftp2fn(efTRX
, NFILE
, fnm
), &t
, &x
, box
);
199 srenew(time
, nalloc
);
200 srenew(sbox
, nalloc
);
203 GMX_RELEASE_ASSERT(time
!= nullptr, "Memory allocation failure; time array is NULL");
204 GMX_RELEASE_ASSERT(sbox
!= nullptr, "Memory allocation failure; sbox array is NULL");
207 copy_mat(box
, sbox
[nfr
]);
208 /* This assumes that the off-diagonal box elements
209 * are not affected by jumps across the periodic boundaries.
211 m_add(avbox
, box
, avbox
);
212 snew(sx
[nfr
], isize
);
213 for (i
= 0; i
< isize
; i
++)
215 copy_rvec(x
[index
[i
]], sx
[nfr
][i
]);
219 } while (read_next_x(oenv
, status
, &t
, x
, box
));
225 fprintf(stderr
, "Read %d frames\n", nfr
);
227 dt
= (time
[nfr
- 1] - time
[0]) / (nfr
- 1);
228 /* Some ugly rounding to get nice nice times in the output */
229 dt
= std::round(10000.0 * dt
) / 10000.0;
235 if (fmmax
<= 0 || fmmax
>= nfr
)
240 nbin
= gmx::roundToInt(rmax
* invbin
);
247 invsbin
= 1.0 / sbin
;
248 mat_nx
= static_cast<int>(std::sqrt(fmmax
* dt
) * invsbin
+ 1);
251 for (f
= 0; f
< mat_nx
; f
++)
255 rmax2
= gmx::square(nbin
* rbin
);
256 /* Initialize time zero */
257 mat
[0][0] = nfr
* isize
;
281 /* Initialize time zero */
290 msmul(avbox
, 1.0 / nfr
, avbox
);
291 for (f
= 0; f
< nfr
; f
++)
295 fprintf(stderr
, "\rProcessing frame %d", f
);
298 if (ePBC
!= epbcNONE
)
300 /* Scale all the configuration to the average box */
301 gmx::invertBoxMatrix(sbox
[f
], corr
);
302 mmul_ur0(avbox
, corr
, corr
);
303 for (i
= 0; i
< isize
; i
++)
305 mvmul_ur0(corr
, sx
[f
][i
], sx
[f
][i
]);
308 /* Correct for periodic jumps */
309 for (m
= DIM
- 1; m
>= 0; m
--)
311 while (sx
[f
][i
][m
] - sx
[f
- 1][i
][m
] > 0.5 * avbox
[m
][m
])
313 rvec_dec(sx
[f
][i
], avbox
[m
]);
315 while (sx
[f
][i
][m
] - sx
[f
- 1][i
][m
] <= -0.5 * avbox
[m
][m
])
317 rvec_inc(sx
[f
][i
], avbox
[m
]);
323 for (ff
= 0; ff
< f
; ff
++)
326 if (fbin
<= fmmax
|| fbin
<= ftmax
)
334 mbin
= gmx::roundToInt(std::sqrt(fbin
* dt
) * invsbin
);
336 for (i
= 0; i
< isize
; i
++)
338 d2
= distance2(sx
[f
][i
], sx
[ff
][i
]);
339 if (mbin
< mat_nx
&& d2
< rmax2
)
341 bin
= gmx::roundToInt(std::sqrt(d2
) * invbin
);
347 if (fbin
<= ftmax
&& d2
<= rint2
)
364 for (fbin
= 0; fbin
< nr
; fbin
++)
366 ff
= f
- (fbin
+ 1) * fshift
;
369 for (i
= 0; i
< isize
; i
++)
371 d2
= distance2(sx
[f
][i
], sx
[ff
][i
]);
372 bin
= gmx::roundToInt(std::sqrt(d2
) * invbin
);
375 nallocn
= 10 * (bin
/ 10) + 11;
376 for (m
= 0; m
< nr
; m
++)
378 srenew(pr
[m
], nallocn
);
379 for (i
= nalloc
; i
< nallocn
; i
++)
393 fprintf(stderr
, "\n");
398 for (f
= 0; f
< mat_nx
; f
++)
400 normfac
= 1.0 / (mcount
[f
] * isize
* rbin
);
401 for (i
= 0; i
< nbin
; i
++)
403 mat
[f
][i
] *= normfac
;
404 if (mat
[f
][i
] > matmax
&& (f
!= 0 || i
!= 0))
410 fprintf(stdout
, "Value at (0,0): %.3f, maximum of the rest %.3f\n", mat
[0][0], matmax
);
416 for (f
= 0; f
< mat_nx
; f
++)
427 snew(ticky
, nbin
+ 1);
428 for (i
= 0; i
<= nbin
; i
++)
432 fp
= gmx_ffopen(matfile
, "w");
433 write_xpm(fp
, MAT_SPATIAL_Y
, "Van Hove function", "G (1/nm)",
434 sbin
== 0 ? "time (ps)" : "sqrt(time) (ps^1/2)", "r (nm)", mat_nx
, nbin
, tickx
,
435 ticky
, mat
, 0, matmax
, rlo
, rhi
, &nlev
);
441 fp
= xvgropen(orfile
, "Van Hove function", "r (nm)", "G (nm\\S-1\\N)", oenv
);
442 if (output_env_get_print_xvgr_codes(oenv
))
444 fprintf(fp
, "@ subtitle \"for particles in group %s\"\n", grpname
);
447 for (fbin
= 0; fbin
< nr
; fbin
++)
449 sprintf(buf
, "%g ps", (fbin
+ 1) * fshift
* dt
);
450 legend
[fbin
] = gmx_strdup(buf
);
452 xvgr_legend(fp
, nr
, legend
, oenv
);
453 for (i
= 0; i
< nalloc
; i
++)
455 fprintf(fp
, "%g", i
* rbin
);
456 for (fbin
= 0; fbin
< nr
; fbin
++)
459 static_cast<real
>(pr
[fbin
][i
]
460 / (rcount
[fbin
] * isize
* rbin
* (i
== 0 ? 0.5 : 1.0))));
469 sprintf(buf
, "Probability of moving less than %g nm", rint
);
470 fp
= xvgropen(otfile
, buf
, "t (ps)", "", oenv
);
471 if (output_env_get_print_xvgr_codes(oenv
))
473 fprintf(fp
, "@ subtitle \"for particles in group %s\"\n", grpname
);
475 for (f
= 0; f
<= ftmax
; f
++)
477 fprintf(fp
, "%g %g\n", f
* dt
, static_cast<real
>(pt
[f
]) / (tcount
[f
] * isize
));
482 do_view(oenv
, matfile
, nullptr);
483 do_view(oenv
, orfile
, nullptr);
484 do_view(oenv
, otfile
, nullptr);