src/gromacs/gmxana/gmx_polystat.c

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