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