Bug Summary

File:gromacs/tools/check.c
Location:line 233, column 5
Description:Value stored to 'devtot' 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-2013, 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 <stdio.h>
43#include <string.h>
44
45#include "macros.h"
46#include "gromacs/commandline/pargs.h"
47#include "txtdump.h"
48#include "gromacs/utility/fatalerror.h"
49#include "atomprop.h"
50#include "gromacs/math/vec.h"
51#include "pbc.h"
52#include "physics.h"
53#include "index.h"
54#include "gromacs/utility/smalloc.h"
55#include "names.h"
56#include "mtop_util.h"
57#include "gromacs/utility/futil.h"
58#include "gromacs/fileio/gmxfio.h"
59#include "gromacs/fileio/trnio.h"
60#include "gromacs/fileio/xtcio.h"
61#include "gromacs/fileio/confio.h"
62#include "gromacs/fileio/enxio.h"
63#include "gromacs/fileio/tpxio.h"
64#include "gromacs/fileio/trxio.h"
65
66#include "compare.h"
67
68typedef struct {
69 int bStep;
70 int bTime;
71 int bLambda;
72 int bX;
73 int bV;
74 int bF;
75 int bBox;
76} t_count;
77
78typedef struct {
79 float bStep;
80 float bTime;
81 float bLambda;
82 float bX;
83 float bV;
84 float bF;
85 float bBox;
86} t_fr_time;
87
88static void tpx2system(FILE *fp, gmx_mtop_t *mtop)
89{
90 int i, nmol, nvsite = 0;
91 gmx_mtop_atomloop_block_t aloop;
92 t_atom *atom;
93
94 fprintf(fp, "\\subsection{Simulation system}\n");
95 aloop = gmx_mtop_atomloop_block_init(mtop);
96 while (gmx_mtop_atomloop_block_next(aloop, &atom, &nmol))
97 {
98 if (atom->ptype == eptVSite)
99 {
100 nvsite += nmol;
101 }
102 }
103 fprintf(fp, "A system of %d molecules (%d atoms) was simulated.\n",
104 mtop->mols.nr, mtop->natoms-nvsite);
105 if (nvsite)
106 {
107 fprintf(fp, "Virtual sites were used in some of the molecules.\n");
108 }
109 fprintf(fp, "\n\n");
110}
111
112static void tpx2params(FILE *fp, t_inputrec *ir)
113{
114 fprintf(fp, "\\subsection{Simulation settings}\n");
115 fprintf(fp, "A total of %g ns were simulated with a time step of %g fs.\n",
116 ir->nsteps*ir->delta_t*0.001, 1000*ir->delta_t);
117 fprintf(fp, "Neighbor searching was performed every %d steps.\n", ir->nstlist);
118 fprintf(fp, "The %s algorithm was used for electrostatic interactions.\n",
119 EELTYPE(ir->coulombtype)((((ir->coulombtype) < 0) || ((ir->coulombtype) >=
(eelNR))) ? "UNDEFINED" : (eel_names)[ir->coulombtype]));
120 fprintf(fp, "with a cut-off of %g nm.\n", ir->rcoulomb);
121 if (ir->coulombtype == eelPME)
122 {
123 fprintf(fp, "A reciprocal grid of %d x %d x %d cells was used with %dth order B-spline interpolation.\n", ir->nkx, ir->nky, ir->nkz, ir->pme_order);
124 }
125 if (ir->rvdw > ir->rlist)
126 {
127 fprintf(fp, "A twin-range Van der Waals cut-off (%g/%g nm) was used, where the long range forces were updated during neighborsearching.\n", ir->rlist, ir->rvdw);
128 }
129 else
130 {
131 fprintf(fp, "A single cut-off of %g was used for Van der Waals interactions.\n", ir->rlist);
132 }
133 if (ir->etc != 0)
134 {
135 fprintf(fp, "Temperature coupling was done with the %s algorithm.\n",
136 etcoupl_names[ir->etc]);
137 }
138 if (ir->epc != 0)
139 {
140 fprintf(fp, "Pressure coupling was done with the %s algorithm.\n",
141 epcoupl_names[ir->epc]);
142 }
143 fprintf(fp, "\n\n");
144}
145
146static void tpx2methods(const char *tpx, const char *tex)
147{
148 FILE *fp;
149 t_tpxheader sh;
150 t_inputrec ir;
151 t_state state;
152 gmx_mtop_t mtop;
153
154 read_tpx_state(tpx, &ir, &state, NULL((void*)0), &mtop);
155 fp = gmx_fio_fopen(tex, "w");
156 fprintf(fp, "\\section{Methods}\n");
157 tpx2system(fp, &mtop);
158 tpx2params(fp, &ir);
159 gmx_fio_fclose(fp);
160}
161
162static void chk_coords(int frame, int natoms, rvec *x, matrix box, real fac, real tol)
163{
164 int i, j;
165 int nNul = 0;
166 real vol = det(box);
167
168 for (i = 0; (i < natoms); i++)
169 {
170 for (j = 0; (j < DIM3); j++)
171 {
172 if ((vol > 0) && (fabs(x[i][j]) > fac*box[j][j]))
173 {
174 printf("Warning at frame %d: coordinates for atom %d are large (%g)\n",
175 frame, i, x[i][j]);
176 }
177 }
178 if ((fabs(x[i][XX0]) < tol) &&
179 (fabs(x[i][YY1]) < tol) &&
180 (fabs(x[i][ZZ2]) < tol))
181 {
182 nNul++;
183 }
184 }
185 if (nNul > 0)
186 {
187 printf("Warning at frame %d: there are %d particles with all coordinates zero\n",
188 frame, nNul);
189 }
190}
191
192static void chk_vels(int frame, int natoms, rvec *v)
193{
194 int i, j;
195
196 for (i = 0; (i < natoms); i++)
197 {
198 for (j = 0; (j < DIM3); j++)
199 {
200 if (fabs(v[i][j]) > 500)
201 {
202 printf("Warning at frame %d. Velocities for atom %d are large (%g)\n",
203 frame, i, v[i][j]);
204 }
205 }
206 }
207}
208
209static void chk_forces(int frame, int natoms, rvec *f)
210{
211 int i, j;
212
213 for (i = 0; (i < natoms); i++)
214 {
215 for (j = 0; (j < DIM3); j++)
216 {
217 if (fabs(f[i][j]) > 10000)
218 {
219 printf("Warning at frame %d. Forces for atom %d are large (%g)\n",
220 frame, i, f[i][j]);
221 }
222 }
223 }
224}
225
226static void chk_bonds(t_idef *idef, int ePBC, rvec *x, matrix box, real tol)
227{
228 int ftype, i, k, ai, aj, type;
229 real b0, blen, deviation, devtot;
230 t_pbc pbc;
231 rvec dx;
232
233 devtot = 0;
Value stored to 'devtot' is never read
234 set_pbc(&pbc, ePBC, box);
235 for (ftype = 0; (ftype < F_NRE); ftype++)
236 {
237 if ((interaction_function[ftype].flags & IF_CHEMBOND1<<3) == IF_CHEMBOND1<<3)
238 {
239 for (k = 0; (k < idef->il[ftype].nr); )
240 {
241 type = idef->il[ftype].iatoms[k++];
242 ai = idef->il[ftype].iatoms[k++];
243 aj = idef->il[ftype].iatoms[k++];
244 b0 = 0;
245 switch (ftype)
246 {
247 case F_BONDS:
248 b0 = idef->iparams[type].harmonic.rA;
249 break;
250 case F_G96BONDS:
251 b0 = sqrt(idef->iparams[type].harmonic.rA);
252 break;
253 case F_MORSE:
254 b0 = idef->iparams[type].morse.b0A;
255 break;
256 case F_CUBICBONDS:
257 b0 = idef->iparams[type].cubic.b0;
258 break;
259 case F_CONSTR:
260 b0 = idef->iparams[type].constr.dA;
261 break;
262 default:
263 break;
264 }
265 if (b0 != 0)
266 {
267 pbc_dx(&pbc, x[ai], x[aj], dx);
268 blen = norm(dx);
269 deviation = sqr(blen-b0);
270 if (sqrt(deviation/sqr(b0) > tol))
271 {
272 fprintf(stderrstderr, "Distance between atoms %d and %d is %.3f, should be %.3f\n", ai+1, aj+1, blen, b0);
273 }
274 }
275 }
276 }
277 }
278}
279
280void chk_trj(const output_env_t oenv, const char *fn, const char *tpr, real tol)
281{
282 t_trxframe fr;
283 t_count count;
284 t_fr_time first, last;
285 int j = -1, new_natoms, natoms;
286 real rdum, tt, old_t1, old_t2, prec;
287 gmx_bool bShowTimestep = TRUE1, bOK, newline = FALSE0;
288 t_trxstatus *status;
289 gmx_mtop_t mtop;
290 gmx_localtop_t *top = NULL((void*)0);
291 t_state state;
292 t_inputrec ir;
293
294 if (tpr)
295 {
296 read_tpx_state(tpr, &ir, &state, NULL((void*)0), &mtop);
297 top = gmx_mtop_generate_local_top(&mtop, &ir);
298 }
299 new_natoms = -1;
300 natoms = -1;
301
302 printf("Checking file %s\n", fn);
303
304 j = 0;
305 old_t2 = -2.0;
306 old_t1 = -1.0;
307
308 count.bStep = 0;
309 count.bTime = 0;
310 count.bLambda = 0;
311 count.bX = 0;
312 count.bV = 0;
313 count.bF = 0;
314 count.bBox = 0;
315
316 first.bStep = 0;
317 first.bTime = 0;
318 first.bLambda = 0;
319 first.bX = 0;
320 first.bV = 0;
321 first.bF = 0;
322 first.bBox = 0;
323
324 last.bStep = 0;
325 last.bTime = 0;
326 last.bLambda = 0;
327 last.bX = 0;
328 last.bV = 0;
329 last.bF = 0;
330 last.bBox = 0;
331
332 read_first_frame(oenv, &status, fn, &fr, TRX_READ_X(1<<0) | TRX_READ_V(1<<2) | TRX_READ_F(1<<4));
333
334 do
335 {
336 if (j == 0)
337 {
338 fprintf(stderrstderr, "\n# Atoms %d\n", fr.natoms);
339 if (fr.bPrec)
340 {
341 fprintf(stderrstderr, "Precision %g (nm)\n", 1/fr.prec);
342 }
343 }
344 newline = TRUE1;
345 if ((natoms > 0) && (new_natoms != natoms))
346 {
347 fprintf(stderrstderr, "\nNumber of atoms at t=%g don't match (%d, %d)\n",
348 old_t1, natoms, new_natoms);
349 newline = FALSE0;
350 }
351 if (j >= 2)
352 {
353 if (fabs((fr.time-old_t1)-(old_t1-old_t2)) >
354 0.1*(fabs(fr.time-old_t1)+fabs(old_t1-old_t2)) )
355 {
356 bShowTimestep = FALSE0;
357 fprintf(stderrstderr, "%sTimesteps at t=%g don't match (%g, %g)\n",
358 newline ? "\n" : "", old_t1, old_t1-old_t2, fr.time-old_t1);
359 }
360 }
361 natoms = new_natoms;
362 if (tpr)
363 {
364 chk_bonds(&top->idef, ir.ePBC, fr.x, fr.box, tol);
365 }
366 if (fr.bX)
367 {
368 chk_coords(j, natoms, fr.x, fr.box, 1e5, tol);
369 }
370 if (fr.bV)
371 {
372 chk_vels(j, natoms, fr.v);
373 }
374 if (fr.bF)
375 {
376 chk_forces(j, natoms, fr.f);
377 }
378
379 old_t2 = old_t1;
380 old_t1 = fr.time;
381 j++;
382 new_natoms = fr.natoms;
383#define INC(s, n, f, l, item) if (s.item != 0) { if (n.item == 0) { first.item = fr.time; } last.item = fr.time; n.item++; \
384}
385 INC(fr, count, first, last, bStep);
386 INC(fr, count, first, last, bTime);
387 INC(fr, count, first, last, bLambda);
388 INC(fr, count, first, last, bX);
389 INC(fr, count, first, last, bV);
390 INC(fr, count, first, last, bF);
391 INC(fr, count, first, last, bBox);
392#undef INC
393 }
394 while (read_next_frame(oenv, status, &fr));
395
396 fprintf(stderrstderr, "\n");
397
398 close_trj(status);
399
400 fprintf(stderrstderr, "\nItem #frames");
401 if (bShowTimestep)
402 {
403 fprintf(stderrstderr, " Timestep (ps)");
404 }
405 fprintf(stderrstderr, "\n");
406#define PRINTITEM(label, item)fprintf(stderr, "%-10s %6d", label, count.item); if ((bShowTimestep
) && (count.item > 1)) {fprintf(stderr, " %g\n"
, (last.item-first.item)/(count.item-1)); }else fprintf(stderr
, "\n") fprintf(stderrstderr, "%-10s %6d", label, count.item); if ((bShowTimestep) && (count.item > 1)) {fprintf(stderrstderr, " %g\n", (last.item-first.item)/(count.item-1)); }else fprintf(stderrstderr, "\n")
407 PRINTITEM ( "Step", bStep )fprintf(stderr, "%-10s %6d", "Step", count.bStep); if ((bShowTimestep
) && (count.bStep > 1)) {fprintf(stderr, " %g\n"
, (last.bStep-first.bStep)/(count.bStep-1)); }else fprintf(stderr
, "\n");
408 PRINTITEM ( "Time", bTime )fprintf(stderr, "%-10s %6d", "Time", count.bTime); if ((bShowTimestep
) && (count.bTime > 1)) {fprintf(stderr, " %g\n"
, (last.bTime-first.bTime)/(count.bTime-1)); }else fprintf(stderr
, "\n");
409 PRINTITEM ( "Lambda", bLambda )fprintf(stderr, "%-10s %6d", "Lambda", count.bLambda); if ((
bShowTimestep) && (count.bLambda > 1)) {fprintf(stderr
, " %g\n", (last.bLambda-first.bLambda)/(count.bLambda-1))
; }else fprintf(stderr, "\n");
410 PRINTITEM ( "Coords", bX )fprintf(stderr, "%-10s %6d", "Coords", count.bX); if ((bShowTimestep
) && (count.bX > 1)) {fprintf(stderr, " %g\n", (
last.bX-first.bX)/(count.bX-1)); }else fprintf(stderr, "\n");
411 PRINTITEM ( "Velocities", bV )fprintf(stderr, "%-10s %6d", "Velocities", count.bV); if ((bShowTimestep
) && (count.bV > 1)) {fprintf(stderr, " %g\n", (
last.bV-first.bV)/(count.bV-1)); }else fprintf(stderr, "\n");
412 PRINTITEM ( "Forces", bF )fprintf(stderr, "%-10s %6d", "Forces", count.bF); if ((bShowTimestep
) && (count.bF > 1)) {fprintf(stderr, " %g\n", (
last.bF-first.bF)/(count.bF-1)); }else fprintf(stderr, "\n");
413 PRINTITEM ( "Box", bBox )fprintf(stderr, "%-10s %6d", "Box", count.bBox); if ((bShowTimestep
) && (count.bBox > 1)) {fprintf(stderr, " %g\n"
, (last.bBox-first.bBox)/(count.bBox-1)); }else fprintf(stderr
, "\n");
414}
415
416void chk_tps(const char *fn, real vdw_fac, real bon_lo, real bon_hi)
417{
418 int natom, i, j, k;
419 char title[STRLEN4096];
420 t_topology top;
421 int ePBC;
422 t_atoms *atoms;
423 rvec *x, *v;
424 rvec dx;
425 matrix box;
426 t_pbc pbc;
427 gmx_bool bV, bX, bB, bFirst, bOut;
428 real r2, ekin, temp1, temp2, dist2, vdwfac2, bonlo2, bonhi2;
429 real *atom_vdw;
430 gmx_atomprop_t aps;
431
432 fprintf(stderrstderr, "Checking coordinate file %s\n", fn);
433 read_tps_conf(fn, title, &top, &ePBC, &x, &v, box, TRUE1);
434 atoms = &top.atoms;
435 natom = atoms->nr;
436 fprintf(stderrstderr, "%d atoms in file\n", atoms->nr);
437
438 /* check coordinates and box */
439 bV = FALSE0;
440 bX = FALSE0;
441 for (i = 0; (i < natom) && !(bV && bX); i++)
442 {
443 for (j = 0; (j < DIM3) && !(bV && bX); j++)
444 {
445 bV = bV || (v[i][j] != 0);
446 bX = bX || (x[i][j] != 0);
447 }
448 }
449 bB = FALSE0;
450 for (i = 0; (i < DIM3) && !bB; i++)
451 {
452 for (j = 0; (j < DIM3) && !bB; j++)
453 {
454 bB = bB || (box[i][j] != 0);
455 }
456 }
457
458 fprintf(stderrstderr, "coordinates %s\n", bX ? "found" : "absent");
459 fprintf(stderrstderr, "box %s\n", bB ? "found" : "absent");
460 fprintf(stderrstderr, "velocities %s\n", bV ? "found" : "absent");
461 fprintf(stderrstderr, "\n");
462
463 /* check velocities */
464 if (bV)
465 {
466 ekin = 0.0;
467 for (i = 0; (i < natom); i++)
468 {
469 for (j = 0; (j < DIM3); j++)
470 {
471 ekin += 0.5*atoms->atom[i].m*v[i][j]*v[i][j];
472 }
473 }
474 temp1 = (2.0*ekin)/(natom*DIM3*BOLTZ(((1.380658e-23)*(6.0221367e23))/(1e3)));
475 temp2 = (2.0*ekin)/(natom*(DIM3-1)*BOLTZ(((1.380658e-23)*(6.0221367e23))/(1e3)));
476 fprintf(stderrstderr, "Kinetic energy: %g (kJ/mol)\n", ekin);
477 fprintf(stderrstderr, "Assuming the number of degrees of freedom to be "
478 "Natoms * %d or Natoms * %d,\n"
479 "the velocities correspond to a temperature of the system\n"
480 "of %g K or %g K respectively.\n\n", DIM3, DIM3-1, temp1, temp2);
481 }
482
483 /* check coordinates */
484 if (bX)
485 {
486 vdwfac2 = sqr(vdw_fac);
487 bonlo2 = sqr(bon_lo);
488 bonhi2 = sqr(bon_hi);
489
490 fprintf(stderrstderr,
491 "Checking for atoms closer than %g and not between %g and %g,\n"
492 "relative to sum of Van der Waals distance:\n",
493 vdw_fac, bon_lo, bon_hi);
494 snew(atom_vdw, natom)(atom_vdw) = save_calloc("atom_vdw", "/home/alexxy/Develop/gromacs/src/gromacs/tools/check.c"
, 494, (natom), sizeof(*(atom_vdw)));
495 aps = gmx_atomprop_init();
496 for (i = 0; (i < natom); i++)
497 {
498 gmx_atomprop_query(aps, epropVDW,
499 *(atoms->resinfo[atoms->atom[i].resind].name),
500 *(atoms->atomname[i]), &(atom_vdw[i]));
501 if (debug)
502 {
503 fprintf(debug, "%5d %4s %4s %7g\n", i+1,
504 *(atoms->resinfo[atoms->atom[i].resind].name),
505 *(atoms->atomname[i]),
506 atom_vdw[i]);
507 }
508 }
509 gmx_atomprop_destroy(aps);
510 if (bB)
511 {
512 set_pbc(&pbc, ePBC, box);
513 }
514
515 bFirst = TRUE1;
516 for (i = 0; (i < natom); i++)
517 {
518 if (((i+1)%10) == 0)
519 {
520 fprintf(stderrstderr, "\r%5d", i+1);
521 }
522 for (j = i+1; (j < natom); j++)
523 {
524 if (bB)
525 {
526 pbc_dx(&pbc, x[i], x[j], dx);
527 }
528 else
529 {
530 rvec_sub(x[i], x[j], dx);
531 }
532 r2 = iprod(dx, dx);
533 dist2 = sqr(atom_vdw[i]+atom_vdw[j]);
534 if ( (r2 <= dist2*bonlo2) ||
535 ( (r2 >= dist2*bonhi2) && (r2 <= dist2*vdwfac2) ) )
536 {
537 if (bFirst)
538 {
539 fprintf(stderrstderr, "\r%5s %4s %8s %5s %5s %4s %8s %5s %6s\n",
540 "atom#", "name", "residue", "r_vdw",
541 "atom#", "name", "residue", "r_vdw", "distance");
542 bFirst = FALSE0;
543 }
544 fprintf(stderrstderr,
545 "\r%5d %4s %4s%4d %-5.3g %5d %4s %4s%4d %-5.3g %-6.4g\n",
546 i+1, *(atoms->atomname[i]),
547 *(atoms->resinfo[atoms->atom[i].resind].name),
548 atoms->resinfo[atoms->atom[i].resind].nr,
549 atom_vdw[i],
550 j+1, *(atoms->atomname[j]),
551 *(atoms->resinfo[atoms->atom[j].resind].name),
552 atoms->resinfo[atoms->atom[j].resind].nr,
553 atom_vdw[j],
554 sqrt(r2) );
555 }
556 }
557 }
558 if (bFirst)
559 {
560 fprintf(stderrstderr, "\rno close atoms found\n");
561 }
562 fprintf(stderrstderr, "\r \n");
563
564 if (bB)
565 {
566 /* check box */
567 bFirst = TRUE1;
568 k = 0;
569 for (i = 0; (i < natom) && (k < 10); i++)
570 {
571 bOut = FALSE0;
572 for (j = 0; (j < DIM3) && !bOut; j++)
573 {
574 bOut = bOut || (x[i][j] < 0) || (x[i][j] > box[j][j]);
575 }
576 if (bOut)
577 {
578 k++;
579 if (bFirst)
580 {
581 fprintf(stderrstderr, "Atoms outside box ( ");
582 for (j = 0; (j < DIM3); j++)
583 {
584 fprintf(stderrstderr, "%g ", box[j][j]);
585 }
586 fprintf(stderrstderr, "):\n"
587 "(These may occur often and are normally not a problem)\n"
588 "%5s %4s %8s %5s %s\n",
589 "atom#", "name", "residue", "r_vdw", "coordinate");
590 bFirst = FALSE0;
591 }
592 fprintf(stderrstderr,
593 "%5d %4s %4s%4d %-5.3g",
594 i, *(atoms->atomname[i]),
595 *(atoms->resinfo[atoms->atom[i].resind].name),
596 atoms->resinfo[atoms->atom[i].resind].nr, atom_vdw[i]);
597 for (j = 0; (j < DIM3); j++)
598 {
599 fprintf(stderrstderr, " %6.3g", x[i][j]);
600 }
601 fprintf(stderrstderr, "\n");
602 }
603 }
604 if (k == 10)
605 {
606 fprintf(stderrstderr, "(maybe more)\n");
607 }
608 if (bFirst)
609 {
610 fprintf(stderrstderr, "no atoms found outside box\n");
611 }
612 fprintf(stderrstderr, "\n");
613 }
614 }
615}
616
617void chk_ndx(const char *fn)
618{
619 t_blocka *grps;
620 char **grpname;
621 int i, j;
622
623 grps = init_index(fn, &grpname);
624 if (debug)
625 {
626 pr_blocka(debug, 0, fn, grps, FALSE0);
627 }
628 else
629 {
630 printf("Contents of index file %s\n", fn);
631 printf("--------------------------------------------------\n");
632 printf("Nr. Group #Entries First Last\n");
633 for (i = 0; (i < grps->nr); i++)
634 {
635 printf("%4d %-20s%8d%8d%8d\n", i, grpname[i],
636 grps->index[i+1]-grps->index[i],
637 grps->a[grps->index[i]]+1,
638 grps->a[grps->index[i+1]-1]+1);
639 }
640 }
641 for (i = 0; (i < grps->nr); i++)
642 {
643 sfree(grpname[i])save_free("grpname[i]", "/home/alexxy/Develop/gromacs/src/gromacs/tools/check.c"
, 643, (grpname[i]));
644 }
645 sfree(grpname)save_free("grpname", "/home/alexxy/Develop/gromacs/src/gromacs/tools/check.c"
, 645, (grpname));
646 done_blocka(grps);
647}
648
649void chk_enx(const char *fn)
650{
651 int nre, fnr, ndr;
652 ener_file_t in;
653 gmx_enxnm_t *enm = NULL((void*)0);
654 t_enxframe *fr;
655 gmx_bool bShowTStep;
656 real t0, old_t1, old_t2;
657 char buf[22];
658
659 fprintf(stderrstderr, "Checking energy file %s\n\n", fn);
660
661 in = open_enx(fn, "r");
662 do_enxnms(in, &nre, &enm);
663 fprintf(stderrstderr, "%d groups in energy file", nre);
664 snew(fr, 1)(fr) = save_calloc("fr", "/home/alexxy/Develop/gromacs/src/gromacs/tools/check.c"
, 664, (1), sizeof(*(fr)));
665 old_t2 = -2.0;
666 old_t1 = -1.0;
667 fnr = 0;
668 t0 = NOTSET-12345;
669 bShowTStep = TRUE1;
670
671 while (do_enx(in, fr))
672 {
673 if (fnr >= 2)
674 {
675 if (fabs((fr->t-old_t1)-(old_t1-old_t2)) >
676 0.1*(fabs(fr->t-old_t1)+fabs(old_t1-old_t2)) )
677 {
678 bShowTStep = FALSE0;
679 fprintf(stderrstderr, "\nTimesteps at t=%g don't match (%g, %g)\n",
680 old_t1, old_t1-old_t2, fr->t-old_t1);
681 }
682 }
683 old_t2 = old_t1;
684 old_t1 = fr->t;
685 if (t0 == NOTSET-12345)
686 {
687 t0 = fr->t;
688 }
689 if (fnr == 0)
690 {
691 fprintf(stderrstderr, "\rframe: %6s (index %6d), t: %10.3f\n",
692 gmx_step_str(fr->step, buf), fnr, fr->t);
693 }
694 fnr++;
695 }
696 fprintf(stderrstderr, "\n\nFound %d frames", fnr);
697 if (bShowTStep && fnr > 1)
698 {
699 fprintf(stderrstderr, " with a timestep of %g ps", (old_t1-t0)/(fnr-1));
700 }
701 fprintf(stderrstderr, ".\n");
702
703 free_enxframe(fr);
704 free_enxnms(nre, enm);
705 sfree(fr)save_free("fr", "/home/alexxy/Develop/gromacs/src/gromacs/tools/check.c"
, 705, (fr));
706}
707
708int gmx_check(int argc, char *argv[])
709{
710 const char *desc[] = {
711 "[THISMODULE] reads a trajectory ([TT].trj[tt], [TT].trr[tt] or ",
712 "[TT].xtc[tt]), an energy file ([TT].ene[tt] or [TT].edr[tt])",
713 "or an index file ([TT].ndx[tt])",
714 "and prints out useful information about them.[PAR]",
715 "Option [TT]-c[tt] checks for presence of coordinates,",
716 "velocities and box in the file, for close contacts (smaller than",
717 "[TT]-vdwfac[tt] and not bonded, i.e. not between [TT]-bonlo[tt]",
718 "and [TT]-bonhi[tt], all relative to the sum of both Van der Waals",
719 "radii) and atoms outside the box (these may occur often and are",
720 "no problem). If velocities are present, an estimated temperature",
721 "will be calculated from them.[PAR]",
722 "If an index file, is given its contents will be summarized.[PAR]",
723 "If both a trajectory and a [TT].tpr[tt] file are given (with [TT]-s1[tt])",
724 "the program will check whether the bond lengths defined in the tpr",
725 "file are indeed correct in the trajectory. If not you may have",
726 "non-matching files due to e.g. deshuffling or due to problems with",
727 "virtual sites. With these flags, [TT]gmx check[tt] provides a quick check for such problems.[PAR]",
728 "The program can compare two run input ([TT].tpr[tt], [TT].tpb[tt] or",
729 "[TT].tpa[tt]) files",
730 "when both [TT]-s1[tt] and [TT]-s2[tt] are supplied.",
731 "Similarly a pair of trajectory files can be compared (using the [TT]-f2[tt]",
732 "option), or a pair of energy files (using the [TT]-e2[tt] option).[PAR]",
733 "For free energy simulations the A and B state topology from one",
734 "run input file can be compared with options [TT]-s1[tt] and [TT]-ab[tt].[PAR]",
735 "In case the [TT]-m[tt] flag is given a LaTeX file will be written",
736 "consisting of a rough outline for a methods section for a paper."
737 };
738 t_filenm fnm[] = {
739 { efTRX, "-f", NULL((void*)0), ffOPTRD(1<<1 | 1<<3) },
740 { efTRX, "-f2", NULL((void*)0), ffOPTRD(1<<1 | 1<<3) },
741 { efTPX, "-s1", "top1", ffOPTRD(1<<1 | 1<<3) },
742 { efTPX, "-s2", "top2", ffOPTRD(1<<1 | 1<<3) },
743 { efTPS, "-c", NULL((void*)0), ffOPTRD(1<<1 | 1<<3) },
744 { efEDR, "-e", NULL((void*)0), ffOPTRD(1<<1 | 1<<3) },
745 { efEDR, "-e2", "ener2", ffOPTRD(1<<1 | 1<<3) },
746 { efNDX, "-n", NULL((void*)0), ffOPTRD(1<<1 | 1<<3) },
747 { efTEX, "-m", NULL((void*)0), ffOPTWR(1<<2| 1<<3) }
748 };
749#define NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))) asize(fnm)((int)(sizeof(fnm)/sizeof((fnm)[0])))
750 const char *fn1 = NULL((void*)0), *fn2 = NULL((void*)0), *tex = NULL((void*)0);
751
752 output_env_t oenv;
753 static real vdw_fac = 0.8;
754 static real bon_lo = 0.4;
755 static real bon_hi = 0.7;
756 static gmx_bool bRMSD = FALSE0;
757 static real ftol = 0.001;
758 static real abstol = 0.001;
759 static gmx_bool bCompAB = FALSE0;
760 static char *lastener = NULL((void*)0);
761 static t_pargs pa[] = {
762 { "-vdwfac", FALSE0, etREAL, {&vdw_fac},
763 "Fraction of sum of VdW radii used as warning cutoff" },
764 { "-bonlo", FALSE0, etREAL, {&bon_lo},
765 "Min. fract. of sum of VdW radii for bonded atoms" },
766 { "-bonhi", FALSE0, etREAL, {&bon_hi},
767 "Max. fract. of sum of VdW radii for bonded atoms" },
768 { "-rmsd", FALSE0, etBOOL, {&bRMSD},
769 "Print RMSD for x, v and f" },
770 { "-tol", FALSE0, etREAL, {&ftol},
771 "Relative tolerance for comparing real values defined as [MATH]2*(a-b)/([MAG]a[mag]+[MAG]b[mag])[math]" },
772 { "-abstol", FALSE0, etREAL, {&abstol},
773 "Absolute tolerance, useful when sums are close to zero." },
774 { "-ab", FALSE0, etBOOL, {&bCompAB},
775 "Compare the A and B topology from one file" },
776 { "-lastener", FALSE0, etSTR, {&lastener},
777 "Last energy term to compare (if not given all are tested). It makes sense to go up until the Pressure." }
778 };
779
780 if (!parse_common_args(&argc, argv, 0, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm, asize(pa)((int)(sizeof(pa)/sizeof((pa)[0]))), pa,
781 asize(desc)((int)(sizeof(desc)/sizeof((desc)[0]))), desc, 0, NULL((void*)0), &oenv))
782 {
783 return 0;
784 }
785
786 fn1 = opt2fn_null("-f", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm);
787 fn2 = opt2fn_null("-f2", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm);
788 tex = opt2fn_null("-m", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm);
789 if (fn1 && fn2)
790 {
791 comp_trx(oenv, fn1, fn2, bRMSD, ftol, abstol);
792 }
793 else if (fn1)
794 {
795 chk_trj(oenv, fn1, opt2fn_null("-s1", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), ftol);
796 }
797 else if (fn2)
798 {
799 fprintf(stderrstderr, "Please give me TWO trajectory (.xtc/.trr/.trj) files!\n");
800 }
801
802 fn1 = opt2fn_null("-s1", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm);
803 fn2 = opt2fn_null("-s2", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm);
804 if ((fn1 && fn2) || bCompAB)
805 {
806 if (bCompAB)
807 {
808 if (fn1 == NULL((void*)0))
809 {
810 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/tools/check.c", 810, "With -ab you need to set the -s1 option");
811 }
812 fn2 = NULL((void*)0);
813 }
814 comp_tpx(fn1, fn2, bRMSD, ftol, abstol);
815 }
816 else if (fn1 && tex)
817 {
818 tpx2methods(fn1, tex);
819 }
820 else if ((fn1 && !opt2fn_null("-f", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm)) || (!fn1 && fn2))
821 {
822 fprintf(stderrstderr, "Please give me TWO run input (.tpr/.tpa/.tpb) files\n"
823 "or specify the -m flag to generate a methods.tex file\n");
824 }
825
826 fn1 = opt2fn_null("-e", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm);
827 fn2 = opt2fn_null("-e2", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm);
828 if (fn1 && fn2)
829 {
830 comp_enx(fn1, fn2, ftol, abstol, lastener);
831 }
832 else if (fn1)
833 {
834 chk_enx(ftp2fn(efEDR, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm));
835 }
836 else if (fn2)
837 {
838 fprintf(stderrstderr, "Please give me TWO energy (.edr/.ene) files!\n");
839 }
840
841 if (ftp2bSet(efTPS, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm))
842 {
843 chk_tps(ftp2fn(efTPS, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), vdw_fac, bon_lo, bon_hi);
844 }
845
846 if (ftp2bSet(efNDX, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm))
847 {
848 chk_ndx(ftp2fn(efNDX, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm));
849 }
850
851 return 0;
852}