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