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 gmx_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   t_trxstatus *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   const 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   gmx_rmpbc_t  gpbc=NULL;
 166
 167   CopyRight(stderr,argv[0]);
 168   parse_common_args(&argc,argv,
 169                     PCA_CAN_VIEW | PCA_CAN_TIME | PCA_TIME_UNIT | PCA_BE_NICE,
 170                     NFILE,fnm,asize(pa),pa,asize(desc),desc,0,NULL,&oenv);
 171
 172   snew(top,1);
 173   ePBC = read_tpx_top(ftp2fn(efTPX,NFILE,fnm),
 174                       NULL,box,&natoms,NULL,NULL,NULL,top);
 175
 176   fprintf(stderr,"Select a group of polymer mainchain atoms:\n");
 177   get_index(&top->atoms,ftp2fn_null(efNDX,NFILE,fnm),
 178             1,&isize,&index,&grpname);
 179
 180   snew(molind,top->mols.nr+1);
 181   nmol = 0;
 182   mol = -1;
 183   for(i=0; i<isize; i++) {
 184     if (i == 0 || index[i] >= top->mols.index[mol+1]) {
 185       molind[nmol++] = i;
 186       do {
 187         mol++;
 188       } while (index[i] >= top->mols.index[mol+1]);
 189     }
 190   }
 191   molind[nmol] = i;
 192   nat_min = top->atoms.nr;
 193   nat_max = 0;
 194   for(mol=0; mol<nmol; mol++) {
 195     nat_min = min(nat_min,molind[mol+1]-molind[mol]);
 196     nat_max = max(nat_max,molind[mol+1]-molind[mol]);
 197   }
 198   fprintf(stderr,"Group %s consists of %d molecules\n",grpname,nmol);
 199   fprintf(stderr,"Group size per molecule, min: %d atoms, max %d atoms\n",
 200           nat_min,nat_max);
 201
 202   sprintf(title,"Size of %d polymers",nmol);
 203   out = xvgropen(opt2fn("-o",NFILE,fnm),title,output_env_get_xvgr_tlabel(oenv),"(nm)",
 204                 oenv);
 205   xvgr_legend(out,bPC ? 8 : 5,leg,oenv);
 206
 207   if (opt2bSet("-v",NFILE,fnm)) {
 208     outv = xvgropen(opt2fn("-v",NFILE,fnm),"Principal components",
 209                     output_env_get_xvgr_tlabel(oenv),"(nm)",oenv);
 210     snew(legp,DIM*DIM);
 211     for(d=0; d<DIM; d++) {
 212       for(d2=0; d2<DIM; d2++) {
 213         sprintf(buf,"eig%d %c",d+1,'x'+d2);
 214         legp[d*DIM+d2] = strdup(buf);
 215       }
 216     }
 217     xvgr_legend(outv,DIM*DIM,(const char**)legp,oenv);
 218   } else {
 219     outv = NULL;
 220   }
 221
 222   if (opt2bSet("-p",NFILE,fnm)) {
 223     outp = xvgropen(opt2fn("-p",NFILE,fnm),"Persistence length",
 224                     output_env_get_xvgr_tlabel(oenv),"bonds",oenv);
 225     snew(bond,nat_max-1);
 226     snew(sum_inp,nat_min/2);
 227     snew(ninp,nat_min/2);
 228   } else {
 229     outp = NULL;
 230   }
 231
 232   if (opt2bSet("-i", NFILE, fnm)) {
 233     outi = xvgropen(opt2fn("-i", NFILE, fnm), "Internal distances",
 234                     "n", "<R\\S2\\N(n)>/n (nm\\S2\\N)",oenv);
 235     i = index[molind[1]-1] - index[molind[0]]; /* Length of polymer -1 */
 236     snew(intd, i);
 237   } else {
 238     intd = NULL;
 239     outi = NULL;
 240   }
 241
 242   natoms = read_first_x(oenv,&status,ftp2fn(efTRX,NFILE,fnm),&t,&x,box);
 243
 244   snew(gyr,DIM);
 245   snew(gyr_all,DIM);
 246   snew(eigv,DIM);
 247   for(d=0; d<DIM; d++) {
 248     snew(gyr[d],DIM);
 249     snew(gyr_all[d],DIM);
 250     snew(eigv[d],DIM);
 251   }
 252
 253   frame        = 0;
 254   sum_eed2_tot = 0;
 255   sum_gyro_tot = 0;
 256   sum_pers_tot = 0;
 257
 258   gpbc = gmx_rmpbc_init(&top->idef,ePBC,natoms,box);
 259
 260   do {
 261     gmx_rmpbc(gpbc,natoms,box,x);
 262
 263     sum_eed2 = 0;
 264     for(d=0; d<DIM; d++)
 265       clear_dvec(gyr_all[d]);
 266
 267     if (bPC)
 268       clear_dvec(sum_eig);
 269
 270     if (outp) {
 271       for(i=0; i<nat_min/2; i++) {
 272         sum_inp[i] = 0;
 273         ninp[i] = 0;
 274       }
 275     }
 276
 277     for(mol=0; mol<nmol; mol++) {
 278       ind0 = molind[mol];
 279       ind1 = molind[mol+1];
 280
 281       /* Determine end to end distance */
 282       sum_eed2 += distance2(x[index[ind0]],x[index[ind1-1]]);
 283
 284       /* Determine internal distances */
 285       if (outi) {
 286         calc_int_dist(intd, x, index[ind0], index[ind1-1]);
 287       }
 288
 289       /* Determine the radius of gyration */
 290       clear_dvec(cm);
 291       for(d=0; d<DIM; d++)
 292         clear_dvec(gyr[d]);
 293       mmol = 0;
 294
 295       for(i=ind0; i<ind1; i++) {
 296         a = index[i];
 297         if (bMW)
 298           m = top->atoms.atom[a].m;
 299         else
 300           m = 1;
 301         mmol += m;
 302         for(d=0; d<DIM; d++) {
 303           cm[d] += m*x[a][d];
 304           for(d2=0; d2<DIM; d2++)
 305             gyr[d][d2] += m*x[a][d]*x[a][d2];
 306         }
 307       }
 308       dsvmul(1/mmol,cm,cm);
 309       for(d=0; d<DIM; d++) {
 310         for(d2=0; d2<DIM; d2++) {
 311           gyr[d][d2] = gyr[d][d2]/mmol - cm[d]*cm[d2];
 312           gyr_all[d][d2] += gyr[d][d2];
 313         }
 314       }
 315       if (bPC) {
 316         gyro_eigen(gyr,eig,eigv,ord);
 317         for(d=0; d<DIM; d++)
 318           sum_eig[d] += eig[ord[d]];
 319       }
 320       if (outp) {
 321         for(i=ind0; i<ind1-1; i++) {
 322           rvec_sub(x[index[i+1]],x[index[i]],bond[i-ind0]);
 323           unitv(bond[i-ind0],bond[i-ind0]);
 324         }
 325         for(i=ind0; i<ind1-1; i++) {
 326           for(j=0; (i+j<ind1-1 && j<nat_min/2); j+=2) {
 327             sum_inp[j] += iprod(bond[i-ind0],bond[i-ind0+j]);
 328             ninp[j]++;
 329           }
 330         }
 331       }
 332     }
 333     sum_eed2 /= nmol;
 334
 335     sum_gyro = 0;
 336     for(d=0; d<DIM; d++) {
 337       for(d2=0; d2<DIM; d2++)
 338         gyr_all[d][d2] /= nmol;
 339       sum_gyro += gyr_all[d][d];
 340     }
 341
 342     gyro_eigen(gyr_all,eig,eigv,ord);
 343
 344     fprintf(out,"%10.3f %8.4f %8.4f %8.4f %8.4f %8.4f",
 345             t*output_env_get_time_factor(oenv),
 346             sqrt(sum_eed2),sqrt(sum_gyro),
 347             sqrt(eig[ord[0]]),sqrt(eig[ord[1]]),sqrt(eig[ord[2]]));
 348     if (bPC) {
 349       for(d=0; d<DIM; d++)
 350         fprintf(out," %8.4f",sqrt(sum_eig[d]/nmol));
 351     }
 352     fprintf(out,"\n");
 353
 354     if (outv) {
 355       fprintf(outv,"%10.3f",t*output_env_get_time_factor(oenv));
 356       for(d=0; d<DIM; d++) {
 357         for(d2=0; d2<DIM; d2++)
 358           fprintf(outv," %6.3f",eigv[ord[d]][d2]);
 359       }
 360       fprintf(outv,"\n");
 361     }
 362
 363     sum_eed2_tot += sum_eed2;
 364     sum_gyro_tot += sum_gyro;
 365
 366     if (outp) {
 367       i = -1;
 368       for(j=0; j<nat_min/2; j+=2) {
 369         sum_inp[j] /= ninp[j];
 370         if (i == -1 && sum_inp[j] <= exp(-1.0))
 371           i = j;
 372       }
 373       if (i == -1) {
 374         pers = j;
 375       } else {
 376         /* Do linear interpolation on a log scale */
 377         pers = i - 2
 378           + 2*(log(sum_inp[i-2]) + 1)/(log(sum_inp[i-2]) - log(sum_inp[i]));
 379       }
 380       fprintf(outp,"%10.3f %8.4f\n",t*output_env_get_time_factor(oenv),pers);
 381       sum_pers_tot += pers;
 382     }
 383
 384     frame++;
 385   } while (read_next_x(oenv,status,&t,natoms,x,box));
 386
 387   gmx_rmpbc_done(gpbc);
 388
 389   close_trx(status);
 390
 391   ffclose(out);
 392   if (outv)
 393     ffclose(outv);
 394   if (outp)
 395     ffclose(outp);
 396
 397   sum_eed2_tot /= frame;
 398   sum_gyro_tot /= frame;
 399   sum_pers_tot /= frame;
 400   fprintf(stdout,"\nAverage end to end distance: %.3f (nm)\n",
 401           sqrt(sum_eed2_tot));
 402   fprintf(stdout,"\nAverage radius of gyration:  %.3f (nm)\n",
 403           sqrt(sum_gyro_tot));
 404   if (opt2bSet("-p",NFILE,fnm))
 405     fprintf(stdout,"\nAverage persistence length:  %.2f bonds\n",
 406             sum_pers_tot);
 407
 408   /* Handle printing of internal distances. */
 409   if (outi) {
 410     fprintf(outi, "@    xaxes scale Logarithmic\n");
 411     ymax = -1;
 412     ymin = 1e300;
 413     j = index[molind[1]-1] - index[molind[0]]; /* Polymer length -1. */
 414     for (i = 0; i < j; i++) {
 415       intd[i] /= (i + 1) * frame * nmol;
 416       if (intd[i] > ymax)
 417         ymax = intd[i];
 418       if (intd[i] < ymin)
 419         ymin = intd[i];
 420     }
 421     xvgr_world(outi, 1, ymin, j, ymax,oenv);
 422     for (i = 0; i < j; i++) {
 423       fprintf(outi, "%d  %8.4f\n", i+1, intd[i]);
 424     }
 425     ffclose(outi);
 426   }
 427
 428   do_view(oenv,opt2fn("-o",NFILE,fnm),"-nxy");
 429   if (opt2bSet("-v",NFILE,fnm))
 430     do_view(oenv,opt2fn("-v",NFILE,fnm),"-nxy");
 431   if (opt2bSet("-p",NFILE,fnm))
 432     do_view(oenv,opt2fn("-p",NFILE,fnm),"-nxy");
 433
 434   thanx(stderr);
 435
 436   return 0;
 437 }