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.
41 #include "gromacs/math/vec.h"
43 #include "gromacs/pbcutil/pbc.h"
44 #include "gromacs/topology/index.h"
45 #include "gromacs/fileio/tpxio.h"
46 #include "gromacs/fileio/trxio.h"
49 #include "gromacs/commandline/pargs.h"
50 #include "gromacs/fileio/xvgr.h"
51 #include "gromacs/pbcutil/rmpbc.h"
52 #include "gromacs/utility/fatalerror.h"
53 #include "gromacs/utility/smalloc.h"
55 static void calc_com_pbc(int nrefat
, t_topology
*top
, rvec x
[], t_pbc
*pbc
,
56 atom_id index
[], rvec xref
, gmx_bool bPBC
)
58 const real tol
= 1e-4;
64 /* First simple calculation */
67 for (m
= 0; (m
< nrefat
); m
++)
70 mass
= top
->atoms
.atom
[ai
].m
;
71 for (j
= 0; (j
< DIM
); j
++)
73 xref
[j
] += mass
*x
[ai
][j
];
77 svmul(1/mtot
, xref
, xref
);
78 /* Now check if any atom is more than half the box from the COM */
85 for (m
= 0; (m
< nrefat
); m
++)
88 mass
= top
->atoms
.atom
[ai
].m
/mtot
;
89 pbc_dx(pbc
, x
[ai
], xref
, dx
);
90 rvec_add(xref
, dx
, xtest
);
91 for (j
= 0; (j
< DIM
); j
++)
93 if (fabs(xtest
[j
]-x
[ai
][j
]) > tol
)
95 /* Here we have used the wrong image for contributing to the COM */
96 xref
[j
] += mass
*(xtest
[j
]-x
[ai
][j
]);
104 printf("COM: %8.3f %8.3f %8.3f iter = %d\n", xref
[XX
], xref
[YY
], xref
[ZZ
], iter
);
112 int gmx_sorient(int argc
, char *argv
[])
124 int i
, j
, p
, sa0
, sa1
, sa2
, n
, ntot
, nf
, m
, *hist1
, *hist2
, *histn
, nbin1
, nbin2
, nrbin
;
125 real
*histi1
, *histi2
, invbw
, invrbw
;
127 int *isize
, nrefgrp
, nrefat
;
130 real inp
, outp
, two_pi
, nav
, normfac
, rmin2
, rmax2
, rcut
, rcut2
, r2
, r
, mass
, mtot
;
134 rvec xref
, dx
, dxh1
, dxh2
, outer
;
135 gmx_rmpbc_t gpbc
= NULL
;
137 const char *legr
[] = {
138 "<cos(\\8q\\4\\s1\\N)>",
139 "<3cos\\S2\\N(\\8q\\4\\s2\\N)-1>"
141 const char *legc
[] = {
142 "cos(\\8q\\4\\s1\\N)",
143 "3cos\\S2\\N(\\8q\\4\\s2\\N)-1"
146 const char *desc
[] = {
147 "[THISMODULE] analyzes solvent orientation around solutes.",
148 "It calculates two angles between the vector from one or more",
149 "reference positions to the first atom of each solvent molecule:[PAR]",
150 "[GRK]theta[grk][SUB]1[sub]: the angle with the vector from the first atom of the solvent",
151 "molecule to the midpoint between atoms 2 and 3.[BR]",
152 "[GRK]theta[grk][SUB]2[sub]: the angle with the normal of the solvent plane, defined by the",
153 "same three atoms, or, when the option [TT]-v23[tt] is set, ",
154 "the angle with the vector between atoms 2 and 3.[PAR]",
155 "The reference can be a set of atoms or",
156 "the center of mass of a set of atoms. The group of solvent atoms should",
157 "consist of 3 atoms per solvent molecule.",
158 "Only solvent molecules between [TT]-rmin[tt] and [TT]-rmax[tt] are",
159 "considered for [TT]-o[tt] and [TT]-no[tt] each frame.[PAR]",
160 "[TT]-o[tt]: distribtion of [MATH][COS][GRK]theta[grk][SUB]1[sub][cos][math] for rmin<=r<=rmax.[PAR]",
161 "[TT]-no[tt]: distribution of [MATH][COS][GRK]theta[grk][SUB]2[sub][cos][math] for rmin<=r<=rmax.[PAR]",
162 "[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",
164 "[TT]-co[tt]: the sum over all solvent molecules within distance r",
165 "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]",
166 "[TT]-rc[tt]: the distribution of the solvent molecules as a function of r"
170 static gmx_bool bCom
= FALSE
, bVec23
= FALSE
, bPBC
= FALSE
;
171 static real rmin
= 0.0, rmax
= 0.5, binwidth
= 0.02, rbinw
= 0.02;
173 { "-com", FALSE
, etBOOL
, {&bCom
},
174 "Use the center of mass as the reference postion" },
175 { "-v23", FALSE
, etBOOL
, {&bVec23
},
176 "Use the vector between atoms 2 and 3" },
177 { "-rmin", FALSE
, etREAL
, {&rmin
}, "Minimum distance (nm)" },
178 { "-rmax", FALSE
, etREAL
, {&rmax
}, "Maximum distance (nm)" },
179 { "-cbin", FALSE
, etREAL
, {&binwidth
}, "Binwidth for the cosine" },
180 { "-rbin", FALSE
, etREAL
, {&rbinw
}, "Binwidth for r (nm)" },
181 { "-pbc", FALSE
, etBOOL
, {&bPBC
}, "Check PBC for the center of mass calculation. Only necessary when your reference group consists of several molecules." }
185 { efTRX
, NULL
, NULL
, ffREAD
},
186 { efTPS
, NULL
, NULL
, ffREAD
},
187 { efNDX
, NULL
, NULL
, ffOPTRD
},
188 { efXVG
, NULL
, "sori.xvg", ffWRITE
},
189 { efXVG
, "-no", "snor.xvg", ffWRITE
},
190 { efXVG
, "-ro", "sord.xvg", ffWRITE
},
191 { efXVG
, "-co", "scum.xvg", ffWRITE
},
192 { efXVG
, "-rc", "scount.xvg", ffWRITE
}
194 #define NFILE asize(fnm)
196 if (!parse_common_args(&argc
, argv
, PCA_CAN_TIME
| PCA_CAN_VIEW
| PCA_BE_NICE
,
197 NFILE
, fnm
, asize(pa
), pa
, asize(desc
), desc
, 0, NULL
, &oenv
))
204 bTPS
= (opt2bSet("-s", NFILE
, fnm
) || !opt2bSet("-n", NFILE
, fnm
) || bCom
);
207 read_tps_conf(ftp2fn(efTPS
, NFILE
, fnm
), title
, &top
, &ePBC
, &xtop
, NULL
, box
,
211 /* get index groups */
212 printf("Select a group of reference particles and a solvent group:\n");
218 get_index(&top
.atoms
, ftp2fn_null(efNDX
, NFILE
, fnm
), 2, isize
, index
, grpname
);
222 get_index(NULL
, ftp2fn(efNDX
, NFILE
, fnm
), 2, isize
, index
, grpname
);
238 gmx_fatal(FARGS
, "The number of solvent atoms (%d) is not a multiple of 3",
242 /* initialize reading trajectory: */
243 natoms
= read_first_x(oenv
, &status
, ftp2fn(efTRX
, NFILE
, fnm
), &t
, &x
, box
);
247 rcut
= 0.99*sqrt(max_cutoff2(guess_ePBC(box
), box
));
255 nbin1
= 1+(int)(2*invbw
+ 0.5);
256 nbin2
= 1+(int)(invbw
+ 0.5);
262 nrbin
= 1+(int)(rcut
/rbinw
);
278 /* make molecules whole again */
279 gpbc
= gmx_rmpbc_init(&top
.idef
, ePBC
, natoms
);
281 /* start analysis of trajectory */
286 /* make molecules whole again */
287 gmx_rmpbc(gpbc
, natoms
, box
, x
);
290 set_pbc(&pbc
, ePBC
, box
);
294 for (p
= 0; (p
< nrefgrp
); p
++)
298 calc_com_pbc(nrefat
, &top
, x
, &pbc
, index
[0], xref
, bPBC
);
302 copy_rvec(x
[index
[0][p
]], xref
);
305 for (m
= 0; m
< isize
[1]; m
+= 3)
310 range_check(sa0
, 0, natoms
);
311 range_check(sa1
, 0, natoms
);
312 range_check(sa2
, 0, natoms
);
313 pbc_dx(&pbc
, x
[sa0
], xref
, dx
);
320 /* Determine the normal to the plain */
321 rvec_sub(x
[sa1
], x
[sa0
], dxh1
);
322 rvec_sub(x
[sa2
], x
[sa0
], dxh2
);
323 rvec_inc(dxh1
, dxh2
);
326 inp
= iprod(dx
, dxh1
);
327 cprod(dxh1
, dxh2
, outer
);
329 outp
= iprod(dx
, outer
);
333 /* Use the vector between the 2nd and 3rd atom */
334 rvec_sub(x
[sa2
], x
[sa1
], dxh2
);
336 outp
= iprod(dx
, dxh2
)/r
;
339 int ii
= (int)(invrbw
*r
);
340 range_check(ii
, 0, nrbin
);
342 histi2
[ii
] += 3*sqr(outp
) - 1;
345 if ((r2
>= rmin2
) && (r2
< rmax2
))
347 int ii1
= (int)(invbw
*(inp
+ 1));
348 int ii2
= (int)(invbw
*fabs(outp
));
350 range_check(ii1
, 0, nbin1
);
351 range_check(ii2
, 0, nbin2
);
365 while (read_next_x(oenv
, status
, &t
, x
, box
));
370 gmx_rmpbc_done(gpbc
);
372 /* Add the bin for the exact maximum to the previous bin */
373 hist1
[nbin1
-1] += hist1
[nbin1
];
374 hist2
[nbin2
-1] += hist2
[nbin2
];
376 nav
= (real
)ntot
/(nrefgrp
*nf
);
377 normfac
= invbw
/ntot
;
379 fprintf(stderr
, "Average nr of molecules between %g and %g nm: %.1f\n",
385 fprintf(stderr
, "Average cos(theta1) between %g and %g nm: %6.3f\n",
387 fprintf(stderr
, "Average 3cos2(theta2)-1 between %g and %g nm: %6.3f\n",
391 sprintf(str
, "Solvent orientation between %g and %g nm", rmin
, rmax
);
392 fp
= xvgropen(opt2fn("-o", NFILE
, fnm
), str
, "cos(\\8q\\4\\s1\\N)", "", oenv
);
393 if (output_env_get_print_xvgr_codes(oenv
))
395 fprintf(fp
, "@ subtitle \"average shell size %.1f molecules\"\n", nav
);
397 for (i
= 0; i
< nbin1
; i
++)
399 fprintf(fp
, "%g %g\n", (i
+0.5)*binwidth
-1, 2*normfac
*hist1
[i
]);
403 sprintf(str
, "Solvent normal orientation between %g and %g nm", rmin
, rmax
);
404 fp
= xvgropen(opt2fn("-no", NFILE
, fnm
), str
, "cos(\\8q\\4\\s2\\N)", "", oenv
);
405 if (output_env_get_print_xvgr_codes(oenv
))
407 fprintf(fp
, "@ subtitle \"average shell size %.1f molecules\"\n", nav
);
409 for (i
= 0; i
< nbin2
; i
++)
411 fprintf(fp
, "%g %g\n", (i
+0.5)*binwidth
, normfac
*hist2
[i
]);
416 sprintf(str
, "Solvent orientation");
417 fp
= xvgropen(opt2fn("-ro", NFILE
, fnm
), str
, "r (nm)", "", oenv
);
418 if (output_env_get_print_xvgr_codes(oenv
))
420 fprintf(fp
, "@ subtitle \"as a function of distance\"\n");
422 xvgr_legend(fp
, 2, legr
, oenv
);
423 for (i
= 0; i
< nrbin
; i
++)
425 fprintf(fp
, "%g %g %g\n", (i
+0.5)*rbinw
,
426 histn
[i
] ? histi1
[i
]/histn
[i
] : 0,
427 histn
[i
] ? histi2
[i
]/histn
[i
] : 0);
431 sprintf(str
, "Cumulative solvent orientation");
432 fp
= xvgropen(opt2fn("-co", NFILE
, fnm
), str
, "r (nm)", "", oenv
);
433 if (output_env_get_print_xvgr_codes(oenv
))
435 fprintf(fp
, "@ subtitle \"as a function of distance\"\n");
437 xvgr_legend(fp
, 2, legc
, oenv
);
438 normfac
= 1.0/(nrefgrp
*nf
);
441 fprintf(fp
, "%g %g %g\n", 0.0, c1
, c2
);
442 for (i
= 0; i
< nrbin
; i
++)
444 c1
+= histi1
[i
]*normfac
;
445 c2
+= histi2
[i
]*normfac
;
446 fprintf(fp
, "%g %g %g\n", (i
+1)*rbinw
, c1
, c2
);
450 sprintf(str
, "Solvent distribution");
451 fp
= xvgropen(opt2fn("-rc", NFILE
, fnm
), str
, "r (nm)", "molecules/nm", oenv
);
452 if (output_env_get_print_xvgr_codes(oenv
))
454 fprintf(fp
, "@ subtitle \"as a function of distance\"\n");
456 normfac
= 1.0/(rbinw
*nf
);
457 for (i
= 0; i
< nrbin
; i
++)
459 fprintf(fp
, "%g %g\n", (i
+0.5)*rbinw
, histn
[i
]*normfac
);
463 do_view(oenv
, opt2fn("-o", NFILE
, fnm
), NULL
);
464 do_view(oenv
, opt2fn("-no", NFILE
, fnm
), NULL
);
465 do_view(oenv
, opt2fn("-ro", NFILE
, fnm
), "-nxy");
466 do_view(oenv
, opt2fn("-co", NFILE
, fnm
), "-nxy");