src/tools/g_h2order.c

   1 /*
   2  * $Id$
   3  *
   4  *       This source code is part of
   5  *
   6  *        G   R   O   M   A   C   S
   7  *
   8  * GROningen MAchine for Chemical Simulations
   9  *
  10  *               VERSION 2.0
  11  *
  12  * Copyright (c) 1991-1999
  13  * BIOSON Research Institute, Dept. of Biophysical Chemistry
  14  * University of Groningen, The Netherlands
  15  *
  16  * Please refer to:
  17  * GROMACS: A message-passing parallel molecular dynamics implementation
  18  * H.J.C. Berendsen, D. van der Spoel and R. van Drunen
  19  * Comp. Phys. Comm. 91, 43-56 (1995)
  20  *
  21  * Also check out our WWW page:
  22  * http://md.chem.rug.nl/~gmx
  23  * or e-mail to:
  24  * gromacs@chem.rug.nl
  25  *
  26  * And Hey:
  27  * GRowing Old MAkes el Chrono Sweat
  28  */
  29 static char *SRCID_g_h2order_c = "$Id$";
  30
  31 #include <math.h>
  32 #include "sysstuff.h"
  33 #include "string.h"
  34 #include "typedefs.h"
  35 #include "smalloc.h"
  36 #include "macros.h"
  37 #include "princ.h"
  38 #include "rmpbc.h"
  39 #include "vec.h"
  40 #include "xvgr.h"
  41 #include "pbc.h"
  42 #include "copyrite.h"
  43 #include "futil.h"
  44 #include "statutil.h"
  45 #include "rdgroup.h"
  46
  47 /****************************************************************************/
  48 /* This program calculates the ordering of water molecules across a box, as */
  49 /* function of the z-coordinate. This implies averaging over slices and over*/
  50 /* time. Output is the average cos of the angle of the dipole with the      */
  51 /* normal, per slice.                                                       */
  52 /* In addition, it calculates the average dipole moment itself in three     */
  53 /* directions.                                                              */
  54 /****************************************************************************/
  55
  56 void calc_order(char *fn, atom_id index[], int ngx, rvec **slDipole,
  57                 real **slOrder, real *slWidth, int *nslices,
  58                 t_topology *top, int axis, bool bMicel, atom_id micel[],
  59                 int nmic)
  60 {
  61   rvec *x0,              /* coordinates with pbc */
  62        dipole,           /* dipole moment due to one molecules */
  63        normal,
  64        com;              /* center of mass of micel, with bMicel */
  65   rvec *dip;             /* sum of dipoles, unnormalized */
  66   matrix box;            /* box (3x3) */
  67   int   status;
  68   real  t,               /* time from trajectory */
  69        *sum,             /* sum of all cosines of dipoles, per slice */
  70        *frame;           /* order over one frame */
  71   int natoms,            /* nr. atoms in trj */
  72       i,j,teller = 0,
  73       slice=0,           /* current slice number */
  74       *count;            /* nr. of atoms in one slice */
  75
  76   if ((natoms = read_first_x(&status,fn,&t,&x0,box)) == 0)
  77     fatal_error(0,"Could not read coordinates from statusfile\n");
  78
  79   if (! *nslices)
  80     *nslices = (int)(box[axis][axis] * 10); /* default value */
  81
  82   switch(axis)
  83   {
  84   case 0:
  85     normal[0] = 1; normal[1] = 0; normal[2] = 0;
  86     break;
  87   case 1:
  88     normal[0] = 0; normal[1] = 1; normal[2] = 0;
  89     break;
  90   case 2:
  91     normal[0] = 0; normal[1] = 0; normal[2] = 1;
  92     break;
  93   default:
  94     fatal_error(0,"No valid value for -axis-. Exiting.\n");
  95   }
  96
  97   clear_rvec(dipole);
  98   snew(count, *nslices);
  99   snew(sum, *nslices);
 100   snew(dip, *nslices);
 101   snew(frame, *nslices);
 102
 103   *slWidth = box[axis][axis]/(*nslices);
 104   fprintf(stderr,"Box divided in %d slices. Initial width of slice: %f\n",
 105           *nslices, *slWidth);
 106
 107   teller = 0;
 108
 109   /*********** Start processing trajectory ***********/
 110   do
 111   {
 112     *slWidth = box[axis][axis]/(*nslices);
 113     teller++;
 114
 115     rm_pbc(&(top->idef),top->atoms.nr,box,x0,x0);
 116
 117     if (bMicel)
 118       calc_xcm(x0, nmic, micel, top->atoms.atom, com, FALSE);
 119
 120     for (i = 0; i < ngx/3; i++)
 121     {
 122       /* put all waters in box */
 123       for (j = 0; j < DIM; j++)
 124       {
 125         if (x0[index[3*i]][j] < 0)
 126         {
 127           x0[index[3*i]][j] += box[j][j];
 128           x0[index[3*i+1]][j] += box[j][j];
 129           x0[index[3*i+2]][j] += box[j][j];
 130         }
 131         if (x0[index[3*i]][j] > box[j][j])
 132         {
 133           x0[index[3*i]][j] -= box[j][j];
 134           x0[index[3*i+1]][j] -= box[j][j];
 135           x0[index[3*i+2]][j] -= box[j][j];
 136         }
 137       }
 138
 139       for (j = 0; j < DIM; j++)
 140         dipole[j] =
 141           x0[index[3*i]][j] * top->atoms.atom[index[3*i]].q +
 142           x0[index[3*i+1]][j] * top->atoms.atom[index[3*i+1]].q +
 143           x0[index[3*i+2]][j] * top->atoms.atom[index[3*i+2]].q;
 144
 145       /* now we have a dipole vector. Might as well safe it. Then the
 146          rest depends on whether we're dealing with
 147          a flat or a spherical interface.
 148        */
 149
 150       if (bMicel)
 151       { /* this is for spherical interfaces */
 152         rvec_sub(com, x0[index[3*i]], normal); /* vector from Oxygen to COM */
 153         slice = norm(normal)/(*slWidth);         /* spherical slice           */
 154
 155         sum[slice] += iprod(dipole, normal) / (norm(dipole) * norm(normal));
 156         frame[slice] += iprod(dipole, normal) / (norm(dipole) * norm(normal));
 157         count[slice]++;
 158
 159       } else {
 160         /* this is for flat interfaces      */
 161
 162         /* determine which slice atom is in */
 163         slice = (x0[index[3*i]][axis] / (*slWidth));
 164         if (slice < 0 || slice >= *nslices)
 165         {
 166           fprintf(stderr,"Coordinate: %f ",x0[index[3*i]][axis]);
 167           fprintf(stderr,"HELP PANIC! slice = %d, OUT OF RANGE!\n",slice);
 168         }
 169         else
 170         {
 171           rvec_add(dipole, dip[slice], dip[slice]);
 172           /* Add dipole to total. mag[slice] is total dipole in axis direction */
 173           sum[slice] += iprod(dipole,normal)/norm(dipole);
 174           frame[slice] += iprod(dipole,normal)/norm(dipole);
 175           /* increase count for that slice */
 176           count[slice]++;
 177         }
 178       }
 179     }
 180
 181   } while (read_next_x(status,&t,natoms,x0,box));
 182   /*********** done with status file **********/
 183
 184   fprintf(stderr,"\nRead trajectory. Printing parameters to file\n");
 185
 186   for (i = 0; i < *nslices; i++)  /* average over frames */
 187   {
 188     fprintf(stderr,"%d waters in slice %d\n",count[i],i);
 189     if (count[i] > 0)   /* divide by number of molecules in each slice */
 190       {
 191         sum[i] = sum[i] / count[i];
 192         dip[i][XX] = dip[i][XX] / count[i];
 193         dip[i][YY] = dip[i][YY] / count[i];
 194         dip[i][ZZ] = dip[i][ZZ] / count[i];
 195       }
 196     else
 197       fprintf(stderr,"No water in slice %d\n",i);
 198   }
 199
 200   *slOrder = sum;  /* copy a pointer, I hope */
 201   *slDipole = dip;
 202   sfree(x0);  /* free memory used by coordinate arrays */
 203 }
 204
 205 void order_plot(rvec dipole[], real order[], char *afile,
 206                 int nslices, real slWidth)
 207 {
 208   FILE       *ord;                /* xvgr files with order parameters  */
 209   int        slice;               /* loop index     */
 210   char       buf[256];            /* for xvgr title */
 211   real       factor;              /* conversion to Debye from electron*nm */
 212
 213   /*  factor = 1e-9*1.60217733e-19/3.336e-30 */
 214   factor = 1.60217733/3.336e-2;
 215   fprintf(stderr,"%d slices\n",nslices);
 216   sprintf(buf,"Water orientation with respect to normal");
 217   ord = xvgropen(afile,buf,"box (nm)","mu_x, mu_y, mu_z (D), cosine with normal");
 218
 219   for (slice = 0; slice < nslices; slice++)
 220     fprintf(ord,"%8.3f %8.3f %8.3f %8.3f %e\n", slWidth*slice,
 221             factor*dipole[slice][XX], factor*dipole[slice][YY],
 222             factor*dipole[slice][ZZ], order[slice]);
 223
 224   fclose(ord);
 225 }
 226
 227 int main(int argc,char *argv[])
 228 {
 229   static char *desc[] = {
 230     "Compute the orientation of water molecules with respect to the normal",
 231     "of the box. The program determines the average cosine of the angle",
 232     "between de dipole moment of water and an axis of the box. The box is",
 233     "divided in slices and the average orientation per slice is printed.",
 234     "Each water molecule is assigned to a slice, per time frame, based on the",
 235     "position of the oxygen. When -nm  is used the angle between the water",
 236     "dipole and the axis from the center of mass to the oxygen is calculated",
 237     "instead of the angle between the dipole and a box axis."
 238   };
 239   static int  axis = 2;                       /* normal to memb. default z  */
 240   static char *axtitle="Z";
 241   static int  nslices = 0;                    /* nr of slices defined       */
 242   t_pargs pa[] = {
 243     { "-d",   FALSE, etSTR, {&axtitle},
 244       "Take the normal on the membrane in direction X, Y or Z." },
 245     { "-sl",  FALSE, etINT, {&nslices},
 246       "Calculate order parameter as function of boxlength, dividing the box"
 247       " in #nr slices."}
 248   };
 249   static char *bugs[] = {
 250     "The program assigns whole water molecules to a slice, based on the first"
 251     "atom of three in the index file group. It assumes an order O,H,H."
 252     "Name is not important, but the order is. If this demand is not met,"
 253     "assigning molecules to slices is different."
 254   };
 255
 256   real      *slOrder,                       /* av. cosine, per slice      */
 257             slWidth = 0.0;                  /* width of a slice           */
 258   rvec      *slDipole;
 259   char      *grpname,                       /* groupnames                 */
 260             *micname;
 261   int       ngx,                            /* nr. of atomsin sol group   */
 262             nmic;                           /* nr. of atoms in micelle    */
 263   t_topology *top;                          /* topology                   */
 264   atom_id    *index,                        /* indices for solvent group  */
 265              *micelle;
 266   bool       bMicel =  FALSE;               /* think we're a micel        */
 267   t_filenm  fnm[] = {                       /* files for g_order          */
 268     { efTRX, "-f", NULL,  ffREAD },         /* trajectory file            */
 269     { efNDX, NULL, NULL,  ffREAD },         /* index file                 */
 270     { efNDX, "-nm", NULL, ffOPTRD },        /* index with micelle atoms   */
 271     { efTPX, NULL, NULL,  ffREAD },         /* topology file              */
 272     { efXVG,"-o",  "order", ffWRITE },      /* xvgr output file           */
 273   };
 274
 275 #define NFILE asize(fnm)
 276
 277   CopyRight(stderr,argv[0]);
 278   parse_common_args(&argc, argv, PCA_CAN_VIEW | PCA_CAN_TIME, TRUE, NFILE,
 279                     fnm, asize(pa),pa,asize(desc),desc,asize(bugs),bugs);
 280   bMicel = opt2bSet("-nm",NFILE,fnm);
 281
 282   top = read_top(ftp2fn(efTPX,NFILE,fnm));     /* read topology file */
 283
 284   rd_index(ftp2fn(efNDX,NFILE,fnm),1,&ngx,&index,&grpname);
 285
 286   if (bMicel)
 287     rd_index(opt2fn("-nm",NFILE,fnm), 1, &nmic, &micelle, &micname);
 288
 289   calc_order(ftp2fn(efTRX,NFILE,fnm), index, ngx, &slDipole, &slOrder,
 290              &slWidth, &nslices, top, axis, bMicel, micelle, nmic);
 291
 292   order_plot(slDipole, slOrder, opt2fn("-o",NFILE,fnm), nslices,
 293              slWidth);
 294
 295   xvgr_file(opt2fn("-o",NFILE,fnm), NULL);      /* view xvgr file */
 296   thanx(stdout);
 297
 298   return 0;
 299 }
 300
 301
 302
 303
 304
 305
 306
 307
 308