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