src/gmxlib/mtop_util.c

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