src/gromacs/topology/mtop_util.c

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