Bug Summary

File:gromacs/gmxana/sfactor.c
Location:line 639, column 9
Description:Value stored to 'success' 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 "gromacs/utility/smalloc.h"
42#include "gromacs/utility/futil.h"
43#include "gromacs/math/utilities.h"
44#include "gromacs/utility/fatalerror.h"
45#include "gromacs/math/vec.h"
46#include "macros.h"
47#include "index.h"
48#include "gromacs/fileio/strdb.h"
49#include "gromacs/fileio/tpxio.h"
50#include "gromacs/fileio/trxio.h"
51#include "typedefs.h"
52#include "oenv.h"
53#include "gromacs/fileio/gmxfio.h"
54#include "gromacs/fileio/xvgr.h"
55#include "gromacs/fileio/matio.h"
56#include "names.h"
57#include "sfactor.h"
58
59
60typedef struct gmx_structurefactors {
61 int nratoms;
62 int *p; /* proton number */
63 int *n; /* neutron number */
64 /* Parameters for the Cromer Mann fit */
65 real **a; /* parameter a */
66 real **b; /* parameter b */
67 real *c; /* parameter c */
68 char **atomnm; /* atomname */
69
70} gmx_structurefactors;
71
72typedef struct reduced_atom{
73 rvec x;
74 int t;
75} reduced_atom;
76
77
78typedef struct structure_factor
79{
80 int n_angles;
81 int n_groups;
82 double lambda;
83 double energy;
84 double momentum;
85 double ref_k;
86 double **F;
87 int nSteps;
88 int total_n_atoms;
89} structure_factor;
90
91
92extern int * create_indexed_atom_type (reduced_atom_t * atm, int size)
93{
94/*
95 * create an index of the atom types found in a group
96 * i.e.: for water index_atp[0]=type_number_of_O and
97 * index_atp[1]=type_number_of_H
98 *
99 * the last element is set to 0
100 */
101 int *index_atp, i, i_tmp, j;
102
103 reduced_atom *att = (reduced_atom *)atm;
104
105 snew (index_atp, 1)(index_atp) = save_calloc("index_atp", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/sfactor.c"
, 105, (1), sizeof(*(index_atp)));
106 i_tmp = 1;
107 index_atp[0] = att[0].t;
108 for (i = 1; i < size; i++)
109 {
110 for (j = 0; j < i_tmp; j++)
111 {
112 if (att[i].t == index_atp[j])
113 {
114 break;
115 }
116 }
117 if (j == i_tmp) /* i.e. no indexed atom type is == to atm[i].t */
118 {
119 i_tmp++;
120 srenew (index_atp, i_tmp * sizeof (int))(index_atp) = save_realloc("index_atp", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/sfactor.c"
, 120, (index_atp), (i_tmp * sizeof (int)), sizeof(*(index_atp
)));
121 index_atp[i_tmp - 1] = att[i].t;
122 }
123 }
124 i_tmp++;
125 srenew (index_atp, i_tmp * sizeof (int))(index_atp) = save_realloc("index_atp", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/sfactor.c"
, 125, (index_atp), (i_tmp * sizeof (int)), sizeof(*(index_atp
)));
126 index_atp[i_tmp - 1] = 0;
127 return index_atp;
128}
129
130
131
132extern t_complex *** rc_tensor_allocation(int x, int y, int z)
133{
134 t_complex ***t;
135 int i, j;
136
137 snew(t, x)(t) = save_calloc("t", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/sfactor.c"
, 137, (x), sizeof(*(t)));
138 snew(t[0], x*y)(t[0]) = save_calloc("t[0]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/sfactor.c"
, 138, (x*y), sizeof(*(t[0])));
139 snew(t[0][0], x*y*z)(t[0][0]) = save_calloc("t[0][0]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/sfactor.c"
, 139, (x*y*z), sizeof(*(t[0][0])));
140
141 for (j = 1; j < y; j++)
142 {
143 t[0][j] = t[0][j-1] + z;
144 }
145 for (i = 1; i < x; i++)
146 {
147 t[i] = t[i-1] + y;
148 t[i][0] = t[i-1][0] + y*z;
149 for (j = 1; j < y; j++)
150 {
151 t[i][j] = t[i][j-1] + z;
152 }
153 }
154 return t;
155}
156
157
158extern void compute_structure_factor (structure_factor_t * sft, matrix box,
159 reduced_atom_t * red, int isize, real start_q,
160 real end_q, int group, real **sf_table)
161{
162 structure_factor *sf = (structure_factor *)sft;
163 reduced_atom *redt = (reduced_atom *)red;
164
165 t_complex ***tmpSF;
166 rvec k_factor;
167 real kdotx, asf, kx, ky, kz, krr;
168 int kr, maxkx, maxky, maxkz, i, j, k, p, *counter;
169
170
171 k_factor[XX0] = 2 * M_PI3.14159265358979323846 / box[XX0][XX0];
172 k_factor[YY1] = 2 * M_PI3.14159265358979323846 / box[YY1][YY1];
173 k_factor[ZZ2] = 2 * M_PI3.14159265358979323846 / box[ZZ2][ZZ2];
174
175 maxkx = (int) (end_q / k_factor[XX0] + 0.5);
176 maxky = (int) (end_q / k_factor[YY1] + 0.5);
177 maxkz = (int) (end_q / k_factor[ZZ2] + 0.5);
178
179 snew (counter, sf->n_angles)(counter) = save_calloc("counter", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/sfactor.c"
, 179, (sf->n_angles), sizeof(*(counter)));
180
181 tmpSF = rc_tensor_allocation(maxkx, maxky, maxkz);
182/*
183 * The big loop...
184 * compute real and imaginary part of the structure factor for every
185 * (kx,ky,kz))
186 */
187 fprintf(stderrstderr, "\n");
188 for (i = 0; i < maxkx; i++)
189 {
190 fprintf (stderrstderr, "\rdone %3.1f%% ", (double)(100.0*(i+1))/maxkx);
191 kx = i * k_factor[XX0];
192 for (j = 0; j < maxky; j++)
193 {
194 ky = j * k_factor[YY1];
195 for (k = 0; k < maxkz; k++)
196 {
197 if (i != 0 || j != 0 || k != 0)
198 {
199 kz = k * k_factor[ZZ2];
200 krr = sqrt (sqr (kx) + sqr (ky) + sqr (kz));
201 if (krr >= start_q && krr <= end_q)
202 {
203 kr = (int) (krr/sf->ref_k + 0.5);
204 if (kr < sf->n_angles)
205 {
206 counter[kr]++; /* will be used for the copmutation
207 of the average*/
208 for (p = 0; p < isize; p++)
209 {
210 asf = sf_table[redt[p].t][kr];
211
212 kdotx = kx * redt[p].x[XX0] +
213 ky * redt[p].x[YY1] + kz * redt[p].x[ZZ2];
214
215 tmpSF[i][j][k].re += cos (kdotx) * asf;
216 tmpSF[i][j][k].im += sin (kdotx) * asf;
217 }
218 }
219 }
220 }
221 }
222 }
223 } /* end loop on i */
224/*
225 * compute the square modulus of the structure factor, averaging on the surface
226 * kx*kx + ky*ky + kz*kz = krr*krr
227 * note that this is correct only for a (on the macroscopic scale)
228 * isotropic system.
229 */
230 for (i = 0; i < maxkx; i++)
231 {
232 kx = i * k_factor[XX0]; for (j = 0; j < maxky; j++)
233 {
234 ky = j * k_factor[YY1]; for (k = 0; k < maxkz; k++)
235 {
236 kz = k * k_factor[ZZ2]; krr = sqrt (sqr (kx) + sqr (ky)
237 + sqr (kz)); if (krr >= start_q && krr <= end_q)
238 {
239 kr = (int) (krr / sf->ref_k + 0.5);
240 if (kr < sf->n_angles && counter[kr] != 0)
241 {
242 sf->F[group][kr] +=
243 (sqr (tmpSF[i][j][k].re) +
244 sqr (tmpSF[i][j][k].im))/ counter[kr];
245 }
246 }
247 }
248 }
249 }
250 sfree(counter)save_free("counter", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/sfactor.c"
, 250, (counter));
251 sfree(tmpSF[0][0])save_free("tmpSF[0][0]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/sfactor.c"
, 251, (tmpSF[0][0]));
252 sfree(tmpSF[0])save_free("tmpSF[0]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/sfactor.c"
, 252, (tmpSF[0]));
253 sfree(tmpSF)save_free("tmpSF", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/sfactor.c"
, 253, (tmpSF));
254}
255
256
257extern gmx_structurefactors_t *gmx_structurefactors_init(const char *datfn)
258{
259
260 /* Read the database for the structure factor of the different atoms */
261
262 FILE *fp;
263 char line[STRLEN4096];
264 gmx_structurefactors *gsf;
265 double a1, a2, a3, a4, b1, b2, b3, b4, c;
266 int p;
267 int i;
268 int nralloc = 10;
269 int line_no;
270 char atomn[32];
271 fp = libopen(datfn);
272 line_no = 0;
273 snew(gsf, 1)(gsf) = save_calloc("gsf", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/sfactor.c"
, 273, (1), sizeof(*(gsf)));
274
275 snew(gsf->atomnm, nralloc)(gsf->atomnm) = save_calloc("gsf->atomnm", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/sfactor.c"
, 275, (nralloc), sizeof(*(gsf->atomnm)));
276 snew(gsf->a, nralloc)(gsf->a) = save_calloc("gsf->a", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/sfactor.c"
, 276, (nralloc), sizeof(*(gsf->a)));
277 snew(gsf->b, nralloc)(gsf->b) = save_calloc("gsf->b", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/sfactor.c"
, 277, (nralloc), sizeof(*(gsf->b)));
278 snew(gsf->c, nralloc)(gsf->c) = save_calloc("gsf->c", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/sfactor.c"
, 278, (nralloc), sizeof(*(gsf->c)));
279 snew(gsf->p, nralloc)(gsf->p) = save_calloc("gsf->p", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/sfactor.c"
, 279, (nralloc), sizeof(*(gsf->p)));
280 gsf->n = NULL((void*)0);
281 gsf->nratoms = line_no;
282 while (get_a_line(fp, line, STRLEN4096))
283 {
284 i = line_no;
285 if (sscanf(line, "%s %d %lf %lf %lf %lf %lf %lf %lf %lf %lf",
286 atomn, &p, &a1, &a2, &a3, &a4, &b1, &b2, &b3, &b4, &c) == 11)
287 {
288 gsf->atomnm[i] = strdup(atomn)(__extension__ (__builtin_constant_p (atomn) && ((size_t
)(const void *)((atomn) + 1) - (size_t)(const void *)(atomn) ==
1) ? (((const char *) (atomn))[0] == '\0' ? (char *) calloc (
(size_t) 1, (size_t) 1) : ({ size_t __len = strlen (atomn) + 1
; char *__retval = (char *) malloc (__len); if (__retval != (
(void*)0)) __retval = (char *) memcpy (__retval, atomn, __len
); __retval; })) : __strdup (atomn)));
289 gsf->p[i] = p;
290 snew(gsf->a[i], 4)(gsf->a[i]) = save_calloc("gsf->a[i]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/sfactor.c"
, 290, (4), sizeof(*(gsf->a[i])));
291 snew(gsf->b[i], 4)(gsf->b[i]) = save_calloc("gsf->b[i]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/sfactor.c"
, 291, (4), sizeof(*(gsf->b[i])));
292 gsf->a[i][0] = a1;
293 gsf->a[i][1] = a2;
294 gsf->a[i][2] = a3;
295 gsf->a[i][3] = a4;
296 gsf->b[i][0] = b1;
297 gsf->b[i][1] = b2;
298 gsf->b[i][2] = b3;
299 gsf->b[i][3] = b4;
300 gsf->c[i] = c;
301 line_no++;
302 gsf->nratoms = line_no;
303 if (line_no == nralloc)
304 {
305 nralloc += 10;
306 srenew(gsf->atomnm, nralloc)(gsf->atomnm) = save_realloc("gsf->atomnm", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/sfactor.c"
, 306, (gsf->atomnm), (nralloc), sizeof(*(gsf->atomnm))
);
307 srenew(gsf->a, nralloc)(gsf->a) = save_realloc("gsf->a", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/sfactor.c"
, 307, (gsf->a), (nralloc), sizeof(*(gsf->a)));
308 srenew(gsf->b, nralloc)(gsf->b) = save_realloc("gsf->b", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/sfactor.c"
, 308, (gsf->b), (nralloc), sizeof(*(gsf->b)));
309 srenew(gsf->c, nralloc)(gsf->c) = save_realloc("gsf->c", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/sfactor.c"
, 309, (gsf->c), (nralloc), sizeof(*(gsf->c)));
310 srenew(gsf->p, nralloc)(gsf->p) = save_realloc("gsf->p", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/sfactor.c"
, 310, (gsf->p), (nralloc), sizeof(*(gsf->p)));
311 }
312 }
313 else
314 {
315 fprintf(stderrstderr, "WARNING: Error in file %s at line %d ignored\n",
316 datfn, line_no);
317 }
318 }
319
320 srenew(gsf->atomnm, gsf->nratoms)(gsf->atomnm) = save_realloc("gsf->atomnm", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/sfactor.c"
, 320, (gsf->atomnm), (gsf->nratoms), sizeof(*(gsf->
atomnm)));
321 srenew(gsf->a, gsf->nratoms)(gsf->a) = save_realloc("gsf->a", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/sfactor.c"
, 321, (gsf->a), (gsf->nratoms), sizeof(*(gsf->a)));
322 srenew(gsf->b, gsf->nratoms)(gsf->b) = save_realloc("gsf->b", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/sfactor.c"
, 322, (gsf->b), (gsf->nratoms), sizeof(*(gsf->b)));
323 srenew(gsf->c, gsf->nratoms)(gsf->c) = save_realloc("gsf->c", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/sfactor.c"
, 323, (gsf->c), (gsf->nratoms), sizeof(*(gsf->c)));
324 srenew(gsf->p, gsf->nratoms)(gsf->p) = save_realloc("gsf->p", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/sfactor.c"
, 324, (gsf->p), (gsf->nratoms), sizeof(*(gsf->p)));
325
326 fclose(fp);
327
328 return (gmx_structurefactors_t *) gsf;
329
330}
331
332
333extern void rearrange_atoms (reduced_atom_t * positions, t_trxframe *fr, atom_id * index,
334 int isize, t_topology * top, gmx_bool flag, gmx_structurefactors_t *gsf)
335/* given the group's index, return the (continuous) array of atoms */
336{
337 int i;
338
339 reduced_atom *pos = (reduced_atom *)positions;
340
341 if (flag)
342 {
343 for (i = 0; i < isize; i++)
344 {
345 pos[i].t =
346 return_atom_type (*(top->atoms.atomname[index[i]]), gsf);
347 }
348 }
349 for (i = 0; i < isize; i++)
350 {
351 copy_rvec (fr->x[index[i]], pos[i].x);
352 }
353
354 positions = (reduced_atom_t *)pos;
355}
356
357
358extern int return_atom_type (const char *name, gmx_structurefactors_t *gsf)
359{
360 typedef struct {
361 const char *name;
362 int nh;
363 } t_united_h;
364 t_united_h uh[] = {
365 { "CH1", 1 }, { "CH2", 2 }, { "CH3", 3 },
366 { "CS1", 1 }, { "CS2", 2 }, { "CS3", 3 },
367 { "CP1", 1 }, { "CP2", 2 }, { "CP3", 3 }
368 };
369 int i, cnt = 0;
370 int *tndx;
371 int nrc;
372 int fndx = 0;
373 int NCMT;
374
375 gmx_structurefactors *gsft = (gmx_structurefactors *)gsf;
376
377 NCMT = gsft->nratoms;
378
379 snew(tndx, NCMT)(tndx) = save_calloc("tndx", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/sfactor.c"
, 379, (NCMT), sizeof(*(tndx)));
380
381 for (i = 0; (i < asize(uh)((int)(sizeof(uh)/sizeof((uh)[0])))); i++)
382 {
383 if (strcmp(name, uh[i].name)__extension__ ({ size_t __s1_len, __s2_len; (__builtin_constant_p
(name) && __builtin_constant_p (uh[i].name) &&
(__s1_len = strlen (name), __s2_len = strlen (uh[i].name), (
!((size_t)(const void *)((name) + 1) - (size_t)(const void *)
(name) == 1) || __s1_len >= 4) && (!((size_t)(const
void *)((uh[i].name) + 1) - (size_t)(const void *)(uh[i].name
) == 1) || __s2_len >= 4)) ? __builtin_strcmp (name, uh[i]
.name) : (__builtin_constant_p (name) && ((size_t)(const
void *)((name) + 1) - (size_t)(const void *)(name) == 1) &&
(__s1_len = strlen (name), __s1_len < 4) ? (__builtin_constant_p
(uh[i].name) && ((size_t)(const void *)((uh[i].name)
+ 1) - (size_t)(const void *)(uh[i].name) == 1) ? __builtin_strcmp
(name, uh[i].name) : (__extension__ ({ const unsigned char *
__s2 = (const unsigned char *) (const char *) (uh[i].name); int
__result = (((const unsigned char *) (const char *) (name))[
0] - __s2[0]); if (__s1_len > 0 && __result == 0) {
__result = (((const unsigned char *) (const char *) (name))[
1] - __s2[1]); if (__s1_len > 1 && __result == 0) {
__result = (((const unsigned char *) (const char *) (name))[
2] - __s2[2]); if (__s1_len > 2 && __result == 0) __result
= (((const unsigned char *) (const char *) (name))[3] - __s2
[3]); } } __result; }))) : (__builtin_constant_p (uh[i].name)
&& ((size_t)(const void *)((uh[i].name) + 1) - (size_t
)(const void *)(uh[i].name) == 1) && (__s2_len = strlen
(uh[i].name), __s2_len < 4) ? (__builtin_constant_p (name
) && ((size_t)(const void *)((name) + 1) - (size_t)(const
void *)(name) == 1) ? __builtin_strcmp (name, uh[i].name) : (
- (__extension__ ({ const unsigned char *__s2 = (const unsigned
char *) (const char *) (name); int __result = (((const unsigned
char *) (const char *) (uh[i].name))[0] - __s2[0]); if (__s2_len
> 0 && __result == 0) { __result = (((const unsigned
char *) (const char *) (uh[i].name))[1] - __s2[1]); if (__s2_len
> 1 && __result == 0) { __result = (((const unsigned
char *) (const char *) (uh[i].name))[2] - __s2[2]); if (__s2_len
> 2 && __result == 0) __result = (((const unsigned
char *) (const char *) (uh[i].name))[3] - __s2[3]); } } __result
; })))) : __builtin_strcmp (name, uh[i].name)))); }) == 0)
384 {
385 return NCMT-1+uh[i].nh;
386 }
387 }
388
389 for (i = 0; (i < NCMT); i++)
390 {
391 if (strncmp (name, gsft->atomnm[i], strlen(gsft->atomnm[i]))(__extension__ (__builtin_constant_p (strlen(gsft->atomnm[
i])) && ((__builtin_constant_p (name) && strlen
(name) < ((size_t) (strlen(gsft->atomnm[i])))) || (__builtin_constant_p
(gsft->atomnm[i]) && strlen (gsft->atomnm[i]) <
((size_t) (strlen(gsft->atomnm[i]))))) ? __extension__ ({
size_t __s1_len, __s2_len; (__builtin_constant_p (name) &&
__builtin_constant_p (gsft->atomnm[i]) && (__s1_len
= strlen (name), __s2_len = strlen (gsft->atomnm[i]), (!(
(size_t)(const void *)((name) + 1) - (size_t)(const void *)(name
) == 1) || __s1_len >= 4) && (!((size_t)(const void
*)((gsft->atomnm[i]) + 1) - (size_t)(const void *)(gsft->
atomnm[i]) == 1) || __s2_len >= 4)) ? __builtin_strcmp (name
, gsft->atomnm[i]) : (__builtin_constant_p (name) &&
((size_t)(const void *)((name) + 1) - (size_t)(const void *)
(name) == 1) && (__s1_len = strlen (name), __s1_len <
4) ? (__builtin_constant_p (gsft->atomnm[i]) && (
(size_t)(const void *)((gsft->atomnm[i]) + 1) - (size_t)(const
void *)(gsft->atomnm[i]) == 1) ? __builtin_strcmp (name, gsft
->atomnm[i]) : (__extension__ ({ const unsigned char *__s2
= (const unsigned char *) (const char *) (gsft->atomnm[i]
); int __result = (((const unsigned char *) (const char *) (name
))[0] - __s2[0]); if (__s1_len > 0 && __result == 0
) { __result = (((const unsigned char *) (const char *) (name
))[1] - __s2[1]); if (__s1_len > 1 && __result == 0
) { __result = (((const unsigned char *) (const char *) (name
))[2] - __s2[2]); if (__s1_len > 2 && __result == 0
) __result = (((const unsigned char *) (const char *) (name))
[3] - __s2[3]); } } __result; }))) : (__builtin_constant_p (gsft
->atomnm[i]) && ((size_t)(const void *)((gsft->
atomnm[i]) + 1) - (size_t)(const void *)(gsft->atomnm[i]) ==
1) && (__s2_len = strlen (gsft->atomnm[i]), __s2_len
< 4) ? (__builtin_constant_p (name) && ((size_t)(
const void *)((name) + 1) - (size_t)(const void *)(name) == 1
) ? __builtin_strcmp (name, gsft->atomnm[i]) : (- (__extension__
({ const unsigned char *__s2 = (const unsigned char *) (const
char *) (name); int __result = (((const unsigned char *) (const
char *) (gsft->atomnm[i]))[0] - __s2[0]); if (__s2_len >
0 && __result == 0) { __result = (((const unsigned char
*) (const char *) (gsft->atomnm[i]))[1] - __s2[1]); if (__s2_len
> 1 && __result == 0) { __result = (((const unsigned
char *) (const char *) (gsft->atomnm[i]))[2] - __s2[2]); if
(__s2_len > 2 && __result == 0) __result = (((const
unsigned char *) (const char *) (gsft->atomnm[i]))[3] - __s2
[3]); } } __result; })))) : __builtin_strcmp (name, gsft->
atomnm[i])))); }) : strncmp (name, gsft->atomnm[i], strlen
(gsft->atomnm[i])))) == 0)
392 {
393 tndx[cnt] = i;
394 cnt++;
395 }
396 }
397
398 if (cnt == 0)
399 {
400 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/sfactor.c"
, 400, "\nError: atom (%s) not in list (%d types checked)!\n",
401 name, i);
402 }
403 else
404 {
405 nrc = 0;
406 for (i = 0; i < cnt; i++)
407 {
408 if (strlen(gsft->atomnm[tndx[i]]) > (size_t)nrc)
409 {
410 nrc = strlen(gsft->atomnm[tndx[i]]);
411 fndx = tndx[i];
412 }
413 }
414
415 return fndx;
416 }
417
418 return 0;
419}
420
421extern int gmx_structurefactors_get_sf(gmx_structurefactors_t *gsf, int elem, real a[4], real b[4], real *c)
422{
423
424 int success;
425 int i;
426 gmx_structurefactors *gsft = (gmx_structurefactors *)gsf;
427 success = 0;
428
429 for (i = 0; i < 4; i++)
430 {
431 a[i] = gsft->a[elem][i];
432 b[i] = gsft->b[elem][i];
433 *c = gsft->c[elem];
434 }
435
436 success += 1;
437 return success;
438}
439
440extern int do_scattering_intensity (const char* fnTPS, const char* fnNDX,
441 const char* fnXVG, const char *fnTRX,
442 const char* fnDAT,
443 real start_q, real end_q,
444 real energy, int ng, const output_env_t oenv)
445{
446 int i, *isize, flags = TRX_READ_X(1<<0), **index_atp;
447 t_trxstatus *status;
448 char **grpname, title[STRLEN4096];
449 atom_id **index;
450 t_topology top;
451 int ePBC;
452 t_trxframe fr;
453 reduced_atom_t **red;
454 structure_factor *sf;
455 rvec *xtop;
456 real **sf_table;
457 int nsftable;
458 matrix box;
459 double r_tmp;
460
461 gmx_structurefactors_t *gmx_sf;
462 real *a, *b, c;
463 int success;
464
465 snew(a, 4)(a) = save_calloc("a", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/sfactor.c"
, 465, (4), sizeof(*(a)));
466 snew(b, 4)(b) = save_calloc("b", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/sfactor.c"
, 466, (4), sizeof(*(b)));
467
468
469 gmx_sf = gmx_structurefactors_init(fnDAT);
470
471 success = gmx_structurefactors_get_sf(gmx_sf, 0, a, b, &c);
472
473 snew (sf, 1)(sf) = save_calloc("sf", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/sfactor.c"
, 473, (1), sizeof(*(sf)));
474 sf->energy = energy;
475
476 /* Read the topology informations */
477 read_tps_conf (fnTPS, title, &top, &ePBC, &xtop, NULL((void*)0), box, TRUE1);
478 sfree (xtop)save_free("xtop", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/sfactor.c"
, 478, (xtop));
479
480 /* groups stuff... */
481 snew (isize, ng)(isize) = save_calloc("isize", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/sfactor.c"
, 481, (ng), sizeof(*(isize)));
482 snew (index, ng)(index) = save_calloc("index", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/sfactor.c"
, 482, (ng), sizeof(*(index)));
483 snew (grpname, ng)(grpname) = save_calloc("grpname", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/sfactor.c"
, 483, (ng), sizeof(*(grpname)));
484
485 fprintf (stderrstderr, "\nSelect %d group%s\n", ng,
486 ng == 1 ? "" : "s");
487 if (fnTPS)
488 {
489 get_index (&top.atoms, fnNDX, ng, isize, index, grpname);
490 }
491 else
492 {
493 rd_index (fnNDX, ng, isize, index, grpname);
494 }
495
496 /* The first time we read data is a little special */
497 read_first_frame (oenv, &status, fnTRX, &fr, flags);
498
499 sf->total_n_atoms = fr.natoms;
500
501 snew (red, ng)(red) = save_calloc("red", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/sfactor.c"
, 501, (ng), sizeof(*(red)));
502 snew (index_atp, ng)(index_atp) = save_calloc("index_atp", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/sfactor.c"
, 502, (ng), sizeof(*(index_atp)));
503
504 r_tmp = max (box[XX][XX], box[YY][YY])(((box[0][0]) > (box[1][1])) ? (box[0][0]) : (box[1][1]) );
505 r_tmp = (double) max (box[ZZ][ZZ], r_tmp)(((box[2][2]) > (r_tmp)) ? (box[2][2]) : (r_tmp) );
506
507 sf->ref_k = (2.0 * M_PI3.14159265358979323846) / (r_tmp);
508 /* ref_k will be the reference momentum unit */
509 sf->n_angles = (int) (end_q / sf->ref_k + 0.5);
510
511 snew (sf->F, ng)(sf->F) = save_calloc("sf->F", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/sfactor.c"
, 511, (ng), sizeof(*(sf->F)));
512 for (i = 0; i < ng; i++)
513 {
514 snew (sf->F[i], sf->n_angles)(sf->F[i]) = save_calloc("sf->F[i]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/sfactor.c"
, 514, (sf->n_angles), sizeof(*(sf->F[i])));
515 }
516 for (i = 0; i < ng; i++)
517 {
518 snew (red[i], isize[i])(red[i]) = save_calloc("red[i]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/sfactor.c"
, 518, (isize[i]), sizeof(*(red[i])));
519 rearrange_atoms (red[i], &fr, index[i], isize[i], &top, TRUE1, gmx_sf);
520 index_atp[i] = create_indexed_atom_type (red[i], isize[i]);
521 }
522
523 sf_table = compute_scattering_factor_table (gmx_sf, (structure_factor_t *)sf);
524
525
526 /* This is the main loop over frames */
527
528 do
529 {
530 sf->nSteps++;
531 for (i = 0; i < ng; i++)
532 {
533 rearrange_atoms (red[i], &fr, index[i], isize[i], &top, FALSE0, gmx_sf);
534
535 compute_structure_factor ((structure_factor_t *)sf, box, red[i], isize[i],
536 start_q, end_q, i, sf_table);
537 }
538 }
539
540 while (read_next_frame (oenv, status, &fr));
541
542 save_data ((structure_factor_t *)sf, fnXVG, ng, start_q, end_q, oenv);
543
544
545 sfree(a)save_free("a", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/sfactor.c"
, 545, (a));
546 sfree(b)save_free("b", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/sfactor.c"
, 546, (b));
547
548 gmx_structurefactors_done(gmx_sf);
549
550 return 0;
551}
552
553
554extern void save_data (structure_factor_t *sft, const char *file, int ngrps,
555 real start_q, real end_q, const output_env_t oenv)
556{
557
558 FILE *fp;
559 int i, g = 0;
560 double *tmp, polarization_factor, A;
561
562 structure_factor *sf = (structure_factor *)sft;
563
564 fp = xvgropen (file, "Scattering Intensity", "q (1/nm)",
565 "Intensity (a.u.)", oenv);
566
567 snew (tmp, ngrps)(tmp) = save_calloc("tmp", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/sfactor.c"
, 567, (ngrps), sizeof(*(tmp)));
568
569 for (g = 0; g < ngrps; g++)
570 {
571 for (i = 0; i < sf->n_angles; i++)
572 {
573
574/*
575 * theta is half the angle between incoming and scattered vectors.
576 *
577 * polar. fact. = 0.5*(1+cos^2(2*theta)) = 1 - 0.5 * sin^2(2*theta)
578 *
579 * sin(theta) = q/(2k) := A -> sin^2(theta) = 4*A^2 (1-A^2) ->
580 * -> 0.5*(1+cos^2(2*theta)) = 1 - 2 A^2 (1-A^2)
581 */
582 A = (double) (i * sf->ref_k) / (2.0 * sf->momentum);
583 polarization_factor = 1 - 2.0 * sqr (A) * (1 - sqr (A));
584 sf->F[g][i] *= polarization_factor;
585 }
586 }
587 for (i = 0; i < sf->n_angles; i++)
588 {
589 if (i * sf->ref_k >= start_q && i * sf->ref_k <= end_q)
590 {
591 fprintf (fp, "%10.5f ", i * sf->ref_k);
592 for (g = 0; g < ngrps; g++)
593 {
594 fprintf (fp, " %10.5f ", (sf->F[g][i]) /( sf->total_n_atoms*
595 sf->nSteps));
596 }
597 fprintf (fp, "\n");
598 }
599 }
600
601 gmx_ffclose (fp);
602}
603
604
605extern double CMSF (gmx_structurefactors_t *gsf, int type, int nh, double lambda, double sin_theta)
606/*
607 * return Cromer-Mann fit for the atomic scattering factor:
608 * sin_theta is the sine of half the angle between incoming and scattered
609 * vectors. See g_sq.h for a short description of CM fit.
610 */
611{
612 int i, success;
613 double tmp = 0.0, k2;
614 real *a, *b;
615 real c;
616
617 snew(a, 4)(a) = save_calloc("a", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/sfactor.c"
, 617, (4), sizeof(*(a)));
618 snew(b, 4)(b) = save_calloc("b", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/sfactor.c"
, 618, (4), sizeof(*(b)));
619
620 /*
621 *
622 * f0[k] = c + [SUM a_i*EXP(-b_i*(k^2)) ]
623 * i=1,4
624 */
625
626 /*
627 * united atoms case
628 * CH2 / CH3 groups
629 */
630 if (nh > 0)
631 {
632 tmp = (CMSF (gsf, return_atom_type ("C", gsf), 0, lambda, sin_theta) +
633 nh*CMSF (gsf, return_atom_type ("H", gsf), 0, lambda, sin_theta));
634 }
635 /* all atom case */
636 else
637 {
638 k2 = (sqr (sin_theta) / sqr (10.0 * lambda));
639 success = gmx_structurefactors_get_sf(gsf, type, a, b, &c);
Value stored to 'success' is never read
640 tmp = c;
641 for (i = 0; (i < 4); i++)
642 {
643 tmp += a[i] * exp (-b[i] * k2);
644 }
645 }
646 return tmp;
647}
648
649
650
651extern real **gmx_structurefactors_table(gmx_structurefactors_t *gsf, real momentum, real ref_k, real lambda, int n_angles)
652{
653
654 int NCMT;
655 int nsftable;
656 int i, j;
657 double q, sin_theta;
658 real **sf_table;
659 gmx_structurefactors *gsft = (gmx_structurefactors *)gsf;
660
661 NCMT = gsft->nratoms;
662 nsftable = NCMT+3;
663
664 snew (sf_table, nsftable)(sf_table) = save_calloc("sf_table", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/sfactor.c"
, 664, (nsftable), sizeof(*(sf_table)));
665 for (i = 0; (i < nsftable); i++)
666 {
667 snew (sf_table[i], n_angles)(sf_table[i]) = save_calloc("sf_table[i]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/sfactor.c"
, 667, (n_angles), sizeof(*(sf_table[i])));
668 for (j = 0; j < n_angles; j++)
669 {
670 q = ((double) j * ref_k);
671 /* theta is half the angle between incoming
672 and scattered wavevectors. */
673 sin_theta = q / (2.0 * momentum);
674 if (i < NCMT)
675 {
676 sf_table[i][j] = CMSF (gsf, i, 0, lambda, sin_theta);
677 }
678 else
679 {
680 sf_table[i][j] = CMSF (gsf, i, i-NCMT+1, lambda, sin_theta);
681 }
682 }
683 }
684 return sf_table;
685}
686
687extern void gmx_structurefactors_done(gmx_structurefactors_t *gsf)
688{
689
690 int i;
691 gmx_structurefactors *sf;
692 sf = (gmx_structurefactors *) gsf;
693
694 for (i = 0; i < sf->nratoms; i++)
695 {
696 sfree(sf->a[i])save_free("sf->a[i]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/sfactor.c"
, 696, (sf->a[i]));
697 sfree(sf->b[i])save_free("sf->b[i]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/sfactor.c"
, 697, (sf->b[i]));
698 sfree(sf->atomnm[i])save_free("sf->atomnm[i]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/sfactor.c"
, 698, (sf->atomnm[i]));
699 }
700
701 sfree(sf->a)save_free("sf->a", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/sfactor.c"
, 701, (sf->a));
702 sfree(sf->b)save_free("sf->b", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/sfactor.c"
, 702, (sf->b));
703 sfree(sf->atomnm)save_free("sf->atomnm", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/sfactor.c"
, 703, (sf->atomnm));
704 sfree(sf->p)save_free("sf->p", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/sfactor.c"
, 704, (sf->p));
705 sfree(sf->c)save_free("sf->c", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/sfactor.c"
, 705, (sf->c));
706
707 sfree(sf)save_free("sf", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/sfactor.c"
, 707, (sf));
708
709}
710
711extern real **compute_scattering_factor_table (gmx_structurefactors_t *gsf, structure_factor_t *sft)
712{
713/*
714 * this function build up a table of scattering factors for every atom
715 * type and for every scattering angle.
716 */
717
718 double hc = 1239.842;
719 real ** sf_table;
720
721 structure_factor *sf = (structure_factor *)sft;
722
723
724 /* \hbar \omega \lambda = hc = 1239.842 eV * nm */
725 sf->momentum = ((double) (2. * 1000.0 * M_PI3.14159265358979323846 * sf->energy) / hc);
726 sf->lambda = hc / (1000.0 * sf->energy);
727 fprintf (stderrstderr, "\nwavelenght = %f nm\n", sf->lambda);
728
729 sf_table = gmx_structurefactors_table(gsf, sf->momentum, sf->ref_k, sf->lambda, sf->n_angles);
730
731 return sf_table;
732}