src/tools/gmx_bundle.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 "statutil.h"
  42 #include "sysstuff.h"
  43 #include "typedefs.h"
  44 #include "smalloc.h"
  45 #include "macros.h"
  46 #include "vec.h"
  47 #include "copyrite.h"
  48 #include "futil.h"
  49 #include "statutil.h"
  50 #include "index.h"
  51 #include "xvgr.h"
  52 #include "rmpbc.h"
  53 #include "tpxio.h"
  54 #include "physics.h"
  55 #include "gmx_ana.h"
  56
  57
  58 #define MAX_ENDS 3
  59
  60 typedef struct{
  61   int n;
  62   int nend;
  63   rvec *end[MAX_ENDS];
  64   rvec *mid;
  65   rvec *dir;
  66   real *len;
  67 } t_bundle;
  68
  69 static void rotate_ends(t_bundle *bun,rvec axis,int c0,int c1)
  70 {
  71   int  end,i;
  72   rvec ax,tmp;
  73
  74   unitv(axis,ax);
  75   for(end=0; end<bun->nend; end++)
  76     for(i=0; i<bun->n; i++) {
  77       copy_rvec(bun->end[end][i],tmp);
  78       bun->end[end][i][c0] = ax[c1]*tmp[c0] - ax[c0]*tmp[c1];
  79       bun->end[end][i][c1] = ax[c0]*tmp[c0] + ax[c1]*tmp[c1];
  80     }
  81   copy_rvec(axis,tmp);
  82   axis[c0] = ax[c1]*tmp[c0] - ax[c0]*tmp[c1];
  83   axis[c1] = ax[c0]*tmp[c0] + ax[c1]*tmp[c1];
  84 }
  85
  86 static void calc_axes(rvec x[],t_atom atom[],int gnx[],atom_id *index[],
  87                       gmx_bool bRot,t_bundle *bun)
  88 {
  89   int  end,i,div,d;
  90   real *mtot,m;
  91   rvec axis[MAX_ENDS],cent;
  92
  93   snew(mtot,bun->n);
  94
  95   for(end=0; end<bun->nend; end++) {
  96     for(i=0; i<bun->n; i++) {
  97       clear_rvec(bun->end[end][i]);
  98       mtot[i] = 0;
  99     }
 100     div = gnx[end]/bun->n;
 101     for(i=0; i<gnx[end]; i++) {
 102       m = atom[index[end][i]].m;
 103       for(d=0; d<DIM; d++)
 104         bun->end[end][i/div][d] += m*x[index[end][i]][d];
 105       mtot[i/div] += m;
 106     }
 107     clear_rvec(axis[end]);
 108     for(i=0; i<bun->n; i++) {
 109       svmul(1.0/mtot[i],bun->end[end][i],bun->end[end][i]);
 110       rvec_inc(axis[end],bun->end[end][i]);
 111     }
 112     svmul(1.0/bun->n,axis[end],axis[end]);
 113   }
 114   sfree(mtot);
 115
 116   rvec_add(axis[0],axis[1],cent);
 117   svmul(0.5,cent,cent);
 118   /* center the bundle on the origin */
 119   for(end=0; end<bun->nend; end++) {
 120     rvec_dec(axis[end],cent);
 121     for(i=0; i<bun->n; i++)
 122       rvec_dec(bun->end[end][i],cent);
 123   }
 124   if (bRot) {
 125     /* rotate the axis parallel to the z-axis */
 126     rotate_ends(bun,axis[0],YY,ZZ);
 127     rotate_ends(bun,axis[0],XX,ZZ);
 128   }
 129   for(i=0; i<bun->n; i++) {
 130     rvec_add(bun->end[0][i],bun->end[1][i],bun->mid[i]);
 131     svmul(0.5,bun->mid[i],bun->mid[i]);
 132     rvec_sub(bun->end[0][i],bun->end[1][i],bun->dir[i]);
 133     bun->len[i] = norm(bun->dir[i]);
 134     unitv(bun->dir[i],bun->dir[i]);
 135   }
 136 }
 137
 138 static void dump_axes(t_trxstatus *status,t_trxframe *fr,t_atoms *outat,
 139                       t_bundle *bun)
 140 {
 141   t_trxframe  frout;
 142   static rvec *xout=NULL;
 143   int         i,end;
 144
 145   if (xout==NULL)
 146     snew(xout,outat->nr);
 147
 148   for(i=0; i<bun->n; i++) {
 149     copy_rvec(bun->end[0][i],xout[3*i]);
 150     if (bun->nend >= 3)
 151       copy_rvec(bun->end[2][i],xout[3*i+1]);
 152     else
 153       copy_rvec(bun->mid[i],xout[3*i+1]);
 154     copy_rvec(bun->end[1][i],xout[3*i+2]);
 155   }
 156   frout = *fr;
 157   frout.bV     = FALSE;
 158   frout.bF     = FALSE;
 159   frout.bBox   = FALSE;
 160   frout.bAtoms = TRUE;
 161   frout.natoms = outat->nr;
 162   frout.atoms  = outat;
 163   frout.x      = xout;
 164   write_trxframe(status,&frout,NULL);
 165 }
 166
 167 int gmx_bundle(int argc,char *argv[])
 168 {
 169   const char *desc[] = {
 170     "g_bundle analyzes bundles of axes. The axes can be for instance",
 171     "helix axes. The program reads two index groups and divides both",
 172     "of them in [TT]-na[tt] parts. The centers of mass of these parts",
 173     "define the tops and bottoms of the axes.",
 174     "Several quantities are written to file:",
 175     "the axis length, the distance and the z-shift of the axis mid-points",
 176     "with respect to the average center of all axes, the total tilt,",
 177     "the radial tilt and the lateral tilt with respect to the average axis.",
 178     "[PAR]",
 179     "With options [TT]-ok[tt], [TT]-okr[tt] and [TT]-okl[tt] the total,",
 180     "radial and lateral kinks of the axes are plotted. An extra index",
 181     "group of kink atoms is required, which is also divided into [TT]-na[tt]",
 182     "parts. The kink angle is defined as the angle between the kink-top and",
 183     "the bottom-kink vectors.",
 184     "[PAR]",
 185     "With option [TT]-oa[tt] the top, mid (or kink when [TT]-ok[tt] is set)",
 186     "and bottom points of each axis",
 187     "are written to a pdb file each frame. The residue numbers correspond",
 188     "to the axis numbers. When viewing this file with [TT]rasmol[tt], use the",
 189     "command line option [TT]-nmrpdb[tt], and type [TT]set axis true[tt] to",
 190     "display the reference axis."
 191   };
 192   static int  n=0;
 193   static gmx_bool bZ=FALSE;
 194   t_pargs pa[] = {
 195     { "-na", FALSE, etINT, {&n},
 196         "Number of axes" },
 197     { "-z", FALSE, etBOOL, {&bZ},
 198         "Use the Z-axis as reference iso the average axis" }
 199   };
 200   FILE       *out,*flen,*fdist,*fz,*ftilt,*ftiltr,*ftiltl;
 201   FILE       *fkink=NULL,*fkinkr=NULL,*fkinkl=NULL;
 202   t_trxstatus *status;
 203   t_trxstatus *fpdb;
 204   t_topology top;
 205   int        ePBC;
 206   rvec       *xtop;
 207   matrix     box;
 208   t_trxframe fr;
 209   t_atoms    outatoms;
 210   real       t,comp;
 211   int        natoms;
 212   char       *grpname[MAX_ENDS],title[256];
 213   /* FIXME: The constness should not be cast away */
 214   char       *anm=(char *)"CA",*rnm=(char *)"GLY";
 215   int        i,j,gnx[MAX_ENDS];
 216   atom_id    *index[MAX_ENDS];
 217   t_bundle   bun;
 218   gmx_bool       bKink;
 219   rvec       va,vb,vc,vr,vl;
 220   output_env_t oenv;
 221   gmx_rmpbc_t  gpbc=NULL;
 222
 223 #define NLEG asize(leg)
 224   t_filenm fnm[] = {
 225     { efTRX, "-f", NULL, ffREAD },
 226     { efTPS, NULL, NULL, ffREAD },
 227     { efNDX, NULL, NULL, ffOPTRD },
 228     { efXVG, "-ol", "bun_len", ffWRITE },
 229     { efXVG, "-od", "bun_dist", ffWRITE },
 230     { efXVG, "-oz", "bun_z", ffWRITE },
 231     { efXVG, "-ot", "bun_tilt", ffWRITE },
 232     { efXVG, "-otr", "bun_tiltr", ffWRITE },
 233     { efXVG, "-otl", "bun_tiltl", ffWRITE },
 234     { efXVG, "-ok", "bun_kink", ffOPTWR },
 235     { efXVG, "-okr", "bun_kinkr", ffOPTWR },
 236     { efXVG, "-okl", "bun_kinkl", ffOPTWR },
 237     { efPDB, "-oa", "axes", ffOPTWR }
 238   };
 239 #define NFILE asize(fnm)
 240
 241   CopyRight(stderr,argv[0]);
 242   parse_common_args(&argc,argv,PCA_CAN_TIME | PCA_TIME_UNIT | PCA_BE_NICE,
 243                     NFILE,fnm,asize(pa),pa,asize(desc),desc,0,NULL,&oenv);
 244
 245   read_tps_conf(ftp2fn(efTPS,NFILE,fnm),title,&top,&ePBC,&xtop,NULL,box,TRUE);
 246
 247   bKink = opt2bSet("-ok",NFILE,fnm) || opt2bSet("-okr",NFILE,fnm)
 248     || opt2bSet("-okl",NFILE,fnm);
 249   if (bKink)
 250     bun.nend = 3;
 251   else
 252     bun.nend = 2;
 253
 254   fprintf(stderr,"Select a group of top and a group of bottom ");
 255   if (bKink)
 256     fprintf(stderr,"and a group of kink ");
 257   fprintf(stderr,"atoms\n");
 258   get_index(&top.atoms,ftp2fn_null(efNDX,NFILE,fnm),bun.nend,
 259             gnx,index,grpname);
 260
 261   if (n<=0 || gnx[0] % n || gnx[1] % n || (bKink && gnx[2] % n))
 262     gmx_fatal(FARGS,
 263                 "The size of one of your index groups is not a multiple of n");
 264   bun.n = n;
 265   snew(bun.end[0],n);
 266   snew(bun.end[1],n);
 267   if (bKink)
 268     snew(bun.end[2],n);
 269   snew(bun.mid,n);
 270   snew(bun.dir,n);
 271   snew(bun.len,n);
 272
 273   flen   = xvgropen(opt2fn("-ol",NFILE,fnm),"Axis lengths",
 274                     output_env_get_xvgr_tlabel(oenv),"(nm)",oenv);
 275   fdist  = xvgropen(opt2fn("-od",NFILE,fnm),"Distance of axis centers",
 276                     output_env_get_xvgr_tlabel(oenv),"(nm)",oenv);
 277   fz     = xvgropen(opt2fn("-oz",NFILE,fnm),"Z-shift of axis centers",
 278                     output_env_get_xvgr_tlabel(oenv),"(nm)",oenv);
 279   ftilt  = xvgropen(opt2fn("-ot",NFILE,fnm),"Axis tilts",
 280                     output_env_get_xvgr_tlabel(oenv),"(degrees)",oenv);
 281   ftiltr = xvgropen(opt2fn("-otr",NFILE,fnm),"Radial axis tilts",
 282                     output_env_get_xvgr_tlabel(oenv),"(degrees)",oenv);
 283   ftiltl = xvgropen(opt2fn("-otl",NFILE,fnm),"Lateral axis tilts",
 284                     output_env_get_xvgr_tlabel(oenv),"(degrees)",oenv);
 285
 286   if (bKink) {
 287     fkink  = xvgropen(opt2fn("-ok",NFILE,fnm),"Kink angles",
 288                       output_env_get_xvgr_tlabel(oenv),"(degrees)",oenv);
 289     fkinkr = xvgropen(opt2fn("-okr",NFILE,fnm),"Radial kink angles",
 290                       output_env_get_xvgr_tlabel(oenv),"(degrees)",oenv);
 291     if (output_env_get_print_xvgr_codes(oenv))
 292       fprintf(fkinkr,"@ subtitle \"+ = ) (   - = ( )\"\n");
 293     fkinkl = xvgropen(opt2fn("-okl",NFILE,fnm),"Lateral kink angles",
 294                       output_env_get_xvgr_tlabel(oenv),"(degrees)",oenv);
 295   }
 296
 297   if (opt2bSet("-oa",NFILE,fnm)) {
 298     init_t_atoms(&outatoms,3*n,FALSE);
 299     outatoms.nr = 3*n;
 300     for(i=0; i<3*n; i++) {
 301       outatoms.atomname[i] = &anm;
 302       outatoms.atom[i].resind = i/3;
 303       outatoms.resinfo[i/3].name = &rnm;
 304       outatoms.resinfo[i/3].nr   = i/3 + 1;
 305       outatoms.resinfo[i/3].ic   = ' ';
 306     }
 307     fpdb = open_trx(opt2fn("-oa",NFILE,fnm),"w");
 308   } else
 309     fpdb = NULL;
 310
 311   read_first_frame(oenv,&status,ftp2fn(efTRX,NFILE,fnm),&fr,TRX_NEED_X);
 312   gpbc = gmx_rmpbc_init(&top.idef,ePBC,fr.natoms,fr.box);
 313
 314   do {
 315     gmx_rmpbc_trxfr(gpbc,&fr);
 316     calc_axes(fr.x,top.atoms.atom,gnx,index,!bZ,&bun);
 317     t = output_env_conv_time(oenv,fr.time);
 318     fprintf(flen," %10g",t);
 319     fprintf(fdist," %10g",t);
 320     fprintf(fz," %10g",t);
 321     fprintf(ftilt," %10g",t);
 322     fprintf(ftiltr," %10g",t);
 323     fprintf(ftiltl," %10g",t);
 324     if (bKink) {
 325       fprintf(fkink," %10g",t);
 326       fprintf(fkinkr," %10g",t);
 327       fprintf(fkinkl," %10g",t);
 328     }
 329
 330     for(i=0; i<bun.n; i++) {
 331       fprintf(flen," %6g",bun.len[i]);
 332       fprintf(fdist," %6g",norm(bun.mid[i]));
 333       fprintf(fz," %6g",bun.mid[i][ZZ]);
 334       fprintf(ftilt," %6g",RAD2DEG*acos(bun.dir[i][ZZ]));
 335       comp = bun.mid[i][XX]*bun.dir[i][XX]+bun.mid[i][YY]*bun.dir[i][YY];
 336       fprintf(ftiltr," %6g",RAD2DEG*
 337               asin(comp/sqrt(sqr(comp)+sqr(bun.dir[i][ZZ]))));
 338       comp = bun.mid[i][YY]*bun.dir[i][XX]-bun.mid[i][XX]*bun.dir[i][YY];
 339       fprintf(ftiltl," %6g",RAD2DEG*
 340               asin(comp/sqrt(sqr(comp)+sqr(bun.dir[i][ZZ]))));
 341       if (bKink) {
 342         rvec_sub(bun.end[0][i],bun.end[2][i],va);
 343         rvec_sub(bun.end[2][i],bun.end[1][i],vb);
 344         unitv_no_table(va,va);
 345         unitv_no_table(vb,vb);
 346         fprintf(fkink," %6g",RAD2DEG*acos(iprod(va,vb)));
 347         cprod(va,vb,vc);
 348         copy_rvec(bun.mid[i],vr);
 349         vr[ZZ] = 0;
 350         unitv_no_table(vr,vr);
 351         fprintf(fkinkr," %6g",RAD2DEG*asin(iprod(vc,vr)));
 352         vl[XX] = vr[YY];
 353         vl[YY] = -vr[XX];
 354         vl[ZZ] = 0;
 355         fprintf(fkinkl," %6g",RAD2DEG*asin(iprod(vc,vl)));
 356       }
 357     }
 358     fprintf(flen,"\n");
 359     fprintf(fdist,"\n");
 360     fprintf(fz,"\n");
 361     fprintf(ftilt,"\n");
 362     fprintf(ftiltr,"\n");
 363     fprintf(ftiltl,"\n");
 364     if (bKink) {
 365       fprintf(fkink,"\n");
 366       fprintf(fkinkr,"\n");
 367       fprintf(fkinkl,"\n");
 368     }
 369     if (fpdb )
 370       dump_axes(fpdb,&fr,&outatoms,&bun);
 371   } while(read_next_frame(oenv,status,&fr));
 372   gmx_rmpbc_done(gpbc);
 373
 374   close_trx(status);
 375
 376   if (fpdb )
 377     close_trx(fpdb);
 378   ffclose(flen);
 379   ffclose(fdist);
 380   ffclose(fz);
 381   ffclose(ftilt);
 382   ffclose(ftiltr);
 383   ffclose(ftiltl);
 384   if (bKink) {
 385     ffclose(fkink);
 386     ffclose(fkinkr);
 387     ffclose(fkinkl);
 388   }
 389
 390   thanx(stderr);
 391
 392   return 0;
 393 }