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