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