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
49 #include "gromacs/fileio/futil.h"
51 #include "gromacs/fileio/tpxio.h"
52 #include "gromacs/fileio/trxio.h"
60 #include "gromacs/fileio/matio.h"
64 static void check_box_c(matrix box)
66 if (fabs(box[ZZ][XX]) > GMX_REAL_EPS*box[ZZ][ZZ] ||
67 fabs(box[ZZ][YY]) > GMX_REAL_EPS*box[ZZ][ZZ])
70 "The last box vector is not parallel to the z-axis: %f %f %f",
71 box[ZZ][XX], box[ZZ][YY], box[ZZ][ZZ]);
75 static void calc_comg(int is, int *coi, int *index, gmx_bool bMass, t_atom *atom,
76 rvec *x, rvec *x_comg)
82 if (bMass && atom == NULL)
84 gmx_fatal(FARGS, "No masses available while mass weighting was requested");
87 for (c = 0; c < is; c++)
91 for (i = coi[c]; i < coi[c+1]; i++)
96 for (d = 0; d < DIM; d++)
98 xc[d] += m*x[index[i]][d];
104 rvec_inc(xc, x[index[i]]);
108 svmul(1/mtot, xc, x_comg[c]);
112 static void split_group(int isize, int *index, char *grpname,
113 t_topology *top, char type,
114 int *is_out, int **coi_out)
116 t_block *mols = NULL;
119 int cur, mol, res, i, a, i1;
121 /* Split up the group in molecules or residues */
128 atom = top->atoms.atom;
131 gmx_fatal(FARGS, "Unknown rdf option '%s'", type);
137 for (i = 0; i < isize; i++)
142 /* Check if the molecule number has changed */
143 i1 = mols->index[mol+1];
147 i1 = mols->index[mol+1];
155 else if (type == 'r')
157 /* Check if the residue index has changed */
158 res = atom[a].resind;
168 printf("Group '%s' of %d atoms consists of %d %s\n",
170 (type == 'm' ? "molecules" : "residues"));
176 static void do_rdf(const char *fnNDX, const char *fnTPS, const char *fnTRX,
177 const char *fnRDF, const char *fnCNRDF, const char *fnHQ,
178 gmx_bool bCM, const char *close,
179 const char **rdft, gmx_bool bXY, gmx_bool bPBC, gmx_bool bNormalize,
180 real cutoff, real binwidth, real fade, int ng,
181 const output_env_t oenv)
185 char outf1[STRLEN], outf2[STRLEN];
186 char title[STRLEN], gtitle[STRLEN], refgt[30];
187 int g, natoms, i, ii, j, k, nbin, j0, j1, n, nframes;
190 int *isize, isize_cm = 0, nrdf = 0, max_i, isize0, isize_g;
191 atom_id **index, *index_cm = NULL;
192 #if (defined SIZEOF_LONG_LONG_INT) && (SIZEOF_LONG_LONG_INT >= 8)
197 real t, rmax2, cut2, r, r2, r2ii, invhbinw, normfac;
198 real segvol, spherevol, prev_spherevol, **rdf;
199 rvec *x, dx, *x0 = NULL, *x_i1, xi;
200 real *inv_segvol, invvol, invvol_sum, rho;
201 gmx_bool bClose, *bExcl, bTop, bNonSelfExcl;
204 atom_id ix, jx, ***pairs;
205 t_topology *top = NULL;
206 int ePBC = -1, ePBCrdf = -1;
207 t_block *mols = NULL;
211 gmx_rmpbc_t gpbc = NULL;
212 int *is = NULL, **coi = NULL, cur, mol, i1, res, a;
216 bClose = (close[0] != 'n');
221 bTop = read_tps_conf(fnTPS, title, top, &ePBC, &x, NULL, box, TRUE);
222 if (bTop && !(bCM || bClose))
224 /* get exclusions from topology */
225 excl = &(top->excls);
231 fprintf(stderr, "\nSelect a reference group and %d group%s\n",
232 ng, ng == 1 ? "" : "s");
235 get_index(&(top->atoms), fnNDX, ng+1, isize, index, grpname);
236 atom = top->atoms.atom;
240 rd_index(fnNDX, ng+1, isize, index, grpname);
249 split_group(isize[0], index[0], grpname[0], top, close[0], &is[0], &coi[0]);
252 if (rdft[0][0] != 'a')
254 /* Split up all the groups in molecules or residues */
257 for (g = ((bCM || bClose) ? 1 : 0); g < ng+1; g++)
259 split_group(isize[g], index[g], grpname[g], top, rdft[0][0], &is[g], &coi[g]);
266 snew(coi[0], is[0]+1);
268 coi[0][1] = isize[0];
272 else if (bClose || rdft[0][0] != 'a')
282 natoms = read_first_x(oenv, &status, fnTRX, &t, &x, box);
285 gmx_fatal(FARGS, "Could not read coordinates from statusfile\n");
289 /* check with topology */
290 if (natoms > top->atoms.nr)
292 gmx_fatal(FARGS, "Trajectory (%d atoms) does not match topology (%d atoms)",
293 natoms, top->atoms.nr);
296 /* check with index groups */
297 for (i = 0; i < ng+1; i++)
299 for (j = 0; j < isize[i]; j++)
301 if (index[i][j] >= natoms)
303 gmx_fatal(FARGS, "Atom index (%d) in index group %s (%d atoms) larger "
304 "than number of atoms in trajectory (%d atoms)",
305 index[i][j], grpname[i], isize[i], natoms);
310 /* initialize some handy things */
313 ePBC = guess_ePBC(box);
315 copy_mat(box, box_pbc);
322 case epbcXY: ePBCrdf = epbcXY; break;
323 case epbcNONE: ePBCrdf = epbcNONE; break;
325 gmx_fatal(FARGS, "xy-rdf's are not supported for pbc type'%s'",
329 /* Make sure the z-height does not influence the cut-off */
330 box_pbc[ZZ][ZZ] = 2*max(box[XX][XX], box[YY][YY]);
338 rmax2 = 0.99*0.99*max_cutoff2(bXY ? epbcXY : epbcXYZ, box_pbc);
342 rmax2 = sqr(3*max(box[XX][XX], max(box[YY][YY], box[ZZ][ZZ])));
346 fprintf(debug, "rmax2 = %g\n", rmax2);
349 /* We use the double amount of bins, so we can correctly
350 * write the rdf and rdf_cn output at i*binwidth values.
352 nbin = (int)(sqrt(rmax2) * 2 / binwidth);
353 invhbinw = 2.0 / binwidth;
362 for (g = 0; g < ng; g++)
364 if (isize[g+1] > max_i)
369 /* this is THE array */
370 snew(count[g], nbin+1);
372 /* make pairlist array for groups and exclusions */
373 snew(pairs[g], isize[0]);
374 snew(npairs[g], isize[0]);
375 for (i = 0; i < isize[0]; i++)
377 /* We can only have exclusions with atomic rdfs */
378 if (!(bCM || bClose || rdft[0][0] != 'a'))
381 for (j = 0; j < natoms; j++)
388 for (j = excl->index[ix]; j < excl->index[ix+1]; j++)
390 bExcl[excl->a[j]] = TRUE;
394 snew(pairs[g][i], isize[g+1]);
395 bNonSelfExcl = FALSE;
396 for (j = 0; j < isize[g+1]; j++)
401 pairs[g][i][k++] = jx;
405 /* Check if we have exclusions other than self exclusions */
412 srenew(pairs[g][i], npairs[g][i]);
416 /* Save a LOT of memory and some cpu cycles */
432 if (bPBC && (NULL != top))
434 gpbc = gmx_rmpbc_init(&top->idef, ePBC, natoms);
439 /* Must init pbc every step because of pressure coupling */
440 copy_mat(box, box_pbc);
445 gmx_rmpbc(gpbc, natoms, box, x);
450 clear_rvec(box_pbc[ZZ]);
452 set_pbc(&pbc, ePBCrdf, box_pbc);
456 /* Set z-size to 1 so we get the surface iso the volume */
460 invvol = 1/det(box_pbc);
461 invvol_sum += invvol;
465 /* Calculate center of mass of the whole group */
466 calc_comg(is[0], coi[0], index[0], TRUE, atom, x, x0);
468 else if (!bClose && rdft[0][0] != 'a')
470 calc_comg(is[0], coi[0], index[0], rdft[0][6] == 'm', atom, x, x0);
473 for (g = 0; g < ng; g++)
475 if (rdft[0][0] == 'a')
477 /* Copy the indexed coordinates to a continuous array */
478 for (i = 0; i < isize[g+1]; i++)
480 copy_rvec(x[index[g+1][i]], x_i1[i]);
485 /* Calculate the COMs/COGs and store in x_i1 */
486 calc_comg(is[g+1], coi[g+1], index[g+1], rdft[0][6] == 'm', atom, x, x_i1);
489 for (i = 0; i < isize0; i++)
493 /* Special loop, since we need to determine the minimum distance
494 * over all selected atoms in the reference molecule/residue.
496 if (rdft[0][0] == 'a')
498 isize_g = isize[g+1];
504 for (j = 0; j < isize_g; j++)
507 /* Loop over the selected atoms in the reference molecule */
508 for (ii = coi[0][i]; ii < coi[0][i+1]; ii++)
512 pbc_dx(&pbc, x[index[0][ii]], x_i1[j], dx);
516 rvec_sub(x[index[0][ii]], x_i1[j], dx);
520 r2ii = dx[XX]*dx[XX] + dx[YY]*dx[YY];
524 r2ii = iprod(dx, dx);
531 if (r2 > cut2 && r2 <= rmax2)
533 count[g][(int)(sqrt(r2)*invhbinw)]++;
539 /* Real rdf between points in space */
540 if (bCM || rdft[0][0] != 'a')
542 copy_rvec(x0[i], xi);
546 copy_rvec(x[index[0][i]], xi);
548 if (rdft[0][0] == 'a' && npairs[g][i] >= 0)
550 /* Expensive loop, because of indexing */
551 for (j = 0; j < npairs[g][i]; j++)
556 pbc_dx(&pbc, xi, x[jx], dx);
560 rvec_sub(xi, x[jx], dx);
565 r2 = dx[XX]*dx[XX] + dx[YY]*dx[YY];
571 if (r2 > cut2 && r2 <= rmax2)
573 count[g][(int)(sqrt(r2)*invhbinw)]++;
579 /* Cheaper loop, no exclusions */
580 if (rdft[0][0] == 'a')
582 isize_g = isize[g+1];
588 for (j = 0; j < isize_g; j++)
592 pbc_dx(&pbc, xi, x_i1[j], dx);
596 rvec_sub(xi, x_i1[j], dx);
600 r2 = dx[XX]*dx[XX] + dx[YY]*dx[YY];
606 if (r2 > cut2 && r2 <= rmax2)
608 count[g][(int)(sqrt(r2)*invhbinw)]++;
617 while (read_next_x(oenv, status, &t, x, box));
618 fprintf(stderr, "\n");
620 if (bPBC && (NULL != top))
622 gmx_rmpbc_done(gpbc);
630 invvol = invvol_sum/nframes;
632 /* Calculate volume of sphere segments or length of circle segments */
633 snew(inv_segvol, (nbin+1)/2);
635 for (i = 0; (i < (nbin+1)/2); i++)
637 r = (i + 0.5)*binwidth;
640 spherevol = M_PI*r*r;
644 spherevol = (4.0/3.0)*M_PI*r*r*r;
646 segvol = spherevol-prev_spherevol;
647 inv_segvol[i] = 1.0/segvol;
648 prev_spherevol = spherevol;
652 for (g = 0; g < ng; g++)
654 /* We have to normalize by dividing by the number of frames */
655 if (rdft[0][0] == 'a')
657 normfac = 1.0/(nframes*invvol*isize0*isize[g+1]);
661 normfac = 1.0/(nframes*invvol*isize0*is[g+1]);
664 /* Do the normalization */
665 nrdf = max((nbin+1)/2, 1+2*fade/binwidth);
667 for (i = 0; i < (nbin+1)/2; i++)
676 j = count[g][i*2-1] + count[g][i*2];
678 if ((fade > 0) && (r >= fade))
680 rdf[g][i] = 1 + (j*inv_segvol[i]*normfac-1)*exp(-16*sqr(r/fade-1));
686 rdf[g][i] = j*inv_segvol[i]*normfac;
690 rdf[g][i] = j/(binwidth*isize0*nframes);
694 for (; (i < nrdf); i++)
700 if (rdft[0][0] == 'a')
702 sprintf(gtitle, "Radial distribution");
706 sprintf(gtitle, "Radial distribution of %s %s",
707 rdft[0][0] == 'm' ? "molecule" : "residue",
708 rdft[0][6] == 'm' ? "COM" : "COG");
710 fp = xvgropen(fnRDF, gtitle, "r", "", oenv);
713 sprintf(refgt, " %s", "COM");
717 sprintf(refgt, " closest atom in %s.", close);
721 sprintf(refgt, "%s", "");
725 if (output_env_get_print_xvgr_codes(oenv))
727 fprintf(fp, "@ subtitle \"%s%s - %s\"\n", grpname[0], refgt, grpname[1]);
732 if (output_env_get_print_xvgr_codes(oenv))
734 fprintf(fp, "@ subtitle \"reference %s%s\"\n", grpname[0], refgt);
736 xvgr_legend(fp, ng, (const char**)(grpname+1), oenv);
738 for (i = 0; (i < nrdf); i++)
740 fprintf(fp, "%10g", i*binwidth);
741 for (g = 0; g < ng; g++)
743 fprintf(fp, " %10g", rdf[g][i]);
749 do_view(oenv, fnRDF, NULL);
751 /* h(Q) function: fourier transform of rdf */
755 real *hq, *integrand, Q;
757 /* Get a better number density later! */
758 rho = isize[1]*invvol;
760 snew(integrand, nrdf);
761 for (i = 0; (i < nhq); i++)
765 for (j = 1; (j < nrdf); j++)
768 integrand[j] = (Q == 0) ? 1.0 : sin(Q*r)/(Q*r);
769 integrand[j] *= 4.0*M_PI*rho*r*r*(rdf[0][j]-1.0);
771 hq[i] = print_and_integrate(debug, nrdf, binwidth, integrand, NULL, 0);
773 fp = xvgropen(fnHQ, "h(Q)", "Q(/nm)", "h(Q)", oenv);
774 for (i = 0; (i < nhq); i++)
776 fprintf(fp, "%10g %10g\n", i*0.5, hq[i]);
779 do_view(oenv, fnHQ, NULL);
786 normfac = 1.0/(isize0*nframes);
787 fp = xvgropen(fnCNRDF, "Cumulative Number RDF", "r", "number", oenv);
790 if (output_env_get_print_xvgr_codes(oenv))
792 fprintf(fp, "@ subtitle \"%s-%s\"\n", grpname[0], grpname[1]);
797 if (output_env_get_print_xvgr_codes(oenv))
799 fprintf(fp, "@ subtitle \"reference %s\"\n", grpname[0]);
801 xvgr_legend(fp, ng, (const char**)(grpname+1), oenv);
804 for (i = 0; (i <= nbin/2); i++)
806 fprintf(fp, "%10g", i*binwidth);
807 for (g = 0; g < ng; g++)
809 fprintf(fp, " %10g", (real)((double)sum[g]*normfac));
812 sum[g] += count[g][i*2] + count[g][i*2+1];
820 do_view(oenv, fnCNRDF, NULL);
823 for (g = 0; g < ng; g++)
831 int gmx_rdf(int argc, char *argv[])
833 const char *desc[] = {
834 "The structure of liquids can be studied by either neutron or X-ray",
835 "scattering. The most common way to describe liquid structure is by a",
836 "radial distribution function. However, this is not easy to obtain from",
837 "a scattering experiment.[PAR]",
838 "[TT]g_rdf[tt] calculates radial distribution functions in different ways.",
839 "The normal method is around a (set of) particle(s), the other methods",
840 "are around the center of mass of a set of particles ([TT]-com[tt])",
841 "or to the closest particle in a set ([TT]-surf[tt]).",
842 "With all methods, the RDF can also be calculated around axes parallel",
843 "to the [IT]z[it]-axis with option [TT]-xy[tt].",
844 "With option [TT]-surf[tt] normalization can not be used.[PAR]",
845 "The option [TT]-rdf[tt] sets the type of RDF to be computed.",
846 "Default is for atoms or particles, but one can also select center",
847 "of mass or geometry of molecules or residues. In all cases, only",
848 "the atoms in the index groups are taken into account.",
849 "For molecules and/or the center of mass option, a run input file",
851 "Weighting other than COM or COG can currently only be achieved",
852 "by providing a run input file with different masses.",
853 "Options [TT]-com[tt] and [TT]-surf[tt] also work in conjunction",
854 "with [TT]-rdf[tt].[PAR]",
855 "If a run input file is supplied ([TT]-s[tt]) and [TT]-rdf[tt] is set",
856 "to [TT]atom[tt], exclusions defined",
857 "in that file are taken into account when calculating the RDF.",
858 "The option [TT]-cut[tt] is meant as an alternative way to avoid",
859 "intramolecular peaks in the RDF plot.",
860 "It is however better to supply a run input file with a higher number of",
861 "exclusions. For e.g. benzene a topology, setting nrexcl to 5",
862 "would eliminate all intramolecular contributions to the RDF.",
863 "Note that all atoms in the selected groups are used, also the ones",
864 "that don't have Lennard-Jones interactions.[PAR]",
865 "Option [TT]-cn[tt] produces the cumulative number RDF,",
866 "i.e. the average number of particles within a distance r.[PAR]"
868 static gmx_bool bCM = FALSE, bXY = FALSE, bPBC = TRUE, bNormalize = TRUE;
869 static real cutoff = 0, binwidth = 0.002, fade = 0.0;
870 static int ngroups = 1;
872 static const char *closet[] = { NULL, "no", "mol", "res", NULL };
873 static const char *rdft[] = { NULL, "atom", "mol_com", "mol_cog", "res_com", "res_cog", NULL };
876 { "-bin", FALSE, etREAL, {&binwidth},
878 { "-com", FALSE, etBOOL, {&bCM},
879 "RDF with respect to the center of mass of first group" },
880 { "-surf", FALSE, etENUM, {closet},
881 "RDF with respect to the surface of the first group" },
882 { "-rdf", FALSE, etENUM, {rdft},
884 { "-pbc", FALSE, etBOOL, {&bPBC},
885 "Use periodic boundary conditions for computing distances. Without PBC the maximum range will be three times the largest box edge." },
886 { "-norm", FALSE, etBOOL, {&bNormalize},
887 "Normalize for volume and density" },
888 { "-xy", FALSE, etBOOL, {&bXY},
889 "Use only the x and y components of the distance" },
890 { "-cut", FALSE, etREAL, {&cutoff},
891 "Shortest distance (nm) to be considered"},
892 { "-ng", FALSE, etINT, {&ngroups},
893 "Number of secondary groups to compute RDFs around a central group" },
894 { "-fade", FALSE, etREAL, {&fade},
895 "From this distance onwards the RDF is tranformed by g'(r) = 1 + [g(r)-1] exp(-(r/fade-1)^2 to make it go to 1 smoothly. If fade is 0.0 nothing is done." }
897 #define NPA asize(pa)
898 const char *fnTPS, *fnNDX;
902 { efTRX, "-f", NULL, ffREAD },
903 { efTPS, NULL, NULL, ffOPTRD },
904 { efNDX, NULL, NULL, ffOPTRD },
905 { efXVG, "-o", "rdf", ffWRITE },
906 { efXVG, "-cn", "rdf_cn", ffOPTWR },
907 { efXVG, "-hq", "hq", ffOPTWR },
909 #define NFILE asize(fnm)
910 if (!parse_common_args(&argc, argv, PCA_CAN_VIEW | PCA_CAN_TIME | PCA_BE_NICE,
911 NFILE, fnm, NPA, pa, asize(desc), desc, 0, NULL, &oenv))
916 if (bCM || closet[0][0] != 'n' || rdft[0][0] == 'm' || rdft[0][6] == 'm')
918 fnTPS = ftp2fn(efTPS, NFILE, fnm);
922 fnTPS = ftp2fn_null(efTPS, NFILE, fnm);
924 fnNDX = ftp2fn_null(efNDX, NFILE, fnm);
926 if (!fnTPS && !fnNDX)
928 gmx_fatal(FARGS, "Neither index file nor topology file specified\n"
932 if (closet[0][0] != 'n')
936 gmx_fatal(FARGS, "Can not have both -com and -surf");
940 fprintf(stderr, "Turning of normalization because of option -surf\n");
945 do_rdf(fnNDX, fnTPS, ftp2fn(efTRX, NFILE, fnm),
946 opt2fn("-o", NFILE, fnm), opt2fn_null("-cn", NFILE, fnm),
947 opt2fn_null("-hq", NFILE, fnm),
948 bCM, closet[0], rdft, bXY, bPBC, bNormalize, cutoff, binwidth, fade, ngroups,