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