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