src/tools/gmx_potential.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 <ctype.h>
  41
  42 #include "sysstuff.h"
  43 #include "string.h"
  44 #include "typedefs.h"
  45 #include "smalloc.h"
  46 #include "macros.h"
  47 #include "princ.h"
  48 #include "rmpbc.h"
  49 #include "vec.h"
  50 #include "xvgr.h"
  51 #include "pbc.h"
  52 #include "copyrite.h"
  53 #include "futil.h"
  54 #include "statutil.h"
  55 #include "index.h"
  56 #include "gmx_ana.h"
  57
  58
  59 #define EPS0 8.85419E-12
  60 #define ELC 1.60219E-19
  61
  62 /****************************************************************************/
  63 /* This program calculates the electrostatic potential across the box by    */
  64 /* determining the charge density in slices of the box and integrating these*/
  65 /* twice.                                                                   */
  66 /* Peter Tieleman, April 1995                                               */
  67 /* It now also calculates electrostatic potential in spherical micelles,    */
  68 /* using \frac{1}{r}\frac{d^2r\Psi}{r^2} = - \frac{\rho}{\epsilon_0}        */
  69 /* This probably sucks but it seems to work.                                */
  70 /****************************************************************************/
  71
  72 static int ce=0, cb=0;
  73
  74 /* this routine integrates the array data and returns the resulting array */
  75 /* routine uses simple trapezoid rule                                     */
  76 void p_integrate(double *result, double data[], int ndata, double slWidth)
  77 {
  78   int i, slice;
  79   double sum;
  80
  81   if (ndata <= 2)
  82     fprintf(stderr,"Warning: nr of slices very small. This will result"
  83             "in nonsense.\n");
  84
  85   fprintf(stderr,"Integrating from slice %d to slice %d\n",cb, ndata-ce);
  86
  87   for (slice = cb; slice < (ndata-ce); slice ++)
  88   {
  89     sum = 0;
  90     for (i = cb; i < slice; i++)
  91       sum += slWidth * (data[i] + 0.5 * (data[i+1] - data[i]));
  92     result[slice] = sum;
  93   }
  94   return;
  95 }
  96
  97 void calc_potential(const char *fn, atom_id **index, int gnx[],
  98                     double ***slPotential, double ***slCharge,
  99                     double ***slField, int *nslices,
 100                     t_topology *top, int ePBC,
 101                     int axis, int nr_grps, double *slWidth,
 102                     double fudge_z, gmx_bool bSpherical, gmx_bool bCorrect,
 103                     const output_env_t oenv)
 104 {
 105   rvec *x0;              /* coordinates without pbc */
 106   matrix box;            /* box (3x3) */
 107   int natoms;            /* nr. atoms in trj */
 108   t_trxstatus *status;
 109   int **slCount,         /* nr. of atoms in one slice for a group */
 110       i,j,n,             /* loop indices */
 111       teller = 0,
 112       ax1=0, ax2=0,
 113       nr_frames = 0,     /* number of frames */
 114       slice;             /* current slice */
 115   double slVolume;         /* volume of slice for spherical averaging */
 116   double qsum,nn;
 117   real t;
 118   double z;
 119   rvec xcm;
 120   gmx_rmpbc_t  gpbc=NULL;
 121
 122   switch(axis)
 123   {
 124   case 0:
 125     ax1 = 1; ax2 = 2;
 126     break;
 127   case 1:
 128     ax1 = 0; ax2 = 2;
 129     break;
 130   case 2:
 131     ax1 = 0; ax2 = 1;
 132     break;
 133   default:
 134     gmx_fatal(FARGS,"Invalid axes. Terminating\n");
 135   }
 136
 137   if ((natoms = read_first_x(oenv,&status,fn,&t,&x0,box)) == 0)
 138     gmx_fatal(FARGS,"Could not read coordinates from statusfile\n");
 139
 140   if (! *nslices)
 141     *nslices = (int)(box[axis][axis] * 10); /* default value */
 142
 143   fprintf(stderr,"\nDividing the box in %d slices\n",*nslices);
 144
 145   snew(*slField, nr_grps);
 146   snew(*slCharge, nr_grps);
 147   snew(*slPotential, nr_grps);
 148
 149   for (i = 0; i < nr_grps; i++)
 150   {
 151     snew((*slField)[i], *nslices);
 152     snew((*slCharge)[i], *nslices);
 153     snew((*slPotential)[i], *nslices);
 154   }
 155
 156
 157   gpbc = gmx_rmpbc_init(&top->idef,ePBC,natoms,box);
 158
 159   /*********** Start processing trajectory ***********/
 160   do
 161   {
 162     *slWidth = box[axis][axis]/(*nslices);
 163     teller++;
 164     gmx_rmpbc(gpbc,natoms,box,x0);
 165
 166     /* calculate position of center of mass based on group 1 */
 167     calc_xcm(x0, gnx[0], index[0], top->atoms.atom, xcm, FALSE);
 168     svmul(-1,xcm,xcm);
 169
 170     for (n = 0; n < nr_grps; n++)
 171     {
 172         /* Check whether we actually have all positions of the requested index
 173          * group in the trajectory file */
 174         if (gnx[n] > natoms)
 175         {
 176             gmx_fatal(FARGS, "You selected a group with %d atoms, but only %d atoms\n"
 177                              "were found in the trajectory.\n", gnx[n], natoms);
 178         }
 179       for (i = 0; i < gnx[n]; i++)   /* loop over all atoms in index file */
 180       {
 181         if (bSpherical)
 182         {
 183           rvec_add(x0[index[n][i]], xcm, x0[index[n][i]]);
 184           /* only distance from origin counts, not sign */
 185           slice = norm(x0[index[n][i]])/(*slWidth);
 186
 187           /* this is a nice check for spherical groups but not for
 188              all water in a cubic box since a lot will fall outside
 189              the sphere
 190             if (slice > (*nslices))
 191             {
 192              fprintf(stderr,"Warning: slice = %d\n",slice);
 193             }
 194           */
 195           (*slCharge)[n][slice] += top->atoms.atom[index[n][i]].q;
 196         }
 197         else
 198         {
 199           z = x0[index[n][i]][axis];
 200           z = z + fudge_z;
 201           if (z < 0)
 202             z += box[axis][axis];
 203           if (z > box[axis][axis])
 204             z -= box[axis][axis];
 205           /* determine which slice atom is in */
 206           slice = (z / (*slWidth));
 207           (*slCharge)[n][slice] += top->atoms.atom[index[n][i]].q;
 208         }
 209       }
 210     }
 211     nr_frames++;
 212   } while (read_next_x(oenv,status,&t,natoms,x0,box));
 213
 214   gmx_rmpbc_done(gpbc);
 215
 216   /*********** done with status file **********/
 217   close_trj(status);
 218
 219   /* slCharge now contains the total charge per slice, summed over all
 220      frames. Now divide by nr_frames and integrate twice
 221    */
 222
 223
 224   if (bSpherical)
 225     fprintf(stderr,"\n\nRead %d frames from trajectory. Calculating potential"
 226             "in spherical coordinates\n", nr_frames);
 227   else
 228     fprintf(stderr,"\n\nRead %d frames from trajectory. Calculating potential\n",
 229             nr_frames);
 230
 231   for (n =0; n < nr_grps; n++)
 232   {
 233       for (i = 0; i < *nslices; i++)
 234       {
 235           if (bSpherical)
 236           {
 237               /* charge per volume is now the summed charge, divided by the nr
 238               of frames and by the volume of the slice it's in, 4pi r^2 dr
 239               */
 240               slVolume = 4*M_PI * sqr(i) * sqr(*slWidth) * *slWidth;
 241               if (slVolume == 0)
 242               {
 243                   (*slCharge)[n][i] = 0;
 244               }
 245               else
 246               {
 247                   (*slCharge)[n][i] = (*slCharge)[n][i] / (nr_frames * slVolume);
 248               }
 249           }
 250           else
 251           {
 252               /* get charge per volume */
 253               (*slCharge)[n][i] = (*slCharge)[n][i] * (*nslices) /
 254               (nr_frames * box[axis][axis] * box[ax1][ax1] * box[ax2][ax2]);
 255           }
 256       }
 257       /* Now we have charge densities */
 258   }
 259
 260   if(bCorrect && !bSpherical)
 261   {
 262       for(n =0; n < nr_grps; n++)
 263       {
 264           nn = 0;
 265           qsum = 0;
 266           for (i = 0; i < *nslices; i++)
 267           {
 268               if( fabs((*slCharge)[n][i]) >= GMX_DOUBLE_MIN )
 269               {
 270                   nn++;
 271                   qsum += (*slCharge)[n][i];
 272               }
 273           }
 274           qsum /= nn;
 275           for (i = 0; i < *nslices; i++)
 276           {
 277               if( fabs((*slCharge)[n][i]) >= GMX_DOUBLE_MIN )
 278               {
 279                   (*slCharge)[n][i] -= qsum;
 280               }
 281           }
 282       }
 283   }
 284
 285   for(n =0; n < nr_grps; n++)
 286   {
 287       /* integrate twice to get field and potential */
 288       p_integrate((*slField)[n], (*slCharge)[n], *nslices, *slWidth);
 289   }
 290
 291
 292   if(bCorrect && !bSpherical)
 293   {
 294       for(n =0; n < nr_grps; n++)
 295       {
 296           nn = 0;
 297           qsum = 0;
 298           for (i = 0; i < *nslices; i++)
 299           {
 300               if( fabs((*slCharge)[n][i]) >= GMX_DOUBLE_MIN )
 301               {
 302                   nn++;
 303                   qsum += (*slField)[n][i];
 304               }
 305           }
 306           qsum /= nn;
 307           for (i = 0; i < *nslices; i++)
 308           {
 309               if( fabs((*slCharge)[n][i]) >= GMX_DOUBLE_MIN )
 310               {
 311                   (*slField)[n][i] -= qsum;
 312               }
 313           }
 314       }
 315   }
 316
 317   for(n =0; n < nr_grps; n++)
 318   {
 319       p_integrate((*slPotential)[n],(*slField)[n], *nslices, *slWidth);
 320   }
 321
 322   /* Now correct for eps0 and in spherical case for r*/
 323   for (n = 0; n < nr_grps; n++)
 324     for (i = 0; i < *nslices; i++)
 325     {
 326       if (bSpherical)
 327       {
 328         (*slPotential)[n][i] = ELC * (*slPotential)[n][i] * -1.0E9 /
 329           (EPS0 * i * (*slWidth));
 330         (*slField)[n][i] = ELC * (*slField)[n][i] * 1E18 /
 331           (EPS0 * i * (*slWidth));
 332       }
 333       else
 334       {
 335         (*slPotential)[n][i] = ELC * (*slPotential)[n][i] * -1.0E9 / EPS0  ;
 336         (*slField)[n][i] = ELC * (*slField)[n][i] * 1E18 / EPS0;
 337       }
 338     }
 339
 340   sfree(x0);  /* free memory used by coordinate array */
 341 }
 342
 343 void plot_potential(double *potential[], double *charge[], double *field[],
 344                     const char *afile, const char *bfile, const char *cfile,
 345                     int nslices, int nr_grps, const char *grpname[], double slWidth,
 346                     const output_env_t oenv)
 347 {
 348   FILE       *pot,     /* xvgr file with potential */
 349              *cha,     /* xvgr file with charges   */
 350              *fie;     /* xvgr files with fields   */
 351   char       buf[256]; /* for xvgr title */
 352   int        slice, n;
 353
 354   sprintf(buf,"Electrostatic Potential");
 355   pot = xvgropen(afile, buf, "Box (nm)","Potential (V)",oenv);
 356   xvgr_legend(pot,nr_grps,grpname,oenv);
 357
 358   sprintf(buf,"Charge Distribution");
 359   cha = xvgropen(bfile, buf, "Box (nm)", "Charge density (q/nm\\S3\\N)",oenv);
 360   xvgr_legend(cha,nr_grps,grpname,oenv);
 361
 362   sprintf(buf, "Electric Field");
 363   fie = xvgropen(cfile, buf, "Box (nm)", "Field (V/nm)",oenv);
 364   xvgr_legend(fie,nr_grps,grpname,oenv);
 365
 366   for (slice = cb; slice < (nslices - ce); slice++)
 367   {
 368     fprintf(pot,"%20.16g  ", slice * slWidth);
 369     fprintf(cha,"%20.16g  ", slice * slWidth);
 370     fprintf(fie,"%20.16g  ", slice * slWidth);
 371     for (n = 0; n < nr_grps; n++)
 372     {
 373       fprintf(pot,"   %20.16g", potential[n][slice]);
 374       fprintf(fie,"   %20.16g", field[n][slice]/1e9);  /* convert to V/nm */
 375       fprintf(cha,"   %20.16g", charge[n][slice]);
 376     }
 377     fprintf(pot,"\n");
 378     fprintf(cha,"\n");
 379     fprintf(fie,"\n");
 380   }
 381
 382   ffclose(pot);
 383   ffclose(cha);
 384   ffclose(fie);
 385 }
 386
 387 int gmx_potential(int argc,char *argv[])
 388 {
 389   const char *desc[] = {
 390     "[TT]g_potential[tt] computes the electrostatical potential across the box. The potential is",
 391     "calculated by first summing the charges per slice and then integrating",
 392     "twice of this charge distribution. Periodic boundaries are not taken",
 393     "into account. Reference of potential is taken to be the left side of",
 394     "the box. It is also possible to calculate the potential in spherical",
 395     "coordinates as function of r by calculating a charge distribution in",
 396     "spherical slices and twice integrating them. epsilon_r is taken as 1,",
 397     "but 2 is more appropriate in many cases."
 398   };
 399   output_env_t oenv;
 400   static int  axis = 2;                      /* normal to memb. default z  */
 401   static const char *axtitle="Z";
 402   static int  nslices = 10;                  /* nr of slices defined       */
 403   static int  ngrps   = 1;
 404   static gmx_bool bSpherical = FALSE;            /* default is bilayer types   */
 405   static real fudge_z = 0;                    /* translate coordinates      */
 406   static gmx_bool bCorrect = 0;
 407   t_pargs pa [] = {
 408     { "-d",   FALSE, etSTR, {&axtitle},
 409       "Take the normal on the membrane in direction X, Y or Z." },
 410     { "-sl",  FALSE, etINT, {&nslices},
 411       "Calculate potential as function of boxlength, dividing the box"
 412       " in #nr slices." } ,
 413     { "-cb",  FALSE, etINT, {&cb},
 414       "Discard first #nr slices of box for integration" },
 415     { "-ce",  FALSE, etINT, {&ce},
 416       "Discard last #nr slices of box for integration" },
 417     { "-tz",  FALSE, etREAL, {&fudge_z},
 418       "Translate all coordinates <distance> in the direction of the box" },
 419     { "-spherical", FALSE, etBOOL, {&bSpherical},
 420       "Calculate spherical thingie" },
 421     { "-ng",       FALSE, etINT, {&ngrps},
 422       "Number of groups to consider" },
 423     { "-correct",  FALSE, etBOOL, {&bCorrect},
 424       "Assume net zero charge of groups to improve accuracy" }
 425   };
 426   const char *bugs[] = {
 427     "Discarding slices for integration should not be necessary."
 428   };
 429
 430   double    **potential,                    /* potential per slice        */
 431             **charge,                       /* total charge per slice     */
 432             **field,                        /* field per slice            */
 433             slWidth;                        /* width of one slice         */
 434   char      **grpname;                      /* groupnames                 */
 435   int       *ngx;                           /* sizes of groups            */
 436   t_topology *top;                          /* topology                   */
 437   int       ePBC;
 438   atom_id   **index;                        /* indices for all groups     */
 439   t_filenm  fnm[] = {                       /* files for g_order          */
 440     { efTRX, "-f", NULL,  ffREAD },         /* trajectory file            */
 441     { efNDX, NULL, NULL,  ffREAD },         /* index file                 */
 442     { efTPX, NULL, NULL,  ffREAD },         /* topology file              */
 443     { efXVG,"-o","potential", ffWRITE },    /* xvgr output file           */
 444     { efXVG,"-oc","charge", ffWRITE },      /* xvgr output file           */
 445     { efXVG,"-of","field", ffWRITE },       /* xvgr output file           */
 446   };
 447
 448 #define NFILE asize(fnm)
 449
 450   CopyRight(stderr,argv[0]);
 451   parse_common_args(&argc, argv, PCA_CAN_VIEW | PCA_CAN_TIME | PCA_BE_NICE,
 452                     NFILE,fnm,asize(pa),pa,asize(desc),desc,asize(bugs),bugs,
 453                     &oenv);
 454
 455   /* Calculate axis */
 456   axis = toupper(axtitle[0]) - 'X';
 457
 458   top = read_top(ftp2fn(efTPX,NFILE,fnm),&ePBC);     /* read topology file */
 459
 460   snew(grpname,ngrps);
 461   snew(index,ngrps);
 462   snew(ngx,ngrps);
 463
 464   rd_index(ftp2fn(efNDX,NFILE,fnm),ngrps,ngx,index,grpname);
 465
 466
 467   calc_potential(ftp2fn(efTRX,NFILE,fnm), index, ngx,
 468                  &potential, &charge, &field,
 469                  &nslices, top, ePBC, axis, ngrps, &slWidth, fudge_z,
 470                  bSpherical, bCorrect,oenv);
 471
 472   plot_potential(potential, charge, field, opt2fn("-o",NFILE,fnm),
 473                  opt2fn("-oc",NFILE,fnm), opt2fn("-of",NFILE,fnm),
 474                  nslices, ngrps, (const char**)grpname, slWidth,oenv);
 475
 476   do_view(oenv,opt2fn("-o",NFILE,fnm), NULL);       /* view xvgr file */
 477   do_view(oenv,opt2fn("-oc",NFILE,fnm), NULL);      /* view xvgr file */
 478   do_view(oenv,opt2fn("-of",NFILE,fnm), NULL);      /* view xvgr file */
 479
 480   thanx(stderr);
 481   return 0;
 482 }
 483
 484
 485
 486
 487
 488
 489
 490
 491