src/tools/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  * check out http://www.gromacs.org for more information.
   7  * Copyright (c) 2012, by the GROMACS development team, led by
   8  * David van der Spoel, Berk Hess, Erik Lindahl, and including many
   9  * others, as listed in the AUTHORS file in the top-level source
  10  * directory and at http://www.gromacs.org.
  11  *
  12  * GROMACS is free software; you can redistribute it and/or
  13  * modify it under the terms of the GNU Lesser General Public License
  14  * as published by the Free Software Foundation; either version 2.1
  15  * of the License, or (at your option) any later version.
  16  *
  17  * GROMACS is distributed in the hope that it will be useful,
  18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  20  * Lesser General Public License for more details.
  21  *
  22  * You should have received a copy of the GNU Lesser General Public
  23  * License along with GROMACS; if not, see
  24  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  25  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  26  *
  27  * If you want to redistribute modifications to GROMACS, please
  28  * consider that scientific software is very special. Version
  29  * control is crucial - bugs must be traceable. We will be happy to
  30  * consider code for inclusion in the official distribution, but
  31  * derived work must not be called official GROMACS. Details are found
  32  * in the README & COPYING files - if they are missing, get the
  33  * official version at http://www.gromacs.org.
  34  *
  35  * To help us fund GROMACS development, we humbly ask that you cite
  36  * the research papers on the package. Check out http://www.gromacs.org.
  37  */
  38 #ifdef HAVE_CONFIG_H
  39 #include <config.h>
  40 #endif
  41 #include <math.h>
  42 #include <string.h>
  43
  44 #include "sysstuff.h"
  45 #include "physics.h"
  46 #include "typedefs.h"
  47 #include "smalloc.h"
  48 #include "futil.h"
  49 #include "statutil.h"
  50 #include "copyrite.h"
  51 #include "vec.h"
  52 #include "index.h"
  53 #include "macros.h"
  54 #include "gmx_fatal.h"
  55 #include "xvgr.h"
  56 #include "rmpbc.h"
  57 #include "tpxio.h"
  58 #include "nrjac.h"
  59 #include "gmx_ana.h"
  60
  61
  62 static void gyro_eigen(double **gyr,double *eig,double **eigv,int *ord)
  63 {
  64   int nrot,d;
  65
  66   jacobi(gyr,DIM,eig,eigv,&nrot);
  67   /* Order the eigenvalues */
  68   ord[0] = 0;
  69   ord[2] = 2;
  70   for(d=0; d<DIM; d++) {
  71     if (eig[d] > eig[ord[0]])
  72       ord[0] = d;
  73     if (eig[d] < eig[ord[2]])
  74       ord[2] = d;
  75   }
  76   for(d=0; d<DIM; d++) {
  77     if (ord[0] != d && ord[2] != d)
  78       ord[1] = d;
  79   }
  80 }
  81
  82 /* Calculate mean square internal distances (Auhl et al., JCP 119, 12718) */
  83 static void calc_int_dist(double *intd, rvec *x, int i0, int i1)
  84 {
  85   int ii, jj;
  86   const int ml = i1 - i0 + 1; /* Number of beads in molecule. */
  87   int bd; /* Distance between beads */
  88   double d;
  89
  90   for (bd = 1; bd < ml; bd++) {
  91     d = 0.;
  92     for (ii = i0; ii <= i1 - bd; ii++)
  93       d += distance2(x[ii], x[ii+bd]);
  94     d /= ml - bd;
  95     intd[bd - 1] += d;
  96   }
  97 }
  98
  99 int gmx_polystat(int argc,char *argv[])
 100 {
 101   const char *desc[] = {
 102     "[TT]g_polystat[tt] plots static properties of polymers as a function of time",
 103     "and prints the average.[PAR]",
 104     "By default it determines the average end-to-end distance and radii",
 105     "of gyration of polymers. It asks for an index group and split this",
 106     "into molecules. The end-to-end distance is then determined using",
 107     "the first and the last atom in the index group for each molecules.",
 108     "For the radius of gyration the total and the three principal components",
 109     "for the average gyration tensor are written.",
 110     "With option [TT]-v[tt] the eigenvectors are written.",
 111     "With option [TT]-pc[tt] also the average eigenvalues of the individual",
 112     "gyration tensors are written.",
 113     "With option [TT]-i[tt] the mean square internal distances are",
 114     "written.[PAR]",
 115     "With option [TT]-p[tt] the persistence length is determined.",
 116     "The chosen index group should consist of atoms that are",
 117     "consecutively bonded in the polymer mainchains.",
 118     "The persistence length is then determined from the cosine of",
 119     "the angles between bonds with an index difference that is even,",
 120     "the odd pairs are not used, because straight polymer backbones",
 121     "are usually all trans and therefore only every second bond aligns.",
 122     "The persistence length is defined as number of bonds where",
 123     "the average cos reaches a value of 1/e. This point is determined",
 124     "by a linear interpolation of log(<cos>)."
 125   };
 126   static gmx_bool bMW = TRUE, bPC = FALSE;
 127   t_pargs pa[] = {
 128     { "-mw", FALSE, etBOOL, {&bMW},
 129       "Use the mass weighting for radii of gyration" },
 130     { "-pc", FALSE, etBOOL, {&bPC},
 131       "Plot average eigenvalues" }
 132   };
 133
 134   t_filenm   fnm[] = {
 135     { efTPX, NULL, NULL,  ffREAD  },
 136     { efTRX, "-f", NULL,  ffREAD  },
 137     { efNDX, NULL, NULL,  ffOPTRD },
 138     { efXVG, "-o", "polystat",  ffWRITE },
 139     { efXVG, "-v", "polyvec", ffOPTWR },
 140     { efXVG, "-p", "persist",  ffOPTWR },
 141     { efXVG, "-i", "intdist", ffOPTWR }
 142   };
 143 #define NFILE asize(fnm)
 144
 145   t_topology *top;
 146   output_env_t oenv;
 147   int    ePBC;
 148   int    isize,*index,nmol,*molind,mol,nat_min=0,nat_max=0;
 149   char   *grpname;
 150   t_trxstatus *status;
 151   real   t;
 152   rvec   *x,*bond=NULL;
 153   matrix box;
 154   int    natoms,i,j,frame,ind0,ind1,a,d,d2,ord[DIM]={0};
 155   dvec   cm,sum_eig={0,0,0};
 156   double **gyr,**gyr_all,eig[DIM],**eigv;
 157   double sum_eed2,sum_eed2_tot,sum_gyro,sum_gyro_tot,sum_pers_tot;
 158   int    *ninp=NULL;
 159   double *sum_inp=NULL,pers;
 160   double *intd, ymax, ymin;
 161   double mmol,m;
 162   char   title[STRLEN];
 163   FILE   *out,*outv,*outp, *outi;
 164   const char *leg[8] = { "end to end", "<R\\sg\\N>",
 165                          "<R\\sg\\N> eig1", "<R\\sg\\N> eig2", "<R\\sg\\N> eig3",
 166                          "<R\\sg\\N eig1>", "<R\\sg\\N eig2>", "<R\\sg\\N eig3>" };
 167   char   **legp,buf[STRLEN];
 168   gmx_rmpbc_t  gpbc=NULL;
 169
 170   CopyRight(stderr,argv[0]);
 171   parse_common_args(&argc,argv,
 172                     PCA_CAN_VIEW | PCA_CAN_TIME | PCA_TIME_UNIT | PCA_BE_NICE,
 173                     NFILE,fnm,asize(pa),pa,asize(desc),desc,0,NULL,&oenv);
 174
 175   snew(top,1);
 176   ePBC = read_tpx_top(ftp2fn(efTPX,NFILE,fnm),
 177                       NULL,box,&natoms,NULL,NULL,NULL,top);
 178
 179   fprintf(stderr,"Select a group of polymer mainchain atoms:\n");
 180   get_index(&top->atoms,ftp2fn_null(efNDX,NFILE,fnm),
 181             1,&isize,&index,&grpname);
 182
 183   snew(molind,top->mols.nr+1);
 184   nmol = 0;
 185   mol = -1;
 186   for(i=0; i<isize; i++) {
 187     if (i == 0 || index[i] >= top->mols.index[mol+1]) {
 188       molind[nmol++] = i;
 189       do {
 190         mol++;
 191       } while (index[i] >= top->mols.index[mol+1]);
 192     }
 193   }
 194   molind[nmol] = i;
 195   nat_min = top->atoms.nr;
 196   nat_max = 0;
 197   for(mol=0; mol<nmol; mol++) {
 198     nat_min = min(nat_min,molind[mol+1]-molind[mol]);
 199     nat_max = max(nat_max,molind[mol+1]-molind[mol]);
 200   }
 201   fprintf(stderr,"Group %s consists of %d molecules\n",grpname,nmol);
 202   fprintf(stderr,"Group size per molecule, min: %d atoms, max %d atoms\n",
 203           nat_min,nat_max);
 204
 205   sprintf(title,"Size of %d polymers",nmol);
 206   out = xvgropen(opt2fn("-o",NFILE,fnm),title,output_env_get_xvgr_tlabel(oenv),"(nm)",
 207                 oenv);
 208   xvgr_legend(out,bPC ? 8 : 5,leg,oenv);
 209
 210   if (opt2bSet("-v",NFILE,fnm)) {
 211     outv = xvgropen(opt2fn("-v",NFILE,fnm),"Principal components",
 212                     output_env_get_xvgr_tlabel(oenv),"(nm)",oenv);
 213     snew(legp,DIM*DIM);
 214     for(d=0; d<DIM; d++) {
 215       for(d2=0; d2<DIM; d2++) {
 216         sprintf(buf,"eig%d %c",d+1,'x'+d2);
 217         legp[d*DIM+d2] = strdup(buf);
 218       }
 219     }
 220     xvgr_legend(outv,DIM*DIM,(const char**)legp,oenv);
 221   } else {
 222     outv = NULL;
 223   }
 224
 225   if (opt2bSet("-p",NFILE,fnm)) {
 226     outp = xvgropen(opt2fn("-p",NFILE,fnm),"Persistence length",
 227                     output_env_get_xvgr_tlabel(oenv),"bonds",oenv);
 228     snew(bond,nat_max-1);
 229     snew(sum_inp,nat_min/2);
 230     snew(ninp,nat_min/2);
 231   } else {
 232     outp = NULL;
 233   }
 234
 235   if (opt2bSet("-i", NFILE, fnm)) {
 236     outi = xvgropen(opt2fn("-i", NFILE, fnm), "Internal distances",
 237                     "n", "<R\\S2\\N(n)>/n (nm\\S2\\N)",oenv);
 238     i = index[molind[1]-1] - index[molind[0]]; /* Length of polymer -1 */
 239     snew(intd, i);
 240   } else {
 241     intd = NULL;
 242     outi = NULL;
 243   }
 244
 245   natoms = read_first_x(oenv,&status,ftp2fn(efTRX,NFILE,fnm),&t,&x,box);
 246
 247   snew(gyr,DIM);
 248   snew(gyr_all,DIM);
 249   snew(eigv,DIM);
 250   for(d=0; d<DIM; d++) {
 251     snew(gyr[d],DIM);
 252     snew(gyr_all[d],DIM);
 253     snew(eigv[d],DIM);
 254   }
 255
 256   frame        = 0;
 257   sum_eed2_tot = 0;
 258   sum_gyro_tot = 0;
 259   sum_pers_tot = 0;
 260
 261   gpbc = gmx_rmpbc_init(&top->idef,ePBC,natoms,box);
 262
 263   do {
 264     gmx_rmpbc(gpbc,natoms,box,x);
 265
 266     sum_eed2 = 0;
 267     for(d=0; d<DIM; d++)
 268       clear_dvec(gyr_all[d]);
 269
 270     if (bPC)
 271       clear_dvec(sum_eig);
 272
 273     if (outp) {
 274       for(i=0; i<nat_min/2; i++) {
 275         sum_inp[i] = 0;
 276         ninp[i] = 0;
 277       }
 278     }
 279
 280     for(mol=0; mol<nmol; mol++) {
 281       ind0 = molind[mol];
 282       ind1 = molind[mol+1];
 283
 284       /* Determine end to end distance */
 285       sum_eed2 += distance2(x[index[ind0]],x[index[ind1-1]]);
 286
 287       /* Determine internal distances */
 288       if (outi) {
 289         calc_int_dist(intd, x, index[ind0], index[ind1-1]);
 290       }
 291
 292       /* Determine the radius of gyration */
 293       clear_dvec(cm);
 294       for(d=0; d<DIM; d++)
 295         clear_dvec(gyr[d]);
 296       mmol = 0;
 297
 298       for(i=ind0; i<ind1; i++) {
 299         a = index[i];
 300         if (bMW)
 301           m = top->atoms.atom[a].m;
 302         else
 303           m = 1;
 304         mmol += m;
 305         for(d=0; d<DIM; d++) {
 306           cm[d] += m*x[a][d];
 307           for(d2=0; d2<DIM; d2++)
 308             gyr[d][d2] += m*x[a][d]*x[a][d2];
 309         }
 310       }
 311       dsvmul(1/mmol,cm,cm);
 312       for(d=0; d<DIM; d++) {
 313         for(d2=0; d2<DIM; d2++) {
 314           gyr[d][d2] = gyr[d][d2]/mmol - cm[d]*cm[d2];
 315           gyr_all[d][d2] += gyr[d][d2];
 316         }
 317       }
 318       if (bPC) {
 319         gyro_eigen(gyr,eig,eigv,ord);
 320         for(d=0; d<DIM; d++)
 321           sum_eig[d] += eig[ord[d]];
 322       }
 323       if (outp) {
 324         for(i=ind0; i<ind1-1; i++) {
 325           rvec_sub(x[index[i+1]],x[index[i]],bond[i-ind0]);
 326           unitv(bond[i-ind0],bond[i-ind0]);
 327         }
 328         for(i=ind0; i<ind1-1; i++) {
 329           for(j=0; (i+j<ind1-1 && j<nat_min/2); j+=2) {
 330             sum_inp[j] += iprod(bond[i-ind0],bond[i-ind0+j]);
 331             ninp[j]++;
 332           }
 333         }
 334       }
 335     }
 336     sum_eed2 /= nmol;
 337
 338     sum_gyro = 0;
 339     for(d=0; d<DIM; d++) {
 340       for(d2=0; d2<DIM; d2++)
 341         gyr_all[d][d2] /= nmol;
 342       sum_gyro += gyr_all[d][d];
 343     }
 344
 345     gyro_eigen(gyr_all,eig,eigv,ord);
 346
 347     fprintf(out,"%10.3f %8.4f %8.4f %8.4f %8.4f %8.4f",
 348             t*output_env_get_time_factor(oenv),
 349             sqrt(sum_eed2),sqrt(sum_gyro),
 350             sqrt(eig[ord[0]]),sqrt(eig[ord[1]]),sqrt(eig[ord[2]]));
 351     if (bPC) {
 352       for(d=0; d<DIM; d++)
 353         fprintf(out," %8.4f",sqrt(sum_eig[d]/nmol));
 354     }
 355     fprintf(out,"\n");
 356
 357     if (outv) {
 358       fprintf(outv,"%10.3f",t*output_env_get_time_factor(oenv));
 359       for(d=0; d<DIM; d++) {
 360         for(d2=0; d2<DIM; d2++)
 361           fprintf(outv," %6.3f",eigv[ord[d]][d2]);
 362       }
 363       fprintf(outv,"\n");
 364     }
 365
 366     sum_eed2_tot += sum_eed2;
 367     sum_gyro_tot += sum_gyro;
 368
 369     if (outp) {
 370       i = -1;
 371       for(j=0; j<nat_min/2; j+=2) {
 372         sum_inp[j] /= ninp[j];
 373         if (i == -1 && sum_inp[j] <= exp(-1.0))
 374           i = j;
 375       }
 376       if (i == -1) {
 377         pers = j;
 378       } else {
 379         /* Do linear interpolation on a log scale */
 380         pers = i - 2
 381           + 2*(log(sum_inp[i-2]) + 1)/(log(sum_inp[i-2]) - log(sum_inp[i]));
 382       }
 383       fprintf(outp,"%10.3f %8.4f\n",t*output_env_get_time_factor(oenv),pers);
 384       sum_pers_tot += pers;
 385     }
 386
 387     frame++;
 388   } while (read_next_x(oenv,status,&t,natoms,x,box));
 389
 390   gmx_rmpbc_done(gpbc);
 391
 392   close_trx(status);
 393
 394   ffclose(out);
 395   if (outv)
 396     ffclose(outv);
 397   if (outp)
 398     ffclose(outp);
 399
 400   sum_eed2_tot /= frame;
 401   sum_gyro_tot /= frame;
 402   sum_pers_tot /= frame;
 403   fprintf(stdout,"\nAverage end to end distance: %.3f (nm)\n",
 404           sqrt(sum_eed2_tot));
 405   fprintf(stdout,"\nAverage radius of gyration:  %.3f (nm)\n",
 406           sqrt(sum_gyro_tot));
 407   if (opt2bSet("-p",NFILE,fnm))
 408     fprintf(stdout,"\nAverage persistence length:  %.2f bonds\n",
 409             sum_pers_tot);
 410
 411   /* Handle printing of internal distances. */
 412   if (outi) {
 413     fprintf(outi, "@    xaxes scale Logarithmic\n");
 414     ymax = -1;
 415     ymin = 1e300;
 416     j = index[molind[1]-1] - index[molind[0]]; /* Polymer length -1. */
 417     for (i = 0; i < j; i++) {
 418       intd[i] /= (i + 1) * frame * nmol;
 419       if (intd[i] > ymax)
 420         ymax = intd[i];
 421       if (intd[i] < ymin)
 422         ymin = intd[i];
 423     }
 424     xvgr_world(outi, 1, ymin, j, ymax,oenv);
 425     for (i = 0; i < j; i++) {
 426       fprintf(outi, "%d  %8.4f\n", i+1, intd[i]);
 427     }
 428     ffclose(outi);
 429   }
 430
 431   do_view(oenv,opt2fn("-o",NFILE,fnm),"-nxy");
 432   if (opt2bSet("-v",NFILE,fnm))
 433     do_view(oenv,opt2fn("-v",NFILE,fnm),"-nxy");
 434   if (opt2bSet("-p",NFILE,fnm))
 435     do_view(oenv,opt2fn("-p",NFILE,fnm),"-nxy");
 436
 437   thanx(stderr);
 438
 439   return 0;
 440 }