File: | gromacs/gmxana/gmx_hydorder.c |
Location: | line 387, column 29 |
Description: | Array access (from variable 'sk_4d') results in a null pointer dereference |
1 | /* | |||
2 | * This file is part of the GROMACS molecular simulation package. | |||
3 | * | |||
4 | * Copyright (c) 2010,2011,2012,2013,2014, by the GROMACS development team, led by | |||
5 | * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl, | |||
6 | * and including many others, as listed in the AUTHORS file in the | |||
7 | * top-level source directory and at http://www.gromacs.org. | |||
8 | * | |||
9 | * GROMACS is free software; you can redistribute it and/or | |||
10 | * modify it under the terms of the GNU Lesser General Public License | |||
11 | * as published by the Free Software Foundation; either version 2.1 | |||
12 | * of the License, or (at your option) any later version. | |||
13 | * | |||
14 | * GROMACS is distributed in the hope that it will be useful, | |||
15 | * but WITHOUT ANY WARRANTY; without even the implied warranty of | |||
16 | * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU | |||
17 | * Lesser General Public License for more details. | |||
18 | * | |||
19 | * You should have received a copy of the GNU Lesser General Public | |||
20 | * License along with GROMACS; if not, see | |||
21 | * http://www.gnu.org/licenses, or write to the Free Software Foundation, | |||
22 | * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. | |||
23 | * | |||
24 | * If you want to redistribute modifications to GROMACS, please | |||
25 | * consider that scientific software is very special. Version | |||
26 | * control is crucial - bugs must be traceable. We will be happy to | |||
27 | * consider code for inclusion in the official distribution, but | |||
28 | * derived work must not be called official GROMACS. Details are found | |||
29 | * in the README & COPYING files - if they are missing, get the | |||
30 | * official version at http://www.gromacs.org. | |||
31 | * | |||
32 | * To help us fund GROMACS development, we humbly ask that you cite | |||
33 | * the research papers on the package. Check out http://www.gromacs.org. | |||
34 | */ | |||
35 | ||||
36 | #ifdef HAVE_CONFIG_H1 | |||
37 | #include <config.h> | |||
38 | #endif | |||
39 | ||||
40 | #include <math.h> | |||
41 | #include <string.h> | |||
42 | ||||
43 | #include "typedefs.h" | |||
44 | #include "gromacs/commandline/pargs.h" | |||
45 | #include "gromacs/utility/smalloc.h" | |||
46 | #include "macros.h" | |||
47 | #include "gstat.h" | |||
48 | #include "gromacs/math/vec.h" | |||
49 | #include "gromacs/fileio/xvgr.h" | |||
50 | #include "pbc.h" | |||
51 | #include "gromacs/utility/futil.h" | |||
52 | #include "index.h" | |||
53 | #include "gromacs/fileio/tpxio.h" | |||
54 | #include "gromacs/fileio/trxio.h" | |||
55 | #include "gromacs/fileio/matio.h" | |||
56 | #include "binsearch.h" | |||
57 | #include "powerspect.h" | |||
58 | ||||
59 | #include "gromacs/utility/fatalerror.h" | |||
60 | ||||
61 | /* Print name of first atom in all groups in index file */ | |||
62 | static void print_types(atom_id index[], atom_id a[], int ngrps, | |||
63 | char *groups[], t_topology *top) | |||
64 | { | |||
65 | int i; | |||
66 | ||||
67 | fprintf(stderrstderr, "Using following groups: \n"); | |||
68 | for (i = 0; i < ngrps; i++) | |||
69 | { | |||
70 | fprintf(stderrstderr, "Groupname: %s First atomname: %s First atomnr %u\n", | |||
71 | groups[i], *(top->atoms.atomname[a[index[i]]]), a[index[i]]); | |||
72 | } | |||
73 | fprintf(stderrstderr, "\n"); | |||
74 | } | |||
75 | ||||
76 | static void check_length(real length, int a, int b) | |||
77 | { | |||
78 | if (length > 0.3) | |||
79 | { | |||
80 | fprintf(stderrstderr, "WARNING: distance between atoms %d and " | |||
81 | "%d > 0.3 nm (%f). Index file might be corrupt.\n", | |||
82 | a, b, length); | |||
83 | } | |||
84 | } | |||
85 | ||||
86 | static void find_tetra_order_grid(t_topology top, int ePBC, | |||
87 | int natoms, matrix box, | |||
88 | rvec x[], int maxidx, atom_id index[], | |||
89 | real *sgmean, real *skmean, | |||
90 | int nslicex, int nslicey, int nslicez, | |||
91 | real ***sggrid, real ***skgrid) | |||
92 | { | |||
93 | int ix, jx, i, j, k, l, n, *nn[4]; | |||
94 | rvec dx, rj, rk, urk, urj; | |||
95 | real cost, cost2, *sgmol, *skmol, rmean, rmean2, r2, box2, *r_nn[4]; | |||
96 | t_pbc pbc; | |||
97 | int slindex_x, slindex_y, slindex_z; | |||
98 | int ***sl_count; | |||
99 | real onethird = 1.0/3.0; | |||
100 | gmx_rmpbc_t gpbc; | |||
101 | ||||
102 | /* dmat = init_mat(maxidx, FALSE); */ | |||
103 | ||||
104 | box2 = box[XX0][XX0] * box[XX0][XX0]; | |||
105 | ||||
106 | /* Initialize expanded sl_count array */ | |||
107 | snew(sl_count, nslicex)(sl_count) = save_calloc("sl_count", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hydorder.c" , 107, (nslicex), sizeof(*(sl_count))); | |||
108 | for (i = 0; i < nslicex; i++) | |||
109 | { | |||
110 | snew(sl_count[i], nslicey)(sl_count[i]) = save_calloc("sl_count[i]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hydorder.c" , 110, (nslicey), sizeof(*(sl_count[i]))); | |||
111 | for (j = 0; j < nslicey; j++) | |||
112 | { | |||
113 | snew(sl_count[i][j], nslicez)(sl_count[i][j]) = save_calloc("sl_count[i][j]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hydorder.c" , 113, (nslicez), sizeof(*(sl_count[i][j]))); | |||
114 | } | |||
115 | } | |||
116 | ||||
117 | ||||
118 | for (i = 0; (i < 4); i++) | |||
119 | { | |||
120 | snew(r_nn[i], natoms)(r_nn[i]) = save_calloc("r_nn[i]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hydorder.c" , 120, (natoms), sizeof(*(r_nn[i]))); | |||
121 | snew(nn[i], natoms)(nn[i]) = save_calloc("nn[i]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hydorder.c" , 121, (natoms), sizeof(*(nn[i]))); | |||
122 | ||||
123 | for (j = 0; (j < natoms); j++) | |||
124 | { | |||
125 | r_nn[i][j] = box2; | |||
126 | } | |||
127 | } | |||
128 | ||||
129 | snew(sgmol, maxidx)(sgmol) = save_calloc("sgmol", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hydorder.c" , 129, (maxidx), sizeof(*(sgmol))); | |||
130 | snew(skmol, maxidx)(skmol) = save_calloc("skmol", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hydorder.c" , 130, (maxidx), sizeof(*(skmol))); | |||
131 | ||||
132 | /* Must init pbc every step because of pressure coupling */ | |||
133 | set_pbc(&pbc, ePBC, box); | |||
134 | gpbc = gmx_rmpbc_init(&top.idef, ePBC, natoms); | |||
135 | gmx_rmpbc(gpbc, natoms, box, x); | |||
136 | ||||
137 | *sgmean = 0.0; | |||
138 | *skmean = 0.0; | |||
139 | l = 0; | |||
140 | for (i = 0; (i < maxidx); i++) | |||
141 | { /* loop over index file */ | |||
142 | ix = index[i]; | |||
143 | for (j = 0; (j < maxidx); j++) | |||
144 | { | |||
145 | ||||
146 | if (i == j) | |||
147 | { | |||
148 | continue; | |||
149 | } | |||
150 | ||||
151 | jx = index[j]; | |||
152 | ||||
153 | pbc_dx(&pbc, x[ix], x[jx], dx); | |||
154 | r2 = iprod(dx, dx); | |||
155 | ||||
156 | /* set_mat_entry(dmat,i,j,r2); */ | |||
157 | ||||
158 | /* determine the nearest neighbours */ | |||
159 | if (r2 < r_nn[0][i]) | |||
160 | { | |||
161 | r_nn[3][i] = r_nn[2][i]; nn[3][i] = nn[2][i]; | |||
162 | r_nn[2][i] = r_nn[1][i]; nn[2][i] = nn[1][i]; | |||
163 | r_nn[1][i] = r_nn[0][i]; nn[1][i] = nn[0][i]; | |||
164 | r_nn[0][i] = r2; nn[0][i] = j; | |||
165 | } | |||
166 | else if (r2 < r_nn[1][i]) | |||
167 | { | |||
168 | r_nn[3][i] = r_nn[2][i]; nn[3][i] = nn[2][i]; | |||
169 | r_nn[2][i] = r_nn[1][i]; nn[2][i] = nn[1][i]; | |||
170 | r_nn[1][i] = r2; nn[1][i] = j; | |||
171 | } | |||
172 | else if (r2 < r_nn[2][i]) | |||
173 | { | |||
174 | r_nn[3][i] = r_nn[2][i]; nn[3][i] = nn[2][i]; | |||
175 | r_nn[2][i] = r2; nn[2][i] = j; | |||
176 | } | |||
177 | else if (r2 < r_nn[3][i]) | |||
178 | { | |||
179 | r_nn[3][i] = r2; nn[3][i] = j; | |||
180 | } | |||
181 | } | |||
182 | ||||
183 | ||||
184 | /* calculate mean distance between nearest neighbours */ | |||
185 | rmean = 0; | |||
186 | for (j = 0; (j < 4); j++) | |||
187 | { | |||
188 | r_nn[j][i] = sqrt(r_nn[j][i]); | |||
189 | rmean += r_nn[j][i]; | |||
190 | } | |||
191 | rmean /= 4; | |||
192 | ||||
193 | n = 0; | |||
194 | sgmol[i] = 0.0; | |||
195 | skmol[i] = 0.0; | |||
196 | ||||
197 | /* Chau1998a eqn 3 */ | |||
198 | /* angular part tetrahedrality order parameter per atom */ | |||
199 | for (j = 0; (j < 3); j++) | |||
200 | { | |||
201 | for (k = j+1; (k < 4); k++) | |||
202 | { | |||
203 | pbc_dx(&pbc, x[ix], x[index[nn[k][i]]], rk); | |||
204 | pbc_dx(&pbc, x[ix], x[index[nn[j][i]]], rj); | |||
205 | ||||
206 | unitv(rk, urk); | |||
207 | unitv(rj, urj); | |||
208 | ||||
209 | cost = iprod(urk, urj) + onethird; | |||
210 | cost2 = cost * cost; | |||
211 | ||||
212 | sgmol[i] += cost2; | |||
213 | l++; | |||
214 | n++; | |||
215 | } | |||
216 | } | |||
217 | /* normalize sgmol between 0.0 and 1.0 */ | |||
218 | sgmol[i] = 3*sgmol[i]/32; | |||
219 | *sgmean += sgmol[i]; | |||
220 | ||||
221 | /* distance part tetrahedrality order parameter per atom */ | |||
222 | rmean2 = 4 * 3 * rmean * rmean; | |||
223 | for (j = 0; (j < 4); j++) | |||
224 | { | |||
225 | skmol[i] += (rmean - r_nn[j][i]) * (rmean - r_nn[j][i]) / rmean2; | |||
226 | /* printf("%d %f (%f %f %f %f) \n", | |||
227 | i, skmol[i], rmean, rmean2, r_nn[j][i], (rmean - r_nn[j][i]) ); | |||
228 | */ | |||
229 | } | |||
230 | ||||
231 | *skmean += skmol[i]; | |||
232 | ||||
233 | /* Compute sliced stuff in x y z*/ | |||
234 | slindex_x = gmx_nint((1+x[i][XX0]/box[XX0][XX0])*nslicex) % nslicex; | |||
235 | slindex_y = gmx_nint((1+x[i][YY1]/box[YY1][YY1])*nslicey) % nslicey; | |||
236 | slindex_z = gmx_nint((1+x[i][ZZ2]/box[ZZ2][ZZ2])*nslicez) % nslicez; | |||
237 | sggrid[slindex_x][slindex_y][slindex_z] += sgmol[i]; | |||
238 | skgrid[slindex_x][slindex_y][slindex_z] += skmol[i]; | |||
239 | (sl_count[slindex_x][slindex_y][slindex_z])++; | |||
240 | } /* loop over entries in index file */ | |||
241 | ||||
242 | *sgmean /= maxidx; | |||
243 | *skmean /= maxidx; | |||
244 | ||||
245 | for (i = 0; (i < nslicex); i++) | |||
246 | { | |||
247 | for (j = 0; j < nslicey; j++) | |||
248 | { | |||
249 | for (k = 0; k < nslicez; k++) | |||
250 | { | |||
251 | if (sl_count[i][j][k] > 0) | |||
252 | { | |||
253 | sggrid[i][j][k] /= sl_count[i][j][k]; | |||
254 | skgrid[i][j][k] /= sl_count[i][j][k]; | |||
255 | } | |||
256 | } | |||
257 | } | |||
258 | } | |||
259 | ||||
260 | sfree(sl_count)save_free("sl_count", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hydorder.c" , 260, (sl_count)); | |||
261 | sfree(sgmol)save_free("sgmol", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hydorder.c" , 261, (sgmol)); | |||
262 | sfree(skmol)save_free("skmol", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hydorder.c" , 262, (skmol)); | |||
263 | for (i = 0; (i < 4); i++) | |||
264 | { | |||
265 | sfree(r_nn[i])save_free("r_nn[i]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hydorder.c" , 265, (r_nn[i])); | |||
266 | sfree(nn[i])save_free("nn[i]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hydorder.c" , 266, (nn[i])); | |||
267 | } | |||
268 | } | |||
269 | ||||
270 | /*Determines interface from tetrahedral order parameter in box with specified binwidth. */ | |||
271 | /*Outputs interface positions(bins), the number of timeframes, and the number of surface-mesh points in xy*/ | |||
272 | ||||
273 | static void calc_tetra_order_interface(const char *fnNDX, const char *fnTPS, const char *fnTRX, real binw, int tblock, | |||
274 | int *nframes, int *nslicex, int *nslicey, | |||
275 | real sgang1, real sgang2, real ****intfpos, | |||
276 | output_env_t oenv) | |||
277 | { | |||
278 | FILE *fpsg = NULL((void*)0), *fpsk = NULL((void*)0); | |||
279 | char *sgslfn = "sg_ang_mesh"; /* Hardcoded filenames for debugging*/ | |||
280 | char *skslfn = "sk_dist_mesh"; | |||
281 | t_topology top; | |||
282 | int ePBC; | |||
283 | char title[STRLEN4096], subtitle[STRLEN4096]; | |||
284 | t_trxstatus *status; | |||
285 | int natoms; | |||
286 | real t; | |||
287 | rvec *xtop, *x; | |||
288 | matrix box; | |||
289 | real sg, sk, sgintf, pos; | |||
290 | atom_id **index = NULL((void*)0); | |||
291 | char **grpname = NULL((void*)0); | |||
292 | int i, j, k, n, *isize, ng, nslicez, framenr; | |||
293 | real ***sg_grid = NULL((void*)0), ***sk_grid = NULL((void*)0), ***sg_fravg = NULL((void*)0), ***sk_fravg = NULL((void*)0), ****sk_4d = NULL((void*)0), ****sg_4d = NULL((void*)0); | |||
294 | int *perm; | |||
295 | int ndx1, ndx2; | |||
296 | int bins; | |||
297 | const real onehalf = 1.0/2.0; | |||
298 | /* real ***intfpos[2]; pointers to arrays of two interface positions zcoord(framenr,xbin,ybin): intfpos[interface_index][t][nslicey*x+y] | |||
299 | * i.e 1D Row-major order in (t,x,y) */ | |||
300 | ||||
301 | ||||
302 | read_tps_conf(fnTPS, title, &top, &ePBC, &xtop, NULL((void*)0), box, FALSE0); | |||
303 | ||||
304 | *nslicex = (int)(box[XX0][XX0]/binw + onehalf); /*Calculate slicenr from binwidth*/ | |||
305 | *nslicey = (int)(box[YY1][YY1]/binw + onehalf); | |||
306 | nslicez = (int)(box[ZZ2][ZZ2]/binw + onehalf); | |||
307 | ||||
308 | ||||
309 | ||||
310 | ng = 1; | |||
311 | /* get index groups */ | |||
312 | printf("Select the group that contains the atoms you want to use for the tetrahedrality order parameter calculation:\n"); | |||
313 | snew(grpname, ng)(grpname) = save_calloc("grpname", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hydorder.c" , 313, (ng), sizeof(*(grpname))); | |||
314 | snew(index, ng)(index) = save_calloc("index", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hydorder.c" , 314, (ng), sizeof(*(index))); | |||
315 | snew(isize, ng)(isize) = save_calloc("isize", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hydorder.c" , 315, (ng), sizeof(*(isize))); | |||
316 | get_index(&top.atoms, fnNDX, ng, isize, index, grpname); | |||
317 | ||||
318 | /* Analyze trajectory */ | |||
319 | natoms = read_first_x(oenv, &status, fnTRX, &t, &x, box); | |||
320 | if (natoms > top.atoms.nr) | |||
| ||||
321 | { | |||
322 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hydorder.c" , 322, "Topology (%d atoms) does not match trajectory (%d atoms)", | |||
323 | top.atoms.nr, natoms); | |||
324 | } | |||
325 | check_index(NULL((void*)0), ng, index[0], NULL((void*)0), natoms); | |||
326 | ||||
327 | ||||
328 | /*Prepare structures for temporary storage of frame info*/ | |||
329 | snew(sg_grid, *nslicex)(sg_grid) = save_calloc("sg_grid", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hydorder.c" , 329, (*nslicex), sizeof(*(sg_grid))); | |||
330 | snew(sk_grid, *nslicex)(sk_grid) = save_calloc("sk_grid", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hydorder.c" , 330, (*nslicex), sizeof(*(sk_grid))); | |||
331 | for (i = 0; i < *nslicex; i++) | |||
332 | { | |||
333 | snew(sg_grid[i], *nslicey)(sg_grid[i]) = save_calloc("sg_grid[i]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hydorder.c" , 333, (*nslicey), sizeof(*(sg_grid[i]))); | |||
334 | snew(sk_grid[i], *nslicey)(sk_grid[i]) = save_calloc("sk_grid[i]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hydorder.c" , 334, (*nslicey), sizeof(*(sk_grid[i]))); | |||
335 | for (j = 0; j < *nslicey; j++) | |||
336 | { | |||
337 | snew(sg_grid[i][j], nslicez)(sg_grid[i][j]) = save_calloc("sg_grid[i][j]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hydorder.c" , 337, (nslicez), sizeof(*(sg_grid[i][j]))); | |||
338 | snew(sk_grid[i][j], nslicez)(sk_grid[i][j]) = save_calloc("sk_grid[i][j]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hydorder.c" , 338, (nslicez), sizeof(*(sk_grid[i][j]))); | |||
339 | } | |||
340 | } | |||
341 | ||||
342 | sg_4d = NULL((void*)0); | |||
343 | sk_4d = NULL((void*)0); | |||
344 | *nframes = 0; | |||
345 | framenr = 0; | |||
346 | ||||
347 | /* Loop over frames*/ | |||
348 | do | |||
349 | { | |||
350 | /*Initialize box meshes (temporary storage for each tblock frame -reinitialise every tblock steps */ | |||
351 | if (framenr%tblock == 0) | |||
352 | { | |||
353 | srenew(sk_4d, *nframes+1)(sk_4d) = save_realloc("sk_4d", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hydorder.c" , 353, (sk_4d), (*nframes+1), sizeof(*(sk_4d))); | |||
354 | srenew(sg_4d, *nframes+1)(sg_4d) = save_realloc("sg_4d", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hydorder.c" , 354, (sg_4d), (*nframes+1), sizeof(*(sg_4d))); | |||
355 | snew(sg_fravg, *nslicex)(sg_fravg) = save_calloc("sg_fravg", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hydorder.c" , 355, (*nslicex), sizeof(*(sg_fravg))); | |||
356 | snew(sk_fravg, *nslicex)(sk_fravg) = save_calloc("sk_fravg", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hydorder.c" , 356, (*nslicex), sizeof(*(sk_fravg))); | |||
357 | for (i = 0; i < *nslicex; i++) | |||
358 | { | |||
359 | snew(sg_fravg[i], *nslicey)(sg_fravg[i]) = save_calloc("sg_fravg[i]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hydorder.c" , 359, (*nslicey), sizeof(*(sg_fravg[i]))); | |||
360 | snew(sk_fravg[i], *nslicey)(sk_fravg[i]) = save_calloc("sk_fravg[i]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hydorder.c" , 360, (*nslicey), sizeof(*(sk_fravg[i]))); | |||
361 | for (j = 0; j < *nslicey; j++) | |||
362 | { | |||
363 | snew(sg_fravg[i][j], nslicez)(sg_fravg[i][j]) = save_calloc("sg_fravg[i][j]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hydorder.c" , 363, (nslicez), sizeof(*(sg_fravg[i][j]))); | |||
364 | snew(sk_fravg[i][j], nslicez)(sk_fravg[i][j]) = save_calloc("sk_fravg[i][j]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hydorder.c" , 364, (nslicez), sizeof(*(sk_fravg[i][j]))); | |||
365 | } | |||
366 | } | |||
367 | } | |||
368 | ||||
369 | find_tetra_order_grid(top, ePBC, natoms, box, x, isize[0], index[0], | |||
370 | &sg, &sk, *nslicex, *nslicey, nslicez, sg_grid, sk_grid); | |||
371 | for (i = 0; i < *nslicex; i++) | |||
372 | { | |||
373 | for (j = 0; j < *nslicey; j++) | |||
374 | { | |||
375 | for (k = 0; k < nslicez; k++) | |||
376 | { | |||
377 | sk_fravg[i][j][k] += sk_grid[i][j][k]/tblock; | |||
378 | sg_fravg[i][j][k] += sg_grid[i][j][k]/tblock; | |||
379 | } | |||
380 | } | |||
381 | } | |||
382 | ||||
383 | framenr++; | |||
384 | ||||
385 | if (framenr%tblock == 0) | |||
386 | { | |||
387 | sk_4d[*nframes] = sk_fravg; | |||
| ||||
388 | sg_4d[*nframes] = sg_fravg; | |||
389 | (*nframes)++; | |||
390 | } | |||
391 | ||||
392 | } | |||
393 | while (read_next_x(oenv, status, &t, x, box)); | |||
394 | close_trj(status); | |||
395 | ||||
396 | sfree(grpname)save_free("grpname", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hydorder.c" , 396, (grpname)); | |||
397 | sfree(index)save_free("index", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hydorder.c" , 397, (index)); | |||
398 | sfree(isize)save_free("isize", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hydorder.c" , 398, (isize)); | |||
399 | ||||
400 | /*Debugging for printing out the entire order parameter meshes.*/ | |||
401 | if (debug) | |||
402 | { | |||
403 | fpsg = xvgropen(sgslfn, "S\\sg\\N Angle Order Parameter / Meshpoint", "(nm)", "S\\sg\\N", oenv); | |||
404 | fpsk = xvgropen(skslfn, "S\\sk\\N Distance Order Parameter / Meshpoint", "(nm)", "S\\sk\\N", oenv); | |||
405 | for (n = 0; n < (*nframes); n++) | |||
406 | { | |||
407 | fprintf(fpsg, "%i\n", n); | |||
408 | fprintf(fpsk, "%i\n", n); | |||
409 | for (i = 0; (i < *nslicex); i++) | |||
410 | { | |||
411 | for (j = 0; j < *nslicey; j++) | |||
412 | { | |||
413 | for (k = 0; k < nslicez; k++) | |||
414 | { | |||
415 | fprintf(fpsg, "%4f %4f %4f %8f\n", (i+0.5)*box[XX0][XX0]/(*nslicex), (j+0.5)*box[YY1][YY1]/(*nslicey), (k+0.5)*box[ZZ2][ZZ2]/nslicez, sg_4d[n][i][j][k]); | |||
416 | fprintf(fpsk, "%4f %4f %4f %8f\n", (i+0.5)*box[XX0][XX0]/(*nslicex), (j+0.5)*box[YY1][YY1]/(*nslicey), (k+0.5)*box[ZZ2][ZZ2]/nslicez, sk_4d[n][i][j][k]); | |||
417 | } | |||
418 | } | |||
419 | } | |||
420 | } | |||
421 | xvgrclose(fpsg); | |||
422 | xvgrclose(fpsk); | |||
423 | } | |||
424 | ||||
425 | ||||
426 | /* Find positions of interface z by scanning orderparam for each frame and for each xy-mesh cylinder along z*/ | |||
427 | ||||
428 | /*Simple trial: assume interface is in the middle of -sgang1 and sgang2*/ | |||
429 | sgintf = 0.5*(sgang1+sgang2); | |||
430 | ||||
431 | ||||
432 | /*Allocate memory for interface arrays; */ | |||
433 | snew((*intfpos), 2)((*intfpos)) = save_calloc("(*intfpos)", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hydorder.c" , 433, (2), sizeof(*((*intfpos)))); | |||
434 | snew((*intfpos)[0], *nframes)((*intfpos)[0]) = save_calloc("(*intfpos)[0]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hydorder.c" , 434, (*nframes), sizeof(*((*intfpos)[0]))); | |||
435 | snew((*intfpos)[1], *nframes)((*intfpos)[1]) = save_calloc("(*intfpos)[1]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hydorder.c" , 435, (*nframes), sizeof(*((*intfpos)[1]))); | |||
436 | ||||
437 | bins = (*nslicex)*(*nslicey); | |||
438 | ||||
439 | ||||
440 | snew(perm, nslicez)(perm) = save_calloc("perm", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hydorder.c" , 440, (nslicez), sizeof(*(perm))); /*permutation array for sorting along normal coordinate*/ | |||
441 | ||||
442 | ||||
443 | for (n = 0; n < *nframes; n++) | |||
444 | { | |||
445 | snew((*intfpos)[0][n], bins)((*intfpos)[0][n]) = save_calloc("(*intfpos)[0][n]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hydorder.c" , 445, (bins), sizeof(*((*intfpos)[0][n]))); | |||
446 | snew((*intfpos)[1][n], bins)((*intfpos)[1][n]) = save_calloc("(*intfpos)[1][n]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hydorder.c" , 446, (bins), sizeof(*((*intfpos)[1][n]))); | |||
447 | for (i = 0; i < *nslicex; i++) | |||
448 | { | |||
449 | for (j = 0; j < *nslicey; j++) | |||
450 | { | |||
451 | rangeArray(perm, nslicez); /*reset permutation array to identity*/ | |||
452 | /*Binsearch returns 2 bin-numbers where the order param is <= setpoint sgintf*/ | |||
453 | ndx1 = start_binsearch(sg_4d[n][i][j], perm, 0, nslicez/2-1, sgintf, 1); | |||
454 | ndx2 = start_binsearch(sg_4d[n][i][j], perm, nslicez/2, nslicez-1, sgintf, -1); | |||
455 | /*Use linear interpolation to smooth out the interface position*/ | |||
456 | ||||
457 | /*left interface (0)*/ | |||
458 | /*if((sg_4d[n][i][j][perm[ndx1+1]]-sg_4d[n][i][j][perm[ndx1]])/sg_4d[n][i][j][perm[ndx1]] > 0.01){ | |||
459 | pos=( (sgintf-sg_4d[n][i][j][perm[ndx1]])*perm[ndx1+1]+(sg_4d[n][i][j][perm[ndx1+1]]-sgintf)*perm[ndx1 ])*/ | |||
460 | (*intfpos)[0][n][j+*nslicey*i] = (perm[ndx1]+onehalf)*binw; | |||
461 | /*right interface (1)*/ | |||
462 | /*alpha=(sgintf-sg_4d[n][i][j][perm[ndx2]])/(sg_4d[n][i][j][perm[ndx2]+1]-sg_4d[n][i][j][perm[ndx2]]);*/ | |||
463 | /*(*intfpos)[1][n][j+*nslicey*i]=((1-alpha)*perm[ndx2]+alpha*(perm[ndx2]+1)+onehalf)*box[ZZ][ZZ]/nslicez;*/ | |||
464 | (*intfpos)[1][n][j+*nslicey*i] = (perm[ndx2]+onehalf)*binw; | |||
465 | } | |||
466 | } | |||
467 | } | |||
468 | ||||
469 | ||||
470 | /*sfree(perm);*/ | |||
471 | sfree(sk_4d)save_free("sk_4d", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hydorder.c" , 471, (sk_4d)); | |||
472 | sfree(sg_4d)save_free("sg_4d", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hydorder.c" , 472, (sg_4d)); | |||
473 | /*sfree(sg_grid);*/ | |||
474 | /*sfree(sk_grid);*/ | |||
475 | ||||
476 | ||||
477 | } | |||
478 | ||||
479 | static void writesurftoxpms(real ***surf, int tblocks, int xbins, int ybins, real bw, char **outfiles, int maplevels ) | |||
480 | { | |||
481 | ||||
482 | char numbuf[8]; | |||
483 | int n, i, j; | |||
484 | real **profile1, **profile2; | |||
485 | real max1, max2, min1, min2, *xticks, *yticks; | |||
486 | t_rgb lo = {1, 1, 1}; | |||
487 | t_rgb hi = {0, 0, 0}; | |||
488 | FILE *xpmfile1, *xpmfile2; | |||
489 | ||||
490 | /*Prepare xpm structures for output*/ | |||
491 | ||||
492 | /*Allocate memory to tick's and matrices*/ | |||
493 | snew (xticks, xbins+1)(xticks) = save_calloc("xticks", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hydorder.c" , 493, (xbins+1), sizeof(*(xticks))); | |||
494 | snew (yticks, ybins+1)(yticks) = save_calloc("yticks", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hydorder.c" , 494, (ybins+1), sizeof(*(yticks))); | |||
495 | ||||
496 | profile1 = mk_matrix(xbins, ybins, FALSE0); | |||
497 | profile2 = mk_matrix(xbins, ybins, FALSE0); | |||
498 | ||||
499 | for (i = 0; i < xbins+1; i++) | |||
500 | { | |||
501 | xticks[i] += bw; | |||
502 | } | |||
503 | for (j = 0; j < ybins+1; j++) | |||
504 | { | |||
505 | yticks[j] += bw; | |||
506 | } | |||
507 | ||||
508 | xpmfile1 = gmx_ffopen(outfiles[0], "w"); | |||
509 | xpmfile2 = gmx_ffopen(outfiles[1], "w"); | |||
510 | ||||
511 | max1 = max2 = 0.0; | |||
512 | min1 = min2 = 1000.00; | |||
513 | ||||
514 | for (n = 0; n < tblocks; n++) | |||
515 | { | |||
516 | sprintf(numbuf, "%5d", n); | |||
517 | /*Filling matrices for inclusion in xpm-files*/ | |||
518 | for (i = 0; i < xbins; i++) | |||
519 | { | |||
520 | for (j = 0; j < ybins; j++) | |||
521 | { | |||
522 | profile1[i][j] = (surf[0][n][j+ybins*i]); | |||
523 | profile2[i][j] = (surf[1][n][j+ybins*i]); | |||
524 | /*Finding max and min values*/ | |||
525 | if (profile1[i][j] > max1) | |||
526 | { | |||
527 | max1 = profile1[i][j]; | |||
528 | } | |||
529 | if (profile1[i][j] < min1) | |||
530 | { | |||
531 | min1 = profile1[i][j]; | |||
532 | } | |||
533 | if (profile2[i][j] > max2) | |||
534 | { | |||
535 | max2 = profile2[i][j]; | |||
536 | } | |||
537 | if (profile2[i][j] < min2) | |||
538 | { | |||
539 | min2 = profile2[i][j]; | |||
540 | } | |||
541 | } | |||
542 | } | |||
543 | ||||
544 | write_xpm(xpmfile1, 3, numbuf, "Height", "x[nm]", "y[nm]", xbins, ybins, xticks, yticks, profile1, min1, max1, lo, hi, &maplevels); | |||
545 | write_xpm(xpmfile2, 3, numbuf, "Height", "x[nm]", "y[nm]", xbins, ybins, xticks, yticks, profile2, min2, max2, lo, hi, &maplevels); | |||
546 | } | |||
547 | ||||
548 | gmx_ffclose(xpmfile1); | |||
549 | gmx_ffclose(xpmfile2); | |||
550 | ||||
551 | ||||
552 | ||||
553 | sfree(profile1)save_free("profile1", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hydorder.c" , 553, (profile1)); | |||
554 | sfree(profile2)save_free("profile2", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hydorder.c" , 554, (profile2)); | |||
555 | sfree(xticks)save_free("xticks", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hydorder.c" , 555, (xticks)); | |||
556 | sfree(yticks)save_free("yticks", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hydorder.c" , 556, (yticks)); | |||
557 | } | |||
558 | ||||
559 | ||||
560 | static void writeraw(real ***surf, int tblocks, int xbins, int ybins, char **fnms) | |||
561 | { | |||
562 | FILE *raw1, *raw2; | |||
563 | int i, j, n; | |||
564 | ||||
565 | raw1 = gmx_ffopen(fnms[0], "w"); | |||
566 | raw2 = gmx_ffopen(fnms[1], "w"); | |||
567 | fprintf(raw1, "#Legend\n#TBlock\n#Xbin Ybin Z t\n"); | |||
568 | fprintf(raw2, "#Legend\n#TBlock\n#Xbin Ybin Z t\n"); | |||
569 | for (n = 0; n < tblocks; n++) | |||
570 | { | |||
571 | fprintf(raw1, "%5d\n", n); | |||
572 | fprintf(raw2, "%5d\n", n); | |||
573 | for (i = 0; i < xbins; i++) | |||
574 | { | |||
575 | for (j = 0; j < ybins; j++) | |||
576 | { | |||
577 | fprintf(raw1, "%i %i %8.5f\n", i, j, (surf[0][n][j+ybins*i])); | |||
578 | fprintf(raw2, "%i %i %8.5f\n", i, j, (surf[1][n][j+ybins*i])); | |||
579 | } | |||
580 | } | |||
581 | } | |||
582 | ||||
583 | gmx_ffclose(raw1); | |||
584 | gmx_ffclose(raw2); | |||
585 | } | |||
586 | ||||
587 | ||||
588 | ||||
589 | int gmx_hydorder(int argc, char *argv[]) | |||
590 | { | |||
591 | static const char *desc[] = { | |||
592 | "[THISMODULE] computes the tetrahedrality order parameters around a ", | |||
593 | "given atom. Both angle an distance order parameters are calculated. See", | |||
594 | "P.-L. Chau and A.J. Hardwick, Mol. Phys., 93, (1998), 511-518.", | |||
595 | "for more details.[PAR]" | |||
596 | "[THISMODULE] calculates the order parameter in a 3d-mesh in the box, and", | |||
597 | "with 2 phases in the box gives the user the option to define a 2D interface in time", | |||
598 | "separating the faces by specifying parameters [TT]-sgang1[tt] and", | |||
599 | "[TT]-sgang2[tt] (it is important to select these judiciously)." | |||
600 | }; | |||
601 | ||||
602 | int axis = 0; | |||
603 | static int nsttblock = 1; | |||
604 | static int nlevels = 100; | |||
605 | static real binwidth = 1.0; /* binwidth in mesh */ | |||
606 | static real sg1 = 1; | |||
607 | static real sg2 = 1; /* order parameters for bulk phases */ | |||
608 | static gmx_bool bFourier = FALSE0; | |||
609 | static gmx_bool bRawOut = FALSE0; | |||
610 | int frames, xslices, yslices; /* Dimensions of interface arrays*/ | |||
611 | real ***intfpos; /* Interface arrays (intfnr,t,xy) -potentially large */ | |||
612 | static char *normal_axis[] = { NULL((void*)0), "z", "x", "y", NULL((void*)0) }; | |||
613 | ||||
614 | t_pargs pa[] = { | |||
615 | { "-d", FALSE0, etENUM, {normal_axis}, | |||
616 | "Direction of the normal on the membrane" }, | |||
617 | { "-bw", FALSE0, etREAL, {&binwidth}, | |||
618 | "Binwidth of box mesh" }, | |||
619 | { "-sgang1", FALSE0, etREAL, {&sg1}, | |||
620 | "tetrahedral angle parameter in Phase 1 (bulk)" }, | |||
621 | { "-sgang2", FALSE0, etREAL, {&sg2}, | |||
622 | "tetrahedral angle parameter in Phase 2 (bulk)" }, | |||
623 | { "-tblock", FALSE0, etINT, {&nsttblock}, | |||
624 | "Number of frames in one time-block average"}, | |||
625 | { "-nlevel", FALSE0, etINT, {&nlevels}, | |||
626 | "Number of Height levels in 2D - XPixMaps"} | |||
627 | }; | |||
628 | ||||
629 | t_filenm fnm[] = { /* files for g_order */ | |||
630 | { efTRX, "-f", NULL((void*)0), ffREAD1<<1 }, /* trajectory file */ | |||
631 | { efNDX, "-n", NULL((void*)0), ffREAD1<<1 }, /* index file */ | |||
632 | { efTPX, "-s", NULL((void*)0), ffREAD1<<1 }, /* topology file */ | |||
633 | { efXPM, "-o", "intf", ffWRMULT(1<<2 | 1<<5)}, /* XPM- surface maps */ | |||
634 | { efOUT, "-or", "raw", ffOPTWRMULT((1<<2 | 1<<5) | 1<<3) }, /* xvgr output file */ | |||
635 | { efOUT, "-Spect", "intfspect", ffOPTWRMULT((1<<2 | 1<<5) | 1<<3)}, /* Fourier spectrum interfaces */ | |||
636 | }; | |||
637 | #define NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))) asize(fnm)((int)(sizeof(fnm)/sizeof((fnm)[0]))) | |||
638 | ||||
639 | /*Filenames*/ | |||
640 | const char *ndxfnm, *tpsfnm, *trxfnm; | |||
641 | char **spectra, **intfn, **raw; | |||
642 | int nfspect, nfxpm, nfraw; | |||
643 | output_env_t oenv; | |||
644 | ||||
645 | if (!parse_common_args(&argc, argv, PCA_CAN_VIEW(1<<5) | PCA_CAN_TIME((1<<6) | (1<<7) | (1<<14)) | PCA_BE_NICE(1<<13), | |||
646 | 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)) | |||
647 | { | |||
648 | return 0; | |||
649 | } | |||
650 | bFourier = opt2bSet("-Spect", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm); | |||
651 | bRawOut = opt2bSet("-or", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm); | |||
652 | ||||
653 | if (binwidth < 0.0) | |||
654 | { | |||
655 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hydorder.c" , 655, "Can not have binwidth < 0"); | |||
656 | } | |||
657 | ||||
658 | ndxfnm = ftp2fn(efNDX, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm); | |||
659 | tpsfnm = ftp2fn(efTPX, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm); | |||
660 | trxfnm = ftp2fn(efTRX, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm); | |||
661 | ||||
662 | /* Calculate axis */ | |||
663 | if (strcmp(normal_axis[0], "x")__extension__ ({ size_t __s1_len, __s2_len; (__builtin_constant_p (normal_axis[0]) && __builtin_constant_p ("x") && (__s1_len = strlen (normal_axis[0]), __s2_len = strlen ("x") , (!((size_t)(const void *)((normal_axis[0]) + 1) - (size_t)( const void *)(normal_axis[0]) == 1) || __s1_len >= 4) && (!((size_t)(const void *)(("x") + 1) - (size_t)(const void * )("x") == 1) || __s2_len >= 4)) ? __builtin_strcmp (normal_axis [0], "x") : (__builtin_constant_p (normal_axis[0]) && ((size_t)(const void *)((normal_axis[0]) + 1) - (size_t)(const void *)(normal_axis[0]) == 1) && (__s1_len = strlen ( normal_axis[0]), __s1_len < 4) ? (__builtin_constant_p ("x" ) && ((size_t)(const void *)(("x") + 1) - (size_t)(const void *)("x") == 1) ? __builtin_strcmp (normal_axis[0], "x") : (__extension__ ({ const unsigned char *__s2 = (const unsigned char *) (const char *) ("x"); int __result = (((const unsigned char *) (const char *) (normal_axis[0]))[0] - __s2[0]); if ( __s1_len > 0 && __result == 0) { __result = (((const unsigned char *) (const char *) (normal_axis[0]))[1] - __s2[ 1]); if (__s1_len > 1 && __result == 0) { __result = (((const unsigned char *) (const char *) (normal_axis[0])) [2] - __s2[2]); if (__s1_len > 2 && __result == 0) __result = (((const unsigned char *) (const char *) (normal_axis [0]))[3] - __s2[3]); } } __result; }))) : (__builtin_constant_p ("x") && ((size_t)(const void *)(("x") + 1) - (size_t )(const void *)("x") == 1) && (__s2_len = strlen ("x" ), __s2_len < 4) ? (__builtin_constant_p (normal_axis[0]) && ((size_t)(const void *)((normal_axis[0]) + 1) - (size_t)(const void *)(normal_axis[0]) == 1) ? __builtin_strcmp (normal_axis [0], "x") : (- (__extension__ ({ const unsigned char *__s2 = ( const unsigned char *) (const char *) (normal_axis[0]); int __result = (((const unsigned char *) (const char *) ("x"))[0] - __s2[ 0]); if (__s2_len > 0 && __result == 0) { __result = (((const unsigned char *) (const char *) ("x"))[1] - __s2[ 1]); if (__s2_len > 1 && __result == 0) { __result = (((const unsigned char *) (const char *) ("x"))[2] - __s2[ 2]); if (__s2_len > 2 && __result == 0) __result = (((const unsigned char *) (const char *) ("x"))[3] - __s2[3] ); } } __result; })))) : __builtin_strcmp (normal_axis[0], "x" )))); }) == 0) | |||
664 | { | |||
665 | axis = XX0; | |||
666 | } | |||
667 | else if (strcmp(normal_axis[0], "y")__extension__ ({ size_t __s1_len, __s2_len; (__builtin_constant_p (normal_axis[0]) && __builtin_constant_p ("y") && (__s1_len = strlen (normal_axis[0]), __s2_len = strlen ("y") , (!((size_t)(const void *)((normal_axis[0]) + 1) - (size_t)( const void *)(normal_axis[0]) == 1) || __s1_len >= 4) && (!((size_t)(const void *)(("y") + 1) - (size_t)(const void * )("y") == 1) || __s2_len >= 4)) ? __builtin_strcmp (normal_axis [0], "y") : (__builtin_constant_p (normal_axis[0]) && ((size_t)(const void *)((normal_axis[0]) + 1) - (size_t)(const void *)(normal_axis[0]) == 1) && (__s1_len = strlen ( normal_axis[0]), __s1_len < 4) ? (__builtin_constant_p ("y" ) && ((size_t)(const void *)(("y") + 1) - (size_t)(const void *)("y") == 1) ? __builtin_strcmp (normal_axis[0], "y") : (__extension__ ({ const unsigned char *__s2 = (const unsigned char *) (const char *) ("y"); int __result = (((const unsigned char *) (const char *) (normal_axis[0]))[0] - __s2[0]); if ( __s1_len > 0 && __result == 0) { __result = (((const unsigned char *) (const char *) (normal_axis[0]))[1] - __s2[ 1]); if (__s1_len > 1 && __result == 0) { __result = (((const unsigned char *) (const char *) (normal_axis[0])) [2] - __s2[2]); if (__s1_len > 2 && __result == 0) __result = (((const unsigned char *) (const char *) (normal_axis [0]))[3] - __s2[3]); } } __result; }))) : (__builtin_constant_p ("y") && ((size_t)(const void *)(("y") + 1) - (size_t )(const void *)("y") == 1) && (__s2_len = strlen ("y" ), __s2_len < 4) ? (__builtin_constant_p (normal_axis[0]) && ((size_t)(const void *)((normal_axis[0]) + 1) - (size_t)(const void *)(normal_axis[0]) == 1) ? __builtin_strcmp (normal_axis [0], "y") : (- (__extension__ ({ const unsigned char *__s2 = ( const unsigned char *) (const char *) (normal_axis[0]); int __result = (((const unsigned char *) (const char *) ("y"))[0] - __s2[ 0]); if (__s2_len > 0 && __result == 0) { __result = (((const unsigned char *) (const char *) ("y"))[1] - __s2[ 1]); if (__s2_len > 1 && __result == 0) { __result = (((const unsigned char *) (const char *) ("y"))[2] - __s2[ 2]); if (__s2_len > 2 && __result == 0) __result = (((const unsigned char *) (const char *) ("y"))[3] - __s2[3] ); } } __result; })))) : __builtin_strcmp (normal_axis[0], "y" )))); }) == 0) | |||
668 | { | |||
669 | axis = YY1; | |||
670 | } | |||
671 | else if (strcmp(normal_axis[0], "z")__extension__ ({ size_t __s1_len, __s2_len; (__builtin_constant_p (normal_axis[0]) && __builtin_constant_p ("z") && (__s1_len = strlen (normal_axis[0]), __s2_len = strlen ("z") , (!((size_t)(const void *)((normal_axis[0]) + 1) - (size_t)( const void *)(normal_axis[0]) == 1) || __s1_len >= 4) && (!((size_t)(const void *)(("z") + 1) - (size_t)(const void * )("z") == 1) || __s2_len >= 4)) ? __builtin_strcmp (normal_axis [0], "z") : (__builtin_constant_p (normal_axis[0]) && ((size_t)(const void *)((normal_axis[0]) + 1) - (size_t)(const void *)(normal_axis[0]) == 1) && (__s1_len = strlen ( normal_axis[0]), __s1_len < 4) ? (__builtin_constant_p ("z" ) && ((size_t)(const void *)(("z") + 1) - (size_t)(const void *)("z") == 1) ? __builtin_strcmp (normal_axis[0], "z") : (__extension__ ({ const unsigned char *__s2 = (const unsigned char *) (const char *) ("z"); int __result = (((const unsigned char *) (const char *) (normal_axis[0]))[0] - __s2[0]); if ( __s1_len > 0 && __result == 0) { __result = (((const unsigned char *) (const char *) (normal_axis[0]))[1] - __s2[ 1]); if (__s1_len > 1 && __result == 0) { __result = (((const unsigned char *) (const char *) (normal_axis[0])) [2] - __s2[2]); if (__s1_len > 2 && __result == 0) __result = (((const unsigned char *) (const char *) (normal_axis [0]))[3] - __s2[3]); } } __result; }))) : (__builtin_constant_p ("z") && ((size_t)(const void *)(("z") + 1) - (size_t )(const void *)("z") == 1) && (__s2_len = strlen ("z" ), __s2_len < 4) ? (__builtin_constant_p (normal_axis[0]) && ((size_t)(const void *)((normal_axis[0]) + 1) - (size_t)(const void *)(normal_axis[0]) == 1) ? __builtin_strcmp (normal_axis [0], "z") : (- (__extension__ ({ const unsigned char *__s2 = ( const unsigned char *) (const char *) (normal_axis[0]); int __result = (((const unsigned char *) (const char *) ("z"))[0] - __s2[ 0]); if (__s2_len > 0 && __result == 0) { __result = (((const unsigned char *) (const char *) ("z"))[1] - __s2[ 1]); if (__s2_len > 1 && __result == 0) { __result = (((const unsigned char *) (const char *) ("z"))[2] - __s2[ 2]); if (__s2_len > 2 && __result == 0) __result = (((const unsigned char *) (const char *) ("z"))[3] - __s2[3] ); } } __result; })))) : __builtin_strcmp (normal_axis[0], "z" )))); }) == 0) | |||
672 | { | |||
673 | axis = ZZ2; | |||
674 | } | |||
675 | else | |||
676 | { | |||
677 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hydorder.c" , 677, "Invalid axis, use x, y or z"); | |||
678 | } | |||
679 | ||||
680 | switch (axis) | |||
681 | { | |||
682 | case 0: | |||
683 | fprintf(stderrstderr, "Taking x axis as normal to the membrane\n"); | |||
684 | break; | |||
685 | case 1: | |||
686 | fprintf(stderrstderr, "Taking y axis as normal to the membrane\n"); | |||
687 | break; | |||
688 | case 2: | |||
689 | fprintf(stderrstderr, "Taking z axis as normal to the membrane\n"); | |||
690 | break; | |||
691 | } | |||
692 | ||||
693 | /* tetraheder order parameter */ | |||
694 | /* If either of the options is set we compute both */ | |||
695 | nfxpm = opt2fns(&intfn, "-o", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm); | |||
696 | if (nfxpm != 2) | |||
697 | { | |||
698 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hydorder.c" , 698, "No or not correct number (2) of output-files: %d", nfxpm); | |||
699 | } | |||
700 | calc_tetra_order_interface(ndxfnm, tpsfnm, trxfnm, binwidth, nsttblock, &frames, &xslices, &yslices, sg1, sg2, &intfpos, oenv); | |||
701 | writesurftoxpms(intfpos, frames, xslices, yslices, binwidth, intfn, nlevels); | |||
702 | ||||
703 | if (bFourier) | |||
704 | { | |||
705 | nfspect = opt2fns(&spectra, "-Spect", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm); | |||
706 | if (nfspect != 2) | |||
707 | { | |||
708 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hydorder.c" , 708, "No or not correct number (2) of output-files: %d", nfspect); | |||
709 | } | |||
710 | powerspectavg(intfpos, frames, xslices, yslices, spectra); | |||
711 | } | |||
712 | ||||
713 | if (bRawOut) | |||
714 | { | |||
715 | nfraw = opt2fns(&raw, "-or", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm); | |||
716 | if (nfraw != 2) | |||
717 | { | |||
718 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hydorder.c" , 718, "No or not correct number (2) of output-files: %d", nfraw); | |||
719 | } | |||
720 | writeraw(intfpos, frames, xslices, yslices, raw); | |||
721 | } | |||
722 | ||||
723 | return 0; | |||
724 | } |