src/gromacs/topology/topology.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-2004, 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 "topology.h"
  40
  41 #include <cstdio>
  42
  43 #include <algorithm>
  44
  45 #include "gromacs/math/vecdump.h"
  46 #include "gromacs/topology/atoms.h"
  47 #include "gromacs/topology/idef.h"
  48 #include "gromacs/topology/ifunc.h"
  49 #include "gromacs/topology/symtab.h"
  50 #include "gromacs/utility/compare.h"
  51 #include "gromacs/utility/gmxassert.h"
  52 #include "gromacs/utility/smalloc.h"
  53 #include "gromacs/utility/strconvert.h"
  54 #include "gromacs/utility/txtdump.h"
  55
  56 const char *gtypes[egcNR+1] = {
  57     "T-Coupling", "Energy Mon.", "Acceleration", "Freeze",
  58     "User1", "User2", "VCM", "Compressed X", "Or. Res. Fit", "QMMM", nullptr
  59 };
  60
  61 static void init_groups(gmx_groups_t *groups)
  62 {
  63     groups->ngrpname = 0;
  64     groups->grpname  = nullptr;
  65     for (int g = 0; g < egcNR; g++)
  66     {
  67         groups->grps[g].nr     = 0;
  68         groups->grps[g].nm_ind = nullptr;
  69         groups->ngrpnr[g]      = 0;
  70         groups->grpnr[g]       = nullptr;
  71     }
  72
  73 }
  74
  75 void init_mtop(gmx_mtop_t *mtop)
  76 {
  77     mtop->name = nullptr;
  78
  79     // TODO: Move to ffparams when that is converted to C++
  80     mtop->ffparams.functype.clear();
  81     mtop->ffparams.iparams.clear();
  82     mtop->ffparams.cmap_grid.grid_spacing = 0;
  83     mtop->ffparams.cmap_grid.cmapdata.clear();
  84
  85     mtop->moltype.clear();
  86     mtop->molblock.clear();
  87     mtop->bIntermolecularInteractions = FALSE;
  88     mtop->intermolecular_ilist        = nullptr;
  89
  90     mtop->natoms       = 0;
  91     mtop->maxres_renum = 0;
  92     mtop->maxresnr     = -1;
  93     init_atomtypes(&mtop->atomtypes);
  94     init_groups(&mtop->groups);
  95     open_symtab(&mtop->symtab);
  96 }
  97
  98 void init_top(t_topology *top)
  99 {
 100     top->name = nullptr;
 101     init_idef(&top->idef);
 102     init_atom(&(top->atoms));
 103     init_atomtypes(&(top->atomtypes));
 104     init_block(&top->cgs);
 105     init_block(&top->mols);
 106     init_blocka(&top->excls);
 107     open_symtab(&top->symtab);
 108 }
 109
 110
 111 gmx_moltype_t::gmx_moltype_t() :
 112     name(nullptr),
 113     cgs(),
 114     excls()
 115 {
 116     init_t_atoms(&atoms, 0, FALSE);
 117 }
 118
 119 gmx_moltype_t::~gmx_moltype_t()
 120 {
 121     done_atom(&atoms);
 122     done_block(&cgs);
 123     done_blocka(&excls);
 124 }
 125
 126 void done_gmx_groups_t(gmx_groups_t *g)
 127 {
 128     int i;
 129
 130     for (i = 0; (i < egcNR); i++)
 131     {
 132         if (nullptr != g->grps[i].nm_ind)
 133         {
 134             sfree(g->grps[i].nm_ind);
 135             g->grps[i].nm_ind = nullptr;
 136         }
 137         if (nullptr != g->grpnr[i])
 138         {
 139             sfree(g->grpnr[i]);
 140             g->grpnr[i] = nullptr;
 141         }
 142     }
 143     /* The contents of this array is in symtab, don't free it here */
 144     sfree(g->grpname);
 145 }
 146
 147 gmx_mtop_t::gmx_mtop_t()
 148 {
 149     init_mtop(this);
 150 }
 151
 152 gmx_mtop_t::~gmx_mtop_t()
 153 {
 154     done_symtab(&symtab);
 155
 156     moltype.clear();
 157     molblock.clear();
 158     done_atomtypes(&atomtypes);
 159     done_gmx_groups_t(&groups);
 160 }
 161
 162 void done_top(t_topology *top)
 163 {
 164     done_idef(&top->idef);
 165     done_atom(&(top->atoms));
 166
 167     /* For GB */
 168     done_atomtypes(&(top->atomtypes));
 169
 170     done_symtab(&(top->symtab));
 171     done_block(&(top->cgs));
 172     done_block(&(top->mols));
 173     done_blocka(&(top->excls));
 174 }
 175
 176 void done_top_mtop(t_topology *top, gmx_mtop_t *mtop)
 177 {
 178     if (mtop != nullptr)
 179     {
 180         if (top != nullptr)
 181         {
 182             done_idef(&top->idef);
 183             done_atom(&top->atoms);
 184             done_block(&top->cgs);
 185             done_blocka(&top->excls);
 186             done_block(&top->mols);
 187             done_symtab(&top->symtab);
 188             open_symtab(&mtop->symtab);
 189             done_atomtypes(&(top->atomtypes));
 190         }
 191
 192         // Note that the rest of mtop will be freed by the destructor
 193     }
 194 }
 195
 196 void init_localtop(gmx_localtop_t *top)
 197 {
 198     init_block(&top->cgs);
 199     init_blocka(&top->excls);
 200     init_idef(&top->idef);
 201     init_atomtypes(&top->atomtypes);
 202 }
 203
 204 void done_localtop(gmx_localtop_t *top)
 205 {
 206     if (top == nullptr)
 207     {
 208         return;
 209     }
 210     done_idef(&top->idef);
 211     done_block(&top->cgs);
 212     done_blocka(&top->excls);
 213     done_atomtypes(&top->atomtypes);
 214 }
 215
 216 void done_and_sfree_localtop(gmx_localtop_t *top)
 217 {
 218     done_localtop(top);
 219     sfree(top);
 220 }
 221
 222 bool gmx_mtop_has_masses(const gmx_mtop_t *mtop)
 223 {
 224     if (mtop == nullptr)
 225     {
 226         return false;
 227     }
 228     return mtop->moltype.empty() || mtop->moltype[0].atoms.haveMass;
 229 }
 230
 231 bool gmx_mtop_has_charges(const gmx_mtop_t *mtop)
 232 {
 233     if (mtop == nullptr)
 234     {
 235         return false;
 236     }
 237     return mtop->moltype.empty() || mtop->moltype[0].atoms.haveCharge;
 238 }
 239
 240 bool gmx_mtop_has_perturbed_charges(const gmx_mtop_t &mtop)
 241 {
 242     for (const gmx_moltype_t &moltype : mtop.moltype)
 243     {
 244         const t_atoms &atoms = moltype.atoms;
 245         if (atoms.haveBState)
 246         {
 247             for (int a = 0; a < atoms.nr; a++)
 248             {
 249                 if (atoms.atom[a].q != atoms.atom[a].qB)
 250                 {
 251                     return true;
 252                 }
 253             }
 254         }
 255     }
 256     return false;
 257 }
 258
 259 bool gmx_mtop_has_atomtypes(const gmx_mtop_t *mtop)
 260 {
 261     if (mtop == nullptr)
 262     {
 263         return false;
 264     }
 265     return mtop->moltype.empty() || mtop->moltype[0].atoms.haveType;
 266 }
 267
 268 bool gmx_mtop_has_pdbinfo(const gmx_mtop_t *mtop)
 269 {
 270     if (mtop == nullptr)
 271     {
 272         return false;
 273     }
 274     return mtop->moltype.empty() || mtop->moltype[0].atoms.havePdbInfo;
 275 }
 276
 277 static void pr_grps(FILE *fp, const char *title, const t_grps grps[], char **grpname[])
 278 {
 279     int i, j;
 280
 281     for (i = 0; (i < egcNR); i++)
 282     {
 283         fprintf(fp, "%s[%-12s] nr=%d, name=[", title, gtypes[i], grps[i].nr);
 284         for (j = 0; (j < grps[i].nr); j++)
 285         {
 286             fprintf(fp, " %s", *(grpname[grps[i].nm_ind[j]]));
 287         }
 288         fprintf(fp, "]\n");
 289     }
 290 }
 291
 292 static void pr_groups(FILE *fp, int indent,
 293                       const gmx_groups_t *groups,
 294                       gmx_bool bShowNumbers)
 295 {
 296     int nat_max, i, g;
 297
 298     pr_grps(fp, "grp", groups->grps, groups->grpname);
 299     pr_strings(fp, indent, "grpname", groups->grpname, groups->ngrpname, bShowNumbers);
 300
 301     pr_indent(fp, indent);
 302     fprintf(fp, "groups          ");
 303     for (g = 0; g < egcNR; g++)
 304     {
 305         printf(" %5.5s", gtypes[g]);
 306     }
 307     printf("\n");
 308
 309     pr_indent(fp, indent);
 310     fprintf(fp, "allocated       ");
 311     nat_max = 0;
 312     for (g = 0; g < egcNR; g++)
 313     {
 314         printf(" %5d", groups->ngrpnr[g]);
 315         nat_max = std::max(nat_max, groups->ngrpnr[g]);
 316     }
 317     printf("\n");
 318
 319     if (nat_max == 0)
 320     {
 321         pr_indent(fp, indent);
 322         fprintf(fp, "groupnr[%5s] =", "*");
 323         for (g = 0; g < egcNR; g++)
 324         {
 325             fprintf(fp, "  %3d ", 0);
 326         }
 327         fprintf(fp, "\n");
 328     }
 329     else
 330     {
 331         for (i = 0; i < nat_max; i++)
 332         {
 333             pr_indent(fp, indent);
 334             fprintf(fp, "groupnr[%5d] =", i);
 335             for (g = 0; g < egcNR; g++)
 336             {
 337                 fprintf(fp, "  %3d ",
 338                         groups->grpnr[g] ? groups->grpnr[g][i] : 0);
 339             }
 340             fprintf(fp, "\n");
 341         }
 342     }
 343 }
 344
 345 static void pr_moltype(FILE *fp, int indent, const char *title,
 346                        const gmx_moltype_t *molt, int n,
 347                        const gmx_ffparams_t *ffparams,
 348                        gmx_bool bShowNumbers, gmx_bool bShowParameters)
 349 {
 350     int j;
 351
 352     indent = pr_title_n(fp, indent, title, n);
 353     pr_indent(fp, indent);
 354     fprintf(fp, "name=\"%s\"\n", *(molt->name));
 355     pr_atoms(fp, indent, "atoms", &(molt->atoms), bShowNumbers);
 356     pr_block(fp, indent, "cgs", &molt->cgs, bShowNumbers);
 357     pr_blocka(fp, indent, "excls", &molt->excls, bShowNumbers);
 358     for (j = 0; (j < F_NRE); j++)
 359     {
 360         pr_ilist(fp, indent, interaction_function[j].longname,
 361                  ffparams->functype.data(), molt->ilist[j],
 362                  bShowNumbers, bShowParameters, ffparams->iparams.data());
 363     }
 364 }
 365
 366 static void pr_molblock(FILE *fp, int indent, const char *title,
 367                         const gmx_molblock_t *molb, int n,
 368                         const std::vector<gmx_moltype_t> &molt)
 369 {
 370     indent = pr_title_n(fp, indent, title, n);
 371     pr_indent(fp, indent);
 372     fprintf(fp, "%-20s = %d \"%s\"\n",
 373             "moltype", molb->type, *(molt[molb->type].name));
 374     pr_int(fp, indent, "#molecules", molb->nmol);
 375     pr_int(fp, indent, "#posres_xA", molb->posres_xA.size());
 376     if (!molb->posres_xA.empty())
 377     {
 378         pr_rvecs(fp, indent, "posres_xA", as_rvec_array(molb->posres_xA.data()), molb->posres_xA.size());
 379     }
 380     pr_int(fp, indent, "#posres_xB", molb->posres_xB.size());
 381     if (!molb->posres_xB.empty())
 382     {
 383         pr_rvecs(fp, indent, "posres_xB", as_rvec_array(molb->posres_xB.data()), molb->posres_xB.size());
 384     }
 385 }
 386
 387 void pr_mtop(FILE *fp, int indent, const char *title, const gmx_mtop_t *mtop,
 388              gmx_bool bShowNumbers, gmx_bool bShowParameters)
 389 {
 390     if (available(fp, mtop, indent, title))
 391     {
 392         indent = pr_title(fp, indent, title);
 393         pr_indent(fp, indent);
 394         fprintf(fp, "name=\"%s\"\n", *(mtop->name));
 395         pr_int(fp, indent, "#atoms", mtop->natoms);
 396         pr_int(fp, indent, "#molblock", mtop->molblock.size());
 397         for (size_t mb = 0; mb < mtop->molblock.size(); mb++)
 398         {
 399             pr_molblock(fp, indent, "molblock", &mtop->molblock[mb], mb, mtop->moltype);
 400         }
 401         pr_str(fp, indent, "bIntermolecularInteractions",
 402                gmx::boolToString(mtop->bIntermolecularInteractions));
 403         if (mtop->bIntermolecularInteractions)
 404         {
 405             for (int j = 0; j < F_NRE; j++)
 406             {
 407                 pr_ilist(fp, indent, interaction_function[j].longname,
 408                          mtop->ffparams.functype.data(),
 409                          (*mtop->intermolecular_ilist)[j],
 410                          bShowNumbers, bShowParameters, mtop->ffparams.iparams.data());
 411             }
 412         }
 413         pr_ffparams(fp, indent, "ffparams", &(mtop->ffparams), bShowNumbers);
 414         pr_atomtypes(fp, indent, "atomtypes", &(mtop->atomtypes), bShowNumbers);
 415         for (size_t mt = 0; mt < mtop->moltype.size(); mt++)
 416         {
 417             pr_moltype(fp, indent, "moltype", &mtop->moltype[mt], mt,
 418                        &mtop->ffparams, bShowNumbers, bShowParameters);
 419         }
 420         pr_groups(fp, indent, &mtop->groups, bShowNumbers);
 421     }
 422 }
 423
 424 void pr_top(FILE *fp, int indent, const char *title, const t_topology *top,
 425             gmx_bool bShowNumbers, gmx_bool bShowParameters)
 426 {
 427     if (available(fp, top, indent, title))
 428     {
 429         indent = pr_title(fp, indent, title);
 430         pr_indent(fp, indent);
 431         fprintf(fp, "name=\"%s\"\n", *(top->name));
 432         pr_atoms(fp, indent, "atoms", &(top->atoms), bShowNumbers);
 433         pr_atomtypes(fp, indent, "atomtypes", &(top->atomtypes), bShowNumbers);
 434         pr_block(fp, indent, "cgs", &top->cgs, bShowNumbers);
 435         pr_block(fp, indent, "mols", &top->mols, bShowNumbers);
 436         pr_str(fp, indent, "bIntermolecularInteractions",
 437                gmx::boolToString(top->bIntermolecularInteractions));
 438         pr_blocka(fp, indent, "excls", &top->excls, bShowNumbers);
 439         pr_idef(fp, indent, "idef", &top->idef, bShowNumbers, bShowParameters);
 440     }
 441 }
 442
 443 static void cmp_ilist(FILE *fp, int ftype, const t_ilist *il1, const t_ilist *il2)
 444 {
 445     int  i;
 446     char buf[256];
 447
 448     fprintf(fp, "comparing ilist %s\n", interaction_function[ftype].name);
 449     sprintf(buf, "%s->nr", interaction_function[ftype].name);
 450     cmp_int(fp, buf, -1, il1->nr, il2->nr);
 451     sprintf(buf, "%s->iatoms", interaction_function[ftype].name);
 452     if (((il1->nr > 0) && (!il1->iatoms)) ||
 453         ((il2->nr > 0) && (!il2->iatoms)) ||
 454         ((il1->nr != il2->nr)))
 455     {
 456         fprintf(fp, "Comparing radically different topologies - %s is different\n",
 457                 buf);
 458     }
 459     else
 460     {
 461         for (i = 0; (i < il1->nr); i++)
 462         {
 463             cmp_int(fp, buf, i, il1->iatoms[i], il2->iatoms[i]);
 464         }
 465     }
 466 }
 467
 468 static void cmp_iparm(FILE *fp, const char *s, t_functype ft,
 469                       const t_iparams &ip1, const t_iparams &ip2, real ftol, real abstol)
 470 {
 471     int      i;
 472     gmx_bool bDiff;
 473
 474     bDiff = FALSE;
 475     for (i = 0; i < MAXFORCEPARAM && !bDiff; i++)
 476     {
 477         bDiff = !equal_real(ip1.generic.buf[i], ip2.generic.buf[i], ftol, abstol);
 478     }
 479     if (bDiff)
 480     {
 481         fprintf(fp, "%s1: ", s);
 482         pr_iparams(fp, ft, &ip1);
 483         fprintf(fp, "%s2: ", s);
 484         pr_iparams(fp, ft, &ip2);
 485     }
 486 }
 487
 488 static void cmp_iparm_AB(FILE *fp, const char *s, t_functype ft,
 489                          const t_iparams &ip1, real ftol, real abstol)
 490 {
 491     int      nrfpA, nrfpB, p0, i;
 492     gmx_bool bDiff;
 493
 494     /* Normally the first parameter is perturbable */
 495     p0    = 0;
 496     nrfpA = interaction_function[ft].nrfpA;
 497     nrfpB = interaction_function[ft].nrfpB;
 498     if (ft == F_PDIHS)
 499     {
 500         nrfpB = 2;
 501     }
 502     else if (interaction_function[ft].flags & IF_TABULATED)
 503     {
 504         /* For tabulated interactions only the second parameter is perturbable */
 505         p0    = 1;
 506         nrfpB = 1;
 507     }
 508     bDiff = FALSE;
 509     for (i = 0; i < nrfpB && !bDiff; i++)
 510     {
 511         bDiff = !equal_real(ip1.generic.buf[p0+i], ip1.generic.buf[nrfpA+i], ftol, abstol);
 512     }
 513     if (bDiff)
 514     {
 515         fprintf(fp, "%s: ", s);
 516         pr_iparams(fp, ft, &ip1);
 517     }
 518 }
 519
 520 static void cmp_cmap(FILE *fp, const gmx_cmap_t *cmap1, const gmx_cmap_t *cmap2, real ftol, real abstol)
 521 {
 522     int cmap1_ngrid = (cmap1 ? cmap1->cmapdata.size() : 0);
 523     int cmap2_ngrid = (cmap2 ? cmap2->cmapdata.size() : 0);
 524
 525     cmp_int(fp, "cmap ngrid", -1, cmap1_ngrid, cmap2_ngrid);
 526
 527     if (cmap1 == nullptr || cmap2 == nullptr)
 528     {
 529         return;
 530     }
 531
 532     cmp_int(fp, "cmap grid_spacing", -1, cmap1->grid_spacing, cmap2->grid_spacing);
 533     if (cmap1->cmapdata.size() == cmap2->cmapdata.size() &&
 534         cmap1->grid_spacing == cmap2->grid_spacing)
 535     {
 536         for (size_t g = 0; g < cmap1->cmapdata.size(); g++)
 537         {
 538             int i;
 539
 540             fprintf(fp, "comparing cmap %zu\n", g);
 541
 542             for (i = 0; i < 4*cmap1->grid_spacing*cmap1->grid_spacing; i++)
 543             {
 544                 cmp_real(fp, "", i, cmap1->cmapdata[g].cmap[i], cmap2->cmapdata[g].cmap[i], ftol, abstol);
 545             }
 546         }
 547     }
 548 }
 549
 550 static void cmp_idef(FILE *fp, const t_idef *id1, const t_idef *id2, real ftol, real abstol)
 551 {
 552     int  i;
 553     char buf1[64], buf2[64];
 554
 555     fprintf(fp, "comparing idef\n");
 556     if (id2)
 557     {
 558         cmp_int(fp, "idef->ntypes", -1, id1->ntypes, id2->ntypes);
 559         cmp_int(fp, "idef->atnr",  -1, id1->atnr, id2->atnr);
 560         for (i = 0; (i < std::min(id1->ntypes, id2->ntypes)); i++)
 561         {
 562             sprintf(buf1, "idef->functype[%d]", i);
 563             sprintf(buf2, "idef->iparam[%d]", i);
 564             cmp_int(fp, buf1, i, static_cast<int>(id1->functype[i]), static_cast<int>(id2->functype[i]));
 565             cmp_iparm(fp, buf2, id1->functype[i],
 566                       id1->iparams[i], id2->iparams[i], ftol, abstol);
 567         }
 568         cmp_real(fp, "fudgeQQ", -1, id1->fudgeQQ, id2->fudgeQQ, ftol, abstol);
 569         cmp_cmap(fp, id1->cmap_grid, id2->cmap_grid, ftol, abstol);
 570         for (i = 0; (i < F_NRE); i++)
 571         {
 572             cmp_ilist(fp, i, &(id1->il[i]), &(id2->il[i]));
 573         }
 574     }
 575     else
 576     {
 577         for (i = 0; (i < id1->ntypes); i++)
 578         {
 579             cmp_iparm_AB(fp, "idef->iparam", id1->functype[i], id1->iparams[i], ftol, abstol);
 580         }
 581     }
 582 }
 583
 584 static void cmp_block(FILE *fp, const t_block *b1, const t_block *b2, const char *s)
 585 {
 586     char buf[32];
 587
 588     fprintf(fp, "comparing block %s\n", s);
 589     sprintf(buf, "%s.nr", s);
 590     cmp_int(fp, buf, -1, b1->nr, b2->nr);
 591 }
 592
 593 static void cmp_blocka(FILE *fp, const t_blocka *b1, const t_blocka *b2, const char *s)
 594 {
 595     char buf[32];
 596
 597     fprintf(fp, "comparing blocka %s\n", s);
 598     sprintf(buf, "%s.nr", s);
 599     cmp_int(fp, buf, -1, b1->nr, b2->nr);
 600     sprintf(buf, "%s.nra", s);
 601     cmp_int(fp, buf, -1, b1->nra, b2->nra);
 602 }
 603
 604 void cmp_top(FILE *fp, const t_topology *t1, const t_topology *t2, real ftol, real abstol)
 605 {
 606     fprintf(fp, "comparing top\n");
 607     if (t2)
 608     {
 609         cmp_idef(fp, &(t1->idef), &(t2->idef), ftol, abstol);
 610         cmp_atoms(fp, &(t1->atoms), &(t2->atoms), ftol, abstol);
 611         cmp_block(fp, &t1->cgs, &t2->cgs, "cgs");
 612         cmp_block(fp, &t1->mols, &t2->mols, "mols");
 613         cmp_bool(fp, "bIntermolecularInteractions", -1, t1->bIntermolecularInteractions, t2->bIntermolecularInteractions);
 614         cmp_blocka(fp, &t1->excls, &t2->excls, "excls");
 615     }
 616     else
 617     {
 618         cmp_idef(fp, &(t1->idef), nullptr, ftol, abstol);
 619         cmp_atoms(fp, &(t1->atoms), nullptr, ftol, abstol);
 620     }
 621 }
 622
 623 void cmp_groups(FILE *fp, const gmx_groups_t *g0, const gmx_groups_t *g1,
 624                 int natoms0, int natoms1)
 625 {
 626     char buf[32];
 627
 628     fprintf(fp, "comparing groups\n");
 629
 630     for (int i = 0; i < egcNR; i++)
 631     {
 632         sprintf(buf, "grps[%d].nr", i);
 633         cmp_int(fp, buf, -1, g0->grps[i].nr, g1->grps[i].nr);
 634         if (g0->grps[i].nr == g1->grps[i].nr)
 635         {
 636             for (int j = 0; j < g0->grps[i].nr; j++)
 637             {
 638                 sprintf(buf, "grps[%d].name[%d]", i, j);
 639                 cmp_str(fp, buf, -1,
 640                         *g0->grpname[g0->grps[i].nm_ind[j]],
 641                         *g1->grpname[g1->grps[i].nm_ind[j]]);
 642             }
 643         }
 644         cmp_int(fp, "ngrpnr", i, g0->ngrpnr[i], g1->ngrpnr[i]);
 645         if (g0->ngrpnr[i] == g1->ngrpnr[i] && natoms0 == natoms1 &&
 646             (g0->grpnr[i] != nullptr || g1->grpnr[i] != nullptr))
 647         {
 648             for (int j = 0; j < natoms0; j++)
 649             {
 650                 cmp_int(fp, gtypes[i], j, getGroupType(g0, i, j), getGroupType(g1, i, j));
 651             }
 652         }
 653     }
 654     /* We have compared the names in the groups lists,
 655      * so we can skip the grpname list comparison.
 656      */
 657 }
 658
 659 int getGroupType(const gmx_groups_t *group, int type, int atom)
 660 {
 661     return (group->grpnr[type] ? group->grpnr[type][atom] : 0);
 662 }