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, 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},
 138           "Use the mass weighting for radii of gyration" },
 139         { "-pc", FALSE, etBOOL, {&bPC},
 140           "Plot average eigenvalues" }
 141     };
 142
 143     t_filenm        fnm[] = {
 144         { efTPR, NULL, NULL,  ffREAD  },
 145         { efTRX, "-f", NULL,  ffREAD  },
 146         { efNDX, NULL, NULL,  ffOPTRD },
 147         { efXVG, "-o", "polystat",  ffWRITE },
 148         { efXVG, "-v", "polyvec", ffOPTWR },
 149         { efXVG, "-p", "persist",  ffOPTWR },
 150         { efXVG, "-i", "intdist", ffOPTWR }
 151     };
 152 #define NFILE asize(fnm)
 153
 154     t_topology       *top;
 155     gmx_output_env_t *oenv;
 156     int               ePBC;
 157     int               isize, *index, nmol, *molind, mol, nat_min = 0, nat_max = 0;
 158     char             *grpname;
 159     t_trxstatus      *status;
 160     real              t;
 161     rvec             *x, *bond = NULL;
 162     matrix            box;
 163     int               natoms, i, j, frame, ind0, ind1, a, d, d2, ord[DIM] = {0};
 164     dvec              cm, sum_eig = {0, 0, 0};
 165     double          **gyr, **gyr_all, eig[DIM], **eigv;
 166     double            sum_eed2, sum_eed2_tot, sum_gyro, sum_gyro_tot, sum_pers_tot;
 167     int              *ninp    = NULL;
 168     double           *sum_inp = NULL, pers;
 169     double           *intd, ymax, ymin;
 170     double            mmol, m;
 171     char              title[STRLEN];
 172     FILE             *out, *outv, *outp, *outi;
 173     const char       *leg[8] = {
 174         "end to end", "<R\\sg\\N>",
 175         "<R\\sg\\N> eig1", "<R\\sg\\N> eig2", "<R\\sg\\N> eig3",
 176         "<R\\sg\\N eig1>", "<R\\sg\\N eig2>", "<R\\sg\\N eig3>"
 177     };
 178     char            **legp, buf[STRLEN];
 179     gmx_rmpbc_t       gpbc = NULL;
 180
 181     if (!parse_common_args(&argc, argv,
 182                            PCA_CAN_VIEW | PCA_CAN_TIME | PCA_TIME_UNIT,
 183                            NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, NULL, &oenv))
 184     {
 185         return 0;
 186     }
 187
 188     snew(top, 1);
 189     ePBC = read_tpx_top(ftp2fn(efTPR, NFILE, fnm),
 190                         NULL, box, &natoms, NULL, NULL, top);
 191
 192     fprintf(stderr, "Select a group of polymer mainchain atoms:\n");
 193     get_index(&top->atoms, ftp2fn_null(efNDX, NFILE, fnm),
 194               1, &isize, &index, &grpname);
 195
 196     snew(molind, top->mols.nr+1);
 197     nmol = 0;
 198     mol  = -1;
 199     for (i = 0; i < isize; i++)
 200     {
 201         if (i == 0 || index[i] >= top->mols.index[mol+1])
 202         {
 203             molind[nmol++] = i;
 204             do
 205             {
 206                 mol++;
 207             }
 208             while (index[i] >= top->mols.index[mol+1]);
 209         }
 210     }
 211     molind[nmol] = i;
 212     nat_min      = top->atoms.nr;
 213     nat_max      = 0;
 214     for (mol = 0; mol < nmol; mol++)
 215     {
 216         nat_min = std::min(nat_min, molind[mol+1]-molind[mol]);
 217         nat_max = std::max(nat_max, molind[mol+1]-molind[mol]);
 218     }
 219     fprintf(stderr, "Group %s consists of %d molecules\n", grpname, nmol);
 220     fprintf(stderr, "Group size per molecule, min: %d atoms, max %d atoms\n",
 221             nat_min, nat_max);
 222
 223     sprintf(title, "Size of %d polymers", nmol);
 224     out = xvgropen(opt2fn("-o", NFILE, fnm), title, output_env_get_xvgr_tlabel(oenv), "(nm)",
 225                    oenv);
 226     xvgr_legend(out, bPC ? 8 : 5, leg, oenv);
 227
 228     if (opt2bSet("-v", NFILE, fnm))
 229     {
 230         outv = xvgropen(opt2fn("-v", NFILE, fnm), "Principal components",
 231                         output_env_get_xvgr_tlabel(oenv), "(nm)", oenv);
 232         snew(legp, DIM*DIM);
 233         for (d = 0; d < DIM; d++)
 234         {
 235             for (d2 = 0; d2 < DIM; d2++)
 236             {
 237                 sprintf(buf, "eig%d %c", d+1, 'x'+d2);
 238                 legp[d*DIM+d2] = gmx_strdup(buf);
 239             }
 240         }
 241         xvgr_legend(outv, DIM*DIM, (const char**)legp, oenv);
 242     }
 243     else
 244     {
 245         outv = NULL;
 246     }
 247
 248     if (opt2bSet("-p", NFILE, fnm))
 249     {
 250         outp = xvgropen(opt2fn("-p", NFILE, fnm), "Persistence length",
 251                         output_env_get_xvgr_tlabel(oenv), "bonds", oenv);
 252         snew(bond, nat_max-1);
 253         snew(sum_inp, nat_min/2);
 254         snew(ninp, nat_min/2);
 255     }
 256     else
 257     {
 258         outp = NULL;
 259     }
 260
 261     if (opt2bSet("-i", NFILE, fnm))
 262     {
 263         outi = xvgropen(opt2fn("-i", NFILE, fnm), "Internal distances",
 264                         "n", "<R\\S2\\N(n)>/n (nm\\S2\\N)", oenv);
 265         i = index[molind[1]-1] - index[molind[0]]; /* Length of polymer -1 */
 266         snew(intd, i);
 267     }
 268     else
 269     {
 270         intd = NULL;
 271         outi = NULL;
 272     }
 273
 274     natoms = read_first_x(oenv, &status, ftp2fn(efTRX, NFILE, fnm), &t, &x, box);
 275
 276     snew(gyr, DIM);
 277     snew(gyr_all, DIM);
 278     snew(eigv, DIM);
 279     for (d = 0; d < DIM; d++)
 280     {
 281         snew(gyr[d], DIM);
 282         snew(gyr_all[d], DIM);
 283         snew(eigv[d], DIM);
 284     }
 285
 286     frame        = 0;
 287     sum_eed2_tot = 0;
 288     sum_gyro_tot = 0;
 289     sum_pers_tot = 0;
 290
 291     gpbc = gmx_rmpbc_init(&top->idef, ePBC, natoms);
 292
 293     do
 294     {
 295         gmx_rmpbc(gpbc, natoms, box, x);
 296
 297         sum_eed2 = 0;
 298         for (d = 0; d < DIM; d++)
 299         {
 300             clear_dvec(gyr_all[d]);
 301         }
 302
 303         if (bPC)
 304         {
 305             clear_dvec(sum_eig);
 306         }
 307
 308         if (outp)
 309         {
 310             for (i = 0; i < nat_min/2; i++)
 311             {
 312                 sum_inp[i] = 0;
 313                 ninp[i]    = 0;
 314             }
 315         }
 316
 317         for (mol = 0; mol < nmol; mol++)
 318         {
 319             ind0 = molind[mol];
 320             ind1 = molind[mol+1];
 321
 322             /* Determine end to end distance */
 323             sum_eed2 += distance2(x[index[ind0]], x[index[ind1-1]]);
 324
 325             /* Determine internal distances */
 326             if (outi)
 327             {
 328                 calc_int_dist(intd, x, index[ind0], index[ind1-1]);
 329             }
 330
 331             /* Determine the radius of gyration */
 332             clear_dvec(cm);
 333             for (d = 0; d < DIM; d++)
 334             {
 335                 clear_dvec(gyr[d]);
 336             }
 337             mmol = 0;
 338
 339             for (i = ind0; i < ind1; i++)
 340             {
 341                 a = index[i];
 342                 if (bMW)
 343                 {
 344                     m = top->atoms.atom[a].m;
 345                 }
 346                 else
 347                 {
 348                     m = 1;
 349                 }
 350                 mmol += m;
 351                 for (d = 0; d < DIM; d++)
 352                 {
 353                     cm[d] += m*x[a][d];
 354                     for (d2 = 0; d2 < DIM; d2++)
 355                     {
 356                         gyr[d][d2] += m*x[a][d]*x[a][d2];
 357                     }
 358                 }
 359             }
 360             dsvmul(1/mmol, cm, cm);
 361             for (d = 0; d < DIM; d++)
 362             {
 363                 for (d2 = 0; d2 < DIM; d2++)
 364                 {
 365                     gyr[d][d2]      = gyr[d][d2]/mmol - cm[d]*cm[d2];
 366                     gyr_all[d][d2] += gyr[d][d2];
 367                 }
 368             }
 369             if (bPC)
 370             {
 371                 gyro_eigen(gyr, eig, eigv, ord);
 372                 for (d = 0; d < DIM; d++)
 373                 {
 374                     sum_eig[d] += eig[ord[d]];
 375                 }
 376             }
 377             if (outp)
 378             {
 379                 for (i = ind0; i < ind1-1; i++)
 380                 {
 381                     rvec_sub(x[index[i+1]], x[index[i]], bond[i-ind0]);
 382                     unitv(bond[i-ind0], bond[i-ind0]);
 383                 }
 384                 for (i = ind0; i < ind1-1; i++)
 385                 {
 386                     for (j = 0; (i+j < ind1-1 && j < nat_min/2); j += 2)
 387                     {
 388                         sum_inp[j] += iprod(bond[i-ind0], bond[i-ind0+j]);
 389                         ninp[j]++;
 390                     }
 391                 }
 392             }
 393         }
 394         sum_eed2 /= nmol;
 395
 396         sum_gyro = 0;
 397         for (d = 0; d < DIM; d++)
 398         {
 399             for (d2 = 0; d2 < DIM; d2++)
 400             {
 401                 gyr_all[d][d2] /= nmol;
 402             }
 403             sum_gyro += gyr_all[d][d];
 404         }
 405
 406         gyro_eigen(gyr_all, eig, eigv, ord);
 407
 408         fprintf(out, "%10.3f %8.4f %8.4f %8.4f %8.4f %8.4f",
 409                 t*output_env_get_time_factor(oenv),
 410                 std::sqrt(sum_eed2), sqrt(sum_gyro),
 411                 std::sqrt(eig[ord[0]]), std::sqrt(eig[ord[1]]), std::sqrt(eig[ord[2]]));
 412         if (bPC)
 413         {
 414             for (d = 0; d < DIM; d++)
 415             {
 416                 fprintf(out, " %8.4f", std::sqrt(sum_eig[d]/nmol));
 417             }
 418         }
 419         fprintf(out, "\n");
 420
 421         if (outv)
 422         {
 423             fprintf(outv, "%10.3f", t*output_env_get_time_factor(oenv));
 424             for (d = 0; d < DIM; d++)
 425             {
 426                 for (d2 = 0; d2 < DIM; d2++)
 427                 {
 428                     fprintf(outv, " %6.3f", eigv[ord[d]][d2]);
 429                 }
 430             }
 431             fprintf(outv, "\n");
 432         }
 433
 434         sum_eed2_tot += sum_eed2;
 435         sum_gyro_tot += sum_gyro;
 436
 437         if (outp)
 438         {
 439             i = -1;
 440             for (j = 0; j < nat_min/2; j += 2)
 441             {
 442                 sum_inp[j] /= ninp[j];
 443                 if (i == -1 && sum_inp[j] <= std::exp(-1.0))
 444                 {
 445                     i = j;
 446                 }
 447             }
 448             if (i == -1)
 449             {
 450                 pers = j;
 451             }
 452             else
 453             {
 454                 /* Do linear interpolation on a log scale */
 455                 pers = i - 2.0
 456                     + 2.0*(std::log(sum_inp[i-2]) + 1.0)/(std::log(sum_inp[i-2]) - std::log(sum_inp[i]));
 457             }
 458             fprintf(outp, "%10.3f %8.4f\n", t*output_env_get_time_factor(oenv), pers);
 459             sum_pers_tot += pers;
 460         }
 461
 462         frame++;
 463     }
 464     while (read_next_x(oenv, status, &t, x, box));
 465
 466     gmx_rmpbc_done(gpbc);
 467
 468     close_trx(status);
 469
 470     xvgrclose(out);
 471     if (outv)
 472     {
 473         xvgrclose(outv);
 474     }
 475     if (outp)
 476     {
 477         xvgrclose(outp);
 478     }
 479
 480     sum_eed2_tot /= frame;
 481     sum_gyro_tot /= frame;
 482     sum_pers_tot /= frame;
 483     fprintf(stdout, "\nAverage end to end distance: %.3f (nm)\n",
 484             std::sqrt(sum_eed2_tot));
 485     fprintf(stdout, "\nAverage radius of gyration:  %.3f (nm)\n",
 486             std::sqrt(sum_gyro_tot));
 487     if (opt2bSet("-p", NFILE, fnm))
 488     {
 489         fprintf(stdout, "\nAverage persistence length:  %.2f bonds\n",
 490                 sum_pers_tot);
 491     }
 492
 493     /* Handle printing of internal distances. */
 494     if (outi)
 495     {
 496         if (output_env_get_print_xvgr_codes(oenv))
 497         {
 498             fprintf(outi, "@    xaxes scale Logarithmic\n");
 499         }
 500         ymax = -1;
 501         ymin = 1e300;
 502         j    = index[molind[1]-1] - index[molind[0]]; /* Polymer length -1. */
 503         for (i = 0; i < j; i++)
 504         {
 505             intd[i] /= (i + 1) * frame * nmol;
 506             if (intd[i] > ymax)
 507             {
 508                 ymax = intd[i];
 509             }
 510             if (intd[i] < ymin)
 511             {
 512                 ymin = intd[i];
 513             }
 514         }
 515         xvgr_world(outi, 1, ymin, j, ymax, oenv);
 516         for (i = 0; i < j; i++)
 517         {
 518             fprintf(outi, "%d  %8.4f\n", i+1, intd[i]);
 519         }
 520         xvgrclose(outi);
 521     }
 522
 523     do_view(oenv, opt2fn("-o", NFILE, fnm), "-nxy");
 524     if (opt2bSet("-v", NFILE, fnm))
 525     {
 526         do_view(oenv, opt2fn("-v", NFILE, fnm), "-nxy");
 527     }
 528     if (opt2bSet("-p", NFILE, fnm))
 529     {
 530         do_view(oenv, opt2fn("-p", NFILE, fnm), "-nxy");
 531     }
 532
 533     return 0;
 534 }