2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5 * Copyright (c) 2001-2004, The GROMACS development team.
6 * Copyright (c) 2013,2014, by the GROMACS development team, led by
7 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
8 * and including many others, as listed in the AUTHORS file in the
9 * top-level source directory and at http://www.gromacs.org.
11 * GROMACS is free software; you can redistribute it and/or
12 * modify it under the terms of the GNU Lesser General Public License
13 * as published by the Free Software Foundation; either version 2.1
14 * of the License, or (at your option) any later version.
16 * GROMACS is distributed in the hope that it will be useful,
17 * but WITHOUT ANY WARRANTY; without even the implied warranty of
18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19 * Lesser General Public License for more details.
21 * You should have received a copy of the GNU Lesser General Public
22 * License along with GROMACS; if not, see
23 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
24 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
26 * If you want to redistribute modifications to GROMACS, please
27 * consider that scientific software is very special. Version
28 * control is crucial - bugs must be traceable. We will be happy to
29 * consider code for inclusion in the official distribution, but
30 * derived work must not be called official GROMACS. Details are found
31 * in the README & COPYING files - if they are missing, get the
32 * official version at http://www.gromacs.org.
34 * To help us fund GROMACS development, we humbly ask that you cite
35 * the research papers on the package. Check out http://www.gromacs.org.
46 #include "gromacs/math/vec.h"
47 #include "gromacs/pbcutil/pbc.h"
49 #include "gromacs/utility/futil.h"
50 #include "gromacs/fileio/tpxio.h"
51 #include "gromacs/fileio/trxio.h"
60 #include "gromacs/commandline/pargs.h"
61 #include "gromacs/fileio/xvgr.h"
62 #include "gromacs/pbcutil/rmpbc.h"
63 #include "gromacs/utility/fatalerror.h"
64 #include "gromacs/utility/smalloc.h"
66 static void check_box_c(matrix box)
68 if (fabs(box[ZZ][XX]) > GMX_REAL_EPS*box[ZZ][ZZ] ||
69 fabs(box[ZZ][YY]) > GMX_REAL_EPS*box[ZZ][ZZ])
72 "The last box vector is not parallel to the z-axis: %f %f %f",
73 box[ZZ][XX], box[ZZ][YY], box[ZZ][ZZ]);
77 static void calc_comg(int is, int *coi, int *index, gmx_bool bMass, t_atom *atom,
78 rvec *x, rvec *x_comg)
84 if (bMass && atom == NULL)
86 gmx_fatal(FARGS, "No masses available while mass weighting was requested");
89 for (c = 0; c < is; c++)
93 for (i = coi[c]; i < coi[c+1]; i++)
98 for (d = 0; d < DIM; d++)
100 xc[d] += m*x[index[i]][d];
106 rvec_inc(xc, x[index[i]]);
110 svmul(1/mtot, xc, x_comg[c]);
114 static void split_group(int isize, int *index, char *grpname,
115 t_topology *top, char type,
116 int *is_out, int **coi_out)
118 t_block *mols = NULL;
121 int cur, mol, res, i, a, i1;
123 /* Split up the group in molecules or residues */
130 atom = top->atoms.atom;
133 gmx_fatal(FARGS, "Unknown rdf option '%s'", type);
139 for (i = 0; i < isize; i++)
144 /* Check if the molecule number has changed */
145 i1 = mols->index[mol+1];
149 i1 = mols->index[mol+1];
157 else if (type == 'r')
159 /* Check if the residue index has changed */
160 res = atom[a].resind;
170 printf("Group '%s' of %d atoms consists of %d %s\n",
172 (type == 'm' ? "molecules" : "residues"));
178 static void do_rdf(const char *fnNDX, const char *fnTPS, const char *fnTRX,
179 const char *fnRDF, const char *fnCNRDF, const char *fnHQ,
180 gmx_bool bCM, const char *close,
181 const char **rdft, gmx_bool bXY, gmx_bool bPBC, gmx_bool bNormalize,
182 real cutoff, real binwidth, real fade, int ng,
183 const output_env_t oenv)
187 char outf1[STRLEN], outf2[STRLEN];
188 char title[STRLEN], gtitle[STRLEN], refgt[30];
189 int g, natoms, i, ii, j, k, nbin, j0, j1, n, nframes;
192 int *isize, isize_cm = 0, nrdf = 0, max_i, isize0, isize_g;
193 atom_id **index, *index_cm = NULL;
195 real t, rmax2, cut2, r, r2, r2ii, invhbinw, normfac;
196 real segvol, spherevol, prev_spherevol, **rdf;
197 rvec *x, dx, *x0 = NULL, *x_i1, xi;
198 real *inv_segvol, invvol, invvol_sum, rho;
199 gmx_bool bClose, *bExcl, bTop, bNonSelfExcl;
202 atom_id ix, jx, ***pairs;
203 t_topology *top = NULL;
204 int ePBC = -1, ePBCrdf = -1;
205 t_block *mols = NULL;
209 gmx_rmpbc_t gpbc = NULL;
210 int *is = NULL, **coi = NULL, cur, mol, i1, res, a;
214 bClose = (close[0] != 'n');
219 bTop = read_tps_conf(fnTPS, title, top, &ePBC, &x, NULL, box, TRUE);
220 if (bTop && !(bCM || bClose))
222 /* get exclusions from topology */
223 excl = &(top->excls);
229 fprintf(stderr, "\nSelect a reference group and %d group%s\n",
230 ng, ng == 1 ? "" : "s");
233 get_index(&(top->atoms), fnNDX, ng+1, isize, index, grpname);
234 atom = top->atoms.atom;
238 rd_index(fnNDX, ng+1, isize, index, grpname);
247 split_group(isize[0], index[0], grpname[0], top, close[0], &is[0], &coi[0]);
250 if (rdft[0][0] != 'a')
252 /* Split up all the groups in molecules or residues */
255 for (g = ((bCM || bClose) ? 1 : 0); g < ng+1; g++)
257 split_group(isize[g], index[g], grpname[g], top, rdft[0][0], &is[g], &coi[g]);
264 snew(coi[0], is[0]+1);
266 coi[0][1] = isize[0];
270 else if (bClose || rdft[0][0] != 'a')
280 natoms = read_first_x(oenv, &status, fnTRX, &t, &x, box);
283 gmx_fatal(FARGS, "Could not read coordinates from statusfile\n");
287 /* check with topology */
288 if (natoms > top->atoms.nr)
290 gmx_fatal(FARGS, "Trajectory (%d atoms) does not match topology (%d atoms)",
291 natoms, top->atoms.nr);
294 /* check with index groups */
295 for (i = 0; i < ng+1; i++)
297 for (j = 0; j < isize[i]; j++)
299 if (index[i][j] >= natoms)
301 gmx_fatal(FARGS, "Atom index (%d) in index group %s (%d atoms) larger "
302 "than number of atoms in trajectory (%d atoms)",
303 index[i][j], grpname[i], isize[i], natoms);
308 /* initialize some handy things */
311 ePBC = guess_ePBC(box);
313 copy_mat(box, box_pbc);
320 case epbcXY: ePBCrdf = epbcXY; break;
321 case epbcNONE: ePBCrdf = epbcNONE; break;
323 gmx_fatal(FARGS, "xy-rdf's are not supported for pbc type'%s'",
327 /* Make sure the z-height does not influence the cut-off */
328 box_pbc[ZZ][ZZ] = 2*max(box[XX][XX], box[YY][YY]);
336 rmax2 = 0.99*0.99*max_cutoff2(bXY ? epbcXY : epbcXYZ, box_pbc);
340 rmax2 = sqr(3*max(box[XX][XX], max(box[YY][YY], box[ZZ][ZZ])));
344 fprintf(debug, "rmax2 = %g\n", rmax2);
347 /* We use the double amount of bins, so we can correctly
348 * write the rdf and rdf_cn output at i*binwidth values.
350 nbin = (int)(sqrt(rmax2) * 2 / binwidth);
351 invhbinw = 2.0 / binwidth;
360 for (g = 0; g < ng; g++)
362 if (isize[g+1] > max_i)
367 /* this is THE array */
368 snew(count[g], nbin+1);
370 /* make pairlist array for groups and exclusions */
371 snew(pairs[g], isize[0]);
372 snew(npairs[g], isize[0]);
373 for (i = 0; i < isize[0]; i++)
375 /* We can only have exclusions with atomic rdfs */
376 if (!(bCM || bClose || rdft[0][0] != 'a'))
379 for (j = 0; j < natoms; j++)
386 for (j = excl->index[ix]; j < excl->index[ix+1]; j++)
388 bExcl[excl->a[j]] = TRUE;
392 snew(pairs[g][i], isize[g+1]);
393 bNonSelfExcl = FALSE;
394 for (j = 0; j < isize[g+1]; j++)
399 pairs[g][i][k++] = jx;
403 /* Check if we have exclusions other than self exclusions */
410 srenew(pairs[g][i], npairs[g][i]);
414 /* Save a LOT of memory and some cpu cycles */
430 if (bPBC && (NULL != top))
432 gpbc = gmx_rmpbc_init(&top->idef, ePBC, natoms);
437 /* Must init pbc every step because of pressure coupling */
438 copy_mat(box, box_pbc);
443 gmx_rmpbc(gpbc, natoms, box, x);
448 clear_rvec(box_pbc[ZZ]);
450 set_pbc(&pbc, ePBCrdf, box_pbc);
454 /* Set z-size to 1 so we get the surface iso the volume */
458 invvol = 1/det(box_pbc);
459 invvol_sum += invvol;
463 /* Calculate center of mass of the whole group */
464 calc_comg(is[0], coi[0], index[0], TRUE, atom, x, x0);
466 else if (!bClose && rdft[0][0] != 'a')
468 calc_comg(is[0], coi[0], index[0], rdft[0][6] == 'm', atom, x, x0);
471 for (g = 0; g < ng; g++)
473 if (rdft[0][0] == 'a')
475 /* Copy the indexed coordinates to a continuous array */
476 for (i = 0; i < isize[g+1]; i++)
478 copy_rvec(x[index[g+1][i]], x_i1[i]);
483 /* Calculate the COMs/COGs and store in x_i1 */
484 calc_comg(is[g+1], coi[g+1], index[g+1], rdft[0][6] == 'm', atom, x, x_i1);
487 for (i = 0; i < isize0; i++)
491 /* Special loop, since we need to determine the minimum distance
492 * over all selected atoms in the reference molecule/residue.
494 if (rdft[0][0] == 'a')
496 isize_g = isize[g+1];
502 for (j = 0; j < isize_g; j++)
505 /* Loop over the selected atoms in the reference molecule */
506 for (ii = coi[0][i]; ii < coi[0][i+1]; ii++)
510 pbc_dx(&pbc, x[index[0][ii]], x_i1[j], dx);
514 rvec_sub(x[index[0][ii]], x_i1[j], dx);
518 r2ii = dx[XX]*dx[XX] + dx[YY]*dx[YY];
522 r2ii = iprod(dx, dx);
529 if (r2 > cut2 && r2 <= rmax2)
531 count[g][(int)(sqrt(r2)*invhbinw)]++;
537 /* Real rdf between points in space */
538 if (bCM || rdft[0][0] != 'a')
540 copy_rvec(x0[i], xi);
544 copy_rvec(x[index[0][i]], xi);
546 if (rdft[0][0] == 'a' && npairs[g][i] >= 0)
548 /* Expensive loop, because of indexing */
549 for (j = 0; j < npairs[g][i]; j++)
554 pbc_dx(&pbc, xi, x[jx], dx);
558 rvec_sub(xi, x[jx], dx);
563 r2 = dx[XX]*dx[XX] + dx[YY]*dx[YY];
569 if (r2 > cut2 && r2 <= rmax2)
571 count[g][(int)(sqrt(r2)*invhbinw)]++;
577 /* Cheaper loop, no exclusions */
578 if (rdft[0][0] == 'a')
580 isize_g = isize[g+1];
586 for (j = 0; j < isize_g; j++)
590 pbc_dx(&pbc, xi, x_i1[j], dx);
594 rvec_sub(xi, x_i1[j], dx);
598 r2 = dx[XX]*dx[XX] + dx[YY]*dx[YY];
604 if (r2 > cut2 && r2 <= rmax2)
606 count[g][(int)(sqrt(r2)*invhbinw)]++;
615 while (read_next_x(oenv, status, &t, x, box));
616 fprintf(stderr, "\n");
618 if (bPBC && (NULL != top))
620 gmx_rmpbc_done(gpbc);
628 invvol = invvol_sum/nframes;
630 /* Calculate volume of sphere segments or length of circle segments */
631 snew(inv_segvol, (nbin+1)/2);
633 for (i = 0; (i < (nbin+1)/2); i++)
635 r = (i + 0.5)*binwidth;
638 spherevol = M_PI*r*r;
642 spherevol = (4.0/3.0)*M_PI*r*r*r;
644 segvol = spherevol-prev_spherevol;
645 inv_segvol[i] = 1.0/segvol;
646 prev_spherevol = spherevol;
650 for (g = 0; g < ng; g++)
652 /* We have to normalize by dividing by the number of frames */
653 if (rdft[0][0] == 'a')
655 normfac = 1.0/(nframes*invvol*isize0*isize[g+1]);
659 normfac = 1.0/(nframes*invvol*isize0*is[g+1]);
662 /* Do the normalization */
663 nrdf = max((nbin+1)/2, 1+2*fade/binwidth);
665 for (i = 0; i < (nbin+1)/2; i++)
674 j = count[g][i*2-1] + count[g][i*2];
676 if ((fade > 0) && (r >= fade))
678 rdf[g][i] = 1 + (j*inv_segvol[i]*normfac-1)*exp(-16*sqr(r/fade-1));
684 rdf[g][i] = j*inv_segvol[i]*normfac;
688 rdf[g][i] = j/(binwidth*isize0*nframes);
692 for (; (i < nrdf); i++)
698 if (rdft[0][0] == 'a')
700 sprintf(gtitle, "Radial distribution");
704 sprintf(gtitle, "Radial distribution of %s %s",
705 rdft[0][0] == 'm' ? "molecule" : "residue",
706 rdft[0][6] == 'm' ? "COM" : "COG");
708 fp = xvgropen(fnRDF, gtitle, "r", "", oenv);
711 sprintf(refgt, " %s", "COM");
715 sprintf(refgt, " closest atom in %s.", close);
719 sprintf(refgt, "%s", "");
723 if (output_env_get_print_xvgr_codes(oenv))
725 fprintf(fp, "@ subtitle \"%s%s - %s\"\n", grpname[0], refgt, grpname[1]);
730 if (output_env_get_print_xvgr_codes(oenv))
732 fprintf(fp, "@ subtitle \"reference %s%s\"\n", grpname[0], refgt);
734 xvgr_legend(fp, ng, (const char**)(grpname+1), oenv);
736 for (i = 0; (i < nrdf); i++)
738 fprintf(fp, "%10g", i*binwidth);
739 for (g = 0; g < ng; g++)
741 fprintf(fp, " %10g", rdf[g][i]);
747 do_view(oenv, fnRDF, NULL);
749 /* h(Q) function: fourier transform of rdf */
753 real *hq, *integrand, Q;
755 /* Get a better number density later! */
756 rho = isize[1]*invvol;
758 snew(integrand, nrdf);
759 for (i = 0; (i < nhq); i++)
763 for (j = 1; (j < nrdf); j++)
766 integrand[j] = (Q == 0) ? 1.0 : sin(Q*r)/(Q*r);
767 integrand[j] *= 4.0*M_PI*rho*r*r*(rdf[0][j]-1.0);
769 hq[i] = print_and_integrate(debug, nrdf, binwidth, integrand, NULL, 0);
771 fp = xvgropen(fnHQ, "h(Q)", "Q(/nm)", "h(Q)", oenv);
772 for (i = 0; (i < nhq); i++)
774 fprintf(fp, "%10g %10g\n", i*0.5, hq[i]);
777 do_view(oenv, fnHQ, NULL);
784 normfac = 1.0/(isize0*nframes);
785 fp = xvgropen(fnCNRDF, "Cumulative Number RDF", "r", "number", oenv);
788 if (output_env_get_print_xvgr_codes(oenv))
790 fprintf(fp, "@ subtitle \"%s-%s\"\n", grpname[0], grpname[1]);
795 if (output_env_get_print_xvgr_codes(oenv))
797 fprintf(fp, "@ subtitle \"reference %s\"\n", grpname[0]);
799 xvgr_legend(fp, ng, (const char**)(grpname+1), oenv);
802 for (i = 0; (i <= nbin/2); i++)
804 fprintf(fp, "%10g", i*binwidth);
805 for (g = 0; g < ng; g++)
807 fprintf(fp, " %10g", (real)((double)sum[g]*normfac));
810 sum[g] += count[g][i*2] + count[g][i*2+1];
818 do_view(oenv, fnCNRDF, NULL);
821 for (g = 0; g < ng; g++)
829 int gmx_rdf(int argc, char *argv[])
831 const char *desc[] = {
832 "The structure of liquids can be studied by either neutron or X-ray",
833 "scattering. The most common way to describe liquid structure is by a",
834 "radial distribution function. However, this is not easy to obtain from",
835 "a scattering experiment.[PAR]",
836 "[THISMODULE] calculates radial distribution functions in different ways.",
837 "The normal method is around a (set of) particle(s), the other methods",
838 "are around the center of mass of a set of particles ([TT]-com[tt])",
839 "or to the closest particle in a set ([TT]-surf[tt]).",
840 "With all methods, the RDF can also be calculated around axes parallel",
841 "to the [IT]z[it]-axis with option [TT]-xy[tt].",
842 "With option [TT]-surf[tt] normalization can not be used.[PAR]",
843 "The option [TT]-rdf[tt] sets the type of RDF to be computed.",
844 "Default is for atoms or particles, but one can also select center",
845 "of mass or geometry of molecules or residues. In all cases, only",
846 "the atoms in the index groups are taken into account.",
847 "For molecules and/or the center of mass option, a run input file",
849 "Weighting other than COM or COG can currently only be achieved",
850 "by providing a run input file with different masses.",
851 "Options [TT]-com[tt] and [TT]-surf[tt] also work in conjunction",
852 "with [TT]-rdf[tt].[PAR]",
853 "If a run input file is supplied ([TT]-s[tt]) and [TT]-rdf[tt] is set",
854 "to [TT]atom[tt], exclusions defined",
855 "in that file are taken into account when calculating the RDF.",
856 "The option [TT]-cut[tt] is meant as an alternative way to avoid",
857 "intramolecular peaks in the RDF plot.",
858 "It is however better to supply a run input file with a higher number of",
859 "exclusions. For e.g. benzene a topology, setting nrexcl to 5",
860 "would eliminate all intramolecular contributions to the RDF.",
861 "Note that all atoms in the selected groups are used, also the ones",
862 "that don't have Lennard-Jones interactions.[PAR]",
863 "Option [TT]-cn[tt] produces the cumulative number RDF,",
864 "i.e. the average number of particles within a distance r.[PAR]"
866 static gmx_bool bCM = FALSE, bXY = FALSE, bPBC = TRUE, bNormalize = TRUE;
867 static real cutoff = 0, binwidth = 0.002, fade = 0.0;
868 static int ngroups = 1;
870 static const char *closet[] = { NULL, "no", "mol", "res", NULL };
871 static const char *rdft[] = { NULL, "atom", "mol_com", "mol_cog", "res_com", "res_cog", NULL };
874 { "-bin", FALSE, etREAL, {&binwidth},
876 { "-com", FALSE, etBOOL, {&bCM},
877 "RDF with respect to the center of mass of first group" },
878 { "-surf", FALSE, etENUM, {closet},
879 "RDF with respect to the surface of the first group" },
880 { "-rdf", FALSE, etENUM, {rdft},
882 { "-pbc", FALSE, etBOOL, {&bPBC},
883 "Use periodic boundary conditions for computing distances. Without PBC the maximum range will be three times the largest box edge." },
884 { "-norm", FALSE, etBOOL, {&bNormalize},
885 "Normalize for volume and density" },
886 { "-xy", FALSE, etBOOL, {&bXY},
887 "Use only the x and y components of the distance" },
888 { "-cut", FALSE, etREAL, {&cutoff},
889 "Shortest distance (nm) to be considered"},
890 { "-ng", FALSE, etINT, {&ngroups},
891 "Number of secondary groups to compute RDFs around a central group" },
892 { "-fade", FALSE, etREAL, {&fade},
893 "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." }
895 #define NPA asize(pa)
896 const char *fnTPS, *fnNDX;
900 { efTRX, "-f", NULL, ffREAD },
901 { efTPS, NULL, NULL, ffOPTRD },
902 { efNDX, NULL, NULL, ffOPTRD },
903 { efXVG, "-o", "rdf", ffWRITE },
904 { efXVG, "-cn", "rdf_cn", ffOPTWR },
905 { efXVG, "-hq", "hq", ffOPTWR },
907 #define NFILE asize(fnm)
908 if (!parse_common_args(&argc, argv, PCA_CAN_VIEW | PCA_CAN_TIME | PCA_BE_NICE,
909 NFILE, fnm, NPA, pa, asize(desc), desc, 0, NULL, &oenv))
914 if (bCM || closet[0][0] != 'n' || rdft[0][0] == 'm' || rdft[0][6] == 'm')
916 fnTPS = ftp2fn(efTPS, NFILE, fnm);
920 fnTPS = ftp2fn_null(efTPS, NFILE, fnm);
922 fnNDX = ftp2fn_null(efNDX, NFILE, fnm);
924 if (!fnTPS && !fnNDX)
926 gmx_fatal(FARGS, "Neither index file nor topology file specified\n"
930 if (closet[0][0] != 'n')
934 gmx_fatal(FARGS, "Can not have both -com and -surf");
938 fprintf(stderr, "Turning of normalization because of option -surf\n");
943 do_rdf(fnNDX, fnTPS, ftp2fn(efTRX, NFILE, fnm),
944 opt2fn("-o", NFILE, fnm), opt2fn_null("-cn", NFILE, fnm),
945 opt2fn_null("-hq", NFILE, fnm),
946 bCM, closet[0], rdft, bXY, bPBC, bNormalize, cutoff, binwidth, fade, ngroups,