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