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