src/tools/g_gyrate.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 2.0
  11  *
  12  * Copyright (c) 1991-1999
  13  * BIOSON Research Institute, Dept. of Biophysical Chemistry
  14  * University of Groningen, The Netherlands
  15  *
  16  * Please refer to:
  17  * GROMACS: A message-passing parallel molecular dynamics implementation
  18  * H.J.C. Berendsen, D. van der Spoel and R. van Drunen
  19  * Comp. Phys. Comm. 91, 43-56 (1995)
  20  *
  21  * Also check out our WWW page:
  22  * http://md.chem.rug.nl/~gmx
  23  * or e-mail to:
  24  * gromacs@chem.rug.nl
  25  *
  26  * And Hey:
  27  * GRowing Old MAkes el Chrono Sweat
  28  */
  29 static char *SRCID_g_gyrate_c = "$Id$";
  30
  31 #include <math.h>
  32 #include <string.h>
  33 #include "statutil.h"
  34 #include "sysstuff.h"
  35 #include "typedefs.h"
  36 #include "smalloc.h"
  37 #include "macros.h"
  38 #include "vec.h"
  39 #include "pbc.h"
  40 #include "copyrite.h"
  41 #include "futil.h"
  42 #include "statutil.h"
  43 #include "rdgroup.h"
  44 #include "mshift.h"
  45 #include "xvgr.h"
  46 #include "princ.h"
  47 #include "rmpbc.h"
  48 #include "txtdump.h"
  49 #include "tpxio.h"
  50
  51 real calc_gyro(rvec x[],int gnx,atom_id index[],t_atom atom[],real tm,
  52                rvec gvec,rvec d,bool bQ,bool bRot)
  53 {
  54   int    i,ii,m;
  55   real   gyro,dx2,m0;
  56   matrix trans;
  57   rvec   comp;
  58
  59   if (bRot) {
  60     principal_comp(gnx,index,atom,x,trans,d);
  61     for(m=0; (m<DIM); m++)
  62       d[m]=sqrt(d[m]/tm);
  63 #ifdef DEBUG
  64     pr_rvecs(stderr,0,"trans",trans,DIM);
  65 #endif
  66     /* rotate_atoms(gnx,index,x,trans); */
  67   }
  68   clear_rvec(comp);
  69   for(i=0; (i<gnx); i++) {
  70     ii=index[i];
  71     if (bQ)
  72       m0=fabs(atom[ii].q);
  73     else
  74       m0=atom[ii].m;
  75     for(m=0; (m<DIM); m++) {
  76       dx2=x[ii][m]*x[ii][m];
  77       comp[m]+=dx2*m0;
  78     }
  79   }
  80   gyro=comp[XX]+comp[YY]+comp[ZZ];
  81
  82   for(m=0; (m<DIM); m++)
  83     gvec[m]=sqrt((gyro-comp[m])/tm);
  84
  85   return sqrt(gyro/tm);
  86 }
  87
  88 int main(int argc,char *argv[])
  89 {
  90   static char *desc[] = {
  91     "g_gyrate computes the radius of gyration of a group of atoms",
  92     "and the radii of gyration about the x, y and z axes,"
  93     "as a function of time. The atoms are explicitly mass weighted."
  94   };
  95   static bool bQ=FALSE,bRot=FALSE;
  96   t_pargs pa[] = {
  97     { "-q", FALSE, etBOOL, {&bQ},
  98       "Use absolute value of the charge of an atom as weighting factor instead of mass" },
  99     { "-p", FALSE, etBOOL, {&bRot},
 100       "Calculate the radii of gyration about the principal axes." }
 101   };
 102   FILE       *out;
 103   int        status;
 104   t_topology top;
 105   rvec       *x,*x_s;
 106   rvec       xcm,gvec;
 107   rvec       d;         /* eigenvalues of inertia tensor */
 108   matrix     box;
 109   real       t,tm,gyro;
 110   int        natoms;
 111   char       *grpname,title[256];
 112   int        i,j,gnx;
 113   atom_id    *index;
 114   char       *leg[] = { "Rg", "RgX", "RgY", "RgZ" };
 115 #define NLEG asize(leg)
 116   t_filenm fnm[] = {
 117     { efTRX, "-f", NULL, ffREAD },
 118     { efTPS, NULL, NULL, ffREAD },
 119     { efXVG, NULL, "gyrate", ffWRITE },
 120     { efNDX, NULL, NULL, ffOPTRD }
 121   };
 122 #define NFILE asize(fnm)
 123
 124   CopyRight(stderr,argv[0]);
 125   parse_common_args(&argc,argv,PCA_CAN_TIME | PCA_CAN_VIEW,TRUE,
 126                     NFILE,fnm,asize(pa),pa,asize(desc),desc,0,NULL);
 127   clear_rvec(d);
 128
 129   for(i=1; (i<argc); i++) {
 130     if (strcmp(argv[i],"-q") == 0) {
 131       bQ=TRUE;
 132       fprintf(stderr,"Will print radius normalised by charge\n");
 133     }
 134     else if (strcmp(argv[i],"-r") == 0) {
 135       bRot=TRUE;
 136       fprintf(stderr,"Will rotate system along principal axes\n");
 137     }
 138   }
 139   read_tps_conf(ftp2fn(efTPS,NFILE,fnm),title,&top,&x,NULL,box,TRUE);
 140   get_index(&top.atoms,ftp2fn_null(efNDX,NFILE,fnm),1,&gnx,&index,&grpname);
 141
 142   natoms=read_first_x(&status,ftp2fn(efTRX,NFILE,fnm),&t,&x,box);
 143   snew(x_s,natoms);
 144
 145   j=0;
 146   if (bQ)
 147     out=xvgropen(ftp2fn(efXVG,NFILE,fnm),
 148                  "Radius of Charge","Time (ps)","Rg (nm)");
 149   else
 150     out=xvgropen(ftp2fn(efXVG,NFILE,fnm),
 151                  "Radius of gyration","Time (ps)","Rg (nm)");
 152   if (bRot)
 153     fprintf(out,"@ subtitle \"Axes are principal component axes\"\n");
 154   xvgr_legend(out,NLEG,leg);
 155   do {
 156     rm_pbc(&top.idef,natoms,box,x,x_s);
 157     tm=sub_xcm(x_s,gnx,index,top.atoms.atom,xcm,bQ);
 158     gyro=calc_gyro(x_s,gnx,index,top.atoms.atom,tm,gvec,d,bQ,bRot);
 159
 160     if (bRot) {
 161       fprintf(out,"%10g  %10g  %10g  %10g  %10g\n",
 162               t,gyro,d[XX],d[YY],d[ZZ]); }
 163     else {
 164       fprintf(out,"%10g  %10g  %10g  %10g  %10g\n",
 165               t,gyro,gvec[XX],gvec[YY],gvec[ZZ]); }
 166     j++;
 167   } while(read_next_x(status,&t,natoms,x,box));
 168   close_trj(status);
 169
 170   fclose(out);
 171
 172   xvgr_file(ftp2fn(efXVG,NFILE,fnm),"-nxy");
 173
 174   thanx(stdout);
 175
 176   return 0;
 177 }