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.
45 #include "gromacs/commandline/pargs.h"
46 #include "gromacs/legacyheaders/typedefs.h"
47 #include "gromacs/utility/smalloc.h"
48 #include "gromacs/legacyheaders/macros.h"
49 #include "gromacs/math/vec.h"
50 #include "gromacs/pbcutil/pbc.h"
51 #include "gromacs/utility/futil.h"
52 #include "gromacs/topology/index.h"
53 #include "gromacs/fileio/xvgr.h"
55 #include "gromacs/pbcutil/rmpbc.h"
56 #include "gromacs/legacyheaders/txtdump.h"
57 #include "gromacs/fileio/tpxio.h"
58 #include "gromacs/fileio/trxio.h"
61 #include "gromacs/utility/fatalerror.h"
70 static int ocomp(const void *a, const void *b)
87 int gmx_trjorder(int argc, char *argv[])
89 const char *desc[] = {
90 "[THISMODULE] orders molecules according to the smallest distance",
91 "to atoms in a reference group",
92 "or on z-coordinate (with option [TT]-z[tt]).",
93 "With distance ordering, it will ask for a group of reference",
94 "atoms and a group of molecules. For each frame of the trajectory",
95 "the selected molecules will be reordered according to the shortest",
96 "distance between atom number [TT]-da[tt] in the molecule and all the",
97 "atoms in the reference group. The center of mass of the molecules can",
98 "be used instead of a reference atom by setting [TT]-da[tt] to 0.",
99 "All atoms in the trajectory are written",
100 "to the output trajectory.[PAR]",
101 "[THISMODULE] can be useful for e.g. analyzing the n waters closest to a",
103 "In that case the reference group would be the protein and the group",
104 "of molecules would consist of all the water atoms. When an index group",
105 "of the first n waters is made, the ordered trajectory can be used",
106 "with any Gromacs program to analyze the n closest waters.",
108 "If the output file is a [TT].pdb[tt] file, the distance to the reference target",
109 "will be stored in the B-factor field in order to color with e.g. Rasmol.",
111 "With option [TT]-nshell[tt] the number of molecules within a shell",
112 "of radius [TT]-r[tt] around the reference group are printed."
114 static int na = 3, ref_a = 1;
115 static real rcut = 0;
116 static gmx_bool bCOM = FALSE, bZ = FALSE;
118 { "-na", FALSE, etINT, {&na},
119 "Number of atoms in a molecule" },
120 { "-da", FALSE, etINT, {&ref_a},
121 "Atom used for the distance calculation, 0 is COM" },
122 { "-com", FALSE, etBOOL, {&bCOM},
123 "Use the distance to the center of mass of the reference group" },
124 { "-r", FALSE, etREAL, {&rcut},
125 "Cutoff used for the distance calculation when computing the number of molecules in a shell around e.g. a protein" },
126 { "-z", FALSE, etBOOL, {&bZ},
127 "Order molecules on z-coordinate" }
132 gmx_bool bNShell, bPDBout;
135 rvec *x, *xsol, xcom, dx;
139 real t, totmass, mass, rcut2 = 0, n2;
140 int natoms, nwat, ncut;
141 char **grpname, title[256];
142 int i, j, d, *isize, isize_ref = 0, isize_sol;
143 atom_id sa, sr, *swi, **index, *ind_ref = NULL, *ind_sol;
146 { efTRX, "-f", NULL, ffREAD },
147 { efTPS, NULL, NULL, ffREAD },
148 { efNDX, NULL, NULL, ffOPTRD },
149 { efTRO, "-o", "ordered", ffOPTWR },
150 { efXVG, "-nshell", "nshell", ffOPTWR }
152 #define NFILE asize(fnm)
154 if (!parse_common_args(&argc, argv, PCA_CAN_TIME | PCA_BE_NICE,
155 NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, NULL, &oenv))
160 read_tps_conf(ftp2fn(efTPS, NFILE, fnm), title, &top, &ePBC, &x, NULL, box, TRUE);
163 /* get index groups */
164 printf("Select %sa group of molecules to be ordered:\n",
165 bZ ? "" : "a group of reference atoms and ");
169 get_index(&top.atoms, ftp2fn_null(efNDX, NFILE, fnm), bZ ? 1 : 2,
170 isize, index, grpname);
174 isize_ref = isize[0];
175 isize_sol = isize[1];
181 isize_sol = isize[0];
185 natoms = read_first_x(oenv, &status, ftp2fn(efTRX, NFILE, fnm), &t, &x, box);
186 if (natoms > top.atoms.nr)
188 gmx_fatal(FARGS, "Number of atoms in the run input file is larger than in the trjactory");
190 for (i = 0; (i < 2); i++)
192 for (j = 0; (j < isize[i]); j++)
194 if (index[i][j] > natoms)
196 gmx_fatal(FARGS, "An atom number in group %s is larger than the number of atoms in the trajectory");
201 if ((isize_sol % na) != 0)
203 gmx_fatal(FARGS, "Number of atoms in the molecule group (%d) is not a multiple of na (%d)",
210 gmx_fatal(FARGS, "The reference atom can not be larger than the number of atoms in a molecule");
216 for (i = 0; (i < natoms); i++)
223 bNShell = ((opt2bSet("-nshell", NFILE, fnm)) ||
224 (opt2parg_bSet("-r", asize(pa), pa)));
229 fp = xvgropen(opt2fn("-nshell", NFILE, fnm), "Number of molecules",
230 "Time (ps)", "N", oenv);
231 printf("Will compute the number of molecules within a radius of %g\n",
234 if (!bNShell || opt2bSet("-o", NFILE, fnm))
236 bPDBout = (fn2ftp(opt2fn("-o", NFILE, fnm)) == efPDB);
237 if (bPDBout && !top.atoms.pdbinfo)
239 fprintf(stderr, "Creating pdbfino records\n");
240 snew(top.atoms.pdbinfo, top.atoms.nr);
242 out = open_trx(opt2fn("-o", NFILE, fnm), "w");
244 gpbc = gmx_rmpbc_init(&top.idef, ePBC, natoms);
247 gmx_rmpbc(gpbc, natoms, box, x);
248 set_pbc(&pbc, ePBC, box);
252 /* Calculate the COM of all solvent molecules */
253 for (i = 0; i < nwat; i++)
257 for (j = 0; j < na; j++)
259 sa = ind_sol[i*na+j];
260 mass = top.atoms.atom[sa].m;
262 for (d = 0; d < DIM; d++)
264 xsol[i][d] += mass*x[sa][d];
267 svmul(1/totmass, xsol[i], xsol[i]);
272 /* Copy the reference atom of all solvent molecules */
273 for (i = 0; i < nwat; i++)
275 copy_rvec(x[ind_sol[i*na+ref_a]], xsol[i]);
281 for (i = 0; (i < nwat); i++)
285 order[i].d2 = xsol[i][ZZ];
292 for (i = 0; i < isize_ref; i++)
294 mass = top.atoms.atom[ind_ref[i]].m;
296 for (j = 0; j < DIM; j++)
298 xcom[j] += mass*x[ind_ref[i]][j];
301 svmul(1/totmass, xcom, xcom);
302 for (i = 0; (i < nwat); i++)
305 pbc_dx(&pbc, xcom, xsol[i], dx);
307 order[i].d2 = norm2(dx);
312 /* Set distance to first atom */
313 for (i = 0; (i < nwat); i++)
316 pbc_dx(&pbc, x[ind_ref[0]], xsol[i], dx);
318 order[i].d2 = norm2(dx);
320 for (j = 1; (j < isize_ref); j++)
323 for (i = 0; (i < nwat); i++)
326 pbc_dx(&pbc, x[sr], xsol[i], dx);
328 if (n2 < order[i].d2)
339 for (i = 0; (i < nwat); i++)
341 if (order[i].d2 <= rcut2)
346 fprintf(fp, "%10.3f %8d\n", t, ncut);
350 qsort(order, nwat, sizeof(*order), ocomp);
351 for (i = 0; (i < nwat); i++)
353 for (j = 0; (j < na); j++)
355 swi[ind_sol[na*i]+j] = order[i].i+j;
359 /* Store the distance as the B-factor */
362 for (i = 0; (i < nwat); i++)
364 for (j = 0; (j < na); j++)
366 top.atoms.pdbinfo[order[i].i+j].bfac = sqrt(order[i].d2);
370 write_trx(out, natoms, swi, &top.atoms, 0, t, box, x, NULL, NULL);
373 while (read_next_x(oenv, status, &t, x, box));
383 gmx_rmpbc_done(gpbc);