src/gromacs/topology/mtop_util.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2008,2009,2010, The GROMACS development team.
   5  * Copyright (c) 2012,2013,2014,2015,2016 The GROMACS development team.
   6  * Copyright (c) 2017,2018,2019,2020,2021, 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 "mtop_util.h"
  40
  41 #include <climits>
  42 #include <cstddef>
  43 #include <cstdio>
  44 #include <cstdlib>
  45 #include <cstring>
  46
  47 #include "gromacs/math/vectypes.h"
  48 #include "gromacs/topology/atoms.h"
  49 #include "gromacs/topology/block.h"
  50 #include "gromacs/topology/exclusionblocks.h"
  51 #include "gromacs/topology/idef.h"
  52 #include "gromacs/topology/ifunc.h"
  53 #include "gromacs/topology/topology.h"
  54 #include "gromacs/topology/topsort.h"
  55 #include "gromacs/utility/arrayref.h"
  56 #include "gromacs/utility/enumerationhelpers.h"
  57 #include "gromacs/utility/fatalerror.h"
  58 #include "gromacs/utility/real.h"
  59 #include "gromacs/utility/smalloc.h"
  60
  61 void gmx_mtop_count_atomtypes(const gmx_mtop_t& mtop, int state, int typecount[])
  62 {
  63     for (int i = 0; i < mtop.ffparams.atnr; ++i)
  64     {
  65         typecount[i] = 0;
  66     }
  67     for (const gmx_molblock_t& molb : mtop.molblock)
  68     {
  69         const t_atoms& atoms = mtop.moltype[molb.type].atoms;
  70         for (int i = 0; i < atoms.nr; ++i)
  71         {
  72             const int tpi = (state == 0) ? atoms.atom[i].type : atoms.atom[i].typeB;
  73             typecount[tpi] += molb.nmol;
  74         }
  75     }
  76 }
  77
  78 int gmx_mtop_num_molecules(const gmx_mtop_t& mtop)
  79 {
  80     int numMolecules = 0;
  81     for (const gmx_molblock_t& molb : mtop.molblock)
  82     {
  83         numMolecules += molb.nmol;
  84     }
  85     return numMolecules;
  86 }
  87
  88 int gmx_mtop_nres(const gmx_mtop_t& mtop)
  89 {
  90     int nres = 0;
  91     for (const gmx_molblock_t& molb : mtop.molblock)
  92     {
  93         nres += molb.nmol * mtop.moltype[molb.type].atoms.nres;
  94     }
  95     return nres;
  96 }
  97
  98 AtomIterator::AtomIterator(const gmx_mtop_t& mtop, int globalAtomNumber) :
  99     mtop_(&mtop),
 100     mblock_(0),
 101     atoms_(&mtop.moltype[mtop.molblock[0].type].atoms),
 102     currentMolecule_(0),
 103     highestResidueNumber_(mtop.maxResNumberNotRenumbered()),
 104     localAtomNumber_(0),
 105     globalAtomNumber_(globalAtomNumber)
 106 {
 107     GMX_ASSERT(globalAtomNumber == 0 || globalAtomNumber == mtop.natoms,
 108                "Starting at other atoms not implemented yet");
 109 }
 110
 111 AtomIterator& AtomIterator::operator++()
 112 {
 113     localAtomNumber_++;
 114     globalAtomNumber_++;
 115
 116     if (localAtomNumber_ >= atoms_->nr)
 117     {
 118         if (atoms_->nres <= mtop_->maxResiduesPerMoleculeToTriggerRenumber())
 119         {
 120             /* Single residue molecule, increase the count with one */
 121             highestResidueNumber_ += atoms_->nres;
 122         }
 123         currentMolecule_++;
 124         localAtomNumber_ = 0;
 125         if (currentMolecule_ >= mtop_->molblock[mblock_].nmol)
 126         {
 127             mblock_++;
 128             if (mblock_ >= mtop_->molblock.size())
 129             {
 130                 return *this;
 131             }
 132             atoms_           = &mtop_->moltype[mtop_->molblock[mblock_].type].atoms;
 133             currentMolecule_ = 0;
 134         }
 135     }
 136     return *this;
 137 }
 138
 139 bool AtomIterator::operator==(const AtomIterator& o) const
 140 {
 141     return mtop_ == o.mtop_ && globalAtomNumber_ == o.globalAtomNumber_;
 142 }
 143
 144 const t_atom& AtomProxy::atom() const
 145 {
 146     return it_->atoms_->atom[it_->localAtomNumber_];
 147 }
 148
 149 int AtomProxy::globalAtomNumber() const
 150 {
 151     return it_->globalAtomNumber_;
 152 }
 153
 154 const char* AtomProxy::atomName() const
 155 {
 156     return *(it_->atoms_->atomname[it_->localAtomNumber_]);
 157 }
 158
 159 const char* AtomProxy::residueName() const
 160 {
 161     int residueIndexInMolecule = it_->atoms_->atom[it_->localAtomNumber_].resind;
 162     return *(it_->atoms_->resinfo[residueIndexInMolecule].name);
 163 }
 164
 165 int AtomProxy::residueNumber() const
 166 {
 167     int residueIndexInMolecule = it_->atoms_->atom[it_->localAtomNumber_].resind;
 168     if (it_->atoms_->nres <= it_->mtop_->maxResiduesPerMoleculeToTriggerRenumber())
 169     {
 170         return it_->highestResidueNumber_ + 1 + residueIndexInMolecule;
 171     }
 172     else
 173     {
 174         return it_->atoms_->resinfo[residueIndexInMolecule].nr;
 175     }
 176 }
 177
 178 const gmx_moltype_t& AtomProxy::moleculeType() const
 179 {
 180     return it_->mtop_->moltype[it_->mtop_->molblock[it_->mblock_].type];
 181 }
 182
 183 int AtomProxy::atomNumberInMol() const
 184 {
 185     return it_->localAtomNumber_;
 186 }
 187
 188 struct gmx_mtop_atomloop_block
 189 {
 190     const gmx_mtop_t* mtop;
 191     size_t            mblock;
 192     const t_atoms*    atoms;
 193     int               at_local;
 194 };
 195
 196 gmx_mtop_atomloop_block_t gmx_mtop_atomloop_block_init(const gmx_mtop_t& mtop)
 197 {
 198     struct gmx_mtop_atomloop_block* aloop = nullptr;
 199
 200     snew(aloop, 1);
 201
 202     aloop->mtop     = &mtop;
 203     aloop->mblock   = 0;
 204     aloop->atoms    = &mtop.moltype[mtop.molblock[aloop->mblock].type].atoms;
 205     aloop->at_local = -1;
 206
 207     return aloop;
 208 }
 209
 210 static void gmx_mtop_atomloop_block_destroy(gmx_mtop_atomloop_block_t aloop)
 211 {
 212     sfree(aloop);
 213 }
 214
 215 gmx_bool gmx_mtop_atomloop_block_next(gmx_mtop_atomloop_block_t aloop, const t_atom** atom, int* nmol)
 216 {
 217     if (aloop == nullptr)
 218     {
 219         gmx_incons("gmx_mtop_atomloop_all_next called without calling gmx_mtop_atomloop_all_init");
 220     }
 221
 222     aloop->at_local++;
 223
 224     if (aloop->at_local >= aloop->atoms->nr)
 225     {
 226         aloop->mblock++;
 227         if (aloop->mblock >= aloop->mtop->molblock.size())
 228         {
 229             gmx_mtop_atomloop_block_destroy(aloop);
 230             return FALSE;
 231         }
 232         aloop->atoms    = &aloop->mtop->moltype[aloop->mtop->molblock[aloop->mblock].type].atoms;
 233         aloop->at_local = 0;
 234     }
 235
 236     *atom = &aloop->atoms->atom[aloop->at_local];
 237     *nmol = aloop->mtop->molblock[aloop->mblock].nmol;
 238
 239     return TRUE;
 240 }
 241
 242 IListIterator::IListIterator(const gmx_mtop_t& mtop, size_t molblock) :
 243     mtop_(&mtop), mblock_(molblock)
 244 {
 245 }
 246
 247 IListIterator& IListIterator::operator++()
 248 {
 249     mblock_++;
 250     return *this;
 251 }
 252
 253 bool IListIterator::operator==(const IListIterator& o) const
 254 {
 255     return mtop_ == o.mtop_ && mblock_ == o.mblock_;
 256 }
 257
 258 const InteractionLists& IListProxy::list() const
 259 {
 260     // one past the end means we want to take the
 261     // intermolecular list instead.
 262     if (it_->mblock_ == it_->mtop_->molblock.size())
 263     {
 264         return *it_->mtop_->intermolecular_ilist;
 265     }
 266     else
 267     {
 268         return it_->mtop_->moltype[it_->mtop_->molblock[it_->mblock_].type].ilist;
 269     }
 270 }
 271
 272 int IListProxy::nmol() const
 273 {
 274     // one past the end means we want to take the
 275     // intermolecular list instead.
 276     if (it_->mblock_ == it_->mtop_->molblock.size())
 277     {
 278         return 1;
 279     }
 280     else
 281     {
 282         return it_->mtop_->molblock[it_->mblock_].nmol;
 283     }
 284 }
 285
 286 IListRange::IListRange(const gmx_mtop_t& mtop) : begin_(mtop), end_(mtop, mtop.molblock.size())
 287 {
 288     if (mtop.bIntermolecularInteractions)
 289     {
 290         end_ = IListIterator(mtop, mtop.molblock.size() + 1);
 291     }
 292 }
 293
 294 int gmx_mtop_ftype_count(const gmx_mtop_t& mtop, int ftype)
 295 {
 296     int n = 0;
 297
 298     for (const IListProxy il : IListRange(mtop))
 299     {
 300         n += il.nmol() * il.list()[ftype].size() / (1 + NRAL(ftype));
 301     }
 302
 303     if (mtop.bIntermolecularInteractions)
 304     {
 305         n += (*mtop.intermolecular_ilist)[ftype].size() / (1 + NRAL(ftype));
 306     }
 307
 308     return n;
 309 }
 310
 311 int gmx_mtop_interaction_count(const gmx_mtop_t& mtop, const int unsigned if_flags)
 312 {
 313     int n = 0;
 314
 315     for (const IListProxy il : IListRange(mtop))
 316     {
 317         for (int ftype = 0; ftype < F_NRE; ftype++)
 318         {
 319             if ((interaction_function[ftype].flags & if_flags) == if_flags)
 320             {
 321                 n += il.nmol() * il.list()[ftype].size() / (1 + NRAL(ftype));
 322             }
 323         }
 324     }
 325
 326     if (mtop.bIntermolecularInteractions)
 327     {
 328         for (int ftype = 0; ftype < F_NRE; ftype++)
 329         {
 330             if ((interaction_function[ftype].flags & if_flags) == if_flags)
 331             {
 332                 n += (*mtop.intermolecular_ilist)[ftype].size() / (1 + NRAL(ftype));
 333             }
 334         }
 335     }
 336
 337     return n;
 338 }
 339
 340 gmx::EnumerationArray<ParticleType, int> gmx_mtop_particletype_count(const gmx_mtop_t& mtop)
 341 {
 342     gmx::EnumerationArray<ParticleType, int> count = { { 0 } };
 343
 344     for (const auto& molblock : mtop.molblock)
 345     {
 346         const t_atoms& atoms = mtop.moltype[molblock.type].atoms;
 347         for (int a = 0; a < atoms.nr; a++)
 348         {
 349             count[atoms.atom[a].ptype] += molblock.nmol;
 350         }
 351     }
 352
 353     return count;
 354 }
 355
 356 static void atomcat(t_atoms* dest, const t_atoms* src, int copies, int maxres_renum, int* maxresnr)
 357 {
 358     int i = 0, j = 0, l = 0, size = 0;
 359     int srcnr  = src->nr;
 360     int destnr = dest->nr;
 361
 362     if (dest->nr == 0)
 363     {
 364         dest->haveMass    = src->haveMass;
 365         dest->haveType    = src->haveType;
 366         dest->haveCharge  = src->haveCharge;
 367         dest->haveBState  = src->haveBState;
 368         dest->havePdbInfo = src->havePdbInfo;
 369     }
 370     else
 371     {
 372         dest->haveMass    = dest->haveMass && src->haveMass;
 373         dest->haveType    = dest->haveType && src->haveType;
 374         dest->haveCharge  = dest->haveCharge && src->haveCharge;
 375         dest->haveBState  = dest->haveBState && src->haveBState;
 376         dest->havePdbInfo = dest->havePdbInfo && src->havePdbInfo;
 377     }
 378
 379     if (srcnr)
 380     {
 381         size = destnr + copies * srcnr;
 382         srenew(dest->atom, size);
 383         srenew(dest->atomname, size);
 384         if (dest->haveType)
 385         {
 386             srenew(dest->atomtype, size);
 387             if (dest->haveBState)
 388             {
 389                 srenew(dest->atomtypeB, size);
 390             }
 391         }
 392         if (dest->havePdbInfo)
 393         {
 394             srenew(dest->pdbinfo, size);
 395         }
 396     }
 397     if (src->nres)
 398     {
 399         size = dest->nres + copies * src->nres;
 400         srenew(dest->resinfo, size);
 401     }
 402
 403     /* residue information */
 404     for (l = dest->nres, j = 0; (j < copies); j++, l += src->nres)
 405     {
 406         memcpy(reinterpret_cast<char*>(&(dest->resinfo[l])),
 407                reinterpret_cast<char*>(&(src->resinfo[0])),
 408                static_cast<size_t>(src->nres * sizeof(src->resinfo[0])));
 409     }
 410
 411     for (l = destnr, j = 0; (j < copies); j++, l += srcnr)
 412     {
 413         memcpy(reinterpret_cast<char*>(&(dest->atom[l])),
 414                reinterpret_cast<char*>(&(src->atom[0])),
 415                static_cast<size_t>(srcnr * sizeof(src->atom[0])));
 416         memcpy(reinterpret_cast<char*>(&(dest->atomname[l])),
 417                reinterpret_cast<char*>(&(src->atomname[0])),
 418                static_cast<size_t>(srcnr * sizeof(src->atomname[0])));
 419         if (dest->haveType)
 420         {
 421             memcpy(reinterpret_cast<char*>(&(dest->atomtype[l])),
 422                    reinterpret_cast<char*>(&(src->atomtype[0])),
 423                    static_cast<size_t>(srcnr * sizeof(src->atomtype[0])));
 424             if (dest->haveBState)
 425             {
 426                 memcpy(reinterpret_cast<char*>(&(dest->atomtypeB[l])),
 427                        reinterpret_cast<char*>(&(src->atomtypeB[0])),
 428                        static_cast<size_t>(srcnr * sizeof(src->atomtypeB[0])));
 429             }
 430         }
 431         if (dest->havePdbInfo)
 432         {
 433             memcpy(reinterpret_cast<char*>(&(dest->pdbinfo[l])),
 434                    reinterpret_cast<char*>(&(src->pdbinfo[0])),
 435                    static_cast<size_t>(srcnr * sizeof(src->pdbinfo[0])));
 436         }
 437     }
 438
 439     /* Increment residue indices */
 440     for (l = destnr, j = 0; (j < copies); j++)
 441     {
 442         for (i = 0; (i < srcnr); i++, l++)
 443         {
 444             dest->atom[l].resind = dest->nres + j * src->nres + src->atom[i].resind;
 445         }
 446     }
 447
 448     if (src->nres <= maxres_renum)
 449     {
 450         /* Single residue molecule, continue counting residues */
 451         for (j = 0; (j < copies); j++)
 452         {
 453             for (l = 0; l < src->nres; l++)
 454             {
 455                 (*maxresnr)++;
 456                 dest->resinfo[dest->nres + j * src->nres + l].nr = *maxresnr;
 457             }
 458         }
 459     }
 460
 461     dest->nres += copies * src->nres;
 462     dest->nr += copies * src->nr;
 463 }
 464
 465 t_atoms gmx_mtop_global_atoms(const gmx_mtop_t& mtop)
 466 {
 467     t_atoms atoms;
 468
 469     init_t_atoms(&atoms, 0, FALSE);
 470
 471     int maxresnr = mtop.maxResNumberNotRenumbered();
 472     for (const gmx_molblock_t& molb : mtop.molblock)
 473     {
 474         atomcat(&atoms,
 475                 &mtop.moltype[molb.type].atoms,
 476                 molb.nmol,
 477                 mtop.maxResiduesPerMoleculeToTriggerRenumber(),
 478                 &maxresnr);
 479     }
 480
 481     return atoms;
 482 }
 483
 484 /*
 485  * The cat routines below are old code from src/kernel/topcat.c
 486  */
 487
 488 static void ilistcat(int ftype, InteractionList* dest, const InteractionList& src, int copies, int dnum, int snum)
 489 {
 490     const int nral = NRAL(ftype);
 491
 492     size_t destIndex = dest->iatoms.size();
 493     dest->iatoms.resize(dest->iatoms.size() + copies * src.size());
 494
 495     for (int c = 0; c < copies; c++)
 496     {
 497         for (int i = 0; i < src.size();)
 498         {
 499             dest->iatoms[destIndex++] = src.iatoms[i++];
 500             for (int a = 0; a < nral; a++)
 501             {
 502                 dest->iatoms[destIndex++] = dnum + src.iatoms[i++];
 503             }
 504         }
 505         dnum += snum;
 506     }
 507 }
 508
 509 static void ilistcat(int ftype, t_ilist* dest, const InteractionList& src, int copies, int dnum, int snum)
 510 {
 511     const int nral = NRAL(ftype);
 512
 513     dest->nalloc = dest->nr + copies * src.size();
 514     srenew(dest->iatoms, dest->nalloc);
 515
 516     for (int c = 0; c < copies; c++)
 517     {
 518         for (int i = 0; i < src.size();)
 519         {
 520             dest->iatoms[dest->nr++] = src.iatoms[i++];
 521             for (int a = 0; a < nral; a++)
 522             {
 523                 dest->iatoms[dest->nr++] = dnum + src.iatoms[i++];
 524             }
 525         }
 526         dnum += snum;
 527     }
 528 }
 529
 530 static const t_iparams& getIparams(const InteractionDefinitions& idef, const int index)
 531 {
 532     return idef.iparams[index];
 533 }
 534
 535 static const t_iparams& getIparams(const t_idef& idef, const int index)
 536 {
 537     return idef.iparams[index];
 538 }
 539
 540 static void resizeIParams(std::vector<t_iparams>* iparams, const int newSize)
 541 {
 542     iparams->resize(newSize);
 543 }
 544
 545 static void resizeIParams(t_iparams** iparams, const int newSize)
 546 {
 547     srenew(*iparams, newSize);
 548 }
 549
 550 template<typename IdefType>
 551 static void set_posres_params(IdefType* idef, const gmx_molblock_t* molb, int i0, int a_offset)
 552 {
 553     auto* il = &idef->il[F_POSRES];
 554     int   i1 = il->size() / 2;
 555     resizeIParams(&idef->iparams_posres, i1);
 556     for (int i = i0; i < i1; i++)
 557     {
 558         t_iparams* ip = &idef->iparams_posres[i];
 559         /* Copy the force constants */
 560         *ip        = getIparams(*idef, il->iatoms[i * 2]);
 561         int a_molb = il->iatoms[i * 2 + 1] - a_offset;
 562         if (molb->posres_xA.empty())
 563         {
 564             gmx_incons("Position restraint coordinates are missing");
 565         }
 566         ip->posres.pos0A[XX] = molb->posres_xA[a_molb][XX];
 567         ip->posres.pos0A[YY] = molb->posres_xA[a_molb][YY];
 568         ip->posres.pos0A[ZZ] = molb->posres_xA[a_molb][ZZ];
 569         if (!molb->posres_xB.empty())
 570         {
 571             ip->posres.pos0B[XX] = molb->posres_xB[a_molb][XX];
 572             ip->posres.pos0B[YY] = molb->posres_xB[a_molb][YY];
 573             ip->posres.pos0B[ZZ] = molb->posres_xB[a_molb][ZZ];
 574         }
 575         else
 576         {
 577             ip->posres.pos0B[XX] = ip->posres.pos0A[XX];
 578             ip->posres.pos0B[YY] = ip->posres.pos0A[YY];
 579             ip->posres.pos0B[ZZ] = ip->posres.pos0A[ZZ];
 580         }
 581         /* Set the parameter index for idef->iparams_posre */
 582         il->iatoms[i * 2] = i;
 583     }
 584 }
 585
 586 template<typename IdefType>
 587 static void set_fbposres_params(IdefType* idef, const gmx_molblock_t* molb, int i0, int a_offset)
 588 {
 589     auto* il = &idef->il[F_FBPOSRES];
 590     int   i1 = il->size() / 2;
 591     resizeIParams(&idef->iparams_fbposres, i1);
 592     for (int i = i0; i < i1; i++)
 593     {
 594         t_iparams* ip = &idef->iparams_fbposres[i];
 595         /* Copy the force constants */
 596         *ip        = getIparams(*idef, il->iatoms[i * 2]);
 597         int a_molb = il->iatoms[i * 2 + 1] - a_offset;
 598         if (molb->posres_xA.empty())
 599         {
 600             gmx_incons("Position restraint coordinates are missing");
 601         }
 602         /* Take flat-bottom posres reference from normal position restraints */
 603         ip->fbposres.pos0[XX] = molb->posres_xA[a_molb][XX];
 604         ip->fbposres.pos0[YY] = molb->posres_xA[a_molb][YY];
 605         ip->fbposres.pos0[ZZ] = molb->posres_xA[a_molb][ZZ];
 606         /* Note: no B-type for flat-bottom posres */
 607
 608         /* Set the parameter index for idef->iparams_posre */
 609         il->iatoms[i * 2] = i;
 610     }
 611 }
 612
 613 /*! \brief Copy parameters to idef structure from mtop.
 614  *
 615  * Makes a deep copy of the force field parameters data structure from a gmx_mtop_t.
 616  * Used to initialize legacy topology types.
 617  *
 618  * \param[in] mtop Reference to input mtop.
 619  * \param[in] idef Pointer to idef to populate.
 620  */
 621 static void copyFFParametersFromMtop(const gmx_mtop_t& mtop, t_idef* idef)
 622 {
 623     const gmx_ffparams_t* ffp = &mtop.ffparams;
 624
 625     idef->ntypes = ffp->numTypes();
 626     idef->atnr   = ffp->atnr;
 627     /* we can no longer copy the pointers to the mtop members,
 628      * because they will become invalid as soon as mtop gets free'd.
 629      * We also need to make sure to only operate on valid data!
 630      */
 631
 632     if (!ffp->functype.empty())
 633     {
 634         snew(idef->functype, ffp->functype.size());
 635         std::copy(ffp->functype.data(), ffp->functype.data() + ffp->functype.size(), idef->functype);
 636     }
 637     else
 638     {
 639         idef->functype = nullptr;
 640     }
 641     if (!ffp->iparams.empty())
 642     {
 643         snew(idef->iparams, ffp->iparams.size());
 644         std::copy(ffp->iparams.data(), ffp->iparams.data() + ffp->iparams.size(), idef->iparams);
 645     }
 646     else
 647     {
 648         idef->iparams = nullptr;
 649     }
 650     idef->iparams_posres   = nullptr;
 651     idef->iparams_fbposres = nullptr;
 652     idef->fudgeQQ          = ffp->fudgeQQ;
 653     idef->ilsort           = ilsortUNKNOWN;
 654 }
 655
 656 /*! \brief Copy idef structure from mtop.
 657  *
 658  * Makes a deep copy of an idef data structure from a gmx_mtop_t.
 659  * Used to initialize legacy topology types.
 660  *
 661  * \param[in] mtop Reference to input mtop.
 662  * \param[in] idef Pointer to idef to populate.
 663  * \param[in] mergeConstr Decide if constraints will be merged.
 664  */
 665 template<typename IdefType>
 666 static void copyIListsFromMtop(const gmx_mtop_t& mtop, IdefType* idef, bool mergeConstr)
 667 {
 668     int natoms = 0;
 669     for (const gmx_molblock_t& molb : mtop.molblock)
 670     {
 671         const gmx_moltype_t& molt = mtop.moltype[molb.type];
 672
 673         int srcnr  = molt.atoms.nr;
 674         int destnr = natoms;
 675
 676         int nposre_old   = idef->il[F_POSRES].size();
 677         int nfbposre_old = idef->il[F_FBPOSRES].size();
 678         for (int ftype = 0; ftype < F_NRE; ftype++)
 679         {
 680             if (mergeConstr && ftype == F_CONSTR && !molt.ilist[F_CONSTRNC].empty())
 681             {
 682                 /* Merge all constrains into one ilist.
 683                  * This simplifies the constraint code.
 684                  */
 685                 for (int mol = 0; mol < molb.nmol; mol++)
 686                 {
 687                     ilistcat(ftype, &idef->il[F_CONSTR], molt.ilist[F_CONSTR], 1, destnr + mol * srcnr, srcnr);
 688                     ilistcat(ftype, &idef->il[F_CONSTR], molt.ilist[F_CONSTRNC], 1, destnr + mol * srcnr, srcnr);
 689                 }
 690             }
 691             else if (!(mergeConstr && ftype == F_CONSTRNC))
 692             {
 693                 ilistcat(ftype, &idef->il[ftype], molt.ilist[ftype], molb.nmol, destnr, srcnr);
 694             }
 695         }
 696         if (idef->il[F_POSRES].size() > nposre_old)
 697         {
 698             /* Executing this line line stops gmxdump -sys working
 699              * correctly. I'm not aware there's an elegant fix. */
 700             set_posres_params(idef, &molb, nposre_old / 2, natoms);
 701         }
 702         if (idef->il[F_FBPOSRES].size() > nfbposre_old)
 703         {
 704             set_fbposres_params(idef, &molb, nfbposre_old / 2, natoms);
 705         }
 706
 707         natoms += molb.nmol * srcnr;
 708     }
 709
 710     if (mtop.bIntermolecularInteractions)
 711     {
 712         for (int ftype = 0; ftype < F_NRE; ftype++)
 713         {
 714             ilistcat(ftype, &idef->il[ftype], (*mtop.intermolecular_ilist)[ftype], 1, 0, mtop.natoms);
 715         }
 716     }
 717
 718     // We have not (yet) sorted free-energy interactions to the end of the ilists
 719     idef->ilsort = ilsortNO_FE;
 720 }
 721
 722 /*! \brief Copy atomtypes from mtop
 723  *
 724  * Makes a deep copy of t_atomtypes from gmx_mtop_t.
 725  * Used to initialize legacy topology types.
 726  *
 727  * \param[in] mtop Reference to input mtop.
 728  * \param[in] atomtypes Pointer to atomtypes to populate.
 729  */
 730 static void copyAtomtypesFromMtop(const gmx_mtop_t& mtop, t_atomtypes* atomtypes)
 731 {
 732     atomtypes->nr = mtop.atomtypes.nr;
 733     if (mtop.atomtypes.atomnumber)
 734     {
 735         snew(atomtypes->atomnumber, mtop.atomtypes.nr);
 736         std::copy(mtop.atomtypes.atomnumber,
 737                   mtop.atomtypes.atomnumber + mtop.atomtypes.nr,
 738                   atomtypes->atomnumber);
 739     }
 740     else
 741     {
 742         atomtypes->atomnumber = nullptr;
 743     }
 744 }
 745
 746 /*! \brief Generate a single list of lists of exclusions for the whole system
 747  *
 748  * \param[in] mtop  Reference to input mtop.
 749  */
 750 static gmx::ListOfLists<int> globalExclusionLists(const gmx_mtop_t& mtop)
 751 {
 752     gmx::ListOfLists<int> excls;
 753
 754     int atomIndex = 0;
 755     for (const gmx_molblock_t& molb : mtop.molblock)
 756     {
 757         const gmx_moltype_t& molt = mtop.moltype[molb.type];
 758
 759         for (int mol = 0; mol < molb.nmol; mol++)
 760         {
 761             excls.appendListOfLists(molt.excls, atomIndex);
 762
 763             atomIndex += molt.atoms.nr;
 764         }
 765     }
 766
 767     return excls;
 768 }
 769
 770 /*! \brief Updates inter-molecular exclusion lists
 771  *
 772  * This function updates inter-molecular exclusions to exclude all
 773  * non-bonded interactions between a given list of atoms
 774  *
 775  * \param[inout]    excls   existing exclusions in local topology
 776  * \param[in]       ids     list of global IDs of atoms
 777  */
 778 static void addMimicExclusions(gmx::ListOfLists<int>* excls, const gmx::ArrayRef<const int> ids)
 779 {
 780     t_blocka inter_excl{};
 781     init_blocka(&inter_excl);
 782     size_t n_q = ids.size();
 783
 784     inter_excl.nr  = excls->ssize();
 785     inter_excl.nra = n_q * n_q;
 786
 787     size_t total_nra = n_q * n_q;
 788
 789     snew(inter_excl.index, excls->ssize() + 1);
 790     snew(inter_excl.a, total_nra);
 791
 792     for (int i = 0; i < inter_excl.nr; ++i)
 793     {
 794         inter_excl.index[i] = 0;
 795     }
 796
 797     /* Here we loop over the list of QM atom ids
 798      *  and create exclusions between all of them resulting in
 799      *  n_q * n_q sized exclusion list
 800      */
 801     int prev_index = 0;
 802     for (int k = 0; k < inter_excl.nr; ++k)
 803     {
 804         inter_excl.index[k] = prev_index;
 805         for (long i = 0; i < ids.ssize(); i++)
 806         {
 807             if (k != ids[i])
 808             {
 809                 continue;
 810             }
 811             size_t index             = n_q * i;
 812             inter_excl.index[ids[i]] = index;
 813             prev_index               = index + n_q;
 814             for (size_t j = 0; j < n_q; ++j)
 815             {
 816                 inter_excl.a[n_q * i + j] = ids[j];
 817             }
 818         }
 819     }
 820     inter_excl.index[ids[n_q - 1] + 1] = n_q * n_q;
 821
 822     inter_excl.index[inter_excl.nr] = n_q * n_q;
 823
 824     std::vector<gmx::ExclusionBlock> qmexcl2(excls->size());
 825     gmx::blockaToExclusionBlocks(&inter_excl, qmexcl2);
 826
 827     // Merge the created exclusion list with the existing one
 828     gmx::mergeExclusions(excls, qmexcl2);
 829 }
 830
 831 static void sortFreeEnergyInteractionsAtEnd(const gmx_mtop_t& mtop, InteractionDefinitions* idef)
 832 {
 833     std::vector<real> qA(mtop.natoms);
 834     std::vector<real> qB(mtop.natoms);
 835     for (const AtomProxy atomP : AtomRange(mtop))
 836     {
 837         const t_atom& local = atomP.atom();
 838         int           index = atomP.globalAtomNumber();
 839         qA[index]           = local.q;
 840         qB[index]           = local.qB;
 841     }
 842     gmx_sort_ilist_fe(idef, qA.data(), qB.data());
 843 }
 844
 845 static void gen_local_top(const gmx_mtop_t& mtop,
 846                           bool              freeEnergyInteractionsAtEnd,
 847                           bool              bMergeConstr,
 848                           gmx_localtop_t*   top)
 849 {
 850     copyIListsFromMtop(mtop, &top->idef, bMergeConstr);
 851     if (freeEnergyInteractionsAtEnd)
 852     {
 853         sortFreeEnergyInteractionsAtEnd(mtop, &top->idef);
 854     }
 855     top->excls = globalExclusionLists(mtop);
 856     if (!mtop.intermolecularExclusionGroup.empty())
 857     {
 858         addMimicExclusions(&top->excls, mtop.intermolecularExclusionGroup);
 859     }
 860 }
 861
 862 void gmx_mtop_generate_local_top(const gmx_mtop_t& mtop, gmx_localtop_t* top, bool freeEnergyInteractionsAtEnd)
 863 {
 864     gen_local_top(mtop, freeEnergyInteractionsAtEnd, true, top);
 865 }
 866
 867 /*! \brief Fills an array with molecule begin/end atom indices
 868  *
 869  * \param[in]  mtop   The global topology
 870  * \param[out] index  Array of size nr. of molecules + 1 to be filled with molecule begin/end indices
 871  */
 872 static void fillMoleculeIndices(const gmx_mtop_t& mtop, gmx::ArrayRef<int> index)
 873 {
 874     int globalAtomIndex   = 0;
 875     int globalMolIndex    = 0;
 876     index[globalMolIndex] = globalAtomIndex;
 877     for (const gmx_molblock_t& molb : mtop.molblock)
 878     {
 879         int numAtomsPerMolecule = mtop.moltype[molb.type].atoms.nr;
 880         for (int mol = 0; mol < molb.nmol; mol++)
 881         {
 882             globalAtomIndex += numAtomsPerMolecule;
 883             globalMolIndex += 1;
 884             index[globalMolIndex] = globalAtomIndex;
 885         }
 886     }
 887 }
 888
 889 gmx::RangePartitioning gmx_mtop_molecules(const gmx_mtop_t& mtop)
 890 {
 891     gmx::RangePartitioning mols;
 892
 893     for (const gmx_molblock_t& molb : mtop.molblock)
 894     {
 895         int numAtomsPerMolecule = mtop.moltype[molb.type].atoms.nr;
 896         for (int mol = 0; mol < molb.nmol; mol++)
 897         {
 898             mols.appendBlock(numAtomsPerMolecule);
 899         }
 900     }
 901
 902     return mols;
 903 }
 904
 905 std::vector<gmx::Range<int>> atomRangeOfEachResidue(const gmx_moltype_t& moltype)
 906 {
 907     std::vector<gmx::Range<int>> atomRanges;
 908     int                          currentResidueNumber = moltype.atoms.atom[0].resind;
 909     int                          startAtom            = 0;
 910     // Go through all atoms in a molecule to store first and last atoms in each residue.
 911     for (int i = 0; i < moltype.atoms.nr; i++)
 912     {
 913         int residueOfThisAtom = moltype.atoms.atom[i].resind;
 914         if (residueOfThisAtom != currentResidueNumber)
 915         {
 916             // This atom belongs to the next residue, so record the range for the previous residue,
 917             // remembering that end points to one place past the last atom.
 918             int endAtom = i;
 919             atomRanges.emplace_back(startAtom, endAtom);
 920             // Prepare for the current residue
 921             startAtom            = endAtom;
 922             currentResidueNumber = residueOfThisAtom;
 923         }
 924     }
 925     // special treatment for last residue in this molecule.
 926     atomRanges.emplace_back(startAtom, moltype.atoms.nr);
 927
 928     return atomRanges;
 929 }
 930
 931 /*! \brief Creates and returns a deprecated t_block struct with molecule indices
 932  *
 933  * \param[in] mtop  The global topology
 934  */
 935 static t_block gmx_mtop_molecules_t_block(const gmx_mtop_t& mtop)
 936 {
 937     t_block mols;
 938
 939     mols.nr           = gmx_mtop_num_molecules(mtop);
 940     mols.nalloc_index = mols.nr + 1;
 941     snew(mols.index, mols.nalloc_index);
 942
 943     fillMoleculeIndices(mtop, gmx::arrayRefFromArray(mols.index, mols.nr + 1));
 944
 945     return mols;
 946 }
 947
 948 static void gen_t_topology(const gmx_mtop_t& mtop, bool bMergeConstr, t_topology* top)
 949 {
 950     copyAtomtypesFromMtop(mtop, &top->atomtypes);
 951     for (int ftype = 0; ftype < F_NRE; ftype++)
 952     {
 953         top->idef.il[ftype].nr     = 0;
 954         top->idef.il[ftype].nalloc = 0;
 955         top->idef.il[ftype].iatoms = nullptr;
 956     }
 957     copyFFParametersFromMtop(mtop, &top->idef);
 958     copyIListsFromMtop(mtop, &top->idef, bMergeConstr);
 959
 960     top->name                        = mtop.name;
 961     top->atoms                       = gmx_mtop_global_atoms(mtop);
 962     top->mols                        = gmx_mtop_molecules_t_block(mtop);
 963     top->bIntermolecularInteractions = mtop.bIntermolecularInteractions;
 964     top->symtab                      = mtop.symtab;
 965 }
 966
 967 t_topology gmx_mtop_t_to_t_topology(gmx_mtop_t* mtop, bool freeMTop)
 968 {
 969     t_topology top;
 970
 971     gen_t_topology(*mtop, false, &top);
 972
 973     if (freeMTop)
 974     {
 975         // Clear pointers and counts, such that the pointers copied to top
 976         // keep pointing to valid data after destroying mtop.
 977         mtop->symtab.symbuf = nullptr;
 978         mtop->symtab.nr     = 0;
 979     }
 980     return top;
 981 }
 982
 983 std::vector<int> get_atom_index(const gmx_mtop_t& mtop)
 984 {
 985
 986     std::vector<int> atom_index;
 987     for (const AtomProxy atomP : AtomRange(mtop))
 988     {
 989         const t_atom& local = atomP.atom();
 990         int           index = atomP.globalAtomNumber();
 991         if (local.ptype == ParticleType::Atom)
 992         {
 993             atom_index.push_back(index);
 994         }
 995     }
 996     return atom_index;
 997 }
 998
 999 void convertAtomsToMtop(t_symtab* symtab, char** name, t_atoms* atoms, gmx_mtop_t* mtop)
