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