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