gmx_vanhove.cpp 16.1 KB
Newer Older
hess's avatar
hess committed
1
/*
2
 * This file is part of the GROMACS molecular simulation package.
3
 *
4
5
6
 * Copyright 1991- The GROMACS Authors
 * and the project initiators Erik Lindahl, Berk Hess and David van der Spoel.
 * Consult the AUTHORS/COPYING files and https://www.gromacs.org for details.
7
8
9
10
 *
 * GROMACS is free software; you can redistribute it and/or
 * modify it under the terms of the GNU Lesser General Public License
 * as published by the Free Software Foundation; either version 2.1
hess's avatar
hess committed
11
 * of the License, or (at your option) any later version.
12
 *
13
14
15
16
 * GROMACS is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * Lesser General Public License for more details.
17
 *
18
19
 * You should have received a copy of the GNU Lesser General Public
 * License along with GROMACS; if not, see
20
 * https://www.gnu.org/licenses, or write to the Free Software Foundation,
21
 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
22
 *
23
24
25
26
27
28
 * If you want to redistribute modifications to GROMACS, please
 * consider that scientific software is very special. Version
 * control is crucial - bugs must be traceable. We will be happy to
 * consider code for inclusion in the official distribution, but
 * derived work must not be called official GROMACS. Details are found
 * in the README & COPYING files - if they are missing, get the
29
 * official version at https://www.gromacs.org.
30
 *
31
 * To help us fund GROMACS development, we humbly ask that you cite
32
 * the research papers on the package. Check out https://www.gromacs.org.
hess's avatar
hess committed
33
 */
34
35
#include "gmxpre.h"

36
37
38
#include <cmath>
#include <cstdlib>
#include <cstring>
hess's avatar
hess committed
39

40
#include "gromacs/commandline/pargs.h"
Teemu Murtola's avatar
Teemu Murtola committed
41
#include "gromacs/commandline/viewit.h"
42
#include "gromacs/fileio/confio.h"
43
#include "gromacs/fileio/matio.h"
Mark Abraham's avatar
Mark Abraham committed
44
#include "gromacs/fileio/trxio.h"
45
46
#include "gromacs/fileio/xvgr.h"
#include "gromacs/gmxana/gmx_ana.h"
Mark Abraham's avatar
Mark Abraham committed
47
#include "gromacs/math/boxmatrix.h"
Erik Lindahl's avatar
Erik Lindahl committed
48
#include "gromacs/math/functions.h"
Teemu Murtola's avatar
Teemu Murtola committed
49
#include "gromacs/math/vec.h"
Berk Hess's avatar
Berk Hess committed
50
#include "gromacs/pbcutil/pbc.h"
51
#include "gromacs/topology/index.h"
52
#include "gromacs/topology/topology.h"
53
#include "gromacs/utility/arraysize.h"
Teemu Murtola's avatar
Teemu Murtola committed
54
#include "gromacs/utility/cstringutil.h"
55
#include "gromacs/utility/futil.h"
56
#include "gromacs/utility/gmxassert.h"
57
#include "gromacs/utility/smalloc.h"
58
#include "gromacs/utility/stringutil.h"
59

hess's avatar
hess committed
60

