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