src/gromacs/gmxana/gmx_rotacf.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   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.
  11  *
  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.
  16  *
  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.
  21  *
  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.
  26  *
  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.
  34  *
  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.
  37  */
  38 #include "gmxpre.h"
  39
  40 #include <cmath>
  41 #include <cstring>
  42
  43 #include "gromacs/commandline/pargs.h"
  44 #include "gromacs/commandline/viewit.h"
  45 #include "gromacs/correlationfunctions/autocorr.h"
  46 #include "gromacs/fileio/trxio.h"
  47 #include "gromacs/gmxana/gmx_ana.h"
  48 #include "gromacs/math/vec.h"
  49 #include "gromacs/pbcutil/rmpbc.h"
  50 #include "gromacs/topology/index.h"
  51 #include "gromacs/topology/topology.h"
  52 #include "gromacs/utility/arraysize.h"
  53 #include "gromacs/utility/fatalerror.h"
  54 #include "gromacs/utility/futil.h"
  55 #include "gromacs/utility/smalloc.h"
  56
  57 int gmx_rotacf(int argc, char* argv[])
  58 {
  59     const char* desc[] = {
  60         "[THISMODULE] calculates the rotational correlation function",
  61         "for molecules. Atom triplets (i,j,k) must be given in the index",
  62         "file, defining two vectors ij and jk. The rotational ACF",
  63         "is calculated as the autocorrelation function of the vector",
  64         "n = ij x jk, i.e. the cross product of the two vectors.",
  65         "Since three atoms span a plane, the order of the three atoms",
  66         "does not matter. Optionally, by invoking the [TT]-d[tt] switch, you can",
  67         "calculate the rotational correlation function for linear molecules",
  68         "by specifying atom pairs (i,j) in the index file.",
  69         "[PAR]",
  70         "EXAMPLES[PAR]",
  71         "[TT]gmx rotacf -P 1 -nparm 2 -fft -n index -o rotacf-x-P1",
  72         "-fa expfit-x-P1 -beginfit 2.5 -endfit 20.0[tt][PAR]",
  73         "This will calculate the rotational correlation function using a first",
  74         "order Legendre polynomial of the angle of a vector defined by the index",
  75         "file. The correlation function will be fitted from 2.5 ps until 20.0 ps",
  76         "to a two-parameter exponential."
  77     };
  78     static gmx_bool bVec = FALSE, bAver = TRUE;
  79
  80     t_pargs pa[] = {
  81         { "-d",
  82           FALSE,
  83           etBOOL,
  84           { &bVec },
  85           "Use index doublets (vectors) for correlation function instead of triplets (planes)" },
  86         { "-aver", FALSE, etBOOL, { &bAver }, "Average over molecules" }
  87     };
  88
  89     t_trxstatus*  status;
  90     int           isize;
  91     int*          index;
  92     char*         grpname;
  93     rvec *        x, *x_s;
  94     matrix        box;
  95     real**        c1;
  96     rvec          xij, xjk, n;
  97     int           i, m, teller, n_alloc, natoms, nvec, ai, aj, ak;
  98     unsigned long mode;
  99     real          t, t0, t1, dt;
 100     gmx_rmpbc_t   gpbc = nullptr;
 101     t_topology*   top;
 102     int           ePBC;
 103     t_filenm      fnm[] = { { efTRX, "-f", nullptr, ffREAD },
 104                        { efTPR, nullptr, nullptr, ffREAD },
 105                        { efNDX, nullptr, nullptr, ffREAD },
 106                        { efXVG, "-o", "rotacf", ffWRITE } };
 107 #define NFILE asize(fnm)
 108     int      npargs;
 109     t_pargs* ppa;
 110
 111     gmx_output_env_t* oenv;
 112
 113     npargs = asize(pa);
 114     ppa    = add_acf_pargs(&npargs, pa);
 115
 116     if (!parse_common_args(&argc, argv, PCA_CAN_VIEW | PCA_CAN_TIME, NFILE, fnm, npargs, ppa,
 117                            asize(desc), desc, 0, nullptr, &oenv))
 118     {
 119         sfree(ppa);
 120         return 0;
 121     }
 122
 123     rd_index(ftp2fn(efNDX, NFILE, fnm), 1, &isize, &index, &grpname);
 124
 125     if (bVec)
 126     {
 127         nvec = isize / 2;
 128     }
 129     else
 130     {
 131         nvec = isize / 3;
 132     }
 133
 134     if (((isize % 3) != 0) && !bVec)
 135     {
 136         gmx_fatal(FARGS,
 137                   "number of index elements not multiple of 3, "
 138                   "these can not be atom triplets\n");
 139     }
 140     if (((isize % 2) != 0) && bVec)
 141     {
 142         gmx_fatal(FARGS,
 143                   "number of index elements not multiple of 2, "
 144                   "these can not be atom doublets\n");
 145     }
 146
 147     top = read_top(ftp2fn(efTPR, NFILE, fnm), &ePBC);
 148
 149     snew(c1, nvec);
 150     for (i = 0; (i < nvec); i++)
 151     {
 152         c1[i] = nullptr;
 153     }
 154     n_alloc = 0;
 155
 156     natoms = read_first_x(oenv, &status, ftp2fn(efTRX, NFILE, fnm), &t, &x, box);
 157     snew(x_s, natoms);
 158
 159     gpbc = gmx_rmpbc_init(&(top->idef), ePBC, natoms);
 160
 161     /* Start the loop over frames */
 162     t0     = t;
 163     teller = 0;
 164     do
 165     {
 166         if (teller >= n_alloc)
 167         {
 168             n_alloc += 100;
 169             for (i = 0; (i < nvec); i++)
 170             {
 171                 srenew(c1[i], DIM * n_alloc);
 172             }
 173         }
 174         t1 = t;
 175
 176         /* Remove periodicity */
 177         gmx_rmpbc_copy(gpbc, natoms, box, x, x_s);
 178
 179         /* Compute crossproducts for all vectors, if triplets.
 180          * else, just get the vectors in case of doublets.
 181          */
 182         if (!bVec)
 183         {
 184             for (i = 0; (i < nvec); i++)
 185             {
 186                 ai = index[3 * i];
 187                 aj = index[3 * i + 1];
 188                 ak = index[3 * i + 2];
 189                 rvec_sub(x_s[ai], x_s[aj], xij);
 190                 rvec_sub(x_s[aj], x_s[ak], xjk);
 191                 cprod(xij, xjk, n);
 192                 for (m = 0; (m < DIM); m++)
 193                 {
 194                     c1[i][DIM * teller + m] = n[m];
 195                 }
 196             }
 197         }
 198         else
 199         {
 200             for (i = 0; (i < nvec); i++)
 201             {
 202                 ai = index[2 * i];
 203                 aj = index[2 * i + 1];
 204                 rvec_sub(x_s[ai], x_s[aj], n);
 205                 for (m = 0; (m < DIM); m++)
 206                 {
 207                     c1[i][DIM * teller + m] = n[m];
 208                 }
 209             }
 210         }
 211         /* Increment loop counter */
 212         teller++;
 213     } while (read_next_x(oenv, status, &t, x, box));
 214     close_trx(status);
 215     fprintf(stderr, "\nDone with trajectory\n");
 216
 217     gmx_rmpbc_done(gpbc);
 218
 219
 220     /* Autocorrelation function */
 221     if (teller < 2)
 222     {
 223         fprintf(stderr, "Not enough frames for correlation function\n");
 224     }
 225     else
 226     {
 227         dt = (t1 - t0) / (static_cast<real>(teller - 1));
 228
 229         mode = eacVector;
 230
 231         do_autocorr(ftp2fn(efXVG, NFILE, fnm), oenv, "Rotational Correlation Function", teller,
 232                     nvec, c1, dt, mode, bAver);
 233     }
 234
 235     do_view(oenv, ftp2fn(efXVG, NFILE, fnm), nullptr);
 236
 237     return 0;
 238 }