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