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 <cmath>
  40 #include <cstdio>
  41 #include <cstring>
  42
  43 #include "gromacs/commandline/pargs.h"
  44 #include "gromacs/fileio/confio.h"
  45 #include "gromacs/fileio/enxio.h"
  46 #include "gromacs/fileio/gmxfio.h"
  47 #include "gromacs/fileio/tpxio.h"
  48 #include "gromacs/fileio/trxio.h"
  49 #include "gromacs/fileio/xtcio.h"
  50 #include "gromacs/math/functions.h"
  51 #include "gromacs/math/units.h"
  52 #include "gromacs/math/vec.h"
  53 #include "gromacs/mdrunutility/mdmodules.h"
  54 #include "gromacs/mdtypes/inputrec.h"
  55 #include "gromacs/mdtypes/md_enums.h"
  56 #include "gromacs/pbcutil/pbc.h"
  57 #include "gromacs/topology/atomprop.h"
  58 #include "gromacs/topology/block.h"
  59 #include "gromacs/topology/ifunc.h"
  60 #include "gromacs/topology/index.h"
  61 #include "gromacs/topology/mtop_util.h"
  62 #include "gromacs/topology/topology.h"
  63 #include "gromacs/trajectory/trajectoryframe.h"
  64 #include "gromacs/utility/arraysize.h"
  65 #include "gromacs/utility/fatalerror.h"
  66 #include "gromacs/utility/futil.h"
  67 #include "gromacs/utility/smalloc.h"
  68
  69 typedef struct {
  70     int bStep;
  71     int bTime;
  72     int bLambda;
  73     int bX;
  74     int bV;
  75     int bF;
  76     int bBox;
  77 } t_count;
  78
  79 typedef struct {
  80     float bStep;
  81     float bTime;
  82     float bLambda;
  83     float bX;
  84     float bV;
  85     float bF;
  86     float bBox;
  87 } t_fr_time;
  88
  89 static void comp_tpx(const char *fn1, const char *fn2,
  90                      gmx_bool bRMSD, real ftol, real abstol)
  91 {
  92     const char    *ff[2];
  93     gmx::MDModules mdModules[2];
  94     t_inputrec    *ir[2] = { mdModules[0].inputrec(), mdModules[1].inputrec() };
  95     t_state        state[2];
  96     gmx_mtop_t     mtop[2];
  97     t_topology     top[2];
  98     int            i;
  99
 100     ff[0] = fn1;
 101     ff[1] = fn2;
 102     for (i = 0; i < (fn2 ? 2 : 1); i++)
 103     {
 104         read_tpx_state(ff[i], ir[i], &state[i], &(mtop[i]));
 105     }
 106     if (fn2)
 107     {
 108         cmp_inputrec(stdout, ir[0], ir[1], ftol, abstol);
 109         /* Convert gmx_mtop_t to t_topology.
 110          * We should implement direct mtop comparison,
 111          * but it might be useful to keep t_topology comparison as an option.
 112          */
 113         top[0] = gmx_mtop_t_to_t_topology(&mtop[0], false);
 114         top[1] = gmx_mtop_t_to_t_topology(&mtop[1], false);
 115         cmp_top(stdout, &top[0], &top[1], ftol, abstol);
 116         cmp_groups(stdout, &mtop[0].groups, &mtop[1].groups,
 117                    mtop[0].natoms, mtop[1].natoms);
 118         comp_state(&state[0], &state[1], bRMSD, ftol, abstol);
 119     }
 120     else
 121     {
 122         if (ir[0]->efep == efepNO)
 123         {
 124             fprintf(stdout, "inputrec->efep = %s\n", efep_names[ir[0]->efep]);
 125         }
 126         else
 127         {
 128             if (ir[0]->bPull)
 129             {
 130                 comp_pull_AB(stdout, ir[0]->pull, ftol, abstol);
 131             }
 132             /* Convert gmx_mtop_t to t_topology.
 133              * We should implement direct mtop comparison,
 134              * but it might be useful to keep t_topology comparison as an option.
 135              */
 136             top[0] = gmx_mtop_t_to_t_topology(&mtop[0], true);
 137             cmp_top(stdout, &top[0], nullptr, ftol, abstol);
 138         }
 139     }
 140 }
 141
 142 static void comp_trx(const gmx_output_env_t *oenv, const char *fn1, const char *fn2,
 143                      gmx_bool bRMSD, real ftol, real abstol)
 144 {
 145     int          i;
 146     const char  *fn[2];
 147     t_trxframe   fr[2];
 148     t_trxstatus *status[2];
 149     gmx_bool     b[2];
 150
 151     fn[0] = fn1;
 152     fn[1] = fn2;
 153     fprintf(stderr, "Comparing trajectory files %s and %s\n", fn1, fn2);
 154     for (i = 0; i < 2; i++)
 155     {
 156         b[i] = read_first_frame(oenv, &status[i], fn[i], &fr[i], TRX_READ_X|TRX_READ_V|TRX_READ_F);
 157     }
 158
 159     if (b[0] && b[1])
 160     {
 161         do
 162         {
 163             comp_frame(stdout, &(fr[0]), &(fr[1]), bRMSD, ftol, abstol);
 164
 165             for (i = 0; i < 2; i++)
 166             {
 167                 b[i] = read_next_frame(oenv, status[i], &fr[i]);
 168             }
 169         }
 170         while (b[0] && b[1]);
 171
 172         for (i = 0; i < 2; i++)
 173         {
 174             if (b[i] && !b[1-i])
 175             {
 176                 fprintf(stdout, "\nEnd of file on %s but not on %s\n", fn[1-i], fn[i]);
 177             }
 178             close_trj(status[i]);
 179         }
 180     }
 181     if (!b[0] && !b[1])
 182     {
 183         fprintf(stdout, "\nBoth files read correctly\n");
 184     }
 185 }
 186
 187 static void tpx2system(FILE *fp, const gmx_mtop_t *mtop)
 188 {
 189     int                       nmol, nvsite = 0;
 190     gmx_mtop_atomloop_block_t aloop;
 191     const t_atom             *atom;
 192
 193     fprintf(fp, "\\subsection{Simulation system}\n");
 194     aloop = gmx_mtop_atomloop_block_init(mtop);
 195     while (gmx_mtop_atomloop_block_next(aloop, &atom, &nmol))
 196     {
 197         if (atom->ptype == eptVSite)
 198         {
 199             nvsite += nmol;
 200         }
 201     }
 202     fprintf(fp, "A system of %d molecules (%d atoms) was simulated.\n",
 203             mtop->mols.nr, mtop->natoms-nvsite);
 204     if (nvsite)
 205     {
 206         fprintf(fp, "Virtual sites were used in some of the molecules.\n");
 207     }
 208     fprintf(fp, "\n\n");
 209 }
 210
 211 static void tpx2params(FILE *fp, const t_inputrec *ir)
 212 {
 213     fprintf(fp, "\\subsection{Simulation settings}\n");
 214     fprintf(fp, "A total of %g ns were simulated with a time step of %g fs.\n",
 215             ir->nsteps*ir->delta_t*0.001, 1000*ir->delta_t);
 216     fprintf(fp, "Neighbor searching was performed every %d steps.\n", ir->nstlist);
 217     fprintf(fp, "The %s algorithm was used for electrostatic interactions.\n",
 218             EELTYPE(ir->coulombtype));
 219     fprintf(fp, "with a cut-off of %g nm.\n", ir->rcoulomb);
 220     if (ir->coulombtype == eelPME)
 221     {
 222         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);
 223     }
 224     if (ir->rvdw > ir->rlist)
 225     {
 226         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);
 227     }
 228     else
 229     {
 230         fprintf(fp, "A single cut-off of %g was used for Van der Waals interactions.\n", ir->rlist);
 231     }
 232     if (ir->etc != 0)
 233     {
 234         fprintf(fp, "Temperature coupling was done with the %s algorithm.\n",
 235                 etcoupl_names[ir->etc]);
 236     }
 237     if (ir->epc != 0)
 238     {
 239         fprintf(fp, "Pressure coupling was done with the %s algorithm.\n",
 240                 epcoupl_names[ir->epc]);
 241     }
 242     fprintf(fp, "\n\n");
 243 }
 244
 245 static void tpx2methods(const char *tpx, const char *tex)
 246 {
 247     FILE          *fp;
 248     t_inputrec    *ir;
 249     t_state        state;
 250     gmx_mtop_t     mtop;
 251
 252     gmx::MDModules mdModules;
 253     ir = mdModules.inputrec();
 254     read_tpx_state(tpx, ir, &state, &mtop);
 255     fp = gmx_fio_fopen(tex, "w");
 256     fprintf(fp, "\\section{Methods}\n");
 257     tpx2system(fp, &mtop);
 258     tpx2params(fp, ir);
 259     gmx_fio_fclose(fp);
 260 }
 261
 262 static void chk_coords(int frame, int natoms, rvec *x, matrix box, real fac, real tol)
 263 {
 264     int  i, j;
 265     int  nNul = 0;
 266     real vol  = det(box);
 267
 268     for (i = 0; (i < natoms); i++)
 269     {
 270         for (j = 0; (j < DIM); j++)
 271         {
 272             if ((vol > 0) && (fabs(x[i][j]) > fac*box[j][j]))
 273             {
 274                 printf("Warning at frame %d: coordinates for atom %d are large (%g)\n",
 275                        frame, i, x[i][j]);
 276             }
 277         }
 278         if ((fabs(x[i][XX]) < tol) &&
 279             (fabs(x[i][YY]) < tol) &&
 280             (fabs(x[i][ZZ]) < tol))
 281         {
 282             nNul++;
 283         }
 284     }
 285     if (nNul > 0)
 286     {
 287         printf("Warning at frame %d: there are %d particles with all coordinates zero\n",
 288                frame, nNul);
 289     }
 290 }
 291
 292 static void chk_vels(int frame, int natoms, rvec *v)
 293 {
 294     int i, j;
 295
 296     for (i = 0; (i < natoms); i++)
 297     {
 298         for (j = 0; (j < DIM); j++)
 299         {
 300             if (fabs(v[i][j]) > 500)
 301             {
 302                 printf("Warning at frame %d. Velocities for atom %d are large (%g)\n",
 303                        frame, i, v[i][j]);
 304             }
 305         }
 306     }
 307 }
 308
 309 static void chk_forces(int frame, int natoms, rvec *f)
 310 {
 311     int i, j;
 312
 313     for (i = 0; (i < natoms); i++)
 314     {
 315         for (j = 0; (j < DIM); j++)
 316         {
 317             if (fabs(f[i][j]) > 10000)
 318             {
 319                 printf("Warning at frame %d. Forces for atom %d are large (%g)\n",
 320                        frame, i, f[i][j]);
 321             }
 322         }
 323     }
 324 }
 325
 326 static void chk_bonds(t_idef *idef, int ePBC, rvec *x, matrix box, real tol)
 327 {
 328     int   ftype, k, ai, aj, type;
 329     real  b0, blen, deviation;
 330     t_pbc pbc;
 331     rvec  dx;
 332
 333     set_pbc(&pbc, ePBC, box);
 334     for (ftype = 0; (ftype < F_NRE); ftype++)
 335     {
 336         if ((interaction_function[ftype].flags & IF_CHEMBOND) == IF_CHEMBOND)
 337         {
 338             for (k = 0; (k < idef->il[ftype].nr); )
 339             {
 340                 type = idef->il[ftype].iatoms[k++];
 341                 ai   = idef->il[ftype].iatoms[k++];
 342                 aj   = idef->il[ftype].iatoms[k++];
 343                 b0   = 0;
 344                 switch (ftype)
 345                 {
 346                     case F_BONDS:
 347                         b0 = idef->iparams[type].harmonic.rA;
 348                         break;
 349                     case F_G96BONDS:
 350                         b0 = std::sqrt(idef->iparams[type].harmonic.rA);
 351                         break;
 352                     case F_MORSE:
 353                         b0 = idef->iparams[type].morse.b0A;
 354                         break;
 355                     case F_CUBICBONDS:
 356                         b0 = idef->iparams[type].cubic.b0;
 357                         break;
 358                     case F_CONSTR:
 359                         b0 = idef->iparams[type].constr.dA;
 360                         break;
 361                     default:
 362                         break;
 363                 }
 364                 if (b0 != 0)
 365                 {
 366                     pbc_dx(&pbc, x[ai], x[aj], dx);
 367                     blen      = norm(dx);
 368                     deviation = gmx::square(blen-b0);
 369                     if (std::sqrt(deviation/gmx::square(b0)) > tol)
 370                     {
 371                         fprintf(stderr, "Distance between atoms %d and %d is %.3f, should be %.3f\n", ai+1, aj+1, blen, b0);
 372                     }
 373                 }
 374             }
 375         }
 376     }
 377 }
 378
 379 void chk_trj(const gmx_output_env_t *oenv, const char *fn, const char *tpr, real tol)
 380 {
 381     t_trxframe       fr;
 382     t_count          count;
 383     t_fr_time        first, last;
 384     int              j = -1, new_natoms, natoms;
 385     real             old_t1, old_t2;
 386     gmx_bool         bShowTimestep = TRUE, newline = FALSE;
 387     t_trxstatus     *status;
 388     gmx_mtop_t       mtop;
 389     gmx_localtop_t  *top = nullptr;
 390     t_state          state;
 391     t_inputrec      *ir;
 392
 393     gmx::MDModules   mdModules;
 394     ir = mdModules.inputrec();
 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_trj(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 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 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 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.",
 835         "Similarly a pair of trajectory files can be compared (using the [TT]-f2[tt]",
 836         "option), or a pair of energy files (using the [TT]-e2[tt] option).[PAR]",
 837         "For free energy simulations the A and B state topology from one",
 838         "run input file can be compared with options [TT]-s1[tt] and [TT]-ab[tt].[PAR]",
 839         "In case the [TT]-m[tt] flag is given a LaTeX file will be written",
 840         "consisting of a rough outline for a methods section for a paper."
 841     };
 842     t_filenm          fnm[] = {
 843         { efTRX, "-f",  nullptr, ffOPTRD },
 844         { efTRX, "-f2",  nullptr, ffOPTRD },
 845         { efTPR, "-s1", "top1", ffOPTRD },
 846         { efTPR, "-s2", "top2", ffOPTRD },
 847         { efTPS, "-c",  nullptr, ffOPTRD },
 848         { efEDR, "-e",  nullptr, ffOPTRD },
 849         { efEDR, "-e2", "ener2", ffOPTRD },
 850         { efNDX, "-n",  nullptr, ffOPTRD },
 851         { efTEX, "-m",  nullptr, ffOPTWR }
 852     };
 853 #define NFILE asize(fnm)
 854     const char       *fn1 = nullptr, *fn2 = nullptr, *tex = nullptr;
 855
 856     gmx_output_env_t *oenv;
 857     static real       vdw_fac  = 0.8;
 858     static real       bon_lo   = 0.4;
 859     static real       bon_hi   = 0.7;
 860     static gmx_bool   bRMSD    = FALSE;
 861     static real       ftol     = 0.001;
 862     static real       abstol   = 0.001;
 863     static gmx_bool   bCompAB  = FALSE;
 864     static char      *lastener = nullptr;
 865     static t_pargs    pa[]     = {
 866         { "-vdwfac", FALSE, etREAL, {&vdw_fac},
 867           "Fraction of sum of VdW radii used as warning cutoff" },
 868         { "-bonlo",  FALSE, etREAL, {&bon_lo},
 869           "Min. fract. of sum of VdW radii for bonded atoms" },
 870         { "-bonhi",  FALSE, etREAL, {&bon_hi},
 871           "Max. fract. of sum of VdW radii for bonded atoms" },
 872         { "-rmsd",   FALSE, etBOOL, {&bRMSD},
 873           "Print RMSD for x, v and f" },
 874         { "-tol",    FALSE, etREAL, {&ftol},
 875           "Relative tolerance for comparing real values defined as [MATH]2*(a-b)/([MAG]a[mag]+[MAG]b[mag])[math]" },
 876         { "-abstol",    FALSE, etREAL, {&abstol},
 877           "Absolute tolerance, useful when sums are close to zero." },
 878         { "-ab",     FALSE, etBOOL, {&bCompAB},
 879           "Compare the A and B topology from one file" },
 880         { "-lastener", FALSE, etSTR,  {&lastener},
 881           "Last energy term to compare (if not given all are tested). It makes sense to go up until the Pressure." }
 882     };
 883
 884     if (!parse_common_args(&argc, argv, 0, NFILE, fnm, asize(pa), pa,
 885                            asize(desc), desc, 0, nullptr, &oenv))
 886     {
 887         return 0;
 888     }
 889
 890     fn1 = opt2fn_null("-f", NFILE, fnm);
 891     fn2 = opt2fn_null("-f2", NFILE, fnm);
 892     tex = opt2fn_null("-m", NFILE, fnm);
 893     if (fn1 && fn2)
 894     {
 895         comp_trx(oenv, fn1, fn2, bRMSD, ftol, abstol);
 896     }
 897     else if (fn1)
 898     {
 899         chk_trj(oenv, fn1, opt2fn_null("-s1", NFILE, fnm), ftol);
 900     }
 901     else if (fn2)
 902     {
 903         fprintf(stderr, "Please give me TWO trajectory (.xtc/.trr/.tng) files!\n");
 904     }
 905
 906     fn1 = opt2fn_null("-s1", NFILE, fnm);
 907     fn2 = opt2fn_null("-s2", NFILE, fnm);
 908     if ((fn1 && fn2) || bCompAB)
 909     {
 910         if (bCompAB)
 911         {
 912             if (fn1 == nullptr)
 913             {
 914                 gmx_fatal(FARGS, "With -ab you need to set the -s1 option");
 915             }
 916             fn2 = nullptr;
 917         }
 918         comp_tpx(fn1, fn2, bRMSD, ftol, abstol);
 919     }
 920     else if (fn1 && tex)
 921     {
 922         tpx2methods(fn1, tex);
 923     }
 924     else if ((fn1 && !opt2fn_null("-f", NFILE, fnm)) || (!fn1 && fn2))
 925     {
 926         fprintf(stderr, "Please give me TWO run input (.tpr) files\n"
 927                 "or specify the -m flag to generate a methods.tex file\n");
 928     }
 929
 930     fn1 = opt2fn_null("-e", NFILE, fnm);
 931     fn2 = opt2fn_null("-e2", NFILE, fnm);
 932     if (fn1 && fn2)
 933     {
 934         comp_enx(fn1, fn2, ftol, abstol, lastener);
 935     }
 936     else if (fn1)
 937     {
 938         chk_enx(ftp2fn(efEDR, NFILE, fnm));
 939     }
 940     else if (fn2)
 941     {
 942         fprintf(stderr, "Please give me TWO energy (.edr/.ene) files!\n");
 943     }
 944
 945     if (ftp2bSet(efTPS, NFILE, fnm))
 946     {
 947         chk_tps(ftp2fn(efTPS, NFILE, fnm), vdw_fac, bon_lo, bon_hi);
 948     }
 949
 950     if (ftp2bSet(efNDX, NFILE, fnm))
 951     {
 952         chk_ndx(ftp2fn(efNDX, NFILE, fnm));
 953     }
 954
 955     return 0;
 956 }