1000 {
1001     mtop->symtab = *symtab;
1002
1003     mtop->name = name;
1004
1005     mtop->moltype.clear();
1006     mtop->moltype.resize(1);
1007     mtop->moltype.back().atoms = *atoms;
1008
1009     mtop->molblock.resize(1);
1010     mtop->molblock[0].type = 0;
1011     mtop->molblock[0].nmol = 1;
1012
1013     mtop->bIntermolecularInteractions = FALSE;
1014
1015     mtop->natoms = atoms->nr;
1016
1017     mtop->haveMoleculeIndices = false;
1018
1019     mtop->finalize();
1020 }
1021
1022 bool haveFepPerturbedNBInteractions(const gmx_mtop_t& mtop)
1023 {
1024     for (const gmx_moltype_t& molt : mtop.moltype)
1025     {
1026         for (int a = 0; a < molt.atoms.nr; a++)
1027         {
1028             if (PERTURBED(molt.atoms.atom[a]))
1029             {
1030                 return true;
1031             }
1032         }
1033     }
1034
1035     return false;
1036 }
1037
1038 bool haveFepPerturbedMasses(const gmx_mtop_t& mtop)
1039 {
1040     for (const gmx_moltype_t& molt : mtop.moltype)
1041     {
1042         for (int a = 0; a < molt.atoms.nr; a++)
1043         {
1044             const t_atom& atom = molt.atoms.atom[a];
1045             if (atom.m != atom.mB)
1046             {
1047                 return true;
1048             }
1049         }
1050     }
1051
1052     return false;
1053 }
1054
1055 bool haveFepPerturbedMassesInSettles(const gmx_mtop_t& mtop)
1056 {
1057     for (const gmx_moltype_t& molt : mtop.moltype)
1058     {
1059         if (!molt.ilist[F_SETTLE].empty())
1060         {
1061             for (int a = 0; a < molt.atoms.nr; a++)
1062             {
1063                 const t_atom& atom = molt.atoms.atom[a];
1064                 if (atom.m != atom.mB)
1065                 {
1066                     return true;
1067                 }
1068             }
1069         }
1070     }
1071     return false;
1072 }
1073
1074 bool havePerturbedConstraints(const gmx_mtop_t& mtop)
1075 {
1076     // This code assumes that all perturbed constraints parameters are actually used
1077     const auto& ffparams = mtop.ffparams;
1078
1079     for (gmx::index i = 0; i < gmx::ssize(ffparams.functype); i++)
1080     {
1081         if (ffparams.functype[i] == F_CONSTR || ffparams.functype[i] == F_CONSTRNC)
1082         {
1083             const auto& iparams = ffparams.iparams[i];
1084             if (iparams.constr.dA != iparams.constr.dB)
1085             {
1086                 return true;
1087             }
1088         }
1089     }
1090
1091     return false;
1092 }