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