src/tools/gmx_dyndom.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
  39 #include "3dview.h"
  40 #include "statutil.h"
  41 #include "smalloc.h"
  42 #include "copyrite.h"
  43 #include "index.h"
  44 #include "confio.h"
  45 #include "gmx_fatal.h"
  46 #include "vec.h"
  47 #include "physics.h"
  48 #include "random.h"
  49 #include "gmx_ana.h"
  50
  51
  52 static void rot_conf(t_atoms *atoms,rvec x[],rvec v[],real trans,real angle,
  53                      rvec head,rvec tail,matrix box,int isize,atom_id index[],
  54                      rvec xout[],rvec vout[])
  55 {
  56   rvec     arrow,center,xcm;
  57   real     theta,phi,arrow_len;
  58   mat4     Rx,Ry,Rz,Rinvy,Rinvz,Mtot,Tcm,Tinvcm,Tx;
  59   mat4     temp1,temp2,temp3,temp4,temp21,temp43;
  60   vec4     xv;
  61   int      i,j,ai;
  62
  63   rvec_sub(tail,head,arrow);
  64   arrow_len = norm(arrow);
  65   if (debug) {
  66     fprintf(debug,"Arrow vector:   %10.4f  %10.4f  %10.4f\n",
  67             arrow[XX],arrow[YY],arrow[ZZ]);
  68     fprintf(debug,"Effective translation %g nm\n",trans);
  69   }
  70   if (arrow_len == 0.0)
  71     gmx_fatal(FARGS,"Arrow vector not given");
  72
  73   /* Copy all aoms to output */
  74   for(i=0; (i<atoms->nr);i ++) {
  75     copy_rvec(x[i],xout[i]);
  76     copy_rvec(v[i],vout[i]);
  77   }
  78
  79   /* Compute center of mass and move atoms there */
  80   clear_rvec(xcm);
  81   for(i=0; (i<isize); i++)
  82     rvec_inc(xcm,x[index[i]]);
  83   for(i=0; (i<DIM); i++)
  84     xcm[i] /= isize;
  85   if (debug)
  86     fprintf(debug,"Center of mass: %10.4f  %10.4f  %10.4f\n",
  87             xcm[XX],xcm[YY],xcm[ZZ]);
  88   for(i=0; (i<isize); i++)
  89     rvec_sub(x[index[i]],xcm,xout[index[i]]);
  90
  91   /* Compute theta and phi that describe the arrow */
  92   theta = acos(arrow[ZZ]/arrow_len);
  93   phi   = atan2(arrow[YY]/arrow_len,arrow[XX]/arrow_len);
  94   if (debug)
  95     fprintf(debug,"Phi = %.1f, Theta = %.1f\n",RAD2DEG*phi,RAD2DEG*theta);
  96
  97   /* Now the total rotation matrix: */
  98   /* Rotate a couple of times */
  99   rotate(ZZ,-phi,Rz);
 100   rotate(YY,M_PI/2-theta,Ry);
 101   rotate(XX,angle*DEG2RAD,Rx);
 102   Rx[WW][XX] = trans;
 103   rotate(YY,theta-M_PI/2,Rinvy);
 104   rotate(ZZ,phi,Rinvz);
 105
 106   mult_matrix(temp1,Ry,Rz);
 107   mult_matrix(temp2,Rinvy,Rx);
 108   mult_matrix(temp3,temp2,temp1);
 109   mult_matrix(Mtot,Rinvz,temp3);
 110
 111   print_m4(debug,"Rz",Rz);
 112   print_m4(debug,"Ry",Ry);
 113   print_m4(debug,"Rx",Rx);
 114   print_m4(debug,"Rinvy",Rinvy);
 115   print_m4(debug,"Rinvz",Rinvz);
 116   print_m4(debug,"Mtot",Mtot);
 117
 118   for(i=0; (i<isize); i++) {
 119     ai = index[i];
 120     m4_op(Mtot,xout[ai],xv);
 121     rvec_add(xv,xcm,xout[ai]);
 122     m4_op(Mtot,v[ai],xv);
 123     copy_rvec(xv,vout[ai]);
 124   }
 125 }
 126
 127 int gmx_dyndom(int argc,char *argv[])
 128 {
 129   const char *desc[] = {
 130     "g_dyndom reads a pdb file output from DynDom",
 131     "(http://www.cmp.uea.ac.uk/dyndom/).",
 132     "It reads the coordinates, the coordinates of the rotation axis,",
 133     "and an index file containing the domains.",
 134     "Furthermore it takes the first and last atom of the arrow file",
 135     "as command line arguments (head and tail) and",
 136     "finally it takes the translation vector (given in DynDom info file)",
 137     "and the angle of rotation (also as command line arguments). If the angle",
 138     "determined by DynDom is given, one should be able to recover the",
 139     "second structure used for generating the DynDom output.",
 140     "Because of limited numerical accuracy this should be verified by",
 141     "computing an all-atom RMSD (using [TT]g_confrms[tt]) rather than by file",
 142     "comparison (using diff).[PAR]",
 143     "The purpose of this program is to interpolate and extrapolate the",
 144     "rotation as found by DynDom. As a result unphysical structures with",
 145     "long or short bonds, or overlapping atoms may be produced. Visual",
 146     "inspection, and energy minimization may be necessary to",
 147     "validate the structure."
 148   };
 149   static real trans0 = 0;
 150   static rvec head  = { 0,0,0 };
 151   static rvec tail  = { 0,0,0 };
 152   static real angle0 = 0,angle1 = 0, maxangle = 0;
 153   static int  label = 0,nframes=11;
 154   t_pargs pa[] = {
 155     { "-firstangle",    FALSE, etREAL, {&angle0},
 156       "Angle of rotation about rotation vector" },
 157     { "-lastangle",    FALSE, etREAL, {&angle1},
 158       "Angle of rotation about rotation vector" },
 159     { "-nframe",   FALSE, etINT,  {&nframes},
 160       "Number of steps on the pathway" },
 161     { "-maxangle", FALSE, etREAL, {&maxangle},
 162       "DymDom dtermined angle of rotation about rotation vector" },
 163     { "-trans",    FALSE, etREAL, {&trans0},
 164       "Translation (Aangstroem) along rotation vector (see DynDom info file)" },
 165     { "-head",     FALSE, etRVEC, {head},
 166       "First atom of the arrow vector" },
 167     { "-tail",     FALSE, etRVEC, {tail},
 168       "Last atom of the arrow vector" }
 169   };
 170   int     i,j,natoms,isize;
 171   t_trxstatus *status;
 172   atom_id *index=NULL,*index_all;
 173   char    title[256],*grpname;
 174   t_atoms atoms;
 175   real    angle,trans;
 176   rvec    *x,*v,*xout,*vout;
 177   matrix  box;
 178   output_env_t oenv;
 179
 180   t_filenm fnm[] = {
 181     { efPDB, "-f", "dyndom",  ffREAD },
 182     { efTRO, "-o", "rotated", ffWRITE },
 183     { efNDX, "-n", "domains", ffREAD }
 184   };
 185 #define NFILE asize(fnm)
 186
 187   CopyRight(stderr,argv[0]);
 188
 189   parse_common_args(&argc,argv,0,NFILE,fnm,asize(pa),pa,
 190                     asize(desc),desc,0,NULL,&oenv);
 191
 192   get_stx_coordnum (opt2fn("-f",NFILE,fnm),&natoms);
 193   init_t_atoms(&atoms,natoms,TRUE);
 194   snew(x,natoms);
 195   snew(v,natoms);
 196   read_stx_conf(opt2fn("-f",NFILE,fnm),title,&atoms,x,v,NULL,box);
 197   snew(xout,natoms);
 198   snew(vout,natoms);
 199
 200   printf("Select group to rotate:\n");
 201   rd_index(ftp2fn(efNDX,NFILE,fnm),1,&isize,&index,&grpname);
 202   printf("Going to rotate %s containg %d atoms\n",grpname,isize);
 203
 204   snew(index_all,atoms.nr);
 205   for(i=0; (i<atoms.nr); i++)
 206     index_all[i] = i;
 207
 208   status = open_trx(opt2fn("-o",NFILE,fnm),"w");
 209
 210   label = 'A';
 211   for(i=0; (i<nframes); i++,label++) {
 212     angle = angle0 + (i*(angle1-angle0))/(nframes-1);
 213     trans = trans0*0.1*angle/maxangle;
 214     printf("Frame: %2d (label %c), angle: %8.3f deg., trans: %8.3f nm\n",
 215            i,label,angle,trans);
 216     rot_conf(&atoms,x,v,trans,angle,head,tail,box,isize,index,xout,vout);
 217
 218     if (label > 'Z')
 219       label-=26;
 220     for(j=0; (j<atoms.nr); j++)
 221       atoms.resinfo[atoms.atom[j].resind].chainid = label;
 222
 223     write_trx(status,atoms.nr,index_all,&atoms,i,angle,box,xout,vout,NULL);
 224   }
 225   close_trx(status);
 226
 227   thanx(stderr);
 228
 229   return 0;
 230 }