src/gromacs/gmxlib/mtop_util.c

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