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
64 /****************************************************************************/
65 /* This program calculates the partial density across the box. */
66 /* Peter Tieleman, Mei 1995 */
67 /****************************************************************************/
69 /* used for sorting the list */
70 int compare(void *a, void *b)
72 t_electron *tmp1,*tmp2;
73 tmp1 = (t_electron *)a; tmp2 = (t_electron *)b;
75 return strcmp(tmp1->atomname,tmp2->atomname);
78 int get_electrons(t_electron **eltab, const char *fn)
80 char buffer[256]; /* to read in a line */
81 char tempname[80]; /* buffer to hold name */
85 int nr; /* number of atomstypes to read */
88 if ( !(in = ffopen(fn,"r")))
89 gmx_fatal(FARGS,"Couldn't open %s. Exiting.\n",fn);
91 if(NULL==fgets(buffer, 255, in))
93 gmx_fatal(FARGS,"Error reading from file %s",fn);
96 if (sscanf(buffer, "%d", &nr) != 1)
97 gmx_fatal(FARGS,"Invalid number of atomtypes in datafile\n");
102 if (fgets(buffer, 255, in) == NULL)
103 gmx_fatal(FARGS,"reading datafile. Check your datafile.\n");
104 if (sscanf(buffer, "%s = %d", tempname, &tempnr) != 2)
105 gmx_fatal(FARGS,"Invalid line in datafile at line %d\n",i+1);
106 (*eltab)[i].nr_el = tempnr;
107 (*eltab)[i].atomname = strdup(tempname);
112 fprintf(stderr,"Sorting list..\n");
113 qsort ((void*)*eltab, nr, sizeof(t_electron),
114 (int(*)(const void*, const void*))compare);
119 void center_coords(t_atoms *atoms,matrix box,rvec x0[],int axis)
123 rvec com,shift,box_center;
127 for(i=0; (i<atoms->nr); i++) {
128 mm = atoms->atom[i].m;
130 for(m=0; (m<DIM); m++)
131 com[m] += mm*x0[i][m];
133 for(m=0; (m<DIM); m++)
135 calc_box_center(ecenterDEF,box,box_center);
136 rvec_sub(box_center,com,shift);
137 shift[axis] -= box_center[axis];
139 for(i=0; (i<atoms->nr); i++)
140 rvec_dec(x0[i],shift);
143 void calc_electron_density(const char *fn, atom_id **index, int gnx[],
144 real ***slDensity, int *nslices, t_topology *top,
146 int axis, int nr_grps, real *slWidth,
147 t_electron eltab[], int nr,gmx_bool bCenter,
148 const output_env_t oenv)
150 rvec *x0; /* coordinates without pbc */
151 matrix box; /* box (3x3) */
152 int natoms; /* nr. atoms in trj */
154 int i,n, /* loop indices */
156 nr_frames = 0, /* number of frames */
157 slice; /* current slice */
158 t_electron *found; /* found by bsearch */
159 t_electron sought; /* thingie thought by bsearch */
160 gmx_rmpbc_t gpbc=NULL;
176 gmx_fatal(FARGS,"Invalid axes. Terminating\n");
179 if ((natoms = read_first_x(oenv,&status,fn,&t,&x0,box)) == 0)
180 gmx_fatal(FARGS,"Could not read coordinates from statusfile\n");
183 *nslices = (int)(box[axis][axis] * 10); /* default value */
184 fprintf(stderr,"\nDividing the box in %d slices\n",*nslices);
186 snew(*slDensity, nr_grps);
187 for (i = 0; i < nr_grps; i++)
188 snew((*slDensity)[i], *nslices);
190 gpbc = gmx_rmpbc_init(&top->idef,ePBC,top->atoms.nr,box);
191 /*********** Start processing trajectory ***********/
193 gmx_rmpbc(gpbc,natoms,box,x0);
196 center_coords(&top->atoms,box,x0,axis);
198 *slWidth = box[axis][axis]/(*nslices);
199 for (n = 0; n < nr_grps; n++) {
200 for (i = 0; i < gnx[n]; i++) { /* loop over all atoms in index file */
201 z = x0[index[n][i]][axis];
203 z += box[axis][axis];
204 while (z > box[axis][axis])
205 z -= box[axis][axis];
207 /* determine which slice atom is in */
208 slice = (z / (*slWidth));
210 sought.atomname = strdup(*(top->atoms.atomname[index[n][i]]));
212 /* now find the number of electrons. This is not efficient. */
213 found = (t_electron *)
214 bsearch((const void *)&sought,
215 (const void *)eltab, nr, sizeof(t_electron),
216 (int(*)(const void*, const void*))compare);
219 fprintf(stderr,"Couldn't find %s. Add it to the .dat file\n",
220 *(top->atoms.atomname[index[n][i]]));
222 (*slDensity)[n][slice] += found->nr_el -
223 top->atoms.atom[index[n][i]].q;
224 free(sought.atomname);
228 } while (read_next_x(oenv,status,&t,natoms,x0,box));
229 gmx_rmpbc_done(gpbc);
231 /*********** done with status file **********/
234 /* slDensity now contains the total number of electrons per slice, summed
235 over all frames. Now divide by nr_frames and volume of slice
238 fprintf(stderr,"\nRead %d frames from trajectory. Counting electrons\n",
241 for (n =0; n < nr_grps; n++) {
242 for (i = 0; i < *nslices; i++)
243 (*slDensity)[n][i] = (*slDensity)[n][i] * (*nslices) /
244 ( nr_frames * box[axis][axis] * box[ax1][ax1] * box[ax2][ax2]);
247 sfree(x0); /* free memory used by coordinate array */
250 void calc_density(const char *fn, atom_id **index, int gnx[],
251 real ***slDensity, int *nslices, t_topology *top, int ePBC,
252 int axis, int nr_grps, real *slWidth, gmx_bool bCenter,
253 const output_env_t oenv)
255 rvec *x0; /* coordinates without pbc */
256 matrix box; /* box (3x3) */
257 int natoms; /* nr. atoms in trj */
259 int **slCount, /* nr. of atoms in one slice for a group */
260 i,j,n, /* loop indices */
263 nr_frames = 0, /* number of frames */
264 slice; /* current slice */
267 char *buf; /* for tmp. keeping atomname */
268 gmx_rmpbc_t gpbc=NULL;
281 gmx_fatal(FARGS,"Invalid axes. Terminating\n");
284 if ((natoms = read_first_x(oenv,&status,fn,&t,&x0,box)) == 0)
285 gmx_fatal(FARGS,"Could not read coordinates from statusfile\n");
288 *nslices = (int)(box[axis][axis] * 10); /* default value */
289 fprintf(stderr,"\nDividing the box in %d slices\n",*nslices);
292 snew(*slDensity, nr_grps);
293 for (i = 0; i < nr_grps; i++)
294 snew((*slDensity)[i], *nslices);
296 gpbc = gmx_rmpbc_init(&top->idef,ePBC,top->atoms.nr,box);
297 /*********** Start processing trajectory ***********/
299 gmx_rmpbc(gpbc,natoms,box,x0);
302 center_coords(&top->atoms,box,x0,axis);
304 *slWidth = box[axis][axis]/(*nslices);
307 for (n = 0; n < nr_grps; n++) {
308 for (i = 0; i < gnx[n]; i++) { /* loop over all atoms in index file */
309 z = x0[index[n][i]][axis];
311 z += box[axis][axis];
312 while (z > box[axis][axis])
313 z -= box[axis][axis];
315 /* determine which slice atom is in */
316 slice = (int)(z / (*slWidth));
317 (*slDensity)[n][slice] += top->atoms.atom[index[n][i]].m;
322 } while (read_next_x(oenv,status,&t,natoms,x0,box));
323 gmx_rmpbc_done(gpbc);
325 /*********** done with status file **********/
328 /* slDensity now contains the total mass per slice, summed over all
329 frames. Now divide by nr_frames and volume of slice
332 fprintf(stderr,"\nRead %d frames from trajectory. Calculating density\n",
335 for (n =0; n < nr_grps; n++) {
336 for (i = 0; i < *nslices; i++) {
337 (*slDensity)[n][i] = (*slDensity)[n][i] * (*nslices) /
338 (nr_frames * box[axis][axis] * box[ax1][ax1] * box[ax2][ax2]);
342 sfree(x0); /* free memory used by coordinate array */
345 void plot_density(real *slDensity[], const char *afile, int nslices,
346 int nr_grps, char *grpname[], real slWidth,
347 const char **dens_opt,
348 gmx_bool bSymmetrize, const output_env_t oenv)
351 const char *ylabel=NULL;
355 switch (dens_opt[0][0]) {
356 case 'm': ylabel = "Density (kg m\\S-3\\N)"; break;
357 case 'n': ylabel = "Number density (nm\\S-3\\N)"; break;
358 case 'c': ylabel = "Charge density (e nm\\S-3\\N)"; break;
359 case 'e': ylabel = "Electron density (e nm\\S-3\\N)"; break;
362 den = xvgropen(afile, "Partial densities", "Box (nm)", ylabel,oenv);
364 xvgr_legend(den,nr_grps,(const char**)grpname,oenv);
366 for (slice = 0; (slice < nslices); slice++) {
367 fprintf(den,"%12g ", slice * slWidth);
368 for (n = 0; (n < nr_grps); n++) {
370 ddd = (slDensity[n][slice]+slDensity[n][nslices-slice-1])*0.5;
372 ddd = slDensity[n][slice];
373 if (dens_opt[0][0] == 'm')
374 fprintf(den," %12g", ddd*AMU/(NANO*NANO*NANO));
376 fprintf(den," %12g", ddd);
384 int gmx_density(int argc,char *argv[])
386 const char *desc[] = {
387 "Compute partial densities across the box, using an index file. Densities",
388 "in kg/m^3, number densities or electron densities can be",
389 "calculated. For electron densities, a file describing the number of",
390 "electrons for each type of atom should be provided using [TT]-ei[tt].",
391 "It should look like:[BR]",
393 " atomname = nrelectrons[BR]",
394 " atomname = nrelectrons[BR]",
395 "The first line contains the number of lines to read from the file.",
396 "There should be one line for each unique atom name in your system.",
397 "The number of electrons for each atom is modified by its atomic",
402 static const char *dens_opt[] =
403 { NULL, "mass", "number", "charge", "electron", NULL };
404 static int axis = 2; /* normal to memb. default z */
405 static const char *axtitle="Z";
406 static int nslices = 50; /* nr of slices defined */
407 static int ngrps = 1; /* nr. of groups */
408 static gmx_bool bSymmetrize=FALSE;
409 static gmx_bool bCenter=FALSE;
411 { "-d", FALSE, etSTR, {&axtitle},
412 "Take the normal on the membrane in direction X, Y or Z." },
413 { "-sl", FALSE, etINT, {&nslices},
414 "Divide the box in #nr slices." },
415 { "-dens", FALSE, etENUM, {dens_opt},
417 { "-ng", FALSE, etINT, {&ngrps},
418 "Number of groups to compute densities of" },
419 { "-symm", FALSE, etBOOL, {&bSymmetrize},
420 "Symmetrize the density along the axis, with respect to the center. Useful for bilayers." },
421 { "-center", FALSE, etBOOL, {&bCenter},
422 "Shift the center of mass along the axis to zero. This means if your axis is Z and your box is bX, bY, bZ, the center of mass will be at bX/2, bY/2, 0."}
425 const char *bugs[] = {
426 "When calculating electron densities, atomnames are used instead of types. This is bad.",
429 real **density; /* density per slice */
430 real slWidth; /* width of one slice */
431 char **grpname; /* groupnames */
432 int nr_electrons; /* nr. electrons */
433 int *ngx; /* sizes of groups */
434 t_electron *el_tab; /* tabel with nr. of electrons*/
435 t_topology *top; /* topology */
437 atom_id **index; /* indices for all groups */
440 t_filenm fnm[] = { /* files for g_density */
441 { efTRX, "-f", NULL, ffREAD },
442 { efNDX, NULL, NULL, ffOPTRD },
443 { efTPX, NULL, NULL, ffREAD },
444 { efDAT, "-ei", "electrons", ffOPTRD }, /* file with nr. of electrons */
445 { efXVG,"-o","density",ffWRITE },
448 #define NFILE asize(fnm)
450 CopyRight(stderr,argv[0]);
452 parse_common_args(&argc,argv,PCA_CAN_VIEW | PCA_CAN_TIME | PCA_BE_NICE,
453 NFILE,fnm,asize(pa),pa,asize(desc),desc,asize(bugs),bugs,
456 if (bSymmetrize && !bCenter) {
457 fprintf(stderr,"Can not symmetrize without centering. Turning on -center\n");
461 axis = toupper(axtitle[0]) - 'X';
463 top = read_top(ftp2fn(efTPX,NFILE,fnm),&ePBC); /* read topology file */
464 if (dens_opt[0][0] == 'n') {
465 for(i=0; (i<top->atoms.nr); i++)
466 top->atoms.atom[i].m = 1;
467 } else if (dens_opt[0][0] == 'c') {
468 for(i=0; (i<top->atoms.nr); i++)
469 top->atoms.atom[i].m = top->atoms.atom[i].q;
476 get_index(&top->atoms,ftp2fn_null(efNDX,NFILE,fnm),ngrps,ngx,index,grpname);
478 if (dens_opt[0][0] == 'e') {
479 nr_electrons = get_electrons(&el_tab,ftp2fn(efDAT,NFILE,fnm));
480 fprintf(stderr,"Read %d atomtypes from datafile\n", nr_electrons);
482 calc_electron_density(ftp2fn(efTRX,NFILE,fnm),index, ngx, &density,
483 &nslices, top, ePBC, axis, ngrps, &slWidth, el_tab,
484 nr_electrons,bCenter,oenv);
486 calc_density(ftp2fn(efTRX,NFILE,fnm),index, ngx, &density, &nslices, top,
487 ePBC, axis, ngrps, &slWidth, bCenter,oenv);
489 plot_density(density, opt2fn("-o",NFILE,fnm),
490 nslices, ngrps, grpname, slWidth, dens_opt,
493 do_view(oenv,opt2fn("-o",NFILE,fnm), "-nxy"); /* view xvgr file */