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.
42 #include "gromacs/legacyheaders/typedefs.h"
43 #include "gromacs/legacyheaders/macros.h"
44 #include "gromacs/math/vec.h"
45 #include "gromacs/pbcutil/pbc.h"
46 #include "gromacs/legacyheaders/viewit.h"
47 #include "gromacs/utility/futil.h"
48 #include "gromacs/fileio/tpxio.h"
49 #include "gromacs/fileio/trxio.h"
50 #include "gromacs/topology/index.h"
51 #include "gromacs/legacyheaders/nrnb.h"
52 #include "gromacs/legacyheaders/coulomb.h"
55 #include "gromacs/legacyheaders/names.h"
57 #include "gromacs/commandline/pargs.h"
58 #include "gromacs/fileio/xvgr.h"
59 #include "gromacs/pbcutil/rmpbc.h"
60 #include "gromacs/utility/fatalerror.h"
61 #include "gromacs/utility/smalloc.h"
63 static void check_box_c(matrix box)
65 if (fabs(box[ZZ][XX]) > GMX_REAL_EPS*box[ZZ][ZZ] ||
66 fabs(box[ZZ][YY]) > GMX_REAL_EPS*box[ZZ][ZZ])
69 "The last box vector is not parallel to the z-axis: %f %f %f",
70 box[ZZ][XX], box[ZZ][YY], box[ZZ][ZZ]);
74 static void calc_comg(int is, int *coi, int *index, gmx_bool bMass, t_atom *atom,
75 rvec *x, rvec *x_comg)
81 if (bMass && atom == NULL)
83 gmx_fatal(FARGS, "No masses available while mass weighting was requested");
86 for (c = 0; c < is; c++)
90 for (i = coi[c]; i < coi[c+1]; i++)
95 for (d = 0; d < DIM; d++)
97 xc[d] += m*x[index[i]][d];
103 rvec_inc(xc, x[index[i]]);
107 svmul(1/mtot, xc, x_comg[c]);
111 static void split_group(int isize, int *index, char *grpname,
112 t_topology *top, char type,
113 int *is_out, int **coi_out)
115 t_block *mols = NULL;
118 int cur, mol, res, i, a, i1;
120 /* Split up the group in molecules or residues */
127 atom = top->atoms.atom;
130 gmx_fatal(FARGS, "Unknown rdf option '%s'", type);
136 for (i = 0; i < isize; i++)
141 /* Check if the molecule number has changed */
142 i1 = mols->index[mol+1];
146 i1 = mols->index[mol+1];
154 else if (type == 'r')
156 /* Check if the residue index has changed */
157 res = atom[a].resind;
167 printf("Group '%s' of %d atoms consists of %d %s\n",
169 (type == 'm' ? "molecules" : "residues"));
175 static void do_rdf(const char *fnNDX, const char *fnTPS, const char *fnTRX,
176 const char *fnRDF, const char *fnCNRDF, const char *fnHQ,
177 gmx_bool bCM, const char *close,
178 const char **rdft, gmx_bool bXY, gmx_bool bPBC, gmx_bool bNormalize,
179 real cutoff, real binwidth, real fade, int ng,
180 const output_env_t oenv)
184 char outf1[STRLEN], outf2[STRLEN];
185 char title[STRLEN], gtitle[STRLEN], refgt[30];
186 int g, natoms, i, ii, j, k, nbin, j0, j1, n, nframes;
189 int *isize, isize_cm = 0, nrdf = 0, max_i, isize0, isize_g;
190 atom_id **index, *index_cm = NULL;
192 real t, rmax2, cut2, r, r2, r2ii, invhbinw, normfac;
193 real segvol, spherevol, prev_spherevol, **rdf;
194 rvec *x, dx, *x0 = NULL, *x_i1, xi;
195 real *inv_segvol, invvol, invvol_sum, rho;
196 gmx_bool bClose, *bExcl, bTop, bNonSelfExcl;
199 atom_id ix, jx, ***pairs;
200 t_topology *top = NULL;
201 int ePBC = -1, ePBCrdf = -1;
202 t_block *mols = NULL;
206 gmx_rmpbc_t gpbc = NULL;
207 int *is = NULL, **coi = NULL, cur, mol, i1, res, a;
211 bClose = (close[0] != 'n');
216 bTop = read_tps_conf(fnTPS, title, top, &ePBC, &x, NULL, box, TRUE);
217 if (bTop && !(bCM || bClose))
219 /* get exclusions from topology */
220 excl = &(top->excls);
226 fprintf(stderr, "\nSelect a reference group and %d group%s\n",
227 ng, ng == 1 ? "" : "s");
230 get_index(&(top->atoms), fnNDX, ng+1, isize, index, grpname);
231 atom = top->atoms.atom;
235 rd_index(fnNDX, ng+1, isize, index, grpname);
244 split_group(isize[0], index[0], grpname[0], top, close[0], &is[0], &coi[0]);
247 if (rdft[0][0] != 'a')
249 /* Split up all the groups in molecules or residues */
252 for (g = ((bCM || bClose) ? 1 : 0); g < ng+1; g++)
254 split_group(isize[g], index[g], grpname[g], top, rdft[0][0], &is[g], &coi[g]);
261 snew(coi[0], is[0]+1);
263 coi[0][1] = isize[0];
267 else if (bClose || rdft[0][0] != 'a')
277 natoms = read_first_x(oenv, &status, fnTRX, &t, &x, box);
280 gmx_fatal(FARGS, "Could not read coordinates from statusfile\n");
284 /* check with topology */
285 if (natoms > top->atoms.nr)
287 gmx_fatal(FARGS, "Trajectory (%d atoms) does not match topology (%d atoms)",
288 natoms, top->atoms.nr);
291 /* check with index groups */
292 for (i = 0; i < ng+1; i++)
294 for (j = 0; j < isize[i]; j++)
296 if (index[i][j] >= natoms)
298 gmx_fatal(FARGS, "Atom index (%d) in index group %s (%d atoms) larger "
299 "than number of atoms in trajectory (%d atoms)",
300 index[i][j], grpname[i], isize[i], natoms);
305 /* initialize some handy things */
308 ePBC = guess_ePBC(box);
310 copy_mat(box, box_pbc);
317 case epbcXY: ePBCrdf = epbcXY; break;
318 case epbcNONE: ePBCrdf = epbcNONE; break;
320 gmx_fatal(FARGS, "xy-rdf's are not supported for pbc type'%s'",
324 /* Make sure the z-height does not influence the cut-off */
325 box_pbc[ZZ][ZZ] = 2*max(box[XX][XX], box[YY][YY]);
333 rmax2 = 0.99*0.99*max_cutoff2(bXY ? epbcXY : epbcXYZ, box_pbc);
337 rmax2 = sqr(3*max(box[XX][XX], max(box[YY][YY], box[ZZ][ZZ])));
341 fprintf(debug, "rmax2 = %g\n", rmax2);
344 /* We use the double amount of bins, so we can correctly
345 * write the rdf and rdf_cn output at i*binwidth values.
347 nbin = (int)(sqrt(rmax2) * 2 / binwidth);
348 invhbinw = 2.0 / binwidth;
357 for (g = 0; g < ng; g++)
359 if (isize[g+1] > max_i)
364 /* this is THE array */
365 snew(count[g], nbin+1);
367 /* make pairlist array for groups and exclusions */
368 snew(pairs[g], isize[0]);
369 snew(npairs[g], isize[0]);
370 for (i = 0; i < isize[0]; i++)
372 /* We can only have exclusions with atomic rdfs */
373 if (!(bCM || bClose || rdft[0][0] != 'a'))
376 for (j = 0; j < natoms; j++)
383 for (j = excl->index[ix]; j < excl->index[ix+1]; j++)
385 bExcl[excl->a[j]] = TRUE;
389 snew(pairs[g][i], isize[g+1]);
390 bNonSelfExcl = FALSE;
391 for (j = 0; j < isize[g+1]; j++)
396 pairs[g][i][k++] = jx;
400 /* Check if we have exclusions other than self exclusions */
407 srenew(pairs[g][i], npairs[g][i]);
411 /* Save a LOT of memory and some cpu cycles */
427 if (bPBC && (NULL != top))
429 gpbc = gmx_rmpbc_init(&top->idef, ePBC, natoms);
434 /* Must init pbc every step because of pressure coupling */
435 copy_mat(box, box_pbc);
440 gmx_rmpbc(gpbc, natoms, box, x);
445 clear_rvec(box_pbc[ZZ]);
447 set_pbc(&pbc, ePBCrdf, box_pbc);
451 /* Set z-size to 1 so we get the surface iso the volume */
455 invvol = 1/det(box_pbc);
456 invvol_sum += invvol;
460 /* Calculate center of mass of the whole group */
461 calc_comg(is[0], coi[0], index[0], TRUE, atom, x, x0);
463 else if (!bClose && rdft[0][0] != 'a')
465 calc_comg(is[0], coi[0], index[0], rdft[0][6] == 'm', atom, x, x0);
468 for (g = 0; g < ng; g++)
470 if (rdft[0][0] == 'a')
472 /* Copy the indexed coordinates to a continuous array */
473 for (i = 0; i < isize[g+1]; i++)
475 copy_rvec(x[index[g+1][i]], x_i1[i]);
480 /* Calculate the COMs/COGs and store in x_i1 */
481 calc_comg(is[g+1], coi[g+1], index[g+1], rdft[0][6] == 'm', atom, x, x_i1);
484 for (i = 0; i < isize0; i++)
488 /* Special loop, since we need to determine the minimum distance
489 * over all selected atoms in the reference molecule/residue.
491 if (rdft[0][0] == 'a')
493 isize_g = isize[g+1];
499 for (j = 0; j < isize_g; j++)
502 /* Loop over the selected atoms in the reference molecule */
503 for (ii = coi[0][i]; ii < coi[0][i+1]; ii++)
507 pbc_dx(&pbc, x[index[0][ii]], x_i1[j], dx);
511 rvec_sub(x[index[0][ii]], x_i1[j], dx);
515 r2ii = dx[XX]*dx[XX] + dx[YY]*dx[YY];
519 r2ii = iprod(dx, dx);
526 if (r2 > cut2 && r2 <= rmax2)
528 count[g][(int)(sqrt(r2)*invhbinw)]++;
534 /* Real rdf between points in space */
535 if (bCM || rdft[0][0] != 'a')
537 copy_rvec(x0[i], xi);
541 copy_rvec(x[index[0][i]], xi);
543 if (rdft[0][0] == 'a' && npairs[g][i] >= 0)
545 /* Expensive loop, because of indexing */
546 for (j = 0; j < npairs[g][i]; j++)
551 pbc_dx(&pbc, xi, x[jx], dx);
555 rvec_sub(xi, x[jx], dx);
560 r2 = dx[XX]*dx[XX] + dx[YY]*dx[YY];
566 if (r2 > cut2 && r2 <= rmax2)
568 count[g][(int)(sqrt(r2)*invhbinw)]++;
574 /* Cheaper loop, no exclusions */
575 if (rdft[0][0] == 'a')
577 isize_g = isize[g+1];
583 for (j = 0; j < isize_g; j++)
587 pbc_dx(&pbc, xi, x_i1[j], dx);
591 rvec_sub(xi, x_i1[j], dx);
595 r2 = dx[XX]*dx[XX] + dx[YY]*dx[YY];
601 if (r2 > cut2 && r2 <= rmax2)
603 count[g][(int)(sqrt(r2)*invhbinw)]++;
612 while (read_next_x(oenv, status, &t, x, box));
613 fprintf(stderr, "\n");
615 if (bPBC && (NULL != top))
617 gmx_rmpbc_done(gpbc);
625 invvol = invvol_sum/nframes;
627 /* Calculate volume of sphere segments or length of circle segments */
628 snew(inv_segvol, (nbin+1)/2);
630 for (i = 0; (i < (nbin+1)/2); i++)
632 r = (i + 0.5)*binwidth;
635 spherevol = M_PI*r*r;
639 spherevol = (4.0/3.0)*M_PI*r*r*r;
641 segvol = spherevol-prev_spherevol;
642 inv_segvol[i] = 1.0/segvol;
643 prev_spherevol = spherevol;
647 for (g = 0; g < ng; g++)
649 /* We have to normalize by dividing by the number of frames */
650 if (rdft[0][0] == 'a')
652 normfac = 1.0/(nframes*invvol*isize0*isize[g+1]);
656 normfac = 1.0/(nframes*invvol*isize0*is[g+1]);
659 /* Do the normalization */
660 nrdf = max((nbin+1)/2, 1+2*fade/binwidth);
662 for (i = 0; i < (nbin+1)/2; i++)
671 j = count[g][i*2-1] + count[g][i*2];
673 if ((fade > 0) && (r >= fade))
675 rdf[g][i] = 1 + (j*inv_segvol[i]*normfac-1)*exp(-16*sqr(r/fade-1));
681 rdf[g][i] = j*inv_segvol[i]*normfac;
685 rdf[g][i] = j/(binwidth*isize0*nframes);
689 for (; (i < nrdf); i++)
695 if (rdft[0][0] == 'a')
697 sprintf(gtitle, "Radial distribution");
701 sprintf(gtitle, "Radial distribution of %s %s",
702 rdft[0][0] == 'm' ? "molecule" : "residue",
703 rdft[0][6] == 'm' ? "COM" : "COG");
705 fp = xvgropen(fnRDF, gtitle, "r", "", oenv);
708 sprintf(refgt, " %s", "COM");
712 sprintf(refgt, " closest atom in %s.", close);
716 sprintf(refgt, "%s", "");
720 if (output_env_get_print_xvgr_codes(oenv))
722 fprintf(fp, "@ subtitle \"%s%s - %s\"\n", grpname[0], refgt, grpname[1]);
727 if (output_env_get_print_xvgr_codes(oenv))
729 fprintf(fp, "@ subtitle \"reference %s%s\"\n", grpname[0], refgt);
731 xvgr_legend(fp, ng, (const char**)(grpname+1), oenv);
733 for (i = 0; (i < nrdf); i++)
735 fprintf(fp, "%10g", i*binwidth);
736 for (g = 0; g < ng; g++)
738 fprintf(fp, " %10g", rdf[g][i]);
744 do_view(oenv, fnRDF, NULL);
746 /* h(Q) function: fourier transform of rdf */
750 real *hq, *integrand, Q;
752 /* Get a better number density later! */
753 rho = isize[1]*invvol;
755 snew(integrand, nrdf);
756 for (i = 0; (i < nhq); i++)
760 for (j = 1; (j < nrdf); j++)
763 integrand[j] = (Q == 0) ? 1.0 : sin(Q*r)/(Q*r);
764 integrand[j] *= 4.0*M_PI*rho*r*r*(rdf[0][j]-1.0);
766 hq[i] = print_and_integrate(debug, nrdf, binwidth, integrand, NULL, 0);
768 fp = xvgropen(fnHQ, "h(Q)", "Q(/nm)", "h(Q)", oenv);
769 for (i = 0; (i < nhq); i++)
771 fprintf(fp, "%10g %10g\n", i*0.5, hq[i]);
774 do_view(oenv, fnHQ, NULL);
781 normfac = 1.0/(isize0*nframes);
782 fp = xvgropen(fnCNRDF, "Cumulative Number RDF", "r", "number", oenv);
785 if (output_env_get_print_xvgr_codes(oenv))
787 fprintf(fp, "@ subtitle \"%s-%s\"\n", grpname[0], grpname[1]);
792 if (output_env_get_print_xvgr_codes(oenv))
794 fprintf(fp, "@ subtitle \"reference %s\"\n", grpname[0]);
796 xvgr_legend(fp, ng, (const char**)(grpname+1), oenv);
799 for (i = 0; (i <= nbin/2); i++)
801 fprintf(fp, "%10g", i*binwidth);
802 for (g = 0; g < ng; g++)
804 fprintf(fp, " %10g", (real)((double)sum[g]*normfac));
807 sum[g] += count[g][i*2] + count[g][i*2+1];
815 do_view(oenv, fnCNRDF, NULL);
818 for (g = 0; g < ng; g++)
826 int gmx_rdf(int argc, char *argv[])
828 const char *desc[] = {
829 "The structure of liquids can be studied by either neutron or X-ray",
830 "scattering. The most common way to describe liquid structure is by a",
831 "radial distribution function. However, this is not easy to obtain from",
832 "a scattering experiment.[PAR]",
833 "[THISMODULE] calculates radial distribution functions in different ways.",
834 "The normal method is around a (set of) particle(s), the other methods",
835 "are around the center of mass of a set of particles ([TT]-com[tt])",
836 "or to the closest particle in a set ([TT]-surf[tt]).",
837 "With all methods, the RDF can also be calculated around axes parallel",
838 "to the [IT]z[it]-axis with option [TT]-xy[tt].",
839 "With option [TT]-surf[tt] normalization can not be used.[PAR]",
840 "The option [TT]-rdf[tt] sets the type of RDF to be computed.",
841 "Default is for atoms or particles, but one can also select center",
842 "of mass or geometry of molecules or residues. In all cases, only",
843 "the atoms in the index groups are taken into account.",
844 "For molecules and/or the center of mass option, a run input file",
846 "Weighting other than COM or COG can currently only be achieved",
847 "by providing a run input file with different masses.",
848 "Options [TT]-com[tt] and [TT]-surf[tt] also work in conjunction",
849 "with [TT]-rdf[tt].[PAR]",
850 "If a run input file is supplied ([TT]-s[tt]) and [TT]-rdf[tt] is set",
851 "to [TT]atom[tt], exclusions defined",
852 "in that file are taken into account when calculating the RDF.",
853 "The option [TT]-cut[tt] is meant as an alternative way to avoid",
854 "intramolecular peaks in the RDF plot.",
855 "It is however better to supply a run input file with a higher number of",
856 "exclusions. For e.g. benzene a topology, setting nrexcl to 5",
857 "would eliminate all intramolecular contributions to the RDF.",
858 "Note that all atoms in the selected groups are used, also the ones",
859 "that don't have Lennard-Jones interactions.[PAR]",
860 "Option [TT]-cn[tt] produces the cumulative number RDF,",
861 "i.e. the average number of particles within a distance r.[PAR]"
863 static gmx_bool bCM = FALSE, bXY = FALSE, bPBC = TRUE, bNormalize = TRUE;
864 static real cutoff = 0, binwidth = 0.002, fade = 0.0;
865 static int ngroups = 1;
867 static const char *closet[] = { NULL, "no", "mol", "res", NULL };
868 static const char *rdft[] = { NULL, "atom", "mol_com", "mol_cog", "res_com", "res_cog", NULL };
871 { "-bin", FALSE, etREAL, {&binwidth},
873 { "-com", FALSE, etBOOL, {&bCM},
874 "RDF with respect to the center of mass of first group" },
875 { "-surf", FALSE, etENUM, {closet},
876 "RDF with respect to the surface of the first group" },
877 { "-rdf", FALSE, etENUM, {rdft},
879 { "-pbc", FALSE, etBOOL, {&bPBC},
880 "Use periodic boundary conditions for computing distances. Without PBC the maximum range will be three times the largest box edge." },
881 { "-norm", FALSE, etBOOL, {&bNormalize},
882 "Normalize for volume and density" },
883 { "-xy", FALSE, etBOOL, {&bXY},
884 "Use only the x and y components of the distance" },
885 { "-cut", FALSE, etREAL, {&cutoff},
886 "Shortest distance (nm) to be considered"},
887 { "-ng", FALSE, etINT, {&ngroups},
888 "Number of secondary groups to compute RDFs around a central group" },
889 { "-fade", FALSE, etREAL, {&fade},
890 "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." }
892 #define NPA asize(pa)
893 const char *fnTPS, *fnNDX;
897 { efTRX, "-f", NULL, ffREAD },
898 { efTPS, NULL, NULL, ffOPTRD },
899 { efNDX, NULL, NULL, ffOPTRD },
900 { efXVG, "-o", "rdf", ffWRITE },
901 { efXVG, "-cn", "rdf_cn", ffOPTWR },
902 { efXVG, "-hq", "hq", ffOPTWR },
904 #define NFILE asize(fnm)
905 if (!parse_common_args(&argc, argv, PCA_CAN_VIEW | PCA_CAN_TIME,
906 NFILE, fnm, NPA, pa, asize(desc), desc, 0, NULL, &oenv))
911 if (bCM || closet[0][0] != 'n' || rdft[0][0] == 'm' || rdft[0][6] == 'm')
913 fnTPS = ftp2fn(efTPS, NFILE, fnm);
917 fnTPS = ftp2fn_null(efTPS, NFILE, fnm);
919 fnNDX = ftp2fn_null(efNDX, NFILE, fnm);
921 if (!fnTPS && !fnNDX)
923 gmx_fatal(FARGS, "Neither index file nor topology file specified\n"
927 if (closet[0][0] != 'n')
931 gmx_fatal(FARGS, "Can not have both -com and -surf");
935 fprintf(stderr, "Turning of normalization because of option -surf\n");
940 do_rdf(fnNDX, fnTPS, ftp2fn(efTRX, NFILE, fnm),
941 opt2fn("-o", NFILE, fnm), opt2fn_null("-cn", NFILE, fnm),
942 opt2fn_null("-hq", NFILE, fnm),
943 bCM, closet[0], rdft, bXY, bPBC, bNormalize, cutoff, binwidth, fade, ngroups,