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