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       for (i = 0; i < gnx[n]; i++)   /* loop over all atoms in index file */
 173       {
 174         if (bSpherical)
 175         {
 176           rvec_add(x0[index[n][i]], xcm, x0[index[n][i]]);
 177           /* only distance from origin counts, not sign */
 178           slice = norm(x0[index[n][i]])/(*slWidth);
 179
 180           /* this is a nice check for spherical groups but not for
 181              all water in a cubic box since a lot will fall outside
 182              the sphere
 183             if (slice > (*nslices))
 184             {
 185              fprintf(stderr,"Warning: slice = %d\n",slice);
 186             }
 187           */
 188           (*slCharge)[n][slice] += top->atoms.atom[index[n][i]].q;
 189         }
 190         else
 191         {
 192           z = x0[index[n][i]][axis];
 193           z = z + fudge_z;
 194           if (z < 0)
 195             z += box[axis][axis];
 196           if (z > box[axis][axis])
 197             z -= box[axis][axis];
 198           /* determine which slice atom is in */
 199           slice = (z / (*slWidth));
 200           (*slCharge)[n][slice] += top->atoms.atom[index[n][i]].q;
 201         }
 202       }
 203     }
 204     nr_frames++;
 205   } while (read_next_x(oenv,status,&t,natoms,x0,box));
 206
 207   gmx_rmpbc_done(gpbc);
 208
 209   /*********** done with status file **********/
 210   close_trj(status);
 211
 212   /* slCharge now contains the total charge per slice, summed over all
 213      frames. Now divide by nr_frames and integrate twice
 214    */
 215
 216
 217   if (bSpherical)
 218     fprintf(stderr,"\n\nRead %d frames from trajectory. Calculating potential"
 219             "in spherical coordinates\n", nr_frames);
 220   else
 221     fprintf(stderr,"\n\nRead %d frames from trajectory. Calculating potential\n",
 222             nr_frames);
 223
 224   for (n =0; n < nr_grps; n++)
 225   {
 226       for (i = 0; i < *nslices; i++)
 227       {
 228           if (bSpherical)
 229           {
 230               /* charge per volume is now the summed charge, divided by the nr
 231               of frames and by the volume of the slice it's in, 4pi r^2 dr
 232               */
 233               slVolume = 4*M_PI * sqr(i) * sqr(*slWidth) * *slWidth;
 234               if (slVolume == 0)
 235               {
 236                   (*slCharge)[n][i] = 0;
 237               }
 238               else
 239               {
 240                   (*slCharge)[n][i] = (*slCharge)[n][i] / (nr_frames * slVolume);
 241               }
 242           }
 243           else
 244           {
 245               /* get charge per volume */
 246               (*slCharge)[n][i] = (*slCharge)[n][i] * (*nslices) /
 247               (nr_frames * box[axis][axis] * box[ax1][ax1] * box[ax2][ax2]);
 248           }
 249       }
 250       /* Now we have charge densities */
 251   }
 252
 253   if(bCorrect && !bSpherical)
 254   {
 255       for(n =0; n < nr_grps; n++)
 256       {
 257           nn = 0;
 258           qsum = 0;
 259           for (i = 0; i < *nslices; i++)
 260           {
 261               if( fabs((*slCharge)[n][i]) >= GMX_DOUBLE_MIN )
 262               {
 263                   nn++;
 264                   qsum += (*slCharge)[n][i];
 265               }
 266           }
 267           qsum /= nn;
 268           for (i = 0; i < *nslices; i++)
 269           {
 270               if( fabs((*slCharge)[n][i]) >= GMX_DOUBLE_MIN )
 271               {
 272                   (*slCharge)[n][i] -= qsum;
 273               }
 274           }
 275       }
 276   }
 277
 278   for(n =0; n < nr_grps; n++)
 279   {
 280       /* integrate twice to get field and potential */
 281       p_integrate((*slField)[n], (*slCharge)[n], *nslices, *slWidth);
 282   }
 283
 284
 285   if(bCorrect && !bSpherical)
 286   {
 287       for(n =0; n < nr_grps; n++)
 288       {
 289           nn = 0;
 290           qsum = 0;
 291           for (i = 0; i < *nslices; i++)
 292           {
 293               if( fabs((*slCharge)[n][i]) >= GMX_DOUBLE_MIN )
 294               {
 295                   nn++;
 296                   qsum += (*slField)[n][i];
 297               }
 298           }
 299           qsum /= nn;
 300           for (i = 0; i < *nslices; i++)
 301           {
 302               if( fabs((*slCharge)[n][i]) >= GMX_DOUBLE_MIN )
 303               {
 304                   (*slField)[n][i] -= qsum;
 305               }
 306           }
 307       }
 308   }
 309
 310   for(n =0; n < nr_grps; n++)
 311   {
 312       p_integrate((*slPotential)[n],(*slField)[n], *nslices, *slWidth);
 313   }
 314
 315   /* Now correct for eps0 and in spherical case for r*/
 316   for (n = 0; n < nr_grps; n++)
 317     for (i = 0; i < *nslices; i++)
 318     {
 319       if (bSpherical)
 320       {
 321         (*slPotential)[n][i] = ELC * (*slPotential)[n][i] * -1.0E9 /
 322           (EPS0 * i * (*slWidth));
 323         (*slField)[n][i] = ELC * (*slField)[n][i] * 1E18 /
 324           (EPS0 * i * (*slWidth));
 325       }
 326       else
 327       {
 328         (*slPotential)[n][i] = ELC * (*slPotential)[n][i] * -1.0E9 / EPS0  ;
 329         (*slField)[n][i] = ELC * (*slField)[n][i] * 1E18 / EPS0;
 330       }
 331     }
 332
 333   sfree(x0);  /* free memory used by coordinate array */
 334 }
 335
 336 void plot_potential(double *potential[], double *charge[], double *field[],
 337                     const char *afile, const char *bfile, const char *cfile,
 338                     int nslices, int nr_grps, const char *grpname[], double slWidth,
 339                     const output_env_t oenv)
 340 {
 341   FILE       *pot,     /* xvgr file with potential */
 342              *cha,     /* xvgr file with charges   */
 343              *fie;     /* xvgr files with fields   */
 344   char       buf[256]; /* for xvgr title */
 345   int        slice, n;
 346
 347   sprintf(buf,"Electrostatic Potential");
 348   pot = xvgropen(afile, buf, "Box (nm)","Potential (V)",oenv);
 349   xvgr_legend(pot,nr_grps,grpname,oenv);
 350
 351   sprintf(buf,"Charge Distribution");
 352   cha = xvgropen(bfile, buf, "Box (nm)", "Charge density (q/nm\\S3\\N)",oenv);
 353   xvgr_legend(cha,nr_grps,grpname,oenv);
 354
 355   sprintf(buf, "Electric Field");
 356   fie = xvgropen(cfile, buf, "Box (nm)", "Field (V/nm)",oenv);
 357   xvgr_legend(fie,nr_grps,grpname,oenv);
 358
 359   for (slice = cb; slice < (nslices - ce); slice++)
 360   {
 361     fprintf(pot,"%20.16g  ", slice * slWidth);
 362     fprintf(cha,"%20.16g  ", slice * slWidth);
 363     fprintf(fie,"%20.16g  ", slice * slWidth);
 364     for (n = 0; n < nr_grps; n++)
 365     {
 366       fprintf(pot,"   %20.16g", potential[n][slice]);
 367       fprintf(fie,"   %20.16g", field[n][slice]);
 368       fprintf(cha,"   %20.16g", charge[n][slice]);
 369     }
 370     fprintf(pot,"\n");
 371     fprintf(cha,"\n");
 372     fprintf(fie,"\n");
 373   }
 374
 375   ffclose(pot);
 376   ffclose(cha);
 377   ffclose(fie);
 378 }
 379
 380 int gmx_potential(int argc,char *argv[])
 381 {
 382   const char *desc[] = {
 383     "Compute the electrostatical potential across the box. The potential is",
 384     "calculated by first summing the charges per slice and then integrating",
 385     "twice of this charge distribution. Periodic boundaries are not taken",
 386     "into account. Reference of potential is taken to be the left side of",
 387     "the box. It's also possible to calculate the potential in spherical",
 388     "coordinates as function of r by calculating a charge distribution in",
 389     "spherical slices and twice integrating them. epsilon_r is taken as 1,",
 390     "2 is more appropriate in many cases."
 391   };
 392   output_env_t oenv;
 393   static int  axis = 2;                      /* normal to memb. default z  */
 394   static const char *axtitle="Z";
 395   static int  nslices = 10;                  /* nr of slices defined       */
 396   static int  ngrps   = 1;
 397   static gmx_bool bSpherical = FALSE;            /* default is bilayer types   */
 398   static real fudge_z = 0;                    /* translate coordinates      */
 399   static gmx_bool bCorrect = 0;
 400   t_pargs pa [] = {
 401     { "-d",   FALSE, etSTR, {&axtitle},
 402       "Take the normal on the membrane in direction X, Y or Z." },
 403     { "-sl",  FALSE, etINT, {&nslices},
 404       "Calculate potential as function of boxlength, dividing the box"
 405       " in #nr slices." } ,
 406     { "-cb",  FALSE, etINT, {&cb},
 407       "Discard first #nr slices of box for integration" },
 408     { "-ce",  FALSE, etINT, {&ce},
 409       "Discard last #nr slices of box for integration" },
 410     { "-tz",  FALSE, etREAL, {&fudge_z},
 411       "Translate all coordinates <distance> in the direction of the box" },
 412     { "-spherical", FALSE, etBOOL, {&bSpherical},
 413       "Calculate spherical thingie" },
 414     { "-ng",       FALSE, etINT, {&ngrps},
 415       "Number of groups to consider" },
 416     { "-correct",  FALSE, etBOOL, {&bCorrect},
 417       "Assume net zero charge of groups to improve accuracy" }
 418   };
 419   const char *bugs[] = {
 420     "Discarding slices for integration should not be necessary."
 421   };
 422
 423   double    **potential,                    /* potential per slice        */
 424             **charge,                       /* total charge per slice     */
 425             **field,                        /* field per slice            */
 426             slWidth;                        /* width of one slice         */
 427   char      **grpname;                      /* groupnames                 */
 428   int       *ngx;                           /* sizes of groups            */
 429   t_topology *top;                          /* topology                   */
 430   int       ePBC;
 431   atom_id   **index;                        /* indices for all groups     */
 432   t_filenm  fnm[] = {                       /* files for g_order          */
 433     { efTRX, "-f", NULL,  ffREAD },         /* trajectory file            */
 434     { efNDX, NULL, NULL,  ffREAD },         /* index file                 */
 435     { efTPX, NULL, NULL,  ffREAD },         /* topology file              */
 436     { efXVG,"-o","potential", ffWRITE },    /* xvgr output file           */
 437     { efXVG,"-oc","charge", ffWRITE },      /* xvgr output file           */
 438     { efXVG,"-of","field", ffWRITE },       /* xvgr output file           */
 439   };
 440
 441 #define NFILE asize(fnm)
 442
 443   CopyRight(stderr,argv[0]);
 444   parse_common_args(&argc, argv, PCA_CAN_VIEW | PCA_CAN_TIME | PCA_BE_NICE,
 445                     NFILE,fnm,asize(pa),pa,asize(desc),desc,asize(bugs),bugs,
 446                     &oenv);
 447
 448   /* Calculate axis */
 449   axis = toupper(axtitle[0]) - 'X';
 450
 451   top = read_top(ftp2fn(efTPX,NFILE,fnm),&ePBC);     /* read topology file */
 452
 453   snew(grpname,ngrps);
 454   snew(index,ngrps);
 455   snew(ngx,ngrps);
 456
 457   rd_index(ftp2fn(efNDX,NFILE,fnm),ngrps,ngx,index,grpname);
 458
 459
 460   calc_potential(ftp2fn(efTRX,NFILE,fnm), index, ngx,
 461                  &potential, &charge, &field,
 462                  &nslices, top, ePBC, axis, ngrps, &slWidth, fudge_z,
 463                  bSpherical, bCorrect,oenv);
 464
 465   plot_potential(potential, charge, field, opt2fn("-o",NFILE,fnm),
 466                  opt2fn("-oc",NFILE,fnm), opt2fn("-of",NFILE,fnm),
 467                  nslices, ngrps, (const char**)grpname, slWidth,oenv);
 468
 469   do_view(oenv,opt2fn("-o",NFILE,fnm), NULL);       /* view xvgr file */
 470   do_view(oenv,opt2fn("-oc",NFILE,fnm), NULL);      /* view xvgr file */
 471   do_view(oenv,opt2fn("-of",NFILE,fnm), NULL);      /* view xvgr file */
 472
 473   thanx(stderr);
 474   return 0;
 475 }
 476
 477
 478
 479
 480
 481
 482
 483
 484