src/gromacs/topology/mtop_util.cpp

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