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