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