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