Bug Summary

File:gromacs/gmxana/gmx_do_dssp.c
Location:line 714, column 29
Description:Dereference of null pointer

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) 2012,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 <stdlib.h>
42
43#include "typedefs.h"
44#include "gromacs/utility/cstringutil.h"
45#include "gromacs/fileio/strdb.h"
46#include "macros.h"
47#include "gromacs/utility/smalloc.h"
48#include "mshift.h"
49#include "gromacs/commandline/pargs.h"
50#include "gromacs/fileio/pdbio.h"
51#include "gromacs/utility/fatalerror.h"
52#include "gromacs/fileio/xvgr.h"
53#include "gromacs/fileio/matio.h"
54#include "index.h"
55#include "gstat.h"
56#include "gromacs/fileio/tpxio.h"
57#include "gromacs/fileio/trxio.h"
58#include "viewit.h"
59
60static int strip_dssp(char *dsspfile, int nres,
61 gmx_bool bPhobres[], real t,
62 real *acc, FILE *fTArea,
63 t_matrix *mat, int average_area[],
64 const output_env_t oenv)
65{
66 static gmx_bool bFirst = TRUE1;
67 static char *ssbuf;
68 FILE *tapeout;
69 static int xsize, frame;
70 char buf[STRLEN4096+1];
71 char SSTP;
72 int i, nr, iacc, nresidues;
73 int naccf, naccb; /* Count hydrophobic and hydrophilic residues */
74 real iaccf, iaccb;
75 t_xpmelmt c;
76
77 tapeout = gmx_ffopen(dsspfile, "r");
78
79 /* Skip header */
80 do
81 {
82 fgets2(buf, STRLEN4096, tapeout);
83 }
84 while (strstr(buf, "KAPPA") == NULL((void*)0));
85 if (bFirst)
86 {
87 /* Since we also have empty lines in the dssp output (temp) file,
88 * and some line content is saved to the ssbuf variable,
89 * we need more memory than just nres elements. To be shure,
90 * we allocate 2*nres-1, since for each chain there is a
91 * separating line in the temp file. (At most each residue
92 * could have been defined as a separate chain.) */
93 snew(ssbuf, 2*nres-1)(ssbuf) = save_calloc("ssbuf", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_do_dssp.c"
, 93, (2*nres-1), sizeof(*(ssbuf)));
94 }
95
96 iaccb = iaccf = 0;
97 nresidues = 0;
98 naccf = 0;
99 naccb = 0;
100 for (nr = 0; (fgets2(buf, STRLEN4096, tapeout) != NULL((void*)0)); nr++)
101 {
102 if (buf[13] == '!') /* Chain separator line has '!' at pos. 13 */
103 {
104 SSTP = '='; /* Chain separator sign '=' */
105 }
106 else
107 {
108 SSTP = buf[16] == ' ' ? '~' : buf[16];
109 }
110 ssbuf[nr] = SSTP;
111
112 buf[39] = '\0';
113
114 /* Only calculate solvent accessible area if needed */
115 if ((NULL((void*)0) != acc) && (buf[13] != '!'))
116 {
117 sscanf(&(buf[34]), "%d", &iacc);
118 acc[nr] = iacc;
119 /* average_area and bPhobres are counted from 0...nres-1 */
120 average_area[nresidues] += iacc;
121 if (bPhobres[nresidues])
122 {
123 naccb++;
124 iaccb += iacc;
125 }
126 else
127 {
128 naccf++;
129 iaccf += iacc;
130 }
131 /* Keep track of the residue number (do not count chain separator lines '!') */
132 nresidues++;
133 }
134 }
135 ssbuf[nr] = '\0';
136
137 if (bFirst)
138 {
139 if (0 != acc)
140 {
141 fprintf(stderrstderr, "%d residues were classified as hydrophobic and %d as hydrophilic.\n", naccb, naccf);
142 }
143
144 sprintf(mat->title, "Secondary structure");
145 mat->legend[0] = 0;
146 sprintf(mat->label_x, "%s", output_env_get_time_label(oenv));
147 sprintf(mat->label_y, "Residue");
148 mat->bDiscrete = TRUE1;
149 mat->ny = nr;
150 snew(mat->axis_y, nr)(mat->axis_y) = save_calloc("mat->axis_y", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_do_dssp.c"
, 150, (nr), sizeof(*(mat->axis_y)));
151 for (i = 0; i < nr; i++)
152 {
153 mat->axis_y[i] = i+1;
154 }
155 mat->axis_x = NULL((void*)0);
156 mat->matrix = NULL((void*)0);
157 xsize = 0;
158 frame = 0;
159 bFirst = FALSE0;
160 }
161 if (frame >= xsize)
162 {
163 xsize += 10;
164 srenew(mat->axis_x, xsize)(mat->axis_x) = save_realloc("mat->axis_x", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_do_dssp.c"
, 164, (mat->axis_x), (xsize), sizeof(*(mat->axis_x)));
165 srenew(mat->matrix, xsize)(mat->matrix) = save_realloc("mat->matrix", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_do_dssp.c"
, 165, (mat->matrix), (xsize), sizeof(*(mat->matrix)));
166 }
167 mat->axis_x[frame] = t;
168 snew(mat->matrix[frame], nr)(mat->matrix[frame]) = save_calloc("mat->matrix[frame]"
, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_do_dssp.c"
, 168, (nr), sizeof(*(mat->matrix[frame])));
169 c.c2 = 0;
170 for (i = 0; i < nr; i++)
171 {
172 c.c1 = ssbuf[i];
173 mat->matrix[frame][i] = max(0, searchcmap(mat->nmap, mat->map, c))(((0) > (searchcmap(mat->nmap, mat->map, c))) ? (0) :
(searchcmap(mat->nmap, mat->map, c)) );
174 }
175 frame++;
176 mat->nx = frame;
177
178 if (fTArea)
179 {
180 fprintf(fTArea, "%10g %10g %10g\n", t, 0.01*iaccb, 0.01*iaccf);
181 }
182 gmx_ffclose(tapeout);
183
184 /* Return the number of lines found in the dssp file (i.e. number
185 * of redidues plus chain separator lines).
186 * This is the number of y elements needed for the area xpm file */
187 return nr;
188}
189
190static gmx_bool *bPhobics(t_atoms *atoms)
191{
192 int i, nb;
193 char **cb;
194 gmx_bool *bb;
195
196
197 nb = get_strings("phbres.dat", &cb);
198 snew(bb, atoms->nres)(bb) = save_calloc("bb", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_do_dssp.c"
, 198, (atoms->nres), sizeof(*(bb)));
199
200 for (i = 0; (i < atoms->nres); i++)
201 {
202 if (-1 != search_str(nb, cb, *atoms->resinfo[i].name) )
203 {
204 bb[i] = TRUE1;
205 }
206 }
207 return bb;
208}
209
210static void check_oo(t_atoms *atoms)
211{
212 char *OOO;
213
214 int i;
215
216 OOO = strdup("O")(__extension__ (__builtin_constant_p ("O") && ((size_t
)(const void *)(("O") + 1) - (size_t)(const void *)("O") == 1
) ? (((const char *) ("O"))[0] == '\0' ? (char *) calloc ((size_t
) 1, (size_t) 1) : ({ size_t __len = strlen ("O") + 1; char *
__retval = (char *) malloc (__len); if (__retval != ((void*)0
)) __retval = (char *) memcpy (__retval, "O", __len); __retval
; })) : __strdup ("O")));
217
218 for (i = 0; (i < atoms->nr); i++)
219 {
220 if (strcmp(*(atoms->atomname[i]), "OXT")__extension__ ({ size_t __s1_len, __s2_len; (__builtin_constant_p
(*(atoms->atomname[i])) && __builtin_constant_p (
"OXT") && (__s1_len = strlen (*(atoms->atomname[i]
)), __s2_len = strlen ("OXT"), (!((size_t)(const void *)((*(atoms
->atomname[i])) + 1) - (size_t)(const void *)(*(atoms->
atomname[i])) == 1) || __s1_len >= 4) && (!((size_t
)(const void *)(("OXT") + 1) - (size_t)(const void *)("OXT") ==
1) || __s2_len >= 4)) ? __builtin_strcmp (*(atoms->atomname
[i]), "OXT") : (__builtin_constant_p (*(atoms->atomname[i]
)) && ((size_t)(const void *)((*(atoms->atomname[i
])) + 1) - (size_t)(const void *)(*(atoms->atomname[i])) ==
1) && (__s1_len = strlen (*(atoms->atomname[i])),
__s1_len < 4) ? (__builtin_constant_p ("OXT") && (
(size_t)(const void *)(("OXT") + 1) - (size_t)(const void *)(
"OXT") == 1) ? __builtin_strcmp (*(atoms->atomname[i]), "OXT"
) : (__extension__ ({ const unsigned char *__s2 = (const unsigned
char *) (const char *) ("OXT"); int __result = (((const unsigned
char *) (const char *) (*(atoms->atomname[i])))[0] - __s2
[0]); if (__s1_len > 0 && __result == 0) { __result
= (((const unsigned char *) (const char *) (*(atoms->atomname
[i])))[1] - __s2[1]); if (__s1_len > 1 && __result
== 0) { __result = (((const unsigned char *) (const char *) (
*(atoms->atomname[i])))[2] - __s2[2]); if (__s1_len > 2
&& __result == 0) __result = (((const unsigned char *
) (const char *) (*(atoms->atomname[i])))[3] - __s2[3]); }
} __result; }))) : (__builtin_constant_p ("OXT") && (
(size_t)(const void *)(("OXT") + 1) - (size_t)(const void *)(
"OXT") == 1) && (__s2_len = strlen ("OXT"), __s2_len <
4) ? (__builtin_constant_p (*(atoms->atomname[i])) &&
((size_t)(const void *)((*(atoms->atomname[i])) + 1) - (size_t
)(const void *)(*(atoms->atomname[i])) == 1) ? __builtin_strcmp
(*(atoms->atomname[i]), "OXT") : (- (__extension__ ({ const
unsigned char *__s2 = (const unsigned char *) (const char *)
(*(atoms->atomname[i])); int __result = (((const unsigned
char *) (const char *) ("OXT"))[0] - __s2[0]); if (__s2_len >
0 && __result == 0) { __result = (((const unsigned char
*) (const char *) ("OXT"))[1] - __s2[1]); if (__s2_len > 1
&& __result == 0) { __result = (((const unsigned char
*) (const char *) ("OXT"))[2] - __s2[2]); if (__s2_len > 2
&& __result == 0) __result = (((const unsigned char *
) (const char *) ("OXT"))[3] - __s2[3]); } } __result; })))) :
__builtin_strcmp (*(atoms->atomname[i]), "OXT")))); }) == 0)
221 {
222 *atoms->atomname[i] = OOO;
223 }
224 else if (strcmp(*(atoms->atomname[i]), "O1")__extension__ ({ size_t __s1_len, __s2_len; (__builtin_constant_p
(*(atoms->atomname[i])) && __builtin_constant_p (
"O1") && (__s1_len = strlen (*(atoms->atomname[i])
), __s2_len = strlen ("O1"), (!((size_t)(const void *)((*(atoms
->atomname[i])) + 1) - (size_t)(const void *)(*(atoms->
atomname[i])) == 1) || __s1_len >= 4) && (!((size_t
)(const void *)(("O1") + 1) - (size_t)(const void *)("O1") ==
1) || __s2_len >= 4)) ? __builtin_strcmp (*(atoms->atomname
[i]), "O1") : (__builtin_constant_p (*(atoms->atomname[i])
) && ((size_t)(const void *)((*(atoms->atomname[i]
)) + 1) - (size_t)(const void *)(*(atoms->atomname[i])) ==
1) && (__s1_len = strlen (*(atoms->atomname[i])),
__s1_len < 4) ? (__builtin_constant_p ("O1") && (
(size_t)(const void *)(("O1") + 1) - (size_t)(const void *)("O1"
) == 1) ? __builtin_strcmp (*(atoms->atomname[i]), "O1") :
(__extension__ ({ const unsigned char *__s2 = (const unsigned
char *) (const char *) ("O1"); int __result = (((const unsigned
char *) (const char *) (*(atoms->atomname[i])))[0] - __s2
[0]); if (__s1_len > 0 && __result == 0) { __result
= (((const unsigned char *) (const char *) (*(atoms->atomname
[i])))[1] - __s2[1]); if (__s1_len > 1 && __result
== 0) { __result = (((const unsigned char *) (const char *) (
*(atoms->atomname[i])))[2] - __s2[2]); if (__s1_len > 2
&& __result == 0) __result = (((const unsigned char *
) (const char *) (*(atoms->atomname[i])))[3] - __s2[3]); }
} __result; }))) : (__builtin_constant_p ("O1") && (
(size_t)(const void *)(("O1") + 1) - (size_t)(const void *)("O1"
) == 1) && (__s2_len = strlen ("O1"), __s2_len < 4
) ? (__builtin_constant_p (*(atoms->atomname[i])) &&
((size_t)(const void *)((*(atoms->atomname[i])) + 1) - (size_t
)(const void *)(*(atoms->atomname[i])) == 1) ? __builtin_strcmp
(*(atoms->atomname[i]), "O1") : (- (__extension__ ({ const
unsigned char *__s2 = (const unsigned char *) (const char *)
(*(atoms->atomname[i])); int __result = (((const unsigned
char *) (const char *) ("O1"))[0] - __s2[0]); if (__s2_len >
0 && __result == 0) { __result = (((const unsigned char
*) (const char *) ("O1"))[1] - __s2[1]); if (__s2_len > 1
&& __result == 0) { __result = (((const unsigned char
*) (const char *) ("O1"))[2] - __s2[2]); if (__s2_len > 2
&& __result == 0) __result = (((const unsigned char *
) (const char *) ("O1"))[3] - __s2[3]); } } __result; })))) :
__builtin_strcmp (*(atoms->atomname[i]), "O1")))); }) == 0)
225 {
226 *atoms->atomname[i] = OOO;
227 }
228 else if (strcmp(*(atoms->atomname[i]), "OC1")__extension__ ({ size_t __s1_len, __s2_len; (__builtin_constant_p
(*(atoms->atomname[i])) && __builtin_constant_p (
"OC1") && (__s1_len = strlen (*(atoms->atomname[i]
)), __s2_len = strlen ("OC1"), (!((size_t)(const void *)((*(atoms
->atomname[i])) + 1) - (size_t)(const void *)(*(atoms->
atomname[i])) == 1) || __s1_len >= 4) && (!((size_t
)(const void *)(("OC1") + 1) - (size_t)(const void *)("OC1") ==
1) || __s2_len >= 4)) ? __builtin_strcmp (*(atoms->atomname
[i]), "OC1") : (__builtin_constant_p (*(atoms->atomname[i]
)) && ((size_t)(const void *)((*(atoms->atomname[i
])) + 1) - (size_t)(const void *)(*(atoms->atomname[i])) ==
1) && (__s1_len = strlen (*(atoms->atomname[i])),
__s1_len < 4) ? (__builtin_constant_p ("OC1") && (
(size_t)(const void *)(("OC1") + 1) - (size_t)(const void *)(
"OC1") == 1) ? __builtin_strcmp (*(atoms->atomname[i]), "OC1"
) : (__extension__ ({ const unsigned char *__s2 = (const unsigned
char *) (const char *) ("OC1"); int __result = (((const unsigned
char *) (const char *) (*(atoms->atomname[i])))[0] - __s2
[0]); if (__s1_len > 0 && __result == 0) { __result
= (((const unsigned char *) (const char *) (*(atoms->atomname
[i])))[1] - __s2[1]); if (__s1_len > 1 && __result
== 0) { __result = (((const unsigned char *) (const char *) (
*(atoms->atomname[i])))[2] - __s2[2]); if (__s1_len > 2
&& __result == 0) __result = (((const unsigned char *
) (const char *) (*(atoms->atomname[i])))[3] - __s2[3]); }
} __result; }))) : (__builtin_constant_p ("OC1") && (
(size_t)(const void *)(("OC1") + 1) - (size_t)(const void *)(
"OC1") == 1) && (__s2_len = strlen ("OC1"), __s2_len <
4) ? (__builtin_constant_p (*(atoms->atomname[i])) &&
((size_t)(const void *)((*(atoms->atomname[i])) + 1) - (size_t
)(const void *)(*(atoms->atomname[i])) == 1) ? __builtin_strcmp
(*(atoms->atomname[i]), "OC1") : (- (__extension__ ({ const
unsigned char *__s2 = (const unsigned char *) (const char *)
(*(atoms->atomname[i])); int __result = (((const unsigned
char *) (const char *) ("OC1"))[0] - __s2[0]); if (__s2_len >
0 && __result == 0) { __result = (((const unsigned char
*) (const char *) ("OC1"))[1] - __s2[1]); if (__s2_len > 1
&& __result == 0) { __result = (((const unsigned char
*) (const char *) ("OC1"))[2] - __s2[2]); if (__s2_len > 2
&& __result == 0) __result = (((const unsigned char *
) (const char *) ("OC1"))[3] - __s2[3]); } } __result; })))) :
__builtin_strcmp (*(atoms->atomname[i]), "OC1")))); }) == 0)
229 {
230 *atoms->atomname[i] = OOO;
231 }
232 }
233}
234
235static void norm_acc(t_atoms *atoms, int nres,
236 real av_area[], real norm_av_area[])
237{
238 int i, n, n_surf;
239
240 char surffn[] = "surface.dat";
241 char **surf_res, **surf_lines;
242 double *surf;
243
244 n_surf = get_lines(surffn, &surf_lines);
245 snew(surf, n_surf)(surf) = save_calloc("surf", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_do_dssp.c"
, 245, (n_surf), sizeof(*(surf)));
246 snew(surf_res, n_surf)(surf_res) = save_calloc("surf_res", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_do_dssp.c"
, 246, (n_surf), sizeof(*(surf_res)));
247 for (i = 0; (i < n_surf); i++)
248 {
249 snew(surf_res[i], 5)(surf_res[i]) = save_calloc("surf_res[i]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_do_dssp.c"
, 249, (5), sizeof(*(surf_res[i])));
250 sscanf(surf_lines[i], "%s %lf", surf_res[i], &surf[i]);
251 }
252
253 for (i = 0; (i < nres); i++)
254 {
255 n = search_str(n_surf, surf_res, *atoms->resinfo[i].name);
256 if (n != -1)
257 {
258 norm_av_area[i] = av_area[i] / surf[n];
259 }
260 else
261 {
262 fprintf(stderrstderr, "Residue %s not found in surface database (%s)\n",
263 *atoms->resinfo[i].name, surffn);
264 }
265 }
266}
267
268void prune_ss_legend(t_matrix *mat)
269{
270 gmx_bool *present;
271 int *newnum;
272 int i, r, f, newnmap;
273 t_mapping *newmap;
274
275 snew(present, mat->nmap)(present) = save_calloc("present", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_do_dssp.c"
, 275, (mat->nmap), sizeof(*(present)));
276 snew(newnum, mat->nmap)(newnum) = save_calloc("newnum", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_do_dssp.c"
, 276, (mat->nmap), sizeof(*(newnum)));
277
278 for (f = 0; f < mat->nx; f++)
19
←
Loop condition is true. Entering loop body
→
22
←
Loop condition is false. Execution continues on line 286
→
279 {
280 for (r = 0; r < mat->ny; r++)
20
←
Loop condition is true. Entering loop body
→
21
←
Loop condition is false. Execution continues on line 278
→
281 {
282 present[mat->matrix[f][r]] = TRUE1;
283 }
284 }
285
286 newnmap = 0;
287 newmap = NULL((void*)0);
23
←
Null pointer value stored to 'newmap'
→
288 for (i = 0; i < mat->nmap; i++)
24
←
Loop condition is false. Execution continues on line 299
→
289 {
290 newnum[i] = -1;
291 if (present[i])
292 {
293 newnum[i] = newnmap;
294 newnmap++;
295 srenew(newmap, newnmap)(newmap) = save_realloc("newmap", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_do_dssp.c"
, 295, (newmap), (newnmap), sizeof(*(newmap)));
296 newmap[newnmap-1] = mat->map[i];
297 }
298 }
299 if (newnmap != mat->nmap)
25
←
Taking true branch
→
300 {
301 mat->nmap = newnmap;
302 mat->map = newmap;
26
←
Null pointer value stored to 'mat.map'
→
303 for (f = 0; f < mat->nx; f++)
27
←
Loop condition is true. Entering loop body
→
30
←
Loop condition is false. Execution continues on line 303
→
304 {
305 for (r = 0; r < mat->ny; r++)
28
←
Loop condition is true. Entering loop body
→
29
←
Loop condition is false. Execution continues on line 303
→
306 {
307 mat->matrix[f][r] = newnum[mat->matrix[f][r]];
308 }
309 }
310 }
311}
312
313void write_sas_mat(const char *fn, real **accr, int nframe, int nres, t_matrix *mat)
314{
315 real lo, hi;
316 int i, j, nlev;
317 t_rgb rlo = {1, 1, 1}, rhi = {0, 0, 0};
318 FILE *fp;
319
320 if (fn)
321 {
322 hi = lo = accr[0][0];
323 for (i = 0; i < nframe; i++)
324 {
325 for (j = 0; j < nres; j++)
326 {
327 lo = min(lo, accr[i][j])(((lo) < (accr[i][j])) ? (lo) : (accr[i][j]) );
328 hi = max(hi, accr[i][j])(((hi) > (accr[i][j])) ? (hi) : (accr[i][j]) );
329 }
330 }
331 fp = gmx_ffopen(fn, "w");
332 nlev = hi-lo+1;
333 write_xpm(fp, 0, "Solvent Accessible Surface", "Surface (A^2)",
334 "Time", "Residue Index", nframe, nres,
335 mat->axis_x, mat->axis_y, accr, lo, hi, rlo, rhi, &nlev);
336 gmx_ffclose(fp);
337 }
338}
339
340void analyse_ss(const char *outfile, t_matrix *mat, const char *ss_string,
341 const output_env_t oenv)
342{
343 FILE *fp;
344 t_mapping *map;
345 int f, r, *count, *total, ss_count, total_count;
346 size_t s;
347 const char** leg;
348
349 map = mat->map;
350 snew(count, mat->nmap)(count) = save_calloc("count", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_do_dssp.c"
, 350, (mat->nmap), sizeof(*(count)));
351 snew(total, mat->nmap)(total) = save_calloc("total", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_do_dssp.c"
, 351, (mat->nmap), sizeof(*(total)));
352 snew(leg, mat->nmap+1)(leg) = save_calloc("leg", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_do_dssp.c"
, 352, (mat->nmap+1), sizeof(*(leg)));
353 leg[0] = "Structure";
354 for (s = 0; s < (size_t)mat->nmap; s++)
355 {
356 leg[s+1] = strdup(map[s].desc)(__extension__ (__builtin_constant_p (map[s].desc) &&
((size_t)(const void *)((map[s].desc) + 1) - (size_t)(const void
*)(map[s].desc) == 1) ? (((const char *) (map[s].desc))[0] ==
'\0' ? (char *) calloc ((size_t) 1, (size_t) 1) : ({ size_t __len
= strlen (map[s].desc) + 1; char *__retval = (char *) malloc
(__len); if (__retval != ((void*)0)) __retval = (char *) memcpy
(__retval, map[s].desc, __len); __retval; })) : __strdup (map
[s].desc)));
357 }
358
359 fp = xvgropen(outfile, "Secondary Structure",
360 output_env_get_xvgr_tlabel(oenv), "Number of Residues", oenv);
361 if (output_env_get_print_xvgr_codes(oenv))
362 {
363 fprintf(fp, "@ subtitle \"Structure = ");
364 }
365 for (s = 0; s < strlen(ss_string); s++)
366 {
367 if (s > 0)
368 {
369 fprintf(fp, " + ");
370 }
371 for (f = 0; f < mat->nmap; f++)
372 {
373 if (ss_string[s] == map[f].code.c1)
374 {
375 fprintf(fp, "%s", map[f].desc);
376 }
377 }
378 }
379 fprintf(fp, "\"\n");
380 xvgr_legend(fp, mat->nmap+1, leg, oenv);
381
382 total_count = 0;
383 for (s = 0; s < (size_t)mat->nmap; s++)
384 {
385 total[s] = 0;
386 }
387 for (f = 0; f < mat->nx; f++)
388 {
389 ss_count = 0;
390 for (s = 0; s < (size_t)mat->nmap; s++)
391 {
392 count[s] = 0;
393 }
394 for (r = 0; r < mat->ny; r++)
395 {
396 count[mat->matrix[f][r]]++;
397 total[mat->matrix[f][r]]++;
398 }
399 for (s = 0; s < (size_t)mat->nmap; s++)
400 {
401 if (strchr(ss_string, map[s].code.c1)(__extension__ (__builtin_constant_p (map[s].code.c1) &&
!__builtin_constant_p (ss_string) && (map[s].code.c1
) == '\0' ? (char *) __rawmemchr (ss_string, map[s].code.c1) :
__builtin_strchr (ss_string, map[s].code.c1))))
402 {
403 ss_count += count[s];
404 total_count += count[s];
405 }
406 }
407 fprintf(fp, "%8g %5d", mat->axis_x[f], ss_count);
408 for (s = 0; s < (size_t)mat->nmap; s++)
409 {
410 fprintf(fp, " %5d", count[s]);
411 }
412 fprintf(fp, "\n");
413 }
414 /* now print column totals */
415 fprintf(fp, "%-8s %5d", "# Totals", total_count);
416 for (s = 0; s < (size_t)mat->nmap; s++)
417 {
418 fprintf(fp, " %5d", total[s]);
419 }
420 fprintf(fp, "\n");
421
422 /* now print percentages */
423 fprintf(fp, "%-8s %5.2f", "# SS %", total_count / (real) (mat->nx * mat->ny));
424 for (s = 0; s < (size_t)mat->nmap; s++)
425 {
426 fprintf(fp, " %5.2f", total[s] / (real) (mat->nx * mat->ny));
427 }
428 fprintf(fp, "\n");
429
430 gmx_ffclose(fp);
431 sfree(leg)save_free("leg", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_do_dssp.c"
, 431, (leg));
432 sfree(count)save_free("count", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_do_dssp.c"
, 432, (count));
433}
434
435int gmx_do_dssp(int argc, char *argv[])
436{
437 const char *desc[] = {
438 "[THISMODULE] ",
439 "reads a trajectory file and computes the secondary structure for",
440 "each time frame ",
441 "calling the dssp program. If you do not have the dssp program,",
442 "get it from http://swift.cmbi.ru.nl/gv/dssp. [THISMODULE] assumes ",
443 "that the dssp executable is located in ",
444 "[TT]/usr/local/bin/dssp[tt]. If this is not the case, then you should",
445 "set an environment variable [TT]DSSP[tt] pointing to the dssp",
446 "executable, e.g.: [PAR]",
447 "[TT]setenv DSSP /opt/dssp/bin/dssp[tt][PAR]",
448 "Since version 2.0.0, dssp is invoked with a syntax that differs",
449 "from earlier versions. If you have an older version of dssp,",
450 "use the [TT]-ver[tt] option to direct do_dssp to use the older syntax.",
451 "By default, do_dssp uses the syntax introduced with version 2.0.0.",
452 "Even newer versions (which at the time of writing are not yet released)",
453 "are assumed to have the same syntax as 2.0.0.[PAR]",
454 "The structure assignment for each residue and time is written to an",
455 "[TT].xpm[tt] matrix file. This file can be visualized with for instance",
456 "[TT]xv[tt] and can be converted to postscript with [TT]xpm2ps[tt].",
457 "Individual chains are separated by light grey lines in the [TT].xpm[tt] and",
458 "postscript files.",
459 "The number of residues with each secondary structure type and the",
460 "total secondary structure ([TT]-sss[tt]) count as a function of",
461 "time are also written to file ([TT]-sc[tt]).[PAR]",
462 "Solvent accessible surface (SAS) per residue can be calculated, both in",
463 "absolute values (A^2) and in fractions of the maximal accessible",
464 "surface of a residue. The maximal accessible surface is defined as",
465 "the accessible surface of a residue in a chain of glycines.",
466 "[BB]Note[bb] that the program [gmx-sas] can also compute SAS",
467 "and that is more efficient.[PAR]",
468 "Finally, this program can dump the secondary structure in a special file",
469 "[TT]ssdump.dat[tt] for usage in the program [gmx-chi]. Together",
470 "these two programs can be used to analyze dihedral properties as a",
471 "function of secondary structure type."
472 };
473 static gmx_bool bVerbose;
474 static const char *ss_string = "HEBT";
475 static int dsspVersion = 2;
476 t_pargs pa[] = {
477 { "-v", FALSE0, etBOOL, {&bVerbose},
478 "HIDDENGenerate miles of useless information" },
479 { "-sss", FALSE0, etSTR, {&ss_string},
480 "Secondary structures for structure count"},
481 { "-ver", FALSE0, etINT, {&dsspVersion},
482 "DSSP major version. Syntax changed with version 2"}
483 };
484
485 t_trxstatus *status;
486 FILE *tapein;
487 FILE *ss, *acc, *fTArea, *tmpf;
488 const char *fnSCount, *fnArea, *fnTArea, *fnAArea;
489 const char *leg[] = { "Phobic", "Phylic" };
490 t_topology top;
491 int ePBC;
492 t_atoms *atoms;
493 t_matrix mat;
494 int nres, nr0, naccr, nres_plus_separators;
495 gmx_bool *bPhbres, bDoAccSurf;
496 real t;
497 int i, j, natoms, nframe = 0;
498 matrix box;
499 int gnx;
500 char *grpnm, *ss_str;
501 atom_id *index;
502 rvec *xp, *x;
503 int *average_area;
504 real **accr, *accr_ptr = NULL((void*)0), *av_area, *norm_av_area;
505 char pdbfile[32], tmpfile[32], title[256];
506 char dssp[256];
507 const char *dptr;
508 output_env_t oenv;
509 gmx_rmpbc_t gpbc = NULL((void*)0);
510
511 t_filenm fnm[] = {
512 { efTRX, "-f", NULL((void*)0), ffREAD1<<1 },
513 { efTPS, NULL((void*)0), NULL((void*)0), ffREAD1<<1 },
514 { efNDX, NULL((void*)0), NULL((void*)0), ffOPTRD(1<<1 | 1<<3) },
515 { efDAT, "-ssdump", "ssdump", ffOPTWR(1<<2| 1<<3) },
516 { efMAP, "-map", "ss", ffLIBRD(1<<1 | 1<<4) },
517 { efXPM, "-o", "ss", ffWRITE1<<2 },
518 { efXVG, "-sc", "scount", ffWRITE1<<2 },
519 { efXPM, "-a", "area", ffOPTWR(1<<2| 1<<3) },
520 { efXVG, "-ta", "totarea", ffOPTWR(1<<2| 1<<3) },
521 { efXVG, "-aa", "averarea", ffOPTWR(1<<2| 1<<3) }
522 };
523#define NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))) asize(fnm)((int)(sizeof(fnm)/sizeof((fnm)[0])))
524
525 if (!parse_common_args(&argc, argv,
1
Taking false branch
→
526 PCA_CAN_TIME((1<<6) | (1<<7) | (1<<14)) | PCA_CAN_VIEW(1<<5) | PCA_TIME_UNIT(1<<15) | PCA_BE_NICE(1<<13),
527 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))
528 {
529 return 0;
530 }
531 fnSCount = opt2fn("-sc", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm);
532 fnArea = opt2fn_null("-a", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm);
533 fnTArea = opt2fn_null("-ta", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm);
534 fnAArea = opt2fn_null("-aa", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm);
535 bDoAccSurf = (fnArea || fnTArea || fnAArea);
536
537 read_tps_conf(ftp2fn(efTPS, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), title, &top, &ePBC, &xp, NULL((void*)0), box, FALSE0);
538 atoms = &(top.atoms);
539 check_oo(atoms);
540 bPhbres = bPhobics(atoms);
541
542 get_index(atoms, ftp2fn_null(efNDX, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), 1, &gnx, &index, &grpnm);
543 nres = 0;
544 nr0 = -1;
545 for (i = 0; (i < gnx); i++)
2
←
Assuming 'i' is >= 'gnx'
→
3
←
Loop condition is false. Execution continues on line 553
→
546 {
547 if (atoms->atom[index[i]].resind != nr0)
548 {
549 nr0 = atoms->atom[index[i]].resind;
550 nres++;
551 }
552 }
553 fprintf(stderrstderr, "There are %d residues in your selected group\n", nres);
554
555 strcpy(pdbfile, "ddXXXXXX");
556 gmx_tmpnam(pdbfile);
557 if ((tmpf = fopen(pdbfile, "w")) == NULL((void*)0))
4
←
Taking false branch
→
558 {
559 sprintf(pdbfile, "%ctmp%cfilterXXXXXX", DIR_SEPARATOR'/', DIR_SEPARATOR'/');
560 gmx_tmpnam(pdbfile);
561 if ((tmpf = fopen(pdbfile, "w")) == NULL((void*)0))
562 {
563 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_do_dssp.c"
, 563, "Can not open tmp file %s", pdbfile);
564 }
565 }
566 else
567 {
568 fclose(tmpf);
569 }
570
571 strcpy(tmpfile, "ddXXXXXX");
572 gmx_tmpnam(tmpfile);
573 if ((tmpf = fopen(tmpfile, "w")) == NULL((void*)0))
5
←
Taking false branch
→
574 {
575 sprintf(tmpfile, "%ctmp%cfilterXXXXXX", DIR_SEPARATOR'/', DIR_SEPARATOR'/');
576 gmx_tmpnam(tmpfile);
577 if ((tmpf = fopen(tmpfile, "w")) == NULL((void*)0))
578 {
579 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_do_dssp.c"
, 579, "Can not open tmp file %s", tmpfile);
580 }
581 }
582 else
583 {
584 fclose(tmpf);
585 }
586
587 if ((dptr = getenv("DSSP")) == NULL((void*)0))
6
←
Taking false branch
→
588 {
589 dptr = "/usr/local/bin/dssp";
590 }
591 if (!gmx_fexist(dptr))
7
←
Taking false branch
→
592 {
593 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_do_dssp.c"
, 593, "DSSP executable (%s) does not exist (use setenv DSSP)",
594 dptr);
595 }
596 if (dsspVersion >= 2)
8
←
Taking true branch
→
597 {
598 if (dsspVersion > 2)
9
←
Taking false branch
→
599 {
600 printf("\nWARNING: You use DSSP version %d, which is not explicitly\nsupported by do_dssp. Assuming version 2 syntax.\n\n", dsspVersion);
601 }
602
603 sprintf(dssp, "%s -i %s -o %s > /dev/null %s",
604 dptr, pdbfile, tmpfile, bVerbose ? "" : "2> /dev/null");
10
←
'?' condition is false
→
605 }
606 else
607 {
608 sprintf(dssp, "%s %s %s %s > /dev/null %s",
609 dptr, bDoAccSurf ? "" : "-na", pdbfile, tmpfile, bVerbose ? "" : "2> /dev/null");
610
611 }
612 fprintf(stderrstderr, "dssp cmd='%s'\n", dssp);
613
614 if (fnTArea)
11
←
Taking false branch
→
615 {
616 fTArea = xvgropen(fnTArea, "Solvent Accessible Surface Area",
617 output_env_get_xvgr_tlabel(oenv), "Area (nm\\S2\\N)", oenv);
618 xvgr_legend(fTArea, 2, leg, oenv);
619 }
620 else
621 {
622 fTArea = NULL((void*)0);
623 }
624
625 mat.map = NULL((void*)0);
626 mat.nmap = getcmap(libopen(opt2fn("-map", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm)),
627 opt2fn("-map", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), &(mat.map));
628
629 natoms = read_first_x(oenv, &status, ftp2fn(efTRX, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), &t, &x, box);
630 if (natoms > atoms->nr)
12
←
Taking false branch
→
631 {
632 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_do_dssp.c"
, 632, "\nTrajectory does not match topology!");
633 }
634 if (gnx > natoms)
13
←
Taking false branch
→
635 {
636 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_do_dssp.c"
, 636, "\nTrajectory does not match selected group!");
637 }
638
639 snew(average_area, atoms->nres)(average_area) = save_calloc("average_area", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_do_dssp.c"
, 639, (atoms->nres), sizeof(*(average_area)));
640 snew(av_area, atoms->nres)(av_area) = save_calloc("av_area", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_do_dssp.c"
, 640, (atoms->nres), sizeof(*(av_area)));
641 snew(norm_av_area, atoms->nres)(norm_av_area) = save_calloc("norm_av_area", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_do_dssp.c"
, 641, (atoms->nres), sizeof(*(norm_av_area)));
642 accr = NULL((void*)0);
643 naccr = 0;
644
645 gpbc = gmx_rmpbc_init(&top.idef, ePBC, natoms);
646 do
16
←
Loop condition is false. Exiting loop
→
647 {
648 t = output_env_conv_time(oenv, t);
649 if (bDoAccSurf && nframe >= naccr)
650 {
651 naccr += 10;
652 srenew(accr, naccr)(accr) = save_realloc("accr", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_do_dssp.c"
, 652, (accr), (naccr), sizeof(*(accr)));
653 for (i = naccr-10; i < naccr; i++)
654 {
655 snew(accr[i], 2*atoms->nres-1)(accr[i]) = save_calloc("accr[i]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_do_dssp.c"
, 655, (2*atoms->nres-1), sizeof(*(accr[i])));
656 }
657 }
658 gmx_rmpbc(gpbc, natoms, box, x);
659 tapein = gmx_ffopen(pdbfile, "w");
660 write_pdbfile_indexed(tapein, NULL((void*)0), atoms, x, ePBC, box, ' ', -1, gnx, index, NULL((void*)0), TRUE1);
661 gmx_ffclose(tapein);
662
663#ifdef GMX_NO_SYSTEM
664 printf("Warning-- No calls to system(3) supported on this platform.");
665 printf("Warning-- Skipping execution of 'system(\"%s\")'.", dssp);
666 exit(1);
667#else
668 if (0 != system(dssp))
14
←
Taking false branch
→
669 {
670 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_do_dssp.c"
, 670, "Failed to execute command: %s\n",
671 "Try specifying your dssp version with the -ver option.", dssp);
672 }
673#endif
674
675 /* strip_dssp returns the number of lines found in the dssp file, i.e.
676 * the number of residues plus the separator lines */
677
678 if (bDoAccSurf)
15
←
Taking false branch
→
679 {
680 accr_ptr = accr[nframe];
681 }
682
683 nres_plus_separators = strip_dssp(tmpfile, nres, bPhbres, t,
684 accr_ptr, fTArea, &mat, average_area, oenv);
685 remove(tmpfile);
686 remove(pdbfile);
687 nframe++;
688 }
689 while (read_next_x(oenv, status, &t, x, box));
690 fprintf(stderrstderr, "\n");
691 close_trj(status);
692 if (fTArea)
17
←
Taking false branch
→
693 {
694 gmx_ffclose(fTArea);
695 }
696 gmx_rmpbc_done(gpbc);
697
698 prune_ss_legend(&mat);
18
←
Calling 'prune_ss_legend'
→
31
←
Returning from 'prune_ss_legend'
→
699
700 ss = opt2FILE("-o", NFILE, fnm, "w")gmx_ffopen(opt2fn("-o", ((int)(sizeof(fnm)/sizeof((fnm)[0])))
, fnm), "w");
701 mat.flags = 0;
702 write_xpm_m(ss, mat);
703 gmx_ffclose(ss);
704
705 if (opt2bSet("-ssdump", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm))
32
←
Taking true branch
→
706 {
707 ss = opt2FILE("-ssdump", NFILE, fnm, "w")gmx_ffopen(opt2fn("-ssdump", ((int)(sizeof(fnm)/sizeof((fnm)[
0]))), fnm), "w");
708 snew(ss_str, nres+1)(ss_str) = save_calloc("ss_str", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_do_dssp.c"
, 708, (nres+1), sizeof(*(ss_str)));
709 fprintf(ss, "%d\n", nres);
710 for (j = 0; j < mat.nx; j++)
33
←
Loop condition is true. Entering loop body
→
711 {
712 for (i = 0; (i < mat.ny); i++)
34
←
Loop condition is true. Entering loop body
→
713 {
714 ss_str[i] = mat.map[mat.matrix[j][i]].code.c1;
35
←
Dereference of null pointer
715 }
716 ss_str[i] = '\0';
717 fprintf(ss, "%s\n", ss_str);
718 }
719 gmx_ffclose(ss);
720 sfree(ss_str)save_free("ss_str", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_do_dssp.c"
, 720, (ss_str));
721 }
722 analyse_ss(fnSCount, &mat, ss_string, oenv);
723
724 if (bDoAccSurf)
725 {
726 write_sas_mat(fnArea, accr, nframe, nres_plus_separators, &mat);
727
728 for (i = 0; i < atoms->nres; i++)
729 {
730 av_area[i] = (average_area[i] / (real)nframe);
731 }
732
733 norm_acc(atoms, nres, av_area, norm_av_area);
734
735 if (fnAArea)
736 {
737 acc = xvgropen(fnAArea, "Average Accessible Area",
738 "Residue", "A\\S2", oenv);
739 for (i = 0; (i < nres); i++)
740 {
741 fprintf(acc, "%5d %10g %10g\n", i+1, av_area[i], norm_av_area[i]);
742 }
743 gmx_ffclose(acc);
744 }
745 }
746
747 view_all(oenv, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm);
748
749 return 0;
750}