File: | gromacs/gmxana/sfactor.c |
Location: | line 354, column 5 |
Description: | Value stored to 'positions' is never read |
1 | /* |
2 | * This file is part of the GROMACS molecular simulation package. |
3 | * |
4 | * Copyright (c) 1991-2000, University of Groningen, The Netherlands. |
5 | * Copyright (c) 2001-2004, The GROMACS development team. |
6 | * Copyright (c) 2013,2014, by the GROMACS development team, led by |
7 | * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl, |
8 | * and including many others, as listed in the AUTHORS file in the |
9 | * top-level source directory and at http://www.gromacs.org. |
10 | * |
11 | * GROMACS is free software; you can redistribute it and/or |
12 | * modify it under the terms of the GNU Lesser General Public License |
13 | * as published by the Free Software Foundation; either version 2.1 |
14 | * of the License, or (at your option) any later version. |
15 | * |
16 | * GROMACS is distributed in the hope that it will be useful, |
17 | * but WITHOUT ANY WARRANTY; without even the implied warranty of |
18 | * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU |
19 | * Lesser General Public License for more details. |
20 | * |
21 | * You should have received a copy of the GNU Lesser General Public |
22 | * License along with GROMACS; if not, see |
23 | * http://www.gnu.org/licenses, or write to the Free Software Foundation, |
24 | * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. |
25 | * |
26 | * If you want to redistribute modifications to GROMACS, please |
27 | * consider that scientific software is very special. Version |
28 | * control is crucial - bugs must be traceable. We will be happy to |
29 | * consider code for inclusion in the official distribution, but |
30 | * derived work must not be called official GROMACS. Details are found |
31 | * in the README & COPYING files - if they are missing, get the |
32 | * official version at http://www.gromacs.org. |
33 | * |
34 | * To help us fund GROMACS development, we humbly ask that you cite |
35 | * the research papers on the package. Check out http://www.gromacs.org. |
36 | */ |
37 | #ifdef HAVE_CONFIG_H1 |
38 | #include <config.h> |
39 | #endif |
40 | |
41 | #include "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 | |
60 | typedef 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 | |
72 | typedef struct reduced_atom{ |
73 | rvec x; |
74 | int t; |
75 | } reduced_atom; |
76 | |
77 | |
78 | typedef 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 | |
92 | extern 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 | |
132 | extern 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 | |
158 | extern 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 | |
257 | extern 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 | |
333 | extern 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; |
Value stored to 'positions' is never read | |
355 | } |
356 | |
357 | |
358 | extern 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 | |
421 | extern 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 | |
440 | extern 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 | |
554 | extern 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 | |
605 | extern 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); |
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 | |
651 | extern 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 | |
687 | extern 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 | |
711 | extern 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 | } |