src/gmxlib/mtop_util.c

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