Bug Summary

File:gromacs/gmxana/gmx_order.c
Location:line 853, column 5
Description:Value stored to 'natoms' is never read

Annotated Source Code

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 <math.h>
42#include <string.h>
43
44#include "typedefs.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 "viewit.h"
51#include "pbc.h"
52#include "gromacs/utility/futil.h"
53#include "gromacs/commandline/pargs.h"
54#include "index.h"
55#include "gromacs/fileio/tpxio.h"
56#include "gromacs/fileio/trxio.h"
57#include "gromacs/fileio/confio.h"
58#include "cmat.h"
59#include "gmx_ana.h"
60
61#include "gromacs/utility/fatalerror.h"
62
63/****************************************************************************/
64/* This program calculates the order parameter per atom for an interface or */
65/* bilayer, averaged over time. */
66/* S = 1/2 * (3 * cos(i)cos(j) - delta(ij)) */
67/* It is assumed that the order parameter with respect to a box-axis */
68/* is calculated. In that case i = j = axis, and delta(ij) = 1. */
69/* */
70/* Peter Tieleman, April 1995 */
71/* P.J. van Maaren, November 2005 Added tetrahedral stuff */
72/****************************************************************************/
73
74static void find_nearest_neighbours(int ePBC,
75 int natoms, matrix box,
76 rvec x[], int maxidx, atom_id index[],
77 real *sgmean, real *skmean,
78 int nslice, int slice_dim,
79 real sgslice[], real skslice[],
80 gmx_rmpbc_t gpbc)
81{
82 FILE *fpoutdist;
83 char fnsgdist[32];
84 int ix, jx, nsgbin, *sgbin;
85 int i1, i2, i, ibin, j, k, l, n, *nn[4];
86 rvec dx, dx1, dx2, rj, rk, urk, urj;
87 real cost, cost2, *sgmol, *skmol, rmean, rmean2, r2, box2, *r_nn[4];
88 t_pbc pbc;
89 t_mat *dmat;
90 t_dist *d;
91 int m1, mm, sl_index;
92 int **nnb, *sl_count;
93 real onethird = 1.0/3.0;
94 /* dmat = init_mat(maxidx, FALSE); */
95 box2 = box[XX0][XX0] * box[XX0][XX0];
96 snew(sl_count, nslice)(sl_count) = save_calloc("sl_count", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_order.c"
, 96, (nslice), sizeof(*(sl_count)));
97 for (i = 0; (i < 4); i++)
98 {
99 snew(r_nn[i], natoms)(r_nn[i]) = save_calloc("r_nn[i]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_order.c"
, 99, (natoms), sizeof(*(r_nn[i])));
100 snew(nn[i], natoms)(nn[i]) = save_calloc("nn[i]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_order.c"
, 100, (natoms), sizeof(*(nn[i])));
101
102 for (j = 0; (j < natoms); j++)
103 {
104 r_nn[i][j] = box2;
105 }
106 }
107
108 snew(sgmol, maxidx)(sgmol) = save_calloc("sgmol", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_order.c"
, 108, (maxidx), sizeof(*(sgmol)));
109 snew(skmol, maxidx)(skmol) = save_calloc("skmol", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_order.c"
, 109, (maxidx), sizeof(*(skmol)));
110
111 /* Must init pbc every step because of pressure coupling */
112 set_pbc(&pbc, ePBC, box);
113
114 gmx_rmpbc(gpbc, natoms, box, x);
115
116 nsgbin = 1 + 1/0.0005;
117 snew(sgbin, nsgbin)(sgbin) = save_calloc("sgbin", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_order.c"
, 117, (nsgbin), sizeof(*(sgbin)));
118
119 *sgmean = 0.0;
120 *skmean = 0.0;
121 l = 0;
122 for (i = 0; (i < maxidx); i++) /* loop over index file */
123 {
124 ix = index[i];
125 for (j = 0; (j < maxidx); j++)
126 {
127 if (i == j)
128 {
129 continue;
130 }
131
132 jx = index[j];
133
134 pbc_dx(&pbc, x[ix], x[jx], dx);
135 r2 = iprod(dx, dx);
136
137 /* set_mat_entry(dmat,i,j,r2); */
138
139 /* determine the nearest neighbours */
140 if (r2 < r_nn[0][i])
141 {
142 r_nn[3][i] = r_nn[2][i]; nn[3][i] = nn[2][i];
143 r_nn[2][i] = r_nn[1][i]; nn[2][i] = nn[1][i];
144 r_nn[1][i] = r_nn[0][i]; nn[1][i] = nn[0][i];
145 r_nn[0][i] = r2; nn[0][i] = j;
146 }
147 else if (r2 < r_nn[1][i])
148 {
149 r_nn[3][i] = r_nn[2][i]; nn[3][i] = nn[2][i];
150 r_nn[2][i] = r_nn[1][i]; nn[2][i] = nn[1][i];
151 r_nn[1][i] = r2; nn[1][i] = j;
152 }
153 else if (r2 < r_nn[2][i])
154 {
155 r_nn[3][i] = r_nn[2][i]; nn[3][i] = nn[2][i];
156 r_nn[2][i] = r2; nn[2][i] = j;
157 }
158 else if (r2 < r_nn[3][i])
159 {
160 r_nn[3][i] = r2; nn[3][i] = j;
161 }
162 }
163
164
165 /* calculate mean distance between nearest neighbours */
166 rmean = 0;
167 for (j = 0; (j < 4); j++)
168 {
169 r_nn[j][i] = sqrt(r_nn[j][i]);
170 rmean += r_nn[j][i];
171 }
172 rmean /= 4;
173
174 n = 0;
175 sgmol[i] = 0.0;
176 skmol[i] = 0.0;
177
178 /* Chau1998a eqn 3 */
179 /* angular part tetrahedrality order parameter per atom */
180 for (j = 0; (j < 3); j++)
181 {
182 for (k = j+1; (k < 4); k++)
183 {
184 pbc_dx(&pbc, x[ix], x[index[nn[k][i]]], rk);
185 pbc_dx(&pbc, x[ix], x[index[nn[j][i]]], rj);
186
187 unitv(rk, urk);
188 unitv(rj, urj);
189
190 cost = iprod(urk, urj) + onethird;
191 cost2 = cost * cost;
192
193 /* sgmol[i] += 3*cost2/32; */
194 sgmol[i] += cost2;
195
196 /* determine distribution */
197 ibin = nsgbin * cost2;
198 if (ibin < nsgbin)
199 {
200 sgbin[ibin]++;
201 }
202 /* printf("%d %d %f %d %d\n", j, k, cost * cost, ibin, sgbin[ibin]);*/
203 l++;
204 n++;
205 }
206 }
207
208 /* normalize sgmol between 0.0 and 1.0 */
209 sgmol[i] = 3*sgmol[i]/32;
210 *sgmean += sgmol[i];
211
212 /* distance part tetrahedrality order parameter per atom */
213 rmean2 = 4 * 3 * rmean * rmean;
214 for (j = 0; (j < 4); j++)
215 {
216 skmol[i] += (rmean - r_nn[j][i]) * (rmean - r_nn[j][i]) / rmean2;
217 /* printf("%d %f (%f %f %f %f) \n",
218 i, skmol[i], rmean, rmean2, r_nn[j][i], (rmean - r_nn[j][i]) );
219 */
220 }
221
222 *skmean += skmol[i];
223
224 /* Compute sliced stuff */
225 sl_index = gmx_nint((1+x[i][slice_dim]/box[slice_dim][slice_dim])*nslice) % nslice;
226 sgslice[sl_index] += sgmol[i];
227 skslice[sl_index] += skmol[i];
228 sl_count[sl_index]++;
229 } /* loop over entries in index file */
230
231 *sgmean /= maxidx;
232 *skmean /= maxidx;
233
234 for (i = 0; (i < nslice); i++)
235 {
236 if (sl_count[i] > 0)
237 {
238 sgslice[i] /= sl_count[i];
239 skslice[i] /= sl_count[i];
240 }
241 }
242 sfree(sl_count)save_free("sl_count", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_order.c"
, 242, (sl_count));
243 sfree(sgbin)save_free("sgbin", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_order.c"
, 243, (sgbin));
244 sfree(sgmol)save_free("sgmol", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_order.c"
, 244, (sgmol));
245 sfree(skmol)save_free("skmol", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_order.c"
, 245, (skmol));
246 for (i = 0; (i < 4); i++)
247 {
248 sfree(r_nn[i])save_free("r_nn[i]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_order.c"
, 248, (r_nn[i]));
249 sfree(nn[i])save_free("nn[i]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_order.c"
, 249, (nn[i]));
250 }
251}
252
253
254static void calc_tetra_order_parm(const char *fnNDX, const char *fnTPS,
255 const char *fnTRX, const char *sgfn,
256 const char *skfn,
257 int nslice, int slice_dim,
258 const char *sgslfn, const char *skslfn,
259 const output_env_t oenv)
260{
261 FILE *fpsg = NULL((void*)0), *fpsk = NULL((void*)0);
262 t_topology top;
263 int ePBC;
264 char title[STRLEN4096], fn[STRLEN4096], subtitle[STRLEN4096];
265 t_trxstatus *status;
266 int natoms;
267 real t;
268 rvec *xtop, *x;
269 matrix box;
270 real sg, sk;
271 atom_id **index;
272 char **grpname;
273 int i, *isize, ng, nframes;
274 real *sg_slice, *sg_slice_tot, *sk_slice, *sk_slice_tot;
275 gmx_rmpbc_t gpbc = NULL((void*)0);
276
277
278 read_tps_conf(fnTPS, title, &top, &ePBC, &xtop, NULL((void*)0), box, FALSE0);
279
280 snew(sg_slice, nslice)(sg_slice) = save_calloc("sg_slice", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_order.c"
, 280, (nslice), sizeof(*(sg_slice)));
281 snew(sk_slice, nslice)(sk_slice) = save_calloc("sk_slice", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_order.c"
, 281, (nslice), sizeof(*(sk_slice)));
282 snew(sg_slice_tot, nslice)(sg_slice_tot) = save_calloc("sg_slice_tot", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_order.c"
, 282, (nslice), sizeof(*(sg_slice_tot)));
283 snew(sk_slice_tot, nslice)(sk_slice_tot) = save_calloc("sk_slice_tot", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_order.c"
, 283, (nslice), sizeof(*(sk_slice_tot)));
284 ng = 1;
285 /* get index groups */
286 printf("Select the group that contains the atoms you want to use for the tetrahedrality order parameter calculation:\n");
287 snew(grpname, ng)(grpname) = save_calloc("grpname", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_order.c"
, 287, (ng), sizeof(*(grpname)));
288 snew(index, ng)(index) = save_calloc("index", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_order.c"
, 288, (ng), sizeof(*(index)));
289 snew(isize, ng)(isize) = save_calloc("isize", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_order.c"
, 289, (ng), sizeof(*(isize)));
290 get_index(&top.atoms, fnNDX, ng, isize, index, grpname);
291
292 /* Analyze trajectory */
293 natoms = read_first_x(oenv, &status, fnTRX, &t, &x, box);
294 if (natoms > top.atoms.nr)
295 {
296 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_order.c"
, 296, "Topology (%d atoms) does not match trajectory (%d atoms)",
297 top.atoms.nr, natoms);
298 }
299 check_index(NULL((void*)0), ng, index[0], NULL((void*)0), natoms);
300
301 fpsg = xvgropen(sgfn, "S\\sg\\N Angle Order Parameter", "Time (ps)", "S\\sg\\N",
302 oenv);
303 fpsk = xvgropen(skfn, "S\\sk\\N Distance Order Parameter", "Time (ps)", "S\\sk\\N",
304 oenv);
305
306 /* loop over frames */
307 gpbc = gmx_rmpbc_init(&top.idef, ePBC, natoms);
308 nframes = 0;
309 do
310 {
311 find_nearest_neighbours(ePBC, natoms, box, x, isize[0], index[0],
312 &sg, &sk, nslice, slice_dim, sg_slice, sk_slice, gpbc);
313 for (i = 0; (i < nslice); i++)
314 {
315 sg_slice_tot[i] += sg_slice[i];
316 sk_slice_tot[i] += sk_slice[i];
317 }
318 fprintf(fpsg, "%f %f\n", t, sg);
319 fprintf(fpsk, "%f %f\n", t, sk);
320 nframes++;
321 }
322 while (read_next_x(oenv, status, &t, x, box));
323 close_trj(status);
324 gmx_rmpbc_done(gpbc);
325
326 sfree(grpname)save_free("grpname", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_order.c"
, 326, (grpname));
327 sfree(index)save_free("index", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_order.c"
, 327, (index));
328 sfree(isize)save_free("isize", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_order.c"
, 328, (isize));
329
330 gmx_ffclose(fpsg);
331 gmx_ffclose(fpsk);
332
333 fpsg = xvgropen(sgslfn,
334 "S\\sg\\N Angle Order Parameter / Slab", "(nm)", "S\\sg\\N",
335 oenv);
336 fpsk = xvgropen(skslfn,
337 "S\\sk\\N Distance Order Parameter / Slab", "(nm)", "S\\sk\\N",
338 oenv);
339 for (i = 0; (i < nslice); i++)
340 {
341 fprintf(fpsg, "%10g %10g\n", (i+0.5)*box[slice_dim][slice_dim]/nslice,
342 sg_slice_tot[i]/nframes);
343 fprintf(fpsk, "%10g %10g\n", (i+0.5)*box[slice_dim][slice_dim]/nslice,
344 sk_slice_tot[i]/nframes);
345 }
346 gmx_ffclose(fpsg);
347 gmx_ffclose(fpsk);
348}
349
350
351/* Print name of first atom in all groups in index file */
352static void print_types(atom_id index[], atom_id a[], int ngrps,
353 char *groups[], t_topology *top)
354{
355 int i;
356
357 fprintf(stderrstderr, "Using following groups: \n");
358 for (i = 0; i < ngrps; i++)
359 {
360 fprintf(stderrstderr, "Groupname: %s First atomname: %s First atomnr %u\n",
361 groups[i], *(top->atoms.atomname[a[index[i]]]), a[index[i]]);
362 }
363 fprintf(stderrstderr, "\n");
364}
365
366static void check_length(real length, int a, int b)
367{
368 if (length > 0.3)
369 {
370 fprintf(stderrstderr, "WARNING: distance between atoms %d and "
371 "%d > 0.3 nm (%f). Index file might be corrupt.\n",
372 a, b, length);
373 }
374}
375
376void calc_order(const char *fn, atom_id *index, atom_id *a, rvec **order,
377 real ***slOrder, real *slWidth, int nslices, gmx_bool bSliced,
378 gmx_bool bUnsat, t_topology *top, int ePBC, int ngrps, int axis,
379 gmx_bool permolecule, gmx_bool radial, gmx_bool distcalc, const char *radfn,
380 real ***distvals,
381 const output_env_t oenv)
382{
383 /* if permolecule = TRUE, order parameters will be calculed per molecule
384 * and stored in slOrder with #slices = # molecules */
385 rvec *x0, /* coordinates with pbc */
386 *x1, /* coordinates without pbc */
387 dist; /* vector between two atoms */
388 matrix box; /* box (3x3) */
389 t_trxstatus *status;
390 rvec cossum, /* sum of vector angles for three axes */
391 Sx, Sy, Sz, /* the three molecular axes */
392 tmp1, tmp2, /* temp. rvecs for calculating dot products */
393 frameorder; /* order parameters for one frame */
394 real *slFrameorder; /* order parameter for one frame, per slice */
395 real length, /* total distance between two atoms */
396 t, /* time from trajectory */
397 z_ave, z1, z2; /* average z, used to det. which slice atom is in */
398 int natoms, /* nr. atoms in trj */
399 nr_tails, /* nr tails, to check if index file is correct */
400 size = 0, /* nr. of atoms in group. same as nr_tails */
401 i, j, m, k, l, teller = 0,
402 slice, /* current slice number */
403 nr_frames = 0;
404 int *slCount; /* nr. of atoms in one slice */
405 real dbangle = 0, /* angle between double bond and axis */
406 sdbangle = 0; /* sum of these angles */
407 gmx_bool use_unitvector = FALSE0; /* use a specified unit vector instead of axis to specify unit normal*/
408 rvec direction, com, dref, dvec;
409 int comsize, distsize;
410 atom_id *comidx = NULL((void*)0), *distidx = NULL((void*)0);
411 char *grpname = NULL((void*)0);
412 t_pbc pbc;
413 real arcdist, tmpdist;
414 gmx_rmpbc_t gpbc = NULL((void*)0);
415
416 /* PBC added for center-of-mass vector*/
417 /* Initiate the pbc structure */
418 memset(&pbc, 0, sizeof(pbc));
419
420 if ((natoms = read_first_x(oenv, &status, fn, &t, &x0, box)) == 0)
421 {
422 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_order.c"
, 422, "Could not read coordinates from statusfile\n");
423 }
424
425 nr_tails = index[1] - index[0];
426 fprintf(stderrstderr, "Number of elements in first group: %d\n", nr_tails);
427 /* take first group as standard. Not rocksolid, but might catch error in index*/
428
429 if (permolecule)
430 {
431 nslices = nr_tails;
432 bSliced = FALSE0; /*force slices off */
433 fprintf(stderrstderr, "Calculating order parameters for each of %d molecules\n",
434 nslices);
435 }
436
437 if (radial)
438 {
439 use_unitvector = TRUE1;
440 fprintf(stderrstderr, "Select an index group to calculate the radial membrane normal\n");
441 get_index(&top->atoms, radfn, 1, &comsize, &comidx, &grpname);
442 }
443 if (distcalc)
444 {
445 if (grpname != NULL((void*)0))
446 {
447 sfree(grpname)save_free("grpname", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_order.c"
, 447, (grpname));
448 }
449 fprintf(stderrstderr, "Select an index group to use as distance reference\n");
450 get_index(&top->atoms, radfn, 1, &distsize, &distidx, &grpname);
451 bSliced = FALSE0; /*force slices off*/
452 }
453
454 if (use_unitvector && bSliced)
455 {
456 fprintf(stderrstderr, "Warning: slicing and specified unit vectors are not currently compatible\n");
457 }
458
459 snew(slCount, nslices)(slCount) = save_calloc("slCount", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_order.c"
, 459, (nslices), sizeof(*(slCount)));
460 snew(*slOrder, nslices)(*slOrder) = save_calloc("*slOrder", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_order.c"
, 460, (nslices), sizeof(*(*slOrder)));
461 for (i = 0; i < nslices; i++)
462 {
463 snew((*slOrder)[i], ngrps)((*slOrder)[i]) = save_calloc("(*slOrder)[i]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_order.c"
, 463, (ngrps), sizeof(*((*slOrder)[i])));
464 }
465 if (distcalc)
466 {
467 snew(*distvals, nslices)(*distvals) = save_calloc("*distvals", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_order.c"
, 467, (nslices), sizeof(*(*distvals)));
468 for (i = 0; i < nslices; i++)
469 {
470 snew((*distvals)[i], ngrps)((*distvals)[i]) = save_calloc("(*distvals)[i]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_order.c"
, 470, (ngrps), sizeof(*((*distvals)[i])));
471 }
472 }
473 snew(*order, ngrps)(*order) = save_calloc("*order", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_order.c"
, 473, (ngrps), sizeof(*(*order)));
474 snew(slFrameorder, nslices)(slFrameorder) = save_calloc("slFrameorder", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_order.c"
, 474, (nslices), sizeof(*(slFrameorder)));
475 snew(x1, natoms)(x1) = save_calloc("x1", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_order.c"
, 475, (natoms), sizeof(*(x1)));
476
477 if (bSliced)
478 {
479 *slWidth = box[axis][axis]/nslices;
480 fprintf(stderrstderr, "Box divided in %d slices. Initial width of slice: %f\n",
481 nslices, *slWidth);
482 }
483
484
485#if 0
486 nr_tails = index[1] - index[0];
487 fprintf(stderrstderr, "Number of elements in first group: %d\n", nr_tails);
488 /* take first group as standard. Not rocksolid, but might catch error
489 in index*/
490#endif
491
492 teller = 0;
493
494 gpbc = gmx_rmpbc_init(&top->idef, ePBC, natoms);
495 /*********** Start processing trajectory ***********/
496 do
497 {
498 if (bSliced)
499 {
500 *slWidth = box[axis][axis]/nslices;
501 }
502 teller++;
503
504 set_pbc(&pbc, ePBC, box);
505 gmx_rmpbc_copy(gpbc, natoms, box, x0, x1);
506
507 /* Now loop over all groups. There are ngrps groups, the order parameter can
508 be calculated for grp 1 to grp ngrps - 1. For each group, loop over all
509 atoms in group, which is index[i] to (index[i+1] - 1) See block.h. Of
510 course, in this case index[i+1] -index[i] has to be the same for all
511 groups, namely the number of tails. i just runs over all atoms in a tail,
512 so for DPPC ngrps = 16 and i runs from 1 to 14, including 14
513 */
514
515
516 if (radial)
517 {
518 /*center-of-mass determination*/
519 com[XX0] = 0.0; com[YY1] = 0.0; com[ZZ2] = 0.0;
520 for (j = 0; j < comsize; j++)
521 {
522 rvec_inc(com, x1[comidx[j]]);
523 }
524 svmul(1.0/comsize, com, com);
525 }
526 if (distcalc)
527 {
528 dref[XX0] = 0.0; dref[YY1] = 0.0; dref[ZZ2] = 0.0;
529 for (j = 0; j < distsize; j++)
530 {
531 rvec_inc(dist, x1[distidx[j]]);
532 }
533 svmul(1.0/distsize, dref, dref);
534 if (radial)
535 {
536 pbc_dx(&pbc, dref, com, dvec);
537 unitv(dvec, dvec);
538 }
539 }
540
541 for (i = 1; i < ngrps - 1; i++)
542 {
543 clear_rvec(frameorder);
544
545 size = index[i+1] - index[i];
546 if (size != nr_tails)
547 {
548 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_order.c"
, 548, "grp %d does not have same number of"
549 " elements as grp 1\n", i);
550 }
551
552 for (j = 0; j < size; j++)
553 {
554 if (radial)
555 /*create unit vector*/
556 {
557 pbc_dx(&pbc, x1[a[index[i]+j]], com, direction);
558 unitv(direction, direction);
559 /*DEBUG*/
560 /*if (j==0)
561 fprintf(stderr,"X %f %f %f\tcom %f %f %f\tdirection %f %f %f\n",x1[a[index[i]+j]][0],x1[a[index[i]+j]][1],x1[a[index[i]+j]][2],com[0],com[1],com[2],
562 direction[0],direction[1],direction[2]);*/
563 }
564
565 if (bUnsat)
566 {
567 /* Using convention for unsaturated carbons */
568 /* first get Sz, the vector from Cn to Cn+1 */
569 rvec_sub(x1[a[index[i+1]+j]], x1[a[index[i]+j]], dist);
570 length = norm(dist);
571 check_length(length, a[index[i]+j], a[index[i+1]+j]);
572 svmul(1/length, dist, Sz);
573
574 /* this is actually the cosine of the angle between the double bond
575 and axis, because Sz is normalized and the two other components of
576 the axis on the bilayer are zero */
577 if (use_unitvector)
578 {
579 sdbangle += gmx_angle(direction, Sz); /*this can probably be optimized*/
580 }
581 else
582 {
583 sdbangle += acos(Sz[axis]);
584 }
585 }
586 else
587 {
588 /* get vector dist(Cn-1,Cn+1) for tail atoms */
589 rvec_sub(x1[a[index[i+1]+j]], x1[a[index[i-1]+j]], dist);
590 length = norm(dist); /* determine distance between two atoms */
591 check_length(length, a[index[i-1]+j], a[index[i+1]+j]);
592
593 svmul(1/length, dist, Sz);
594 /* Sz is now the molecular axis Sz, normalized and all that */
595 }
596
597 /* now get Sx. Sx is normal to the plane of Cn-1, Cn and Cn+1 so
598 we can use the outer product of Cn-1->Cn and Cn+1->Cn, I hope */
599 rvec_sub(x1[a[index[i+1]+j]], x1[a[index[i]+j]], tmp1);
600 rvec_sub(x1[a[index[i-1]+j]], x1[a[index[i]+j]], tmp2);
601 cprod(tmp1, tmp2, Sx);
602 svmul(1/norm(Sx), Sx, Sx);
603
604 /* now we can get Sy from the outer product of Sx and Sz */
605 cprod(Sz, Sx, Sy);
606 svmul(1/norm(Sy), Sy, Sy);
607
608 /* the square of cosine of the angle between dist and the axis.
609 Using the innerproduct, but two of the three elements are zero
610 Determine the sum of the orderparameter of all atoms in group
611 */
612 if (use_unitvector)
613 {
614 cossum[XX0] = sqr(iprod(Sx, direction)); /* this is allowed, since Sa is normalized */
615 cossum[YY1] = sqr(iprod(Sy, direction));
616 cossum[ZZ2] = sqr(iprod(Sz, direction));
617 }
618 else
619 {
620 cossum[XX0] = sqr(Sx[axis]); /* this is allowed, since Sa is normalized */
621 cossum[YY1] = sqr(Sy[axis]);
622 cossum[ZZ2] = sqr(Sz[axis]);
623 }
624
625 for (m = 0; m < DIM3; m++)
626 {
627 frameorder[m] += 0.5 * (3 * cossum[m] - 1);
628 }
629
630 if (bSliced)
631 {
632 /* get average coordinate in box length for slicing,
633 determine which slice atom is in, increase count for that
634 slice. slFrameorder and slOrder are reals, not
635 rvecs. Only the component [axis] of the order tensor is
636 kept, until I find it necessary to know the others too
637 */
638
639 z1 = x1[a[index[i-1]+j]][axis];
640 z2 = x1[a[index[i+1]+j]][axis];
641 z_ave = 0.5 * (z1 + z2);
642 if (z_ave < 0)
643 {
644 z_ave += box[axis][axis];
645 }
646 if (z_ave > box[axis][axis])
647 {
648 z_ave -= box[axis][axis];
649 }
650
651 slice = (int)(0.5 + (z_ave / (*slWidth))) - 1;
652 slCount[slice]++; /* determine slice, increase count */
653
654 slFrameorder[slice] += 0.5 * (3 * cossum[axis] - 1);
655 }
656 else if (permolecule)
657 {
658 /* store per-molecule order parameter
659 * To just track single-axis order: (*slOrder)[j][i] += 0.5 * (3 * iprod(cossum,direction) - 1);
660 * following is for Scd order: */
661 (*slOrder)[j][i] += -1* (0.3333 * (3 * cossum[XX0] - 1) + 0.3333 * 0.5 * (3 * cossum[YY1] - 1));
662 }
663 if (distcalc)
664 {
665 if (radial)
666 {
667 /* bin order parameter by arc distance from reference group*/
668 arcdist = gmx_angle(dvec, direction);
669 (*distvals)[j][i] += arcdist;
670 }
671 else if (i == 1)
672 {
673 /* Want minimum lateral distance to first group calculated */
674 tmpdist = trace(box); /* should be max value */
675 for (k = 0; k < distsize; k++)
676 {
677 pbc_dx(&pbc, x1[distidx[k]], x1[a[index[i]+j]], dvec);
678 /* at the moment, just remove dvec[axis] */
679 dvec[axis] = 0;
680 tmpdist = min(tmpdist, norm2(dvec))(((tmpdist) < (norm2(dvec))) ? (tmpdist) : (norm2(dvec)) );
681 }
682 //fprintf(stderr, "Min dist %f; trace %f\n", tmpdist, trace(box));
683 (*distvals)[j][i] += sqrt(tmpdist);
684 }
685 }
686 } /* end loop j, over all atoms in group */
687
688 for (m = 0; m < DIM3; m++)
689 {
690 (*order)[i][m] += (frameorder[m]/size);
691 }
692
693 if (!permolecule)
694 { /*Skip following if doing per-molecule*/
695 for (k = 0; k < nslices; k++)
696 {
697 if (slCount[k]) /* if no elements, nothing has to be added */
698 {
699 (*slOrder)[k][i] += slFrameorder[k]/slCount[k];
700 slFrameorder[k] = 0; slCount[k] = 0;
701 }
702 }
703 } /* end loop i, over all groups in indexfile */
704 }
705 nr_frames++;
706
707 }
708 while (read_next_x(oenv, status, &t, x0, box));
709 /*********** done with status file **********/
710
711 fprintf(stderrstderr, "\nRead trajectory. Printing parameters to file\n");
712 gmx_rmpbc_done(gpbc);
713
714 /* average over frames */
715 for (i = 1; i < ngrps - 1; i++)
716 {
717 svmul(1.0/nr_frames, (*order)[i], (*order)[i]);
718 fprintf(stderrstderr, "Atom %d Tensor: x=%g , y=%g, z=%g\n", i, (*order)[i][XX0],
719 (*order)[i][YY1], (*order)[i][ZZ2]);
720 if (bSliced || permolecule)
721 {
722 for (k = 0; k < nslices; k++)
723 {
724 (*slOrder)[k][i] /= nr_frames;
725 }
726 }
727 if (distcalc)
728 {
729 for (k = 0; k < nslices; k++)
730 {
731 (*distvals)[k][i] /= nr_frames;
732 }
733 }
734 }
735
736 if (bUnsat)
737 {
738 fprintf(stderrstderr, "Average angle between double bond and normal: %f\n",
739 180*sdbangle/(nr_frames * size*M_PI3.14159265358979323846));
740 }
741
742 sfree(x0)save_free("x0", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_order.c"
, 742, (x0)); /* free memory used by coordinate arrays */
743 sfree(x1)save_free("x1", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_order.c"
, 743, (x1));
744 if (comidx != NULL((void*)0))
745 {
746 sfree(comidx)save_free("comidx", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_order.c"
, 746, (comidx));
747 }
748 if (distidx != NULL((void*)0))
749 {
750 sfree(distidx)save_free("distidx", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_order.c"
, 750, (distidx));
751 }
752 if (grpname != NULL((void*)0))
753 {
754 sfree(grpname)save_free("grpname", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_order.c"
, 754, (grpname));
755 }
756}
757
758
759void order_plot(rvec order[], real *slOrder[], const char *afile, const char *bfile,
760 const char *cfile, int ngrps, int nslices, real slWidth, gmx_bool bSzonly,
761 gmx_bool permolecule, real **distvals, const output_env_t oenv)
762{
763 FILE *ord, *slOrd; /* xvgr files with order parameters */
764 int atom, slice; /* atom corresponding to order para.*/
765 char buf[256]; /* for xvgr title */
766 real S; /* order parameter averaged over all atoms */
767
768 if (permolecule)
769 {
770 sprintf(buf, "Scd order parameters");
771 ord = xvgropen(afile, buf, "Atom", "S", oenv);
772 sprintf(buf, "Orderparameters per atom per slice");
773 slOrd = xvgropen(bfile, buf, "Molecule", "S", oenv);
774 for (atom = 1; atom < ngrps - 1; atom++)
775 {
776 fprintf(ord, "%12d %12g\n", atom, -1 * (0.6667 * order[atom][XX0] +
777 0.333 * order[atom][YY1]));
778 }
779
780 for (slice = 0; slice < nslices; slice++)
781 {
782 fprintf(slOrd, "%12d\t", slice);
783 if (distvals)
784 {
785 fprintf(slOrd, "%12g\t", distvals[slice][1]); /*use distance value at second carbon*/
786 }
787 for (atom = 1; atom < ngrps - 1; atom++)
788 {
789 fprintf(slOrd, "%12g\t", slOrder[slice][atom]);
790 }
791 fprintf(slOrd, "\n");
792 }
793
794 }
795 else if (bSzonly)
796 {
797 sprintf(buf, "Orderparameters Sz per atom");
798 ord = xvgropen(afile, buf, "Atom", "S", oenv);
799 fprintf(stderrstderr, "ngrps = %d, nslices = %d", ngrps, nslices);
800
801 sprintf(buf, "Orderparameters per atom per slice");
802 slOrd = xvgropen(bfile, buf, "Slice", "S", oenv);
803
804 for (atom = 1; atom < ngrps - 1; atom++)
805 {
806 fprintf(ord, "%12d %12g\n", atom, order[atom][ZZ2]);
807 }
808
809 for (slice = 0; slice < nslices; slice++)
810 {
811 S = 0;
812 for (atom = 1; atom < ngrps - 1; atom++)
813 {
814 S += slOrder[slice][atom];
815 }
816 fprintf(slOrd, "%12g %12g\n", slice*slWidth, S/atom);
817 }
818
819 }
820 else
821 {
822 sprintf(buf, "Order tensor diagonal elements");
823 ord = xvgropen(afile, buf, "Atom", "S", oenv);
824 sprintf(buf, "Deuterium order parameters");
825 slOrd = xvgropen(cfile, buf, "Atom", "Scd", oenv);
826
827 for (atom = 1; atom < ngrps - 1; atom++)
828 {
829 fprintf(ord, "%12d %12g %12g %12g\n", atom, order[atom][XX0],
830 order[atom][YY1], order[atom][ZZ2]);
831 fprintf(slOrd, "%12d %12g\n", atom, -1 * (0.6667 * order[atom][XX0] +
832 0.333 * order[atom][YY1]));
833 }
834
835 gmx_ffclose(ord);
836 gmx_ffclose(slOrd);
837 }
838}
839
840void write_bfactors(t_filenm *fnm, int nfile, atom_id *index, atom_id *a, int nslices, int ngrps, real **order, t_topology *top, real **distvals, output_env_t oenv)
841{
842 /*function to write order parameters as B factors in PDB file using
843 first frame of trajectory*/
844 t_trxstatus *status;
845 int natoms;
846 t_trxframe fr, frout;
847 t_atoms useatoms;
848 int i, j, ctr, nout;
849
850 ngrps -= 2; /*we don't have an order parameter for the first or
851 last atom in each chain*/
852 nout = nslices*ngrps;
853 natoms = read_first_frame(oenv, &status, ftp2fn(efTRX, nfile, fnm), &fr,
Value stored to 'natoms' is never read
854 TRX_NEED_X(1<<1));
855 close_trj(status);
856 frout = fr;
857 frout.natoms = nout;
858 frout.bF = FALSE0;
859 frout.bV = FALSE0;
860 frout.x = 0;
861 snew(frout.x, nout)(frout.x) = save_calloc("frout.x", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_order.c"
, 861, (nout), sizeof(*(frout.x)));
862
863 init_t_atoms(&useatoms, nout, TRUE1);
864 useatoms.nr = nout;
865
866 /*initialize PDBinfo*/
867 for (i = 0; i < useatoms.nr; ++i)
868 {
869 useatoms.pdbinfo[i].type = 0;
870 useatoms.pdbinfo[i].occup = 0.0;
871 useatoms.pdbinfo[i].bfac = 0.0;
872 useatoms.pdbinfo[i].bAnisotropic = FALSE0;
873 }
874
875 for (j = 0, ctr = 0; j < nslices; j++)
876 {
877 for (i = 0; i < ngrps; i++, ctr++)
878 {
879 /*iterate along each chain*/
880 useatoms.pdbinfo[ctr].bfac = order[j][i+1];
881 if (distvals)
882 {
883 useatoms.pdbinfo[ctr].occup = distvals[j][i+1];
884 }
885 copy_rvec(fr.x[a[index[i+1]+j]], frout.x[ctr]);
886 useatoms.atomname[ctr] = top->atoms.atomname[a[index[i+1]+j]];
887 useatoms.atom[ctr] = top->atoms.atom[a[index[i+1]+j]];
888 useatoms.nres = max(useatoms.nres, useatoms.atom[ctr].resind+1)(((useatoms.nres) > (useatoms.atom[ctr].resind+1)) ? (useatoms
.nres) : (useatoms.atom[ctr].resind+1) );
889 useatoms.resinfo[useatoms.atom[ctr].resind] = top->atoms.resinfo[useatoms.atom[ctr].resind]; /*copy resinfo*/
890 }
891 }
892
893 write_sto_conf(opt2fn("-ob", nfile, fnm), "Order parameters", &useatoms, frout.x, NULL((void*)0), frout.ePBC, frout.box);
894
895 sfree(frout.x)save_free("frout.x", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_order.c"
, 895, (frout.x));
896 free_t_atoms(&useatoms, FALSE0);
897}
898
899int gmx_order(int argc, char *argv[])
900{
901 const char *desc[] = {
902 "[THISMODULE] computes the order parameter per atom for carbon tails. For atom i the",
903 "vector i-1, i+1 is used together with an axis. ",
904 "The index file should contain only the groups to be used for calculations,",
905 "with each group of equivalent carbons along the relevant acyl chain in its own",
906 "group. There should not be any generic groups (like System, Protein) in the index",
907 "file to avoid confusing the program (this is not relevant to tetrahedral order",
908 "parameters however, which only work for water anyway).[PAR]",
909 "[THISMODULE] can also give all",
910 "diagonal elements of the order tensor and even calculate the deuterium",
911 "order parameter Scd (default). If the option [TT]-szonly[tt] is given, only one",
912 "order tensor component (specified by the [TT]-d[tt] option) is given and the",
913 "order parameter per slice is calculated as well. If [TT]-szonly[tt] is not",
914 "selected, all diagonal elements and the deuterium order parameter is",
915 "given.[PAR]"
916 "The tetrahedrality order parameters can be determined",
917 "around an atom. Both angle an distance order parameters are calculated. See",
918 "P.-L. Chau and A.J. Hardwick, Mol. Phys., 93, (1998), 511-518.",
919 "for more details.[BR]",
920 ""
921 };
922
923 static int nslices = 1; /* nr of slices defined */
924 static gmx_bool bSzonly = FALSE0; /* True if only Sz is wanted */
925 static gmx_bool bUnsat = FALSE0; /* True if carbons are unsat. */
926 static const char *normal_axis[] = { NULL((void*)0), "z", "x", "y", NULL((void*)0) };
927 static gmx_bool permolecule = FALSE0; /*compute on a per-molecule basis */
928 static gmx_bool radial = FALSE0; /*compute a radial membrane normal */
929 static gmx_bool distcalc = FALSE0; /*calculate distance from a reference group */
930 t_pargs pa[] = {
931 { "-d", FALSE0, etENUM, {normal_axis},
932 "Direction of the normal on the membrane" },
933 { "-sl", FALSE0, etINT, {&nslices},
934 "Calculate order parameter as function of box length, dividing the box"
935 " into this number of slices." },
936 { "-szonly", FALSE0, etBOOL, {&bSzonly},
937 "Only give Sz element of order tensor. (axis can be specified with [TT]-d[tt])" },
938 { "-unsat", FALSE0, etBOOL, {&bUnsat},
939 "Calculate order parameters for unsaturated carbons. Note that this can"
940 "not be mixed with normal order parameters." },
941 { "-permolecule", FALSE0, etBOOL, {&permolecule},
942 "Compute per-molecule Scd order parameters" },
943 { "-radial", FALSE0, etBOOL, {&radial},
944 "Compute a radial membrane normal" },
945 { "-calcdist", FALSE0, etBOOL, {&distcalc},
946 "Compute distance from a reference" },
947 };
948
949 rvec *order; /* order par. for each atom */
950 real **slOrder; /* same, per slice */
951 real slWidth = 0.0; /* width of a slice */
952 char **grpname; /* groupnames */
953 int ngrps, /* nr. of groups */
954 i,
955 axis = 0; /* normal axis */
956 t_topology *top; /* topology */
957 int ePBC;
958 atom_id *index, /* indices for a */
959 *a; /* atom numbers in each group */
960 t_blocka *block; /* data from index file */
961 t_filenm fnm[] = { /* files for g_order */
962 { efTRX, "-f", NULL((void*)0), ffREAD1<<1 }, /* trajectory file */
963 { efNDX, "-n", NULL((void*)0), ffREAD1<<1 }, /* index file */
964 { efNDX, "-nr", NULL((void*)0), ffREAD1<<1 }, /* index for radial axis calculation */
965 { efTPX, NULL((void*)0), NULL((void*)0), ffREAD1<<1 }, /* topology file */
966 { efXVG, "-o", "order", ffWRITE1<<2 }, /* xvgr output file */
967 { efXVG, "-od", "deuter", ffWRITE1<<2 }, /* xvgr output file */
968 { efPDB, "-ob", NULL((void*)0), ffWRITE1<<2 }, /* write Scd as B factors to PDB if permolecule */
969 { efXVG, "-os", "sliced", ffWRITE1<<2 }, /* xvgr output file */
970 { efXVG, "-Sg", "sg-ang", ffOPTWR(1<<2| 1<<3) }, /* xvgr output file */
971 { efXVG, "-Sk", "sk-dist", ffOPTWR(1<<2| 1<<3) }, /* xvgr output file */
972 { efXVG, "-Sgsl", "sg-ang-slice", ffOPTWR(1<<2| 1<<3) }, /* xvgr output file */
973 { efXVG, "-Sksl", "sk-dist-slice", ffOPTWR(1<<2| 1<<3) }, /* xvgr output file */
974 };
975 gmx_bool bSliced = FALSE0; /* True if box is sliced */
976#define NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))) asize(fnm)((int)(sizeof(fnm)/sizeof((fnm)[0])))
977 real **distvals = NULL((void*)0);
978 const char *sgfnm, *skfnm, *ndxfnm, *tpsfnm, *trxfnm;
979 output_env_t oenv;
980
981 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),
982 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))
983 {
984 return 0;
985 }
986 if (nslices < 1)
987 {
988 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_order.c"
, 988, "Can not have nslices < 1");
989 }
990 sgfnm = opt2fn_null("-Sg", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm);
991 skfnm = opt2fn_null("-Sk", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm);
992 ndxfnm = opt2fn_null("-n", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm);
993 tpsfnm = ftp2fn(efTPX, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm);
994 trxfnm = ftp2fn(efTRX, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm);
995
996 /* Calculate axis */
997 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)
998 {
999 axis = XX0;
1000 }
1001 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)
1002 {
1003 axis = YY1;
1004 }
1005 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)
1006 {
1007 axis = ZZ2;
1008 }
1009 else
1010 {
1011 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_order.c"
, 1011, "Invalid axis, use x, y or z");
1012 }
1013
1014 switch (axis)
1015 {
1016 case 0:
1017 fprintf(stderrstderr, "Taking x axis as normal to the membrane\n");
1018 break;
1019 case 1:
1020 fprintf(stderrstderr, "Taking y axis as normal to the membrane\n");
1021 break;
1022 case 2:
1023 fprintf(stderrstderr, "Taking z axis as normal to the membrane\n");
1024 break;
1025 }
1026
1027 /* tetraheder order parameter */
1028 if (skfnm || sgfnm)
1029 {
1030 /* If either of theoptions is set we compute both */
1031 sgfnm = opt2fn("-Sg", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm);
1032 skfnm = opt2fn("-Sk", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm);
1033 calc_tetra_order_parm(ndxfnm, tpsfnm, trxfnm, sgfnm, skfnm, nslices, axis,
1034 opt2fn("-Sgsl", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), opt2fn("-Sksl", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm),
1035 oenv);
1036 /* view xvgr files */
1037 do_view(oenv, opt2fn("-Sg", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), NULL((void*)0));
1038 do_view(oenv, opt2fn("-Sk", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), NULL((void*)0));
1039 if (nslices > 1)
1040 {
1041 do_view(oenv, opt2fn("-Sgsl", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), NULL((void*)0));
1042 do_view(oenv, opt2fn("-Sksl", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), NULL((void*)0));
1043 }
1044 }
1045 else
1046 {
1047 /* tail order parameter */
1048
1049 if (nslices > 1)
1050 {
1051 bSliced = TRUE1;
1052 fprintf(stderrstderr, "Dividing box in %d slices.\n\n", nslices);
1053 }
1054
1055 if (bSzonly)
1056 {
1057 fprintf(stderrstderr, "Only calculating Sz\n");
1058 }
1059 if (bUnsat)
1060 {
1061 fprintf(stderrstderr, "Taking carbons as unsaturated!\n");
1062 }
1063
1064 top = read_top(ftp2fn(efTPX, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), &ePBC); /* read topology file */
1065
1066 block = init_index(ftp2fn(efNDX, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), &grpname);
1067 index = block->index; /* get indices from t_block block */
1068 a = block->a; /* see block.h */
1069 ngrps = block->nr;
1070
1071 if (permolecule)
1072 {
1073 nslices = index[1] - index[0]; /*I think this assumes contiguous lipids in topology*/
1074 fprintf(stderrstderr, "Calculating Scd order parameters for each of %d molecules\n", nslices);
1075 }
1076
1077 if (radial)
1078 {
1079 fprintf(stderrstderr, "Calculating radial distances\n");
1080 if (!permolecule)
1081 {
1082 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_order.c"
, 1082, "Cannot yet output radial distances without permolecule\n");
1083 }
1084 }
1085
1086 /* show atomtypes, to check if index file is correct */
1087 print_types(index, a, ngrps, grpname, top);
1088
1089 calc_order(ftp2fn(efTRX, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), index, a, &order,
1090 &slOrder, &slWidth, nslices, bSliced, bUnsat,
1091 top, ePBC, ngrps, axis, permolecule, radial, distcalc, opt2fn_null("-nr", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), &distvals, oenv);
1092
1093 if (radial)
1094 {
1095 ngrps--; /*don't print the last group--was used for
1096 center-of-mass determination*/
1097
1098 }
1099 order_plot(order, slOrder, opt2fn("-o", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), opt2fn("-os", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm),
1100 opt2fn("-od", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), ngrps, nslices, slWidth, bSzonly, permolecule, distvals, oenv);
1101
1102 if (opt2bSet("-ob", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm))
1103 {
1104 if (!permolecule)
1105 {
1106 fprintf(stderrstderr,
1107 "Won't write B-factors with averaged order parameters; use -permolecule\n");
1108 }
1109 else
1110 {
1111 write_bfactors(fnm, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), index, a, nslices, ngrps, slOrder, top, distvals, oenv);
1112 }
1113 }
1114
1115
1116 do_view(oenv, opt2fn("-o", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), NULL((void*)0)); /* view xvgr file */
1117 do_view(oenv, opt2fn("-os", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), NULL((void*)0)); /* view xvgr file */
1118 do_view(oenv, opt2fn("-od", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), NULL((void*)0)); /* view xvgr file */
1119 }
1120
1121 if (distvals != NULL((void*)0))
1122 {
1123 for (i = 0; i < nslices; ++i)
1124 {
1125 sfree(distvals[i])save_free("distvals[i]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_order.c"
, 1125, (distvals[i]));
1126 }
1127 sfree(distvals)save_free("distvals", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_order.c"
, 1127, (distvals));
1128 }
1129
1130 return 0;
1131}