src/tools/gmx_h2order.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
  41 #include "sysstuff.h"
  42 #include <string.h>
  43 #include "typedefs.h"
  44 #include "smalloc.h"
  45 #include "macros.h"
  46 #include "princ.h"
  47 #include "rmpbc.h"
  48 #include "vec.h"
  49 #include "xvgr.h"
  50 #include "pbc.h"
  51 #include "copyrite.h"
  52 #include "futil.h"
  53 #include "statutil.h"
  54 #include "index.h"
  55 #include "gmx_ana.h"
  56
  57
  58 /****************************************************************************/
  59 /* This program calculates the ordering of water molecules across a box, as */
  60 /* function of the z-coordinate. This implies averaging over slices and over*/
  61 /* time. Output is the average cos of the angle of the dipole with the      */
  62 /* normal, per slice.                                                       */
  63 /* In addition, it calculates the average dipole moment itself in three     */
  64 /* directions.                                                              */
  65 /****************************************************************************/
  66
  67 void calc_h2order(const char *fn, atom_id index[], int ngx, rvec **slDipole,
  68                   real **slOrder, real *slWidth, int *nslices,
  69                   t_topology *top, int ePBC,
  70                   int axis, gmx_bool bMicel, atom_id micel[], int nmic,
  71                   const output_env_t oenv)
  72 {
  73   rvec *x0,              /* coordinates with pbc */
  74        dipole,           /* dipole moment due to one molecules */
  75        normal,
  76        com;              /* center of mass of micel, with bMicel */
  77   rvec *dip;             /* sum of dipoles, unnormalized */
  78   matrix box;            /* box (3x3) */
  79   t_trxstatus *status;
  80   real  t,               /* time from trajectory */
  81        *sum,             /* sum of all cosines of dipoles, per slice */
  82        *frame;           /* order over one frame */
  83   int natoms,            /* nr. atoms in trj */
  84       i,j,teller = 0,
  85       slice=0,           /* current slice number */
  86       *count;            /* nr. of atoms in one slice */
  87   gmx_rmpbc_t  gpbc=NULL;
  88
  89   if ((natoms = read_first_x(oenv,&status,fn,&t,&x0,box)) == 0)
  90     gmx_fatal(FARGS,"Could not read coordinates from statusfile\n");
  91
  92   if (! *nslices)
  93     *nslices = (int)(box[axis][axis] * 10); /* default value */
  94
  95
  96   switch(axis)
  97   {
  98   case 0:
  99     normal[0] = 1; normal[1] = 0; normal[2] = 0;
 100     break;
 101   case 1:
 102     normal[0] = 0; normal[1] = 1; normal[2] = 0;
 103     break;
 104   case 2:
 105     normal[0] = 0; normal[1] = 0; normal[2] = 1;
 106     break;
 107   default:
 108       gmx_fatal(FARGS,"No valid value for -axis-. Exiting.\n");
 109       /* make compiler happy */
 110       normal[0] = 1; normal[1] = 0; normal[2] = 0;
 111   }
 112
 113   clear_rvec(dipole);
 114   snew(count, *nslices);
 115   snew(sum, *nslices);
 116   snew(dip, *nslices);
 117   snew(frame, *nslices);
 118
 119   *slWidth = box[axis][axis]/(*nslices);
 120   fprintf(stderr,"Box divided in %d slices. Initial width of slice: %f\n",
 121           *nslices, *slWidth);
 122
 123   teller = 0;
 124
 125   gpbc = gmx_rmpbc_init(&top->idef,ePBC,natoms,box);
 126   /*********** Start processing trajectory ***********/
 127   do
 128   {
 129     *slWidth = box[axis][axis]/(*nslices);
 130     teller++;
 131
 132     gmx_rmpbc(gpbc,natoms,box,x0);
 133
 134     if (bMicel)
 135       calc_xcm(x0, nmic, micel, top->atoms.atom, com, FALSE);
 136
 137     for (i = 0; i < ngx/3; i++)
 138     {
 139       /* put all waters in box */
 140       for (j = 0; j < DIM; j++)
 141       {
 142         if (x0[index[3*i]][j] < 0)
 143         {
 144           x0[index[3*i]][j] += box[j][j];
 145           x0[index[3*i+1]][j] += box[j][j];
 146           x0[index[3*i+2]][j] += box[j][j];
 147         }
 148         if (x0[index[3*i]][j] > box[j][j])
 149         {
 150           x0[index[3*i]][j] -= box[j][j];
 151           x0[index[3*i+1]][j] -= box[j][j];
 152           x0[index[3*i+2]][j] -= box[j][j];
 153         }
 154       }
 155
 156       for (j = 0; j < DIM; j++)
 157         dipole[j] =
 158           x0[index[3*i]][j] * top->atoms.atom[index[3*i]].q +
 159           x0[index[3*i+1]][j] * top->atoms.atom[index[3*i+1]].q +
 160           x0[index[3*i+2]][j] * top->atoms.atom[index[3*i+2]].q;
 161
 162       /* now we have a dipole vector. Might as well safe it. Then the
 163          rest depends on whether we're dealing with
 164          a flat or a spherical interface.
 165        */
 166
 167       if (bMicel)
 168       { /* this is for spherical interfaces */
 169         rvec_sub(com, x0[index[3*i]], normal); /* vector from Oxygen to COM */
 170         slice = norm(normal)/(*slWidth);         /* spherical slice           */
 171
 172         sum[slice] += iprod(dipole, normal) / (norm(dipole) * norm(normal));
 173         frame[slice] += iprod(dipole, normal) / (norm(dipole) * norm(normal));
 174         count[slice]++;
 175
 176       } else {
 177         /* this is for flat interfaces      */
 178
 179         /* determine which slice atom is in */
 180         slice = (x0[index[3*i]][axis] / (*slWidth));
 181         if (slice < 0 || slice >= *nslices)
 182         {
 183           fprintf(stderr,"Coordinate: %f ",x0[index[3*i]][axis]);
 184           fprintf(stderr,"HELP PANIC! slice = %d, OUT OF RANGE!\n",slice);
 185         }
 186         else
 187         {
 188           rvec_add(dipole, dip[slice], dip[slice]);
 189           /* Add dipole to total. mag[slice] is total dipole in axis direction */
 190           sum[slice] += iprod(dipole,normal)/norm(dipole);
 191           frame[slice] += iprod(dipole,normal)/norm(dipole);
 192           /* increase count for that slice */
 193           count[slice]++;
 194         }
 195       }
 196     }
 197
 198   } while (read_next_x(oenv,status,&t,natoms,x0,box));
 199   /*********** done with status file **********/
 200
 201   fprintf(stderr,"\nRead trajectory. Printing parameters to file\n");
 202   gmx_rmpbc_done(gpbc);
 203
 204   for (i = 0; i < *nslices; i++)  /* average over frames */
 205   {
 206     fprintf(stderr,"%d waters in slice %d\n",count[i],i);
 207     if (count[i] > 0)   /* divide by number of molecules in each slice */
 208       {
 209         sum[i] = sum[i] / count[i];
 210         dip[i][XX] = dip[i][XX] / count[i];
 211         dip[i][YY] = dip[i][YY] / count[i];
 212         dip[i][ZZ] = dip[i][ZZ] / count[i];
 213       }
 214     else
 215       fprintf(stderr,"No water in slice %d\n",i);
 216   }
 217
 218   *slOrder = sum;  /* copy a pointer, I hope */
 219   *slDipole = dip;
 220   sfree(x0);  /* free memory used by coordinate arrays */
 221 }
 222
 223 void h2order_plot(rvec dipole[], real order[], const char *afile,
 224                   int nslices, real slWidth, const output_env_t oenv)
 225 {
 226   FILE       *ord;                /* xvgr files with order parameters  */
 227   int        slice;               /* loop index     */
 228   char       buf[256];            /* for xvgr title */
 229   real       factor;              /* conversion to Debye from electron*nm */
 230
 231   /*  factor = 1e-9*1.60217733e-19/3.336e-30 */
 232   factor = 1.60217733/3.336e-2;
 233   fprintf(stderr,"%d slices\n",nslices);
 234   sprintf(buf,"Water orientation with respect to normal");
 235   ord = xvgropen(afile,buf,
 236                  "box (nm)","mu_x, mu_y, mu_z (D), cosine with normal",oenv);
 237
 238   for (slice = 0; slice < nslices; slice++)
 239     fprintf(ord,"%8.3f %8.3f %8.3f %8.3f %e\n", slWidth*slice,
 240             factor*dipole[slice][XX], factor*dipole[slice][YY],
 241             factor*dipole[slice][ZZ], order[slice]);
 242
 243   ffclose(ord);
 244 }
 245
 246 int gmx_h2order(int argc,char *argv[])
 247 {
 248   const char *desc[] = {
 249     "[TT]g_h2order[tt] computes the orientation of water molecules with respect to the normal",
 250     "of the box. The program determines the average cosine of the angle",
 251     "between the dipole moment of water and an axis of the box. The box is",
 252     "divided in slices and the average orientation per slice is printed.",
 253     "Each water molecule is assigned to a slice, per time frame, based on the",
 254     "position of the oxygen. When [TT]-nm[tt] is used, the angle between the water",
 255     "dipole and the axis from the center of mass to the oxygen is calculated",
 256     "instead of the angle between the dipole and a box axis."
 257   };
 258   static int  axis = 2;                       /* normal to memb. default z  */
 259   static const char *axtitle="Z";
 260   static int  nslices = 0;                    /* nr of slices defined       */
 261   t_pargs pa[] = {
 262     { "-d",   FALSE, etSTR, {&axtitle},
 263       "Take the normal on the membrane in direction X, Y or Z." },
 264     { "-sl",  FALSE, etINT, {&nslices},
 265       "Calculate order parameter as function of boxlength, dividing the box"
 266       " in this number of slices."}
 267   };
 268   const char *bugs[] = {
 269     "The program assigns whole water molecules to a slice, based on the first "
 270     "atom of three in the index file group. It assumes an order O,H,H. "
 271     "Name is not important, but the order is. If this demand is not met, "
 272     "assigning molecules to slices is different."
 273   };
 274
 275   output_env_t oenv;
 276   real      *slOrder,                       /* av. cosine, per slice      */
 277             slWidth = 0.0;                  /* width of a slice           */
 278   rvec      *slDipole;
 279   char      *grpname,                       /* groupnames                 */
 280             *micname;
 281   int       ngx,                            /* nr. of atomsin sol group   */
 282             nmic=0;                           /* nr. of atoms in micelle    */
 283   t_topology *top;                          /* topology                   */
 284   int       ePBC;
 285   atom_id    *index,                        /* indices for solvent group  */
 286              *micelle=NULL;
 287   gmx_bool       bMicel =  FALSE;               /* think we're a micel        */
 288   t_filenm  fnm[] = {                       /* files for g_order          */
 289     { efTRX, "-f", NULL,  ffREAD },         /* trajectory file            */
 290     { efNDX, NULL, NULL,  ffREAD },         /* index file                 */
 291     { efNDX, "-nm", NULL, ffOPTRD },        /* index with micelle atoms   */
 292     { efTPX, NULL, NULL,  ffREAD },         /* topology file              */
 293     { efXVG,"-o",  "order", ffWRITE },      /* xvgr output file           */
 294   };
 295
 296 #define NFILE asize(fnm)
 297
 298   CopyRight(stderr,argv[0]);
 299   parse_common_args(&argc, argv,
 300                     PCA_CAN_VIEW | PCA_CAN_TIME | PCA_BE_NICE, NFILE,
 301                     fnm, asize(pa),pa,asize(desc),desc,asize(bugs),bugs,&oenv);
 302   bMicel = opt2bSet("-nm",NFILE,fnm);
 303
 304   top = read_top(ftp2fn(efTPX,NFILE,fnm),&ePBC);   /* read topology file */
 305
 306   rd_index(ftp2fn(efNDX,NFILE,fnm),1,&ngx,&index,&grpname);
 307
 308   if (bMicel)
 309     rd_index(opt2fn("-nm",NFILE,fnm), 1, &nmic, &micelle, &micname);
 310
 311   calc_h2order(ftp2fn(efTRX,NFILE,fnm), index, ngx, &slDipole, &slOrder,
 312                &slWidth, &nslices, top, ePBC, axis, bMicel, micelle, nmic,
 313                oenv);
 314
 315   h2order_plot(slDipole, slOrder, opt2fn("-o",NFILE,fnm), nslices,
 316                slWidth, oenv);
 317
 318   do_view(oenv,opt2fn("-o",NFILE,fnm),"-nxy");      /* view xvgr file */
 319   thanx(stderr);
 320
 321   return 0;
 322 }
 323
 324
 325
 326
 327
 328
 329
 330
 331