src/tools/g_dyndom.c

   1 /*
   2  * $Id$
   3  *
   4  *                This source code is part of
   5  *
   6  *                 G   R   O   M   A   C   S
   7  *
   8  *          GROningen MAchine for Chemical Simulations
   9  *
  10  *                        VERSION 3.1
  11  * Copyright (c) 1991-2001, University of Groningen, The Netherlands
  12  * This program is free software; you can redistribute it and/or
  13  * modify it under the terms of the GNU General Public License
  14  * as published by the Free Software Foundation; either version 2
  15  * of the License, or (at your option) any later version.
  16  *
  17  * If you want to redistribute modifications, please consider that
  18  * scientific software is very special. Version control is crucial -
  19  * bugs must be traceable. We will be happy to consider code for
  20  * inclusion in the official distribution, but derived work must not
  21  * be called official GROMACS. Details are found in the README & COPYING
  22  * files - if they are missing, get the official version at www.gromacs.org.
  23  *
  24  * To help us fund GROMACS development, we humbly ask that you cite
  25  * the papers on the package - you can find them in the top README file.
  26  *
  27  * For more info, check our website at http://www.gromacs.org
  28  *
  29  * And Hey:
  30  * Gromacs Runs One Microsecond At Cannonball Speeds
  31  */
  32
  33 #include "3dview.h"
  34 #include "statutil.h"
  35 #include "smalloc.h"
  36 #include "copyrite.h"
  37 #include "rdgroup.h"
  38 #include "confio.h"
  39 #include "fatal.h"
  40 #include "vec.h"
  41 #include "physics.h"
  42 #include "random.h"
  43
  44 static void rot_conf(t_atoms *atoms,rvec x[],rvec v[],real trans,real angle,
  45                      rvec head,rvec tail,matrix box,int isize,atom_id index[],
  46                      rvec xout[],rvec vout[])
  47 {
  48   rvec     arrow,center,xcm;
  49   real     theta,phi,arrow_len;
  50   mat4     Rx,Ry,Rz,Rinvy,Rinvz,Mtot,Tcm,Tinvcm,Tx;
  51   mat4     temp1,temp2,temp3,temp4,temp21,temp43;
  52   vec4     xv;
  53   int      i,j,ai;
  54
  55   rvec_sub(tail,head,arrow);
  56   arrow_len = norm(arrow);
  57   if (debug) {
  58     fprintf(debug,"Arrow vector:   %10.4f  %10.4f  %10.4f\n",
  59             arrow[XX],arrow[YY],arrow[ZZ]);
  60     fprintf(debug,"Effective translation %g nm\n",trans);
  61   }
  62   if (arrow_len == 0.0)
  63     fatal_error(0,"Arrow vector not given");
  64
  65   /* Copy all aoms to output */
  66   for(i=0; (i<atoms->nr);i ++) {
  67     copy_rvec(x[i],xout[i]);
  68     copy_rvec(v[i],vout[i]);
  69   }
  70
  71   /* Compute center of mass and move atoms there */
  72   clear_rvec(xcm);
  73   for(i=0; (i<isize); i++)
  74     rvec_inc(xcm,x[index[i]]);
  75   for(i=0; (i<DIM); i++)
  76     xcm[i] /= isize;
  77   if (debug)
  78     fprintf(debug,"Center of mass: %10.4f  %10.4f  %10.4f\n",
  79             xcm[XX],xcm[YY],xcm[ZZ]);
  80   for(i=0; (i<isize); i++)
  81     rvec_sub(x[index[i]],xcm,xout[index[i]]);
  82
  83   /* Compute theta and phi that describe the arrow */
  84   theta = acos(arrow[ZZ]/arrow_len);
  85   phi   = atan2(arrow[YY]/arrow_len,arrow[XX]/arrow_len);
  86   if (debug)
  87     fprintf(debug,"Phi = %.1f, Theta = %.1f\n",RAD2DEG*phi,RAD2DEG*theta);
  88
  89   /* Now the total rotation matrix: */
  90   /* Rotate a couple of times */
  91   rotate(ZZ,-phi,Rz);
  92   rotate(YY,M_PI/2-theta,Ry);
  93   rotate(XX,angle*DEG2RAD,Rx);
  94   Rx[WW][XX] = trans;
  95   rotate(YY,theta-M_PI/2,Rinvy);
  96   rotate(ZZ,phi,Rinvz);
  97
  98   mult_matrix(temp1,Ry,Rz);
  99   mult_matrix(temp2,Rinvy,Rx);
 100   mult_matrix(temp3,temp2,temp1);
 101   mult_matrix(Mtot,Rinvz,temp3);
 102
 103   print_m4(debug,"Rz",Rz);
 104   print_m4(debug,"Ry",Ry);
 105   print_m4(debug,"Rx",Rx);
 106   print_m4(debug,"Rinvy",Rinvy);
 107   print_m4(debug,"Rinvz",Rinvz);
 108   print_m4(debug,"Mtot",Mtot);
 109
 110   for(i=0; (i<isize); i++) {
 111     ai = index[i];
 112     m4_op(Mtot,xout[ai],xv);
 113     rvec_add(xv,xcm,xout[ai]);
 114     m4_op(Mtot,v[ai],xv);
 115     copy_rvec(xv,vout[ai]);
 116   }
 117 }
 118
 119 int main(int argc,char *argv[])
 120 {
 121   char *desc[] = {
 122     "g_dyndom reads a pdb file output from DynDom",
 123     "http://md.chem.rug.nl/~steve/DynDom/dyndom.home.html",
 124     "It reads the coordinates, and the coordinates of the rotation axis",
 125     "furthermore it reads an index file containing the domains.",
 126     "Furthermore it takes the first and last atom of the arrow file",
 127     "as command line arguments (head and tail) and",
 128     "finally it takes the translation vector (given in DynDom info file)",
 129     "and the angle of rotation (also as command line arguments). If the angle",
 130     "determined by DynDom is given, one should be able to recover the",
 131     "second structure used for generating the DynDom output.",
 132     "Because of limited numerical accuracy this should be verified by",
 133     "computing an all-atom RMSD (using [TT]g_confrms[tt]) rather than by file",
 134     "comparison (using diff).[PAR]",
 135     "The purpose of this program is to interpolate and extrapolate the",
 136     "rotation as found by DynDom. As a result unphysical structures with",
 137     "long or short bonds, or overlapping atoms may be produced. Visual",
 138     "inspection, and energy minimization may be necessary to",
 139     "validate the structure."
 140   };
 141   static real trans0 = 0;
 142   static rvec head  = { 0,0,0 };
 143   static rvec tail  = { 0,0,0 };
 144   static real angle0 = 0,angle1 = 0, maxangle = 0;
 145   static int  label = 0,nframes=11;
 146   t_pargs pa[] = {
 147     { "-firstangle",    FALSE, etREAL, {&angle0},
 148       "Angle of rotation about rotation vector" },
 149     { "-lastangle",    FALSE, etREAL, {&angle1},
 150       "Angle of rotation about rotation vector" },
 151     { "-nframe",   FALSE, etINT,  {&nframes},
 152       "Number of steps on the pathway" },
 153     { "-maxangle", FALSE, etREAL, {&maxangle},
 154       "DymDom dtermined angle of rotation about rotation vector" },
 155     { "-trans",    FALSE, etREAL, {&trans0},
 156       "Translation (Aangstroem) along rotation vector (see DynDom info file)" },
 157     { "-head",     FALSE, etRVEC, {head},
 158       "First atom of the arrow vector" },
 159     { "-tail",     FALSE, etRVEC, {tail},
 160       "Last atom of the arrow vector" }
 161   };
 162   int     i,j,natoms,isize,status;
 163   atom_id *index=NULL,*index_all;
 164   char    title[256],*grpname;
 165   t_atoms atoms;
 166   real    angle,trans;
 167   rvec    *x,*v,*xout,*vout;
 168   matrix  box;
 169
 170   t_filenm fnm[] = {
 171     { efPDB, "-f", "dyndom",  ffREAD },
 172     { efTRX, "-o", "rotated", ffWRITE },
 173     { efNDX, "-n", "domains", ffREAD }
 174   };
 175 #define NFILE asize(fnm)
 176
 177   CopyRight(stderr,argv[0]);
 178
 179   parse_common_args(&argc,argv,0,NFILE,fnm,asize(pa),pa,
 180                     asize(desc),desc,0,NULL);
 181
 182   get_stx_coordnum (opt2fn("-f",NFILE,fnm),&natoms);
 183   init_t_atoms(&atoms,natoms,TRUE);
 184   snew(x,natoms);
 185   snew(v,natoms);
 186   read_stx_conf(opt2fn("-f",NFILE,fnm),title,&atoms,x,v,box);
 187   snew(xout,natoms);
 188   snew(vout,natoms);
 189
 190   printf("Select group to rotate:\n");
 191   rd_index(ftp2fn(efNDX,NFILE,fnm),1,&isize,&index,&grpname);
 192   printf("Going to rotate %s containg %d atoms\n",grpname,isize);
 193
 194   snew(index_all,atoms.nr);
 195   for(i=0; (i<atoms.nr); i++)
 196     index_all[i] = i;
 197
 198   status = open_trx(opt2fn("-o",NFILE,fnm),"w");
 199
 200   label = 'A';
 201   for(i=0; (i<nframes); i++,label++) {
 202     angle = angle0 + (i*(angle1-angle0))/(nframes-1);
 203     trans = trans0*0.1*angle/maxangle;
 204     printf("Frame: %2d (label %c), angle: %8.3f deg., trans: %8.3f nm\n",
 205            i,label,angle,trans);
 206     rot_conf(&atoms,x,v,trans,angle,head,tail,box,isize,index,xout,vout);
 207
 208     if (label > 'Z')
 209       label-=26;
 210     for(j=0; (j<atoms.nr); j++)
 211       atoms.atom[j].chain = label;
 212
 213     write_trx(status,atoms.nr,index_all,&atoms,i,angle,box,xout,vout);
 214   }
 215   close_trx(status);
 216
 217   thanx(stderr);
 218
 219   return 0;
 220 }