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,gmx_bool bQ,gmx_bool bRot,gmx_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     "[TT]g_gyrate[tt] computes the radius of gyration of a group of atoms",
 155     "and the radii of gyration about the [IT]x[it]-, [IT]y[it]- and [IT]z[it]-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 [IT]x-y[it] plane",
 161     "of slices along the [IT]z[it]-axis are calculated."
 162   };
 163   static int  nmol=1,nz=0;
 164   static gmx_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   t_trxstatus *status;
 179   t_topology top;
 180   int        ePBC;
 181   rvec       *x,*x_s;
 182   rvec       xcm,gvec,gvec1;
 183   matrix     box,trans;
 184   gmx_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   gmx_rmpbc_t  gpbc=NULL;
 195   const char *leg[]  = { "Rg", "RgX", "RgY", "RgZ" };
 196   const char *legI[] = { "Itot", "I1", "I2", "I3" };
 197 #define NLEG asize(leg)
 198   t_filenm fnm[] = {
 199     { efTRX, "-f",   NULL,       ffREAD },
 200     { efTPS, NULL,   NULL,       ffREAD },
 201     { efNDX, NULL,   NULL,       ffOPTRD },
 202     { efXVG, NULL,   "gyrate",   ffWRITE },
 203     { efXVG, "-acf", "moi-acf",  ffOPTWR },
 204   };
 205 #define NFILE asize(fnm)
 206   int     npargs;
 207   t_pargs *ppa;
 208
 209   CopyRight(stderr,argv[0]);
 210   npargs = asize(pa);
 211   ppa    = add_acf_pargs(&npargs,pa);
 212
 213   parse_common_args(&argc,argv,PCA_CAN_TIME | PCA_CAN_VIEW | PCA_BE_NICE,
 214                     NFILE,fnm,npargs,ppa,asize(desc),desc,0,NULL,&oenv);
 215   bACF = opt2bSet("-acf",NFILE,fnm);
 216   if (bACF && nmol!=1)
 217     gmx_fatal(FARGS,"Can only do acf with nmol=1");
 218   bRot = bRot || bMOI || bACF;
 219   /*
 220     if (nz > 0)
 221     bMOI = TRUE;
 222   */
 223   if (bRot) {
 224     printf("Will rotate system along principal axes\n");
 225     snew(moi_trans,DIM);
 226   }
 227   if (bMOI) {
 228     printf("Will print moments of inertia\n");
 229     bQ = FALSE;
 230   }
 231   if (bQ)
 232     printf("Will print radius normalised by charge\n");
 233
 234   read_tps_conf(ftp2fn(efTPS,NFILE,fnm),title,&top,&ePBC,&x,NULL,box,TRUE);
 235   get_index(&top.atoms,ftp2fn_null(efNDX,NFILE,fnm),1,&gnx,&index,&grpname);
 236
 237   if (nmol > gnx || gnx % nmol != 0) {
 238     gmx_fatal(FARGS,"The number of atoms in the group (%d) is not a multiple of nmol (%d)",gnx,nmol);
 239   }
 240   nam = gnx/nmol;
 241
 242   natoms=read_first_x(oenv,&status,ftp2fn(efTRX,NFILE,fnm),&t,&x,box);
 243   snew(x_s,natoms);
 244
 245   j  = 0;
 246   t0 = t;
 247   if (bQ)
 248     out=xvgropen(ftp2fn(efXVG,NFILE,fnm),
 249                  "Radius of Charge","Time (ps)","Rg (nm)",oenv);
 250   else if (bMOI)
 251     out=xvgropen(ftp2fn(efXVG,NFILE,fnm),
 252                  "Moments of inertia","Time (ps)","I (a.m.u. nm\\S2\\N)",oenv);
 253   else
 254     out=xvgropen(ftp2fn(efXVG,NFILE,fnm),
 255                  "Radius of gyration","Time (ps)","Rg (nm)",oenv);
 256   if (bMOI)
 257     xvgr_legend(out,NLEG,legI,oenv);
 258   else {
 259     if (bRot)
 260       if (output_env_get_print_xvgr_codes(oenv))
 261         fprintf(out,"@ subtitle \"Axes are principal component axes\"\n");
 262     xvgr_legend(out,NLEG,leg,oenv);
 263   }
 264   if (nz == 0)
 265     gpbc = gmx_rmpbc_init(&top.idef,ePBC,natoms,box);
 266   do {
 267     if (nz == 0)
 268       gmx_rmpbc_copy(gpbc,natoms,box,x,x_s);
 269     gyro = 0;
 270     clear_rvec(gvec);
 271     clear_rvec(d);
 272     for(mol=0; mol<nmol; mol++) {
 273       tm    = sub_xcm(nz==0?x_s:x,nam,index+mol*nam,top.atoms.atom,xcm,bQ);
 274       if (nz == 0)
 275         gyro += calc_gyro(x_s,nam,index+mol*nam,top.atoms.atom,
 276                           tm,gvec1,d1,bQ,bRot,bMOI,trans);
 277       else
 278         calc_gyro_z(x,box,nam,index+mol*nam,top.atoms.atom,nz,t,out);
 279       rvec_inc(gvec,gvec1);
 280       rvec_inc(d,d1);
 281     }
 282     if (nmol > 0) {
 283       gyro /= nmol;
 284       svmul(1.0/nmol,gvec,gvec);
 285       svmul(1.0/nmol,d,d);
 286     }
 287
 288     if (nz == 0) {
 289       if (bRot) {
 290         if (j >= max_moi) {
 291           max_moi += delta_moi;
 292           for(m=0; (m<DIM); m++)
 293             srenew(moi_trans[m],max_moi*DIM);
 294         }
 295         for(m=0; (m<DIM); m++)
 296           copy_rvec(trans[m],moi_trans[m]+DIM*j);
 297         fprintf(out,"%10g  %10g  %10g  %10g  %10g\n",
 298                 t,gyro,d[XX],d[YY],d[ZZ]); }
 299       else {
 300         fprintf(out,"%10g  %10g  %10g  %10g  %10g\n",
 301                 t,gyro,gvec[XX],gvec[YY],gvec[ZZ]); }
 302     }
 303     j++;
 304   } while(read_next_x(oenv,status,&t,natoms,x,box));
 305   close_trj(status);
 306   if (nz == 0)
 307     gmx_rmpbc_done(gpbc);
 308
 309   ffclose(out);
 310
 311   if (bACF) {
 312     int mode = eacVector;
 313
 314     do_autocorr(opt2fn("-acf",NFILE,fnm),oenv,
 315                 "Moment of inertia vector ACF",
 316                 j,3,moi_trans,(t-t0)/j,mode,FALSE);
 317     do_view(oenv,opt2fn("-acf",NFILE,fnm),"-nxy");
 318   }
 319
 320   do_view(oenv,ftp2fn(efXVG,NFILE,fnm),"-nxy");
 321
 322   thanx(stderr);
 323
 324   return 0;
 325 }