src/gromacs/tools/check.cpp

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