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,2021, 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 "gromacs/utility/enumerationhelpers.h"
  41 #include "topology.h"
  42
  43 #include <cstdio>
  44
  45 #include <algorithm>
  46
  47 #include "gromacs/math/vecdump.h"
  48 #include "gromacs/topology/atoms.h"
  49 #include "gromacs/topology/idef.h"
  50 #include "gromacs/topology/ifunc.h"
  51 #include "gromacs/topology/symtab.h"
  52 #include "gromacs/utility/arrayref.h"
  53 #include "gromacs/utility/compare.h"
  54 #include "gromacs/utility/gmxassert.h"
  55 #include "gromacs/utility/smalloc.h"
  56 #include "gromacs/utility/strconvert.h"
  57 #include "gromacs/utility/txtdump.h"
  58
  59 const char* shortName(SimulationAtomGroupType type)
  60 {
  61     constexpr gmx::EnumerationArray<SimulationAtomGroupType, const char*> sc_simulationAtomGroupTypeShortNames = {
  62         { "T-Coupling",
  63           "Energy Mon.",
  64           "Acc. not used",
  65           "Freeze",
  66           "User1",
  67           "User2",
  68           "VCM",
  69           "Compressed X",
  70           "Or. Res. Fit",
  71           "QMMM" }
  72     };
  73
  74     return sc_simulationAtomGroupTypeShortNames[type];
  75 }
  76
  77 void init_top(t_topology* top)
  78 {
  79     top->name = nullptr;
  80     init_idef(&top->idef);
  81     init_atom(&(top->atoms));
  82     init_atomtypes(&(top->atomtypes));
  83     init_block(&top->mols);
  84     open_symtab(&top->symtab);
  85 }
  86
  87
  88 gmx_moltype_t::gmx_moltype_t() : name(nullptr)
  89 {
  90     init_t_atoms(&atoms, 0, FALSE);
  91 }
  92
  93 gmx_moltype_t::~gmx_moltype_t()
  94 {
  95     done_atom(&atoms);
  96 }
  97
  98 static int gmx_mtop_maxresnr(const gmx::ArrayRef<const gmx_moltype_t> moltypes, int maxres_renum)
  99 {
 100     int maxresnr = 0;
 101
 102     for (const gmx_moltype_t& moltype : moltypes)
 103     {
 104         const t_atoms& atoms = moltype.atoms;
 105         if (atoms.nres > maxres_renum)
 106         {
 107             for (int r = 0; r < atoms.nres; r++)
 108             {
 109                 if (atoms.resinfo[r].nr > maxresnr)
 110                 {
 111                     maxresnr = atoms.resinfo[r].nr;
 112                 }
 113             }
 114         }
 115     }
 116
 117     return maxresnr;
 118 }
 119
 120 gmx_mtop_t::gmx_mtop_t()
 121 {
 122     init_atomtypes(&atomtypes);
 123     open_symtab(&symtab);
 124 }
 125
 126 gmx_mtop_t::~gmx_mtop_t()
 127 {
 128     done_symtab(&symtab);
 129
 130     moltype.clear();
 131     molblock.clear();
 132     done_atomtypes(&atomtypes);
 133 }
 134
 135 void gmx_mtop_t::finalize()
 136 {
 137     if (molblock.size() == 1 && molblock[0].nmol == 1)
 138     {
 139         /* We have a single molecule only, no renumbering needed.
 140          * This case also covers an mtop converted from pdb/gro/... input,
 141          * so we retain the original residue numbering.
 142          */
 143         maxResiduesPerMoleculeToTriggerRenumber_ = 0;
 144     }
 145     else
 146     {
 147         /* We only renumber single residue molecules. Their intra-molecular
 148          * residue numbering is anyhow irrelevant.
 149          */
 150         maxResiduesPerMoleculeToTriggerRenumber_ = 1;
 151     }
 152
 153     const char* env = getenv("GMX_MAXRESRENUM");
 154     if (env != nullptr)
 155     {
 156         sscanf(env, "%d", &maxResiduesPerMoleculeToTriggerRenumber_);
 157     }
 158     if (maxResiduesPerMoleculeToTriggerRenumber_ == -1)
 159     {
 160         /* -1 signals renumber residues in all molecules */
 161         maxResiduesPerMoleculeToTriggerRenumber_ = std::numeric_limits<int>::max();
 162     }
 163
 164     maxResNumberNotRenumbered_ = gmx_mtop_maxresnr(moltype, maxResiduesPerMoleculeToTriggerRenumber_);
 165
 166     buildMolblockIndices();
 167 }
 168
 169 void gmx_mtop_t::buildMolblockIndices()
 170 {
 171     moleculeBlockIndices.resize(molblock.size());
 172
 173     int atomIndex          = 0;
 174     int residueIndex       = 0;
 175     int residueNumberStart = maxResNumberNotRenumbered_ + 1;
 176     int moleculeIndexStart = 0;
 177     for (size_t mb = 0; mb < molblock.size(); mb++)
 178     {
 179         const gmx_molblock_t& molb         = molblock[mb];
 180         MoleculeBlockIndices& indices      = moleculeBlockIndices[mb];
 181         const int             numResPerMol = moltype[molb.type].atoms.nres;
 182
 183         indices.numAtomsPerMolecule = moltype[molb.type].atoms.nr;
 184         indices.globalAtomStart     = atomIndex;
 185         indices.globalResidueStart  = residueIndex;
 186         atomIndex += molb.nmol * indices.numAtomsPerMolecule;
 187         residueIndex += molb.nmol * numResPerMol;
 188         indices.globalAtomEnd      = atomIndex;
 189         indices.residueNumberStart = residueNumberStart;
 190         if (numResPerMol <= maxResiduesPerMoleculeToTriggerRenumber_)
 191         {
 192             residueNumberStart += molb.nmol * numResPerMol;
 193         }
 194         indices.moleculeIndexStart = moleculeIndexStart;
 195         moleculeIndexStart += molb.nmol;
 196     }
 197 }
 198
 199 void done_top(t_topology* top)
 200 {
 201     done_idef(&top->idef);
 202     done_atom(&(top->atoms));
 203
 204     /* For GB */
 205     done_atomtypes(&(top->atomtypes));
 206
 207     done_symtab(&(top->symtab));
 208     done_block(&(top->mols));
 209 }
 210
 211 void done_top_mtop(t_topology* top, gmx_mtop_t* mtop)
 212 {
 213     if (mtop != nullptr)
 214     {
 215         if (top != nullptr)
 216         {
 217             done_idef(&top->idef);
 218             done_atom(&top->atoms);
 219             done_block(&top->mols);
 220             done_symtab(&top->symtab);
 221             open_symtab(&mtop->symtab);
 222             done_atomtypes(&(top->atomtypes));
 223         }
 224
 225         // Note that the rest of mtop will be freed by the destructor
 226     }
 227 }
 228
 229 gmx_localtop_t::gmx_localtop_t(const gmx_ffparams_t& ffparams) : idef(ffparams) {}
 230
 231 bool gmx_mtop_has_masses(const gmx_mtop_t* mtop)
 232 {
 233     if (mtop == nullptr)
 234     {
 235         return false;
 236     }
 237     return mtop->moltype.empty() || mtop->moltype[0].atoms.haveMass;
 238 }
 239
 240 bool gmx_mtop_has_charges(const gmx_mtop_t* mtop)
 241 {
 242     if (mtop == nullptr)
 243     {
 244         return false;
 245     }
 246     return mtop->moltype.empty() || mtop->moltype[0].atoms.haveCharge;
 247 }
 248
 249 bool gmx_mtop_has_perturbed_charges(const gmx_mtop_t& mtop)
 250 {
 251     for (const gmx_moltype_t& moltype : mtop.moltype)
 252     {
 253         const t_atoms& atoms = moltype.atoms;
 254         if (atoms.haveBState)
 255         {
 256             for (int a = 0; a < atoms.nr; a++)
 257             {
 258                 if (atoms.atom[a].q != atoms.atom[a].qB)
 259                 {
 260                     return true;
 261                 }
 262             }
 263         }
 264     }
 265     return false;
 266 }
 267
 268 bool gmx_mtop_has_atomtypes(const gmx_mtop_t* mtop)
 269 {
 270     if (mtop == nullptr)
 271     {
 272         return false;
 273     }
 274     return mtop->moltype.empty() || mtop->moltype[0].atoms.haveType;
 275 }
 276
 277 bool gmx_mtop_has_pdbinfo(const gmx_mtop_t* mtop)
 278 {
 279     if (mtop == nullptr)
 280     {
 281         return false;
 282     }
 283     return mtop->moltype.empty() || mtop->moltype[0].atoms.havePdbInfo;
 284 }
 285
 286 static void pr_grps(FILE* fp, const char* title, gmx::ArrayRef<const AtomGroupIndices> grps, char*** grpname)
 287 {
 288     int index = 0;
 289     for (const auto& group : grps)
 290     {
 291         fprintf(fp,
 292                 "%s[%-12s] nr=%zu, name=[",
 293                 title,
 294                 shortName(static_cast<SimulationAtomGroupType>(index)),
 295                 group.size());
 296         for (const auto& entry : group)
 297         {
 298             fprintf(fp, " %s", *(grpname[entry]));
 299         }
 300         fprintf(fp, "]\n");
 301         index++;
 302     }
 303 }
 304
 305 static void pr_groups(FILE* fp, int indent, const SimulationGroups& groups, gmx_bool bShowNumbers)
 306 {
 307     pr_grps(fp, "grp", groups.groups, const_cast<char***>(groups.groupNames.data()));
 308     pr_strings(fp,
 309                indent,
 310                "grpname",
 311                const_cast<char***>(groups.groupNames.data()),
 312                groups.groupNames.size(),
 313                bShowNumbers);
 314
 315     pr_indent(fp, indent);
 316     fprintf(fp, "groups          ");
 317     for (const auto group : gmx::EnumerationWrapper<SimulationAtomGroupType>{})
 318     {
 319         printf(" %5.5s", shortName(group));
 320     }
 321     printf("\n");
 322
 323     pr_indent(fp, indent);
 324     fprintf(fp, "allocated       ");
 325     int nat_max = 0;
 326     for (auto group : keysOf(groups.groups))
 327     {
 328         printf(" %5d", groups.numberOfGroupNumbers(group));
 329         nat_max = std::max(nat_max, groups.numberOfGroupNumbers(group));
 330     }
 331     printf("\n");
 332
 333     if (nat_max == 0)
 334     {
 335         pr_indent(fp, indent);
 336         fprintf(fp, "groupnr[%5s] =", "*");
 337         for (auto gmx_unused group : keysOf(groups.groups))
 338         {
 339             fprintf(fp, "  %3d ", 0);
 340         }
 341         fprintf(fp, "\n");
 342     }
 343     else
 344     {
 345         for (int i = 0; i < nat_max; i++)
 346         {
 347             pr_indent(fp, indent);
 348             fprintf(fp, "groupnr[%5d] =", i);
 349             for (auto group : keysOf(groups.groups))
 350             {
 351                 fprintf(fp,
 352                         "  %3d ",
 353                         !groups.groupNumbers[group].empty() ? groups.groupNumbers[group][i] : 0);
 354             }
 355             fprintf(fp, "\n");
 356         }
 357     }
 358 }
 359
 360 static void pr_moltype(FILE*                 fp,
 361                        int                   indent,
 362                        const char*           title,
 363                        const gmx_moltype_t*  molt,
 364                        int                   n,
 365                        const gmx_ffparams_t* ffparams,
 366                        gmx_bool              bShowNumbers,
 367                        gmx_bool              bShowParameters)
 368 {
 369     indent = pr_title_n(fp, indent, title, n);
 370     pr_indent(fp, indent);
 371     fprintf(fp, "name=\"%s\"\n", *(molt->name));
 372     pr_atoms(fp, indent, "atoms", &(molt->atoms), bShowNumbers);
 373     pr_listoflists(fp, indent, "excls", &molt->excls, bShowNumbers);
 374     for (int j = 0; (j < F_NRE); j++)
 375     {
 376         pr_ilist(fp,
 377                  indent,
 378                  interaction_function[j].longname,
 379                  ffparams->functype.data(),
 380                  molt->ilist[j],
 381                  bShowNumbers,
 382                  bShowParameters,
 383                  ffparams->iparams.data());
 384     }
 385 }
 386
 387 static void pr_molblock(FILE*                             fp,
 388                         int                               indent,
 389                         const char*                       title,
 390                         const gmx_molblock_t*             molb,
 391                         int                               n,
 392                         const std::vector<gmx_moltype_t>& molt)
 393 {
 394     indent = pr_title_n(fp, indent, title, n);
 395     pr_indent(fp, indent);
 396     fprintf(fp, "%-20s = %d \"%s\"\n", "moltype", molb->type, *(molt[molb->type].name));
 397     pr_int(fp, indent, "#molecules", molb->nmol);
 398     pr_int(fp, indent, "#posres_xA", molb->posres_xA.size());
 399     if (!molb->posres_xA.empty())
 400     {
 401         pr_rvecs(fp, indent, "posres_xA", as_rvec_array(molb->posres_xA.data()), molb->posres_xA.size());
 402     }
 403     pr_int(fp, indent, "#posres_xB", molb->posres_xB.size());
 404     if (!molb->posres_xB.empty())
 405     {
 406         pr_rvecs(fp, indent, "posres_xB", as_rvec_array(molb->posres_xB.data()), molb->posres_xB.size());
 407     }
 408 }
 409
 410 void pr_mtop(FILE* fp, int indent, const char* title, const gmx_mtop_t* mtop, gmx_bool bShowNumbers, gmx_bool bShowParameters)
 411 {
 412     if (available(fp, mtop, indent, title))
 413     {
 414         indent = pr_title(fp, indent, title);
 415         pr_indent(fp, indent);
 416         fprintf(fp, "name=\"%s\"\n", *(mtop->name));
 417         pr_int(fp, indent, "#atoms", mtop->natoms);
 418         pr_int(fp, indent, "#molblock", mtop->molblock.size());
 419         for (size_t mb = 0; mb < mtop->molblock.size(); mb++)
 420         {
 421             pr_molblock(fp, indent, "molblock", &mtop->molblock[mb], mb, mtop->moltype);
 422         }
 423         pr_str(fp, indent, "bIntermolecularInteractions", gmx::boolToString(mtop->bIntermolecularInteractions));
 424         if (mtop->bIntermolecularInteractions)
 425         {
 426             for (int j = 0; j < F_NRE; j++)
 427             {
 428                 pr_ilist(fp,
 429                          indent,
 430                          interaction_function[j].longname,
 431                          mtop->ffparams.functype.data(),
 432                          (*mtop->intermolecular_ilist)[j],
 433                          bShowNumbers,
 434                          bShowParameters,
 435                          mtop->ffparams.iparams.data());
 436             }
 437         }
 438         pr_ffparams(fp, indent, "ffparams", &(mtop->ffparams), bShowNumbers);
 439         pr_atomtypes(fp, indent, "atomtypes", &(mtop->atomtypes), bShowNumbers);
 440         for (size_t mt = 0; mt < mtop->moltype.size(); mt++)
 441         {
 442             pr_moltype(fp, indent, "moltype", &mtop->moltype[mt], mt, &mtop->ffparams, bShowNumbers, bShowParameters);
 443         }
 444         pr_groups(fp, indent, mtop->groups, bShowNumbers);
 445     }
 446 }
 447
 448 void pr_top(FILE* fp, int indent, const char* title, const t_topology* top, gmx_bool bShowNumbers, gmx_bool bShowParameters)
 449 {
 450     if (available(fp, top, indent, title))
 451     {
 452         indent = pr_title(fp, indent, title);
 453         pr_indent(fp, indent);
 454         fprintf(fp, "name=\"%s\"\n", *(top->name));
 455         pr_atoms(fp, indent, "atoms", &(top->atoms), bShowNumbers);
 456         pr_atomtypes(fp, indent, "atomtypes", &(top->atomtypes), bShowNumbers);
 457         pr_block(fp, indent, "mols", &top->mols, bShowNumbers);
 458         pr_str(fp, indent, "bIntermolecularInteractions", gmx::boolToString(top->bIntermolecularInteractions));
 459         pr_idef(fp, indent, "idef", &top->idef, bShowNumbers, bShowParameters);
 460     }
 461 }
 462
 463 static void cmp_iparm(FILE*            fp,
 464                       const char*      s,
 465                       t_functype       ft,
 466                       const t_iparams& ip1,
 467                       const t_iparams& ip2,
 468                       real             relativeTolerance,
 469                       real             absoluteTolerance)
 470 {
 471     bool bDiff = false;
 472     for (int i = 0; i < MAXFORCEPARAM && !bDiff; i++)
 473     {
 474         bDiff = !equal_real(ip1.generic.buf[i], ip2.generic.buf[i], relativeTolerance, absoluteTolerance);
 475     }
 476     if (bDiff)
 477     {
 478         fprintf(fp, "%s1: ", s);
 479         pr_iparams(fp, ft, ip1);
 480         fprintf(fp, "%s2: ", s);
 481         pr_iparams(fp, ft, ip2);
 482     }
 483 }
 484
 485 static void cmp_iparm_AB(FILE* fp, const char* s, t_functype ft, const t_iparams& ip1, real relativeTolerance, real absoluteTolerance)
 486 {
 487     /* Normally the first parameter is perturbable */
 488     int p0    = 0;
 489     int nrfpA = interaction_function[ft].nrfpA;
 490     int nrfpB = interaction_function[ft].nrfpB;
 491     if (ft == F_PDIHS)
 492     {
 493         nrfpB = 2;
 494     }
 495     else if (interaction_function[ft].flags & IF_TABULATED)
 496     {
 497         /* For tabulated interactions only the second parameter is perturbable */
 498         p0    = 1;
 499         nrfpB = 1;
 500     }
 501     bool bDiff = false;
 502     for (int i = 0; i < nrfpB && !bDiff; i++)
 503     {
 504         bDiff = !equal_real(
 505                 ip1.generic.buf[p0 + i], ip1.generic.buf[nrfpA + i], relativeTolerance, absoluteTolerance);
 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 relativeTolerance, real absoluteTolerance)
 515 {
 516     int cmap1_ngrid = (cmap1 ? cmap1->cmapdata.size() : 0);
 517     int cmap2_ngrid = (cmap2 ? cmap2->cmapdata.size() : 0);
 518
 519     cmp_int(fp, "cmap ngrid", -1, cmap1_ngrid, cmap2_ngrid);
 520
 521     if (cmap1 == nullptr || cmap2 == nullptr)
 522     {
 523         return;
 524     }
 525
 526     cmp_int(fp, "cmap grid_spacing", -1, cmap1->grid_spacing, cmap2->grid_spacing);
 527     if (cmap1->cmapdata.size() == cmap2->cmapdata.size() && cmap1->grid_spacing == cmap2->grid_spacing)
 528     {
 529         for (size_t g = 0; g < cmap1->cmapdata.size(); g++)
 530         {
 531             fprintf(fp, "comparing cmap %zu\n", g);
 532
 533             for (int i = 0; i < 4 * cmap1->grid_spacing * cmap1->grid_spacing; i++)
 534             {
 535                 cmp_real(fp, "", i, cmap1->cmapdata[g].cmap[i], cmap2->cmapdata[g].cmap[i], relativeTolerance, absoluteTolerance);
 536             }
 537         }
 538     }
 539 }
 540
 541 static void cmp_listoflists(FILE*                        fp,
 542                             const gmx::ListOfLists<int>& list1,
 543                             const gmx::ListOfLists<int>& list2,
 544                             const char*                  s)
 545 {
 546     char buf[32];
 547
 548     fprintf(fp, "comparing blocka %s\n", s);
 549     sprintf(buf, "%s.numLists", s);
 550     cmp_int(fp, buf, -1, list1.ssize(), list2.ssize());
 551     sprintf(buf, "%s.numElements", s);
 552     cmp_int(fp, buf, -1, list1.numElements(), list2.numElements());
 553 }
 554
 555 static void compareFfparams(FILE*                 fp,
 556                             const gmx_ffparams_t& ff1,
 557                             const gmx_ffparams_t& ff2,
 558                             real                  relativeTolerance,
 559                             real                  absoluteTolerance)
 560 {
 561     fprintf(fp, "comparing force field parameters\n");
 562     cmp_int(fp, "numTypes", -1, ff1.numTypes(), ff2.numTypes());
 563     cmp_int(fp, "atnr", -1, ff1.atnr, ff1.atnr);
 564     cmp_double(fp, "reppow", -1, ff1.reppow, ff2.reppow, relativeTolerance, absoluteTolerance);
 565     cmp_real(fp, "fudgeQQ", -1, ff1.fudgeQQ, ff2.fudgeQQ, relativeTolerance, absoluteTolerance);
 566     cmp_cmap(fp, &ff1.cmap_grid, &ff2.cmap_grid, relativeTolerance, absoluteTolerance);
 567     for (int i = 0; i < std::min(ff1.numTypes(), ff2.numTypes()); i++)
 568     {
 569         std::string buf = gmx::formatString("ffparams->functype[%d]", i);
 570         cmp_int(fp, buf.c_str(), i, ff1.functype[i], ff2.functype[i]);
 571         buf = gmx::formatString("ffparams->iparams[%d]", i);
 572         cmp_iparm(fp, buf.c_str(), ff1.functype[i], ff1.iparams[i], ff2.iparams[i], relativeTolerance, absoluteTolerance);
 573     }
 574 }
 575
 576 static void compareFfparamAB(FILE* fp, const gmx_ffparams_t& ff1, real relativeTolerance, real absoluteTolerance)
 577 {
 578     fprintf(fp, "comparing free energy parameters\n");
 579     for (int i = 0; i < ff1.numTypes(); i++)
 580     {
 581         std::string buf = gmx::formatString("ffparams->iparams[%d]", i);
 582         cmp_iparm_AB(fp, buf.c_str(), ff1.functype[i], ff1.iparams[i], relativeTolerance, absoluteTolerance);
 583     }
 584 }
 585 static void compareInteractionLists(FILE* fp, const InteractionLists* il1, const InteractionLists* il2)
 586 {
 587     fprintf(fp, "comparing InteractionLists\n");
 588     if ((il1 || il2) && (!il1 || !il2))
 589     {
 590         fprintf(fp, "InteractionLists are present in topology %d but not in the other\n", il1 ? 1 : 2);
 591     }
 592     if (il1 && il2)
 593     {
 594         for (int i = 0; i < F_NRE; i++)
 595         {
 596             cmp_int(fp, "InteractionList size", i, il1->at(i).size(), il2->at(i).size());
 597             int nr = std::min(il1->at(i).size(), il2->at(i).size());
 598             for (int j = 0; j < nr; j++)
 599             {
 600                 cmp_int(fp, "InteractionList entry", j, il1->at(i).iatoms.at(j), il2->at(i).iatoms.at(j));
 601             }
 602         }
 603     }
 604 }
 605
 606 static void compareMoltypes(FILE*                              fp,
 607                             gmx::ArrayRef<const gmx_moltype_t> mt1,
 608                             gmx::ArrayRef<const gmx_moltype_t> mt2,
 609                             real                               relativeTolerance,
 610                             real                               absoluteTolerance)
 611 {
 612     fprintf(fp, "comparing molecule types\n");
 613     cmp_int(fp, "moltype size", -1, mt1.size(), mt2.size());
 614     for (int i = 0; i < std::min(mt1.ssize(), mt2.ssize()); i++)
 615     {
 616         cmp_str(fp, "Name", i, *mt1[i].name, *mt2[i].name);
 617         compareAtoms(fp, &mt1[i].atoms, &mt2[i].atoms, relativeTolerance, absoluteTolerance);
 618         compareInteractionLists(fp, &mt1[i].ilist, &mt2[i].ilist);
 619         std::string buf = gmx::formatString("excls[%d]", i);
 620         cmp_listoflists(fp, mt1[i].excls, mt2[i].excls, buf.c_str());
 621     }
 622 }
 623
 624 static void compareMoletypeAB(FILE* fp, gmx::ArrayRef<const gmx_moltype_t> mt1, real relativeTolerance, real absoluteTolerance)
 625 {
 626     fprintf(fp, "comparing free energy molecule types\n");
 627     for (gmx::index i = 0; i < mt1.ssize(); i++)
 628     {
 629         compareAtoms(fp, &mt1[i].atoms, nullptr, relativeTolerance, absoluteTolerance);
 630     }
 631 }
 632 static void compareMolblocks(FILE*                               fp,
 633                              gmx::ArrayRef<const gmx_molblock_t> mb1,
 634                              gmx::ArrayRef<const gmx_molblock_t> mb2)
 635 {
 636     fprintf(fp, "comparing molecule blocks\n");
 637     cmp_int(fp, "molblock size", -1, mb1.size(), mb2.size());
 638     int nr = std::min(mb1.size(), mb2.size());
 639     for (int i = 0; i < nr; i++)
 640     {
 641         cmp_int(fp, "type", i, mb1[i].type, mb2[i].type);
 642         cmp_int(fp, "nmol", i, mb1[i].nmol, mb2[i].nmol);
 643         // Only checking size of restraint vectors for now
 644         cmp_int(fp, "posres_xA size", i, mb1[i].posres_xA.size(), mb2[i].posres_xA.size());
 645         cmp_int(fp, "posres_xB size", i, mb1[i].posres_xB.size(), mb2[i].posres_xB.size());
 646     }
 647 }
 648
 649 static void compareAtomtypes(FILE* fp, const t_atomtypes& at1, const t_atomtypes& at2)
 650 {
 651     fprintf(fp, "comparing atomtypes\n");
 652     cmp_int(fp, "nr", -1, at1.nr, at2.nr);
 653     int nr = std::min(at1.nr, at2.nr);
 654     for (int i = 0; i < nr; i++)
 655     {
 656         cmp_int(fp, "atomtype", i, at1.atomnumber[i], at2.atomnumber[i]);
 657     }
 658 }
 659
 660 static void compareIntermolecularExclusions(FILE*                    fp,
 661                                             gmx::ArrayRef<const int> ime1,
 662                                             gmx::ArrayRef<const int> ime2)
 663 {
 664     fprintf(fp, "comparing intermolecular exclusions\n");
 665     cmp_int(fp, "exclusion number", -1, ime1.size(), ime2.size());
 666     int nr = std::min(ime1.size(), ime2.size());
 667     for (int i = 0; i < nr; i++)
 668     {
 669         cmp_int(fp, "exclusion", i, ime1[i], ime2[i]);
 670     }
 671 }
 672
 673 static void compareBlockIndices(FILE*                                     fp,
 674                                 gmx::ArrayRef<const MoleculeBlockIndices> mbi1,
 675                                 gmx::ArrayRef<const MoleculeBlockIndices> mbi2)
 676 {
 677     fprintf(fp, "comparing moleculeBlockIndices\n");
 678     cmp_int(fp, "size", -1, mbi1.size(), mbi2.size());
 679     int nr = std::min(mbi1.size(), mbi2.size());
 680     for (int i = 0; i < nr; i++)
 681     {
 682         cmp_int(fp, "numAtomsPerMolecule", i, mbi1[i].numAtomsPerMolecule, mbi2[i].numAtomsPerMolecule);
 683         cmp_int(fp, "globalAtomStart", i, mbi1[i].globalAtomStart, mbi2[i].globalAtomStart);
 684         cmp_int(fp, "globalAtomEnd", i, mbi1[i].globalAtomEnd, mbi2[i].globalAtomEnd);
 685         cmp_int(fp, "globalResidueStart", i, mbi1[i].globalResidueStart, mbi2[i].globalResidueStart);
 686         cmp_int(fp, "moleculeIndexStart", i, mbi1[i].moleculeIndexStart, mbi2[i].moleculeIndexStart);
 687     }
 688 }
 689
 690 void compareMtop(FILE* fp, const gmx_mtop_t& mtop1, const gmx_mtop_t& mtop2, real relativeTolerance, real absoluteTolerance)
 691 {
 692     fprintf(fp, "comparing mtop topology\n");
 693     cmp_str(fp, "Name", -1, *mtop1.name, *mtop2.name);
 694     cmp_int(fp, "natoms", -1, mtop1.natoms, mtop2.natoms);
 695     cmp_int(fp,
 696             "maxres_renum",
 697             -1,
 698             mtop1.maxResiduesPerMoleculeToTriggerRenumber(),
 699             mtop2.maxResiduesPerMoleculeToTriggerRenumber());
 700     cmp_int(fp, "maxresnr", -1, mtop1.maxResNumberNotRenumbered(), mtop2.maxResNumberNotRenumbered());
 701     cmp_bool(fp, "bIntermolecularInteractions", -1, mtop1.bIntermolecularInteractions, mtop2.bIntermolecularInteractions);
 702     cmp_bool(fp, "haveMoleculeIndices", -1, mtop1.haveMoleculeIndices, mtop2.haveMoleculeIndices);
 703
 704     compareFfparams(fp, mtop1.ffparams, mtop2.ffparams, relativeTolerance, absoluteTolerance);
 705     compareMoltypes(fp, mtop1.moltype, mtop2.moltype, relativeTolerance, absoluteTolerance);
 706     compareMolblocks(fp, mtop1.molblock, mtop2.molblock);
 707     compareInteractionLists(fp, mtop1.intermolecular_ilist.get(), mtop2.intermolecular_ilist.get());
 708     compareAtomtypes(fp, mtop1.atomtypes, mtop2.atomtypes);
 709     compareAtomGroups(fp, mtop1.groups, mtop2.groups, mtop1.natoms, mtop2.natoms);
 710     compareIntermolecularExclusions(
 711             fp, mtop1.intermolecularExclusionGroup, mtop2.intermolecularExclusionGroup);
 712     compareBlockIndices(fp, mtop1.moleculeBlockIndices, mtop2.moleculeBlockIndices);
 713 }
 714
 715 void compareMtopAB(FILE* fp, const gmx_mtop_t& mtop1, real relativeTolerance, real absoluteTolerance)
 716 {
 717     fprintf(fp, "comparing topAB\n");
 718     compareFfparamAB(fp, mtop1.ffparams, relativeTolerance, absoluteTolerance);
 719     compareMoletypeAB(fp, mtop1.moltype, relativeTolerance, absoluteTolerance);
 720 }
 721
 722 void compareAtomGroups(FILE* fp, const SimulationGroups& g0, const SimulationGroups& g1, int natoms0, int natoms1)
 723 {
 724     fprintf(fp, "comparing groups\n");
 725
 726     for (auto group : keysOf(g0.groups))
 727     {
 728         std::string buf = gmx::formatString("grps[%d].nr", static_cast<int>(group));
 729         cmp_int(fp, buf.c_str(), -1, g0.groups[group].size(), g1.groups[group].size());
 730         if (g0.groups[group].size() == g1.groups[group].size())
 731         {
 732             for (gmx::index j = 0; j < gmx::ssize(g0.groups[group]); j++)
 733             {
 734                 buf = gmx::formatString("grps[%d].name[%zd]", static_cast<int>(group), j);
 735                 cmp_str(fp,
 736                         buf.c_str(),
 737                         -1,
 738                         *g0.groupNames[g0.groups[group][j]],
 739                         *g1.groupNames[g1.groups[group][j]]);
 740             }
 741         }
 742         cmp_int(fp,
 743                 "ngrpnr",
 744                 static_cast<int>(group),
 745                 g0.numberOfGroupNumbers(group),
 746                 g1.numberOfGroupNumbers(group));
 747         if (g0.numberOfGroupNumbers(group) == g1.numberOfGroupNumbers(group) && natoms0 == natoms1
 748             && (!g0.groupNumbers[group].empty() || !g1.groupNumbers[group].empty()))
 749         {
 750             for (int j = 0; j < natoms0; j++)
 751             {
 752                 cmp_int(fp, shortName(group), j, getGroupType(g0, group, j), getGroupType(g1, group, j));
 753             }
 754         }
 755     }
 756     /* We have compared the names in the groups lists,
 757      * so we can skip the grpname list comparison.
 758      */
 759 }
 760
 761 int getGroupType(const SimulationGroups& group, SimulationAtomGroupType type, int atom)
 762 {
 763     return (group.groupNumbers[type].empty() ? 0 : group.groupNumbers[type][atom]);
 764 }
 765
 766 void copy_moltype(const gmx_moltype_t* src, gmx_moltype_t* dst)
 767 {
 768     dst->name          = src->name;
 769     dst->excls         = src->excls;
 770     t_atoms* atomsCopy = copy_t_atoms(&src->atoms);
 771     dst->atoms         = *atomsCopy;
 772     sfree(atomsCopy);
 773
 774     for (int i = 0; i < F_NRE; ++i)
 775     {
 776         dst->ilist[i] = src->ilist[i];
 777     }
 778 }