src/gromacs/topology/mtop_util.c

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2008,2009,2010,2012,2013,2014,2015, 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/legacyheaders/types/enums.h"
  46 #include "gromacs/legacyheaders/types/ifunc.h"
  47 #include "gromacs/legacyheaders/types/inputrec.h"
  48 #include "gromacs/legacyheaders/types/simple.h"
  49 #include "gromacs/math/vectypes.h"
  50 #include "gromacs/topology/atoms.h"
  51 #include "gromacs/topology/block.h"
  52 #include "gromacs/topology/idef.h"
  53 #include "gromacs/topology/topology.h"
  54 #include "gromacs/topology/topsort.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, mt, r;
  62     const t_atoms *atoms;
  63
  64     maxresnr = 0;
  65
  66     for (mt = 0; mt < mtop->nmoltype; mt++)
  67     {
  68         atoms = &mtop->moltype[mt].atoms;
  69         if (atoms->nres > maxres_renum)
  70         {
  71             for (r = 0; r < atoms->nres; r++)
  72             {
  73                 if (atoms->resinfo[r].nr > maxresnr)
  74                 {
  75                     maxresnr = atoms->resinfo[r].nr;
  76                 }
  77             }
  78         }
  79     }
  80
  81     return maxresnr;
  82 }
  83
  84 void gmx_mtop_finalize(gmx_mtop_t *mtop)
  85 {
  86     char *env;
  87
  88     mtop->maxres_renum = 1;
  89
  90     env = getenv("GMX_MAXRESRENUM");
  91     if (env != NULL)
  92     {
  93         sscanf(env, "%d", &mtop->maxres_renum);
  94     }
  95     if (mtop->maxres_renum == -1)
  96     {
  97         /* -1 signals renumber residues in all molecules */
  98         mtop->maxres_renum = INT_MAX;
  99     }
 100
 101     mtop->maxresnr = gmx_mtop_maxresnr(mtop, mtop->maxres_renum);
 102 }
 103
 104 void gmx_mtop_count_atomtypes(const gmx_mtop_t *mtop, int state, int typecount[])
 105 {
 106     int      i, mb, nmol, tpi;
 107     t_atoms *atoms;
 108
 109     for (i = 0; i < mtop->ffparams.atnr; ++i)
 110     {
 111         typecount[i] = 0;
 112     }
 113     for (mb = 0; mb < mtop->nmolblock; ++mb)
 114     {
 115         nmol  = mtop->molblock[mb].nmol;
 116         atoms = &mtop->moltype[mtop->molblock[mb].type].atoms;
 117         for (i = 0; i < atoms->nr; ++i)
 118         {
 119             if (state == 0)
 120             {
 121                 tpi = atoms->atom[i].type;
 122             }
 123             else
 124             {
 125                 tpi = atoms->atom[i].typeB;
 126             }
 127             typecount[tpi] += nmol;
 128         }
 129     }
 130 }
 131
 132 int ncg_mtop(const gmx_mtop_t *mtop)
 133 {
 134     int ncg;
 135     int mb;
 136
 137     ncg = 0;
 138     for (mb = 0; mb < mtop->nmolblock; mb++)
 139     {
 140         ncg +=
 141             mtop->molblock[mb].nmol*
 142             mtop->moltype[mtop->molblock[mb].type].cgs.nr;
 143     }
 144
 145     return ncg;
 146 }
 147
 148 void gmx_mtop_remove_chargegroups(gmx_mtop_t *mtop)
 149 {
 150     int      mt;
 151     t_block *cgs;
 152     int      i;
 153
 154     for (mt = 0; mt < mtop->nmoltype; mt++)
 155     {
 156         cgs = &mtop->moltype[mt].cgs;
 157         if (cgs->nr < mtop->moltype[mt].atoms.nr)
 158         {
 159             cgs->nr = mtop->moltype[mt].atoms.nr;
 160             srenew(cgs->index, cgs->nr+1);
 161             for (i = 0; i < cgs->nr+1; i++)
 162             {
 163                 cgs->index[i] = i;
 164             }
 165         }
 166     }
 167 }
 168
 169
 170 typedef struct
 171 {
 172     int a_start;
 173     int a_end;
 174     int na_mol;
 175 } mb_at_t;
 176
 177 typedef struct gmx_mtop_atomlookup
 178 {
 179     const gmx_mtop_t *mtop;
 180     int               nmb;
 181     int               mb_start;
 182     mb_at_t          *mba;
 183 } t_gmx_mtop_atomlookup;
 184
 185
 186 gmx_mtop_atomlookup_t
 187 gmx_mtop_atomlookup_init(const gmx_mtop_t *mtop)
 188 {
 189     t_gmx_mtop_atomlookup *alook;
 190     int                    mb;
 191     int                    a_start, a_end, na, na_start = -1;
 192
 193     snew(alook, 1);
 194
 195     alook->mtop     = mtop;
 196     alook->nmb      = mtop->nmolblock;
 197     alook->mb_start = 0;
 198     snew(alook->mba, alook->nmb);
 199
 200     a_start = 0;
 201     for (mb = 0; mb < mtop->nmolblock; mb++)
 202     {
 203         na    = mtop->molblock[mb].nmol*mtop->molblock[mb].natoms_mol;
 204         a_end = a_start + na;
 205
 206         alook->mba[mb].a_start = a_start;
 207         alook->mba[mb].a_end   = a_end;
 208         alook->mba[mb].na_mol  = mtop->molblock[mb].natoms_mol;
 209
 210         /* We start the binary search with the largest block */
 211         if (mb == 0 || na > na_start)
 212         {
 213             alook->mb_start = mb;
 214             na_start        = na;
 215         }
 216
 217         a_start = a_end;
 218     }
 219
 220     return alook;
 221 }
 222
 223 gmx_mtop_atomlookup_t
 224 gmx_mtop_atomlookup_settle_init(const gmx_mtop_t *mtop)
 225 {
 226     t_gmx_mtop_atomlookup *alook;
 227     int                    mb;
 228     int                    na, na_start = -1;
 229
 230     alook = gmx_mtop_atomlookup_init(mtop);
 231
 232     /* Check if the starting molblock has settle */
 233     if (mtop->moltype[mtop->molblock[alook->mb_start].type].ilist[F_SETTLE].nr  == 0)
 234     {
 235         /* Search the largest molblock with settle */
 236         alook->mb_start = -1;
 237         for (mb = 0; mb < mtop->nmolblock; mb++)
 238         {
 239             if (mtop->moltype[mtop->molblock[mb].type].ilist[F_SETTLE].nr > 0)
 240             {
 241                 na = alook->mba[mb].a_end - alook->mba[mb].a_start;
 242                 if (alook->mb_start == -1 || na > na_start)
 243                 {
 244                     alook->mb_start = mb;
 245                     na_start        = na;
 246                 }
 247             }
 248         }
 249
 250         if (alook->mb_start == -1)
 251         {
 252             gmx_incons("gmx_mtop_atomlookup_settle_init called without settles");
 253         }
 254     }
 255
 256     return alook;
 257 }
 258
 259 void
 260 gmx_mtop_atomlookup_destroy(gmx_mtop_atomlookup_t alook)
 261 {
 262     sfree(alook->mba);
 263     sfree(alook);
 264 }
 265
 266 void gmx_mtop_atomnr_to_atom(const gmx_mtop_atomlookup_t alook,
 267                              int                         atnr_global,
 268                              t_atom                    **atom)
 269 {
 270     int mb0, mb1, mb;
 271     int a_start, atnr_mol;
 272
 273 #ifdef DEBUG_MTOP
 274     if (atnr_global < 0 || atnr_global >= mtop->natoms)
 275     {
 276         gmx_fatal(FARGS, "gmx_mtop_atomnr_to_moltype was called with atnr_global=%d which is not in the atom range of this system (%d-%d)",
 277                   atnr_global, 0, mtop->natoms-1);
 278     }
 279 #endif
 280
 281     mb0 = -1;
 282     mb1 = alook->nmb;
 283     mb  = alook->mb_start;
 284
 285     while (TRUE)
 286     {
 287         a_start = alook->mba[mb].a_start;
 288         if (atnr_global < a_start)
 289         {
 290             mb1 = mb;
 291         }
 292         else if (atnr_global >= alook->mba[mb].a_end)
 293         {
 294             mb0 = mb;
 295         }
 296         else
 297         {
 298             break;
 299         }
 300         mb = ((mb0 + mb1 + 1)>>1);
 301     }
 302
 303     atnr_mol = (atnr_global - a_start) % alook->mba[mb].na_mol;
 304
 305     *atom = &alook->mtop->moltype[alook->mtop->molblock[mb].type].atoms.atom[atnr_mol];
 306 }
 307
 308 void gmx_mtop_atomnr_to_ilist(const gmx_mtop_atomlookup_t alook,
 309                               int atnr_global,
 310                               t_ilist **ilist_mol, int *atnr_offset)
 311 {
 312     int mb0, mb1, mb;
 313     int a_start, atnr_local;
 314
 315 #ifdef DEBUG_MTOP
 316     if (atnr_global < 0 || atnr_global >= mtop->natoms)
 317     {
 318         gmx_fatal(FARGS, "gmx_mtop_atomnr_to_moltype was called with atnr_global=%d which is not in the atom range of this system (%d-%d)",
 319                   atnr_global, 0, mtop->natoms-1);
 320     }
 321 #endif
 322
 323     mb0 = -1;
 324     mb1 = alook->nmb;
 325     mb  = alook->mb_start;
 326
 327     while (TRUE)
 328     {
 329         a_start = alook->mba[mb].a_start;
 330         if (atnr_global < a_start)
 331         {
 332             mb1 = mb;
 333         }
 334         else if (atnr_global >= alook->mba[mb].a_end)
 335         {
 336             mb0 = mb;
 337         }
 338         else
 339         {
 340             break;
 341         }
 342         mb = ((mb0 + mb1 + 1)>>1);
 343     }
 344
 345     *ilist_mol = alook->mtop->moltype[alook->mtop->molblock[mb].type].ilist;
 346
 347     atnr_local = (atnr_global - a_start) % alook->mba[mb].na_mol;
 348
 349     *atnr_offset = atnr_global - atnr_local;
 350 }
 351
 352 void gmx_mtop_atomnr_to_molblock_ind(const gmx_mtop_atomlookup_t alook,
 353                                      int atnr_global,
 354                                      int *molb, int *molnr, int *atnr_mol)
 355 {
 356     int mb0, mb1, mb;
 357     int a_start;
 358
 359 #ifdef DEBUG_MTOP
 360     if (atnr_global < 0 || atnr_global >= mtop->natoms)
 361     {
 362         gmx_fatal(FARGS, "gmx_mtop_atomnr_to_moltype was called with atnr_global=%d which is not in the atom range of this system (%d-%d)",
 363                   atnr_global, 0, mtop->natoms-1);
 364     }
 365 #endif
 366
 367     mb0 = -1;
 368     mb1 = alook->nmb;
 369     mb  = alook->mb_start;
 370
 371     while (TRUE)
 372     {
 373         a_start = alook->mba[mb].a_start;
 374         if (atnr_global < a_start)
 375         {
 376             mb1 = mb;
 377         }
 378         else if (atnr_global >= alook->mba[mb].a_end)
 379         {
 380             mb0 = mb;
 381         }
 382         else
 383         {
 384             break;
 385         }
 386         mb = ((mb0 + mb1 + 1)>>1);
 387     }
 388
 389     *molb     = mb;
 390     *molnr    = (atnr_global - a_start) / alook->mba[mb].na_mol;
 391     *atnr_mol = atnr_global - a_start - (*molnr)*alook->mba[mb].na_mol;
 392 }
 393
 394 void gmx_mtop_atominfo_global(const gmx_mtop_t *mtop, int atnr_global,
 395                               char **atomname, int *resnr, char **resname)
 396 {
 397     int             mb, a_start, a_end, maxresnr, at_loc;
 398     gmx_molblock_t *molb;
 399     t_atoms        *atoms = NULL;
 400
 401     if (atnr_global < 0 || atnr_global >= mtop->natoms)
 402     {
 403         gmx_fatal(FARGS, "gmx_mtop_atominfo_global was called with atnr_global=%d which is not in the atom range of this system (%d-%d)",
 404                   atnr_global, 0, mtop->natoms-1);
 405     }
 406
 407     mb       = -1;
 408     a_end    = 0;
 409     maxresnr = mtop->maxresnr;
 410     do
 411     {
 412         if (mb >= 0)
 413         {
 414             /* cppcheck-suppress nullPointer #6330*/
 415             if (atoms->nres <= mtop->maxres_renum)
 416             {
 417                 /* Single residue molecule, keep counting */
 418                 /* cppcheck-suppress nullPointer #6330*/
 419                 maxresnr += mtop->molblock[mb].nmol*atoms->nres;
 420             }
 421         }
 422         mb++;
 423         atoms   = &mtop->moltype[mtop->molblock[mb].type].atoms;
 424         a_start = a_end;
 425         a_end   = a_start + mtop->molblock[mb].nmol*atoms->nr;
 426     }
 427     while (atnr_global >= a_end);
 428
 429     at_loc    = (atnr_global - a_start) % atoms->nr;
 430     *atomname = *(atoms->atomname[at_loc]);
 431     if (atoms->nres > mtop->maxres_renum)
 432     {
 433         *resnr = atoms->resinfo[atoms->atom[at_loc].resind].nr;
 434     }
 435     else
 436     {
 437         /* Single residue molecule, keep counting */
 438         *resnr = maxresnr + 1 + (atnr_global - a_start)/atoms->nr*atoms->nres + atoms->atom[at_loc].resind;
 439     }
 440     *resname  = *(atoms->resinfo[atoms->atom[at_loc].resind].name);
 441 }
 442
 443 typedef struct gmx_mtop_atomloop_all
 444 {
 445     const gmx_mtop_t *mtop;
 446     int               mblock;
 447     t_atoms          *atoms;
 448     int               mol;
 449     int               maxresnr;
 450     int               at_local;
 451     int               at_global;
 452 } t_gmx_mtop_atomloop_all;
 453
 454 gmx_mtop_atomloop_all_t
 455 gmx_mtop_atomloop_all_init(const gmx_mtop_t *mtop)
 456 {
 457     struct gmx_mtop_atomloop_all *aloop;
 458
 459     snew(aloop, 1);
 460
 461     aloop->mtop         = mtop;
 462     aloop->mblock       = 0;
 463     aloop->atoms        =
 464         &mtop->moltype[mtop->molblock[aloop->mblock].type].atoms;
 465     aloop->mol          = 0;
 466     aloop->maxresnr     = mtop->maxresnr;
 467     aloop->at_local     = -1;
 468     aloop->at_global    = -1;
 469
 470     return aloop;
 471 }
 472
 473 static void gmx_mtop_atomloop_all_destroy(gmx_mtop_atomloop_all_t aloop)
 474 {
 475     sfree(aloop);
 476 }
 477
 478 gmx_bool gmx_mtop_atomloop_all_next(gmx_mtop_atomloop_all_t aloop,
 479                                     int *at_global, t_atom **atom)
 480 {
 481     if (aloop == NULL)
 482     {
 483         gmx_incons("gmx_mtop_atomloop_all_next called without calling gmx_mtop_atomloop_all_init");
 484     }
 485
 486     aloop->at_local++;
 487     aloop->at_global++;
 488
 489     if (aloop->at_local >= aloop->atoms->nr)
 490     {
 491         if (aloop->atoms->nres <= aloop->mtop->maxres_renum)
 492         {
 493             /* Single residue molecule, increase the count with one */
 494             aloop->maxresnr += aloop->atoms->nres;
 495         }
 496         aloop->mol++;
 497         aloop->at_local = 0;
 498         if (aloop->mol >= aloop->mtop->molblock[aloop->mblock].nmol)
 499         {
 500             aloop->mblock++;
 501             if (aloop->mblock >= aloop->mtop->nmolblock)
 502             {
 503                 gmx_mtop_atomloop_all_destroy(aloop);
 504                 return FALSE;
 505             }
 506             aloop->atoms = &aloop->mtop->moltype[aloop->mtop->molblock[aloop->mblock].type].atoms;
 507             aloop->mol   = 0;
 508         }
 509     }
 510
 511     *at_global = aloop->at_global;
 512     *atom      = &aloop->atoms->atom[aloop->at_local];
 513
 514     return TRUE;
 515 }
 516
 517 void gmx_mtop_atomloop_all_names(gmx_mtop_atomloop_all_t aloop,
 518                                  char **atomname, int *resnr, char **resname)
 519 {
 520     int resind_mol;
 521
 522     *atomname  = *(aloop->atoms->atomname[aloop->at_local]);
 523     resind_mol = aloop->atoms->atom[aloop->at_local].resind;
 524     *resnr     = aloop->atoms->resinfo[resind_mol].nr;
 525     if (aloop->atoms->nres <= aloop->mtop->maxres_renum)
 526     {
 527         *resnr = aloop->maxresnr + 1 + resind_mol;
 528     }
 529     *resname  = *(aloop->atoms->resinfo[resind_mol].name);
 530 }
 531
 532 void gmx_mtop_atomloop_all_moltype(gmx_mtop_atomloop_all_t aloop,
 533                                    gmx_moltype_t **moltype, int *at_mol)
 534 {
 535     *moltype = &aloop->mtop->moltype[aloop->mtop->molblock[aloop->mblock].type];
 536     *at_mol  = aloop->at_local;
 537 }
 538
 539 typedef struct gmx_mtop_atomloop_block
 540 {
 541     const gmx_mtop_t *mtop;
 542     int               mblock;
 543     t_atoms          *atoms;
 544     int               at_local;
 545 } t_gmx_mtop_atomloop_block;
 546
 547 gmx_mtop_atomloop_block_t
 548 gmx_mtop_atomloop_block_init(const gmx_mtop_t *mtop)
 549 {
 550     struct gmx_mtop_atomloop_block *aloop;
 551
 552     snew(aloop, 1);
 553
 554     aloop->mtop      = mtop;
 555     aloop->mblock    = 0;
 556     aloop->atoms     = &mtop->moltype[mtop->molblock[aloop->mblock].type].atoms;
 557     aloop->at_local  = -1;
 558
 559     return aloop;
 560 }
 561
 562 static void gmx_mtop_atomloop_block_destroy(gmx_mtop_atomloop_block_t aloop)
 563 {
 564     sfree(aloop);
 565 }
 566
 567 gmx_bool gmx_mtop_atomloop_block_next(gmx_mtop_atomloop_block_t aloop,
 568                                       t_atom **atom, int *nmol)
 569 {
 570     if (aloop == NULL)
 571     {
 572         gmx_incons("gmx_mtop_atomloop_all_next called without calling gmx_mtop_atomloop_all_init");
 573     }
 574
 575     aloop->at_local++;
 576
 577     if (aloop->at_local >= aloop->atoms->nr)
 578     {
 579         aloop->mblock++;
 580         if (aloop->mblock >= aloop->mtop->nmolblock)
 581         {
 582             gmx_mtop_atomloop_block_destroy(aloop);
 583             return FALSE;
 584         }
 585         aloop->atoms    = &aloop->mtop->moltype[aloop->mtop->molblock[aloop->mblock].type].atoms;
 586         aloop->at_local = 0;
 587     }
 588
 589     *atom = &aloop->atoms->atom[aloop->at_local];
 590     *nmol = aloop->mtop->molblock[aloop->mblock].nmol;
 591
 592     return TRUE;
 593 }
 594
 595 typedef struct gmx_mtop_ilistloop
 596 {
 597     const gmx_mtop_t *mtop;
 598     int               mblock;
 599 } t_gmx_mtop_ilist;
 600
 601 gmx_mtop_ilistloop_t
 602 gmx_mtop_ilistloop_init(const gmx_mtop_t *mtop)
 603 {
 604     struct gmx_mtop_ilistloop *iloop;
 605
 606     snew(iloop, 1);
 607
 608     iloop->mtop      = mtop;
 609     iloop->mblock    = -1;
 610
 611     return iloop;
 612 }
 613
 614 static void gmx_mtop_ilistloop_destroy(gmx_mtop_ilistloop_t iloop)
 615 {
 616     sfree(iloop);
 617 }
 618
 619 gmx_bool gmx_mtop_ilistloop_next(gmx_mtop_ilistloop_t iloop,
 620                                  t_ilist **ilist_mol, int *nmol)
 621 {
 622     if (iloop == NULL)
 623     {
 624         gmx_incons("gmx_mtop_ilistloop_next called without calling gmx_mtop_ilistloop_init");
 625     }
 626
 627     iloop->mblock++;
 628     if (iloop->mblock >= iloop->mtop->nmolblock)
 629     {
 630         if (iloop->mblock == iloop->mtop->nmolblock &&
 631             iloop->mtop->bIntermolecularInteractions)
 632         {
 633             *ilist_mol = iloop->mtop->intermolecular_ilist;
 634             *nmol      = 1;
 635             return TRUE;
 636         }
 637
 638         gmx_mtop_ilistloop_destroy(iloop);
 639         return FALSE;
 640     }
 641
 642     *ilist_mol =
 643         iloop->mtop->moltype[iloop->mtop->molblock[iloop->mblock].type].ilist;
 644
 645     *nmol = iloop->mtop->molblock[iloop->mblock].nmol;
 646
 647     return TRUE;
 648 }
 649 typedef struct gmx_mtop_ilistloop_all
 650 {
 651     const gmx_mtop_t *mtop;
 652     int               mblock;
 653     int               mol;
 654     int               a_offset;
 655 } t_gmx_mtop_ilist_all;
 656
 657 gmx_mtop_ilistloop_all_t
 658 gmx_mtop_ilistloop_all_init(const gmx_mtop_t *mtop)
 659 {
 660     struct gmx_mtop_ilistloop_all *iloop;
 661
 662     snew(iloop, 1);
 663
 664     iloop->mtop      = mtop;
 665     iloop->mblock    = 0;
 666     iloop->mol       = -1;
 667     iloop->a_offset  = 0;
 668
 669     return iloop;
 670 }
 671
 672 static void gmx_mtop_ilistloop_all_destroy(gmx_mtop_ilistloop_all_t iloop)
 673 {
 674     sfree(iloop);
 675 }
 676
 677 gmx_bool gmx_mtop_ilistloop_all_next(gmx_mtop_ilistloop_all_t iloop,
 678                                      t_ilist **ilist_mol, int *atnr_offset)
 679 {
 680     gmx_molblock_t *molb;
 681
 682     if (iloop == NULL)
 683     {
 684         gmx_incons("gmx_mtop_ilistloop_all_next called without calling gmx_mtop_ilistloop_all_init");
 685     }
 686
 687     if (iloop->mol >= 0)
 688     {
 689         iloop->a_offset += iloop->mtop->molblock[iloop->mblock].natoms_mol;
 690     }
 691
 692     iloop->mol++;
 693
 694     /* Inter-molecular interactions, if present, are indexed with
 695      * iloop->mblock == iloop->mtop->nmolblock, thus we should separately
 696      * check for this value in this conditional.
 697      */
 698     if (iloop->mblock == iloop->mtop->nmolblock ||
 699         iloop->mol >= iloop->mtop->molblock[iloop->mblock].nmol)
 700     {
 701         iloop->mblock++;
 702         iloop->mol = 0;
 703         if (iloop->mblock >= iloop->mtop->nmolblock)
 704         {
 705             if (iloop->mblock == iloop->mtop->nmolblock &&
 706                 iloop->mtop->bIntermolecularInteractions)
 707             {
 708                 *ilist_mol   = iloop->mtop->intermolecular_ilist;
 709                 *atnr_offset = 0;
 710                 return TRUE;
 711             }
 712
 713             gmx_mtop_ilistloop_all_destroy(iloop);
 714             return FALSE;
 715         }
 716     }
 717
 718     *ilist_mol =
 719         iloop->mtop->moltype[iloop->mtop->molblock[iloop->mblock].type].ilist;
 720
 721     *atnr_offset = iloop->a_offset;
 722
 723     return TRUE;
 724 }
 725
 726 int gmx_mtop_ftype_count(const gmx_mtop_t *mtop, int ftype)
 727 {
 728     gmx_mtop_ilistloop_t iloop;
 729     t_ilist             *il;
 730     int                  n, nmol;
 731
 732     n = 0;
 733
 734     iloop = gmx_mtop_ilistloop_init(mtop);
 735     while (gmx_mtop_ilistloop_next(iloop, &il, &nmol))
 736     {
 737         n += nmol*il[ftype].nr/(1+NRAL(ftype));
 738     }
 739
 740     if (mtop->bIntermolecularInteractions)
 741     {
 742         n += mtop->intermolecular_ilist[ftype].nr/(1+NRAL(ftype));
 743     }
 744
 745     return n;
 746 }
 747
 748 t_block gmx_mtop_global_cgs(const gmx_mtop_t *mtop)
 749 {
 750     t_block         cgs_gl, *cgs_mol;
 751     int             mb, mol, cg;
 752     gmx_molblock_t *molb;
 753     t_atoms        *atoms;
 754
 755     /* In most cases this is too much, but we realloc at the end */
 756     snew(cgs_gl.index, mtop->natoms+1);
 757
 758     cgs_gl.nr       = 0;
 759     cgs_gl.index[0] = 0;
 760     for (mb = 0; mb < mtop->nmolblock; mb++)
 761     {
 762         molb    = &mtop->molblock[mb];
 763         cgs_mol = &mtop->moltype[molb->type].cgs;
 764         for (mol = 0; mol < molb->nmol; mol++)
 765         {
 766             for (cg = 0; cg < cgs_mol->nr; cg++)
 767             {
 768                 cgs_gl.index[cgs_gl.nr+1] =
 769                     cgs_gl.index[cgs_gl.nr] +
 770                     cgs_mol->index[cg+1] - cgs_mol->index[cg];
 771                 cgs_gl.nr++;
 772             }
 773         }
 774     }
 775     cgs_gl.nalloc_index = cgs_gl.nr + 1;
 776     srenew(cgs_gl.index, cgs_gl.nalloc_index);
 777
 778     return cgs_gl;
 779 }
 780
 781 static void atomcat(t_atoms *dest, t_atoms *src, int copies,
 782                     int maxres_renum, int *maxresnr)
 783 {
 784     int i, j, l, size;
 785     int srcnr  = src->nr;
 786     int destnr = dest->nr;
 787
 788     if (srcnr)
 789     {
 790         size = destnr+copies*srcnr;
 791         srenew(dest->atom, size);
 792         srenew(dest->atomname, size);
 793         srenew(dest->atomtype, size);
 794         srenew(dest->atomtypeB, size);
 795     }
 796     if (src->nres)
 797     {
 798         size = dest->nres+copies*src->nres;
 799         srenew(dest->resinfo, size);
 800     }
 801
 802     /* residue information */
 803     for (l = dest->nres, j = 0; (j < copies); j++, l += src->nres)
 804     {
 805         memcpy((char *) &(dest->resinfo[l]), (char *) &(src->resinfo[0]),
 806                (size_t)(src->nres*sizeof(src->resinfo[0])));
 807     }
 808
 809     for (l = destnr, j = 0; (j < copies); j++, l += srcnr)
 810     {
 811         memcpy((char *) &(dest->atomname[l]), (char *) &(src->atomname[0]),
 812                (size_t)(srcnr*sizeof(src->atomname[0])));
 813         memcpy((char *) &(dest->atomtype[l]), (char *) &(src->atomtype[0]),
 814                (size_t)(srcnr*sizeof(src->atomtype[0])));
 815         memcpy((char *) &(dest->atomtypeB[l]), (char *) &(src->atomtypeB[0]),
 816                (size_t)(srcnr*sizeof(src->atomtypeB[0])));
 817         memcpy((char *) &(dest->atom[l]), (char *) &(src->atom[0]),
 818                (size_t)(srcnr*sizeof(src->atom[0])));
 819     }
 820
 821     /* Increment residue indices */
 822     for (l = destnr, j = 0; (j < copies); j++)
 823     {
 824         for (i = 0; (i < srcnr); i++, l++)
 825         {
 826             dest->atom[l].resind = dest->nres+j*src->nres+src->atom[i].resind;
 827         }
 828     }
 829
 830     if (src->nres <= maxres_renum)
 831     {
 832         /* Single residue molecule, continue counting residues */
 833         for (j = 0; (j < copies); j++)
 834         {
 835             for (l = 0; l < src->nres; l++)
 836             {
 837                 (*maxresnr)++;
 838                 dest->resinfo[dest->nres+j*src->nres+l].nr = *maxresnr;
 839             }
 840         }
 841     }
 842
 843     dest->nres += copies*src->nres;
 844     dest->nr   += copies*src->nr;
 845 }
 846
 847 t_atoms gmx_mtop_global_atoms(const gmx_mtop_t *mtop)
 848 {
 849     t_atoms         atoms;
 850     int             maxresnr, mb;
 851     gmx_molblock_t *molb;
 852
 853     init_t_atoms(&atoms, 0, FALSE);
 854
 855     maxresnr = mtop->maxresnr;
 856     for (mb = 0; mb < mtop->nmolblock; mb++)
 857     {
 858         molb = &mtop->molblock[mb];
 859         atomcat(&atoms, &mtop->moltype[molb->type].atoms, molb->nmol,
 860                 mtop->maxres_renum, &maxresnr);
 861     }
 862
 863     return atoms;
 864 }
 865
 866 void gmx_mtop_make_atomic_charge_groups(gmx_mtop_t *mtop,
 867                                         gmx_bool    bKeepSingleMolCG)
 868 {
 869     int      mb, cg;
 870     t_block *cgs_mol;
 871
 872     for (mb = 0; mb < mtop->nmolblock; mb++)
 873     {
 874         cgs_mol = &mtop->moltype[mtop->molblock[mb].type].cgs;
 875         if (!(bKeepSingleMolCG && cgs_mol->nr == 1))
 876         {
 877             cgs_mol->nr           = mtop->molblock[mb].natoms_mol;
 878             cgs_mol->nalloc_index = cgs_mol->nr + 1;
 879             srenew(cgs_mol->index, cgs_mol->nalloc_index);
 880             for (cg = 0; cg < cgs_mol->nr+1; cg++)
 881             {
 882                 cgs_mol->index[cg] = cg;
 883             }
 884         }
 885     }
 886 }
 887
 888 /*
 889  * The cat routines below are old code from src/kernel/topcat.c
 890  */
 891
 892 static void blockcat(t_block *dest, t_block *src, int copies)
 893 {
 894     int i, j, l, nra, size;
 895
 896     if (src->nr)
 897     {
 898         size = (dest->nr+copies*src->nr+1);
 899         srenew(dest->index, size);
 900     }
 901
 902     nra = dest->index[dest->nr];
 903     for (l = dest->nr, j = 0; (j < copies); j++)
 904     {
 905         for (i = 0; (i < src->nr); i++)
 906         {
 907             dest->index[l++] = nra + src->index[i];
 908         }
 909         nra += src->index[src->nr];
 910     }
 911     dest->nr             += copies*src->nr;
 912     dest->index[dest->nr] = nra;
 913 }
 914
 915 static void blockacat(t_blocka *dest, t_blocka *src, int copies,
 916                       int dnum, int snum)
 917 {
 918     int i, j, l, size;
 919     int destnr  = dest->nr;
 920     int destnra = dest->nra;
 921
 922     if (src->nr)
 923     {
 924         size = (dest->nr+copies*src->nr+1);
 925         srenew(dest->index, size);
 926     }
 927     if (src->nra)
 928     {
 929         size = (dest->nra+copies*src->nra);
 930         srenew(dest->a, size);
 931     }
 932
 933     for (l = destnr, j = 0; (j < copies); j++)
 934     {
 935         for (i = 0; (i < src->nr); i++)
 936         {
 937             dest->index[l++] = dest->nra+src->index[i];
 938         }
 939         dest->nra += src->nra;
 940     }
 941     for (l = destnra, j = 0; (j < copies); j++)
 942     {
 943         for (i = 0; (i < src->nra); i++)
 944         {
 945             dest->a[l++] = dnum+src->a[i];
 946         }
 947         dnum     += snum;
 948         dest->nr += src->nr;
 949     }
 950     dest->index[dest->nr] = dest->nra;
 951 }
 952
 953 static void ilistcat(int ftype, t_ilist *dest, t_ilist *src, int copies,
 954                      int dnum, int snum)
 955 {
 956     int nral, c, i, a;
 957
 958     nral = NRAL(ftype);
 959
 960     dest->nalloc = dest->nr + copies*src->nr;
 961     srenew(dest->iatoms, dest->nalloc);
 962
 963     for (c = 0; c < copies; c++)
 964     {
 965         for (i = 0; i < src->nr; )
 966         {
 967             dest->iatoms[dest->nr++] = src->iatoms[i++];
 968             for (a = 0; a < nral; a++)
 969             {
 970                 dest->iatoms[dest->nr++] = dnum + src->iatoms[i++];
 971             }
 972         }
 973         dnum += snum;
 974     }
 975 }
 976
 977 static void set_posres_params(t_idef *idef, gmx_molblock_t *molb,
 978                               int i0, int a_offset)
 979 {
 980     t_ilist   *il;
 981     int        i1, i, a_molb;
 982     t_iparams *ip;
 983
 984     il = &idef->il[F_POSRES];
 985     i1 = il->nr/2;
 986     idef->iparams_posres_nalloc = i1;
 987     srenew(idef->iparams_posres, idef->iparams_posres_nalloc);
 988     for (i = i0; i < i1; i++)
 989     {
 990         ip = &idef->iparams_posres[i];
 991         /* Copy the force constants */
 992         *ip    = idef->iparams[il->iatoms[i*2]];
 993         a_molb = il->iatoms[i*2+1] - a_offset;
 994         if (molb->nposres_xA == 0)
 995         {
 996             gmx_incons("Position restraint coordinates are missing");
 997         }
 998         ip->posres.pos0A[XX] = molb->posres_xA[a_molb][XX];
 999         ip->posres.pos0A[YY] = molb->posres_xA[a_molb][YY];
1000         ip->posres.pos0A[ZZ] = molb->posres_xA[a_molb][ZZ];
1001         if (molb->nposres_xB > 0)
1002         {
1003             ip->posres.pos0B[XX] = molb->posres_xB[a_molb][XX];
1004             ip->posres.pos0B[YY] = molb->posres_xB[a_molb][YY];
1005             ip->posres.pos0B[ZZ] = molb->posres_xB[a_molb][ZZ];
1006         }
1007         else
1008         {
1009             ip->posres.pos0B[XX] = ip->posres.pos0A[XX];
1010             ip->posres.pos0B[YY] = ip->posres.pos0A[YY];
1011             ip->posres.pos0B[ZZ] = ip->posres.pos0A[ZZ];
1012         }
1013         /* Set the parameter index for idef->iparams_posre */
1014         il->iatoms[i*2] = i;
1015     }
1016 }
1017
1018 static void set_fbposres_params(t_idef *idef, gmx_molblock_t *molb,
1019                                 int i0, int a_offset)
1020 {
1021     t_ilist   *il;
1022     int        i1, i, a_molb;
1023     t_iparams *ip;
1024
1025     il = &idef->il[F_FBPOSRES];
1026     i1 = il->nr/2;
1027     idef->iparams_fbposres_nalloc = i1;
1028     srenew(idef->iparams_fbposres, idef->iparams_fbposres_nalloc);
1029     for (i = i0; i < i1; i++)
1030     {
1031         ip = &idef->iparams_fbposres[i];
1032         /* Copy the force constants */
1033         *ip    = idef->iparams[il->iatoms[i*2]];
1034         a_molb = il->iatoms[i*2+1] - a_offset;
1035         if (molb->nposres_xA == 0)
1036         {
1037             gmx_incons("Position restraint coordinates are missing");
1038         }
1039         /* Take flat-bottom posres reference from normal position restraints */
1040         ip->fbposres.pos0[XX] = molb->posres_xA[a_molb][XX];
1041         ip->fbposres.pos0[YY] = molb->posres_xA[a_molb][YY];
1042         ip->fbposres.pos0[ZZ] = molb->posres_xA[a_molb][ZZ];
1043         /* Note: no B-type for flat-bottom posres */
1044
1045         /* Set the parameter index for idef->iparams_posre */
1046         il->iatoms[i*2] = i;
1047     }
1048 }
1049
1050 static void gen_local_top(const gmx_mtop_t *mtop, const t_inputrec *ir,
1051                           gmx_bool bMergeConstr,
1052                           gmx_localtop_t *top)
1053 {
1054     int                     mb, srcnr, destnr, ftype, ftype_dest, mt, natoms, mol, nposre_old, nfbposre_old;
1055     gmx_molblock_t         *molb;
1056     gmx_moltype_t          *molt;
1057     const gmx_ffparams_t   *ffp;
1058     t_idef                 *idef;
1059     real                   *qA, *qB;
1060     gmx_mtop_atomloop_all_t aloop;
1061     int                     ag;
1062     t_atom                 *atom;
1063
1064     top->atomtypes = mtop->atomtypes;
1065
1066     ffp = &mtop->ffparams;
1067
1068     idef                          = &top->idef;
1069     idef->ntypes                  = ffp->ntypes;
1070     idef->atnr                    = ffp->atnr;
1071     idef->functype                = ffp->functype;
1072     idef->iparams                 = ffp->iparams;
1073     idef->iparams_posres          = NULL;
1074     idef->iparams_posres_nalloc   = 0;
1075     idef->iparams_fbposres        = NULL;
1076     idef->iparams_fbposres_nalloc = 0;
1077     idef->fudgeQQ                 = ffp->fudgeQQ;
1078     idef->cmap_grid               = ffp->cmap_grid;
1079     idef->ilsort                  = ilsortUNKNOWN;
1080
1081     init_block(&top->cgs);
1082     init_blocka(&top->excls);
1083     for (ftype = 0; ftype < F_NRE; ftype++)
1084     {
1085         idef->il[ftype].nr     = 0;
1086         idef->il[ftype].nalloc = 0;
1087         idef->il[ftype].iatoms = NULL;
1088     }
1089
1090     natoms = 0;
1091     for (mb = 0; mb < mtop->nmolblock; mb++)
1092     {
1093         molb = &mtop->molblock[mb];
1094         molt = &mtop->moltype[molb->type];
1095
1096         srcnr  = molt->atoms.nr;
1097         destnr = natoms;
1098
1099         blockcat(&top->cgs, &molt->cgs, molb->nmol);
1100
1101         blockacat(&top->excls, &molt->excls, molb->nmol, destnr, srcnr);
1102
1103         nposre_old   = idef->il[F_POSRES].nr;
1104         nfbposre_old = idef->il[F_FBPOSRES].nr;
1105         for (ftype = 0; ftype < F_NRE; ftype++)
1106         {
1107             if (bMergeConstr &&
1108                 ftype == F_CONSTR && molt->ilist[F_CONSTRNC].nr > 0)
1109             {
1110                 /* Merge all constrains into one ilist.
1111                  * This simplifies the constraint code.
1112                  */
1113                 for (mol = 0; mol < molb->nmol; mol++)
1114                 {
1115                     ilistcat(ftype, &idef->il[F_CONSTR], &molt->ilist[F_CONSTR],
1116                              1, destnr+mol*srcnr, srcnr);
1117                     ilistcat(ftype, &idef->il[F_CONSTR], &molt->ilist[F_CONSTRNC],
1118                              1, destnr+mol*srcnr, srcnr);
1119                 }
1120             }
1121             else if (!(bMergeConstr && ftype == F_CONSTRNC))
1122             {
1123                 ilistcat(ftype, &idef->il[ftype], &molt->ilist[ftype],
1124                          molb->nmol, destnr, srcnr);
1125             }
1126         }
1127         if (idef->il[F_POSRES].nr > nposre_old)
1128         {
1129             /* Executing this line line stops gmxdump -sys working
1130              * correctly. I'm not aware there's an elegant fix. */
1131             set_posres_params(idef, molb, nposre_old/2, natoms);
1132         }
1133         if (idef->il[F_FBPOSRES].nr > nfbposre_old)
1134         {
1135             set_fbposres_params(idef, molb, nfbposre_old/2, natoms);
1136         }
1137
1138         natoms += molb->nmol*srcnr;
1139     }
1140
1141     if (mtop->bIntermolecularInteractions)
1142     {
1143         for (ftype = 0; ftype < F_NRE; ftype++)
1144         {
1145             ilistcat(ftype, &idef->il[ftype], &mtop->intermolecular_ilist[ftype],
1146                      1, 0, mtop->natoms);
1147         }
1148     }
1149
1150     if (ir == NULL)
1151     {
1152         top->idef.ilsort = ilsortUNKNOWN;
1153     }
1154     else
1155     {
1156         if (ir->efep != efepNO && gmx_mtop_bondeds_free_energy(mtop))
1157         {
1158             snew(qA, mtop->natoms);
1159             snew(qB, mtop->natoms);
1160             aloop = gmx_mtop_atomloop_all_init(mtop);
1161             while (gmx_mtop_atomloop_all_next(aloop, &ag, &atom))
1162             {
1163                 qA[ag] = atom->q;
1164                 qB[ag] = atom->qB;
1165             }
1166             gmx_sort_ilist_fe(&top->idef, qA, qB);
1167             sfree(qA);
1168             sfree(qB);
1169         }
1170         else
1171         {
1172             top->idef.ilsort = ilsortNO_FE;
1173         }
1174     }
1175 }
1176
1177 gmx_localtop_t *gmx_mtop_generate_local_top(const gmx_mtop_t *mtop,
1178                                             const t_inputrec *ir)
1179 {
1180     gmx_localtop_t *top;
1181
1182     snew(top, 1);
1183
1184     gen_local_top(mtop, ir, TRUE, top);
1185
1186     return top;
1187 }
1188
1189 t_topology gmx_mtop_t_to_t_topology(gmx_mtop_t *mtop)
1190 {
1191     int            mt, mb;
1192     gmx_localtop_t ltop;
1193     t_topology     top;
1194
1195     gen_local_top(mtop, NULL, FALSE, &ltop);
1196
1197     top.name                        = mtop->name;
1198     top.idef                        = ltop.idef;
1199     top.atomtypes                   = ltop.atomtypes;
1200     top.cgs                         = ltop.cgs;
1201     top.excls                       = ltop.excls;
1202     top.atoms                       = gmx_mtop_global_atoms(mtop);
1203     top.mols                        = mtop->mols;
1204     top.bIntermolecularInteractions = mtop->bIntermolecularInteractions;
1205     top.symtab                      = mtop->symtab;
1206
1207     /* We only need to free the moltype and molblock data,
1208      * all other pointers have been copied to top.
1209      *
1210      * Well, except for the group data, but we can't free those, because they
1211      * are used somewhere even after a call to this function.
1212      */
1213     for (mt = 0; mt < mtop->nmoltype; mt++)
1214     {
1215         done_moltype(&mtop->moltype[mt]);
1216     }
1217     sfree(mtop->moltype);
1218
1219     for (mb = 0; mb < mtop->nmolblock; mb++)
1220     {
1221         done_molblock(&mtop->molblock[mb]);
1222     }
1223     sfree(mtop->molblock);
1224
1225     return top;
1226 }