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 by the GROMACS development team.
   7  * Copyright (c) 2018,2019,2020, by the GROMACS development team, led by
   8  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   9  * and including many others, as listed in the AUTHORS file in the
  10  * top-level source directory and at http://www.gromacs.org.
  11  *
  12  * GROMACS is free software; you can redistribute it and/or
  13  * modify it under the terms of the GNU Lesser General Public License
  14  * as published by the Free Software Foundation; either version 2.1
  15  * of the License, or (at your option) any later version.
  16  *
  17  * GROMACS is distributed in the hope that it will be useful,
  18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  20  * Lesser General Public License for more details.
  21  *
  22  * You should have received a copy of the GNU Lesser General Public
  23  * License along with GROMACS; if not, see
  24  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  25  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  26  *
  27  * If you want to redistribute modifications to GROMACS, please
  28  * consider that scientific software is very special. Version
  29  * control is crucial - bugs must be traceable. We will be happy to
  30  * consider code for inclusion in the official distribution, but
  31  * derived work must not be called official GROMACS. Details are found
  32  * in the README & COPYING files - if they are missing, get the
  33  * official version at http://www.gromacs.org.
  34  *
  35  * To help us fund GROMACS development, we humbly ask that you cite
  36  * the research papers on the package. Check out http://www.gromacs.org.
  37  */
  38 #include "gmxpre.h"
  39
  40 #include "topology.h"
  41
  42 #include <cstdio>
  43
  44 #include <algorithm>
  45
  46 #include "gromacs/math/vecdump.h"
  47 #include "gromacs/topology/atoms.h"
  48 #include "gromacs/topology/idef.h"
  49 #include "gromacs/topology/ifunc.h"
  50 #include "gromacs/topology/symtab.h"
  51 #include "gromacs/utility/arrayref.h"
  52 #include "gromacs/utility/compare.h"
  53 #include "gromacs/utility/gmxassert.h"
  54 #include "gromacs/utility/smalloc.h"
  55 #include "gromacs/utility/strconvert.h"
  56 #include "gromacs/utility/txtdump.h"
  57
  58 static gmx::EnumerationArray<SimulationAtomGroupType, const char*> c_simulationAtomGroupTypeShortNames = {
  59     { "T-Coupling", "Energy Mon.", "Acceleration", "Freeze", "User1", "User2", "VCM",
  60       "Compressed X", "Or. Res. Fit", "QMMM" }
  61 };
  62
  63 const char* shortName(SimulationAtomGroupType type)
  64 {
  65     return c_simulationAtomGroupTypeShortNames[type];
  66 }
  67
  68 void init_top(t_topology* top)
  69 {
  70     top->name = nullptr;
  71     init_idef(&top->idef);
  72     init_atom(&(top->atoms));
  73     init_atomtypes(&(top->atomtypes));
  74     init_block(&top->mols);
  75     open_symtab(&top->symtab);
  76 }
  77
  78
  79 gmx_moltype_t::gmx_moltype_t() : name(nullptr)
  80 {
  81     init_t_atoms(&atoms, 0, FALSE);
  82 }
  83
  84 gmx_moltype_t::~gmx_moltype_t()
  85 {
  86     done_atom(&atoms);
  87 }
  88
  89 gmx_mtop_t::gmx_mtop_t()
  90 {
  91     init_atomtypes(&atomtypes);
  92     open_symtab(&symtab);
  93 }
  94
  95 gmx_mtop_t::~gmx_mtop_t()
  96 {
  97     done_symtab(&symtab);
  98
  99     moltype.clear();
 100     molblock.clear();
 101     done_atomtypes(&atomtypes);
 102 }
 103
 104 void done_top(t_topology* top)
 105 {
 106     done_idef(&top->idef);
 107     done_atom(&(top->atoms));
 108
 109     /* For GB */
 110     done_atomtypes(&(top->atomtypes));
 111
 112     done_symtab(&(top->symtab));
 113     done_block(&(top->mols));
 114 }
 115
 116 void done_top_mtop(t_topology* top, gmx_mtop_t* mtop)
 117 {
 118     if (mtop != nullptr)
 119     {
 120         if (top != nullptr)
 121         {
 122             done_idef(&top->idef);
 123             done_atom(&top->atoms);
 124             done_block(&top->mols);
 125             done_symtab(&top->symtab);
 126             open_symtab(&mtop->symtab);
 127             done_atomtypes(&(top->atomtypes));
 128         }
 129
 130         // Note that the rest of mtop will be freed by the destructor
 131     }
 132 }
 133
 134 gmx_localtop_t::gmx_localtop_t(const gmx_ffparams_t& ffparams) : idef(ffparams) {}
 135
 136 bool gmx_mtop_has_masses(const gmx_mtop_t* mtop)
 137 {
 138     if (mtop == nullptr)
 139     {
 140         return false;
 141     }
 142     return mtop->moltype.empty() || mtop->moltype[0].atoms.haveMass;
 143 }
 144
 145 bool gmx_mtop_has_charges(const gmx_mtop_t* mtop)
 146 {
 147     if (mtop == nullptr)
 148     {
 149         return false;
 150     }
 151     return mtop->moltype.empty() || mtop->moltype[0].atoms.haveCharge;
 152 }
 153
 154 bool gmx_mtop_has_perturbed_charges(const gmx_mtop_t& mtop)
 155 {
 156     for (const gmx_moltype_t& moltype : mtop.moltype)
 157     {
 158         const t_atoms& atoms = moltype.atoms;
 159         if (atoms.haveBState)
 160         {
 161             for (int a = 0; a < atoms.nr; a++)
 162             {
 163                 if (atoms.atom[a].q != atoms.atom[a].qB)
 164                 {
 165                     return true;
 166                 }
 167             }
 168         }
 169     }
 170     return false;
 171 }
 172
 173 bool gmx_mtop_has_atomtypes(const gmx_mtop_t* mtop)
 174 {
 175     if (mtop == nullptr)
 176     {
 177         return false;
 178     }
 179     return mtop->moltype.empty() || mtop->moltype[0].atoms.haveType;
 180 }
 181
 182 bool gmx_mtop_has_pdbinfo(const gmx_mtop_t* mtop)
 183 {
 184     if (mtop == nullptr)
 185     {
 186         return false;
 187     }
 188     return mtop->moltype.empty() || mtop->moltype[0].atoms.havePdbInfo;
 189 }
 190
 191 static void pr_grps(FILE* fp, const char* title, gmx::ArrayRef<const AtomGroupIndices> grps, char*** grpname)
 192 {
 193     int index = 0;
 194     for (const auto& group : grps)
 195     {
 196         fprintf(fp, "%s[%-12s] nr=%zu, name=[", title, c_simulationAtomGroupTypeShortNames[index],
 197                 group.size());
 198         for (const auto& entry : group)
 199         {
 200             fprintf(fp, " %s", *(grpname[entry]));
 201         }
 202         fprintf(fp, "]\n");
 203         index++;
 204     }
 205 }
 206
 207 static void pr_groups(FILE* fp, int indent, const SimulationGroups& groups, gmx_bool bShowNumbers)
 208 {
 209     pr_grps(fp, "grp", groups.groups, const_cast<char***>(groups.groupNames.data()));
 210     pr_strings(fp, indent, "grpname", const_cast<char***>(groups.groupNames.data()),
 211                groups.groupNames.size(), bShowNumbers);
 212
 213     pr_indent(fp, indent);
 214     fprintf(fp, "groups          ");
 215     for (const auto& group : c_simulationAtomGroupTypeShortNames)
 216     {
 217         printf(" %5.5s", group);
 218     }
 219     printf("\n");
 220
 221     pr_indent(fp, indent);
 222     fprintf(fp, "allocated       ");
 223     int nat_max = 0;
 224     for (auto group : keysOf(groups.groups))
 225     {
 226         printf(" %5d", groups.numberOfGroupNumbers(group));
 227         nat_max = std::max(nat_max, groups.numberOfGroupNumbers(group));
 228     }
 229     printf("\n");
 230
 231     if (nat_max == 0)
 232     {
 233         pr_indent(fp, indent);
 234         fprintf(fp, "groupnr[%5s] =", "*");
 235         for (auto gmx_unused group : keysOf(groups.groups))
 236         {
 237             fprintf(fp, "  %3d ", 0);
 238         }
 239         fprintf(fp, "\n");
 240     }
 241     else
 242     {
 243         for (int i = 0; i < nat_max; i++)
 244         {
 245             pr_indent(fp, indent);
 246             fprintf(fp, "groupnr[%5d] =", i);
 247             for (auto group : keysOf(groups.groups))
 248             {
 249                 fprintf(fp, "  %3d ",
 250                         !groups.groupNumbers[group].empty() ? groups.groupNumbers[group][i] : 0);
 251             }
 252             fprintf(fp, "\n");
 253         }
 254     }
 255 }
 256
 257 static void pr_moltype(FILE*                 fp,
 258                        int                   indent,
 259                        const char*           title,
 260                        const gmx_moltype_t*  molt,
 261                        int                   n,
 262                        const gmx_ffparams_t* ffparams,
 263                        gmx_bool              bShowNumbers,
 264                        gmx_bool              bShowParameters)
 265 {
 266     int j;
 267
 268     indent = pr_title_n(fp, indent, title, n);
 269     pr_indent(fp, indent);
 270     fprintf(fp, "name=\"%s\"\n", *(molt->name));
 271     pr_atoms(fp, indent, "atoms", &(molt->atoms), bShowNumbers);
 272     pr_listoflists(fp, indent, "excls", &molt->excls, bShowNumbers);
 273     for (j = 0; (j < F_NRE); j++)
 274     {
 275         pr_ilist(fp, indent, interaction_function[j].longname, ffparams->functype.data(),
 276                  molt->ilist[j], bShowNumbers, bShowParameters, ffparams->iparams.data());
 277     }
 278 }
 279
 280 static void pr_molblock(FILE*                             fp,
 281                         int                               indent,
 282                         const char*                       title,
 283                         const gmx_molblock_t*             molb,
 284                         int                               n,
 285                         const std::vector<gmx_moltype_t>& molt)
 286 {
 287     indent = pr_title_n(fp, indent, title, n);
 288     pr_indent(fp, indent);
 289     fprintf(fp, "%-20s = %d \"%s\"\n", "moltype", molb->type, *(molt[molb->type].name));
 290     pr_int(fp, indent, "#molecules", molb->nmol);
 291     pr_int(fp, indent, "#posres_xA", molb->posres_xA.size());
 292     if (!molb->posres_xA.empty())
 293     {
 294         pr_rvecs(fp, indent, "posres_xA", as_rvec_array(molb->posres_xA.data()), molb->posres_xA.size());
 295     }
 296     pr_int(fp, indent, "#posres_xB", molb->posres_xB.size());
 297     if (!molb->posres_xB.empty())
 298     {
 299         pr_rvecs(fp, indent, "posres_xB", as_rvec_array(molb->posres_xB.data()), molb->posres_xB.size());
 300     }
 301 }
 302
 303 void pr_mtop(FILE* fp, int indent, const char* title, const gmx_mtop_t* mtop, gmx_bool bShowNumbers, gmx_bool bShowParameters)
 304 {
 305     if (available(fp, mtop, indent, title))
 306     {
 307         indent = pr_title(fp, indent, title);
 308         pr_indent(fp, indent);
 309         fprintf(fp, "name=\"%s\"\n", *(mtop->name));
 310         pr_int(fp, indent, "#atoms", mtop->natoms);
 311         pr_int(fp, indent, "#molblock", mtop->molblock.size());
 312         for (size_t mb = 0; mb < mtop->molblock.size(); mb++)
 313         {
 314             pr_molblock(fp, indent, "molblock", &mtop->molblock[mb], mb, mtop->moltype);
 315         }
 316         pr_str(fp, indent, "bIntermolecularInteractions",
 317                gmx::boolToString(mtop->bIntermolecularInteractions));
 318         if (mtop->bIntermolecularInteractions)
 319         {
 320             for (int j = 0; j < F_NRE; j++)
 321             {
 322                 pr_ilist(fp, indent, interaction_function[j].longname,
 323                          mtop->ffparams.functype.data(), (*mtop->intermolecular_ilist)[j],
 324                          bShowNumbers, bShowParameters, mtop->ffparams.iparams.data());
 325             }
 326         }
 327         pr_ffparams(fp, indent, "ffparams", &(mtop->ffparams), bShowNumbers);
 328         pr_atomtypes(fp, indent, "atomtypes", &(mtop->atomtypes), bShowNumbers);
 329         for (size_t mt = 0; mt < mtop->moltype.size(); mt++)
 330         {
 331             pr_moltype(fp, indent, "moltype", &mtop->moltype[mt], mt, &mtop->ffparams, bShowNumbers,
 332                        bShowParameters);
 333         }
 334         pr_groups(fp, indent, mtop->groups, bShowNumbers);
 335     }
 336 }
 337
 338 void pr_top(FILE* fp, int indent, const char* title, const t_topology* top, gmx_bool bShowNumbers, gmx_bool bShowParameters)
 339 {
 340     if (available(fp, top, indent, title))
 341     {
 342         indent = pr_title(fp, indent, title);
 343         pr_indent(fp, indent);
 344         fprintf(fp, "name=\"%s\"\n", *(top->name));
 345         pr_atoms(fp, indent, "atoms", &(top->atoms), bShowNumbers);
 346         pr_atomtypes(fp, indent, "atomtypes", &(top->atomtypes), bShowNumbers);
 347         pr_block(fp, indent, "mols", &top->mols, bShowNumbers);
 348         pr_str(fp, indent, "bIntermolecularInteractions",
 349                gmx::boolToString(top->bIntermolecularInteractions));
 350         pr_idef(fp, indent, "idef", &top->idef, bShowNumbers, bShowParameters);
 351     }
 352 }
 353
 354 static void cmp_iparm(FILE*            fp,
 355                       const char*      s,
 356                       t_functype       ft,
 357                       const t_iparams& ip1,
 358                       const t_iparams& ip2,
 359                       real             relativeTolerance,
 360                       real             absoluteTolerance)
 361 {
 362     int      i;
 363     gmx_bool bDiff;
 364
 365     bDiff = FALSE;
 366     for (i = 0; i < MAXFORCEPARAM && !bDiff; i++)
 367     {
 368         bDiff = !equal_real(ip1.generic.buf[i], ip2.generic.buf[i], relativeTolerance, absoluteTolerance);
 369     }
 370     if (bDiff)
 371     {
 372         fprintf(fp, "%s1: ", s);
 373         pr_iparams(fp, ft, ip1);
 374         fprintf(fp, "%s2: ", s);
 375         pr_iparams(fp, ft, ip2);
 376     }
 377 }
 378
 379 static void cmp_iparm_AB(FILE* fp, const char* s, t_functype ft, const t_iparams& ip1, real relativeTolerance, real absoluteTolerance)
 380 {
 381     int      nrfpA, nrfpB, p0, i;
 382     gmx_bool bDiff;
 383
 384     /* Normally the first parameter is perturbable */
 385     p0    = 0;
 386     nrfpA = interaction_function[ft].nrfpA;
 387     nrfpB = interaction_function[ft].nrfpB;
 388     if (ft == F_PDIHS)
 389     {
 390         nrfpB = 2;
 391     }
 392     else if (interaction_function[ft].flags & IF_TABULATED)
 393     {
 394         /* For tabulated interactions only the second parameter is perturbable */
 395         p0    = 1;
 396         nrfpB = 1;
 397     }
 398     bDiff = FALSE;
 399     for (i = 0; i < nrfpB && !bDiff; i++)
 400     {
 401         bDiff = !equal_real(ip1.generic.buf[p0 + i], ip1.generic.buf[nrfpA + i], relativeTolerance,
 402                             absoluteTolerance);
 403     }
 404     if (bDiff)
 405     {
 406         fprintf(fp, "%s: ", s);
 407         pr_iparams(fp, ft, ip1);
 408     }
 409 }
 410
 411 static void cmp_cmap(FILE* fp, const gmx_cmap_t* cmap1, const gmx_cmap_t* cmap2, real relativeTolerance, real absoluteTolerance)
 412 {
 413     int cmap1_ngrid = (cmap1 ? cmap1->cmapdata.size() : 0);
 414     int cmap2_ngrid = (cmap2 ? cmap2->cmapdata.size() : 0);
 415
 416     cmp_int(fp, "cmap ngrid", -1, cmap1_ngrid, cmap2_ngrid);
 417
 418     if (cmap1 == nullptr || cmap2 == nullptr)
 419     {
 420         return;
 421     }
 422
 423     cmp_int(fp, "cmap grid_spacing", -1, cmap1->grid_spacing, cmap2->grid_spacing);
 424     if (cmap1->cmapdata.size() == cmap2->cmapdata.size() && cmap1->grid_spacing == cmap2->grid_spacing)
 425     {
 426         for (size_t g = 0; g < cmap1->cmapdata.size(); g++)
 427         {
 428             int i;
 429
 430             fprintf(fp, "comparing cmap %zu\n", g);
 431
 432             for (i = 0; i < 4 * cmap1->grid_spacing * cmap1->grid_spacing; i++)
 433             {
 434                 cmp_real(fp, "", i, cmap1->cmapdata[g].cmap[i], cmap2->cmapdata[g].cmap[i],
 435                          relativeTolerance, absoluteTolerance);
 436             }
 437         }
 438     }
 439 }
 440
 441 static void cmp_listoflists(FILE*                        fp,
 442                             const gmx::ListOfLists<int>& list1,
 443                             const gmx::ListOfLists<int>& list2,
 444                             const char*                  s)
 445 {
 446     char buf[32];
 447
 448     fprintf(fp, "comparing blocka %s\n", s);
 449     sprintf(buf, "%s.numLists", s);
 450     cmp_int(fp, buf, -1, list1.ssize(), list2.ssize());
 451     sprintf(buf, "%s.numElements", s);
 452     cmp_int(fp, buf, -1, list1.numElements(), list2.numElements());
 453 }
 454
 455 static void compareFfparams(FILE*                 fp,
 456                             const gmx_ffparams_t& ff1,
 457                             const gmx_ffparams_t& ff2,
 458                             real                  relativeTolerance,
 459                             real                  absoluteTolerance)
 460 {
 461     fprintf(fp, "comparing force field parameters\n");
 462     cmp_int(fp, "numTypes", -1, ff1.numTypes(), ff2.numTypes());
 463     cmp_int(fp, "atnr", -1, ff1.atnr, ff1.atnr);
 464     cmp_double(fp, "reppow", -1, ff1.reppow, ff2.reppow, relativeTolerance, absoluteTolerance);
 465     cmp_real(fp, "fudgeQQ", -1, ff1.fudgeQQ, ff2.fudgeQQ, relativeTolerance, absoluteTolerance);
 466     cmp_cmap(fp, &ff1.cmap_grid, &ff2.cmap_grid, relativeTolerance, absoluteTolerance);
 467     for (int i = 0; i < std::min(ff1.numTypes(), ff2.numTypes()); i++)
 468     {
 469         std::string buf = gmx::formatString("ffparams->functype[%d]", i);
 470         cmp_int(fp, buf.c_str(), i, ff1.functype[i], ff2.functype[i]);
 471         buf = gmx::formatString("ffparams->iparams[%d]", i);
 472         cmp_iparm(fp, buf.c_str(), ff1.functype[i], ff1.iparams[i], ff2.iparams[i],
 473                   relativeTolerance, absoluteTolerance);
 474     }
 475 }
 476
 477 static void compareFfparamAB(FILE* fp, const gmx_ffparams_t& ff1, real relativeTolerance, real absoluteTolerance)
 478 {
 479     fprintf(fp, "comparing free energy parameters\n");
 480     for (int i = 0; i < ff1.numTypes(); i++)
 481     {
 482         std::string buf = gmx::formatString("ffparams->iparams[%d]", i);
 483         cmp_iparm_AB(fp, buf.c_str(), ff1.functype[i], ff1.iparams[i], relativeTolerance, absoluteTolerance);
 484     }
 485 }
 486 static void compareInteractionLists(FILE* fp, const InteractionLists* il1, const InteractionLists* il2)
 487 {
 488     fprintf(fp, "comparing InteractionLists\n");
 489     if ((il1 || il2) && (!il1 || !il2))
 490     {
 491         fprintf(fp, "InteractionLists are present in topology %d but not in the other\n", il1 ? 1 : 2);
 492     }
 493     if (il1 && il2)
 494     {
 495         for (int i = 0; i < F_NRE; i++)
 496         {
 497             cmp_int(fp, "InteractionList size", i, il1->at(i).size(), il2->at(i).size());
 498             int nr = std::min(il1->at(i).size(), il2->at(i).size());
 499             for (int j = 0; j < nr; j++)
 500             {
 501                 cmp_int(fp, "InteractionList entry", j, il1->at(i).iatoms.at(j), il2->at(i).iatoms.at(j));
 502             }
 503         }
 504     }
 505 }
 506
 507 static void compareMoltypes(FILE*                              fp,
 508                             gmx::ArrayRef<const gmx_moltype_t> mt1,
 509                             gmx::ArrayRef<const gmx_moltype_t> mt2,
 510                             real                               relativeTolerance,
 511                             real                               absoluteTolerance)
 512 {
 513     fprintf(fp, "comparing molecule types\n");
 514     cmp_int(fp, "moltype size", -1, mt1.size(), mt2.size());
 515     for (int i = 0; i < std::min(mt1.ssize(), mt2.ssize()); i++)
 516     {
 517         cmp_str(fp, "Name", i, *mt1[i].name, *mt2[i].name);
 518         compareAtoms(fp, &mt1[i].atoms, &mt2[i].atoms, relativeTolerance, absoluteTolerance);
 519         compareInteractionLists(fp, &mt1[i].ilist, &mt2[i].ilist);
 520         std::string buf = gmx::formatString("excls[%d]", i);
 521         cmp_listoflists(fp, mt1[i].excls, mt2[i].excls, buf.c_str());
 522     }
 523 }
 524
 525 static void compareMoletypeAB(FILE* fp, gmx::ArrayRef<const gmx_moltype_t> mt1, real relativeTolerance, real absoluteTolerance)
 526 {
 527     fprintf(fp, "comparing free energy molecule types\n");
 528     for (gmx::index i = 0; i < mt1.ssize(); i++)
 529     {
 530         compareAtoms(fp, &mt1[i].atoms, nullptr, relativeTolerance, absoluteTolerance);
 531     }
 532 }
 533 static void compareMolblocks(FILE*                               fp,
 534                              gmx::ArrayRef<const gmx_molblock_t> mb1,
 535                              gmx::ArrayRef<const gmx_molblock_t> mb2)
 536 {
 537     fprintf(fp, "comparing molecule blocks\n");
 538     cmp_int(fp, "molblock size", -1, mb1.size(), mb2.size());
 539     int nr = std::min(mb1.size(), mb2.size());
 540     for (int i = 0; i < nr; i++)
 541     {
 542         cmp_int(fp, "type", i, mb1[i].type, mb2[i].type);
 543         cmp_int(fp, "nmol", i, mb1[i].nmol, mb2[i].nmol);
 544         // Only checking size of restraint vectors for now
 545         cmp_int(fp, "posres_xA size", i, mb1[i].posres_xA.size(), mb2[i].posres_xA.size());
 546         cmp_int(fp, "posres_xB size", i, mb1[i].posres_xB.size(), mb2[i].posres_xB.size());
 547     }
 548 }
 549
 550 static void compareAtomtypes(FILE* fp, const t_atomtypes& at1, const t_atomtypes& at2)
 551 {
 552     fprintf(fp, "comparing atomtypes\n");
 553     cmp_int(fp, "nr", -1, at1.nr, at2.nr);
 554     int nr = std::min(at1.nr, at2.nr);
 555     for (int i = 0; i < nr; i++)
 556     {
 557         cmp_int(fp, "atomtype", i, at1.atomnumber[i], at2.atomnumber[i]);
 558     }
 559 }
 560
 561 static void compareIntermolecularExclusions(FILE*                    fp,
 562                                             gmx::ArrayRef<const int> ime1,
 563                                             gmx::ArrayRef<const int> ime2)
 564 {
 565     fprintf(fp, "comparing intermolecular exclusions\n");
 566     cmp_int(fp, "exclusion number", -1, ime1.size(), ime2.size());
 567     int nr = std::min(ime1.size(), ime2.size());
 568     for (int i = 0; i < nr; i++)
 569     {
 570         cmp_int(fp, "exclusion", i, ime1[i], ime2[i]);
 571     }
 572 }
 573
 574 static void compareBlockIndices(FILE*                                     fp,
 575                                 gmx::ArrayRef<const MoleculeBlockIndices> mbi1,
 576                                 gmx::ArrayRef<const MoleculeBlockIndices> mbi2)
 577 {
 578     fprintf(fp, "comparing moleculeBlockIndices\n");
 579     cmp_int(fp, "size", -1, mbi1.size(), mbi2.size());
 580     int nr = std::min(mbi1.size(), mbi2.size());
 581     for (int i = 0; i < nr; i++)
 582     {
 583         cmp_int(fp, "numAtomsPerMolecule", i, mbi1[i].numAtomsPerMolecule, mbi2[i].numAtomsPerMolecule);
 584         cmp_int(fp, "globalAtomStart", i, mbi1[i].globalAtomStart, mbi2[i].globalAtomStart);
 585         cmp_int(fp, "globalAtomEnd", i, mbi1[i].globalAtomEnd, mbi2[i].globalAtomEnd);
 586         cmp_int(fp, "globalResidueStart", i, mbi1[i].globalResidueStart, mbi2[i].globalResidueStart);
 587         cmp_int(fp, "moleculeIndexStart", i, mbi1[i].moleculeIndexStart, mbi2[i].moleculeIndexStart);
 588     }
 589 }
 590
 591 void compareMtop(FILE* fp, const gmx_mtop_t& mtop1, const gmx_mtop_t& mtop2, real relativeTolerance, real absoluteTolerance)
 592 {
 593     fprintf(fp, "comparing mtop topology\n");
 594     cmp_str(fp, "Name", -1, *mtop1.name, *mtop2.name);
 595     cmp_int(fp, "natoms", -1, mtop1.natoms, mtop2.natoms);
 596     cmp_int(fp, "maxres_renum", -1, mtop1.maxres_renum, mtop2.maxres_renum);
 597     cmp_int(fp, "maxresnr", -1, mtop1.maxresnr, mtop2.maxresnr);
 598     cmp_bool(fp, "bIntermolecularInteractions", -1, mtop1.bIntermolecularInteractions,
 599              mtop2.bIntermolecularInteractions);
 600     cmp_bool(fp, "haveMoleculeIndices", -1, mtop1.haveMoleculeIndices, mtop2.haveMoleculeIndices);
 601
 602     compareFfparams(fp, mtop1.ffparams, mtop2.ffparams, relativeTolerance, absoluteTolerance);
 603     compareMoltypes(fp, mtop1.moltype, mtop2.moltype, relativeTolerance, absoluteTolerance);
 604     compareMolblocks(fp, mtop1.molblock, mtop2.molblock);
 605     compareInteractionLists(fp, mtop1.intermolecular_ilist.get(), mtop2.intermolecular_ilist.get());
 606     compareAtomtypes(fp, mtop1.atomtypes, mtop2.atomtypes);
 607     compareAtomGroups(fp, mtop1.groups, mtop2.groups, mtop1.natoms, mtop2.natoms);
 608     compareIntermolecularExclusions(fp, mtop1.intermolecularExclusionGroup,
 609                                     mtop2.intermolecularExclusionGroup);
 610     compareBlockIndices(fp, mtop1.moleculeBlockIndices, mtop2.moleculeBlockIndices);
 611 }
 612
 613 void compareMtopAB(FILE* fp, const gmx_mtop_t& mtop1, real relativeTolerance, real absoluteTolerance)
 614 {
 615     fprintf(fp, "comparing topAB\n");
 616     compareFfparamAB(fp, mtop1.ffparams, relativeTolerance, absoluteTolerance);
 617     compareMoletypeAB(fp, mtop1.moltype, relativeTolerance, absoluteTolerance);
 618 }
 619
 620 void compareAtomGroups(FILE* fp, const SimulationGroups& g0, const SimulationGroups& g1, int natoms0, int natoms1)
 621 {
 622     fprintf(fp, "comparing groups\n");
 623
 624     for (auto group : keysOf(g0.groups))
 625     {
 626         std::string buf = gmx::formatString("grps[%d].nr", static_cast<int>(group));
 627         cmp_int(fp, buf.c_str(), -1, g0.groups[group].size(), g1.groups[group].size());
 628         if (g0.groups[group].size() == g1.groups[group].size())
 629         {
 630             for (gmx::index j = 0; j < gmx::ssize(g0.groups[group]); j++)
 631             {
 632                 buf = gmx::formatString("grps[%d].name[%zd]", static_cast<int>(group), j);
 633                 cmp_str(fp, buf.c_str(), -1, *g0.groupNames[g0.groups[group][j]],
 634                         *g1.groupNames[g1.groups[group][j]]);
 635             }
 636         }
 637         cmp_int(fp, "ngrpnr", static_cast<int>(group), g0.numberOfGroupNumbers(group),
 638                 g1.numberOfGroupNumbers(group));
 639         if (g0.numberOfGroupNumbers(group) == g1.numberOfGroupNumbers(group) && natoms0 == natoms1
 640             && (!g0.groupNumbers[group].empty() || !g1.groupNumbers[group].empty()))
 641         {
 642             for (int j = 0; j < natoms0; j++)
 643             {
 644                 cmp_int(fp, c_simulationAtomGroupTypeShortNames[group], j,
 645                         getGroupType(g0, group, j), getGroupType(g1, group, j));
 646             }
 647         }
 648     }
 649     /* We have compared the names in the groups lists,
 650      * so we can skip the grpname list comparison.
 651      */
 652 }
 653
 654 int getGroupType(const SimulationGroups& group, SimulationAtomGroupType type, int atom)
 655 {
 656     return (group.groupNumbers[type].empty() ? 0 : group.groupNumbers[type][atom]);
 657 }
 658
 659 void copy_moltype(const gmx_moltype_t* src, gmx_moltype_t* dst)
 660 {
 661     dst->name          = src->name;
 662     dst->excls         = src->excls;
 663     t_atoms* atomsCopy = copy_t_atoms(&src->atoms);
 664     dst->atoms         = *atomsCopy;
 665     sfree(atomsCopy);
 666
 667     for (int i = 0; i < F_NRE; ++i)
 668     {
 669         dst->ilist[i] = src->ilist[i];
 670     }
 671 }