src/tools/gmx_gyrate.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 <math.h>
  40 #include <string.h>
  41
  42 #include "statutil.h"
  43 #include "sysstuff.h"
  44 #include "typedefs.h"
  45 #include "smalloc.h"
  46 #include "macros.h"
  47 #include "vec.h"
  48 #include "pbc.h"
  49 #include "copyrite.h"
  50 #include "futil.h"
  51 #include "statutil.h"
  52 #include "index.h"
  53 #include "mshift.h"
  54 #include "xvgr.h"
  55 #include "princ.h"
  56 #include "rmpbc.h"
  57 #include "txtdump.h"
  58 #include "tpxio.h"
  59 #include "gstat.h"
  60 #include "gmx_ana.h"
  61
  62
  63 real calc_gyro(rvec x[],int gnx,atom_id index[],t_atom atom[],real tm,
  64                rvec gvec,rvec d,bool bQ,bool bRot,bool bMOI,matrix trans)
  65 {
  66   int    i,ii,m;
  67   real   gyro,dx2,m0,Itot;
  68   rvec   comp;
  69
  70   if (bRot) {
  71     principal_comp(gnx,index,atom,x,trans,d);
  72     Itot = norm(d);
  73     if (bMOI)
  74       return Itot;
  75     for(m=0; (m<DIM); m++)
  76       d[m]=sqrt(d[m]/tm);
  77 #ifdef DEBUG
  78     pr_rvecs(stderr,0,"trans",trans,DIM);
  79 #endif
  80     /* rotate_atoms(gnx,index,x,trans); */
  81   }
  82   clear_rvec(comp);
  83   for(i=0; (i<gnx); i++) {
  84     ii=index[i];
  85     if (bQ)
  86       m0=fabs(atom[ii].q);
  87     else
  88       m0=atom[ii].m;
  89     for(m=0; (m<DIM); m++) {
  90       dx2=x[ii][m]*x[ii][m];
  91       comp[m]+=dx2*m0;
  92     }
  93   }
  94   gyro=comp[XX]+comp[YY]+comp[ZZ];
  95
  96   for(m=0; (m<DIM); m++)
  97     gvec[m]=sqrt((gyro-comp[m])/tm);
  98
  99   return sqrt(gyro/tm);
 100 }
 101
 102 void calc_gyro_z(rvec x[],matrix box,
 103                  int gnx,atom_id index[],t_atom atom[],
 104                  int nz,real time,FILE *out)
 105 {
 106   static dvec *inertia=NULL;
 107   static double *tm=NULL;
 108   int    i,ii,j,zi;
 109   real   zf,w,sdet,e1,e2;
 110
 111   if (inertia == NULL) {
 112     snew(inertia,nz);
 113     snew(tm,nz);
 114   }
 115
 116   for(i=0; i<nz; i++) {
 117     clear_dvec(inertia[i]);
 118     tm[i] = 0;
 119   }
 120
 121   for(i=0; (i<gnx); i++) {
 122     ii = index[i];
 123     zf = nz*x[ii][ZZ]/box[ZZ][ZZ];
 124     if (zf >= nz)
 125       zf -= nz;
 126     if (zf < 0)
 127       zf += nz;
 128     for(j=0; j<2; j++) {
 129       zi = zf + j;
 130       if (zi == nz)
 131         zi = 0;
 132       w = atom[ii].m*(1 + cos(M_PI*(zf - zi)));
 133       inertia[zi][0] += w*sqr(x[ii][YY]);
 134       inertia[zi][1] += w*sqr(x[ii][XX]);
 135       inertia[zi][2] -= w*x[ii][XX]*x[ii][YY];
 136       tm[zi] += w;
 137     }
 138   }
 139   fprintf(out,"%10g",time);
 140   for(j=0; j<nz; j++) {
 141     for(i=0; i<3; i++)
 142       inertia[j][i] /= tm[j];
 143     sdet = sqrt(sqr(inertia[j][0] - inertia[j][1]) + 4*sqr(inertia[j][2]));
 144     e1 = sqrt(0.5*(inertia[j][0] + inertia[j][1] + sdet));
 145     e2 = sqrt(0.5*(inertia[j][0] + inertia[j][1] - sdet));
 146     fprintf(out," %5.3f %5.3f",e1,e2);
 147   }
 148   fprintf(out,"\n");
 149 }
 150
 151 int gmx_gyrate(int argc,char *argv[])
 152 {
 153   const char *desc[] = {
 154     "g_gyrate computes the radius of gyration of a group of atoms",
 155     "and the radii of gyration about the x, y and z axes,",
 156     "as a function of time. The atoms are explicitly mass weighted.[PAR]",
 157     "With the [TT]-nmol[tt] option the radius of gyration will be calculated",
 158     "for multiple molecules by splitting the analysis group in equally",
 159     "sized parts.[PAR]",
 160     "With the option [TT]-nz[tt] 2D radii of gyration in the x-y plane",
 161     "of slices along the z-axis are calculated."
 162   };
 163   static int  nmol=1,nz=0;
 164   static bool bQ=FALSE,bRot=FALSE,bMOI=FALSE;
 165   t_pargs pa[] = {
 166     { "-nmol", FALSE, etINT, {&nmol},
 167       "The number of molecules to analyze" },
 168     { "-q", FALSE, etBOOL, {&bQ},
 169       "Use absolute value of the charge of an atom as weighting factor instead of mass" },
 170     { "-p", FALSE, etBOOL, {&bRot},
 171       "Calculate the radii of gyration about the principal axes." },
 172     { "-moi", FALSE, etBOOL, {&bMOI},
 173       "Calculate the moments of inertia (defined by the principal axes)." },
 174     { "-nz", FALSE, etINT, {&nz},
 175       "Calculate the 2D radii of gyration of # slices along the z-axis" },
 176   };
 177   FILE       *out;
 178   int        status;
 179   t_topology top;
 180   int        ePBC;
 181   rvec       *x,*x_s;
 182   rvec       xcm,gvec,gvec1;
 183   matrix     box,trans;
 184   bool       bACF;
 185   real       **moi_trans=NULL;
 186   int        max_moi=0,delta_moi=100;
 187   rvec       d,d1;         /* eigenvalues of inertia tensor */
 188   real       t,t0,tm,gyro;
 189   int        natoms;
 190   char       *grpname,title[256];
 191   int        i,j,m,gnx,nam,mol;
 192   atom_id    *index;
 193   output_env_t oenv;
 194   char *leg[]  = { "Rg", "RgX", "RgY", "RgZ" };
 195   char *legI[] = { "Itot", "I1", "I2", "I3" };
 196 #define NLEG asize(leg)
 197   t_filenm fnm[] = {
 198     { efTRX, "-f",   NULL,       ffREAD },
 199     { efTPS, NULL,   NULL,       ffREAD },
 200     { efNDX, NULL,   NULL,       ffOPTRD },
 201     { efXVG, NULL,   "gyrate",   ffWRITE },
 202     { efXVG, "-acf", "moi-acf",  ffOPTWR },
 203   };
 204 #define NFILE asize(fnm)
 205   int     npargs;
 206   t_pargs *ppa;
 207
 208   CopyRight(stderr,argv[0]);
 209   npargs = asize(pa);
 210   ppa    = add_acf_pargs(&npargs,pa);
 211
 212   parse_common_args(&argc,argv,PCA_CAN_TIME | PCA_CAN_VIEW | PCA_BE_NICE,
 213                     NFILE,fnm,npargs,ppa,asize(desc),desc,0,NULL,&oenv);
 214   bACF = opt2bSet("-acf",NFILE,fnm);
 215   if (bACF && nmol!=1)
 216     gmx_fatal(FARGS,"Can only do acf with nmol=1");
 217   bRot = bRot || bMOI || bACF;
 218   /*
 219     if (nz > 0)
 220     bMOI = TRUE;
 221   */
 222   if (bRot) {
 223     printf("Will rotate system along principal axes\n");
 224     snew(moi_trans,DIM);
 225   }
 226   if (bMOI) {
 227     printf("Will print moments of inertia\n");
 228     bQ = FALSE;
 229   }
 230   if (bQ)
 231     printf("Will print radius normalised by charge\n");
 232
 233   read_tps_conf(ftp2fn(efTPS,NFILE,fnm),title,&top,&ePBC,&x,NULL,box,TRUE);
 234   get_index(&top.atoms,ftp2fn_null(efNDX,NFILE,fnm),1,&gnx,&index,&grpname);
 235
 236   if (nmol > gnx || gnx % nmol != 0) {
 237     gmx_fatal(FARGS,"The number of atoms in the group (%d) is not a multiple of nmol (%d)",gnx,nmol);
 238   }
 239   nam = gnx/nmol;
 240
 241   natoms=read_first_x(oenv,&status,ftp2fn(efTRX,NFILE,fnm),&t,&x,box);
 242   snew(x_s,natoms);
 243
 244   j  = 0;
 245   t0 = t;
 246   if (bQ)
 247     out=xvgropen(ftp2fn(efXVG,NFILE,fnm),
 248                  "Radius of Charge","Time (ps)","Rg (nm)",oenv);
 249   else if (bMOI)
 250     out=xvgropen(ftp2fn(efXVG,NFILE,fnm),
 251                  "Moments of inertia","Time (ps)","I (a.m.u. nm\\S2\\N)",oenv);
 252   else
 253     out=xvgropen(ftp2fn(efXVG,NFILE,fnm),
 254                  "Radius of gyration","Time (ps)","Rg (nm)",oenv);
 255   if (bMOI)
 256     xvgr_legend(out,NLEG,legI,oenv);
 257   else {
 258     if (bRot)
 259       if (output_env_get_print_xvgr_codes(oenv))
 260         fprintf(out,"@ subtitle \"Axes are principal component axes\"\n");
 261     xvgr_legend(out,NLEG,leg,oenv);
 262   }
 263   do {
 264     if (nz == 0)
 265       rm_pbc(&top.idef,ePBC,natoms,box,x,x_s);
 266     gyro = 0;
 267     clear_rvec(gvec);
 268     clear_rvec(d);
 269     for(mol=0; mol<nmol; mol++) {
 270       tm    = sub_xcm(nz==0?x_s:x,nam,index+mol*nam,top.atoms.atom,xcm,bQ);
 271       if (nz == 0)
 272         gyro += calc_gyro(x_s,nam,index+mol*nam,top.atoms.atom,
 273                           tm,gvec1,d1,bQ,bRot,bMOI,trans);
 274       else
 275         calc_gyro_z(x,box,nam,index+mol*nam,top.atoms.atom,nz,t,out);
 276       rvec_inc(gvec,gvec1);
 277       rvec_inc(d,d1);
 278     }
 279     if (nmol > 0) {
 280       gyro /= nmol;
 281       svmul(1.0/nmol,gvec,gvec);
 282       svmul(1.0/nmol,d,d);
 283     }
 284
 285     if (nz == 0) {
 286       if (bRot) {
 287         if (j >= max_moi) {
 288           max_moi += delta_moi;
 289           for(m=0; (m<DIM); m++)
 290             srenew(moi_trans[m],max_moi*DIM);
 291         }
 292         for(m=0; (m<DIM); m++)
 293           copy_rvec(trans[m],moi_trans[m]+DIM*j);
 294         fprintf(out,"%10g  %10g  %10g  %10g  %10g\n",
 295                 t,gyro,d[XX],d[YY],d[ZZ]); }
 296       else {
 297         fprintf(out,"%10g  %10g  %10g  %10g  %10g\n",
 298                 t,gyro,gvec[XX],gvec[YY],gvec[ZZ]); }
 299     }
 300     j++;
 301   } while(read_next_x(oenv,status,&t,natoms,x,box));
 302   close_trj(status);
 303
 304   ffclose(out);
 305
 306   if (bACF) {
 307     int mode = eacVector;
 308
 309     do_autocorr(opt2fn("-acf",NFILE,fnm),oenv,
 310                 "Moment of inertia vector ACF",
 311                 j,3,moi_trans,(t-t0)/j,mode,FALSE);
 312     do_view(oenv,opt2fn("-acf",NFILE,fnm),"-nxy");
 313   }
 314
 315   do_view(oenv,ftp2fn(efXVG,NFILE,fnm),"-nxy");
 316
 317   thanx(stderr);
 318
 319   return 0;
 320 }