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