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, 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 <limits.h>
  40 #include <stddef.h>
  41 #include <stdio.h>
  42 #include <stdlib.h>
  43 #include <string.h>
  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/fatalerror.h"
  53 #include "gromacs/utility/real.h"
  54 #include "gromacs/utility/smalloc.h"
  55
  56 static int gmx_mtop_maxresnr(const gmx_mtop_t *mtop, int maxres_renum)
  57 {
  58     int            maxresnr, mt, r;
  59     const t_atoms *atoms;
  60
  61     maxresnr = 0;
  62
  63     for (mt = 0; mt < mtop->nmoltype; mt++)
  64     {
  65         atoms = &mtop->moltype[mt].atoms;
  66         if (atoms->nres > maxres_renum)
  67         {
  68             for (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 finalizeMolblocks(gmx_mtop_t *mtop)
  82 {
  83     int atomIndex          = 0;
  84     int residueIndex       = 0;
  85     int residueNumberStart = mtop->maxresnr + 1;
  86     for (int mb = 0; mb < mtop->nmolblock; mb++)
  87     {
  88         gmx_molblock_t *molb          = &mtop->molblock[mb];
  89         int             numResPerMol  = mtop->moltype[molb->type].atoms.nres;
  90         molb->globalAtomStart         = atomIndex;
  91         molb->globalResidueStart      = residueIndex;
  92         atomIndex                    += molb->nmol*molb->natoms_mol;
  93         residueIndex                 += molb->nmol*numResPerMol;
  94         molb->globalAtomEnd           = atomIndex;
  95         molb->residueNumberStart      = residueNumberStart;
  96         if (numResPerMol <= mtop->maxres_renum)
  97         {
  98             residueNumberStart       += molb->nmol*numResPerMol;
  99         }
 100     }
 101 }
 102
 103 void gmx_mtop_finalize(gmx_mtop_t *mtop)
 104 {
 105     char *env;
 106
 107     if (mtop->nmolblock == 1 && mtop->molblock[0].nmol == 1)
 108     {
 109         /* We have a single molecule only, no renumbering needed.
 110          * This case also covers an mtop converted from pdb/gro/... input,
 111          * so we retain the original residue numbering.
 112          */
 113         mtop->maxres_renum = 0;
 114     }
 115     else
 116     {
 117         /* We only renumber single residue molecules. Their intra-molecular
 118          * residue numbering is anyhow irrelevant.
 119          */
 120         mtop->maxres_renum = 1;
 121     }
 122
 123     env = getenv("GMX_MAXRESRENUM");
 124     if (env != nullptr)
 125     {
 126         sscanf(env, "%d", &mtop->maxres_renum);
 127     }
 128     if (mtop->maxres_renum == -1)
 129     {
 130         /* -1 signals renumber residues in all molecules */
 131         mtop->maxres_renum = INT_MAX;
 132     }
 133
 134     mtop->maxresnr = gmx_mtop_maxresnr(mtop, mtop->maxres_renum);
 135
 136     finalizeMolblocks(mtop);
 137 }
 138
 139 void gmx_mtop_count_atomtypes(const gmx_mtop_t *mtop, int state, int typecount[])
 140 {
 141     int      i, mb, nmol, tpi;
 142     t_atoms *atoms;
 143
 144     for (i = 0; i < mtop->ffparams.atnr; ++i)
 145     {
 146         typecount[i] = 0;
 147     }
 148     for (mb = 0; mb < mtop->nmolblock; ++mb)
 149     {
 150         nmol  = mtop->molblock[mb].nmol;
 151         atoms = &mtop->moltype[mtop->molblock[mb].type].atoms;
 152         for (i = 0; i < atoms->nr; ++i)
 153         {
 154             if (state == 0)
 155             {
 156                 tpi = atoms->atom[i].type;
 157             }
 158             else
 159             {
 160                 tpi = atoms->atom[i].typeB;
 161             }
 162             typecount[tpi] += nmol;
 163         }
 164     }
 165 }
 166
 167 int ncg_mtop(const gmx_mtop_t *mtop)
 168 {
 169     int ncg;
 170     int mb;
 171
 172     ncg = 0;
 173     for (mb = 0; mb < mtop->nmolblock; mb++)
 174     {
 175         ncg +=
 176             mtop->molblock[mb].nmol*
 177             mtop->moltype[mtop->molblock[mb].type].cgs.nr;
 178     }
 179
 180     return ncg;
 181 }
 182
 183 int gmx_mtop_nres(const gmx_mtop_t *mtop)
 184 {
 185     int nres = 0;
 186     for (int mb = 0; mb < mtop->nmolblock; ++mb)
 187     {
 188         nres +=
 189             mtop->molblock[mb].nmol*
 190             mtop->moltype[mtop->molblock[mb].type].atoms.nres;
 191     }
 192     return nres;
 193 }
 194
 195 void gmx_mtop_remove_chargegroups(gmx_mtop_t *mtop)
 196 {
 197     int      mt;
 198     t_block *cgs;
 199     int      i;
 200
 201     for (mt = 0; mt < mtop->nmoltype; mt++)
 202     {
 203         cgs = &mtop->moltype[mt].cgs;
 204         if (cgs->nr < mtop->moltype[mt].atoms.nr)
 205         {
 206             cgs->nr = mtop->moltype[mt].atoms.nr;
 207             srenew(cgs->index, cgs->nr+1);
 208             for (i = 0; i < cgs->nr+1; i++)
 209             {
 210                 cgs->index[i] = i;
 211             }
 212         }
 213     }
 214 }
 215
 216 typedef struct gmx_mtop_atomloop_all
 217 {
 218     const gmx_mtop_t *mtop;
 219     int               mblock;
 220     t_atoms          *atoms;
 221     int               mol;
 222     int               maxresnr;
 223     int               at_local;
 224     int               at_global;
 225 } t_gmx_mtop_atomloop_all;
 226
 227 gmx_mtop_atomloop_all_t
 228 gmx_mtop_atomloop_all_init(const gmx_mtop_t *mtop)
 229 {
 230     struct gmx_mtop_atomloop_all *aloop;
 231
 232     snew(aloop, 1);
 233
 234     aloop->mtop         = mtop;
 235     aloop->mblock       = 0;
 236     aloop->atoms        =
 237         &mtop->moltype[mtop->molblock[aloop->mblock].type].atoms;
 238     aloop->mol          = 0;
 239     aloop->maxresnr     = mtop->maxresnr;
 240     aloop->at_local     = -1;
 241     aloop->at_global    = -1;
 242
 243     return aloop;
 244 }
 245
 246 static void gmx_mtop_atomloop_all_destroy(gmx_mtop_atomloop_all_t aloop)
 247 {
 248     sfree(aloop);
 249 }
 250
 251 gmx_bool gmx_mtop_atomloop_all_next(gmx_mtop_atomloop_all_t aloop,
 252                                     int *at_global, const t_atom **atom)
 253 {
 254     if (aloop == nullptr)
 255     {
 256         gmx_incons("gmx_mtop_atomloop_all_next called without calling gmx_mtop_atomloop_all_init");
 257     }
 258
 259     aloop->at_local++;
 260     aloop->at_global++;
 261
 262     if (aloop->at_local >= aloop->atoms->nr)
 263     {
 264         if (aloop->atoms->nres <= aloop->mtop->maxres_renum)
 265         {
 266             /* Single residue molecule, increase the count with one */
 267             aloop->maxresnr += aloop->atoms->nres;
 268         }
 269         aloop->mol++;
 270         aloop->at_local = 0;
 271         if (aloop->mol >= aloop->mtop->molblock[aloop->mblock].nmol)
 272         {
 273             aloop->mblock++;
 274             if (aloop->mblock >= aloop->mtop->nmolblock)
 275             {
 276                 gmx_mtop_atomloop_all_destroy(aloop);
 277                 return FALSE;
 278             }
 279             aloop->atoms = &aloop->mtop->moltype[aloop->mtop->molblock[aloop->mblock].type].atoms;
 280             aloop->mol   = 0;
 281         }
 282     }
 283
 284     *at_global = aloop->at_global;
 285     *atom      = &aloop->atoms->atom[aloop->at_local];
 286
 287     return TRUE;
 288 }
 289
 290 void gmx_mtop_atomloop_all_names(gmx_mtop_atomloop_all_t aloop,
 291                                  char **atomname, int *resnr, char **resname)
 292 {
 293     int resind_mol;
 294
 295     *atomname  = *(aloop->atoms->atomname[aloop->at_local]);
 296     resind_mol = aloop->atoms->atom[aloop->at_local].resind;
 297     *resnr     = aloop->atoms->resinfo[resind_mol].nr;
 298     if (aloop->atoms->nres <= aloop->mtop->maxres_renum)
 299     {
 300         *resnr = aloop->maxresnr + 1 + resind_mol;
 301     }
 302     *resname  = *(aloop->atoms->resinfo[resind_mol].name);
 303 }
 304
 305 void gmx_mtop_atomloop_all_moltype(gmx_mtop_atomloop_all_t aloop,
 306                                    gmx_moltype_t **moltype, int *at_mol)
 307 {
 308     *moltype = &aloop->mtop->moltype[aloop->mtop->molblock[aloop->mblock].type];
 309     *at_mol  = aloop->at_local;
 310 }
 311
 312 typedef struct gmx_mtop_atomloop_block
 313 {
 314     const gmx_mtop_t *mtop;
 315     int               mblock;
 316     t_atoms          *atoms;
 317     int               at_local;
 318 } t_gmx_mtop_atomloop_block;
 319
 320 gmx_mtop_atomloop_block_t
 321 gmx_mtop_atomloop_block_init(const gmx_mtop_t *mtop)
 322 {
 323     struct gmx_mtop_atomloop_block *aloop;
 324
 325     snew(aloop, 1);
 326
 327     aloop->mtop      = mtop;
 328     aloop->mblock    = 0;
 329     aloop->atoms     = &mtop->moltype[mtop->molblock[aloop->mblock].type].atoms;
 330     aloop->at_local  = -1;
 331
 332     return aloop;
 333 }
 334
 335 static void gmx_mtop_atomloop_block_destroy(gmx_mtop_atomloop_block_t aloop)
 336 {
 337     sfree(aloop);
 338 }
 339
 340 gmx_bool gmx_mtop_atomloop_block_next(gmx_mtop_atomloop_block_t aloop,
 341                                       const t_atom **atom, int *nmol)
 342 {
 343     if (aloop == nullptr)
 344     {
 345         gmx_incons("gmx_mtop_atomloop_all_next called without calling gmx_mtop_atomloop_all_init");
 346     }
 347
 348     aloop->at_local++;
 349
 350     if (aloop->at_local >= aloop->atoms->nr)
 351     {
 352         aloop->mblock++;
 353         if (aloop->mblock >= aloop->mtop->nmolblock)
 354         {
 355             gmx_mtop_atomloop_block_destroy(aloop);
 356             return FALSE;
 357         }
 358         aloop->atoms    = &aloop->mtop->moltype[aloop->mtop->molblock[aloop->mblock].type].atoms;
 359         aloop->at_local = 0;
 360     }
 361
 362     *atom = &aloop->atoms->atom[aloop->at_local];
 363     *nmol = aloop->mtop->molblock[aloop->mblock].nmol;
 364
 365     return TRUE;
 366 }
 367
 368 typedef struct gmx_mtop_ilistloop
 369 {
 370     const gmx_mtop_t *mtop;
 371     int               mblock;
 372 } t_gmx_mtop_ilist;
 373
 374 gmx_mtop_ilistloop_t
 375 gmx_mtop_ilistloop_init(const gmx_mtop_t *mtop)
 376 {
 377     struct gmx_mtop_ilistloop *iloop;
 378
 379     snew(iloop, 1);
 380
 381     iloop->mtop      = mtop;
 382     iloop->mblock    = -1;
 383
 384     return iloop;
 385 }
 386
 387 static void gmx_mtop_ilistloop_destroy(gmx_mtop_ilistloop_t iloop)
 388 {
 389     sfree(iloop);
 390 }
 391
 392 gmx_bool gmx_mtop_ilistloop_next(gmx_mtop_ilistloop_t iloop,
 393                                  t_ilist **ilist_mol, int *nmol)
 394 {
 395     if (iloop == nullptr)
 396     {
 397         gmx_incons("gmx_mtop_ilistloop_next called without calling gmx_mtop_ilistloop_init");
 398     }
 399
 400     iloop->mblock++;
 401     if (iloop->mblock >= iloop->mtop->nmolblock)
 402     {
 403         if (iloop->mblock == iloop->mtop->nmolblock &&
 404             iloop->mtop->bIntermolecularInteractions)
 405         {
 406             *ilist_mol = iloop->mtop->intermolecular_ilist;
 407             *nmol      = 1;
 408             return TRUE;
 409         }
 410
 411         gmx_mtop_ilistloop_destroy(iloop);
 412         return FALSE;
 413     }
 414
 415     *ilist_mol =
 416         iloop->mtop->moltype[iloop->mtop->molblock[iloop->mblock].type].ilist;
 417
 418     *nmol = iloop->mtop->molblock[iloop->mblock].nmol;
 419
 420     return TRUE;
 421 }
 422 typedef struct gmx_mtop_ilistloop_all
 423 {
 424     const gmx_mtop_t *mtop;
 425     int               mblock;
 426     int               mol;
 427     int               a_offset;
 428 } t_gmx_mtop_ilist_all;
 429
 430 gmx_mtop_ilistloop_all_t
 431 gmx_mtop_ilistloop_all_init(const gmx_mtop_t *mtop)
 432 {
 433     struct gmx_mtop_ilistloop_all *iloop;
 434
 435     snew(iloop, 1);
 436
 437     iloop->mtop      = mtop;
 438     iloop->mblock    = 0;
 439     iloop->mol       = -1;
 440     iloop->a_offset  = 0;
 441
 442     return iloop;
 443 }
 444
 445 static void gmx_mtop_ilistloop_all_destroy(gmx_mtop_ilistloop_all_t iloop)
 446 {
 447     sfree(iloop);
 448 }
 449
 450 gmx_bool gmx_mtop_ilistloop_all_next(gmx_mtop_ilistloop_all_t iloop,
 451                                      t_ilist **ilist_mol, int *atnr_offset)
 452 {
 453
 454     if (iloop == nullptr)
 455     {
 456         gmx_incons("gmx_mtop_ilistloop_all_next called without calling gmx_mtop_ilistloop_all_init");
 457     }
 458
 459     if (iloop->mol >= 0)
 460     {
 461         iloop->a_offset += iloop->mtop->molblock[iloop->mblock].natoms_mol;
 462     }
 463
 464     iloop->mol++;
 465
 466     /* Inter-molecular interactions, if present, are indexed with
 467      * iloop->mblock == iloop->mtop->nmolblock, thus we should separately
 468      * check for this value in this conditional.
 469      */
 470     if (iloop->mblock == iloop->mtop->nmolblock ||
 471         iloop->mol >= iloop->mtop->molblock[iloop->mblock].nmol)
 472     {
 473         iloop->mblock++;
 474         iloop->mol = 0;
 475         if (iloop->mblock >= iloop->mtop->nmolblock)
 476         {
 477             if (iloop->mblock == iloop->mtop->nmolblock &&
 478                 iloop->mtop->bIntermolecularInteractions)
 479             {
 480                 *ilist_mol   = iloop->mtop->intermolecular_ilist;
 481                 *atnr_offset = 0;
 482                 return TRUE;
 483             }
 484
 485             gmx_mtop_ilistloop_all_destroy(iloop);
 486             return FALSE;
 487         }
 488     }
 489
 490     *ilist_mol =
 491         iloop->mtop->moltype[iloop->mtop->molblock[iloop->mblock].type].ilist;
 492
 493     *atnr_offset = iloop->a_offset;
 494
 495     return TRUE;
 496 }
 497
 498 int gmx_mtop_ftype_count(const gmx_mtop_t *mtop, int ftype)
 499 {
 500     gmx_mtop_ilistloop_t iloop;
 501     t_ilist             *il;
 502     int                  n, nmol;
 503
 504     n = 0;
 505
 506     iloop = gmx_mtop_ilistloop_init(mtop);
 507     while (gmx_mtop_ilistloop_next(iloop, &il, &nmol))
 508     {
 509         n += nmol*il[ftype].nr/(1+NRAL(ftype));
 510     }
 511
 512     if (mtop->bIntermolecularInteractions)
 513     {
 514         n += mtop->intermolecular_ilist[ftype].nr/(1+NRAL(ftype));
 515     }
 516
 517     return n;
 518 }
 519
 520 t_block gmx_mtop_global_cgs(const gmx_mtop_t *mtop)
 521 {
 522     t_block         cgs_gl, *cgs_mol;
 523     int             mb, mol, cg;
 524     gmx_molblock_t *molb;
 525
 526     /* In most cases this is too much, but we realloc at the end */
 527     snew(cgs_gl.index, mtop->natoms+1);
 528
 529     cgs_gl.nr       = 0;
 530     cgs_gl.index[0] = 0;
 531     for (mb = 0; mb < mtop->nmolblock; mb++)
 532     {
 533         molb    = &mtop->molblock[mb];
 534         cgs_mol = &mtop->moltype[molb->type].cgs;
 535         for (mol = 0; mol < molb->nmol; mol++)
 536         {
 537             for (cg = 0; cg < cgs_mol->nr; cg++)
 538             {
 539                 cgs_gl.index[cgs_gl.nr+1] =
 540                     cgs_gl.index[cgs_gl.nr] +
 541                     cgs_mol->index[cg+1] - cgs_mol->index[cg];
 542                 cgs_gl.nr++;
 543             }
 544         }
 545     }
 546     cgs_gl.nalloc_index = cgs_gl.nr + 1;
 547     srenew(cgs_gl.index, cgs_gl.nalloc_index);
 548
 549     return cgs_gl;
 550 }
 551
 552 static void atomcat(t_atoms *dest, t_atoms *src, int copies,
 553                     int maxres_renum, int *maxresnr)
 554 {
 555     int i, j, l, size;
 556     int srcnr  = src->nr;
 557     int destnr = dest->nr;
 558
 559     if (dest->nr == 0)
 560     {
 561         dest->haveMass    = src->haveMass;
 562         dest->haveType    = src->haveType;
 563         dest->haveCharge  = src->haveCharge;
 564         dest->haveBState  = src->haveBState;
 565         dest->havePdbInfo = src->havePdbInfo;
 566     }
 567     else
 568     {
 569         dest->haveMass    = dest->haveMass    && src->haveMass;
 570         dest->haveType    = dest->haveType    && src->haveType;
 571         dest->haveCharge  = dest->haveCharge  && src->haveCharge;
 572         dest->haveBState  = dest->haveBState  && src->haveBState;
 573         dest->havePdbInfo = dest->havePdbInfo && src->havePdbInfo;
 574     }
 575
 576     if (srcnr)
 577     {
 578         size = destnr+copies*srcnr;
 579         srenew(dest->atom, size);
 580         srenew(dest->atomname, size);
 581         if (dest->haveType)
 582         {
 583             srenew(dest->atomtype, size);
 584             if (dest->haveBState)
 585             {
 586                 srenew(dest->atomtypeB, size);
 587             }
 588         }
 589         if (dest->havePdbInfo)
 590         {
 591             srenew(dest->pdbinfo, size);
 592         }
 593     }
 594     if (src->nres)
 595     {
 596         size = dest->nres+copies*src->nres;
 597         srenew(dest->resinfo, size);
 598     }
 599
 600     /* residue information */
 601     for (l = dest->nres, j = 0; (j < copies); j++, l += src->nres)
 602     {
 603         memcpy((char *) &(dest->resinfo[l]), (char *) &(src->resinfo[0]),
 604                (size_t)(src->nres*sizeof(src->resinfo[0])));
 605     }
 606
 607     for (l = destnr, j = 0; (j < copies); j++, l += srcnr)
 608     {
 609         memcpy((char *) &(dest->atom[l]), (char *) &(src->atom[0]),
 610                (size_t)(srcnr*sizeof(src->atom[0])));
 611         memcpy((char *) &(dest->atomname[l]), (char *) &(src->atomname[0]),
 612                (size_t)(srcnr*sizeof(src->atomname[0])));
 613         if (dest->haveType)
 614         {
 615             memcpy((char *) &(dest->atomtype[l]), (char *) &(src->atomtype[0]),
 616                    (size_t)(srcnr*sizeof(src->atomtype[0])));
 617             if (dest->haveBState)
 618             {
 619                 memcpy((char *) &(dest->atomtypeB[l]), (char *) &(src->atomtypeB[0]),
 620                        (size_t)(srcnr*sizeof(src->atomtypeB[0])));
 621             }
 622         }
 623         if (dest->havePdbInfo)
 624         {
 625             memcpy((char *) &(dest->pdbinfo[l]), (char *) &(src->pdbinfo[0]),
 626                    (size_t)(srcnr*sizeof(src->pdbinfo[0])));
 627         }
 628     }
 629
 630     /* Increment residue indices */
 631     for (l = destnr, j = 0; (j < copies); j++)
 632     {
 633         for (i = 0; (i < srcnr); i++, l++)
 634         {
 635             dest->atom[l].resind = dest->nres+j*src->nres+src->atom[i].resind;
 636         }
 637     }
 638
 639     if (src->nres <= maxres_renum)
 640     {
 641         /* Single residue molecule, continue counting residues */
 642         for (j = 0; (j < copies); j++)
 643         {
 644             for (l = 0; l < src->nres; l++)
 645             {
 646                 (*maxresnr)++;
 647                 dest->resinfo[dest->nres+j*src->nres+l].nr = *maxresnr;
 648             }
 649         }
 650     }
 651
 652     dest->nres += copies*src->nres;
 653     dest->nr   += copies*src->nr;
 654 }
 655
 656 t_atoms gmx_mtop_global_atoms(const gmx_mtop_t *mtop)
 657 {
 658     t_atoms         atoms;
 659     int             maxresnr, mb;
 660     gmx_molblock_t *molb;
 661
 662     init_t_atoms(&atoms, 0, FALSE);
 663
 664     maxresnr = mtop->maxresnr;
 665     for (mb = 0; mb < mtop->nmolblock; mb++)
 666     {
 667         molb = &mtop->molblock[mb];
 668         atomcat(&atoms, &mtop->moltype[molb->type].atoms, molb->nmol,
 669                 mtop->maxres_renum, &maxresnr);
 670     }
 671
 672     return atoms;
 673 }
 674
 675 /*
 676  * The cat routines below are old code from src/kernel/topcat.c
 677  */
 678
 679 static void blockcat(t_block *dest, t_block *src, int copies)
 680 {
 681     int i, j, l, nra, size;
 682
 683     if (src->nr)
 684     {
 685         size = (dest->nr+copies*src->nr+1);
 686         srenew(dest->index, size);
 687     }
 688
 689     nra = dest->index[dest->nr];
 690     for (l = dest->nr, j = 0; (j < copies); j++)
 691     {
 692         for (i = 0; (i < src->nr); i++)
 693         {
 694             dest->index[l++] = nra + src->index[i];
 695         }
 696         nra += src->index[src->nr];
 697     }
 698     dest->nr             += copies*src->nr;
 699     dest->index[dest->nr] = nra;
 700 }
 701
 702 static void blockacat(t_blocka *dest, t_blocka *src, int copies,
 703                       int dnum, int snum)
 704 {
 705     int i, j, l, size;
 706     int destnr  = dest->nr;
 707     int destnra = dest->nra;
 708
 709     if (src->nr)
 710     {
 711         size = (dest->nr+copies*src->nr+1);
 712         srenew(dest->index, size);
 713     }
 714     if (src->nra)
 715     {
 716         size = (dest->nra+copies*src->nra);
 717         srenew(dest->a, size);
 718     }
 719
 720     for (l = destnr, j = 0; (j < copies); j++)
 721     {
 722         for (i = 0; (i < src->nr); i++)
 723         {
 724             dest->index[l++] = dest->nra+src->index[i];
 725         }
 726         dest->nra += src->nra;
 727     }
 728     for (l = destnra, j = 0; (j < copies); j++)
 729     {
 730         for (i = 0; (i < src->nra); i++)
 731         {
 732             dest->a[l++] = dnum+src->a[i];
 733         }
 734         dnum     += snum;
 735         dest->nr += src->nr;
 736     }
 737     dest->index[dest->nr] = dest->nra;
 738 }
 739
 740 static void ilistcat(int ftype, t_ilist *dest, t_ilist *src, int copies,
 741                      int dnum, int snum)
 742 {
 743     int nral, c, i, a;
 744
 745     nral = NRAL(ftype);
 746
 747     dest->nalloc = dest->nr + copies*src->nr;
 748     srenew(dest->iatoms, dest->nalloc);
 749
 750     for (c = 0; c < copies; c++)
 751     {
 752         for (i = 0; i < src->nr; )
 753         {
 754             dest->iatoms[dest->nr++] = src->iatoms[i++];
 755             for (a = 0; a < nral; a++)
 756             {
 757                 dest->iatoms[dest->nr++] = dnum + src->iatoms[i++];
 758             }
 759         }
 760         dnum += snum;
 761     }
 762 }
 763
 764 static void set_posres_params(t_idef *idef, gmx_molblock_t *molb,
 765                               int i0, int a_offset)
 766 {
 767     t_ilist   *il;
 768     int        i1, i, a_molb;
 769     t_iparams *ip;
 770
 771     il = &idef->il[F_POSRES];
 772     i1 = il->nr/2;
 773     idef->iparams_posres_nalloc = i1;
 774     srenew(idef->iparams_posres, idef->iparams_posres_nalloc);
 775     for (i = i0; i < i1; i++)
 776     {
 777         ip = &idef->iparams_posres[i];
 778         /* Copy the force constants */
 779         *ip    = idef->iparams[il->iatoms[i*2]];
 780         a_molb = il->iatoms[i*2+1] - a_offset;
 781         if (molb->nposres_xA == 0)
 782         {
 783             gmx_incons("Position restraint coordinates are missing");
 784         }
 785         ip->posres.pos0A[XX] = molb->posres_xA[a_molb][XX];
 786         ip->posres.pos0A[YY] = molb->posres_xA[a_molb][YY];
 787         ip->posres.pos0A[ZZ] = molb->posres_xA[a_molb][ZZ];
 788         if (molb->nposres_xB > 0)
 789         {
 790             ip->posres.pos0B[XX] = molb->posres_xB[a_molb][XX];
 791             ip->posres.pos0B[YY] = molb->posres_xB[a_molb][YY];
 792             ip->posres.pos0B[ZZ] = molb->posres_xB[a_molb][ZZ];
 793         }
 794         else
 795         {
 796             ip->posres.pos0B[XX] = ip->posres.pos0A[XX];
 797             ip->posres.pos0B[YY] = ip->posres.pos0A[YY];
 798             ip->posres.pos0B[ZZ] = ip->posres.pos0A[ZZ];
 799         }
 800         /* Set the parameter index for idef->iparams_posre */
 801         il->iatoms[i*2] = i;
 802     }
 803 }
 804
 805 static void set_fbposres_params(t_idef *idef, gmx_molblock_t *molb,
 806                                 int i0, int a_offset)
 807 {
 808     t_ilist   *il;
 809     int        i1, i, a_molb;
 810     t_iparams *ip;
 811
 812     il = &idef->il[F_FBPOSRES];
 813     i1 = il->nr/2;
 814     idef->iparams_fbposres_nalloc = i1;
 815     srenew(idef->iparams_fbposres, idef->iparams_fbposres_nalloc);
 816     for (i = i0; i < i1; i++)
 817     {
 818         ip = &idef->iparams_fbposres[i];
 819         /* Copy the force constants */
 820         *ip    = idef->iparams[il->iatoms[i*2]];
 821         a_molb = il->iatoms[i*2+1] - a_offset;
 822         if (molb->nposres_xA == 0)
 823         {
 824             gmx_incons("Position restraint coordinates are missing");
 825         }
 826         /* Take flat-bottom posres reference from normal position restraints */
 827         ip->fbposres.pos0[XX] = molb->posres_xA[a_molb][XX];
 828         ip->fbposres.pos0[YY] = molb->posres_xA[a_molb][YY];
 829         ip->fbposres.pos0[ZZ] = molb->posres_xA[a_molb][ZZ];
 830         /* Note: no B-type for flat-bottom posres */
 831
 832         /* Set the parameter index for idef->iparams_posre */
 833         il->iatoms[i*2] = i;
 834     }
 835 }
 836
 837 static void gen_local_top(const gmx_mtop_t *mtop,
 838                           bool              freeEnergyInteractionsAtEnd,
 839                           bool              bMergeConstr,
 840                           gmx_localtop_t   *top)
 841 {
 842     int                     mb, srcnr, destnr, ftype, natoms, mol, nposre_old, nfbposre_old;
 843     gmx_molblock_t         *molb;
 844     gmx_moltype_t          *molt;
 845     const gmx_ffparams_t   *ffp;
 846     t_idef                 *idef;
 847     real                   *qA, *qB;
 848     gmx_mtop_atomloop_all_t aloop;
 849     int                     ag;
 850
 851     top->atomtypes = mtop->atomtypes;
 852
 853     ffp = &mtop->ffparams;
 854
 855     idef                          = &top->idef;
 856     idef->ntypes                  = ffp->ntypes;
 857     idef->atnr                    = ffp->atnr;
 858     idef->functype                = ffp->functype;
 859     idef->iparams                 = ffp->iparams;
 860     idef->iparams_posres          = nullptr;
 861     idef->iparams_posres_nalloc   = 0;
 862     idef->iparams_fbposres        = nullptr;
 863     idef->iparams_fbposres_nalloc = 0;
 864     idef->fudgeQQ                 = ffp->fudgeQQ;
 865     idef->cmap_grid               = ffp->cmap_grid;
 866     idef->ilsort                  = ilsortUNKNOWN;
 867
 868     init_block(&top->cgs);
 869     init_blocka(&top->excls);
 870     for (ftype = 0; ftype < F_NRE; ftype++)
 871     {
 872         idef->il[ftype].nr     = 0;
 873         idef->il[ftype].nalloc = 0;
 874         idef->il[ftype].iatoms = nullptr;
 875     }
 876
 877     natoms = 0;
 878     for (mb = 0; mb < mtop->nmolblock; mb++)
 879     {
 880         molb = &mtop->molblock[mb];
 881         molt = &mtop->moltype[molb->type];
 882
 883         srcnr  = molt->atoms.nr;
 884         destnr = natoms;
 885
 886         blockcat(&top->cgs, &molt->cgs, molb->nmol);
 887
 888         blockacat(&top->excls, &molt->excls, molb->nmol, destnr, srcnr);
 889
 890         nposre_old   = idef->il[F_POSRES].nr;
 891         nfbposre_old = idef->il[F_FBPOSRES].nr;
 892         for (ftype = 0; ftype < F_NRE; ftype++)
 893         {
 894             if (bMergeConstr &&
 895                 ftype == F_CONSTR && molt->ilist[F_CONSTRNC].nr > 0)
 896             {
 897                 /* Merge all constrains into one ilist.
 898                  * This simplifies the constraint code.
 899                  */
 900                 for (mol = 0; mol < molb->nmol; mol++)
 901                 {
 902                     ilistcat(ftype, &idef->il[F_CONSTR], &molt->ilist[F_CONSTR],
 903                              1, destnr+mol*srcnr, srcnr);
 904                     ilistcat(ftype, &idef->il[F_CONSTR], &molt->ilist[F_CONSTRNC],
 905                              1, destnr+mol*srcnr, srcnr);
 906                 }
 907             }
 908             else if (!(bMergeConstr && ftype == F_CONSTRNC))
 909             {
 910                 ilistcat(ftype, &idef->il[ftype], &molt->ilist[ftype],
 911                          molb->nmol, destnr, srcnr);
 912             }
 913         }
 914         if (idef->il[F_POSRES].nr > nposre_old)
 915         {
 916             /* Executing this line line stops gmxdump -sys working
 917              * correctly. I'm not aware there's an elegant fix. */
 918             set_posres_params(idef, molb, nposre_old/2, natoms);
 919         }
 920         if (idef->il[F_FBPOSRES].nr > nfbposre_old)
 921         {
 922             set_fbposres_params(idef, molb, nfbposre_old/2, natoms);
 923         }
 924
 925         natoms += molb->nmol*srcnr;
 926     }
 927
 928     if (mtop->bIntermolecularInteractions)
 929     {
 930         for (ftype = 0; ftype < F_NRE; ftype++)
 931         {
 932             ilistcat(ftype, &idef->il[ftype], &mtop->intermolecular_ilist[ftype],
 933                      1, 0, mtop->natoms);
 934         }
 935     }
 936
 937     if (freeEnergyInteractionsAtEnd && gmx_mtop_bondeds_free_energy(mtop))
 938     {
 939         snew(qA, mtop->natoms);
 940         snew(qB, mtop->natoms);
 941         aloop = gmx_mtop_atomloop_all_init(mtop);
 942         const t_atom *atom;
 943         while (gmx_mtop_atomloop_all_next(aloop, &ag, &atom))
 944         {
 945             qA[ag] = atom->q;
 946             qB[ag] = atom->qB;
 947         }
 948         gmx_sort_ilist_fe(&top->idef, qA, qB);
 949         sfree(qA);
 950         sfree(qB);
 951     }
 952     else
 953     {
 954         top->idef.ilsort = ilsortNO_FE;
 955     }
 956 }
 957
 958 gmx_localtop_t *
 959 gmx_mtop_generate_local_top(const gmx_mtop_t *mtop,
 960                             bool              freeEnergyInteractionsAtEnd)
 961 {
 962     gmx_localtop_t *top;
 963
 964     snew(top, 1);
 965
 966     gen_local_top(mtop, freeEnergyInteractionsAtEnd, true, top);
 967
 968     return top;
 969 }
 970
 971 t_topology gmx_mtop_t_to_t_topology(gmx_mtop_t *mtop, bool freeMTop)
 972 {
 973     int            mt, mb;
 974     gmx_localtop_t ltop;
 975     t_topology     top;
 976
 977     gen_local_top(mtop, false, FALSE, &ltop);
 978     ltop.idef.ilsort = ilsortUNKNOWN;
 979
 980     top.name                        = mtop->name;
 981     top.idef                        = ltop.idef;
 982     top.atomtypes                   = ltop.atomtypes;
 983     top.cgs                         = ltop.cgs;
 984     top.excls                       = ltop.excls;
 985     top.atoms                       = gmx_mtop_global_atoms(mtop);
 986     top.mols                        = mtop->mols;
 987     top.bIntermolecularInteractions = mtop->bIntermolecularInteractions;
 988     top.symtab                      = mtop->symtab;
 989
 990     if (freeMTop)
 991     {
 992         // Free pointers that have not been copied to top.
 993         for (mt = 0; mt < mtop->nmoltype; mt++)
 994         {
 995             done_moltype(&mtop->moltype[mt]);
 996         }
 997         sfree(mtop->moltype);
 998
 999         for (mb = 0; mb < mtop->nmolblock; mb++)
1000         {
1001             done_molblock(&mtop->molblock[mb]);
1002         }
1003         sfree(mtop->molblock);
1004
1005         done_gmx_groups_t(&mtop->groups);
1006     }
1007
1008     return top;
1009 }
1010
1011 std::vector<size_t> get_atom_index(const gmx_mtop_t *mtop)
1012 {
1013
1014     std::vector<size_t>       atom_index;
1015     gmx_mtop_atomloop_block_t aloopb = gmx_mtop_atomloop_block_init(mtop);
1016     const t_atom             *atom;
1017     int                       nmol, j = 0;
1018     while (gmx_mtop_atomloop_block_next(aloopb, &atom, &nmol))
1019     {
1020         if (atom->ptype == eptAtom)
1021         {
1022             atom_index.push_back(j);
1023         }
1024         j++;
1025     }
1026     return atom_index;
1027 }
1028
1029 void convertAtomsToMtop(t_symtab    *symtab,
1030                         char       **name,
1031                         t_atoms     *atoms,
1032                         gmx_mtop_t  *mtop)
1033 {
1034     mtop->symtab                 = *symtab;
1035
1036     mtop->name                   = name;
1037
1038     mtop->nmoltype               = 1;
1039     // This snew clears all entries, we should replace it by an initializer
1040     snew(mtop->moltype, mtop->nmoltype);
1041     mtop->moltype[0].atoms       = *atoms;
1042     init_block(&mtop->moltype[0].cgs);
1043     init_blocka(&mtop->moltype[0].excls);
1044
1045     mtop->nmolblock              = 1;
1046     // This snew clears all entries, we should replace it by an initializer
1047     snew(mtop->molblock, mtop->nmolblock);
1048     mtop->molblock[0].type       = 0;
1049     mtop->molblock[0].nmol       = 1;
1050     mtop->molblock[0].natoms_mol = atoms->nr;
1051
1052     mtop->bIntermolecularInteractions = FALSE;
1053
1054     mtop->natoms                 = atoms->nr;
1055
1056     gmx_mtop_finalize(mtop);
1057 }