Bug Summary

File:gromacs/gmxana/gmx_disre.c
Location:line 218, column 9
Description:Value stored to 'ener' is never read

Annotated Source Code

1/*
2 * This file is part of the GROMACS molecular simulation package.
3 *
4 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5 * Copyright (c) 2001-2004, The GROMACS development team.
6 * Copyright (c) 2013,2014, by the GROMACS development team, led by
7 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
8 * and including many others, as listed in the AUTHORS file in the
9 * top-level source directory and at http://www.gromacs.org.
10 *
11 * GROMACS is free software; you can redistribute it and/or
12 * modify it under the terms of the GNU Lesser General Public License
13 * as published by the Free Software Foundation; either version 2.1
14 * of the License, or (at your option) any later version.
15 *
16 * GROMACS is distributed in the hope that it will be useful,
17 * but WITHOUT ANY WARRANTY; without even the implied warranty of
18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19 * Lesser General Public License for more details.
20 *
21 * You should have received a copy of the GNU Lesser General Public
22 * License along with GROMACS; if not, see
23 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
24 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
25 *
26 * If you want to redistribute modifications to GROMACS, please
27 * consider that scientific software is very special. Version
28 * control is crucial - bugs must be traceable. We will be happy to
29 * consider code for inclusion in the official distribution, but
30 * derived work must not be called official GROMACS. Details are found
31 * in the README & COPYING files - if they are missing, get the
32 * official version at http://www.gromacs.org.
33 *
34 * To help us fund GROMACS development, we humbly ask that you cite
35 * the research papers on the package. Check out http://www.gromacs.org.
36 */
37#ifdef HAVE_CONFIG_H1
38#include <config.h>
39#endif
40
41#include <math.h>
42#include <stdlib.h>
43#include <string.h>
44
45#include "typedefs.h"
46#include "macros.h"
47#include "mshift.h"
48#include "gromacs/fileio/xvgr.h"
49#include "viewit.h"
50#include "gromacs/math/vec.h"
51#include "gromacs/fileio/confio.h"
52#include "gromacs/utility/smalloc.h"
53#include "nrnb.h"
54#include "disre.h"
55#include "gromacs/commandline/pargs.h"
56#include "force.h"
57#include "gstat.h"
58#include "main.h"
59#include "gromacs/fileio/pdbio.h"
60#include "index.h"
61#include "mdatoms.h"
62#include "gromacs/fileio/tpxio.h"
63#include "gromacs/fileio/trxio.h"
64#include "mdrun.h"
65#include "names.h"
66#include "gromacs/fileio/matio.h"
67#include "mtop_util.h"
68#include "gmx_ana.h"
69
70#include "gromacs/math/do_fit.h"
71#include "gromacs/utility/fatalerror.h"
72
73typedef struct {
74 int n;
75 real v;
76} t_toppop;
77
78t_toppop *top = NULL((void*)0);
79int ntop = 0;
80
81typedef struct {
82 int nv, nframes;
83 real sumv, averv, maxv;
84 real *aver1, *aver2, *aver_3, *aver_6;
85} t_dr_result;
86
87static void init5(int n)
88{
89 ntop = n;
90 snew(top, ntop)(top) = save_calloc("top", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_disre.c"
, 90, (ntop), sizeof(*(top)));
91}
92
93static void reset5(void)
94{
95 int i;
96
97 for (i = 0; (i < ntop); i++)
98 {
99 top[i].n = -1;
100 top[i].v = 0;
101 }
102}
103
104int tpcomp(const void *a, const void *b)
105{
106 t_toppop *tpa;
107 t_toppop *tpb;
108
109 tpa = (t_toppop *)a;
110 tpb = (t_toppop *)b;
111
112 return (1e7*(tpb->v-tpa->v));
113}
114
115static void add5(int ndr, real viol)
116{
117 int i, mini;
118
119 mini = 0;
120 for (i = 1; (i < ntop); i++)
121 {
122 if (top[i].v < top[mini].v)
123 {
124 mini = i;
125 }
126 }
127 if (viol > top[mini].v)
128 {
129 top[mini].v = viol;
130 top[mini].n = ndr;
131 }
132}
133
134static void print5(FILE *fp)
135{
136 int i;
137
138 qsort(top, ntop, sizeof(top[0]), tpcomp);
139 fprintf(fp, "Index:");
140 for (i = 0; (i < ntop); i++)
141 {
142 fprintf(fp, " %6d", top[i].n);
143 }
144 fprintf(fp, "\nViol: ");
145 for (i = 0; (i < ntop); i++)
146 {
147 fprintf(fp, " %6.3f", top[i].v);
148 }
149 fprintf(fp, "\n");
150}
151
152static void check_viol(FILE *log,
153 t_ilist *disres, t_iparams forceparams[],
154 rvec x[], rvec f[],
155 t_pbc *pbc, t_graph *g, t_dr_result dr[],
156 int clust_id, int isize, atom_id index[], real vvindex[],
157 t_fcdata *fcd)
158{
159 t_iatom *forceatoms;
160 int i, j, nat, n, type, nviol, ndr, label;
161 real ener, rt, mviol, tviol, viol, lam, dvdl, drt;
162 rvec *fshift;
163 static gmx_bool bFirst = TRUE1;
164
165 lam = 0;
166 dvdl = 0;
167 tviol = 0;
168 nviol = 0;
169 mviol = 0;
170 ndr = 0;
171 if (ntop)
172 {
173 reset5();
174 }
175 forceatoms = disres->iatoms;
176 for (j = 0; (j < isize); j++)
177 {
178 vvindex[j] = 0;
179 }
180 nat = interaction_function[F_DISRES].nratoms+1;
181 for (i = 0; (i < disres->nr); )
182 {
183 type = forceatoms[i];
184 n = 0;
185 label = forceparams[type].disres.label;
186 if (debug)
187 {
188 fprintf(debug, "DISRE: ndr = %d, label = %d i=%d, n =%d\n",
189 ndr, label, i, n);
190 }
191 if (ndr != label)
192 {
193 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_disre.c"
, 193, "tpr inconsistency. ndr = %d, label = %d\n", ndr, label);
194 }
195 do
196 {
197 n += nat;
198 }
199 while (((i+n) < disres->nr) &&
200 (forceparams[forceatoms[i+n]].disres.label == label));
201
202 calc_disres_R_6(n, &forceatoms[i], forceparams,
203 (const rvec*)x, pbc, fcd, NULL((void*)0));
204
205 if (fcd->disres.Rt_6[0] <= 0)
206 {
207 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_disre.c"
, 207, "ndr = %d, rt_6 = %f", ndr, fcd->disres.Rt_6[0]);
208 }
209
210 rt = pow(fcd->disres.Rt_6[0], -1.0/6.0);
211 dr[clust_id].aver1[ndr] += rt;
212 dr[clust_id].aver2[ndr] += sqr(rt);
213 drt = pow(rt, -3.0);
214 dr[clust_id].aver_3[ndr] += drt;
215 dr[clust_id].aver_6[ndr] += fcd->disres.Rt_6[0];
216
217 snew(fshift, SHIFTS)(fshift) = save_calloc("fshift", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_disre.c"
, 217, (((2*1 +1)*(2*1 +1)*(2*2 +1))), sizeof(*(fshift)));
218 ener = interaction_function[F_DISRES].ifunc(n, &forceatoms[i], forceparams,
Value stored to 'ener' is never read
219 (const rvec*)x, f, fshift,
220 pbc, g, lam, &dvdl, NULL((void*)0), fcd, NULL((void*)0));
221 sfree(fshift)save_free("fshift", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_disre.c"
, 221, (fshift));
222 viol = fcd->disres.sumviol;
223
224 if (viol > 0)
225 {
226 nviol++;
227 if (ntop)
228 {
229 add5(forceparams[type].disres.label, viol);
230 }
231 if (viol > mviol)
232 {
233 mviol = viol;
234 }
235 tviol += viol;
236 for (j = 0; (j < isize); j++)
237 {
238 if (index[j] == forceparams[type].disres.label)
239 {
240 vvindex[j] = pow(fcd->disres.Rt_6[0], -1.0/6.0);
241 }
242 }
243 }
244 ndr++;
245 i += n;
246 }
247 dr[clust_id].nv = nviol;
248 dr[clust_id].maxv = mviol;
249 dr[clust_id].sumv = tviol;
250 dr[clust_id].averv = tviol/ndr;
251 dr[clust_id].nframes++;
252
253 if (bFirst)
254 {
255 fprintf(stderrstderr, "\nThere are %d restraints and %d pairs\n",
256 ndr, disres->nr/nat);
257 bFirst = FALSE0;
258 }
259 if (ntop)
260 {
261 print5(log);
262 }
263}
264
265typedef struct {
266 int index;
267 gmx_bool bCore;
268 real up1, r, rT3, rT6, viol, violT3, violT6;
269} t_dr_stats;
270
271static int drs_comp(const void *a, const void *b)
272{
273 t_dr_stats *da, *db;
274
275 da = (t_dr_stats *)a;
276 db = (t_dr_stats *)b;
277
278 if (da->viol > db->viol)
279 {
280 return -1;
281 }
282 else if (da->viol < db->viol)
283 {
284 return 1;
285 }
286 else
287 {
288 return 0;
289 }
290}
291
292static void dump_dump(FILE *log, int ndr, t_dr_stats drs[])
293{
294 static const char *core[] = { "All restraints", "Core restraints" };
295 static const char *tp[] = { "linear", "third power", "sixth power" };
296 real viol_tot, viol_max, viol = 0;
297 gmx_bool bCore;
298 int nviol, nrestr;
299 int i, kkk;
300
301 for (bCore = FALSE0; (bCore <= TRUE1); bCore++)
302 {
303 for (kkk = 0; (kkk < 3); kkk++)
304 {
305 viol_tot = 0;
306 viol_max = 0;
307 nviol = 0;
308 nrestr = 0;
309 for (i = 0; (i < ndr); i++)
310 {
311 if (!bCore || (bCore && drs[i].bCore))
312 {
313 switch (kkk)
314 {
315 case 0:
316 viol = drs[i].viol;
317 break;
318 case 1:
319 viol = drs[i].violT3;
320 break;
321 case 2:
322 viol = drs[i].violT6;
323 break;
324 default:
325 gmx_incons("Dumping violations")_gmx_error("incons", "Dumping violations", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_disre.c"
, 325);
326 }
327 viol_max = max(viol_max, viol)(((viol_max) > (viol)) ? (viol_max) : (viol) );
328 if (viol > 0)
329 {
330 nviol++;
331 }
332 viol_tot += viol;
333 nrestr++;
334 }
335 }
336 if ((nrestr > 0) || (bCore && (nrestr < ndr)))
337 {
338 fprintf(log, "\n");
339 fprintf(log, "+++++++ %s ++++++++\n", core[bCore]);
340 fprintf(log, "+++++++ Using %s averaging: ++++++++\n", tp[kkk]);
341 fprintf(log, "Sum of violations: %8.3f nm\n", viol_tot);
342 if (nrestr > 0)
343 {
344 fprintf(log, "Average violation: %8.3f nm\n", viol_tot/nrestr);
345 }
346 fprintf(log, "Largest violation: %8.3f nm\n", viol_max);
347 fprintf(log, "Number of violated restraints: %d/%d\n", nviol, nrestr);
348 }
349 }
350 }
351}
352
353static void dump_viol(FILE *log, int ndr, t_dr_stats *drs, gmx_bool bLinear)
354{
355 int i;
356
357 fprintf(log, "Restr. Core Up1 <r> <rT3> <rT6> <viol><violT3><violT6>\n");
358 for (i = 0; (i < ndr); i++)
359 {
360 if (bLinear && (drs[i].viol == 0))
361 {
362 break;
363 }
364 fprintf(log, "%6d%5s%8.3f%8.3f%8.3f%8.3f%8.3f%8.3f%8.3f\n",
365 drs[i].index, yesno_names[drs[i].bCore],
366 drs[i].up1, drs[i].r, drs[i].rT3, drs[i].rT6,
367 drs[i].viol, drs[i].violT3, drs[i].violT6);
368 }
369}
370
371static gmx_bool is_core(int i, int isize, atom_id index[])
372{
373 int kk;
374 gmx_bool bIC = FALSE0;
375
376 for (kk = 0; !bIC && (kk < isize); kk++)
377 {
378 bIC = (index[kk] == i);
379 }
380
381 return bIC;
382}
383
384static void dump_stats(FILE *log, int nsteps, int ndr, t_ilist *disres,
385 t_iparams ip[], t_dr_result *dr,
386 int isize, atom_id index[], t_atoms *atoms)
387{
388 int i, j, nra;
389 t_dr_stats *drs;
390
391 fprintf(log, "\n");
392 fprintf(log, "++++++++++++++ STATISTICS ++++++++++++++++++++++++\n");
393 snew(drs, ndr)(drs) = save_calloc("drs", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_disre.c"
, 393, (ndr), sizeof(*(drs)));
394 j = 0;
395 nra = interaction_function[F_DISRES].nratoms+1;
396 for (i = 0; (i < ndr); i++)
397 {
398 /* Search for the appropriate restraint */
399 for (; (j < disres->nr); j += nra)
400 {
401 if (ip[disres->iatoms[j]].disres.label == i)
402 {
403 break;
404 }
405 }
406 drs[i].index = i;
407 drs[i].bCore = is_core(i, isize, index);
408 drs[i].up1 = ip[disres->iatoms[j]].disres.up1;
409 drs[i].r = dr->aver1[i]/nsteps;
410 drs[i].rT3 = pow(dr->aver_3[i]/nsteps, -1.0/3.0);
411 drs[i].rT6 = pow(dr->aver_6[i]/nsteps, -1.0/6.0);
412 drs[i].viol = max(0, drs[i].r-drs[i].up1)(((0) > (drs[i].r-drs[i].up1)) ? (0) : (drs[i].r-drs[i].up1
) );
413 drs[i].violT3 = max(0, drs[i].rT3-drs[i].up1)(((0) > (drs[i].rT3-drs[i].up1)) ? (0) : (drs[i].rT3-drs[i
].up1) );
414 drs[i].violT6 = max(0, drs[i].rT6-drs[i].up1)(((0) > (drs[i].rT6-drs[i].up1)) ? (0) : (drs[i].rT6-drs[i
].up1) );
415 if (atoms)
416 {
417 int j1 = disres->iatoms[j+1];
418 int j2 = disres->iatoms[j+2];
419 atoms->pdbinfo[j1].bfac += drs[i].violT3*5;
420 atoms->pdbinfo[j2].bfac += drs[i].violT3*5;
421 }
422 }
423 dump_viol(log, ndr, drs, FALSE0);
424
425 fprintf(log, "+++ Sorted by linear averaged violations: +++\n");
426 qsort(drs, ndr, sizeof(drs[0]), drs_comp);
427 dump_viol(log, ndr, drs, TRUE1);
428
429 dump_dump(log, ndr, drs);
430
431 sfree(drs)save_free("drs", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_disre.c"
, 431, (drs));
432}
433
434static void dump_clust_stats(FILE *fp, int ndr, t_ilist *disres,
435 t_iparams ip[], t_blocka *clust, t_dr_result dr[],
436 char *clust_name[], int isize, atom_id index[])
437{
438 int i, j, k, nra, mmm = 0;
439 double sumV, maxV, sumVT3, sumVT6, maxVT3, maxVT6;
440 t_dr_stats *drs;
441
442 fprintf(fp, "\n");
443 fprintf(fp, "++++++++++++++ STATISTICS ++++++++++++++++++++++\n");
444 fprintf(fp, "Cluster NFrames SumV MaxV SumVT MaxVT SumVS MaxVS\n");
445
446 snew(drs, ndr)(drs) = save_calloc("drs", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_disre.c"
, 446, (ndr), sizeof(*(drs)));
447
448 for (k = 0; (k < clust->nr); k++)
449 {
450 if (dr[k].nframes == 0)
451 {
452 continue;
453 }
454 if (dr[k].nframes != (clust->index[k+1]-clust->index[k]))
455 {
456 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_disre.c"
, 456, "Inconsistency in cluster %s.\n"
457 "Found %d frames in trajectory rather than the expected %d\n",
458 clust_name[k], dr[k].nframes,
459 clust->index[k+1]-clust->index[k]);
460 }
461 if (!clust_name[k])
462 {
463 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_disre.c"
, 463, "Inconsistency with cluster %d. Invalid name", k);
464 }
465 j = 0;
466 nra = interaction_function[F_DISRES].nratoms+1;
467 sumV = sumVT3 = sumVT6 = maxV = maxVT3 = maxVT6 = 0;
468 for (i = 0; (i < ndr); i++)
469 {
470 /* Search for the appropriate restraint */
471 for (; (j < disres->nr); j += nra)
472 {
473 if (ip[disres->iatoms[j]].disres.label == i)
474 {
475 break;
476 }
477 }
478 drs[i].index = i;
479 drs[i].bCore = is_core(i, isize, index);
480 drs[i].up1 = ip[disres->iatoms[j]].disres.up1;
481 drs[i].r = dr[k].aver1[i]/dr[k].nframes;
482 if ((dr[k].aver_3[i] <= 0) || (dr[k].aver_3[i] != dr[k].aver_3[i]))
483 {
484 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_disre.c"
, 484, "dr[%d].aver_3[%d] = %f", k, i, dr[k].aver_3[i]);
485 }
486 drs[i].rT3 = pow(dr[k].aver_3[i]/dr[k].nframes, -1.0/3.0);
487 drs[i].rT6 = pow(dr[k].aver_6[i]/dr[k].nframes, -1.0/6.0);
488 drs[i].viol = max(0, drs[i].r-drs[i].up1)(((0) > (drs[i].r-drs[i].up1)) ? (0) : (drs[i].r-drs[i].up1
) );
489 drs[i].violT3 = max(0, drs[i].rT3-drs[i].up1)(((0) > (drs[i].rT3-drs[i].up1)) ? (0) : (drs[i].rT3-drs[i
].up1) );
490 drs[i].violT6 = max(0, drs[i].rT6-drs[i].up1)(((0) > (drs[i].rT6-drs[i].up1)) ? (0) : (drs[i].rT6-drs[i
].up1) );
491 sumV += drs[i].viol;
492 sumVT3 += drs[i].violT3;
493 sumVT6 += drs[i].violT6;
494 maxV = max(maxV, drs[i].viol)(((maxV) > (drs[i].viol)) ? (maxV) : (drs[i].viol) );
495 maxVT3 = max(maxVT3, drs[i].violT3)(((maxVT3) > (drs[i].violT3)) ? (maxVT3) : (drs[i].violT3)
);
496 maxVT6 = max(maxVT6, drs[i].violT6)(((maxVT6) > (drs[i].violT6)) ? (maxVT6) : (drs[i].violT6)
);
497 }
498 if (strcmp(clust_name[k], "1000")__extension__ ({ size_t __s1_len, __s2_len; (__builtin_constant_p
(clust_name[k]) && __builtin_constant_p ("1000") &&
(__s1_len = strlen (clust_name[k]), __s2_len = strlen ("1000"
), (!((size_t)(const void *)((clust_name[k]) + 1) - (size_t)(
const void *)(clust_name[k]) == 1) || __s1_len >= 4) &&
(!((size_t)(const void *)(("1000") + 1) - (size_t)(const void
*)("1000") == 1) || __s2_len >= 4)) ? __builtin_strcmp (clust_name
[k], "1000") : (__builtin_constant_p (clust_name[k]) &&
((size_t)(const void *)((clust_name[k]) + 1) - (size_t)(const
void *)(clust_name[k]) == 1) && (__s1_len = strlen (
clust_name[k]), __s1_len < 4) ? (__builtin_constant_p ("1000"
) && ((size_t)(const void *)(("1000") + 1) - (size_t)
(const void *)("1000") == 1) ? __builtin_strcmp (clust_name[k
], "1000") : (__extension__ ({ const unsigned char *__s2 = (const
unsigned char *) (const char *) ("1000"); int __result = (((
const unsigned char *) (const char *) (clust_name[k]))[0] - __s2
[0]); if (__s1_len > 0 && __result == 0) { __result
= (((const unsigned char *) (const char *) (clust_name[k]))[
1] - __s2[1]); if (__s1_len > 1 && __result == 0) {
__result = (((const unsigned char *) (const char *) (clust_name
[k]))[2] - __s2[2]); if (__s1_len > 2 && __result ==
0) __result = (((const unsigned char *) (const char *) (clust_name
[k]))[3] - __s2[3]); } } __result; }))) : (__builtin_constant_p
("1000") && ((size_t)(const void *)(("1000") + 1) - (
size_t)(const void *)("1000") == 1) && (__s2_len = strlen
("1000"), __s2_len < 4) ? (__builtin_constant_p (clust_name
[k]) && ((size_t)(const void *)((clust_name[k]) + 1) -
(size_t)(const void *)(clust_name[k]) == 1) ? __builtin_strcmp
(clust_name[k], "1000") : (- (__extension__ ({ const unsigned
char *__s2 = (const unsigned char *) (const char *) (clust_name
[k]); int __result = (((const unsigned char *) (const char *)
("1000"))[0] - __s2[0]); if (__s2_len > 0 && __result
== 0) { __result = (((const unsigned char *) (const char *) (
"1000"))[1] - __s2[1]); if (__s2_len > 1 && __result
== 0) { __result = (((const unsigned char *) (const char *) (
"1000"))[2] - __s2[2]); if (__s2_len > 2 && __result
== 0) __result = (((const unsigned char *) (const char *) ("1000"
))[3] - __s2[3]); } } __result; })))) : __builtin_strcmp (clust_name
[k], "1000")))); }) == 0)
499 {
500 mmm++;
501 }
502 fprintf(fp, "%-10s%6d%8.3f %8.3f %8.3f %8.3f %8.3f %8.3f\n",
503 clust_name[k],
504 dr[k].nframes, sumV, maxV, sumVT3, maxVT3, sumVT6, maxVT6);
505
506 }
507 fflush(fp);
508 sfree(drs)save_free("drs", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_disre.c"
, 508, (drs));
509}
510
511static void init_dr_res(t_dr_result *dr, int ndr)
512{
513 snew(dr->aver1, ndr+1)(dr->aver1) = save_calloc("dr->aver1", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_disre.c"
, 513, (ndr+1), sizeof(*(dr->aver1)));
514 snew(dr->aver2, ndr+1)(dr->aver2) = save_calloc("dr->aver2", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_disre.c"
, 514, (ndr+1), sizeof(*(dr->aver2)));
515 snew(dr->aver_3, ndr+1)(dr->aver_3) = save_calloc("dr->aver_3", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_disre.c"
, 515, (ndr+1), sizeof(*(dr->aver_3)));
516 snew(dr->aver_6, ndr+1)(dr->aver_6) = save_calloc("dr->aver_6", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_disre.c"
, 516, (ndr+1), sizeof(*(dr->aver_6)));
517 dr->nv = 0;
518 dr->nframes = 0;
519 dr->sumv = 0;
520 dr->maxv = 0;
521 dr->averv = 0;
522}
523
524static void dump_disre_matrix(const char *fn, t_dr_result *dr, int ndr,
525 int nsteps, t_idef *idef, gmx_mtop_t *mtop,
526 real max_dr, int nlevels, gmx_bool bThird)
527{
528 FILE *fp;
529 int *resnr;
530 int n_res, a_offset, mb, mol, a;
531 t_atoms *atoms;
532 int iii, i, j, nra, nratoms, tp, ri, rj, index, nlabel, label;
533 atom_id ai, aj, *ptr;
534 real **matrix, *t_res, hi, *w_dr, rav, rviol;
535 t_rgb rlo = { 1, 1, 1 };
536 t_rgb rhi = { 0, 0, 0 };
537 if (fn == NULL((void*)0))
538 {
539 return;
540 }
541
542 snew(resnr, mtop->natoms)(resnr) = save_calloc("resnr", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_disre.c"
, 542, (mtop->natoms), sizeof(*(resnr)));
543 n_res = 0;
544 a_offset = 0;
545 for (mb = 0; mb < mtop->nmolblock; mb++)
546 {
547 atoms = &mtop->moltype[mtop->molblock[mb].type].atoms;
548 for (mol = 0; mol < mtop->molblock[mb].nmol; mol++)
549 {
550 for (a = 0; a < atoms->nr; a++)
551 {
552 resnr[a_offset+a] = n_res + atoms->atom[a].resind;
553 }
554 n_res += atoms->nres;
555 a_offset += atoms->nr;
556 }
557 }
558
559 snew(t_res, n_res)(t_res) = save_calloc("t_res", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_disre.c"
, 559, (n_res), sizeof(*(t_res)));
560 for (i = 0; (i < n_res); i++)
561 {
562 t_res[i] = i+1;
563 }
564 snew(matrix, n_res)(matrix) = save_calloc("matrix", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_disre.c"
, 564, (n_res), sizeof(*(matrix)));
565 for (i = 0; (i < n_res); i++)
566 {
567 snew(matrix[i], n_res)(matrix[i]) = save_calloc("matrix[i]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_disre.c"
, 567, (n_res), sizeof(*(matrix[i])));
568 }
569 nratoms = interaction_function[F_DISRES].nratoms;
570 nra = (idef->il[F_DISRES].nr/(nratoms+1));
571 snew(ptr, nra+1)(ptr) = save_calloc("ptr", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_disre.c"
, 571, (nra+1), sizeof(*(ptr)));
572 index = 0;
573 nlabel = 0;
574 ptr[0] = 0;
575 snew(w_dr, ndr)(w_dr) = save_calloc("w_dr", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_disre.c"
, 575, (ndr), sizeof(*(w_dr)));
576 for (i = 0; (i < idef->il[F_DISRES].nr); i += nratoms+1)
577 {
578 tp = idef->il[F_DISRES].iatoms[i];
579 label = idef->iparams[tp].disres.label;
580
581 if (label != index)
582 {
583 /* Set index pointer */
584 ptr[index+1] = i;
585 if (nlabel <= 0)
586 {
587 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_disre.c"
, 587, "nlabel is %d, label = %d", nlabel, label);
588 }
589 if (index >= ndr)
590 {
591 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_disre.c"
, 591, "ndr = %d, index = %d", ndr, index);
592 }
593 /* Update the weight */
594 w_dr[index] = 1.0/nlabel;
595 index = label;
596 nlabel = 1;
597 }
598 else
599 {
600 nlabel++;
601 }
602 }
603 printf("nlabel = %d, index = %d, ndr = %d\n", nlabel, index, ndr);
604 hi = 0;
605 for (i = 0; (i < ndr); i++)
606 {
607 for (j = ptr[i]; (j < ptr[i+1]); j += nratoms+1)
608 {
609 tp = idef->il[F_DISRES].iatoms[j];
610 ai = idef->il[F_DISRES].iatoms[j+1];
611 aj = idef->il[F_DISRES].iatoms[j+2];
612
613 ri = resnr[ai];
614 rj = resnr[aj];
615 if (bThird)
616 {
617 rav = pow(dr->aver_3[i]/nsteps, -1.0/3.0);
618 }
619 else
620 {
621 rav = dr->aver1[i]/nsteps;
622 }
623 if (debug)
624 {
625 fprintf(debug, "DR %d, atoms %d, %d, distance %g\n", i, ai, aj, rav);
626 }
627 rviol = max(0, rav-idef->iparams[tp].disres.up1)(((0) > (rav-idef->iparams[tp].disres.up1)) ? (0) : (rav
-idef->iparams[tp].disres.up1) );
628 matrix[ri][rj] += w_dr[i]*rviol;
629 matrix[rj][ri] += w_dr[i]*rviol;
630 hi = max(hi, matrix[ri][rj])(((hi) > (matrix[ri][rj])) ? (hi) : (matrix[ri][rj]) );
631 hi = max(hi, matrix[rj][ri])(((hi) > (matrix[rj][ri])) ? (hi) : (matrix[rj][ri]) );
632 }
633 }
634
635 sfree(resnr)save_free("resnr", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_disre.c"
, 635, (resnr));
636
637 if (max_dr > 0)
638 {
639 if (hi > max_dr)
640 {
641 printf("Warning: the maxdr that you have specified (%g) is smaller than\nthe largest value in your simulation (%g)\n", max_dr, hi);
642 }
643 hi = max_dr;
644 }
645 printf("Highest level in the matrix will be %g\n", hi);
646 fp = gmx_ffopen(fn, "w");
647 write_xpm(fp, 0, "Distance Violations", "<V> (nm)", "Residue", "Residue",
648 n_res, n_res, t_res, t_res, matrix, 0, hi, rlo, rhi, &nlevels);
649 gmx_ffclose(fp);
650}
651
652int gmx_disre(int argc, char *argv[])
653{
654 const char *desc[] = {
655 "[THISMODULE] computes violations of distance restraints.",
656 "If necessary, all protons can be added to a protein molecule ",
657 "using the [gmx-protonate] program.[PAR]",
658 "The program always",
659 "computes the instantaneous violations rather than time-averaged,",
660 "because this analysis is done from a trajectory file afterwards",
661 "it does not make sense to use time averaging. However,",
662 "the time averaged values per restraint are given in the log file.[PAR]",
663 "An index file may be used to select specific restraints for",
664 "printing.[PAR]",
665 "When the optional [TT]-q[tt] flag is given a [TT].pdb[tt] file coloured by the",
666 "amount of average violations.[PAR]",
667 "When the [TT]-c[tt] option is given, an index file will be read",
668 "containing the frames in your trajectory corresponding to the clusters",
669 "(defined in another manner) that you want to analyze. For these clusters",
670 "the program will compute average violations using the third power",
671 "averaging algorithm and print them in the log file."
672 };
673 static int ntop = 0;
674 static int nlevels = 20;
675 static real max_dr = 0;
676 static gmx_bool bThird = TRUE1;
677 t_pargs pa[] = {
678 { "-ntop", FALSE0, etINT, {&ntop},
679 "Number of large violations that are stored in the log file every step" },
680 { "-maxdr", FALSE0, etREAL, {&max_dr},
681 "Maximum distance violation in matrix output. If less than or equal to 0 the maximum will be determined by the data." },
682 { "-nlevels", FALSE0, etINT, {&nlevels},
683 "Number of levels in the matrix output" },
684 { "-third", FALSE0, etBOOL, {&bThird},
685 "Use inverse third power averaging or linear for matrix output" }
686 };
687
688 FILE *out = NULL((void*)0), *aver = NULL((void*)0), *numv = NULL((void*)0), *maxxv = NULL((void*)0), *xvg = NULL((void*)0);
689 t_tpxheader header;
690 t_inputrec ir;
691 gmx_mtop_t mtop;
692 rvec *xtop;
693 gmx_localtop_t *top;
694 t_atoms *atoms = NULL((void*)0);
695 t_fcdata fcd;
696 t_nrnb nrnb;
697 t_graph *g;
698 int ntopatoms, natoms, i, j, kkk;
699 t_trxstatus *status;
700 real t;
701 rvec *x, *f, *xav = NULL((void*)0);
702 matrix box;
703 gmx_bool bPDB;
704 int isize;
705 atom_id *index = NULL((void*)0), *ind_fit = NULL((void*)0);
706 char *grpname;
707 t_cluster_ndx *clust = NULL((void*)0);
708 t_dr_result dr, *dr_clust = NULL((void*)0);
709 char **leg;
710 real *vvindex = NULL((void*)0), *w_rls = NULL((void*)0);
711 t_mdatoms *mdatoms;
712 t_pbc pbc, *pbc_null;
713 int my_clust;
714 FILE *fplog;
715 output_env_t oenv;
716 gmx_rmpbc_t gpbc = NULL((void*)0);
717
718 t_filenm fnm[] = {
719 { efTPX, NULL((void*)0), NULL((void*)0), ffREAD1<<1 },
720 { efTRX, "-f", NULL((void*)0), ffREAD1<<1 },
721 { efXVG, "-ds", "drsum", ffWRITE1<<2 },
722 { efXVG, "-da", "draver", ffWRITE1<<2 },
723 { efXVG, "-dn", "drnum", ffWRITE1<<2 },
724 { efXVG, "-dm", "drmax", ffWRITE1<<2 },
725 { efXVG, "-dr", "restr", ffWRITE1<<2 },
726 { efLOG, "-l", "disres", ffWRITE1<<2 },
727 { efNDX, NULL((void*)0), "viol", ffOPTRD(1<<1 | 1<<3) },
728 { efPDB, "-q", "viol", ffOPTWR(1<<2| 1<<3) },
729 { efNDX, "-c", "clust", ffOPTRD(1<<1 | 1<<3) },
730 { efXPM, "-x", "matrix", ffOPTWR(1<<2| 1<<3) }
731 };
732#define NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))) asize(fnm)((int)(sizeof(fnm)/sizeof((fnm)[0])))
733
734 if (!parse_common_args(&argc, argv, PCA_CAN_TIME((1<<6) | (1<<7) | (1<<14)) | PCA_CAN_VIEW(1<<5) | PCA_BE_NICE(1<<13),
735 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))
736 {
737 return 0;
738 }
739
740 fplog = ftp2FILE(efLOG, NFILE, fnm, "w")gmx_ffopen(ftp2fn(efLOG, ((int)(sizeof(fnm)/sizeof((fnm)[0]))
), fnm), "w");
741
742 if (ntop)
743 {
744 init5(ntop);
745 }
746
747 read_tpxheader(ftp2fn(efTPX, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), &header, FALSE0, NULL((void*)0), NULL((void*)0));
748 snew(xtop, header.natoms)(xtop) = save_calloc("xtop", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_disre.c"
, 748, (header.natoms), sizeof(*(xtop)));
749 read_tpx(ftp2fn(efTPX, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), &ir, box, &ntopatoms, xtop, NULL((void*)0), NULL((void*)0), &mtop);
750 bPDB = opt2bSet("-q", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm);
751 if (bPDB)
752 {
753 snew(xav, ntopatoms)(xav) = save_calloc("xav", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_disre.c"
, 753, (ntopatoms), sizeof(*(xav)));
754 snew(ind_fit, ntopatoms)(ind_fit) = save_calloc("ind_fit", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_disre.c"
, 754, (ntopatoms), sizeof(*(ind_fit)));
755 snew(w_rls, ntopatoms)(w_rls) = save_calloc("w_rls", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_disre.c"
, 755, (ntopatoms), sizeof(*(w_rls)));
756 for (kkk = 0; (kkk < ntopatoms); kkk++)
757 {
758 w_rls[kkk] = 1;
759 ind_fit[kkk] = kkk;
760 }
761
762 snew(atoms, 1)(atoms) = save_calloc("atoms", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_disre.c"
, 762, (1), sizeof(*(atoms)));
763 *atoms = gmx_mtop_global_atoms(&mtop);
764
765 if (atoms->pdbinfo == NULL((void*)0))
766 {
767 snew(atoms->pdbinfo, atoms->nr)(atoms->pdbinfo) = save_calloc("atoms->pdbinfo", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_disre.c"
, 767, (atoms->nr), sizeof(*(atoms->pdbinfo)));
768 }
769 }
770
771 top = gmx_mtop_generate_local_top(&mtop, &ir);
772
773 g = NULL((void*)0);
774 pbc_null = NULL((void*)0);
775 if (ir.ePBC != epbcNONE)
776 {
777 if (ir.bPeriodicMols)
778 {
779 pbc_null = &pbc;
780 }
781 else
782 {
783 g = mk_graph(fplog, &top->idef, 0, mtop.natoms, FALSE0, FALSE0);
784 }
785 }
786
787 if (ftp2bSet(efNDX, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm))
788 {
789 rd_index(ftp2fn(efNDX, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), 1, &isize, &index, &grpname);
790 xvg = xvgropen(opt2fn("-dr", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), "Individual Restraints", "Time (ps)",
791 "nm", oenv);
792 snew(vvindex, isize)(vvindex) = save_calloc("vvindex", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_disre.c"
, 792, (isize), sizeof(*(vvindex)));
793 snew(leg, isize)(leg) = save_calloc("leg", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_disre.c"
, 793, (isize), sizeof(*(leg)));
794 for (i = 0; (i < isize); i++)
795 {
796 index[i]++;
797 snew(leg[i], 12)(leg[i]) = save_calloc("leg[i]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_disre.c"
, 797, (12), sizeof(*(leg[i])));
798 sprintf(leg[i], "index %d", index[i]);
799 }
800 xvgr_legend(xvg, isize, (const char**)leg, oenv);
801 }
802 else
803 {
804 isize = 0;
805 }
806
807 ir.dr_tau = 0.0;
808 init_disres(fplog, &mtop, &ir, NULL((void*)0), &fcd, NULL((void*)0), FALSE0);
809
810 natoms = read_first_x(oenv, &status, ftp2fn(efTRX, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), &t, &x, box);
811 snew(f, 5*natoms)(f) = save_calloc("f", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_disre.c"
, 811, (5*natoms), sizeof(*(f)));
812
813 init_dr_res(&dr, fcd.disres.nres);
814 if (opt2bSet("-c", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm))
815 {
816 clust = cluster_index(fplog, opt2fn("-c", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm));
817 snew(dr_clust, clust->clust->nr+1)(dr_clust) = save_calloc("dr_clust", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_disre.c"
, 817, (clust->clust->nr+1), sizeof(*(dr_clust)));
818 for (i = 0; (i <= clust->clust->nr); i++)
819 {
820 init_dr_res(&dr_clust[i], fcd.disres.nres);
821 }
822 }
823 else
824 {
825 out = xvgropen(opt2fn("-ds", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm),
826 "Sum of Violations", "Time (ps)", "nm", oenv);
827 aver = xvgropen(opt2fn("-da", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm),
828 "Average Violation", "Time (ps)", "nm", oenv);
829 numv = xvgropen(opt2fn("-dn", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm),
830 "# Violations", "Time (ps)", "#", oenv);
831 maxxv = xvgropen(opt2fn("-dm", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm),
832 "Largest Violation", "Time (ps)", "nm", oenv);
833 }
834
835 mdatoms = init_mdatoms(fplog, &mtop, ir.efep != efepNO);
836 atoms2md(&mtop, &ir, 0, NULL((void*)0), mtop.natoms, mdatoms);
837 update_mdatoms(mdatoms, ir.fepvals->init_lambda);
838 init_nrnb(&nrnb);
839 if (ir.ePBC != epbcNONE)
840 {
841 gpbc = gmx_rmpbc_init(&top->idef, ir.ePBC, natoms);
842 }
843
844 j = 0;
845 do
846 {
847 if (ir.ePBC != epbcNONE)
848 {
849 if (ir.bPeriodicMols)
850 {
851 set_pbc(&pbc, ir.ePBC, box);
852 }
853 else
854 {
855 gmx_rmpbc(gpbc, natoms, box, x);
856 }
857 }
858
859 if (clust)
860 {
861 if (j > clust->maxframe)
862 {
863 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_disre.c"
, 863, "There are more frames in the trajectory than in the cluster index file. t = %8f\n", t);
864 }
865 my_clust = clust->inv_clust[j];
866 range_check(my_clust, 0, clust->clust->nr)_range_check(my_clust, 0, clust->clust->nr, ((void*)0),
"my_clust", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_disre.c"
, 866);
867 check_viol(fplog, &(top->idef.il[F_DISRES]),
868 top->idef.iparams,
869 x, f, pbc_null, g, dr_clust, my_clust, isize, index, vvindex, &fcd);
870 }
871 else
872 {
873 check_viol(fplog, &(top->idef.il[F_DISRES]),
874 top->idef.iparams,
875 x, f, pbc_null, g, &dr, 0, isize, index, vvindex, &fcd);
876 }
877 if (bPDB)
878 {
879 reset_x(atoms->nr, ind_fit, atoms->nr, NULL((void*)0), x, w_rls);
880 do_fit(atoms->nr, w_rls, x, x);
881 if (j == 0)
882 {
883 /* Store the first frame of the trajectory as 'characteristic'
884 * for colouring with violations.
885 */
886 for (kkk = 0; (kkk < atoms->nr); kkk++)
887 {
888 copy_rvec(x[kkk], xav[kkk]);
889 }
890 }
891 }
892 if (!clust)
893 {
894 if (isize > 0)
895 {
896 fprintf(xvg, "%10g", t);
897 for (i = 0; (i < isize); i++)
898 {
899 fprintf(xvg, " %10g", vvindex[i]);
900 }
901 fprintf(xvg, "\n");
902 }
903 fprintf(out, "%10g %10g\n", t, dr.sumv);
904 fprintf(aver, "%10g %10g\n", t, dr.averv);
905 fprintf(maxxv, "%10g %10g\n", t, dr.maxv);
906 fprintf(numv, "%10g %10d\n", t, dr.nv);
907 }
908 j++;
909 }
910 while (read_next_x(oenv, status, &t, x, box));
911 close_trj(status);
912 if (ir.ePBC != epbcNONE)
913 {
914 gmx_rmpbc_done(gpbc);
915 }
916
917 if (clust)
918 {
919 dump_clust_stats(fplog, fcd.disres.nres, &(top->idef.il[F_DISRES]),
920 top->idef.iparams, clust->clust, dr_clust,
921 clust->grpname, isize, index);
922 }
923 else
924 {
925 dump_stats(fplog, j, fcd.disres.nres, &(top->idef.il[F_DISRES]),
926 top->idef.iparams, &dr, isize, index,
927 bPDB ? atoms : NULL((void*)0));
928 if (bPDB)
929 {
930 write_sto_conf(opt2fn("-q", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm),
931 "Coloured by average violation in Angstrom",
932 atoms, xav, NULL((void*)0), ir.ePBC, box);
933 }
934 dump_disre_matrix(opt2fn_null("-x", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), &dr, fcd.disres.nres,
935 j, &top->idef, &mtop, max_dr, nlevels, bThird);
936 gmx_ffclose(out);
937 gmx_ffclose(aver);
938 gmx_ffclose(numv);
939 gmx_ffclose(maxxv);
940 if (isize > 0)
941 {
942 gmx_ffclose(xvg);
943 do_view(oenv, opt2fn("-dr", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), "-nxy");
944 }
945 do_view(oenv, opt2fn("-dn", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), "-nxy");
946 do_view(oenv, opt2fn("-da", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), "-nxy");
947 do_view(oenv, opt2fn("-ds", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), "-nxy");
948 do_view(oenv, opt2fn("-dm", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), "-nxy");
949 }
950
951 gmx_log_close(fplog);
952
953 return 0;
954}