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 t_block gmx_mtop_global_cgs(const gmx_mtop_t *mtop)
 534 {
 535     t_block cgs_gl;
 536     /* In most cases this is too much, but we realloc at the end */
 537     snew(cgs_gl.index, mtop->natoms+1);
 538
 539     cgs_gl.nr       = 0;
 540     cgs_gl.index[0] = 0;
 541     for (const gmx_molblock_t &molb : mtop->molblock)
 542     {
 543         const t_block &cgs_mol = mtop->moltype[molb.type].cgs;
 544         for (int mol = 0; mol < molb.nmol; mol++)
 545         {
 546             for (int cg = 0; cg < cgs_mol.nr; cg++)
 547             {
 548                 cgs_gl.index[cgs_gl.nr + 1] =
 549                     cgs_gl.index[cgs_gl.nr] +
 550                     cgs_mol.index[cg + 1] - cgs_mol.index[cg];
 551                 cgs_gl.nr++;
 552             }
 553         }
 554     }
 555     cgs_gl.nalloc_index = cgs_gl.nr + 1;
 556     srenew(cgs_gl.index, cgs_gl.nalloc_index);
 557
 558     return cgs_gl;
 559 }
 560
 561 static void atomcat(t_atoms *dest, const t_atoms *src, int copies,
 562                     int maxres_renum, int *maxresnr)
 563 {
 564     int i, j, l, size;
 565     int srcnr  = src->nr;
 566     int destnr = dest->nr;
 567
 568     if (dest->nr == 0)
 569     {
 570         dest->haveMass    = src->haveMass;
 571         dest->haveType    = src->haveType;
 572         dest->haveCharge  = src->haveCharge;
 573         dest->haveBState  = src->haveBState;
 574         dest->havePdbInfo = src->havePdbInfo;
 575     }
 576     else
 577     {
 578         dest->haveMass    = dest->haveMass    && src->haveMass;
 579         dest->haveType    = dest->haveType    && src->haveType;
 580         dest->haveCharge  = dest->haveCharge  && src->haveCharge;
 581         dest->haveBState  = dest->haveBState  && src->haveBState;
 582         dest->havePdbInfo = dest->havePdbInfo && src->havePdbInfo;
 583     }
 584
 585     if (srcnr)
 586     {
 587         size = destnr+copies*srcnr;
 588         srenew(dest->atom, size);
 589         srenew(dest->atomname, size);
 590         if (dest->haveType)
 591         {
 592             srenew(dest->atomtype, size);
 593             if (dest->haveBState)
 594             {
 595                 srenew(dest->atomtypeB, size);
 596             }
 597         }
 598         if (dest->havePdbInfo)
 599         {
 600             srenew(dest->pdbinfo, size);
 601         }
 602     }
 603     if (src->nres)
 604     {
 605         size = dest->nres+copies*src->nres;
 606         srenew(dest->resinfo, size);
 607     }
 608
 609     /* residue information */
 610     for (l = dest->nres, j = 0; (j < copies); j++, l += src->nres)
 611     {
 612         memcpy(reinterpret_cast<char *>(&(dest->resinfo[l])), reinterpret_cast<char *>(&(src->resinfo[0])),
 613                static_cast<size_t>(src->nres*sizeof(src->resinfo[0])));
 614     }
 615
 616     for (l = destnr, j = 0; (j < copies); j++, l += srcnr)
 617     {
 618         memcpy(reinterpret_cast<char *>(&(dest->atom[l])), reinterpret_cast<char *>(&(src->atom[0])),
 619                static_cast<size_t>(srcnr*sizeof(src->atom[0])));
 620         memcpy(reinterpret_cast<char *>(&(dest->atomname[l])), reinterpret_cast<char *>(&(src->atomname[0])),
 621                static_cast<size_t>(srcnr*sizeof(src->atomname[0])));
 622         if (dest->haveType)
 623         {
 624             memcpy(reinterpret_cast<char *>(&(dest->atomtype[l])), reinterpret_cast<char *>(&(src->atomtype[0])),
 625                    static_cast<size_t>(srcnr*sizeof(src->atomtype[0])));
 626             if (dest->haveBState)
 627             {
 628                 memcpy(reinterpret_cast<char *>(&(dest->atomtypeB[l])), reinterpret_cast<char *>(&(src->atomtypeB[0])),
 629                        static_cast<size_t>(srcnr*sizeof(src->atomtypeB[0])));
 630             }
 631         }
 632         if (dest->havePdbInfo)
 633         {
 634             memcpy(reinterpret_cast<char *>(&(dest->pdbinfo[l])), reinterpret_cast<char *>(&(src->pdbinfo[0])),
 635                    static_cast<size_t>(srcnr*sizeof(src->pdbinfo[0])));
 636         }
 637     }
 638
 639     /* Increment residue indices */
 640     for (l = destnr, j = 0; (j < copies); j++)
 641     {
 642         for (i = 0; (i < srcnr); i++, l++)
 643         {
 644             dest->atom[l].resind = dest->nres+j*src->nres+src->atom[i].resind;
 645         }
 646     }
 647
 648     if (src->nres <= maxres_renum)
 649     {
 650         /* Single residue molecule, continue counting residues */
 651         for (j = 0; (j < copies); j++)
 652         {
 653             for (l = 0; l < src->nres; l++)
 654             {
 655                 (*maxresnr)++;
 656                 dest->resinfo[dest->nres+j*src->nres+l].nr = *maxresnr;
 657             }
 658         }
 659     }
 660
 661     dest->nres += copies*src->nres;
 662     dest->nr   += copies*src->nr;
 663 }
 664
 665 t_atoms gmx_mtop_global_atoms(const gmx_mtop_t *mtop)
 666 {
 667     t_atoms         atoms;
 668
 669     init_t_atoms(&atoms, 0, FALSE);
 670
 671     int maxresnr = mtop->maxresnr;
 672     for (const gmx_molblock_t &molb : mtop->molblock)
 673     {
 674         atomcat(&atoms, &mtop->moltype[molb.type].atoms, molb.nmol,
 675                 mtop->maxres_renum, &maxresnr);
 676     }
 677
 678     return atoms;
 679 }
 680
 681 /*
 682  * The cat routines below are old code from src/kernel/topcat.c
 683  */
 684
 685 static void blockcat(t_block *dest, const t_block *src, int copies)
 686 {
 687     int i, j, l, nra, size;
 688
 689     if (src->nr)
 690     {
 691         size = (dest->nr+copies*src->nr+1);
 692         srenew(dest->index, size);
 693     }
 694
 695     nra = dest->index[dest->nr];
 696     for (l = dest->nr, j = 0; (j < copies); j++)
 697     {
 698         for (i = 0; (i < src->nr); i++)
 699         {
 700             dest->index[l++] = nra + src->index[i];
 701         }
 702         if (src->nr > 0)
 703         {
 704             nra += src->index[src->nr];
 705         }
 706     }
 707     dest->nr             += copies*src->nr;
 708     dest->index[dest->nr] = nra;
 709 }
 710
 711 static void blockacat(t_blocka *dest, const t_blocka *src, int copies,
 712                       int dnum, int snum)
 713 {
 714     int i, j, l, size;
 715     int destnr  = dest->nr;
 716     int destnra = dest->nra;
 717
 718     if (src->nr)
 719     {
 720         size = (dest->nr+copies*src->nr+1);
 721         srenew(dest->index, size);
 722     }
 723     if (src->nra)
 724     {
 725         size = (dest->nra+copies*src->nra);
 726         srenew(dest->a, size);
 727     }
 728
 729     for (l = destnr, j = 0; (j < copies); j++)
 730     {
 731         for (i = 0; (i < src->nr); i++)
 732         {
 733             dest->index[l++] = dest->nra+src->index[i];
 734         }
 735         dest->nra += src->nra;
 736     }
 737     for (l = destnra, j = 0; (j < copies); j++)
 738     {
 739         for (i = 0; (i < src->nra); i++)
 740         {
 741             dest->a[l++] = dnum+src->a[i];
 742         }
 743         dnum     += snum;
 744         dest->nr += src->nr;
 745     }
 746     dest->index[dest->nr] = dest->nra;
 747     dest->nalloc_index    = dest->nr;
 748     dest->nalloc_a        = dest->nra;
 749 }
 750
 751 static void ilistcat(int                    ftype,
 752                      t_ilist               *dest,
 753                      const InteractionList &src,
 754                      int                    copies,
 755                      int                    dnum,
 756                      int                    snum)
 757 {
 758     int nral, c, i, a;
 759
 760     nral = NRAL(ftype);
 761
 762     dest->nalloc = dest->nr + copies*src.size();
 763     srenew(dest->iatoms, dest->nalloc);
 764
 765     for (c = 0; c < copies; c++)
 766     {
 767         for (i = 0; i < src.size(); )
 768         {
 769             dest->iatoms[dest->nr++] = src.iatoms[i++];
 770             for (a = 0; a < nral; a++)
 771             {
 772                 dest->iatoms[dest->nr++] = dnum + src.iatoms[i++];
 773             }
 774         }
 775         dnum += snum;
 776     }
 777 }
 778
 779 static void set_posres_params(t_idef *idef, const gmx_molblock_t *molb,
 780                               int i0, int a_offset)
 781 {
 782     t_ilist   *il;
 783     int        i1, i, a_molb;
 784     t_iparams *ip;
 785
 786     il = &idef->il[F_POSRES];
 787     i1 = il->nr/2;
 788     idef->iparams_posres_nalloc = i1;
 789     srenew(idef->iparams_posres, idef->iparams_posres_nalloc);
 790     for (i = i0; i < i1; i++)
 791     {
 792         ip = &idef->iparams_posres[i];
 793         /* Copy the force constants */
 794         *ip    = idef->iparams[il->iatoms[i*2]];
 795         a_molb = il->iatoms[i*2+1] - a_offset;
 796         if (molb->posres_xA.empty())
 797         {
 798             gmx_incons("Position restraint coordinates are missing");
 799         }
 800         ip->posres.pos0A[XX] = molb->posres_xA[a_molb][XX];
 801         ip->posres.pos0A[YY] = molb->posres_xA[a_molb][YY];
 802         ip->posres.pos0A[ZZ] = molb->posres_xA[a_molb][ZZ];
 803         if (!molb->posres_xB.empty())
 804         {
 805             ip->posres.pos0B[XX] = molb->posres_xB[a_molb][XX];
 806             ip->posres.pos0B[YY] = molb->posres_xB[a_molb][YY];
 807             ip->posres.pos0B[ZZ] = molb->posres_xB[a_molb][ZZ];
 808         }
 809         else
 810         {
 811             ip->posres.pos0B[XX] = ip->posres.pos0A[XX];
 812             ip->posres.pos0B[YY] = ip->posres.pos0A[YY];
 813             ip->posres.pos0B[ZZ] = ip->posres.pos0A[ZZ];
 814         }
 815         /* Set the parameter index for idef->iparams_posre */
 816         il->iatoms[i*2] = i;
 817     }
 818 }
 819
 820 static void set_fbposres_params(t_idef *idef, const gmx_molblock_t *molb,
 821                                 int i0, int a_offset)
 822 {
 823     t_ilist   *il;
 824     int        i1, i, a_molb;
 825     t_iparams *ip;
 826
 827     il = &idef->il[F_FBPOSRES];
 828     i1 = il->nr/2;
 829     idef->iparams_fbposres_nalloc = i1;
 830     srenew(idef->iparams_fbposres, idef->iparams_fbposres_nalloc);
 831     for (i = i0; i < i1; i++)
 832     {
 833         ip = &idef->iparams_fbposres[i];
 834         /* Copy the force constants */
 835         *ip    = idef->iparams[il->iatoms[i*2]];
 836         a_molb = il->iatoms[i*2+1] - a_offset;
 837         if (molb->posres_xA.empty())
 838         {
 839             gmx_incons("Position restraint coordinates are missing");
 840         }
 841         /* Take flat-bottom posres reference from normal position restraints */
 842         ip->fbposres.pos0[XX] = molb->posres_xA[a_molb][XX];
 843         ip->fbposres.pos0[YY] = molb->posres_xA[a_molb][YY];
 844         ip->fbposres.pos0[ZZ] = molb->posres_xA[a_molb][ZZ];
 845         /* Note: no B-type for flat-bottom posres */
 846
 847         /* Set the parameter index for idef->iparams_posre */
 848         il->iatoms[i*2] = i;
 849     }
 850 }
 851
 852 /*! \brief Copy idef structure from mtop.
 853  *
 854  * Makes a deep copy of an idef data structure from a gmx_mtop_t.
 855  * Used to initialize legacy topology types.
 856  *
 857  * \param[in] mtop Reference to input mtop.
 858  * \param[in] idef Pointer to idef to populate.
 859  * \param[in] mergeConstr Decide if constraints will be merged.
 860  * \param[in] freeEnergyInteractionsAtEnd Decide if free energy stuff should
 861  *              be added at the end.
 862  */
 863 static void copyIdefFromMtop(const gmx_mtop_t &mtop,
 864                              t_idef           *idef,
 865                              bool              freeEnergyInteractionsAtEnd,
 866                              bool              mergeConstr)
 867 {
 868     const gmx_ffparams_t   *ffp = &mtop.ffparams;
 869
 870     idef->ntypes                  = ffp->numTypes();
 871     idef->atnr                    = ffp->atnr;
 872     /* we can no longer copy the pointers to the mtop members,
 873      * because they will become invalid as soon as mtop gets free'd.
 874      * We also need to make sure to only operate on valid data!
 875      */
 876
 877     if (!ffp->functype.empty())
 878     {
 879         snew(idef->functype, ffp->functype.size());
 880         std::copy(ffp->functype.data(), ffp->functype.data() + ffp->functype.size(), idef->functype);
 881     }
 882     else
 883     {
 884         idef->functype = nullptr;
 885     }
 886     if (!ffp->iparams.empty())
 887     {
 888         snew(idef->iparams, ffp->iparams.size());
 889         std::copy(ffp->iparams.data(), ffp->iparams.data() + ffp->iparams.size(), idef->iparams);
 890     }
 891     else
 892     {
 893         idef->iparams = nullptr;
 894     }
 895     idef->iparams_posres          = nullptr;
 896     idef->iparams_posres_nalloc   = 0;
 897     idef->iparams_fbposres        = nullptr;
 898     idef->iparams_fbposres_nalloc = 0;
 899     idef->fudgeQQ                 = ffp->fudgeQQ;
 900     idef->cmap_grid               = new gmx_cmap_t;
 901     *idef->cmap_grid              = ffp->cmap_grid;
 902     idef->ilsort                  = ilsortUNKNOWN;
 903
 904     for (int ftype = 0; ftype < F_NRE; ftype++)
 905     {
 906         idef->il[ftype].nr     = 0;
 907         idef->il[ftype].nalloc = 0;
 908         idef->il[ftype].iatoms = nullptr;
 909     }
 910
 911     int natoms = 0;
 912     for (const gmx_molblock_t &molb : mtop.molblock)
 913     {
 914         const gmx_moltype_t &molt = mtop.moltype[molb.type];
 915
 916         int                  srcnr  = molt.atoms.nr;
 917         int                  destnr = natoms;
 918
 919         int                  nposre_old   = idef->il[F_POSRES].nr;
 920         int                  nfbposre_old = idef->il[F_FBPOSRES].nr;
 921         for (int ftype = 0; ftype < F_NRE; ftype++)
 922         {
 923             if (mergeConstr &&
 924                 ftype == F_CONSTR && molt.ilist[F_CONSTRNC].size() > 0)
 925             {
 926                 /* Merge all constrains into one ilist.
 927                  * This simplifies the constraint code.
 928                  */
 929                 for (int mol = 0; mol < molb.nmol; mol++)
 930                 {
 931                     ilistcat(ftype, &idef->il[F_CONSTR], molt.ilist[F_CONSTR],
 932                              1, destnr + mol*srcnr, srcnr);
 933                     ilistcat(ftype, &idef->il[F_CONSTR], molt.ilist[F_CONSTRNC],
 934                              1, destnr + mol*srcnr, srcnr);
 935                 }
 936             }
 937             else if (!(mergeConstr && ftype == F_CONSTRNC))
 938             {
 939                 ilistcat(ftype, &idef->il[ftype], molt.ilist[ftype],
 940                          molb.nmol, destnr, srcnr);
 941             }
 942         }
 943         if (idef->il[F_POSRES].nr > nposre_old)
 944         {
 945             /* Executing this line line stops gmxdump -sys working
 946              * correctly. I'm not aware there's an elegant fix. */
 947             set_posres_params(idef, &molb, nposre_old/2, natoms);
 948         }
 949         if (idef->il[F_FBPOSRES].nr > nfbposre_old)
 950         {
 951             set_fbposres_params(idef, &molb, nfbposre_old/2, natoms);
 952         }
 953
 954         natoms += molb.nmol*srcnr;
 955     }
 956
 957     if (mtop.bIntermolecularInteractions)
 958     {
 959         for (int ftype = 0; ftype < F_NRE; ftype++)
 960         {
 961             ilistcat(ftype, &idef->il[ftype], (*mtop.intermolecular_ilist)[ftype],
 962                      1, 0, mtop.natoms);
 963         }
 964     }
 965
 966     if (freeEnergyInteractionsAtEnd && gmx_mtop_bondeds_free_energy(&mtop))
 967     {
 968         std::vector<real>          qA(mtop.natoms);
 969         std::vector<real>          qB(mtop.natoms);
 970         for (const AtomProxy atomP : AtomRange(mtop))
 971         {
 972             const t_atom &local = atomP.atom();
 973             int           index = atomP.globalAtomNumber();
 974             qA[index] = local.q;
 975             qB[index] = local.qB;
 976         }
 977         gmx_sort_ilist_fe(idef, qA.data(), qB.data());
 978     }
 979     else
 980     {
 981         idef->ilsort = ilsortNO_FE;
 982     }
 983 }
 984
 985 /*! \brief Copy atomtypes from mtop
 986  *
 987  * Makes a deep copy of t_atomtypes from gmx_mtop_t.
 988  * Used to initialize legacy topology types.
 989  *
 990  * \param[in] mtop Reference to input mtop.
 991  * \param[in] atomtypes Pointer to atomtypes to populate.
 992  */
 993 static void copyAtomtypesFromMtop(const gmx_mtop_t &mtop,
 994                                   t_atomtypes      *atomtypes)
 995 {
 996     atomtypes->nr = mtop.atomtypes.nr;
 997     if (mtop.atomtypes.atomnumber)
 998     {
 999         snew(atomtypes->atomnumber, mtop.atomtypes.nr);
1000         std::copy(mtop.atomtypes.atomnumber, mtop.atomtypes.atomnumber + mtop.atomtypes.nr, atomtypes->atomnumber);
1001     }
1002     else
1003     {
1004         atomtypes->atomnumber = nullptr;
1005     }
1006 }
1007
1008 /*! \brief Copy cgs from mtop.
1009  *
1010  * Makes a deep copy of cgs(t_block) from gmx_mtop_t.
1011  * Used to initialize legacy topology types.
1012  *
1013  * \param[in] mtop Reference to input mtop.
1014  * \param[in] cgs  Pointer to final cgs data structure.
1015  */
1016 static void copyCgsFromMtop(const gmx_mtop_t &mtop,
1017                             t_block          *cgs)
1018 {
1019     init_block(cgs);
1020     for (const gmx_molblock_t &molb : mtop.molblock)
1021     {
1022         const gmx_moltype_t &molt = mtop.moltype[molb.type];
1023         blockcat(cgs, &molt.cgs, molb.nmol);
1024     }
1025 }
1026
1027 /*! \brief Copy excls from mtop.
1028  *
1029  * Makes a deep copy of excls(t_blocka) from gmx_mtop_t.
1030  * Used to initialize legacy topology types.
1031  *
1032  * \param[in] mtop  Reference to input mtop.
1033  * \param[in] excls Pointer to final excls data structure.
1034  */
1035 static void copyExclsFromMtop(const gmx_mtop_t &mtop,
1036                               t_blocka         *excls)
1037 {
1038     init_blocka(excls);
1039     int natoms = 0;
1040     for (const gmx_molblock_t &molb : mtop.molblock)
1041     {
1042         const gmx_moltype_t &molt = mtop.moltype[molb.type];
1043
1044         int                  srcnr  = molt.atoms.nr;
1045         int                  destnr = natoms;
1046
1047         blockacat(excls, &molt.excls, molb.nmol, destnr, srcnr);
1048
1049         natoms += molb.nmol*srcnr;
1050     }
1051 }
1052
1053 /*! \brief Updates inter-molecular exclusion lists
1054  *
1055  * This function updates inter-molecular exclusions to exclude all
1056  * non-bonded interactions between a given list of atoms
1057  *
1058  * \param[inout]    excls   existing exclusions in local topology
1059  * \param[in]       ids     list of global IDs of atoms
1060  */
1061 static void addMimicExclusions(t_blocka                      *excls,
1062                                const gmx::ArrayRef<const int> ids)
1063 {
1064     t_blocka inter_excl {};
1065     init_blocka(&inter_excl);
1066     size_t   n_q = ids.size();
1067
1068     inter_excl.nr  = excls->nr;
1069     inter_excl.nra = n_q * n_q;
1070
1071     size_t total_nra = n_q * n_q;
1072
1073     snew(inter_excl.index, excls->nr + 1);
1074     snew(inter_excl.a, total_nra);
1075
1076     for (int i = 0; i < excls->nr; ++i)
1077     {
1078         inter_excl.index[i] = 0;
1079     }
1080
1081     /* Here we loop over the list of QM atom ids
1082      *  and create exclusions between all of them resulting in
1083      *  n_q * n_q sized exclusion list
1084      */
1085     int prev_index = 0;
1086     for (int k = 0; k < inter_excl.nr; ++k)
1087     {
1088         inter_excl.index[k] = prev_index;
1089         for (long i = 0; i < ids.ssize(); i++)
1090         {
1091             if (k != ids[i])
1092             {
1093                 continue;
1094             }
1095             size_t index = n_q * i;
1096             inter_excl.index[ids[i]]     = index;
1097             prev_index                   = index + n_q;
1098             for (size_t j = 0; j < n_q; ++j)
1099             {
1100                 inter_excl.a[n_q * i + j] = ids[j];
1101             }
1102         }
1103     }
1104     inter_excl.index[ids[n_q - 1] + 1] = n_q * n_q;
1105
1106     inter_excl.index[inter_excl.nr] = n_q * n_q;
1107
1108     std::vector<gmx::ExclusionBlock> qmexcl2(excls->nr);
1109     gmx::blockaToExclusionBlocks(&inter_excl, qmexcl2);
1110
1111     // Merge the created exclusion list with the existing one
1112     gmx::mergeExclusions(excls, qmexcl2);
1113 }
1114
1115 static void gen_local_top(const gmx_mtop_t  &mtop,
1116                           bool               freeEnergyInteractionsAtEnd,
1117                           bool               bMergeConstr,
1118                           gmx_localtop_t    *top)
1119 {
1120     copyAtomtypesFromMtop(mtop, &top->atomtypes);
1121     copyIdefFromMtop(mtop, &top->idef, freeEnergyInteractionsAtEnd, bMergeConstr);
1122     copyExclsFromMtop(mtop, &top->excls);
1123     if (!mtop.intermolecularExclusionGroup.empty())
1124     {
1125         addMimicExclusions(&top->excls,
1126                            mtop.intermolecularExclusionGroup);
1127     }
1128 }
1129
1130 void
1131 gmx_mtop_generate_local_top(const gmx_mtop_t &mtop,
1132                             gmx_localtop_t   *top,
1133                             bool              freeEnergyInteractionsAtEnd)
1134 {
1135     gen_local_top(mtop, freeEnergyInteractionsAtEnd, true, top);
1136 }
1137
1138 /*! \brief Fills an array with molecule begin/end atom indices
1139  *
1140  * \param[in]  mtop   The global topology
1141  * \param[out] index  Array of size nr. of molecules + 1 to be filled with molecule begin/end indices
1142  */
1143 static void fillMoleculeIndices(const gmx_mtop_t  &mtop,
1144                                 gmx::ArrayRef<int> index)
1145 {
1146     int globalAtomIndex   = 0;
1147     int globalMolIndex    = 0;
1148     index[globalMolIndex] = globalAtomIndex;
1149     for (const gmx_molblock_t &molb : mtop.molblock)
1150     {
1151         int numAtomsPerMolecule = mtop.moltype[molb.type].atoms.nr;
1152         for (int mol = 0; mol < molb.nmol; mol++)
1153         {
1154             globalAtomIndex       += numAtomsPerMolecule;
1155             globalMolIndex        += 1;
1156             index[globalMolIndex]  = globalAtomIndex;
1157         }
1158     }
1159 }
1160
1161 gmx::RangePartitioning gmx_mtop_molecules(const gmx_mtop_t &mtop)
1162 {
1163     gmx::RangePartitioning mols;
1164
1165     for (const gmx_molblock_t &molb : mtop.molblock)
1166     {
1167         int numAtomsPerMolecule = mtop.moltype[molb.type].atoms.nr;
1168         for (int mol = 0; mol < molb.nmol; mol++)
1169         {
1170             mols.appendBlock(numAtomsPerMolecule);
1171         }
1172     }
1173
1174     return mols;
1175 }
1176
1177 /*! \brief Creates and returns a deprecated t_block struct with molecule indices
1178  *
1179  * \param[in] mtop  The global topology
1180  */
1181 static t_block gmx_mtop_molecules_t_block(const gmx_mtop_t &mtop)
1182 {
1183     t_block mols;
1184
1185     mols.nr           = gmx_mtop_num_molecules(mtop);
1186     mols.nalloc_index = mols.nr + 1;
1187     snew(mols.index, mols.nalloc_index);
1188
1189     fillMoleculeIndices(mtop, gmx::arrayRefFromArray(mols.index, mols.nr + 1));
1190
1191     return mols;
1192 }
1193
1194 static void gen_t_topology(const gmx_mtop_t &mtop,
1195                            bool              freeEnergyInteractionsAtEnd,
1196                            bool              bMergeConstr,
1197                            t_topology       *top)
1198 {
1199     copyAtomtypesFromMtop(mtop, &top->atomtypes);
1200     copyIdefFromMtop(mtop, &top->idef, freeEnergyInteractionsAtEnd, bMergeConstr);
1201     copyCgsFromMtop(mtop, &top->cgs);
1202     copyExclsFromMtop(mtop, &top->excls);
1203
1204     top->name                        = mtop.name;
1205     top->atoms                       = gmx_mtop_global_atoms(&mtop);
1206     top->mols                        = gmx_mtop_molecules_t_block(mtop);
1207     top->bIntermolecularInteractions = mtop.bIntermolecularInteractions;
1208     top->symtab                      = mtop.symtab;
1209 }
1210
1211 t_topology gmx_mtop_t_to_t_topology(gmx_mtop_t *mtop, bool freeMTop)
1212 {
1213     t_topology     top;
1214
1215     gen_t_topology(*mtop, false, false, &top);
1216
1217     if (freeMTop)
1218     {
1219         // Clear pointers and counts, such that the pointers copied to top
1220         // keep pointing to valid data after destroying mtop.
1221         mtop->symtab.symbuf     = nullptr;
1222         mtop->symtab.nr         = 0;
1223     }
1224     return top;
1225 }
1226
1227 std::vector<int> get_atom_index(const gmx_mtop_t *mtop)
1228 {
1229
1230     std::vector<int>             atom_index;
1231     for (const AtomProxy atomP : AtomRange(*mtop))
1232     {
1233         const t_atom &local = atomP.atom();
1234         int           index = atomP.globalAtomNumber();
1235         if (local.ptype == eptAtom)
1236         {
1237             atom_index.push_back(index);
1238         }
1239     }
1240     return atom_index;
1241 }
1242
1243 void convertAtomsToMtop(t_symtab    *symtab,
1244                         char       **name,
1245                         t_atoms     *atoms,
1246                         gmx_mtop_t  *mtop)
1247 {
1248     mtop->symtab                 = *symtab;
1249
1250     mtop->name                   = name;
1251
1252     mtop->moltype.clear();
1253     mtop->moltype.resize(1);
1254     mtop->moltype.back().atoms   = *atoms;
1255
1256     mtop->molblock.resize(1);
1257     mtop->molblock[0].type       = 0;
1258     mtop->molblock[0].nmol       = 1;
1259
1260     mtop->bIntermolecularInteractions = FALSE;
1261
1262     mtop->natoms                 = atoms->nr;
1263
1264     mtop->haveMoleculeIndices    = false;
1265
1266     gmx_mtop_finalize(mtop);
1267 }