src/tools/gmx_potential.c

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