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