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