3 * This source code is part of
7 * GROningen MAchine for Chemical Simulations
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.
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.
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.
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.
30 * For more info, check our website at http://www.gromacs.org
33 * Green Red Orange Magenta Azure Cyan Skyblue
59 #define EPS0 8.85419E-12
60 #define ELC 1.60219E-19
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*/
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 /****************************************************************************/
72 static int ce=0, cb=0;
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)
82 fprintf(stderr,"Warning: nr of slices very small. This will result"
85 fprintf(stderr,"Integrating from slice %d to slice %d\n",cb, ndata-ce);
87 for (slice = cb; slice < (ndata-ce); slice ++)
90 for (i = cb; i < slice; i++)
91 sum += slWidth * (data[i] + 0.5 * (data[i+1] - data[i]));
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, bool bSpherical, bool bCorrect,
103 const output_env_t oenv)
105 rvec *x0; /* coordinates without pbc */
106 matrix box; /* box (3x3) */
107 int natoms; /* nr. atoms in trj */
109 int **slCount, /* nr. of atoms in one slice for a group */
110 i,j,n, /* loop indices */
113 nr_frames = 0, /* number of frames */
114 slice; /* current slice */
115 double slVolume; /* volume of slice for spherical averaging */
120 gmx_rmpbc_t gpbc=NULL;
134 gmx_fatal(FARGS,"Invalid axes. Terminating\n");
137 if ((natoms = read_first_x(oenv,&status,fn,&t,&x0,box)) == 0)
138 gmx_fatal(FARGS,"Could not read coordinates from statusfile\n");
141 *nslices = (int)(box[axis][axis] * 10); /* default value */
143 fprintf(stderr,"\nDividing the box in %d slices\n",*nslices);
145 snew(*slField, nr_grps);
146 snew(*slCharge, nr_grps);
147 snew(*slPotential, nr_grps);
149 for (i = 0; i < nr_grps; i++)
151 snew((*slField)[i], *nslices);
152 snew((*slCharge)[i], *nslices);
153 snew((*slPotential)[i], *nslices);
157 gpbc = gmx_rmpbc_init(&top->idef,ePBC,natoms,box);
159 /*********** Start processing trajectory ***********/
162 *slWidth = box[axis][axis]/(*nslices);
164 gmx_rmpbc(gpbc,natoms,box,x0);
166 /* calculate position of center of mass based on group 1 */
167 calc_xcm(x0, gnx[0], index[0], top->atoms.atom, xcm, FALSE);
170 for (n = 0; n < nr_grps; n++)
172 for (i = 0; i < gnx[n]; i++) /* loop over all atoms in index file */
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);
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
183 if (slice > (*nslices))
185 fprintf(stderr,"Warning: slice = %d\n",slice);
188 (*slCharge)[n][slice] += top->atoms.atom[index[n][i]].q;
192 z = x0[index[n][i]][axis];
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;
205 } while (read_next_x(oenv,status,&t,natoms,x0,box));
207 gmx_rmpbc_done(gpbc);
209 /*********** done with status file **********/
212 /* slCharge now contains the total charge per slice, summed over all
213 frames. Now divide by nr_frames and integrate twice
218 fprintf(stderr,"\n\nRead %d frames from trajectory. Calculating potential"
219 "in spherical coordinates\n", nr_frames);
221 fprintf(stderr,"\n\nRead %d frames from trajectory. Calculating potential\n",
224 for (n =0; n < nr_grps; n++)
226 for (i = 0; i < *nslices; i++)
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
233 slVolume = 4*M_PI * sqr(i) * sqr(*slWidth) * *slWidth;
236 (*slCharge)[n][i] = 0;
240 (*slCharge)[n][i] = (*slCharge)[n][i] / (nr_frames * slVolume);
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]);
250 /* Now we have charge densities */
253 if(bCorrect && !bSpherical)
255 for(n =0; n < nr_grps; n++)
259 for (i = 0; i < *nslices; i++)
261 if( fabs((*slCharge)[n][i]) >= GMX_DOUBLE_MIN )
264 qsum += (*slCharge)[n][i];
268 for (i = 0; i < *nslices; i++)
270 if( fabs((*slCharge)[n][i]) >= GMX_DOUBLE_MIN )
272 (*slCharge)[n][i] -= qsum;
278 for(n =0; n < nr_grps; n++)
280 /* integrate twice to get field and potential */
281 p_integrate((*slField)[n], (*slCharge)[n], *nslices, *slWidth);
285 if(bCorrect && !bSpherical)
287 for(n =0; n < nr_grps; n++)
291 for (i = 0; i < *nslices; i++)
293 if( fabs((*slCharge)[n][i]) >= GMX_DOUBLE_MIN )
296 qsum += (*slField)[n][i];
300 for (i = 0; i < *nslices; i++)
302 if( fabs((*slCharge)[n][i]) >= GMX_DOUBLE_MIN )
304 (*slField)[n][i] -= qsum;
310 for(n =0; n < nr_grps; n++)
312 p_integrate((*slPotential)[n],(*slField)[n], *nslices, *slWidth);
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++)
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));
328 (*slPotential)[n][i] = ELC * (*slPotential)[n][i] * -1.0E9 / EPS0 ;
329 (*slField)[n][i] = ELC * (*slField)[n][i] * 1E18 / EPS0;
333 sfree(x0); /* free memory used by coordinate array */
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)
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 */
347 sprintf(buf,"Electrostatic Potential");
348 pot = xvgropen(afile, buf, "Box (nm)","Potential (V)",oenv);
349 xvgr_legend(pot,nr_grps,grpname,oenv);
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);
355 sprintf(buf, "Electric Field");
356 fie = xvgropen(cfile, buf, "Box (nm)", "Field (V/nm)",oenv);
357 xvgr_legend(fie,nr_grps,grpname,oenv);
359 for (slice = cb; slice < (nslices - ce); slice++)
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++)
366 fprintf(pot," %20.16g", potential[n][slice]);
367 fprintf(fie," %20.16g", field[n][slice]);
368 fprintf(cha," %20.16g", charge[n][slice]);
380 int gmx_potential(int argc,char *argv[])
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."
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 bool bSpherical = FALSE; /* default is bilayer types */
398 static real fudge_z = 0; /* translate coordinates */
399 static bool bCorrect = 0;
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" }
419 const char *bugs[] = {
420 "Discarding slices for integration should not be necessary."
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 */
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 */
441 #define NFILE asize(fnm)
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,
449 axis = toupper(axtitle[0]) - 'X';
451 top = read_top(ftp2fn(efTPX,NFILE,fnm),&ePBC); /* read topology file */
457 rd_index(ftp2fn(efNDX,NFILE,fnm),ngrps,ngx,index,grpname);
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);
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);
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 */