File: | gromacs/gmxana/gmx_sorient.c |
Location: | line 203, column 5 |
Description: | Value stored to 'two_pi' is never read |
1 | /* |
2 | * This file is part of the GROMACS molecular simulation package. |
3 | * |
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. |
10 | * |
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. |
15 | * |
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. |
20 | * |
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. |
25 | * |
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. |
33 | * |
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. |
36 | */ |
37 | #ifdef HAVE_CONFIG_H1 |
38 | #include <config.h> |
39 | #endif |
40 | |
41 | #include "macros.h" |
42 | #include "gromacs/commandline/pargs.h" |
43 | #include "gromacs/utility/smalloc.h" |
44 | #include "gstat.h" |
45 | #include "gromacs/math/vec.h" |
46 | #include "gromacs/fileio/xvgr.h" |
47 | #include "viewit.h" |
48 | #include "pbc.h" |
49 | #include "index.h" |
50 | #include "gromacs/fileio/tpxio.h" |
51 | #include "gromacs/fileio/trxio.h" |
52 | #include "gmx_ana.h" |
53 | |
54 | #include "gromacs/utility/fatalerror.h" |
55 | |
56 | static void calc_com_pbc(int nrefat, t_topology *top, rvec x[], t_pbc *pbc, |
57 | atom_id index[], rvec xref, gmx_bool bPBC) |
58 | { |
59 | const real tol = 1e-4; |
60 | gmx_bool bChanged; |
61 | int m, j, ai, iter; |
62 | real mass, mtot; |
63 | rvec dx, xtest; |
64 | |
65 | /* First simple calculation */ |
66 | clear_rvec(xref); |
67 | mtot = 0; |
68 | for (m = 0; (m < nrefat); m++) |
69 | { |
70 | ai = index[m]; |
71 | mass = top->atoms.atom[ai].m; |
72 | for (j = 0; (j < DIM3); j++) |
73 | { |
74 | xref[j] += mass*x[ai][j]; |
75 | } |
76 | mtot += mass; |
77 | } |
78 | svmul(1/mtot, xref, xref); |
79 | /* Now check if any atom is more than half the box from the COM */ |
80 | if (bPBC) |
81 | { |
82 | iter = 0; |
83 | do |
84 | { |
85 | bChanged = FALSE0; |
86 | for (m = 0; (m < nrefat); m++) |
87 | { |
88 | ai = index[m]; |
89 | mass = top->atoms.atom[ai].m/mtot; |
90 | pbc_dx(pbc, x[ai], xref, dx); |
91 | rvec_add(xref, dx, xtest); |
92 | for (j = 0; (j < DIM3); j++) |
93 | { |
94 | if (fabs(xtest[j]-x[ai][j]) > tol) |
95 | { |
96 | /* Here we have used the wrong image for contributing to the COM */ |
97 | xref[j] += mass*(xtest[j]-x[ai][j]); |
98 | x[ai][j] = xtest[j]; |
99 | bChanged = TRUE1; |
100 | } |
101 | } |
102 | } |
103 | if (bChanged) |
104 | { |
105 | printf("COM: %8.3f %8.3f %8.3f iter = %d\n", xref[XX0], xref[YY1], xref[ZZ2], iter); |
106 | } |
107 | iter++; |
108 | } |
109 | while (bChanged); |
110 | } |
111 | } |
112 | |
113 | int gmx_sorient(int argc, char *argv[]) |
114 | { |
115 | t_topology top; |
116 | int ePBC = -1; |
117 | char title[STRLEN4096]; |
118 | t_trxstatus *status; |
119 | int natoms; |
120 | real t; |
121 | rvec *xtop, *x; |
122 | matrix box; |
123 | |
124 | FILE *fp; |
125 | int i, j, p, sa0, sa1, sa2, n, ntot, nf, m, *hist1, *hist2, *histn, nbin1, nbin2, nrbin; |
126 | real *histi1, *histi2, invbw, invrbw; |
127 | double sum1, sum2; |
128 | int *isize, nrefgrp, nrefat; |
129 | atom_id **index; |
130 | char **grpname; |
131 | real inp, outp, two_pi, nav, normfac, rmin2, rmax2, rcut, rcut2, r2, r, mass, mtot; |
132 | real c1, c2; |
133 | char str[STRLEN4096]; |
134 | gmx_bool bTPS; |
135 | rvec xref, dx, dxh1, dxh2, outer; |
136 | gmx_rmpbc_t gpbc = NULL((void*)0); |
137 | t_pbc pbc; |
138 | const char *legr[] = { |
139 | "<cos(\\8q\\4\\s1\\N)>", |
140 | "<3cos\\S2\\N(\\8q\\4\\s2\\N)-1>" |
141 | }; |
142 | const char *legc[] = { |
143 | "cos(\\8q\\4\\s1\\N)", |
144 | "3cos\\S2\\N(\\8q\\4\\s2\\N)-1" |
145 | }; |
146 | |
147 | const char *desc[] = { |
148 | "[THISMODULE] analyzes solvent orientation around solutes.", |
149 | "It calculates two angles between the vector from one or more", |
150 | "reference positions to the first atom of each solvent molecule:[PAR]", |
151 | "[GRK]theta[grk][SUB]1[sub]: the angle with the vector from the first atom of the solvent", |
152 | "molecule to the midpoint between atoms 2 and 3.[BR]", |
153 | "[GRK]theta[grk][SUB]2[sub]: the angle with the normal of the solvent plane, defined by the", |
154 | "same three atoms, or, when the option [TT]-v23[tt] is set, ", |
155 | "the angle with the vector between atoms 2 and 3.[PAR]", |
156 | "The reference can be a set of atoms or", |
157 | "the center of mass of a set of atoms. The group of solvent atoms should", |
158 | "consist of 3 atoms per solvent molecule.", |
159 | "Only solvent molecules between [TT]-rmin[tt] and [TT]-rmax[tt] are", |
160 | "considered for [TT]-o[tt] and [TT]-no[tt] each frame.[PAR]", |
161 | "[TT]-o[tt]: distribtion of [MATH][COS][GRK]theta[grk][SUB]1[sub][cos][math] for rmin<=r<=rmax.[PAR]", |
162 | "[TT]-no[tt]: distribution of [MATH][COS][GRK]theta[grk][SUB]2[sub][cos][math] for rmin<=r<=rmax.[PAR]", |
163 | "[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 | "distance.[PAR]", |
165 | "[TT]-co[tt]: the sum over all solvent molecules within distance r", |
166 | "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]", |
167 | "[TT]-rc[tt]: the distribution of the solvent molecules as a function of r" |
168 | }; |
169 | |
170 | output_env_t oenv; |
171 | static gmx_bool bCom = FALSE0, bVec23 = FALSE0, bPBC = FALSE0; |
172 | static real rmin = 0.0, rmax = 0.5, binwidth = 0.02, rbinw = 0.02; |
173 | t_pargs pa[] = { |
174 | { "-com", FALSE0, etBOOL, {&bCom}, |
175 | "Use the center of mass as the reference postion" }, |
176 | { "-v23", FALSE0, etBOOL, {&bVec23}, |
177 | "Use the vector between atoms 2 and 3" }, |
178 | { "-rmin", FALSE0, etREAL, {&rmin}, "Minimum distance (nm)" }, |
179 | { "-rmax", FALSE0, etREAL, {&rmax}, "Maximum distance (nm)" }, |
180 | { "-cbin", FALSE0, etREAL, {&binwidth}, "Binwidth for the cosine" }, |
181 | { "-rbin", FALSE0, etREAL, {&rbinw}, "Binwidth for r (nm)" }, |
182 | { "-pbc", FALSE0, etBOOL, {&bPBC}, "Check PBC for the center of mass calculation. Only necessary when your reference group consists of several molecules." } |
183 | }; |
184 | |
185 | t_filenm fnm[] = { |
186 | { efTRX, NULL((void*)0), NULL((void*)0), ffREAD1<<1 }, |
187 | { efTPS, NULL((void*)0), NULL((void*)0), ffREAD1<<1 }, |
188 | { efNDX, NULL((void*)0), NULL((void*)0), ffOPTRD(1<<1 | 1<<3) }, |
189 | { efXVG, NULL((void*)0), "sori.xvg", ffWRITE1<<2 }, |
190 | { efXVG, "-no", "snor.xvg", ffWRITE1<<2 }, |
191 | { efXVG, "-ro", "sord.xvg", ffWRITE1<<2 }, |
192 | { efXVG, "-co", "scum.xvg", ffWRITE1<<2 }, |
193 | { efXVG, "-rc", "scount.xvg", ffWRITE1<<2 } |
194 | }; |
195 | #define NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))) asize(fnm)((int)(sizeof(fnm)/sizeof((fnm)[0]))) |
196 | |
197 | if (!parse_common_args(&argc, argv, PCA_CAN_TIME((1<<6) | (1<<7) | (1<<14)) | PCA_CAN_VIEW(1<<5) | PCA_BE_NICE(1<<13), |
198 | NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm, asize(pa)((int)(sizeof(pa)/sizeof((pa)[0]))), pa, asize(desc)((int)(sizeof(desc)/sizeof((desc)[0]))), desc, 0, NULL((void*)0), &oenv)) |
199 | { |
200 | return 0; |
201 | } |
202 | |
203 | two_pi = 2/M_PI3.14159265358979323846; |
Value stored to 'two_pi' is never read | |
204 | |
205 | bTPS = (opt2bSet("-s", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm) || !opt2bSet("-n", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm) || bCom); |
206 | if (bTPS) |
207 | { |
208 | read_tps_conf(ftp2fn(efTPS, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), title, &top, &ePBC, &xtop, NULL((void*)0), box, |
209 | bCom); |
210 | } |
211 | |
212 | /* get index groups */ |
213 | printf("Select a group of reference particles and a solvent group:\n"); |
214 | snew(grpname, 2)(grpname) = save_calloc("grpname", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_sorient.c" , 214, (2), sizeof(*(grpname))); |
215 | snew(index, 2)(index) = save_calloc("index", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_sorient.c" , 215, (2), sizeof(*(index))); |
216 | snew(isize, 2)(isize) = save_calloc("isize", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_sorient.c" , 216, (2), sizeof(*(isize))); |
217 | if (bTPS) |
218 | { |
219 | get_index(&top.atoms, ftp2fn_null(efNDX, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), 2, isize, index, grpname); |
220 | } |
221 | else |
222 | { |
223 | get_index(NULL((void*)0), ftp2fn(efNDX, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), 2, isize, index, grpname); |
224 | } |
225 | |
226 | if (bCom) |
227 | { |
228 | nrefgrp = 1; |
229 | nrefat = isize[0]; |
230 | } |
231 | else |
232 | { |
233 | nrefgrp = isize[0]; |
234 | nrefat = 1; |
235 | } |
236 | |
237 | if (isize[1] % 3) |
238 | { |
239 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_sorient.c" , 239, "The number of solvent atoms (%d) is not a multiple of 3", |
240 | isize[1]); |
241 | } |
242 | |
243 | /* initialize reading trajectory: */ |
244 | natoms = read_first_x(oenv, &status, ftp2fn(efTRX, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), &t, &x, box); |
245 | |
246 | rmin2 = sqr(rmin); |
247 | rmax2 = sqr(rmax); |
248 | rcut = 0.99*sqrt(max_cutoff2(guess_ePBC(box), box)); |
249 | if (rcut == 0) |
250 | { |
251 | rcut = 10*rmax; |
252 | } |
253 | rcut2 = sqr(rcut); |
254 | |
255 | invbw = 1/binwidth; |
256 | nbin1 = 1+(int)(2*invbw + 0.5); |
257 | nbin2 = 1+(int)(invbw + 0.5); |
258 | |
259 | invrbw = 1/rbinw; |
260 | |
261 | snew(hist1, nbin1)(hist1) = save_calloc("hist1", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_sorient.c" , 261, (nbin1), sizeof(*(hist1))); |
262 | snew(hist2, nbin2)(hist2) = save_calloc("hist2", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_sorient.c" , 262, (nbin2), sizeof(*(hist2))); |
263 | nrbin = 1+(int)(rcut/rbinw); |
264 | if (nrbin == 0) |
265 | { |
266 | nrbin = 1; |
267 | } |
268 | snew(histi1, nrbin)(histi1) = save_calloc("histi1", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_sorient.c" , 268, (nrbin), sizeof(*(histi1))); |
269 | snew(histi2, nrbin)(histi2) = save_calloc("histi2", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_sorient.c" , 269, (nrbin), sizeof(*(histi2))); |
270 | snew(histn, nrbin)(histn) = save_calloc("histn", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_sorient.c" , 270, (nrbin), sizeof(*(histn))); |
271 | |
272 | ntot = 0; |
273 | nf = 0; |
274 | sum1 = 0; |
275 | sum2 = 0; |
276 | |
277 | if (bTPS) |
278 | { |
279 | /* make molecules whole again */ |
280 | gpbc = gmx_rmpbc_init(&top.idef, ePBC, natoms); |
281 | } |
282 | /* start analysis of trajectory */ |
283 | do |
284 | { |
285 | if (bTPS) |
286 | { |
287 | /* make molecules whole again */ |
288 | gmx_rmpbc(gpbc, natoms, box, x); |
289 | } |
290 | |
291 | set_pbc(&pbc, ePBC, box); |
292 | n = 0; |
293 | inp = 0; |
294 | outp = 0; |
295 | for (p = 0; (p < nrefgrp); p++) |
296 | { |
297 | if (bCom) |
298 | { |
299 | calc_com_pbc(nrefat, &top, x, &pbc, index[0], xref, bPBC); |
300 | } |
301 | else |
302 | { |
303 | copy_rvec(x[index[0][p]], xref); |
304 | } |
305 | |
306 | for (m = 0; m < isize[1]; m += 3) |
307 | { |
308 | sa0 = index[1][m]; |
309 | sa1 = index[1][m+1]; |
310 | sa2 = index[1][m+2]; |
311 | range_check(sa0, 0, natoms)_range_check(sa0, 0, natoms, ((void*)0),"sa0", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_sorient.c" , 311); |
312 | range_check(sa1, 0, natoms)_range_check(sa1, 0, natoms, ((void*)0),"sa1", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_sorient.c" , 312); |
313 | range_check(sa2, 0, natoms)_range_check(sa2, 0, natoms, ((void*)0),"sa2", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_sorient.c" , 313); |
314 | pbc_dx(&pbc, x[sa0], xref, dx); |
315 | r2 = norm2(dx); |
316 | if (r2 < rcut2) |
317 | { |
318 | r = sqrt(r2); |
319 | if (!bVec23) |
320 | { |
321 | /* Determine the normal to the plain */ |
322 | rvec_sub(x[sa1], x[sa0], dxh1); |
323 | rvec_sub(x[sa2], x[sa0], dxh2); |
324 | rvec_inc(dxh1, dxh2); |
325 | svmul(1/r, dx, dx); |
326 | unitv(dxh1, dxh1); |
327 | inp = iprod(dx, dxh1); |
328 | cprod(dxh1, dxh2, outer); |
329 | unitv(outer, outer); |
330 | outp = iprod(dx, outer); |
331 | } |
332 | else |
333 | { |
334 | /* Use the vector between the 2nd and 3rd atom */ |
335 | rvec_sub(x[sa2], x[sa1], dxh2); |
336 | unitv(dxh2, dxh2); |
337 | outp = iprod(dx, dxh2)/r; |
338 | } |
339 | { |
340 | int ii = (int)(invrbw*r); |
341 | range_check(ii, 0, nrbin)_range_check(ii, 0, nrbin, ((void*)0),"ii", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_sorient.c" , 341); |
342 | histi1[ii] += inp; |
343 | histi2[ii] += 3*sqr(outp) - 1; |
344 | histn[ii]++; |
345 | } |
346 | if ((r2 >= rmin2) && (r2 < rmax2)) |
347 | { |
348 | int ii1 = (int)(invbw*(inp + 1)); |
349 | int ii2 = (int)(invbw*fabs(outp)); |
350 | |
351 | range_check(ii1, 0, nbin1)_range_check(ii1, 0, nbin1, ((void*)0),"ii1", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_sorient.c" , 351); |
352 | range_check(ii2, 0, nbin2)_range_check(ii2, 0, nbin2, ((void*)0),"ii2", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_sorient.c" , 352); |
353 | hist1[ii1]++; |
354 | hist2[ii2]++; |
355 | sum1 += inp; |
356 | sum2 += outp; |
357 | n++; |
358 | } |
359 | } |
360 | } |
361 | } |
362 | ntot += n; |
363 | nf++; |
364 | |
365 | } |
366 | while (read_next_x(oenv, status, &t, x, box)); |
367 | |
368 | /* clean up */ |
369 | sfree(x)save_free("x", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_sorient.c" , 369, (x)); |
370 | close_trj(status); |
371 | gmx_rmpbc_done(gpbc); |
372 | |
373 | /* Add the bin for the exact maximum to the previous bin */ |
374 | hist1[nbin1-1] += hist1[nbin1]; |
375 | hist2[nbin2-1] += hist2[nbin2]; |
376 | |
377 | nav = (real)ntot/(nrefgrp*nf); |
378 | normfac = invbw/ntot; |
379 | |
380 | fprintf(stderrstderr, "Average nr of molecules between %g and %g nm: %.1f\n", |
381 | rmin, rmax, nav); |
382 | if (ntot > 0) |
383 | { |
384 | sum1 /= ntot; |
385 | sum2 /= ntot; |
386 | fprintf(stderrstderr, "Average cos(theta1) between %g and %g nm: %6.3f\n", |
387 | rmin, rmax, sum1); |
388 | fprintf(stderrstderr, "Average 3cos2(theta2)-1 between %g and %g nm: %6.3f\n", |
389 | rmin, rmax, sum2); |
390 | } |
391 | |
392 | sprintf(str, "Solvent orientation between %g and %g nm", rmin, rmax); |
393 | fp = xvgropen(opt2fn("-o", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), str, "cos(\\8q\\4\\s1\\N)", "", oenv); |
394 | if (output_env_get_print_xvgr_codes(oenv)) |
395 | { |
396 | fprintf(fp, "@ subtitle \"average shell size %.1f molecules\"\n", nav); |
397 | } |
398 | for (i = 0; i < nbin1; i++) |
399 | { |
400 | fprintf(fp, "%g %g\n", (i+0.5)*binwidth-1, 2*normfac*hist1[i]); |
401 | } |
402 | gmx_ffclose(fp); |
403 | |
404 | sprintf(str, "Solvent normal orientation between %g and %g nm", rmin, rmax); |
405 | fp = xvgropen(opt2fn("-no", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), str, "cos(\\8q\\4\\s2\\N)", "", oenv); |
406 | if (output_env_get_print_xvgr_codes(oenv)) |
407 | { |
408 | fprintf(fp, "@ subtitle \"average shell size %.1f molecules\"\n", nav); |
409 | } |
410 | for (i = 0; i < nbin2; i++) |
411 | { |
412 | fprintf(fp, "%g %g\n", (i+0.5)*binwidth, normfac*hist2[i]); |
413 | } |
414 | gmx_ffclose(fp); |
415 | |
416 | |
417 | sprintf(str, "Solvent orientation"); |
418 | fp = xvgropen(opt2fn("-ro", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), str, "r (nm)", "", oenv); |
419 | if (output_env_get_print_xvgr_codes(oenv)) |
420 | { |
421 | fprintf(fp, "@ subtitle \"as a function of distance\"\n"); |
422 | } |
423 | xvgr_legend(fp, 2, legr, oenv); |
424 | for (i = 0; i < nrbin; i++) |
425 | { |
426 | fprintf(fp, "%g %g %g\n", (i+0.5)*rbinw, |
427 | histn[i] ? histi1[i]/histn[i] : 0, |
428 | histn[i] ? histi2[i]/histn[i] : 0); |
429 | } |
430 | gmx_ffclose(fp); |
431 | |
432 | sprintf(str, "Cumulative solvent orientation"); |
433 | fp = xvgropen(opt2fn("-co", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), str, "r (nm)", "", oenv); |
434 | if (output_env_get_print_xvgr_codes(oenv)) |
435 | { |
436 | fprintf(fp, "@ subtitle \"as a function of distance\"\n"); |
437 | } |
438 | xvgr_legend(fp, 2, legc, oenv); |
439 | normfac = 1.0/(nrefgrp*nf); |
440 | c1 = 0; |
441 | c2 = 0; |
442 | fprintf(fp, "%g %g %g\n", 0.0, c1, c2); |
443 | for (i = 0; i < nrbin; i++) |
444 | { |
445 | c1 += histi1[i]*normfac; |
446 | c2 += histi2[i]*normfac; |
447 | fprintf(fp, "%g %g %g\n", (i+1)*rbinw, c1, c2); |
448 | } |
449 | gmx_ffclose(fp); |
450 | |
451 | sprintf(str, "Solvent distribution"); |
452 | fp = xvgropen(opt2fn("-rc", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), str, "r (nm)", "molecules/nm", oenv); |
453 | if (output_env_get_print_xvgr_codes(oenv)) |
454 | { |
455 | fprintf(fp, "@ subtitle \"as a function of distance\"\n"); |
456 | } |
457 | normfac = 1.0/(rbinw*nf); |
458 | for (i = 0; i < nrbin; i++) |
459 | { |
460 | fprintf(fp, "%g %g\n", (i+0.5)*rbinw, histn[i]*normfac); |
461 | } |
462 | gmx_ffclose(fp); |
463 | |
464 | do_view(oenv, opt2fn("-o", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), NULL((void*)0)); |
465 | do_view(oenv, opt2fn("-no", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), NULL((void*)0)); |
466 | do_view(oenv, opt2fn("-ro", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), "-nxy"); |
467 | do_view(oenv, opt2fn("-co", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), "-nxy"); |
468 | |
469 | return 0; |
470 | } |