src/gromacs/gmxana/gmx_polystat.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,2017,2018,2019, by the GROMACS development team, led by
   7  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   8  * and including many others, as listed in the AUTHORS file in the
   9  * top-level source directory and at http://www.gromacs.org.
  10  *
  11  * GROMACS is free software; you can redistribute it and/or
  12  * modify it under the terms of the GNU Lesser General Public License
  13  * as published by the Free Software Foundation; either version 2.1
  14  * of the License, or (at your option) any later version.
  15  *
  16  * GROMACS is distributed in the hope that it will be useful,
  17  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  18  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19  * Lesser General Public License for more details.
  20  *
  21  * You should have received a copy of the GNU Lesser General Public
  22  * License along with GROMACS; if not, see
  23  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  24  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  25  *
  26  * If you want to redistribute modifications to GROMACS, please
  27  * consider that scientific software is very special. Version
  28  * control is crucial - bugs must be traceable. We will be happy to
  29  * consider code for inclusion in the official distribution, but
  30  * derived work must not be called official GROMACS. Details are found
  31  * in the README & COPYING files - if they are missing, get the
  32  * official version at http://www.gromacs.org.
  33  *
  34  * To help us fund GROMACS development, we humbly ask that you cite
  35  * the research papers on the package. Check out http://www.gromacs.org.
  36  */
  37 #include "gmxpre.h"
  38
  39 #include <cmath>
  40 #include <cstring>
  41
  42 #include <algorithm>
  43
  44 #include "gromacs/commandline/pargs.h"
  45 #include "gromacs/commandline/viewit.h"
  46 #include "gromacs/fileio/tpxio.h"
  47 #include "gromacs/fileio/trxio.h"
  48 #include "gromacs/fileio/xvgr.h"
  49 #include "gromacs/gmxana/gmx_ana.h"
  50 #include "gromacs/linearalgebra/nrjac.h"
  51 #include "gromacs/math/vec.h"
  52 #include "gromacs/pbcutil/rmpbc.h"
  53 #include "gromacs/topology/index.h"
  54 #include "gromacs/topology/topology.h"
  55 #include "gromacs/utility/arraysize.h"
  56 #include "gromacs/utility/cstringutil.h"
  57 #include "gromacs/utility/futil.h"
  58 #include "gromacs/utility/smalloc.h"
  59
  60 static void gyro_eigen(double** gyr, double* eig, double** eigv, int* ord)
  61 {
  62     int nrot, d;
  63
  64     jacobi(gyr, DIM, eig, eigv, &nrot);
  65     /* Order the eigenvalues */
  66     ord[0] = 0;
  67     ord[2] = 2;
  68     for (d = 0; d < DIM; d++)
  69     {
  70         if (eig[d] > eig[ord[0]])
  71         {
  72             ord[0] = d;
  73         }
  74         if (eig[d] < eig[ord[2]])
  75         {
  76             ord[2] = d;
  77         }
  78     }
  79     for (d = 0; d < DIM; d++)
  80     {
  81         if (ord[0] != d && ord[2] != d)
  82         {
  83             ord[1] = d;
  84         }
  85     }
  86 }
  87
  88 /* Calculate mean square internal distances (Auhl et al., JCP 119, 12718) */
  89 static void calc_int_dist(double* intd, rvec* x, int i0, int i1)
  90 {
  91     int       ii;
  92     const int ml = i1 - i0 + 1; /* Number of beads in molecule. */
  93     int       bd;               /* Distance between beads */
  94     double    d;
  95
  96     for (bd = 1; bd < ml; bd++)
  97     {
  98         d = 0.;
  99         for (ii = i0; ii <= i1 - bd; ii++)
 100         {
 101             d += distance2(x[ii], x[ii + bd]);
 102         }
 103         d /= ml - bd;
 104         intd[bd - 1] += d;
 105     }
 106 }
 107
 108 int gmx_polystat(int argc, char* argv[])
 109 {
 110     const char* desc[] = {
 111         "[THISMODULE] plots static properties of polymers as a function of time",
 112         "and prints the average.[PAR]",
 113         "By default it determines the average end-to-end distance and radii",
 114         "of gyration of polymers. It asks for an index group and split this",
 115         "into molecules. The end-to-end distance is then determined using",
 116         "the first and the last atom in the index group for each molecules.",
 117         "For the radius of gyration the total and the three principal components",
 118         "for the average gyration tensor are written.",
 119         "With option [TT]-v[tt] the eigenvectors are written.",
 120         "With option [TT]-pc[tt] also the average eigenvalues of the individual",
 121         "gyration tensors are written.",
 122         "With option [TT]-i[tt] the mean square internal distances are",
 123         "written.[PAR]",
 124         "With option [TT]-p[tt] the persistence length is determined.",
 125         "The chosen index group should consist of atoms that are",
 126         "consecutively bonded in the polymer mainchains.",
 127         "The persistence length is then determined from the cosine of",
 128         "the angles between bonds with an index difference that is even,",
 129         "the odd pairs are not used, because straight polymer backbones",
 130         "are usually all trans and therefore only every second bond aligns.",
 131         "The persistence length is defined as number of bonds where",
 132         "the average cos reaches a value of 1/e. This point is determined",
 133         "by a linear interpolation of [LOG]<cos>[log]."
 134     };
 135     static gmx_bool bMW = TRUE, bPC = FALSE;
 136     t_pargs         pa[] = {
 137         { "-mw", FALSE, etBOOL, { &bMW }, "Use the mass weighting for radii of gyration" },
 138         { "-pc", FALSE, etBOOL, { &bPC }, "Plot average eigenvalues" }
 139     };
 140
 141     t_filenm fnm[] = { { efTPR, nullptr, nullptr, ffREAD },  { efTRX, "-f", nullptr, ffREAD },
 142                        { efNDX, nullptr, nullptr, ffOPTRD }, { efXVG, "-o", "polystat", ffWRITE },
 143                        { efXVG, "-v", "polyvec", ffOPTWR },  { efXVG, "-p", "persist", ffOPTWR },
 144                        { efXVG, "-i", "intdist", ffOPTWR } };
 145 #define NFILE asize(fnm)
 146
 147     t_topology*       top;
 148     gmx_output_env_t* oenv;
 149     int               ePBC;
 150     int               isize, *index, nmol, *molind, mol, nat_min = 0, nat_max = 0;
 151     char*             grpname;
 152     t_trxstatus*      status;
 153     real              t;
 154     rvec *            x, *bond = nullptr;
 155     matrix            box;
 156     int               natoms, i, j, frame, ind0, ind1, a, d, d2, ord[DIM] = { 0 };
 157     dvec              cm, sum_eig = { 0, 0, 0 };
 158     double **         gyr, **gyr_all, eig[DIM], **eigv;
 159     double            sum_eed2, sum_eed2_tot, sum_gyro, sum_gyro_tot, sum_pers_tot;
 160     int*              ninp    = nullptr;
 161     double *          sum_inp = nullptr, pers;
 162     double *          intd, ymax, ymin;
 163     double            mmol, m;
 164     char              title[STRLEN];
 165     FILE *            out, *outv, *outp, *outi;
 166     const char*       leg[8] = { "end to end",      "<R\\sg\\N>",      "<R\\sg\\N> eig1",
 167                            "<R\\sg\\N> eig2", "<R\\sg\\N> eig3", "<R\\sg\\N eig1>",
 168                            "<R\\sg\\N eig2>", "<R\\sg\\N eig3>" };
 169     char **           legp, buf[STRLEN];
 170     gmx_rmpbc_t       gpbc = nullptr;
 171
 172     if (!parse_common_args(&argc, argv, PCA_CAN_VIEW | PCA_CAN_TIME | PCA_TIME_UNIT, NFILE, fnm,
 173                            asize(pa), pa, asize(desc), desc, 0, nullptr, &oenv))
 174     {
 175         return 0;
 176     }
 177
 178     snew(top, 1);
 179     ePBC = read_tpx_top(ftp2fn(efTPR, NFILE, fnm), nullptr, box, &natoms, nullptr, nullptr, top);
 180
 181     fprintf(stderr, "Select a group of polymer mainchain atoms:\n");
 182     get_index(&top->atoms, ftp2fn_null(efNDX, NFILE, fnm), 1, &isize, &index, &grpname);
 183
 184     snew(molind, top->mols.nr + 1);
 185     nmol = 0;
 186     mol  = -1;
 187     for (i = 0; i < isize; i++)
 188     {
 189         if (i == 0 || index[i] >= top->mols.index[mol + 1])
 190         {
 191             molind[nmol++] = i;
 192             do
 193             {
 194                 mol++;
 195             } while (index[i] >= top->mols.index[mol + 1]);
 196         }
 197     }
 198     molind[nmol] = i;
 199     nat_min      = top->atoms.nr;
 200     nat_max      = 0;
 201     for (mol = 0; mol < nmol; mol++)
 202     {
 203         nat_min = std::min(nat_min, molind[mol + 1] - molind[mol]);
 204         nat_max = std::max(nat_max, molind[mol + 1] - molind[mol]);
 205     }
 206     fprintf(stderr, "Group %s consists of %d molecules\n", grpname, nmol);
 207     fprintf(stderr, "Group size per molecule, min: %d atoms, max %d atoms\n", nat_min, nat_max);
 208
 209     sprintf(title, "Size of %d polymers", nmol);
 210     out = xvgropen(opt2fn("-o", NFILE, fnm), title, output_env_get_xvgr_tlabel(oenv), "(nm)", oenv);
 211     xvgr_legend(out, bPC ? 8 : 5, leg, oenv);
 212
 213     if (opt2bSet("-v", NFILE, fnm))
 214     {
 215         outv = xvgropen(opt2fn("-v", NFILE, fnm), "Principal components",
 216                         output_env_get_xvgr_tlabel(oenv), "(nm)", oenv);
 217         snew(legp, DIM * DIM);
 218         for (d = 0; d < DIM; d++)
 219         {
 220             for (d2 = 0; d2 < DIM; d2++)
 221             {
 222                 sprintf(buf, "eig%d %c", d + 1, 'x' + d2);
 223                 legp[d * DIM + d2] = gmx_strdup(buf);
 224             }
 225         }
 226         xvgr_legend(outv, DIM * DIM, legp, oenv);
 227     }
 228     else
 229     {
 230         outv = nullptr;
 231     }
 232
 233     if (opt2bSet("-p", NFILE, fnm))
 234     {
 235         outp = xvgropen(opt2fn("-p", NFILE, fnm), "Persistence length",
 236                         output_env_get_xvgr_tlabel(oenv), "bonds", oenv);
 237         snew(bond, nat_max - 1);
 238         snew(sum_inp, nat_min / 2);
 239         snew(ninp, nat_min / 2);
 240     }
 241     else
 242     {
 243         outp = nullptr;
 244     }
 245
 246     if (opt2bSet("-i", NFILE, fnm))
 247     {
 248         outi = xvgropen(opt2fn("-i", NFILE, fnm), "Internal distances", "n",
 249                         "<R\\S2\\N(n)>/n (nm\\S2\\N)", oenv);
 250         i    = index[molind[1] - 1] - index[molind[0]]; /* Length of polymer -1 */
 251         snew(intd, i);
 252     }
 253     else
 254     {
 255         intd = nullptr;
 256         outi = nullptr;
 257     }
 258
 259     natoms = read_first_x(oenv, &status, ftp2fn(efTRX, NFILE, fnm), &t, &x, box);
 260
 261     snew(gyr, DIM);
 262     snew(gyr_all, DIM);
 263     snew(eigv, DIM);
 264     for (d = 0; d < DIM; d++)
 265     {
 266         snew(gyr[d], DIM);
 267         snew(gyr_all[d], DIM);
 268         snew(eigv[d], DIM);
 269     }
 270
 271     frame        = 0;
 272     sum_eed2_tot = 0;
 273     sum_gyro_tot = 0;
 274     sum_pers_tot = 0;
 275
 276     gpbc = gmx_rmpbc_init(&top->idef, ePBC, natoms);
 277
 278     do
 279     {
 280         gmx_rmpbc(gpbc, natoms, box, x);
 281
 282         sum_eed2 = 0;
 283         for (d = 0; d < DIM; d++)
 284         {
 285             clear_dvec(gyr_all[d]);
 286         }
 287
 288         if (bPC)
 289         {
 290             clear_dvec(sum_eig);
 291         }
 292
 293         if (outp)
 294         {
 295             for (i = 0; i < nat_min / 2; i++)
 296             {
 297                 sum_inp[i] = 0;
 298                 ninp[i]    = 0;
 299             }
 300         }
 301
 302         for (mol = 0; mol < nmol; mol++)
 303         {
 304             ind0 = molind[mol];
 305             ind1 = molind[mol + 1];
 306
 307             /* Determine end to end distance */
 308             sum_eed2 += distance2(x[index[ind0]], x[index[ind1 - 1]]);
 309
 310             /* Determine internal distances */
 311             if (outi)
 312             {
 313                 calc_int_dist(intd, x, index[ind0], index[ind1 - 1]);
 314             }
 315
 316             /* Determine the radius of gyration */
 317             clear_dvec(cm);
 318             for (d = 0; d < DIM; d++)
 319             {
 320                 clear_dvec(gyr[d]);
 321             }
 322             mmol = 0;
 323
 324             for (i = ind0; i < ind1; i++)
 325             {
 326                 a = index[i];
 327                 if (bMW)
 328                 {
 329                     m = top->atoms.atom[a].m;
 330                 }
 331                 else
 332                 {
 333                     m = 1;
 334                 }
 335                 mmol += m;
 336                 for (d = 0; d < DIM; d++)
 337                 {
 338                     cm[d] += m * x[a][d];
 339                     for (d2 = 0; d2 < DIM; d2++)
 340                     {
 341                         gyr[d][d2] += m * x[a][d] * x[a][d2];
 342                     }
 343                 }
 344             }
 345             dsvmul(1 / mmol, cm, cm);
 346             for (d = 0; d < DIM; d++)
 347             {
 348                 for (d2 = 0; d2 < DIM; d2++)
 349                 {
 350                     gyr[d][d2] = gyr[d][d2] / mmol - cm[d] * cm[d2];
 351                     gyr_all[d][d2] += gyr[d][d2];
 352                 }
 353             }
 354             if (bPC)
 355             {
 356                 gyro_eigen(gyr, eig, eigv, ord);
 357                 for (d = 0; d < DIM; d++)
 358                 {
 359                     sum_eig[d] += eig[ord[d]];
 360                 }
 361             }
 362             if (outp)
 363             {
 364                 for (i = ind0; i < ind1 - 1; i++)
 365                 {
 366                     rvec_sub(x[index[i + 1]], x[index[i]], bond[i - ind0]);
 367                     unitv(bond[i - ind0], bond[i - ind0]);
 368                 }
 369                 for (i = ind0; i < ind1 - 1; i++)
 370                 {
 371                     for (j = 0; (i + j < ind1 - 1 && j < nat_min / 2); j += 2)
 372                     {
 373                         sum_inp[j] += iprod(bond[i - ind0], bond[i - ind0 + j]);
 374                         ninp[j]++;
 375                     }
 376                 }
 377             }
 378         }
 379         sum_eed2 /= nmol;
 380
 381         sum_gyro = 0;
 382         for (d = 0; d < DIM; d++)
 383         {
 384             for (d2 = 0; d2 < DIM; d2++)
 385             {
 386                 gyr_all[d][d2] /= nmol;
 387             }
 388             sum_gyro += gyr_all[d][d];
 389         }
 390
 391         gyro_eigen(gyr_all, eig, eigv, ord);
 392
 393         fprintf(out, "%10.3f %8.4f %8.4f %8.4f %8.4f %8.4f", t * output_env_get_time_factor(oenv),
 394                 std::sqrt(sum_eed2), sqrt(sum_gyro), std::sqrt(eig[ord[0]]), std::sqrt(eig[ord[1]]),
 395                 std::sqrt(eig[ord[2]]));
 396         if (bPC)
 397         {
 398             for (d = 0; d < DIM; d++)
 399             {
 400                 fprintf(out, " %8.4f", std::sqrt(sum_eig[d] / nmol));
 401             }
 402         }
 403         fprintf(out, "\n");
 404
 405         if (outv)
 406         {
 407             fprintf(outv, "%10.3f", t * output_env_get_time_factor(oenv));
 408             for (d = 0; d < DIM; d++)
 409             {
 410                 for (d2 = 0; d2 < DIM; d2++)
 411                 {
 412                     fprintf(outv, " %6.3f", eigv[ord[d]][d2]);
 413                 }
 414             }
 415             fprintf(outv, "\n");
 416         }
 417
 418         sum_eed2_tot += sum_eed2;
 419         sum_gyro_tot += sum_gyro;
 420
 421         if (outp)
 422         {
 423             i = -1;
 424             for (j = 0; j < nat_min / 2; j += 2)
 425             {
 426                 sum_inp[j] /= ninp[j];
 427                 if (i == -1 && sum_inp[j] <= std::exp(-1.0))
 428                 {
 429                     i = j;
 430                 }
 431             }
 432             if (i == -1)
 433             {
 434                 pers = j;
 435             }
 436             else
 437             {
 438                 /* Do linear interpolation on a log scale */
 439                 pers = i - 2.0
 440                        + 2.0 * (std::log(sum_inp[i - 2]) + 1.0)
 441                                  / (std::log(sum_inp[i - 2]) - std::log(sum_inp[i]));
 442             }
 443             fprintf(outp, "%10.3f %8.4f\n", t * output_env_get_time_factor(oenv), pers);
 444             sum_pers_tot += pers;
 445         }
 446
 447         frame++;
 448     } while (read_next_x(oenv, status, &t, x, box));
 449
 450     gmx_rmpbc_done(gpbc);
 451
 452     close_trx(status);
 453
 454     xvgrclose(out);
 455     if (outv)
 456     {
 457         xvgrclose(outv);
 458     }
 459     if (outp)
 460     {
 461         xvgrclose(outp);
 462     }
 463
 464     sum_eed2_tot /= frame;
 465     sum_gyro_tot /= frame;
 466     sum_pers_tot /= frame;
 467     fprintf(stdout, "\nAverage end to end distance: %.3f (nm)\n", std::sqrt(sum_eed2_tot));
 468     fprintf(stdout, "\nAverage radius of gyration:  %.3f (nm)\n", std::sqrt(sum_gyro_tot));
 469     if (opt2bSet("-p", NFILE, fnm))
 470     {
 471         fprintf(stdout, "\nAverage persistence length:  %.2f bonds\n", sum_pers_tot);
 472     }
 473
 474     /* Handle printing of internal distances. */
 475     if (outi)
 476     {
 477         if (output_env_get_print_xvgr_codes(oenv))
 478         {
 479             fprintf(outi, "@    xaxes scale Logarithmic\n");
 480         }
 481         ymax = -1;
 482         ymin = 1e300;
 483         j    = index[molind[1] - 1] - index[molind[0]]; /* Polymer length -1. */
 484         for (i = 0; i < j; i++)
 485         {
 486             intd[i] /= (i + 1) * frame * nmol;
 487             if (intd[i] > ymax)
 488             {
 489                 ymax = intd[i];
 490             }
 491             if (intd[i] < ymin)
 492             {
 493                 ymin = intd[i];
 494             }
 495         }
 496         xvgr_world(outi, 1, ymin, j, ymax, oenv);
 497         for (i = 0; i < j; i++)
 498         {
 499             fprintf(outi, "%d  %8.4f\n", i + 1, intd[i]);
 500         }
 501         xvgrclose(outi);
 502     }
 503
 504     do_view(oenv, opt2fn("-o", NFILE, fnm), "-nxy");
 505     if (opt2bSet("-v", NFILE, fnm))
 506     {
 507         do_view(oenv, opt2fn("-v", NFILE, fnm), "-nxy");
 508     }
 509     if (opt2bSet("-p", NFILE, fnm))
 510     {
 511         do_view(oenv, opt2fn("-p", NFILE, fnm), "-nxy");
 512     }
 513
 514     return 0;
 515 }