3 * This source code is part of
7 * GROningen MAchine for Chemical Simulations
10 * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
11 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
12 * Copyright (c) 2001-2004, The GROMACS development team,
13 * check out http://www.gromacs.org for more information.
15 * This program is free software; you can redistribute it and/or
16 * modify it under the terms of the GNU General Public License
17 * as published by the Free Software Foundation; either version 2
18 * of the License, or (at your option) any later version.
20 * If you want to redistribute modifications, please consider that
21 * scientific software is very special. Version control is crucial -
22 * bugs must be traceable. We will be happy to consider code for
23 * inclusion in the official distribution, but derived work must not
24 * be called official GROMACS. Details are found in the README & COPYING
25 * files - if they are missing, get the official version at www.gromacs.org.
27 * To help us fund GROMACS development, we humbly ask that you cite
28 * the papers on the package - you can find them in the top README file.
30 * For more info, check our website at http://www.gromacs.org
33 * Green Red Orange Magenta Azure Cyan Skyblue
52 static void calc_com_pbc(int nrefat
, t_topology
*top
, rvec x
[], t_pbc
*pbc
,
53 atom_id index
[], rvec xref
, gmx_bool bPBC
, matrix box
)
55 const real tol
= 1e-4;
61 /* First simple calculation */
64 for (m
= 0; (m
< nrefat
); m
++)
67 mass
= top
->atoms
.atom
[ai
].m
;
68 for (j
= 0; (j
< DIM
); j
++)
70 xref
[j
] += mass
*x
[ai
][j
];
74 svmul(1/mtot
, xref
, xref
);
75 /* Now check if any atom is more than half the box from the COM */
82 for (m
= 0; (m
< nrefat
); m
++)
85 mass
= top
->atoms
.atom
[ai
].m
/mtot
;
86 pbc_dx(pbc
, x
[ai
], xref
, dx
);
87 rvec_add(xref
, dx
, xtest
);
88 for (j
= 0; (j
< DIM
); j
++)
90 if (fabs(xtest
[j
]-x
[ai
][j
]) > tol
)
92 /* Here we have used the wrong image for contributing to the COM */
93 xref
[j
] += mass
*(xtest
[j
]-x
[ai
][j
]);
101 printf("COM: %8.3f %8.3f %8.3f iter = %d\n", xref
[XX
], xref
[YY
], xref
[ZZ
], iter
);
109 int gmx_sorient(int argc
, char *argv
[])
121 int i
, j
, p
, sa0
, sa1
, sa2
, n
, ntot
, nf
, m
, *hist1
, *hist2
, *histn
, nbin1
, nbin2
, nrbin
;
122 real
*histi1
, *histi2
, invbw
, invrbw
;
124 int *isize
, nrefgrp
, nrefat
;
127 real inp
, outp
, two_pi
, nav
, normfac
, rmin2
, rmax2
, rcut
, rcut2
, r2
, r
, mass
, mtot
;
131 rvec xref
, dx
, dxh1
, dxh2
, outer
;
132 gmx_rmpbc_t gpbc
= NULL
;
134 const char *legr
[] = {
135 "<cos(\\8q\\4\\s1\\N)>",
136 "<3cos\\S2\\N(\\8q\\4\\s2\\N)-1>"
138 const char *legc
[] = {
139 "cos(\\8q\\4\\s1\\N)",
140 "3cos\\S2\\N(\\8q\\4\\s2\\N)-1"
143 const char *desc
[] = {
144 "[TT]g_sorient[tt] analyzes solvent orientation around solutes.",
145 "It calculates two angles between the vector from one or more",
146 "reference positions to the first atom of each solvent molecule:[PAR]",
147 "[GRK]theta[grk][SUB]1[sub]: the angle with the vector from the first atom of the solvent",
148 "molecule to the midpoint between atoms 2 and 3.[BR]",
149 "[GRK]theta[grk][SUB]2[sub]: the angle with the normal of the solvent plane, defined by the",
150 "same three atoms, or, when the option [TT]-v23[tt] is set, ",
151 "the angle with the vector between atoms 2 and 3.[PAR]",
152 "The reference can be a set of atoms or",
153 "the center of mass of a set of atoms. The group of solvent atoms should",
154 "consist of 3 atoms per solvent molecule.",
155 "Only solvent molecules between [TT]-rmin[tt] and [TT]-rmax[tt] are",
156 "considered for [TT]-o[tt] and [TT]-no[tt] each frame.[PAR]",
157 "[TT]-o[tt]: distribtion of [MATH][COS][GRK]theta[grk][SUB]1[sub][cos][math] for rmin<=r<=rmax.[PAR]",
158 "[TT]-no[tt]: distribution of [MATH][COS][GRK]theta[grk][SUB]2[sub][cos][math] for rmin<=r<=rmax.[PAR]",
159 "[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",
161 "[TT]-co[tt]: the sum over all solvent molecules within distance r",
162 "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]",
163 "[TT]-rc[tt]: the distribution of the solvent molecules as a function of r"
167 static gmx_bool bCom
= FALSE
, bVec23
= FALSE
, bPBC
= FALSE
;
168 static real rmin
= 0.0, rmax
= 0.5, binwidth
= 0.02, rbinw
= 0.02;
170 { "-com", FALSE
, etBOOL
, {&bCom
},
171 "Use the center of mass as the reference postion" },
172 { "-v23", FALSE
, etBOOL
, {&bVec23
},
173 "Use the vector between atoms 2 and 3" },
174 { "-rmin", FALSE
, etREAL
, {&rmin
}, "Minimum distance (nm)" },
175 { "-rmax", FALSE
, etREAL
, {&rmax
}, "Maximum distance (nm)" },
176 { "-cbin", FALSE
, etREAL
, {&binwidth
}, "Binwidth for the cosine" },
177 { "-rbin", FALSE
, etREAL
, {&rbinw
}, "Binwidth for r (nm)" },
178 { "-pbc", FALSE
, etBOOL
, {&bPBC
}, "Check PBC for the center of mass calculation. Only necessary when your reference group consists of several molecules." }
182 { efTRX
, NULL
, NULL
, ffREAD
},
183 { efTPS
, NULL
, NULL
, ffREAD
},
184 { efNDX
, NULL
, NULL
, ffOPTRD
},
185 { efXVG
, NULL
, "sori.xvg", ffWRITE
},
186 { efXVG
, "-no", "snor.xvg", ffWRITE
},
187 { efXVG
, "-ro", "sord.xvg", ffWRITE
},
188 { efXVG
, "-co", "scum.xvg", ffWRITE
},
189 { efXVG
, "-rc", "scount.xvg", ffWRITE
}
191 #define NFILE asize(fnm)
193 parse_common_args(&argc
, argv
, PCA_CAN_TIME
| PCA_CAN_VIEW
| PCA_BE_NICE
,
194 NFILE
, fnm
, asize(pa
), pa
, asize(desc
), desc
, 0, NULL
, &oenv
);
198 bTPS
= (opt2bSet("-s", NFILE
, fnm
) || !opt2bSet("-n", NFILE
, fnm
) || bCom
);
201 read_tps_conf(ftp2fn(efTPS
, NFILE
, fnm
), title
, &top
, &ePBC
, &xtop
, NULL
, box
,
205 /* get index groups */
206 printf("Select a group of reference particles and a solvent group:\n");
212 get_index(&top
.atoms
, ftp2fn_null(efNDX
, NFILE
, fnm
), 2, isize
, index
, grpname
);
216 get_index(NULL
, ftp2fn(efNDX
, NFILE
, fnm
), 2, isize
, index
, grpname
);
232 gmx_fatal(FARGS
, "The number of solvent atoms (%d) is not a multiple of 3",
236 /* initialize reading trajectory: */
237 natoms
= read_first_x(oenv
, &status
, ftp2fn(efTRX
, NFILE
, fnm
), &t
, &x
, box
);
241 rcut
= 0.99*sqrt(max_cutoff2(guess_ePBC(box
), box
));
249 nbin1
= 1+(int)(2*invbw
+ 0.5);
250 nbin2
= 1+(int)(invbw
+ 0.5);
256 nrbin
= 1+(int)(rcut
/rbinw
);
272 /* make molecules whole again */
273 gpbc
= gmx_rmpbc_init(&top
.idef
, ePBC
, natoms
, box
);
275 /* start analysis of trajectory */
280 /* make molecules whole again */
281 gmx_rmpbc(gpbc
, natoms
, box
, x
);
284 set_pbc(&pbc
, ePBC
, box
);
288 for (p
= 0; (p
< nrefgrp
); p
++)
292 calc_com_pbc(nrefat
, &top
, x
, &pbc
, index
[0], xref
, bPBC
, box
);
296 copy_rvec(x
[index
[0][p
]], xref
);
299 for (m
= 0; m
< isize
[1]; m
+= 3)
304 range_check(sa0
, 0, natoms
);
305 range_check(sa1
, 0, natoms
);
306 range_check(sa2
, 0, natoms
);
307 pbc_dx(&pbc
, x
[sa0
], xref
, dx
);
314 /* Determine the normal to the plain */
315 rvec_sub(x
[sa1
], x
[sa0
], dxh1
);
316 rvec_sub(x
[sa2
], x
[sa0
], dxh2
);
317 rvec_inc(dxh1
, dxh2
);
320 inp
= iprod(dx
, dxh1
);
321 cprod(dxh1
, dxh2
, outer
);
323 outp
= iprod(dx
, outer
);
327 /* Use the vector between the 2nd and 3rd atom */
328 rvec_sub(x
[sa2
], x
[sa1
], dxh2
);
330 outp
= iprod(dx
, dxh2
)/r
;
333 int ii
= (int)(invrbw
*r
);
334 range_check(ii
, 0, nrbin
);
336 histi2
[ii
] += 3*sqr(outp
) - 1;
339 if ((r2
>= rmin2
) && (r2
< rmax2
))
341 int ii1
= (int)(invbw
*(inp
+ 1));
342 int ii2
= (int)(invbw
*fabs(outp
));
344 range_check(ii1
, 0, nbin1
);
345 range_check(ii2
, 0, nbin2
);
359 while (read_next_x(oenv
, status
, &t
, natoms
, x
, box
));
364 gmx_rmpbc_done(gpbc
);
366 /* Add the bin for the exact maximum to the previous bin */
367 hist1
[nbin1
-1] += hist1
[nbin1
];
368 hist2
[nbin2
-1] += hist2
[nbin2
];
370 nav
= (real
)ntot
/(nrefgrp
*nf
);
371 normfac
= invbw
/ntot
;
373 fprintf(stderr
, "Average nr of molecules between %g and %g nm: %.1f\n",
379 fprintf(stderr
, "Average cos(theta1) between %g and %g nm: %6.3f\n",
381 fprintf(stderr
, "Average 3cos2(theta2)-1 between %g and %g nm: %6.3f\n",
385 sprintf(str
, "Solvent orientation between %g and %g nm", rmin
, rmax
);
386 fp
= xvgropen(opt2fn("-o", NFILE
, fnm
), str
, "cos(\\8q\\4\\s1\\N)", "", oenv
);
387 if (output_env_get_print_xvgr_codes(oenv
))
389 fprintf(fp
, "@ subtitle \"average shell size %.1f molecules\"\n", nav
);
391 for (i
= 0; i
< nbin1
; i
++)
393 fprintf(fp
, "%g %g\n", (i
+0.5)*binwidth
-1, 2*normfac
*hist1
[i
]);
397 sprintf(str
, "Solvent normal orientation between %g and %g nm", rmin
, rmax
);
398 fp
= xvgropen(opt2fn("-no", NFILE
, fnm
), str
, "cos(\\8q\\4\\s2\\N)", "", oenv
);
399 if (output_env_get_print_xvgr_codes(oenv
))
401 fprintf(fp
, "@ subtitle \"average shell size %.1f molecules\"\n", nav
);
403 for (i
= 0; i
< nbin2
; i
++)
405 fprintf(fp
, "%g %g\n", (i
+0.5)*binwidth
, normfac
*hist2
[i
]);
410 sprintf(str
, "Solvent orientation");
411 fp
= xvgropen(opt2fn("-ro", NFILE
, fnm
), str
, "r (nm)", "", oenv
);
412 if (output_env_get_print_xvgr_codes(oenv
))
414 fprintf(fp
, "@ subtitle \"as a function of distance\"\n");
416 xvgr_legend(fp
, 2, legr
, oenv
);
417 for (i
= 0; i
< nrbin
; i
++)
419 fprintf(fp
, "%g %g %g\n", (i
+0.5)*rbinw
,
420 histn
[i
] ? histi1
[i
]/histn
[i
] : 0,
421 histn
[i
] ? histi2
[i
]/histn
[i
] : 0);
425 sprintf(str
, "Cumulative solvent orientation");
426 fp
= xvgropen(opt2fn("-co", NFILE
, fnm
), str
, "r (nm)", "", oenv
);
427 if (output_env_get_print_xvgr_codes(oenv
))
429 fprintf(fp
, "@ subtitle \"as a function of distance\"\n");
431 xvgr_legend(fp
, 2, legc
, oenv
);
432 normfac
= 1.0/(nrefgrp
*nf
);
435 fprintf(fp
, "%g %g %g\n", 0.0, c1
, c2
);
436 for (i
= 0; i
< nrbin
; i
++)
438 c1
+= histi1
[i
]*normfac
;
439 c2
+= histi2
[i
]*normfac
;
440 fprintf(fp
, "%g %g %g\n", (i
+1)*rbinw
, c1
, c2
);
444 sprintf(str
, "Solvent distribution");
445 fp
= xvgropen(opt2fn("-rc", NFILE
, fnm
), str
, "r (nm)", "molecules/nm", oenv
);
446 if (output_env_get_print_xvgr_codes(oenv
))
448 fprintf(fp
, "@ subtitle \"as a function of distance\"\n");
450 normfac
= 1.0/(rbinw
*nf
);
451 for (i
= 0; i
< nrbin
; i
++)
453 fprintf(fp
, "%g %g\n", (i
+0.5)*rbinw
, histn
[i
]*normfac
);
457 do_view(oenv
, opt2fn("-o", NFILE
, fnm
), NULL
);
458 do_view(oenv
, opt2fn("-no", NFILE
, fnm
), NULL
);
459 do_view(oenv
, opt2fn("-ro", NFILE
, fnm
), "-nxy");
460 do_view(oenv
, opt2fn("-co", NFILE
, fnm
), "-nxy");