src/gromacs/gmxlib/mtop_util.c

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