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