src/gromacs/topology/mtop_util.c

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