src/gromacs/topology/index.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
   5  * Copyright (c) 2001-2004, The GROMACS development team.
   6  * Copyright (c) 2013,2014, by the GROMACS development team, led by
   7  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   8  * and including many others, as listed in the AUTHORS file in the
   9  * top-level source directory and at http://www.gromacs.org.
  10  *
  11  * GROMACS is free software; you can redistribute it and/or
  12  * modify it under the terms of the GNU Lesser General Public License
  13  * as published by the Free Software Foundation; either version 2.1
  14  * of the License, or (at your option) any later version.
  15  *
  16  * GROMACS is distributed in the hope that it will be useful,
  17  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  18  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19  * Lesser General Public License for more details.
  20  *
  21  * You should have received a copy of the GNU Lesser General Public
  22  * License along with GROMACS; if not, see
  23  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  24  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  25  *
  26  * If you want to redistribute modifications to GROMACS, please
  27  * consider that scientific software is very special. Version
  28  * control is crucial - bugs must be traceable. We will be happy to
  29  * consider code for inclusion in the official distribution, but
  30  * derived work must not be called official GROMACS. Details are found
  31  * in the README & COPYING files - if they are missing, get the
  32  * official version at http://www.gromacs.org.
  33  *
  34  * To help us fund GROMACS development, we humbly ask that you cite
  35  * the research papers on the package. Check out http://www.gromacs.org.
  36  */
  37 #include "gmxpre.h"
  38
  39 #include "gromacs/topology/index.h"
  40
  41 #include <assert.h>
  42 #include <ctype.h>
  43 #include <stdlib.h>
  44 #include <string.h>
  45
  46 #include <algorithm>
  47
  48 #include "gromacs/legacyheaders/macros.h"
  49 #include "gromacs/legacyheaders/txtdump.h"
  50
  51 #include "gromacs/fileio/gmxfio.h"
  52 #include "gromacs/fileio/strdb.h"
  53 #include "gromacs/topology/atoms.h"
  54 #include "gromacs/topology/block.h"
  55 #include "gromacs/topology/invblock.h"
  56 #include "gromacs/topology/residuetypes.h"
  57 #include "gromacs/utility/cstringutil.h"
  58 #include "gromacs/utility/fatalerror.h"
  59 #include "gromacs/utility/smalloc.h"
  60
  61 static gmx_bool gmx_ask_yesno(gmx_bool bASK)
  62 {
  63     char c;
  64
  65     if (bASK)
  66     {
  67         do
  68         {
  69             c = toupper(fgetc(stdin));
  70         }
  71         while ((c != 'Y') && (c != 'N'));
  72
  73         return (c == 'Y');
  74     }
  75     else
  76     {
  77         return FALSE;
  78     }
  79 }
  80
  81 void write_index(const char *outf, t_blocka *b, char **gnames, gmx_bool bDuplicate, int natoms)
  82 {
  83     FILE *out;
  84     int   i, j, k;
  85
  86     out = gmx_fio_fopen(outf, "w");
  87     /* fprintf(out,"%5d  %5d\n",b->nr,b->nra); */
  88     for (i = 0; (i < b->nr); i++)
  89     {
  90         fprintf(out, "[ %s ]\n", gnames[i]);
  91         for (k = 0, j = b->index[i]; j < b->index[i+1]; j++, k++)
  92         {
  93             fprintf(out, "%4d ", b->a[j]+1);
  94             if ((k % 15) == 14)
  95             {
  96                 fprintf(out, "\n");
  97             }
  98         }
  99         fprintf(out, "\n");
 100     }
 101
 102     /* Duplicate copy, useful for computational electrophysiology double-layer setups */
 103     if (bDuplicate)
 104     {
 105         fprintf(stderr, "Duplicating the whole system with an atom offset of %d atoms.\n", natoms);
 106         for (i = 0; (i < b->nr); i++)
 107         {
 108             fprintf(out, "[ %s_copy ]\n", gnames[i]);
 109             for (k = 0, j = b->index[i]; j < b->index[i+1]; j++, k++)
 110             {
 111                 fprintf(out, "%4d ", b->a[j]+1 + natoms );
 112                 if ((k % 15) == 14)
 113                 {
 114                     fprintf(out, "\n");
 115                 }
 116             }
 117             fprintf(out, "\n");
 118         }
 119     }
 120
 121     gmx_fio_fclose(out);
 122 }
 123
 124 void add_grp(t_blocka *b, char ***gnames, int nra, atom_id a[], const char *name)
 125 {
 126     int i;
 127
 128     srenew(b->index, b->nr+2);
 129     srenew(*gnames, b->nr+1);
 130     (*gnames)[b->nr] = gmx_strdup(name);
 131
 132     srenew(b->a, b->nra+nra);
 133     for (i = 0; (i < nra); i++)
 134     {
 135         b->a[b->nra++] = a[i];
 136     }
 137     b->nr++;
 138     b->index[b->nr] = b->nra;
 139 }
 140
 141 /* compare index in `a' with group in `b' at `index',
 142    when `index'<0 it is relative to end of `b' */
 143 static gmx_bool grp_cmp(t_blocka *b, int nra, atom_id a[], int index)
 144 {
 145     int i;
 146
 147     if (index < 0)
 148     {
 149         index = b->nr-1+index;
 150     }
 151     if (index >= b->nr)
 152     {
 153         gmx_fatal(FARGS, "no such index group %d in t_blocka (nr=%d)", index, b->nr);
 154     }
 155     /* compare sizes */
 156     if (nra != b->index[index+1] - b->index[index])
 157     {
 158         return FALSE;
 159     }
 160     for (i = 0; i < nra; i++)
 161     {
 162         if (a[i] != b->a[b->index[index]+i])
 163         {
 164             return FALSE;
 165         }
 166     }
 167     return TRUE;
 168 }
 169
 170 static void
 171 p_status(const char **restype, int nres, const char **typenames, int ntypes)
 172 {
 173     int   i, j;
 174     int * counter;
 175
 176     snew(counter, ntypes);
 177     for (i = 0; i < ntypes; i++)
 178     {
 179         counter[i] = 0;
 180     }
 181     for (i = 0; i < nres; i++)
 182     {
 183         for (j = 0; j < ntypes; j++)
 184         {
 185             if (!gmx_strcasecmp(restype[i], typenames[j]))
 186             {
 187                 counter[j]++;
 188             }
 189         }
 190     }
 191
 192     for (i = 0; (i < ntypes); i++)
 193     {
 194         if (counter[i] > 0)
 195         {
 196             printf("There are: %5d %10s residues\n", counter[i], typenames[i]);
 197         }
 198     }
 199
 200     sfree(counter);
 201 }
 202
 203
 204 static atom_id *
 205 mk_aid(t_atoms *atoms, const char ** restype, const char * typestring, int *nra, gmx_bool bMatch)
 206 /* Make an array of atom_ids for all atoms with residuetypes matching typestring, or the opposite if bMatch is false */
 207 {
 208     atom_id *a;
 209     int      i;
 210     int      res;
 211
 212     snew(a, atoms->nr);
 213     *nra = 0;
 214     for (i = 0; (i < atoms->nr); i++)
 215     {
 216         res = !gmx_strcasecmp(restype[atoms->atom[i].resind], typestring);
 217         if (bMatch == FALSE)
 218         {
 219             res = !res;
 220         }
 221         if (res)
 222         {
 223             a[(*nra)++] = i;
 224         }
 225     }
 226
 227     return a;
 228 }
 229
 230 typedef struct {
 231     char    *rname;
 232     gmx_bool bNeg;
 233     char    *gname;
 234 } restp_t;
 235
 236 static void analyse_other(const char ** restype, t_atoms *atoms,
 237                           t_blocka *gb, char ***gn, gmx_bool bASK, gmx_bool bVerb)
 238 {
 239     restp_t *restp = NULL;
 240     char   **attp  = NULL;
 241     char    *rname, *aname;
 242     atom_id *other_ndx, *aid, *aaid;
 243     int      i, j, k, l, resind, naid, naaid, natp, nrestp = 0;
 244
 245     for (i = 0; (i < atoms->nres); i++)
 246     {
 247         if (gmx_strcasecmp(restype[i], "Protein") && gmx_strcasecmp(restype[i], "DNA") && gmx_strcasecmp(restype[i], "RNA") && gmx_strcasecmp(restype[i], "Water"))
 248         {
 249             break;
 250         }
 251     }
 252     if (i < atoms->nres)
 253     {
 254         /* we have others */
 255         if (bVerb)
 256         {
 257             printf("Analysing residues not classified as Protein/DNA/RNA/Water and splitting into groups...\n");
 258         }
 259         snew(other_ndx, atoms->nr);
 260         for (k = 0; (k < atoms->nr); k++)
 261         {
 262             resind = atoms->atom[k].resind;
 263             rname  = *atoms->resinfo[resind].name;
 264             if (gmx_strcasecmp(restype[resind], "Protein") && gmx_strcasecmp(restype[resind], "DNA") &&
 265                 gmx_strcasecmp(restype[resind], "RNA") && gmx_strcasecmp(restype[resind], "Water"))
 266             {
 267
 268                 for (l = 0; (l < nrestp); l++)
 269                 {
 270                     if (strcmp(restp[l].rname, rname) == 0)
 271                     {
 272                         break;
 273                     }
 274                 }
 275                 if (l == nrestp)
 276                 {
 277                     srenew(restp, nrestp+1);
 278                     restp[nrestp].rname = gmx_strdup(rname);
 279                     restp[nrestp].bNeg  = FALSE;
 280                     restp[nrestp].gname = gmx_strdup(rname);
 281                     nrestp++;
 282
 283                 }
 284             }
 285         }
 286         for (i = 0; (i < nrestp); i++)
 287         {
 288             snew(aid, atoms->nr);
 289             naid = 0;
 290             for (j = 0; (j < atoms->nr); j++)
 291             {
 292                 rname = *atoms->resinfo[atoms->atom[j].resind].name;
 293                 if ((strcmp(restp[i].rname, rname) == 0 && !restp[i].bNeg) ||
 294                     (strcmp(restp[i].rname, rname) != 0 &&  restp[i].bNeg))
 295                 {
 296                     aid[naid++] = j;
 297                 }
 298             }
 299             add_grp(gb, gn, naid, aid, restp[i].gname);
 300             if (bASK)
 301             {
 302                 printf("split %s into atoms (y/n) ? ", restp[i].gname);
 303                 fflush(stdout);
 304                 if (gmx_ask_yesno(bASK))
 305                 {
 306                     natp = 0;
 307                     for (k = 0; (k < naid); k++)
 308                     {
 309                         aname = *atoms->atomname[aid[k]];
 310                         for (l = 0; (l < natp); l++)
 311                         {
 312                             if (strcmp(aname, attp[l]) == 0)
 313                             {
 314                                 break;
 315                             }
 316                         }
 317                         if (l == natp)
 318                         {
 319                             srenew(attp, ++natp);
 320                             attp[natp-1] = aname;
 321                         }
 322                     }
 323                     if (natp > 1)
 324                     {
 325                         for (l = 0; (l < natp); l++)
 326                         {
 327                             snew(aaid, naid);
 328                             naaid = 0;
 329                             for (k = 0; (k < naid); k++)
 330                             {
 331                                 aname = *atoms->atomname[aid[k]];
 332                                 if (strcmp(aname, attp[l]) == 0)
 333                                 {
 334                                     aaid[naaid++] = aid[k];
 335                                 }
 336                             }
 337                             add_grp(gb, gn, naaid, aaid, attp[l]);
 338                             sfree(aaid);
 339                         }
 340                     }
 341                     sfree(attp);
 342                     attp = NULL;
 343                 }
 344                 sfree(aid);
 345             }
 346         }
 347         sfree(other_ndx);
 348     }
 349 }
 350
 351 /*! \brief
 352  * Cata necessary to construct a single (protein) index group in
 353  * analyse_prot().
 354  */
 355 typedef struct gmx_help_make_index_group
 356 {
 357     /** The set of atom names that will be used to form this index group */
 358     const char **defining_atomnames;
 359     /** Size of the defining_atomnames array */
 360     int          num_defining_atomnames;
 361     /** Name of this index group */
 362     const char  *group_name;
 363     /** Whether the above atom names name the atoms in the group, or
 364         those not in the group */
 365     gmx_bool bTakeComplement;
 366     /** The index in wholename gives the first item in the arrays of
 367        atomnames that should be tested with 'gmx_strncasecmp' in stead of
 368        gmx_strcasecmp, or -1 if all items should be tested with strcasecmp
 369        This is comparable to using a '*' wildcard at the end of specific
 370        atom names, but that is more involved to implement...
 371      */
 372     int wholename;
 373     /** Only create this index group if it differs from the one specified in compareto,
 374        where -1 means to always create this group. */
 375     int compareto;
 376 } t_gmx_help_make_index_group;
 377
 378 static void analyse_prot(const char ** restype, t_atoms *atoms,
 379                          t_blocka *gb, char ***gn, gmx_bool bASK, gmx_bool bVerb)
 380 {
 381     /* lists of atomnames to be used in constructing index groups: */
 382     static const char *pnoh[]    = { "H", "HN" };
 383     static const char *pnodum[]  = {
 384         "MN1",  "MN2",  "MCB1", "MCB2", "MCG1", "MCG2",
 385         "MCD1", "MCD2", "MCE1", "MCE2", "MNZ1", "MNZ2"
 386     };
 387     static const char *calpha[]  = { "CA" };
 388     static const char *bb[]      = { "N", "CA", "C" };
 389     static const char *mc[]      = { "N", "CA", "C", "O", "O1", "O2", "OC1", "OC2", "OT", "OXT" };
 390     static const char *mcb[]     = { "N", "CA", "CB", "C", "O", "O1", "O2", "OC1", "OC2", "OT", "OXT" };
 391     static const char *mch[]     = {
 392         "N", "CA", "C", "O", "O1", "O2", "OC1", "OC2", "OT", "OXT",
 393         "H1", "H2", "H3", "H", "HN"
 394     };
 395
 396     static const t_gmx_help_make_index_group constructing_data[] =
 397     {{ NULL,   0, "Protein",      TRUE,  -1, -1},
 398      { pnoh,   asize(pnoh),   "Protein-H",    TRUE,  0,  -1},
 399      { calpha, asize(calpha), "C-alpha",      FALSE, -1, -1},
 400      { bb,     asize(bb),     "Backbone",     FALSE, -1, -1},
 401      { mc,     asize(mc),     "MainChain",    FALSE, -1, -1},
 402      { mcb,    asize(mcb),    "MainChain+Cb", FALSE, -1, -1},
 403      { mch,    asize(mch),    "MainChain+H",  FALSE, -1, -1},
 404      { mch,    asize(mch),    "SideChain",    TRUE,  -1, -1},
 405      { mch,    asize(mch),    "SideChain-H",  TRUE,  11, -1},
 406      { pnodum, asize(pnodum), "Prot-Masses",  TRUE,  -1, 0}, };
 407     const int   num_index_groups = asize(constructing_data);
 408
 409     int         n, j;
 410     atom_id    *aid;
 411     int         nra, npres;
 412     gmx_bool    match;
 413     char        ndx_name[STRLEN], *atnm;
 414     int         i;
 415
 416     if (bVerb)
 417     {
 418         printf("Analysing Protein...\n");
 419     }
 420     snew(aid, atoms->nr);
 421
 422     /* calculate the number of protein residues */
 423     npres = 0;
 424     for (i = 0; (i < atoms->nres); i++)
 425     {
 426         if (0 == gmx_strcasecmp(restype[i], "Protein"))
 427         {
 428             npres++;
 429         }
 430     }
 431     /* find matching or complement atoms */
 432     for (i = 0; (i < (int)num_index_groups); i++)
 433     {
 434         nra = 0;
 435         for (n = 0; (n < atoms->nr); n++)
 436         {
 437             if (0 == gmx_strcasecmp(restype[atoms->atom[n].resind], "Protein"))
 438             {
 439                 match = FALSE;
 440                 for (j = 0; (j < constructing_data[i].num_defining_atomnames); j++)
 441                 {
 442                     /* skip digits at beginning of atomname, e.g. 1H */
 443                     atnm = *atoms->atomname[n];
 444                     while (isdigit(atnm[0]))
 445                     {
 446                         atnm++;
 447                     }
 448                     if ( (constructing_data[i].wholename == -1) || (j < constructing_data[i].wholename) )
 449                     {
 450                         if (0 == gmx_strcasecmp(constructing_data[i].defining_atomnames[j], atnm))
 451                         {
 452                             match = TRUE;
 453                         }
 454                     }
 455                     else
 456                     {
 457                         if (0 == gmx_strncasecmp(constructing_data[i].defining_atomnames[j], atnm, strlen(constructing_data[i].defining_atomnames[j])))
 458                         {
 459                             match = TRUE;
 460                         }
 461                     }
 462                 }
 463                 if (constructing_data[i].bTakeComplement != match)
 464                 {
 465                     aid[nra++] = n;
 466                 }
 467             }
 468         }
 469         /* if we want to add this group always or it differs from previous
 470            group, add it: */
 471         if (-1 == constructing_data[i].compareto || !grp_cmp(gb, nra, aid, constructing_data[i].compareto-i) )
 472         {
 473             add_grp(gb, gn, nra, aid, constructing_data[i].group_name);
 474         }
 475     }
 476
 477     if (bASK)
 478     {
 479         for (i = 0; (i < (int)num_index_groups); i++)
 480         {
 481             printf("Split %12s into %5d residues (y/n) ? ", constructing_data[i].group_name, npres);
 482             if (gmx_ask_yesno(bASK))
 483             {
 484                 int resind;
 485                 nra = 0;
 486                 for (n = 0; ((atoms->atom[n].resind < npres) && (n < atoms->nr)); )
 487                 {
 488                     resind = atoms->atom[n].resind;
 489                     for (; ((atoms->atom[n].resind == resind) && (n < atoms->nr)); n++)
 490                     {
 491                         match = FALSE;
 492                         for (j = 0; (j < constructing_data[i].num_defining_atomnames); j++)
 493                         {
 494                             if (0 == gmx_strcasecmp(constructing_data[i].defining_atomnames[j], *atoms->atomname[n]))
 495                             {
 496                                 match = TRUE;
 497                             }
 498                         }
 499                         if (constructing_data[i].bTakeComplement != match)
 500                         {
 501                             aid[nra++] = n;
 502                         }
 503                     }
 504                     /* copy the residuename to the tail of the groupname */
 505                     if (nra > 0)
 506                     {
 507                         t_resinfo *ri;
 508                         ri = &atoms->resinfo[resind];
 509                         sprintf(ndx_name, "%s_%s%d%c",
 510                                 constructing_data[i].group_name, *ri->name, ri->nr, ri->ic == ' ' ? '\0' : ri->ic);
 511                         add_grp(gb, gn, nra, aid, ndx_name);
 512                         nra = 0;
 513                     }
 514                 }
 515             }
 516         }
 517         printf("Make group with sidechain and C=O swapped (y/n) ? ");
 518         if (gmx_ask_yesno(bASK))
 519         {
 520             /* Make swap sidechain C=O index */
 521             int resind, hold;
 522             nra = 0;
 523             for (n = 0; ((atoms->atom[n].resind < npres) && (n < atoms->nr)); )
 524             {
 525                 resind = atoms->atom[n].resind;
 526                 hold   = -1;
 527                 for (; ((atoms->atom[n].resind == resind) && (n < atoms->nr)); n++)
 528                 {
 529                     if (strcmp("CA", *atoms->atomname[n]) == 0)
 530                     {
 531                         aid[nra++] = n;
 532                         hold       = nra;
 533                         nra       += 2;
 534                     }
 535                     else if (strcmp("C", *atoms->atomname[n]) == 0)
 536                     {
 537                         if (hold == -1)
 538                         {
 539                             gmx_incons("Atom naming problem");
 540                         }
 541                         aid[hold] = n;
 542                     }
 543                     else if (strcmp("O", *atoms->atomname[n]) == 0)
 544                     {
 545                         if (hold == -1)
 546                         {
 547                             gmx_incons("Atom naming problem");
 548                         }
 549                         aid[hold+1] = n;
 550                     }
 551                     else if (strcmp("O1", *atoms->atomname[n]) == 0)
 552                     {
 553                         if (hold == -1)
 554                         {
 555                             gmx_incons("Atom naming problem");
 556                         }
 557                         aid[hold+1] = n;
 558                     }
 559                     else
 560                     {
 561                         aid[nra++] = n;
 562                     }
 563                 }
 564             }
 565             /* copy the residuename to the tail of the groupname */
 566             if (nra > 0)
 567             {
 568                 add_grp(gb, gn, nra, aid, "SwapSC-CO");
 569             }
 570         }
 571     }
 572     sfree(aid);
 573 }
 574
 575
 576 void analyse(t_atoms *atoms, t_blocka *gb, char ***gn, gmx_bool bASK, gmx_bool bVerb)
 577 {
 578     gmx_residuetype_t*rt = NULL;
 579     char             *resnm;
 580     atom_id          *aid;
 581     const char    **  restype;
 582     int               nra;
 583     int               i, k;
 584     int               ntypes;
 585     const char     ** p_typename;
 586     int               iwater, iion;
 587     int               nwater, nion;
 588     int               found;
 589
 590     if (bVerb)
 591     {
 592         printf("Analysing residue names:\n");
 593     }
 594     /* Create system group, every single atom */
 595     snew(aid, atoms->nr);
 596     for (i = 0; i < atoms->nr; i++)
 597     {
 598         aid[i] = i;
 599     }
 600     add_grp(gb, gn, atoms->nr, aid, "System");
 601     sfree(aid);
 602
 603     /* For every residue, get a pointer to the residue type name */
 604     gmx_residuetype_init(&rt);
 605     assert(rt);
 606
 607     snew(restype, atoms->nres);
 608     ntypes     = 0;
 609     p_typename = NULL;
 610     for (i = 0; i < atoms->nres; i++)
 611     {
 612         resnm = *atoms->resinfo[i].name;
 613         gmx_residuetype_get_type(rt, resnm, &(restype[i]));
 614
 615         /* Note that this does not lead to a N*N loop, but N*K, where
 616          * K is the number of residue _types_, which is small and independent of N.
 617          */
 618         found = 0;
 619         for (k = 0; k < ntypes && !found; k++)
 620         {
 621             assert(p_typename != NULL);
 622             found = !strcmp(restype[i], p_typename[k]);
 623         }
 624         if (!found)
 625         {
 626             srenew(p_typename, ntypes+1);
 627             p_typename[ntypes++] = gmx_strdup(restype[i]);
 628         }
 629     }
 630
 631     if (bVerb)
 632     {
 633         p_status(restype, atoms->nres, p_typename, ntypes);
 634     }
 635
 636     for (k = 0; k < ntypes; k++)
 637     {
 638         aid = mk_aid(atoms, restype, p_typename[k], &nra, TRUE);
 639
 640         /* Check for special types to do fancy stuff with */
 641
 642         if (!gmx_strcasecmp(p_typename[k], "Protein") && nra > 0)
 643         {
 644             sfree(aid);
 645             /* PROTEIN */
 646             analyse_prot(restype, atoms, gb, gn, bASK, bVerb);
 647
 648             /* Create a Non-Protein group */
 649             aid = mk_aid(atoms, restype, "Protein", &nra, FALSE);
 650             if ((nra > 0) && (nra < atoms->nr))
 651             {
 652                 add_grp(gb, gn, nra, aid, "non-Protein");
 653             }
 654             sfree(aid);
 655         }
 656         else if (!gmx_strcasecmp(p_typename[k], "Water") && nra > 0)
 657         {
 658             add_grp(gb, gn, nra, aid, p_typename[k]);
 659             /* Add this group as 'SOL' too, for backward compatibility with older gromacs versions */
 660             add_grp(gb, gn, nra, aid, "SOL");
 661
 662             sfree(aid);
 663
 664             /* Solvent, create a negated group too */
 665             aid = mk_aid(atoms, restype, "Water", &nra, FALSE);
 666             if ((nra > 0) && (nra < atoms->nr))
 667             {
 668                 add_grp(gb, gn, nra, aid, "non-Water");
 669             }
 670             sfree(aid);
 671         }
 672         else if (nra > 0)
 673         {
 674             /* Other groups */
 675             add_grp(gb, gn, nra, aid, p_typename[k]);
 676             sfree(aid);
 677             analyse_other(restype, atoms, gb, gn, bASK, bVerb);
 678         }
 679     }
 680
 681     sfree(p_typename);
 682     sfree(restype);
 683     gmx_residuetype_destroy(rt);
 684
 685     /* Create a merged water_and_ions group */
 686     iwater = -1;
 687     iion   = -1;
 688     nwater = 0;
 689     nion   = 0;
 690
 691     for (i = 0; i < gb->nr; i++)
 692     {
 693         if (!gmx_strcasecmp((*gn)[i], "Water"))
 694         {
 695             iwater = i;
 696             nwater = gb->index[i+1]-gb->index[i];
 697         }
 698         else if (!gmx_strcasecmp((*gn)[i], "Ion"))
 699         {
 700             iion = i;
 701             nion = gb->index[i+1]-gb->index[i];
 702         }
 703     }
 704
 705     if (nwater > 0 && nion > 0)
 706     {
 707         srenew(gb->index, gb->nr+2);
 708         srenew(*gn, gb->nr+1);
 709         (*gn)[gb->nr] = gmx_strdup("Water_and_ions");
 710         srenew(gb->a, gb->nra+nwater+nion);
 711         if (nwater > 0)
 712         {
 713             for (i = gb->index[iwater]; i < gb->index[iwater+1]; i++)
 714             {
 715                 gb->a[gb->nra++] = gb->a[i];
 716             }
 717         }
 718         if (nion > 0)
 719         {
 720             for (i = gb->index[iion]; i < gb->index[iion+1]; i++)
 721             {
 722                 gb->a[gb->nra++] = gb->a[i];
 723             }
 724         }
 725         gb->nr++;
 726         gb->index[gb->nr] = gb->nra;
 727     }
 728 }
 729
 730
 731 void check_index(char *gname, int n, atom_id index[], char *traj, int natoms)
 732 {
 733     int i;
 734
 735     for (i = 0; i < n; i++)
 736     {
 737         if (index[i] >= natoms)
 738         {
 739             gmx_fatal(FARGS, "%s atom number (index[%d]=%d) is larger than the number of atoms in %s (%d)",
 740                       gname ? gname : "Index", i+1, index[i]+1,
 741                       traj ? traj : "the trajectory", natoms);
 742         }
 743         else if (index[i] < 0)
 744         {
 745             gmx_fatal(FARGS, "%s atom number (index[%d]=%d) is less than zero",
 746                       gname ? gname : "Index", i+1, index[i]+1);
 747         }
 748     }
 749 }
 750
 751 t_blocka *init_index(const char *gfile, char ***grpname)
 752 {
 753     FILE      *in;
 754     t_blocka  *b;
 755     int        maxentries;
 756     int        i, j;
 757     char       line[STRLEN], *pt, str[STRLEN];
 758
 759     in = gmx_fio_fopen(gfile, "r");
 760     snew(b, 1);
 761     b->nr      = 0;
 762     b->index   = NULL;
 763     b->nra     = 0;
 764     b->a       = NULL;
 765     *grpname   = NULL;
 766     maxentries = 0;
 767     while (get_a_line(in, line, STRLEN))
 768     {
 769         if (get_header(line, str))
 770         {
 771             b->nr++;
 772             srenew(b->index, b->nr+1);
 773             srenew(*grpname, b->nr);
 774             if (b->nr == 1)
 775             {
 776                 b->index[0] = 0;
 777             }
 778             b->index[b->nr]     = b->index[b->nr-1];
 779             (*grpname)[b->nr-1] = gmx_strdup(str);
 780         }
 781         else
 782         {
 783             if (b->nr == 0)
 784             {
 785                 gmx_fatal(FARGS, "The first header of your indexfile is invalid");
 786             }
 787             pt = line;
 788             while (sscanf(pt, "%s", str) == 1)
 789             {
 790                 i = b->index[b->nr];
 791                 if (i >= maxentries)
 792                 {
 793                     maxentries += 1024;
 794                     srenew(b->a, maxentries);
 795                 }
 796                 assert(b->a != NULL); // for clang analyzer
 797                 b->a[i] = strtol(str, NULL, 10)-1;
 798                 b->index[b->nr]++;
 799                 (b->nra)++;
 800                 pt = strstr(pt, str)+strlen(str);
 801             }
 802         }
 803     }
 804     gmx_fio_fclose(in);
 805
 806     for (i = 0; (i < b->nr); i++)
 807     {
 808         assert(b->a != NULL); // for clang analyzer
 809         for (j = b->index[i]; (j < b->index[i+1]); j++)
 810         {
 811             if (b->a[j] < 0)
 812             {
 813                 fprintf(stderr, "\nWARNING: negative index %d in group %s\n\n",
 814                         b->a[j], (*grpname)[i]);
 815             }
 816         }
 817     }
 818
 819     return b;
 820 }
 821
 822 static void minstring(char *str)
 823 {
 824     int i;
 825
 826     for (i = 0; (i < (int)strlen(str)); i++)
 827     {
 828         if (str[i] == '-')
 829         {
 830             str[i] = '_';
 831         }
 832     }
 833 }
 834
 835 int find_group(char s[], int ngrps, char **grpname)
 836 {
 837     int      aa, i, n;
 838     char     string[STRLEN];
 839     gmx_bool bMultiple;
 840
 841     bMultiple = FALSE;
 842     n         = strlen(s);
 843     aa        = NOTSET;
 844     /* first look for whole name match */
 845     if (aa == NOTSET)
 846     {
 847         for (i = 0; i < ngrps; i++)
 848         {
 849             if (gmx_strcasecmp_min(s, grpname[i]) == 0)
 850             {
 851                 if (aa != NOTSET)
 852                 {
 853                     bMultiple = TRUE;
 854                 }
 855                 aa = i;
 856             }
 857         }
 858     }
 859     /* second look for first string match */
 860     if (aa == NOTSET)
 861     {
 862         for (i = 0; i < ngrps; i++)
 863         {
 864             if (gmx_strncasecmp_min(s, grpname[i], n) == 0)
 865             {
 866                 if (aa != NOTSET)
 867                 {
 868                     bMultiple = TRUE;
 869                 }
 870                 aa = i;
 871             }
 872         }
 873     }
 874     /* last look for arbitrary substring match */
 875     if (aa == NOTSET)
 876     {
 877         upstring(s);
 878         minstring(s);
 879         for (i = 0; i < ngrps; i++)
 880         {
 881             strcpy(string, grpname[i]);
 882             upstring(string);
 883             minstring(string);
 884             if (strstr(string, s) != NULL)
 885             {
 886                 if (aa != NOTSET)
 887                 {
 888                     bMultiple = TRUE;
 889                 }
 890                 aa = i;
 891             }
 892         }
 893     }
 894     if (bMultiple)
 895     {
 896         printf("Error: Multiple groups '%s' selected\n", s);
 897         aa = NOTSET;
 898     }
 899     return aa;
 900 }
 901
 902 static int qgroup(int *a, int ngrps, char **grpname)
 903 {
 904     char     s[STRLEN];
 905     int      aa;
 906     gmx_bool bInRange;
 907     char    *end;
 908
 909     do
 910     {
 911         fprintf(stderr, "Select a group: ");
 912         do
 913         {
 914             if (scanf("%s", s) != 1)
 915             {
 916                 gmx_fatal(FARGS, "Cannot read from input");
 917             }
 918             trim(s); /* remove spaces */
 919         }
 920         while (strlen(s) == 0);
 921         aa = strtol(s, &end, 10);
 922         if (aa == 0 && end[0] != '\0') /* string entered */
 923         {
 924             aa = find_group(s, ngrps, grpname);
 925         }
 926         bInRange = (aa >= 0 && aa < ngrps);
 927         if (!bInRange)
 928         {
 929             printf("Error: No such group '%s'\n", s);
 930         }
 931     }
 932     while (!bInRange);
 933     printf("Selected %d: '%s'\n", aa, grpname[aa]);
 934     *a = aa;
 935     return aa;
 936 }
 937
 938 static void rd_groups(t_blocka *grps, char **grpname, char *gnames[],
 939                       int ngrps, int isize[], atom_id *index[], int grpnr[])
 940 {
 941     int i, j, gnr1;
 942
 943     if (grps->nr == 0)
 944     {
 945         gmx_fatal(FARGS, "Error: no groups in indexfile");
 946     }
 947     for (i = 0; (i < grps->nr); i++)
 948     {
 949         fprintf(stderr, "Group %5d (%15s) has %5d elements\n", i, grpname[i],
 950                 grps->index[i+1]-grps->index[i]);
 951     }
 952     for (i = 0; (i < ngrps); i++)
 953     {
 954         if (grps->nr > 1)
 955         {
 956             do
 957             {
 958                 gnr1 = qgroup(&grpnr[i], grps->nr, grpname);
 959                 if ((gnr1 < 0) || (gnr1 >= grps->nr))
 960                 {
 961                     fprintf(stderr, "Select between %d and %d.\n", 0, grps->nr-1);
 962                 }
 963             }
 964             while ((gnr1 < 0) || (gnr1 >= grps->nr));
 965         }
 966         else
 967         {
 968             fprintf(stderr, "There is one group in the index\n");
 969             gnr1 = 0;
 970         }
 971         gnames[i] = gmx_strdup(grpname[gnr1]);
 972         isize[i]  = grps->index[gnr1+1]-grps->index[gnr1];
 973         snew(index[i], isize[i]);
 974         for (j = 0; (j < isize[i]); j++)
 975         {
 976             index[i][j] = grps->a[grps->index[gnr1]+j];
 977         }
 978     }
 979 }
 980
 981 void rd_index(const char *statfile, int ngrps, int isize[],
 982               atom_id *index[], char *grpnames[])
 983 {
 984     char    **gnames;
 985     t_blocka *grps;
 986     int      *grpnr;
 987
 988     snew(grpnr, ngrps);
 989     if (!statfile)
 990     {
 991         gmx_fatal(FARGS, "No index file specified");
 992     }
 993     grps = init_index(statfile, &gnames);
 994     rd_groups(grps, gnames, grpnames, ngrps, isize, index, grpnr);
 995 }
 996
 997 void rd_index_nrs(char *statfile, int ngrps, int isize[],
 998                   atom_id *index[], char *grpnames[], int grpnr[])
 999 {
1000     char    **gnames;
1001     t_blocka *grps;
1002
1003     if (!statfile)
1004     {
1005         gmx_fatal(FARGS, "No index file specified");
1006     }
1007     grps = init_index(statfile, &gnames);
1008
1009     rd_groups(grps, gnames, grpnames, ngrps, isize, index, grpnr);
1010 }
1011
1012 void get_index(t_atoms *atoms, const char *fnm, int ngrps,
1013                int isize[], atom_id *index[], char *grpnames[])
1014 {
1015     char    ***gnames;
1016     t_blocka  *grps = NULL;
1017     int       *grpnr;
1018
1019     snew(grpnr, ngrps);
1020     snew(gnames, 1);
1021     if (fnm != NULL)
1022     {
1023         grps = init_index(fnm, gnames);
1024     }
1025     else if (atoms)
1026     {
1027         snew(grps, 1);
1028         snew(grps->index, 1);
1029         analyse(atoms, grps, gnames, FALSE, FALSE);
1030     }
1031     else
1032     {
1033         gmx_incons("You need to supply a valid atoms structure or a valid index file name");
1034     }
1035
1036     rd_groups(grps, *gnames, grpnames, ngrps, isize, index, grpnr);
1037 }
1038
1039 t_cluster_ndx *cluster_index(FILE *fplog, const char *ndx)
1040 {
1041     t_cluster_ndx *c;
1042     int            i;
1043
1044     snew(c, 1);
1045     c->clust     = init_index(ndx, &c->grpname);
1046     c->maxframe  = -1;
1047     for (i = 0; (i < c->clust->nra); i++)
1048     {
1049         c->maxframe = std::max(c->maxframe, c->clust->a[i]);
1050     }
1051     fprintf(fplog ? fplog : stdout,
1052             "There are %d clusters containing %d structures, highest framenr is %d\n",
1053             c->clust->nr, c->clust->nra, c->maxframe);
1054     if (debug)
1055     {
1056         pr_blocka(debug, 0, "clust", c->clust, TRUE);
1057         for (i = 0; (i < c->clust->nra); i++)
1058         {
1059             if ((c->clust->a[i] < 0) || (c->clust->a[i] > c->maxframe))
1060             {
1061                 gmx_fatal(FARGS, "Range check error for c->clust->a[%d] = %d\n"
1062                           "should be within 0 and %d", i, c->clust->a[i], c->maxframe+1);
1063             }
1064         }
1065     }
1066     c->inv_clust = make_invblocka(c->clust, c->maxframe);
1067
1068     return c;
1069 }