Paul Bauer's avatar
Paul Bauer committed
61
int gmx_vanhove(int argc, char* argv[])
hess's avatar
hess committed
62
{
Paul Bauer's avatar
Paul Bauer committed
63
    const char* desc[] = {
64
        "[THISMODULE] computes the Van Hove correlation function.",
65
66
        "The Van Hove G(r,t) is the probability that a particle that is at r[SUB]0[sub]",
        "at time zero can be found at position r[SUB]0[sub]+r at time t.",
67
        "[THISMODULE] determines G not for a vector r, but for the length of r.",
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
        "Thus it gives the probability that a particle moves a distance of r",
        "in time t.",
        "Jumps across the periodic boundaries are removed.",
        "Corrections are made for scaling due to isotropic",
        "or anisotropic pressure coupling.",
        "[PAR]",
        "With option [TT]-om[tt] the whole matrix can be written as a function",
        "of t and r or as a function of [SQRT]t[sqrt] and r (option [TT]-sqrt[tt]).",
        "[PAR]",
        "With option [TT]-or[tt] the Van Hove function is plotted for one",
        "or more values of t. Option [TT]-nr[tt] sets the number of times,",
        "option [TT]-fr[tt] the number spacing between the times.",
        "The binwidth is set with option [TT]-rbin[tt]. The number of bins",
        "is determined automatically.",
        "[PAR]",
        "With option [TT]-ot[tt] the integral up to a certain distance",
        "(option [TT]-rt[tt]) is plotted as a function of time.",
        "[PAR]",
        "For all frames that are read the coordinates of the selected particles",
        "are stored in memory. Therefore the program may use a lot of memory.",
        "For options [TT]-om[tt] and [TT]-ot[tt] the program may be slow.",
        "This is because the calculation scales as the number of frames times",
        "[TT]-fm[tt] or [TT]-ft[tt].",
        "Note that with the [TT]-dt[tt] option the memory usage and calculation",
        "time can be reduced."
    };
    static int  fmmax = 0, ftmax = 0, nlev = 81, nr = 1, fshift = 0;
Paul Bauer's avatar
Paul Bauer committed
95
96
97
98
99
100
    static real sbin = 0, rmax = 2, rbin = 0.01, mmax = 0, rint = 0;
    t_pargs     pa[] = {
        { "-sqrt",
          FALSE,
          etREAL,
          { &sbin },
101
          "Use [SQRT]t[sqrt] on the matrix axis which binspacing # in [SQRT]ps[sqrt]" },
Paul Bauer's avatar
Paul Bauer committed
102
103
104
105
106
107
108
        { "-fm", FALSE, etINT, { &fmmax }, "Number of frames in the matrix, 0 is plot all" },
        { "-rmax", FALSE, etREAL, { &rmax }, "Maximum r in the matrix (nm)" },
        { "-rbin", FALSE, etREAL, { &rbin }, "Binwidth in the matrix and for [TT]-or[tt] (nm)" },
        { "-mmax",
          FALSE,
          etREAL,
          { &mmax },
109
          "Maximum density in the matrix, 0 is calculate (1/nm)" },
Paul Bauer's avatar
Paul Bauer committed
110
111
112
113
114
115
116
117
        { "-nlevels", FALSE, etINT, { &nlev }, "Number of levels in the matrix" },
        { "-nr", FALSE, etINT, { &nr }, "Number of curves for the [TT]-or[tt] output" },
        { "-fr", FALSE, etINT, { &fshift }, "Frame spacing for the [TT]-or[tt] output" },
        { "-rt", FALSE, etREAL, { &rint }, "Integration limit for the [TT]-ot[tt] output (nm)" },
        { "-ft",
          FALSE,
          etINT,
          { &ftmax },
118
119
          "Number of frames in the [TT]-ot[tt] output, 0 is plot all" }
    };
hess's avatar
hess committed
120
121
#define NPA asize(pa)

122
    t_filenm fnm[] = {
Paul Bauer's avatar
Paul Bauer committed
123
124
125
        { efTRX, nullptr, nullptr, ffREAD },    { efTPS, nullptr, nullptr, ffREAD },
        { efNDX, nullptr, nullptr, ffOPTRD },   { efXPM, "-om", "vanhove", ffOPTWR },
        { efXVG, "-or", "vanhove_r", ffOPTWR }, { efXVG, "-ot", "vanhove_t", ffOPTWR }
126
    };
hess's avatar
hess committed
127
128
#define NFILE asize(fnm)

129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
    gmx_output_env_t*        oenv;
    const char *             matfile, *otfile, *orfile;
    t_topology               top;
    PbcType                  pbcType;
    matrix                   boxtop, box, *sbox, avbox, corr;
    rvec *                   xtop, *x, **sx;
    int                      isize, nalloc, nallocn;
    t_trxstatus*             status;
    int*                     index;
    char*                    grpname;
    int                      nfr, f, ff, i, m, mat_nx = 0, nbin = 0, bin, mbin, fbin;
    real *                   time, t, invbin = 0, rmax2 = 0, rint2 = 0, d2;
    real                     invsbin = 0, matmax, normfac, dt, *tickx, *ticky;
    std::vector<std::string> legend;
    real**                   mat = nullptr;
Paul Bauer's avatar
Paul Bauer committed
144
145
146
    int * pt = nullptr, **pr = nullptr, *mcount = nullptr, *tcount = nullptr, *rcount = nullptr;
    FILE* fp;
    t_rgb rlo = { 1, 1, 1 }, rhi = { 0, 0, 0 };
hess's avatar
hess committed
147

148
149
    if (!parse_common_args(
                &argc, argv, PCA_CAN_VIEW | PCA_CAN_TIME, NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, nullptr, &oenv))
150
151
152
    {
        return 0;
    }
hess's avatar
hess committed
153

154
155
156
    matfile = opt2fn_null("-om", NFILE, fnm);
    if (opt2parg_bSet("-fr", NPA, pa))
    {
Paul Bauer's avatar
Paul Bauer committed
157
        orfile = opt2fn("-or", NFILE, fnm);
158
159
160
    }
    else
    {
Paul Bauer's avatar
Paul Bauer committed
161
        orfile = opt2fn_null("-or", NFILE, fnm);
162
163
164
    }
    if (opt2parg_bSet("-rt", NPA, pa))
    {
Paul Bauer's avatar
Paul Bauer committed
165
        otfile = opt2fn("-ot", NFILE, fnm);
166
167
168
    }
    else
    {
Paul Bauer's avatar
Paul Bauer committed
169
        otfile = opt2fn_null("-ot", NFILE, fnm);
hess's avatar
hess committed
170
171
    }

172
173
    if (!matfile && !otfile && !orfile)
    {
Paul Bauer's avatar
Paul Bauer committed
174
        fprintf(stderr, "For output set one (or more) of the output file options\n");
175
176
177
        exit(0);
    }

178
    read_tps_conf(ftp2fn(efTPS, NFILE, fnm), &top, &pbcType, &xtop, nullptr, boxtop, FALSE);
179
    get_index(&top.atoms, ftp2fn_null(efNDX, NFILE, fnm), 1, &isize, &index, &grpname);
hess's avatar
hess committed
180

181
    nalloc = 0;
Roland Schulz's avatar
Roland Schulz committed
182
183
184
    time   = nullptr;
    sbox   = nullptr;
    sx     = nullptr;
185
    clear_mat(avbox);
hess's avatar
hess committed
186

187
    read_first_x(oenv, &status, ftp2fn(efTRX, NFILE, fnm), &t, &x, box);
Paul Bauer's avatar
Paul Bauer committed
188
    nfr = 0;
189
190
191
192
193
194
195
196
197
    do
    {
        if (nfr >= nalloc)
        {
            nalloc += 100;
            srenew(time, nalloc);
            srenew(sbox, nalloc);
            srenew(sx, nalloc);
        }
Roland Schulz's avatar
Roland Schulz committed
198
199
        GMX_RELEASE_ASSERT(time != nullptr, "Memory allocation failure; time array is NULL");
        GMX_RELEASE_ASSERT(sbox != nullptr, "Memory allocation failure; sbox array is NULL");
hess's avatar
hess committed
200

201
202
203
204
205
206
207
208
209
210
211
212
213
        time[nfr] = t;
        copy_mat(box, sbox[nfr]);
        /* This assumes that the off-diagonal box elements
         * are not affected by jumps across the periodic boundaries.
         */
        m_add(avbox, box, avbox);
        snew(sx[nfr], isize);
        for (i = 0; i < isize; i++)
        {
            copy_rvec(x[index[i]], sx[nfr][i]);
        }

        nfr++;
Paul Bauer's avatar
Paul Bauer committed
214
    } while (read_next_x(oenv, status, &t, x, box));
215
216
217

    /* clean up */
    sfree(x);
218
    close_trx(status);
219
220
221

    fprintf(stderr, "Read %d frames\n", nfr);

Paul Bauer's avatar
Paul Bauer committed
222
    dt = (time[nfr - 1] - time[0]) / (nfr - 1);
223
    /* Some ugly rounding to get nice nice times in the output */
Paul Bauer's avatar
Paul Bauer committed
224
    dt = std::round(10000.0 * dt) / 10000.0;
hess's avatar
hess committed
225

Paul Bauer's avatar
Paul Bauer committed
226
    invbin = 1.0 / rbin;
227
228
229
230
231
232
233
234

    if (matfile)
    {
        if (fmmax <= 0 || fmmax >= nfr)
        {
            fmmax = nfr - 1;
        }
        snew(mcount, fmmax);
Paul Bauer's avatar
Paul Bauer committed
235
        nbin = gmx::roundToInt(rmax * invbin);
236
237
238
239
240
241
        if (sbin == 0)
        {
            mat_nx = fmmax + 1;
        }
        else
        {
Paul Bauer's avatar
Paul Bauer committed
242
243
            invsbin = 1.0 / sbin;
            mat_nx  = static_cast<int>(std::sqrt(fmmax * dt) * invsbin + 1);
244
245
246
247
248
249
        }
        snew(mat, mat_nx);
        for (f = 0; f < mat_nx; f++)
        {
            snew(mat[f], nbin);
        }
Paul Bauer's avatar
Paul Bauer committed
250
        rmax2 = gmx::square(nbin * rbin);
251
        /* Initialize time zero */
Paul Bauer's avatar
Paul Bauer committed
252
        mat[0][0] = nfr * isize;
253
254
255
256
257
        mcount[0] += nfr;
    }
    else
    {
        fmmax = 0;
hess's avatar
hess committed
258
    }
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274

    if (orfile)
    {
        snew(pr, nr);
        nalloc = 0;
        snew(rcount, nr);
    }

    if (otfile)
    {
        if (ftmax <= 0)
        {
            ftmax = nfr - 1;
        }
        snew(tcount, ftmax);
        snew(pt, nfr);
Paul Bauer's avatar
Paul Bauer committed
275
        rint2 = rint * rint;
276
        /* Initialize time zero */
Paul Bauer's avatar
Paul Bauer committed
277
        pt[0] = nfr * isize;
278
        tcount[0] += nfr;
hess's avatar
hess committed
279
    }
280
281
282
    else
    {
        ftmax = 0;
hess's avatar
hess committed
283
    }
284

Paul Bauer's avatar
Paul Bauer committed
285
    msmul(avbox, 1.0 / nfr, avbox);
286
287
288
289
290
    for (f = 0; f < nfr; f++)
    {
        if (f % 100 == 0)
        {
            fprintf(stderr, "\rProcessing frame %d", f);
291
            fflush(stderr);
292
        }
293
        if (pbcType != PbcType::No)
294
        {
Berk Hess's avatar
Berk Hess committed
295
296
297
298
            /* Scale all the configuration to the average box */
            gmx::invertBoxMatrix(sbox[f], corr);
            mmul_ur0(avbox, corr, corr);
            for (i = 0; i < isize; i++)
299
            {
Berk Hess's avatar
Berk Hess committed
300
301
                mvmul_ur0(corr, sx[f][i], sx[f][i]);
                if (f > 0)
302
                {
Berk Hess's avatar
Berk Hess committed
303
                    /* Correct for periodic jumps */
Paul Bauer's avatar
Paul Bauer committed
304
                    for (m = DIM - 1; m >= 0; m--)
305
                    {
Paul Bauer's avatar
Paul Bauer committed
306
                        while (sx[f][i][m] - sx[f - 1][i][m] > 0.5 * avbox[m][m])
Berk Hess's avatar
Berk Hess committed
307
308
309
                        {
                            rvec_dec(sx[f][i], avbox[m]);
                        }
Paul Bauer's avatar
Paul Bauer committed
310
                        while (sx[f][i][m] - sx[f - 1][i][m] <= -0.5 * avbox[m][m])
Berk Hess's avatar
Berk Hess committed
311
312
313
                        {
                            rvec_inc(sx[f][i], avbox[m]);
                        }
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
                    }
                }
            }
        }
        for (ff = 0; ff < f; ff++)
        {
            fbin = f - ff;
            if (fbin <= fmmax || fbin <= ftmax)
            {
                if (sbin == 0)
                {
                    mbin = fbin;
                }
                else
                {
Paul Bauer's avatar
Paul Bauer committed
329
                    mbin = gmx::roundToInt(std::sqrt(fbin * dt) * invsbin);
330
331
332
333
334
335
                }
                for (i = 0; i < isize; i++)
                {
                    d2 = distance2(sx[f][i], sx[ff][i]);
                    if (mbin < mat_nx && d2 < rmax2)
                    {
Paul Bauer's avatar
Paul Bauer committed
336
                        bin = gmx::roundToInt(std::sqrt(d2) * invbin);
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
                        if (bin < nbin)
                        {
                            mat[mbin][bin] += 1;
                        }
                    }
                    if (fbin <= ftmax && d2 <= rint2)
                    {
                        pt[fbin]++;
                    }
                }
                if (matfile)
                {
                    mcount[mbin]++;
                }
                if (otfile)
                {
                    tcount[fbin]++;
                }
            }
        }
        if (orfile)
        {
            for (fbin = 0; fbin < nr; fbin++)
            {
Paul Bauer's avatar
Paul Bauer committed
361
                ff = f - (fbin + 1) * fshift;
362
363
364
365
366
                if (ff >= 0)
                {
                    for (i = 0; i < isize; i++)
                    {
                        d2  = distance2(sx[f][i], sx[ff][i]);
Paul Bauer's avatar
Paul Bauer committed
367
                        bin = gmx::roundToInt(std::sqrt(d2) * invbin);
368
369
                        if (bin >= nalloc)
                        {
Paul Bauer's avatar
Paul Bauer committed
370
                            nallocn = 10 * (bin / 10) + 11;
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
                            for (m = 0; m < nr; m++)
                            {
                                srenew(pr[m], nallocn);
                                for (i = nalloc; i < nallocn; i++)
                                {
                                    pr[m][i] = 0;
                                }
                            }
                            nalloc = nallocn;
                        }
                        pr[fbin][bin]++;
                    }
                    rcount[fbin]++;
                }
            }
        }
hess's avatar
hess committed
387
    }
388
389
390
391
392
393
394
    fprintf(stderr, "\n");

    if (matfile)
    {
        matmax = 0;
        for (f = 0; f < mat_nx; f++)
        {
Paul Bauer's avatar
Paul Bauer committed
395
            normfac = 1.0 / (mcount[f] * isize * rbin);
396
397
398
399
400
401
402
403
404
            for (i = 0; i < nbin; i++)
            {
                mat[f][i] *= normfac;
                if (mat[f][i] > matmax && (f != 0 || i != 0))
                {
                    matmax = mat[f][i];
                }
            }
        }
Paul Bauer's avatar
Paul Bauer committed
405
        fprintf(stdout, "Value at (0,0): %.3f, maximum of the rest %.3f\n", mat[0][0], matmax);
406
407
408
409
410
411
412
413
414
        if (mmax > 0)
        {
            matmax = mmax;
        }
        snew(tickx, mat_nx);
        for (f = 0; f < mat_nx; f++)
        {
            if (sbin == 0)
            {
Paul Bauer's avatar
Paul Bauer committed
415
                tickx[f] = f * dt;
416
417
418
            }
            else
            {
Paul Bauer's avatar
Paul Bauer committed
419
                tickx[f] = f * sbin;
420
421
            }
        }
Paul Bauer's avatar
Paul Bauer committed
422
        snew(ticky, nbin + 1);
423
424
        for (i = 0; i <= nbin; i++)
        {
Paul Bauer's avatar
Paul Bauer committed
425
            ticky[i] = i * rbin;
426
        }
427
        fp = gmx_ffopen(matfile, "w");
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
        write_xpm(fp,
                  MAT_SPATIAL_Y,
                  "Van Hove function",
                  "G (1/nm)",
                  sbin == 0 ? "time (ps)" : "sqrt(time) (ps^1/2)",
                  "r (nm)",
                  mat_nx,
                  nbin,
                  tickx,
                  ticky,
                  mat,
                  0,
                  matmax,
                  rlo,
                  rhi,
                  &nlev);
444
        gmx_ffclose(fp);
hess's avatar
hess committed
445
    }
446
447
448
449

    if (orfile)
    {
        fp = xvgropen(orfile, "Van Hove function", "r (nm)", "G (nm\\S-1\\N)", oenv);
450
        if (output_env_get_print_xvgr_codes(oenv))
451
452
453
        {
            fprintf(fp, "@ subtitle \"for particles in group %s\"\n", grpname);
        }
454
455
        for (fbin = 0; fbin < nr; fbin++)
        {
456
            legend.emplace_back(gmx::formatString("%g ps", (fbin + 1) * fshift * dt));
457
        }
458
        xvgrLegend(fp, legend, oenv);
459
460
        for (i = 0; i < nalloc; i++)
        {
Paul Bauer's avatar
Paul Bauer committed
461
            fprintf(fp, "%g", i * rbin);
462
463
            for (fbin = 0; fbin < nr; fbin++)
            {
464
465
                fprintf(fp,
                        " %g",
Paul Bauer's avatar
Paul Bauer committed
466
467
                        static_cast<real>(pr[fbin][i]
                                          / (rcount[fbin] * isize * rbin * (i == 0 ? 0.5 : 1.0))));
468
469
470
            }
            fprintf(fp, "\n");
        }
471
        xvgrclose(fp);
hess's avatar
hess committed
472
    }
473
474
475

    if (otfile)
    {
476
477
        auto buf = gmx::formatString("Probability of moving less than %g nm", rint);
        fp       = xvgropen(otfile, buf.c_str(), "t (ps)", "", oenv);
478
        if (output_env_get_print_xvgr_codes(oenv))
479
480
481
        {
            fprintf(fp, "@ subtitle \"for particles in group %s\"\n", grpname);
        }
482
483
        for (f = 0; f <= ftmax; f++)
        {
Paul Bauer's avatar
Paul Bauer committed
484
            fprintf(fp, "%g %g\n", f * dt, static_cast<real>(pt[f]) / (tcount[f] * isize));
485
        }
486
        xvgrclose(fp);
hess's avatar
hess committed
487
488
    }

Roland Schulz's avatar
Roland Schulz committed
489
490
491
    do_view(oenv, matfile, nullptr);
    do_view(oenv, orfile, nullptr);
    do_view(oenv, otfile, nullptr);
492
493

    return 0;
hess's avatar
hess committed
494
}