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 "futil.h"
  46 #include "statutil.h"
  47 #include "copyrite.h"
  48 #include "vec.h"
  49 #include "index.h"
  50 #include "macros.h"
  51 #include "gmx_fatal.h"
  52 #include "xvgr.h"
  53 #include "rmpbc.h"
  54 #include "tpxio.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         "[TT]g_polystat[tt] 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     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     snew(top, 1);
 185     ePBC = read_tpx_top(ftp2fn(efTPX, NFILE, fnm),
 186                         NULL, box, &natoms, NULL, NULL, NULL, top);
 187
 188     fprintf(stderr, "Select a group of polymer mainchain atoms:\n");
 189     get_index(&top->atoms, ftp2fn_null(efNDX, NFILE, fnm),
 190               1, &isize, &index, &grpname);
 191
 192     snew(molind, top->mols.nr+1);
 193     nmol = 0;
 194     mol  = -1;
 195     for (i = 0; i < isize; i++)
 196     {
 197         if (i == 0 || index[i] >= top->mols.index[mol+1])
 198         {
 199             molind[nmol++] = i;
 200             do
 201             {
 202                 mol++;
 203             }
 204             while (index[i] >= top->mols.index[mol+1]);
 205         }
 206     }
 207     molind[nmol] = i;
 208     nat_min      = top->atoms.nr;
 209     nat_max      = 0;
 210     for (mol = 0; mol < nmol; mol++)
 211     {
 212         nat_min = min(nat_min, molind[mol+1]-molind[mol]);
 213         nat_max = max(nat_max, molind[mol+1]-molind[mol]);
 214     }
 215     fprintf(stderr, "Group %s consists of %d molecules\n", grpname, nmol);
 216     fprintf(stderr, "Group size per molecule, min: %d atoms, max %d atoms\n",
 217             nat_min, nat_max);
 218
 219     sprintf(title, "Size of %d polymers", nmol);
 220     out = xvgropen(opt2fn("-o", NFILE, fnm), title, output_env_get_xvgr_tlabel(oenv), "(nm)",
 221                    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), "Principal components",
 227                         output_env_get_xvgr_tlabel(oenv), "(nm)", oenv);
 228         snew(legp, DIM*DIM);
 229         for (d = 0; d < DIM; d++)
 230         {
 231             for (d2 = 0; d2 < DIM; d2++)
 232             {
 233                 sprintf(buf, "eig%d %c", d+1, 'x'+d2);
 234                 legp[d*DIM+d2] = strdup(buf);
 235             }
 236         }
 237         xvgr_legend(outv, DIM*DIM, (const char**)legp, oenv);
 238     }
 239     else
 240     {
 241         outv = NULL;
 242     }
 243
 244     if (opt2bSet("-p", NFILE, fnm))
 245     {
 246         outp = xvgropen(opt2fn("-p", NFILE, fnm), "Persistence length",
 247                         output_env_get_xvgr_tlabel(oenv), "bonds", oenv);
 248         snew(bond, nat_max-1);
 249         snew(sum_inp, nat_min/2);
 250         snew(ninp, nat_min/2);
 251     }
 252     else
 253     {
 254         outp = NULL;
 255     }
 256
 257     if (opt2bSet("-i", NFILE, fnm))
 258     {
 259         outi = xvgropen(opt2fn("-i", NFILE, fnm), "Internal distances",
 260                         "n", "<R\\S2\\N(n)>/n (nm\\S2\\N)", oenv);
 261         i = index[molind[1]-1] - index[molind[0]]; /* Length of polymer -1 */
 262         snew(intd, i);
 263     }
 264     else
 265     {
 266         intd = NULL;
 267         outi = NULL;
 268     }
 269
 270     natoms = read_first_x(oenv, &status, ftp2fn(efTRX, NFILE, fnm), &t, &x, box);
 271
 272     snew(gyr, DIM);
 273     snew(gyr_all, DIM);
 274     snew(eigv, DIM);
 275     for (d = 0; d < DIM; d++)
 276     {
 277         snew(gyr[d], DIM);
 278         snew(gyr_all[d], DIM);
 279         snew(eigv[d], DIM);
 280     }
 281
 282     frame        = 0;
 283     sum_eed2_tot = 0;
 284     sum_gyro_tot = 0;
 285     sum_pers_tot = 0;
 286
 287     gpbc = gmx_rmpbc_init(&top->idef, ePBC, natoms, box);
 288
 289     do
 290     {
 291         gmx_rmpbc(gpbc, natoms, box, x);
 292
 293         sum_eed2 = 0;
 294         for (d = 0; d < DIM; d++)
 295         {
 296             clear_dvec(gyr_all[d]);
 297         }
 298
 299         if (bPC)
 300         {
 301             clear_dvec(sum_eig);
 302         }
 303
 304         if (outp)
 305         {
 306             for (i = 0; i < nat_min/2; i++)
 307             {
 308                 sum_inp[i] = 0;
 309                 ninp[i]    = 0;
 310             }
 311         }
 312
 313         for (mol = 0; mol < nmol; mol++)
 314         {
 315             ind0 = molind[mol];
 316             ind1 = molind[mol+1];
 317
 318             /* Determine end to end distance */
 319             sum_eed2 += distance2(x[index[ind0]], x[index[ind1-1]]);
 320
 321             /* Determine internal distances */
 322             if (outi)
 323             {
 324                 calc_int_dist(intd, x, index[ind0], index[ind1-1]);
 325             }
 326
 327             /* Determine the radius of gyration */
 328             clear_dvec(cm);
 329             for (d = 0; d < DIM; d++)
 330             {
 331                 clear_dvec(gyr[d]);
 332             }
 333             mmol = 0;
 334
 335             for (i = ind0; i < ind1; i++)
 336             {
 337                 a = index[i];
 338                 if (bMW)
 339                 {
 340                     m = top->atoms.atom[a].m;
 341                 }
 342                 else
 343                 {
 344                     m = 1;
 345                 }
 346                 mmol += m;
 347                 for (d = 0; d < DIM; d++)
 348                 {
 349                     cm[d] += m*x[a][d];
 350                     for (d2 = 0; d2 < DIM; d2++)
 351                     {
 352                         gyr[d][d2] += m*x[a][d]*x[a][d2];
 353                     }
 354                 }
 355             }
 356             dsvmul(1/mmol, cm, cm);
 357             for (d = 0; d < DIM; d++)
 358             {
 359                 for (d2 = 0; d2 < DIM; d2++)
 360                 {
 361                     gyr[d][d2]      = gyr[d][d2]/mmol - cm[d]*cm[d2];
 362                     gyr_all[d][d2] += gyr[d][d2];
 363                 }
 364             }
 365             if (bPC)
 366             {
 367                 gyro_eigen(gyr, eig, eigv, ord);
 368                 for (d = 0; d < DIM; d++)
 369                 {
 370                     sum_eig[d] += eig[ord[d]];
 371                 }
 372             }
 373             if (outp)
 374             {
 375                 for (i = ind0; i < ind1-1; i++)
 376                 {
 377                     rvec_sub(x[index[i+1]], x[index[i]], bond[i-ind0]);
 378                     unitv(bond[i-ind0], bond[i-ind0]);
 379                 }
 380                 for (i = ind0; i < ind1-1; i++)
 381                 {
 382                     for (j = 0; (i+j < ind1-1 && j < nat_min/2); j += 2)
 383                     {
 384                         sum_inp[j] += iprod(bond[i-ind0], bond[i-ind0+j]);
 385                         ninp[j]++;
 386                     }
 387                 }
 388             }
 389         }
 390         sum_eed2 /= nmol;
 391
 392         sum_gyro = 0;
 393         for (d = 0; d < DIM; d++)
 394         {
 395             for (d2 = 0; d2 < DIM; d2++)
 396             {
 397                 gyr_all[d][d2] /= nmol;
 398             }
 399             sum_gyro += gyr_all[d][d];
 400         }
 401
 402         gyro_eigen(gyr_all, eig, eigv, ord);
 403
 404         fprintf(out, "%10.3f %8.4f %8.4f %8.4f %8.4f %8.4f",
 405                 t*output_env_get_time_factor(oenv),
 406                 sqrt(sum_eed2), sqrt(sum_gyro),
 407                 sqrt(eig[ord[0]]), sqrt(eig[ord[1]]), sqrt(eig[ord[2]]));
 408         if (bPC)
 409         {
 410             for (d = 0; d < DIM; d++)
 411             {
 412                 fprintf(out, " %8.4f", sqrt(sum_eig[d]/nmol));
 413             }
 414         }
 415         fprintf(out, "\n");
 416
 417         if (outv)
 418         {
 419             fprintf(outv, "%10.3f", t*output_env_get_time_factor(oenv));
 420             for (d = 0; d < DIM; d++)
 421             {
 422                 for (d2 = 0; d2 < DIM; d2++)
 423                 {
 424                     fprintf(outv, " %6.3f", eigv[ord[d]][d2]);
 425                 }
 426             }
 427             fprintf(outv, "\n");
 428         }
 429
 430         sum_eed2_tot += sum_eed2;
 431         sum_gyro_tot += sum_gyro;
 432
 433         if (outp)
 434         {
 435             i = -1;
 436             for (j = 0; j < nat_min/2; j += 2)
 437             {
 438                 sum_inp[j] /= ninp[j];
 439                 if (i == -1 && sum_inp[j] <= exp(-1.0))
 440                 {
 441                     i = j;
 442                 }
 443             }
 444             if (i == -1)
 445             {
 446                 pers = j;
 447             }
 448             else
 449             {
 450                 /* Do linear interpolation on a log scale */
 451                 pers = i - 2
 452                     + 2*(log(sum_inp[i-2]) + 1)/(log(sum_inp[i-2]) - log(sum_inp[i]));
 453             }
 454             fprintf(outp, "%10.3f %8.4f\n", t*output_env_get_time_factor(oenv), pers);
 455             sum_pers_tot += pers;
 456         }
 457
 458         frame++;
 459     }
 460     while (read_next_x(oenv, status, &t, natoms, x, box));
 461
 462     gmx_rmpbc_done(gpbc);
 463
 464     close_trx(status);
 465
 466     ffclose(out);
 467     if (outv)
 468     {
 469         ffclose(outv);
 470     }
 471     if (outp)
 472     {
 473         ffclose(outp);
 474     }
 475
 476     sum_eed2_tot /= frame;
 477     sum_gyro_tot /= frame;
 478     sum_pers_tot /= frame;
 479     fprintf(stdout, "\nAverage end to end distance: %.3f (nm)\n",
 480             sqrt(sum_eed2_tot));
 481     fprintf(stdout, "\nAverage radius of gyration:  %.3f (nm)\n",
 482             sqrt(sum_gyro_tot));
 483     if (opt2bSet("-p", NFILE, fnm))
 484     {
 485         fprintf(stdout, "\nAverage persistence length:  %.2f bonds\n",
 486                 sum_pers_tot);
 487     }
 488
 489     /* Handle printing of internal distances. */
 490     if (outi)
 491     {
 492         fprintf(outi, "@    xaxes scale Logarithmic\n");
 493         ymax = -1;
 494         ymin = 1e300;
 495         j    = index[molind[1]-1] - index[molind[0]]; /* Polymer length -1. */
 496         for (i = 0; i < j; i++)
 497         {
 498             intd[i] /= (i + 1) * frame * nmol;
 499             if (intd[i] > ymax)
 500             {
 501                 ymax = intd[i];
 502             }
 503             if (intd[i] < ymin)
 504             {
 505                 ymin = intd[i];
 506             }
 507         }
 508         xvgr_world(outi, 1, ymin, j, ymax, oenv);
 509         for (i = 0; i < j; i++)
 510         {
 511             fprintf(outi, "%d  %8.4f\n", i+1, intd[i]);
 512         }
 513         ffclose(outi);
 514     }
 515
 516     do_view(oenv, opt2fn("-o", NFILE, fnm), "-nxy");
 517     if (opt2bSet("-v", NFILE, fnm))
 518     {
 519         do_view(oenv, opt2fn("-v", NFILE, fnm), "-nxy");
 520     }
 521     if (opt2bSet("-p", NFILE, fnm))
 522     {
 523         do_view(oenv, opt2fn("-p", NFILE, fnm), "-nxy");
 524     }
 525
 526     thanx(stderr);
 527
 528     return 0;
 529 }