src/tools/gmx_sorient.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 "macros.h"
  40 #include "statutil.h"
  41 #include "smalloc.h"
  42 #include "copyrite.h"
  43 #include "gstat.h"
  44 #include "vec.h"
  45 #include "xvgr.h"
  46 #include "pbc.h"
  47 #include "index.h"
  48 #include "tpxio.h"
  49 #include "gmx_ana.h"
  50
  51
  52 static void calc_com_pbc(int nrefat,t_topology *top,rvec x[],t_pbc *pbc,
  53                          atom_id index[],rvec xref,gmx_bool bPBC,matrix box)
  54 {
  55   const real tol=1e-4;
  56   gmx_bool  bChanged;
  57   int   m,j,ai,iter;
  58   real  mass,mtot;
  59   rvec  dx,xtest;
  60
  61   /* First simple calculation */
  62   clear_rvec(xref);
  63   mtot = 0;
  64   for(m=0; (m<nrefat); m++) {
  65     ai = index[m];
  66     mass = top->atoms.atom[ai].m;
  67     for(j=0; (j<DIM); j++)
  68       xref[j] += mass*x[ai][j];
  69     mtot += mass;
  70   }
  71   svmul(1/mtot,xref,xref);
  72   /* Now check if any atom is more than half the box from the COM */
  73   if (bPBC) {
  74     iter = 0;
  75     do {
  76       bChanged = FALSE;
  77       for(m=0; (m<nrefat); m++) {
  78         ai   = index[m];
  79         mass = top->atoms.atom[ai].m/mtot;
  80         pbc_dx(pbc,x[ai],xref,dx);
  81         rvec_add(xref,dx,xtest);
  82         for(j=0; (j<DIM); j++)
  83           if (fabs(xtest[j]-x[ai][j]) > tol) {
  84             /* Here we have used the wrong image for contributing to the COM */
  85             xref[j] += mass*(xtest[j]-x[ai][j]);
  86             x[ai][j] = xtest[j];
  87             bChanged = TRUE;
  88           }
  89       }
  90       if (bChanged)
  91         printf("COM: %8.3f  %8.3f  %8.3f  iter = %d\n",xref[XX],xref[YY],xref[ZZ],iter);
  92       iter++;
  93     } while (bChanged);
  94   }
  95 }
  96
  97 int gmx_sorient(int argc,char *argv[])
  98 {
  99   t_topology top;
 100   int      ePBC;
 101   char     title[STRLEN];
 102   t_trxstatus *status;
 103   int      natoms;
 104   real     t;
 105   rvec     *xtop,*x;
 106   matrix   box;
 107
 108   FILE    *fp;
 109   int     i,j,p,sa0,sa1,sa2,n,ntot,nf,m,*hist1,*hist2,*histn,nbin1,nbin2,nrbin;
 110   real    *histi1,*histi2,invbw,invrbw;
 111   double  sum1,sum2;
 112   int     *isize,nrefgrp,nrefat;
 113   atom_id **index;
 114   char    **grpname;
 115   real    inp,outp,two_pi,nav,normfac,rmin2,rmax2,rcut,rcut2,r2,r,mass,mtot;
 116   real    c1,c2;
 117   char    str[STRLEN];
 118   gmx_bool    bTPS;
 119   rvec    xref,dx,dxh1,dxh2,outer;
 120   gmx_rmpbc_t gpbc=NULL;
 121   t_pbc   pbc;
 122   const char *legr[] = { "<cos(\\8q\\4\\s1\\N)>",
 123                          "<3cos\\S2\\N(\\8q\\4\\s2\\N)-1>" };
 124   const char *legc[] = { "cos(\\8q\\4\\s1\\N)",
 125                          "3cos\\S2\\N(\\8q\\4\\s2\\N)-1" };
 126
 127   const char *desc[] = {
 128     "[TT]g_sorient[tt] analyzes solvent orientation around solutes.",
 129     "It calculates two angles between the vector from one or more",
 130     "reference positions to the first atom of each solvent molecule:[PAR]",
 131     "[GRK]theta[grk]1: the angle with the vector from the first atom of the solvent",
 132     "molecule to the midpoint between atoms 2 and 3.[BR]",
 133     "[GRK]theta[grk]2: the angle with the normal of the solvent plane, defined by the",
 134     "same three atoms, or, when the option [TT]-v23[tt] is set, ",
 135     "the angle with the vector between atoms 2 and 3.[PAR]",
 136     "The reference can be a set of atoms or",
 137     "the center of mass of a set of atoms. The group of solvent atoms should",
 138     "consist of 3 atoms per solvent molecule.",
 139     "Only solvent molecules between [TT]-rmin[tt] and [TT]-rmax[tt] are",
 140     "considered for [TT]-o[tt] and [TT]-no[tt] each frame.[PAR]",
 141     "[TT]-o[tt]: distribtion of cos([GRK]theta[grk]1) for rmin<=r<=rmax.[PAR]",
 142     "[TT]-no[tt]: distribution of cos([GRK]theta[grk]2) for rmin<=r<=rmax.[PAR]",
 143     "[TT]-ro[tt]: <cos([GRK]theta[grk]1)> and <3cos^2([GRK]theta[grk]2)-1> as a function of the",
 144     "distance.[PAR]",
 145     "[TT]-co[tt]: the sum over all solvent molecules within distance r",
 146     "of cos([GRK]theta[grk]1) and 3cos^2([GRK]theta[grk]2)-1 as a function of r.[PAR]",
 147     "[TT]-rc[tt]: the distribution of the solvent molecules as a function of r"
 148   };
 149
 150   output_env_t oenv;
 151   static gmx_bool bCom = FALSE,bVec23=FALSE,bPBC = FALSE;
 152   static real rmin=0.0,rmax=0.5,binwidth=0.02,rbinw=0.02;
 153   t_pargs pa[] = {
 154     { "-com",  FALSE, etBOOL,  {&bCom},
 155       "Use the center of mass as the reference postion" },
 156     { "-v23",  FALSE, etBOOL,  {&bVec23},
 157       "Use the vector between atoms 2 and 3" },
 158     { "-rmin",  FALSE, etREAL, {&rmin}, "Minimum distance (nm)" },
 159     { "-rmax",  FALSE, etREAL, {&rmax}, "Maximum distance (nm)" },
 160     { "-cbin",  FALSE, etREAL, {&binwidth}, "Binwidth for the cosine" },
 161     { "-rbin",  FALSE, etREAL, {&rbinw}, "Binwidth for r (nm)" },
 162     { "-pbc",   FALSE, etBOOL, {&bPBC}, "Check PBC for the center of mass calculation. Only necessary when your reference group consists of several molecules." }
 163   };
 164
 165   t_filenm fnm[] = {
 166     { efTRX, NULL,  NULL,  ffREAD },
 167     { efTPS, NULL,  NULL,  ffREAD },
 168     { efNDX, NULL,  NULL,  ffOPTRD },
 169     { efXVG, NULL,  "sori.xvg",  ffWRITE },
 170     { efXVG, "-no", "snor.xvg",  ffWRITE },
 171     { efXVG, "-ro", "sord.xvg",  ffWRITE },
 172     { efXVG, "-co", "scum.xvg",  ffWRITE },
 173     { efXVG, "-rc", "scount.xvg",  ffWRITE }
 174   };
 175 #define NFILE asize(fnm)
 176
 177   CopyRight(stderr,argv[0]);
 178   parse_common_args(&argc,argv,PCA_CAN_TIME | PCA_CAN_VIEW | PCA_BE_NICE,
 179                     NFILE,fnm,asize(pa),pa,asize(desc),desc,0,NULL,&oenv);
 180
 181   two_pi = 2/M_PI;
 182
 183   bTPS = (opt2bSet("-s",NFILE,fnm) || !opt2bSet("-n",NFILE,fnm) || bCom);
 184   if (bTPS) {
 185     read_tps_conf(ftp2fn(efTPS,NFILE,fnm),title,&top,&ePBC,&xtop,NULL,box,
 186                   bCom);
 187   }
 188
 189   /* get index groups */
 190   printf("Select a group of reference particles and a solvent group:\n");
 191   snew(grpname,2);
 192   snew(index,2);
 193   snew(isize,2);
 194   if (bTPS) {
 195     get_index(&top.atoms,ftp2fn_null(efNDX,NFILE,fnm),2,isize,index,grpname);
 196   } else {
 197     get_index(NULL,ftp2fn(efNDX,NFILE,fnm),2,isize,index,grpname);
 198   }
 199
 200   if (bCom) {
 201     nrefgrp = 1;
 202     nrefat  = isize[0];
 203   } else {
 204     nrefgrp = isize[0];
 205     nrefat  = 1;
 206   }
 207
 208   if (isize[1] % 3)
 209     gmx_fatal(FARGS,"The number of solvent atoms (%d) is not a multiple of 3",
 210                 isize[1]);
 211
 212   /* initialize reading trajectory:                         */
 213   natoms=read_first_x(oenv,&status,ftp2fn(efTRX,NFILE,fnm),&t,&x,box);
 214
 215   rmin2 = sqr(rmin);
 216   rmax2 = sqr(rmax);
 217   rcut  = 0.99*sqrt(max_cutoff2(guess_ePBC(box),box));
 218   if (rcut == 0)
 219     rcut = 10*rmax;
 220   rcut2 = sqr(rcut);
 221
 222   invbw = 1/binwidth;
 223   nbin1 = 1+(int)(2*invbw + 0.5);
 224   nbin2 = 1+(int)(invbw + 0.5);
 225
 226   invrbw = 1/rbinw;
 227
 228   snew(hist1,nbin1);
 229   snew(hist2,nbin2);
 230   nrbin = 1+(int)(rcut/rbinw);
 231   if (nrbin == 0)
 232     nrbin = 1;
 233   snew(histi1,nrbin);
 234   snew(histi2,nrbin);
 235   snew(histn,nrbin);
 236
 237   ntot = 0;
 238   nf = 0;
 239   sum1 = 0;
 240   sum2 = 0;
 241
 242   if (bTPS) {
 243     /* make molecules whole again */
 244     gpbc = gmx_rmpbc_init(&top.idef,ePBC,natoms,box);
 245   }
 246   /* start analysis of trajectory */
 247   do {
 248     if (bTPS) {
 249       /* make molecules whole again */
 250       gmx_rmpbc(gpbc,natoms,box,x);
 251     }
 252
 253     set_pbc(&pbc,ePBC,box);
 254     n    = 0;
 255     inp  = 0;
 256     outp = 0;
 257     for(p=0; (p<nrefgrp); p++) {
 258       if (bCom)
 259         calc_com_pbc(nrefat,&top,x,&pbc,index[0],xref,bPBC,box);
 260       else
 261         copy_rvec(x[index[0][p]],xref);
 262
 263       for(m=0; m<isize[1]; m+=3) {
 264         sa0 = index[1][m];
 265         sa1 = index[1][m+1];
 266         sa2 = index[1][m+2];
 267         range_check(sa0,0,natoms);
 268         range_check(sa1,0,natoms);
 269         range_check(sa2,0,natoms);
 270         pbc_dx(&pbc,x[sa0],xref,dx);
 271         r2  = norm2(dx);
 272         if (r2 < rcut2) {
 273           r = sqrt(r2);
 274           if (!bVec23) {
 275             /* Determine the normal to the plain */
 276             rvec_sub(x[sa1],x[sa0],dxh1);
 277             rvec_sub(x[sa2],x[sa0],dxh2);
 278             rvec_inc(dxh1,dxh2);
 279             svmul(1/r,dx,dx);
 280             unitv(dxh1,dxh1);
 281             inp = iprod(dx,dxh1);
 282             cprod(dxh1,dxh2,outer);
 283             unitv(outer,outer);
 284             outp = iprod(dx,outer);
 285           } else {
 286             /* Use the vector between the 2nd and 3rd atom */
 287             rvec_sub(x[sa2],x[sa1],dxh2);
 288             unitv(dxh2,dxh2);
 289             outp = iprod(dx,dxh2)/r;
 290           }
 291           {
 292             int ii = (int)(invrbw*r);
 293             range_check(ii,0,nrbin);
 294             histi1[ii] += inp;
 295             histi2[ii] += 3*sqr(outp) - 1;
 296             histn[ii]++;
 297           }
 298           if ((r2>=rmin2) && (r2<rmax2)) {
 299             int ii1 = (int)(invbw*(inp + 1));
 300             int ii2 = (int)(invbw*fabs(outp));
 301
 302             range_check(ii1,0,nbin1);
 303             range_check(ii2,0,nbin2);
 304             hist1[ii1]++;
 305             hist2[ii2]++;
 306             sum1 += inp;
 307             sum2 += outp;
 308             n++;
 309           }
 310         }
 311       }
 312     }
 313     ntot += n;
 314     nf++;
 315
 316   }  while (read_next_x(oenv,status,&t,natoms,x,box));
 317
 318   /* clean up */
 319   sfree(x);
 320   close_trj(status);
 321   gmx_rmpbc_done(gpbc);
 322
 323   /* Add the bin for the exact maximum to the previous bin */
 324   hist1[nbin1-1] += hist1[nbin1];
 325   hist2[nbin2-1] += hist2[nbin2];
 326
 327   nav     = (real)ntot/(nrefgrp*nf);
 328   normfac = invbw/ntot;
 329
 330   fprintf(stderr,  "Average nr of molecules between %g and %g nm: %.1f\n",
 331           rmin,rmax,nav);
 332   if (ntot > 0) {
 333     sum1 /= ntot;
 334     sum2 /= ntot;
 335     fprintf(stderr,"Average cos(theta1)     between %g and %g nm: %6.3f\n",
 336             rmin,rmax,sum1);
 337     fprintf(stderr,"Average 3cos2(theta2)-1 between %g and %g nm: %6.3f\n",
 338             rmin,rmax,sum2);
 339   }
 340
 341   sprintf(str,"Solvent orientation between %g and %g nm",rmin,rmax);
 342   fp=xvgropen(opt2fn("-o",NFILE,fnm), str,"cos(\\8q\\4\\s1\\N)","",oenv);
 343   if (output_env_get_print_xvgr_codes(oenv))
 344     fprintf(fp,"@ subtitle \"average shell size %.1f molecules\"\n",nav);
 345   for(i=0; i<nbin1; i++) {
 346     fprintf(fp,"%g %g\n",(i+0.5)*binwidth-1,2*normfac*hist1[i]);
 347   }
 348   ffclose(fp);
 349
 350   sprintf(str,"Solvent normal orientation between %g and %g nm",rmin,rmax);
 351   fp=xvgropen(opt2fn("-no",NFILE,fnm), str,"cos(\\8q\\4\\s2\\N)","",oenv);
 352   if (output_env_get_print_xvgr_codes(oenv))
 353     fprintf(fp,"@ subtitle \"average shell size %.1f molecules\"\n",nav);
 354   for(i=0; i<nbin2; i++) {
 355     fprintf(fp,"%g %g\n",(i+0.5)*binwidth,normfac*hist2[i]);
 356   }
 357   ffclose(fp);
 358
 359
 360   sprintf(str,"Solvent orientation");
 361   fp=xvgropen(opt2fn("-ro",NFILE,fnm),str,"r (nm)","",oenv);
 362   if (output_env_get_print_xvgr_codes(oenv))
 363     fprintf(fp,"@ subtitle \"as a function of distance\"\n");
 364   xvgr_legend(fp,2,legr,oenv);
 365   for(i=0; i<nrbin; i++)
 366     fprintf(fp,"%g %g %g\n",(i+0.5)*rbinw,
 367             histn[i] ? histi1[i]/histn[i] : 0,
 368             histn[i] ? histi2[i]/histn[i] : 0);
 369   ffclose(fp);
 370
 371   sprintf(str,"Cumulative solvent orientation");
 372   fp=xvgropen(opt2fn("-co",NFILE,fnm),str,"r (nm)","",oenv);
 373   if (output_env_get_print_xvgr_codes(oenv))
 374     fprintf(fp,"@ subtitle \"as a function of distance\"\n");
 375   xvgr_legend(fp,2,legc,oenv);
 376   normfac = 1.0/(nrefgrp*nf);
 377   c1 = 0;
 378   c2 = 0;
 379   fprintf(fp,"%g %g %g\n",0.0,c1,c2);
 380   for(i=0; i<nrbin; i++) {
 381     c1 += histi1[i]*normfac;
 382     c2 += histi2[i]*normfac;
 383     fprintf(fp,"%g %g %g\n",(i+1)*rbinw,c1,c2);
 384   }
 385   ffclose(fp);
 386
 387   sprintf(str,"Solvent distribution");
 388   fp=xvgropen(opt2fn("-rc",NFILE,fnm),str,"r (nm)","molecules/nm",oenv);
 389   if (output_env_get_print_xvgr_codes(oenv))
 390     fprintf(fp,"@ subtitle \"as a function of distance\"\n");
 391   normfac = 1.0/(rbinw*nf);
 392   for(i=0; i<nrbin; i++) {
 393     fprintf(fp,"%g %g\n",(i+0.5)*rbinw,histn[i]*normfac);
 394   }
 395   ffclose(fp);
 396
 397   do_view(oenv, opt2fn("-o",NFILE,fnm),NULL);
 398   do_view(oenv, opt2fn("-no",NFILE,fnm),NULL);
 399   do_view(oenv, opt2fn("-ro",NFILE,fnm),"-nxy");
 400   do_view(oenv, opt2fn("-co",NFILE,fnm),"-nxy");
 401
 402   thanx(stderr);
 403
 404   return 0;
 405 }