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.
43 #include "gromacs/math/vec.h"
45 #include "gromacs/pbcutil/pbc.h"
46 #include "gromacs/topology/index.h"
47 #include "gromacs/fileio/tpxio.h"
48 #include "gromacs/fileio/trxio.h"
51 #include "gromacs/commandline/pargs.h"
52 #include "gromacs/fileio/xvgr.h"
53 #include "gromacs/pbcutil/rmpbc.h"
54 #include "gromacs/utility/fatalerror.h"
55 #include "gromacs/utility/smalloc.h"
57 static void calc_com_pbc(int nrefat, t_topology *top, rvec x[], t_pbc *pbc,
58 atom_id index[], rvec xref, gmx_bool bPBC)
60 const real tol = 1e-4;
66 /* First simple calculation */
69 for (m = 0; (m < nrefat); m++)
72 mass = top->atoms.atom[ai].m;
73 for (j = 0; (j < DIM); j++)
75 xref[j] += mass*x[ai][j];
79 svmul(1/mtot, xref, xref);
80 /* Now check if any atom is more than half the box from the COM */
87 for (m = 0; (m < nrefat); m++)
90 mass = top->atoms.atom[ai].m/mtot;
91 pbc_dx(pbc, x[ai], xref, dx);
92 rvec_add(xref, dx, xtest);
93 for (j = 0; (j < DIM); j++)
95 if (fabs(xtest[j]-x[ai][j]) > tol)
97 /* Here we have used the wrong image for contributing to the COM */
98 xref[j] += mass*(xtest[j]-x[ai][j]);
106 printf("COM: %8.3f %8.3f %8.3f iter = %d\n", xref[XX], xref[YY], xref[ZZ], iter);
114 int gmx_sorient(int argc, char *argv[])
126 int i, j, p, sa0, sa1, sa2, n, ntot, nf, m, *hist1, *hist2, *histn, nbin1, nbin2, nrbin;
127 real *histi1, *histi2, invbw, invrbw;
129 int *isize, nrefgrp, nrefat;
132 real inp, outp, two_pi, nav, normfac, rmin2, rmax2, rcut, rcut2, r2, r, mass, mtot;
136 rvec xref, dx, dxh1, dxh2, outer;
137 gmx_rmpbc_t gpbc = NULL;
139 const char *legr[] = {
140 "<cos(\\8q\\4\\s1\\N)>",
141 "<3cos\\S2\\N(\\8q\\4\\s2\\N)-1>"
143 const char *legc[] = {
144 "cos(\\8q\\4\\s1\\N)",
145 "3cos\\S2\\N(\\8q\\4\\s2\\N)-1"
148 const char *desc[] = {
149 "[THISMODULE] analyzes solvent orientation around solutes.",
150 "It calculates two angles between the vector from one or more",
151 "reference positions to the first atom of each solvent molecule:[PAR]",
152 "[GRK]theta[grk][SUB]1[sub]: the angle with the vector from the first atom of the solvent",
153 "molecule to the midpoint between atoms 2 and 3.[BR]",
154 "[GRK]theta[grk][SUB]2[sub]: the angle with the normal of the solvent plane, defined by the",
155 "same three atoms, or, when the option [TT]-v23[tt] is set, ",
156 "the angle with the vector between atoms 2 and 3.[PAR]",
157 "The reference can be a set of atoms or",
158 "the center of mass of a set of atoms. The group of solvent atoms should",
159 "consist of 3 atoms per solvent molecule.",
160 "Only solvent molecules between [TT]-rmin[tt] and [TT]-rmax[tt] are",
161 "considered for [TT]-o[tt] and [TT]-no[tt] each frame.[PAR]",
162 "[TT]-o[tt]: distribtion of [MATH][COS][GRK]theta[grk][SUB]1[sub][cos][math] for rmin<=r<=rmax.[PAR]",
163 "[TT]-no[tt]: distribution of [MATH][COS][GRK]theta[grk][SUB]2[sub][cos][math] for rmin<=r<=rmax.[PAR]",
164 "[TT]-ro[tt]: [MATH][CHEVRON][COS][GRK]theta[grk][SUB]1[sub][cos][chevron][math] and [MATH][CHEVRON]3[COS]^2[GRK]theta[grk][SUB]2[sub][cos]-1[chevron][math] as a function of the",
166 "[TT]-co[tt]: the sum over all solvent molecules within distance r",
167 "of [MATH][COS][GRK]theta[grk][SUB]1[sub][cos][math] and [MATH]3[COS]^2([GRK]theta[grk][SUB]2[sub])-1[cos][math] as a function of r.[PAR]",
168 "[TT]-rc[tt]: the distribution of the solvent molecules as a function of r"
172 static gmx_bool bCom = FALSE, bVec23 = FALSE, bPBC = FALSE;
173 static real rmin = 0.0, rmax = 0.5, binwidth = 0.02, rbinw = 0.02;
175 { "-com", FALSE, etBOOL, {&bCom},
176 "Use the center of mass as the reference postion" },
177 { "-v23", FALSE, etBOOL, {&bVec23},
178 "Use the vector between atoms 2 and 3" },
179 { "-rmin", FALSE, etREAL, {&rmin}, "Minimum distance (nm)" },
180 { "-rmax", FALSE, etREAL, {&rmax}, "Maximum distance (nm)" },
181 { "-cbin", FALSE, etREAL, {&binwidth}, "Binwidth for the cosine" },
182 { "-rbin", FALSE, etREAL, {&rbinw}, "Binwidth for r (nm)" },
183 { "-pbc", FALSE, etBOOL, {&bPBC}, "Check PBC for the center of mass calculation. Only necessary when your reference group consists of several molecules." }
187 { efTRX, NULL, NULL, ffREAD },
188 { efTPS, NULL, NULL, ffREAD },
189 { efNDX, NULL, NULL, ffOPTRD },
190 { efXVG, NULL, "sori.xvg", ffWRITE },
191 { efXVG, "-no", "snor.xvg", ffWRITE },
192 { efXVG, "-ro", "sord.xvg", ffWRITE },
193 { efXVG, "-co", "scum.xvg", ffWRITE },
194 { efXVG, "-rc", "scount.xvg", ffWRITE }
196 #define NFILE asize(fnm)
198 if (!parse_common_args(&argc, argv, PCA_CAN_TIME | PCA_CAN_VIEW | PCA_BE_NICE,
199 NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, NULL, &oenv))
206 bTPS = (opt2bSet("-s", NFILE, fnm) || !opt2bSet("-n", NFILE, fnm) || bCom);
209 read_tps_conf(ftp2fn(efTPS, NFILE, fnm), title, &top, &ePBC, &xtop, NULL, box,
213 /* get index groups */
214 printf("Select a group of reference particles and a solvent group:\n");
220 get_index(&top.atoms, ftp2fn_null(efNDX, NFILE, fnm), 2, isize, index, grpname);
224 get_index(NULL, ftp2fn(efNDX, NFILE, fnm), 2, isize, index, grpname);
240 gmx_fatal(FARGS, "The number of solvent atoms (%d) is not a multiple of 3",
244 /* initialize reading trajectory: */
245 natoms = read_first_x(oenv, &status, ftp2fn(efTRX, NFILE, fnm), &t, &x, box);
249 rcut = 0.99*sqrt(max_cutoff2(guess_ePBC(box), box));
257 nbin1 = 1+(int)(2*invbw + 0.5);
258 nbin2 = 1+(int)(invbw + 0.5);
264 nrbin = 1+(int)(rcut/rbinw);
280 /* make molecules whole again */
281 gpbc = gmx_rmpbc_init(&top.idef, ePBC, natoms);
283 /* start analysis of trajectory */
288 /* make molecules whole again */
289 gmx_rmpbc(gpbc, natoms, box, x);
292 set_pbc(&pbc, ePBC, box);
296 for (p = 0; (p < nrefgrp); p++)
300 calc_com_pbc(nrefat, &top, x, &pbc, index[0], xref, bPBC);
304 copy_rvec(x[index[0][p]], xref);
307 for (m = 0; m < isize[1]; m += 3)
312 range_check(sa0, 0, natoms);
313 range_check(sa1, 0, natoms);
314 range_check(sa2, 0, natoms);
315 pbc_dx(&pbc, x[sa0], xref, dx);
322 /* Determine the normal to the plain */
323 rvec_sub(x[sa1], x[sa0], dxh1);
324 rvec_sub(x[sa2], x[sa0], dxh2);
325 rvec_inc(dxh1, dxh2);
328 inp = iprod(dx, dxh1);
329 cprod(dxh1, dxh2, outer);
331 outp = iprod(dx, outer);
335 /* Use the vector between the 2nd and 3rd atom */
336 rvec_sub(x[sa2], x[sa1], dxh2);
338 outp = iprod(dx, dxh2)/r;
341 int ii = (int)(invrbw*r);
342 range_check(ii, 0, nrbin);
344 histi2[ii] += 3*sqr(outp) - 1;
347 if ((r2 >= rmin2) && (r2 < rmax2))
349 int ii1 = (int)(invbw*(inp + 1));
350 int ii2 = (int)(invbw*fabs(outp));
352 range_check(ii1, 0, nbin1);
353 range_check(ii2, 0, nbin2);
367 while (read_next_x(oenv, status, &t, x, box));
372 gmx_rmpbc_done(gpbc);
374 /* Add the bin for the exact maximum to the previous bin */
375 hist1[nbin1-1] += hist1[nbin1];
376 hist2[nbin2-1] += hist2[nbin2];
378 nav = (real)ntot/(nrefgrp*nf);
379 normfac = invbw/ntot;
381 fprintf(stderr, "Average nr of molecules between %g and %g nm: %.1f\n",
387 fprintf(stderr, "Average cos(theta1) between %g and %g nm: %6.3f\n",
389 fprintf(stderr, "Average 3cos2(theta2)-1 between %g and %g nm: %6.3f\n",
393 sprintf(str, "Solvent orientation between %g and %g nm", rmin, rmax);
394 fp = xvgropen(opt2fn("-o", NFILE, fnm), str, "cos(\\8q\\4\\s1\\N)", "", oenv);
395 if (output_env_get_print_xvgr_codes(oenv))
397 fprintf(fp, "@ subtitle \"average shell size %.1f molecules\"\n", nav);
399 for (i = 0; i < nbin1; i++)
401 fprintf(fp, "%g %g\n", (i+0.5)*binwidth-1, 2*normfac*hist1[i]);
405 sprintf(str, "Solvent normal orientation between %g and %g nm", rmin, rmax);
406 fp = xvgropen(opt2fn("-no", NFILE, fnm), str, "cos(\\8q\\4\\s2\\N)", "", oenv);
407 if (output_env_get_print_xvgr_codes(oenv))
409 fprintf(fp, "@ subtitle \"average shell size %.1f molecules\"\n", nav);
411 for (i = 0; i < nbin2; i++)
413 fprintf(fp, "%g %g\n", (i+0.5)*binwidth, normfac*hist2[i]);
418 sprintf(str, "Solvent orientation");
419 fp = xvgropen(opt2fn("-ro", NFILE, fnm), str, "r (nm)", "", oenv);
420 if (output_env_get_print_xvgr_codes(oenv))
422 fprintf(fp, "@ subtitle \"as a function of distance\"\n");
424 xvgr_legend(fp, 2, legr, oenv);
425 for (i = 0; i < nrbin; i++)
427 fprintf(fp, "%g %g %g\n", (i+0.5)*rbinw,
428 histn[i] ? histi1[i]/histn[i] : 0,
429 histn[i] ? histi2[i]/histn[i] : 0);
433 sprintf(str, "Cumulative solvent orientation");
434 fp = xvgropen(opt2fn("-co", NFILE, fnm), str, "r (nm)", "", oenv);
435 if (output_env_get_print_xvgr_codes(oenv))
437 fprintf(fp, "@ subtitle \"as a function of distance\"\n");
439 xvgr_legend(fp, 2, legc, oenv);
440 normfac = 1.0/(nrefgrp*nf);
443 fprintf(fp, "%g %g %g\n", 0.0, c1, c2);
444 for (i = 0; i < nrbin; i++)
446 c1 += histi1[i]*normfac;
447 c2 += histi2[i]*normfac;
448 fprintf(fp, "%g %g %g\n", (i+1)*rbinw, c1, c2);
452 sprintf(str, "Solvent distribution");
453 fp = xvgropen(opt2fn("-rc", NFILE, fnm), str, "r (nm)", "molecules/nm", oenv);
454 if (output_env_get_print_xvgr_codes(oenv))
456 fprintf(fp, "@ subtitle \"as a function of distance\"\n");
458 normfac = 1.0/(rbinw*nf);
459 for (i = 0; i < nrbin; i++)
461 fprintf(fp, "%g %g\n", (i+0.5)*rbinw, histn[i]*normfac);
465 do_view(oenv, opt2fn("-o", NFILE, fnm), NULL);
466 do_view(oenv, opt2fn("-no", NFILE, fnm), NULL);
467 do_view(oenv, opt2fn("-ro", NFILE, fnm), "-nxy");
468 do_view(oenv, opt2fn("-co", NFILE, fnm), "-nxy");