src/gromacs/gmxana/gmx_polystat.c

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