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