src/gromacs/gmxana/gmx_make_ndx.c

   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) 2012,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 "config.h"
  40
  41 #include <ctype.h>
  42 #include <string.h>
  43
  44 #include "gromacs/utility/futil.h"
  45 #include "gromacs/legacyheaders/macros.h"
  46 #include "gromacs/fileio/confio.h"
  47 #include "gromacs/legacyheaders/typedefs.h"
  48 #include "gromacs/topology/index.h"
  49 #include "gromacs/math/vec.h"
  50
  51 #include "gromacs/commandline/pargs.h"
  52 #include "gromacs/topology/block.h"
  53 #include "gromacs/utility/cstringutil.h"
  54 #include "gromacs/utility/fatalerror.h"
  55 #include "gromacs/utility/smalloc.h"
  56
  57 /* It's not nice to have size limits, but we should not spend more time
  58  * on this ancient tool, but instead use the new selection library.
  59  */
  60 #define MAXNAMES 1024
  61 #define NAME_LEN 1024
  62
  63 gmx_bool bCase = FALSE;
  64
  65 static int or_groups(atom_id nr1, atom_id *at1, atom_id nr2, atom_id *at2,
  66                      atom_id *nr, atom_id *at)
  67 {
  68     atom_id  i1, i2, max = 0;
  69     gmx_bool bNotIncr;
  70
  71     *nr = 0;
  72
  73     bNotIncr = FALSE;
  74     for (i1 = 0; i1 < nr1; i1++)
  75     {
  76         if ((i1 > 0) && (at1[i1] <= max))
  77         {
  78             bNotIncr = TRUE;
  79         }
  80         max = at1[i1];
  81     }
  82     for (i1 = 0; i1 < nr2; i1++)
  83     {
  84         if ((i1 > 0) && (at2[i1] <= max))
  85         {
  86             bNotIncr = TRUE;
  87         }
  88         max = at2[i1];
  89     }
  90
  91     if (bNotIncr)
  92     {
  93         printf("One of your groups is not ascending\n");
  94     }
  95     else
  96     {
  97         i1  = 0;
  98         i2  = 0;
  99         *nr = 0;
 100         while ((i1 < nr1) || (i2 < nr2))
 101         {
 102             if ((i2 == nr2) || ((i1 < nr1) && (at1[i1] < at2[i2])))
 103             {
 104                 at[*nr] = at1[i1];
 105                 (*nr)++;
 106                 i1++;
 107             }
 108             else
 109             {
 110                 if ((i2 < nr2) && ((i1 == nr1) || (at1[i1] > at2[i2])))
 111                 {
 112                     at[*nr] = at2[i2];
 113                     (*nr)++;
 114                 }
 115                 i2++;
 116             }
 117         }
 118
 119         printf("Merged two groups with OR: %d %d -> %d\n", nr1, nr2, *nr);
 120     }
 121
 122     return *nr;
 123 }
 124
 125 static int and_groups(atom_id nr1, atom_id *at1, atom_id nr2, atom_id *at2,
 126                       atom_id *nr, atom_id *at)
 127 {
 128     atom_id i1, i2;
 129
 130     *nr = 0;
 131     for (i1 = 0; i1 < nr1; i1++)
 132     {
 133         for (i2 = 0; i2 < nr2; i2++)
 134         {
 135             if (at1[i1] == at2[i2])
 136             {
 137                 at[*nr] = at1[i1];
 138                 (*nr)++;
 139             }
 140         }
 141     }
 142
 143     printf("Merged two groups with AND: %d %d -> %d\n", nr1, nr2, *nr);
 144
 145     return *nr;
 146 }
 147
 148 static gmx_bool is_name_char(char c)
 149 {
 150     /* This string should contain all characters that can not be
 151      * the first letter of a name due to the make_ndx syntax.
 152      */
 153     const char *spec = " !&|";
 154
 155     return (c != '\0' && strchr(spec, c) == NULL);
 156 }
 157
 158 static int parse_names(char **string, int *n_names, char **names)
 159 {
 160     int i;
 161
 162     *n_names = 0;
 163     while (is_name_char((*string)[0]) || (*string)[0] == ' ')
 164     {
 165         if (is_name_char((*string)[0]))
 166         {
 167             if (*n_names >= MAXNAMES)
 168             {
 169                 gmx_fatal(FARGS, "To many names: %d\n", *n_names+1);
 170             }
 171             i = 0;
 172             while (is_name_char((*string)[i]))
 173             {
 174                 names[*n_names][i] = (*string)[i];
 175                 i++;
 176                 if (i > NAME_LEN)
 177                 {
 178                     printf("Name is too long, the maximum is %d characters\n", NAME_LEN);
 179                     return 0;
 180                 }
 181             }
 182             names[*n_names][i] = '\0';
 183             if (!bCase)
 184             {
 185                 upstring(names[*n_names]);
 186             }
 187             *string += i;
 188             (*n_names)++;
 189         }
 190         else
 191         {
 192             (*string)++;
 193         }
 194     }
 195
 196     return *n_names;
 197 }
 198
 199 static gmx_bool parse_int_char(char **string, int *nr, char *c)
 200 {
 201     char    *orig;
 202     gmx_bool bRet;
 203
 204     orig = *string;
 205
 206     while ((*string)[0] == ' ')
 207     {
 208         (*string)++;
 209     }
 210
 211     bRet = FALSE;
 212
 213     *c = ' ';
 214
 215     if (isdigit((*string)[0]))
 216     {
 217         *nr = (*string)[0]-'0';
 218         (*string)++;
 219         while (isdigit((*string)[0]))
 220         {
 221             *nr = (*nr)*10+(*string)[0]-'0';
 222             (*string)++;
 223         }
 224         if (isalpha((*string)[0]))
 225         {
 226             *c = (*string)[0];
 227             (*string)++;
 228         }
 229         /* Check if there is at most one non-digit character */
 230         if (!isalnum((*string)[0]))
 231         {
 232             bRet = TRUE;
 233         }
 234         else
 235         {
 236             *string = orig;
 237         }
 238     }
 239     else
 240     {
 241         *nr = NOTSET;
 242     }
 243
 244     return bRet;
 245 }
 246
 247 static gmx_bool parse_int(char **string, int *nr)
 248 {
 249     char    *orig, c;
 250     gmx_bool bRet;
 251
 252     orig = *string;
 253     bRet = parse_int_char(string, nr, &c);
 254     if (bRet && c != ' ')
 255     {
 256         *string = orig;
 257         bRet    = FALSE;
 258     }
 259
 260     return bRet;
 261 }
 262
 263 static gmx_bool isquote(char c)
 264 {
 265     return (c == '\"');
 266 }
 267
 268 static gmx_bool parse_string(char **string, int *nr, int ngrps, char **grpname)
 269 {
 270     char *s, *sp;
 271     char  c;
 272
 273     while ((*string)[0] == ' ')
 274     {
 275         (*string)++;
 276     }
 277
 278     (*nr) = NOTSET;
 279     if (isquote((*string)[0]))
 280     {
 281         c = (*string)[0];
 282         (*string)++;
 283         s  = gmx_strdup((*string));
 284         sp = strchr(s, c);
 285         if (sp != NULL)
 286         {
 287             (*string) += sp-s + 1;
 288             sp[0]      = '\0';
 289             (*nr)      = find_group(s, ngrps, grpname);
 290         }
 291     }
 292
 293     return (*nr) != NOTSET;
 294 }
 295
 296 static int select_atomnumbers(char **string, t_atoms *atoms, atom_id n1,
 297                               atom_id *nr, atom_id *index, char *gname)
 298 {
 299     char    buf[STRLEN];
 300     int     i, up;
 301
 302     *nr = 0;
 303     while ((*string)[0] == ' ')
 304     {
 305         (*string)++;
 306     }
 307     if ((*string)[0] == '-')
 308     {
 309         (*string)++;
 310         parse_int(string, &up);
 311         if ((n1 < 1) || (n1 > atoms->nr) || (up < 1) || (up > atoms->nr))
 312         {
 313             printf("Invalid atom range\n");
 314         }
 315         else
 316         {
 317             for (i = n1-1; i <= up-1; i++)
 318             {
 319                 index[*nr] = i;
 320                 (*nr)++;
 321             }
 322             printf("Found %d atom%s in range %d-%d\n", *nr, (*nr == 1) ? "" : "s", n1, up);
 323             if (n1 == up)
 324             {
 325                 sprintf(buf, "a_%d", n1);
 326             }
 327             else
 328             {
 329                 sprintf(buf, "a_%d-%d", n1, up);
 330             }
 331             strcpy(gname, buf);
 332         }
 333     }
 334     else
 335     {
 336         i = n1;
 337         sprintf(gname, "a");
 338         do
 339         {
 340             if ((i-1 >= 0) && (i-1 < atoms->nr))
 341             {
 342                 index[*nr] = i-1;
 343                 (*nr)++;
 344                 sprintf(buf, "_%d", i);
 345                 strcat(gname, buf);
 346             }
 347             else
 348             {
 349                 printf("Invalid atom number %d\n", i);
 350                 *nr = 0;
 351             }
 352         }
 353         while ((*nr != 0) && (parse_int(string, &i)));
 354     }
 355
 356     return *nr;
 357 }
 358
 359 static int select_residuenumbers(char **string, t_atoms *atoms,
 360                                  atom_id n1, char c,
 361                                  atom_id *nr, atom_id *index, char *gname)
 362 {
 363     char       buf[STRLEN];
 364     int        i, j, up;
 365     t_resinfo *ri;
 366
 367     *nr = 0;
 368     while ((*string)[0] == ' ')
 369     {
 370         (*string)++;
 371     }
 372     if ((*string)[0] == '-')
 373     {
 374         /* Residue number range selection */
 375         if (c != ' ')
 376         {
 377             printf("Error: residue insertion codes can not be used with residue range selection\n");
 378             return 0;
 379         }
 380         (*string)++;
 381         parse_int(string, &up);
 382
 383         for (i = 0; i < atoms->nr; i++)
 384         {
 385             ri = &atoms->resinfo[atoms->atom[i].resind];
 386             for (j = n1; (j <= up); j++)
 387             {
 388                 if (ri->nr == j && (c == ' ' || ri->ic == c))
 389                 {
 390                     index[*nr] = i;
 391                     (*nr)++;
 392                 }
 393             }
 394         }
 395         printf("Found %d atom%s with res.nr. in range %d-%d\n",
 396                *nr, (*nr == 1) ? "" : "s", n1, up);
 397         if (n1 == up)
 398         {
 399             sprintf(buf, "r_%d", n1);
 400         }
 401         else
 402         {
 403             sprintf(buf, "r_%d-%d", n1, up);
 404         }
 405         strcpy(gname, buf);
 406     }
 407     else
 408     {
 409         /* Individual residue number/insertion code selection */
 410         j = n1;
 411         sprintf(gname, "r");
 412         do
 413         {
 414             for (i = 0; i < atoms->nr; i++)
 415             {
 416                 ri = &atoms->resinfo[atoms->atom[i].resind];
 417                 if (ri->nr == j && ri->ic == c)
 418                 {
 419                     index[*nr] = i;
 420                     (*nr)++;
 421                 }
 422             }
 423             sprintf(buf, "_%d", j);
 424             strcat(gname, buf);
 425         }
 426         while (parse_int_char(string, &j, &c));
 427     }
 428
 429     return *nr;
 430 }
 431
 432 static int select_residueindices(char **string, t_atoms *atoms,
 433                                  atom_id n1, char c,
 434                                  atom_id *nr, atom_id *index, char *gname)
 435 {
 436     /*this should be similar to select_residuenumbers except select by index (sequential numbering in file)*/
 437     /*resind+1 for 1-indexing*/
 438     char       buf[STRLEN];
 439     int        i, j, up;
 440     t_resinfo *ri;
 441
 442     *nr = 0;
 443     while ((*string)[0] == ' ')
 444     {
 445         (*string)++;
 446     }
 447     if ((*string)[0] == '-')
 448     {
 449         /* Residue number range selection */
 450         if (c != ' ')
 451         {
 452             printf("Error: residue insertion codes can not be used with residue range selection\n");
 453             return 0;
 454         }
 455         (*string)++;
 456         parse_int(string, &up);
 457
 458         for (i = 0; i < atoms->nr; i++)
 459         {
 460             ri = &atoms->resinfo[atoms->atom[i].resind];
 461             for (j = n1; (j <= up); j++)
 462             {
 463                 if (atoms->atom[i].resind+1 == j && (c == ' ' || ri->ic == c))
 464                 {
 465                     index[*nr] = i;
 466                     (*nr)++;
 467                 }
 468             }
 469         }
 470         printf("Found %d atom%s with resind.+1 in range %d-%d\n",
 471                *nr, (*nr == 1) ? "" : "s", n1, up);
 472         if (n1 == up)
 473         {
 474             sprintf(buf, "r_%d", n1);
 475         }
 476         else
 477         {
 478             sprintf(buf, "r_%d-%d", n1, up);
 479         }
 480         strcpy(gname, buf);
 481     }
 482     else
 483     {
 484         /* Individual residue number/insertion code selection */
 485         j = n1;
 486         sprintf(gname, "r");
 487         do
 488         {
 489             for (i = 0; i < atoms->nr; i++)
 490             {
 491                 ri = &atoms->resinfo[atoms->atom[i].resind];
 492                 if (atoms->atom[i].resind+1 == j && ri->ic == c)
 493                 {
 494                     index[*nr] = i;
 495                     (*nr)++;
 496                 }
 497             }
 498             sprintf(buf, "_%d", j);
 499             strcat(gname, buf);
 500         }
 501         while (parse_int_char(string, &j, &c));
 502     }
 503
 504     return *nr;
 505 }
 506
 507
 508 static gmx_bool atoms_from_residuenumbers(t_atoms *atoms, int group, t_blocka *block,
 509                                           atom_id *nr, atom_id *index, char *gname)
 510 {
 511     int i, j, j0, j1, resnr, nres;
 512
 513     j0   = block->index[group];
 514     j1   = block->index[group+1];
 515     nres = atoms->nres;
 516     for (j = j0; j < j1; j++)
 517     {
 518         if (block->a[j] >= nres)
 519         {
 520             printf("Index %s contains number>nres (%d>%d)\n",
 521                    gname, block->a[j]+1, nres);
 522             return FALSE;
 523         }
 524     }
 525     for (i = 0; i < atoms->nr; i++)
 526     {
 527         resnr = atoms->resinfo[atoms->atom[i].resind].nr;
 528         for (j = j0; j < j1; j++)
 529         {
 530             if (block->a[j]+1 == resnr)
 531             {
 532                 index[*nr] = i;
 533                 (*nr)++;
 534                 break;
 535             }
 536         }
 537     }
 538     printf("Found %d atom%s in %d residues from group %s\n",
 539            *nr, (*nr == 1) ? "" : "s", j1-j0, gname);
 540     return *nr;
 541 }
 542
 543 static gmx_bool comp_name(char *name, char *search)
 544 {
 545     while (name[0] != '\0' && search[0] != '\0')
 546     {
 547         switch (search[0])
 548         {
 549             case '?':
 550                 /* Always matches */
 551                 break;
 552             case '*':
 553                 if (search[1] != '\0')
 554                 {
 555                     printf("WARNING: Currently '*' is only supported at the end of an expression\n");
 556                     return FALSE;
 557                 }
 558                 else
 559                 {
 560                     return TRUE;
 561                 }
 562                 break;
 563             default:
 564                 /* Compare a single character */
 565                 if (( bCase && strncmp(name, search, 1)) ||
 566                     (!bCase && gmx_strncasecmp(name, search, 1)))
 567                 {
 568                     return FALSE;
 569                 }
 570         }
 571         name++;
 572         search++;
 573     }
 574
 575     return (name[0] == '\0' && (search[0] == '\0' || search[0] == '*'));
 576 }
 577
 578 static int select_chainnames(t_atoms *atoms, int n_names, char **names,
 579                              atom_id *nr, atom_id *index)
 580 {
 581     char    name[2];
 582     int     j;
 583     atom_id i;
 584
 585     name[1] = 0;
 586     *nr     = 0;
 587     for (i = 0; i < atoms->nr; i++)
 588     {
 589         name[0] = atoms->resinfo[atoms->atom[i].resind].chainid;
 590         j       = 0;
 591         while (j < n_names && !comp_name(name, names[j]))
 592         {
 593             j++;
 594         }
 595         if (j < n_names)
 596         {
 597             index[*nr] = i;
 598             (*nr)++;
 599         }
 600     }
 601     printf("Found %d atom%s with chain identifier%s",
 602            *nr, (*nr == 1) ? "" : "s", (n_names == 1) ? "" : "s");
 603     for (j = 0; (j < n_names); j++)
 604     {
 605         printf(" %s", names[j]);
 606     }
 607     printf("\n");
 608
 609     return *nr;
 610 }
 611
 612 static int select_atomnames(t_atoms *atoms, int n_names, char **names,
 613                             atom_id *nr, atom_id *index, gmx_bool bType)
 614 {
 615     char   *name;
 616     int     j;
 617     atom_id i;
 618
 619     *nr = 0;
 620     for (i = 0; i < atoms->nr; i++)
 621     {
 622         if (bType)
 623         {
 624             name = *(atoms->atomtype[i]);
 625         }
 626         else
 627         {
 628             name = *(atoms->atomname[i]);
 629         }
 630         j = 0;
 631         while (j < n_names && !comp_name(name, names[j]))
 632         {
 633             j++;
 634         }
 635         if (j < n_names)
 636         {
 637             index[*nr] = i;
 638             (*nr)++;
 639         }
 640     }
 641     printf("Found %d atoms with %s%s",
 642            *nr, bType ? "type" : "name", (n_names == 1) ? "" : "s");
 643     for (j = 0; (j < n_names); j++)
 644     {
 645         printf(" %s", names[j]);
 646     }
 647     printf("\n");
 648
 649     return *nr;
 650 }
 651
 652 static int select_residuenames(t_atoms *atoms, int n_names, char **names,
 653                                atom_id *nr, atom_id *index)
 654 {
 655     char   *name;
 656     int     j;
 657     atom_id i;
 658
 659     *nr = 0;
 660     for (i = 0; i < atoms->nr; i++)
 661     {
 662         name = *(atoms->resinfo[atoms->atom[i].resind].name);
 663         j    = 0;
 664         while (j < n_names && !comp_name(name, names[j]))
 665         {
 666             j++;
 667         }
 668         if (j < n_names)
 669         {
 670             index[*nr] = i;
 671             (*nr)++;
 672         }
 673     }
 674     printf("Found %d atoms with residue name%s", *nr, (n_names == 1) ? "" : "s");
 675     for (j = 0; (j < n_names); j++)
 676     {
 677         printf(" %s", names[j]);
 678     }
 679     printf("\n");
 680
 681     return *nr;
 682 }
 683
 684 static void copy2block(int n, atom_id *index, t_blocka *block)
 685 {
 686     int i, n0;
 687
 688     block->nr++;
 689     n0         = block->nra;
 690     block->nra = n0+n;
 691     srenew(block->index, block->nr+1);
 692     block->index[block->nr] = n0+n;
 693     srenew(block->a, n0+n);
 694     for (i = 0; (i < n); i++)
 695     {
 696         block->a[n0+i] = index[i];
 697     }
 698 }
 699
 700 static void make_gname(int n, char **names, char *gname)
 701 {
 702     int i;
 703
 704     strcpy(gname, names[0]);
 705     for (i = 1; i < n; i++)
 706     {
 707         strcat(gname, "_");
 708         strcat(gname, names[i]);
 709     }
 710 }
 711
 712 static void copy_group(int g, t_blocka *block, atom_id *nr, atom_id *index)
 713 {
 714     int i, i0;
 715
 716     i0  = block->index[g];
 717     *nr = block->index[g+1]-i0;
 718     for (i = 0; i < *nr; i++)
 719     {
 720         index[i] = block->a[i0+i];
 721     }
 722 }
 723
 724 static void remove_group(int nr, int nr2, t_blocka *block, char ***gn)
 725 {
 726     int   i, j, shift;
 727     char *name;
 728
 729     if (nr2 == NOTSET)
 730     {
 731         nr2 = nr;
 732     }
 733
 734     for (j = 0; j <= nr2-nr; j++)
 735     {
 736         if ((nr < 0) || (nr >= block->nr))
 737         {
 738             printf("Group %d does not exist\n", nr+j);
 739         }
 740         else
 741         {
 742             shift = block->index[nr+1]-block->index[nr];
 743             for (i = block->index[nr+1]; i < block->nra; i++)
 744             {
 745                 block->a[i-shift] = block->a[i];
 746             }
 747
 748             for (i = nr; i < block->nr; i++)
 749             {
 750                 block->index[i] = block->index[i+1]-shift;
 751             }
 752             name = gmx_strdup((*gn)[nr]);
 753             sfree((*gn)[nr]);
 754             for (i = nr; i < block->nr-1; i++)
 755             {
 756                 (*gn)[i] = (*gn)[i+1];
 757             }
 758             block->nr--;
 759             block->nra = block->index[block->nr];
 760             printf("Removed group %d '%s'\n", nr+j, name);
 761             sfree(name);
 762         }
 763     }
 764 }
 765
 766 static void split_group(t_atoms *atoms, int sel_nr, t_blocka *block, char ***gn,
 767                         gmx_bool bAtom)
 768 {
 769     char    buf[STRLEN], *name;
 770     int     i, resind;
 771     atom_id a, n0, n1;
 772
 773     printf("Splitting group %d '%s' into %s\n", sel_nr, (*gn)[sel_nr],
 774            bAtom ? "atoms" : "residues");
 775
 776     n0 = block->index[sel_nr];
 777     n1 = block->index[sel_nr+1];
 778     srenew(block->a, block->nra+n1-n0);
 779     for (i = n0; i < n1; i++)
 780     {
 781         a      = block->a[i];
 782         resind = atoms->atom[a].resind;
 783         name   = *(atoms->resinfo[resind].name);
 784         if (bAtom || (i == n0) || (atoms->atom[block->a[i-1]].resind != resind))
 785         {
 786             if (i > n0)
 787             {
 788                 block->index[block->nr] = block->nra;
 789             }
 790             block->nr++;
 791             srenew(block->index, block->nr+1);
 792             srenew(*gn, block->nr);
 793             if (bAtom)
 794             {
 795                 sprintf(buf, "%s_%s_%d", (*gn)[sel_nr], *atoms->atomname[a], a+1);
 796             }
 797             else
 798             {
 799                 sprintf(buf, "%s_%s_%d", (*gn)[sel_nr], name, atoms->resinfo[resind].nr);
 800             }
 801             (*gn)[block->nr-1] = gmx_strdup(buf);
 802         }
 803         block->a[block->nra] = a;
 804         block->nra++;
 805     }
 806     block->index[block->nr] = block->nra;
 807 }
 808
 809 static int split_chain(t_atoms *atoms, rvec *x,
 810                        int sel_nr, t_blocka *block, char ***gn)
 811 {
 812     char    buf[STRLEN];
 813     int     j, nchain;
 814     atom_id i, a, natoms, *start = NULL, *end = NULL, ca_start, ca_end;
 815     rvec    vec;
 816
 817     natoms   = atoms->nr;
 818     nchain   = 0;
 819     ca_start = 0;
 820
 821     while (ca_start < natoms)
 822     {
 823         while ((ca_start < natoms) && strcmp(*atoms->atomname[ca_start], "CA"))
 824         {
 825             ca_start++;
 826         }
 827         if (ca_start < natoms)
 828         {
 829             srenew(start, nchain+1);
 830             srenew(end, nchain+1);
 831             start[nchain] = ca_start;
 832             while ((start[nchain] > 0) &&
 833                    (atoms->atom[start[nchain]-1].resind ==
 834                     atoms->atom[ca_start].resind))
 835             {
 836                 start[nchain]--;
 837             }
 838
 839             i = ca_start;
 840             do
 841             {
 842                 ca_end = i;
 843                 do
 844                 {
 845                     i++;
 846                 }
 847                 while ((i < natoms) && strcmp(*atoms->atomname[i], "CA"));
 848                 if (i < natoms)
 849                 {
 850                     rvec_sub(x[ca_end], x[i], vec);
 851                 }
 852                 else
 853                 {
 854                     break;
 855                 }
 856             }
 857             while (norm(vec) < 0.45);
 858
 859             end[nchain] = ca_end;
 860             while ((end[nchain]+1 < natoms) &&
 861                    (atoms->atom[end[nchain]+1].resind == atoms->atom[ca_end].resind))
 862             {
 863                 end[nchain]++;
 864             }
 865             ca_start = end[nchain]+1;
 866             nchain++;
 867         }
 868     }
 869     if (nchain == 1)
 870     {
 871         printf("Found 1 chain, will not split\n");
 872     }
 873     else
 874     {
 875         printf("Found %d chains\n", nchain);
 876     }
 877     for (j = 0; j < nchain; j++)
 878     {
 879         printf("%d:%6d atoms (%d to %d)\n",
 880                j+1, end[j]-start[j]+1, start[j]+1, end[j]+1);
 881     }
 882
 883     if (nchain > 1)
 884     {
 885         srenew(block->a, block->nra+block->index[sel_nr+1]-block->index[sel_nr]);
 886         for (j = 0; j < nchain; j++)
 887         {
 888             block->nr++;
 889             srenew(block->index, block->nr+1);
 890             srenew(*gn, block->nr);
 891             sprintf(buf, "%s_chain%d", (*gn)[sel_nr], j+1);
 892             (*gn)[block->nr-1] = gmx_strdup(buf);
 893             for (i = block->index[sel_nr]; i < block->index[sel_nr+1]; i++)
 894             {
 895                 a = block->a[i];
 896                 if ((a >= start[j]) && (a <= end[j]))
 897                 {
 898                     block->a[block->nra] = a;
 899                     block->nra++;
 900                 }
 901             }
 902             block->index[block->nr] = block->nra;
 903             if (block->index[block->nr-1] == block->index[block->nr])
 904             {
 905                 remove_group(block->nr-1, NOTSET, block, gn);
 906             }
 907         }
 908     }
 909     sfree(start);
 910     sfree(end);
 911
 912     return nchain;
 913 }
 914
 915 static gmx_bool check_have_atoms(t_atoms *atoms, char *string)
 916 {
 917     if (atoms == NULL)
 918     {
 919         printf("Can not process '%s' without atom info, use option -f\n", string);
 920         return FALSE;
 921     }
 922     else
 923     {
 924         return TRUE;
 925     }
 926 }
 927
 928 static gmx_bool parse_entry(char **string, int natoms, t_atoms *atoms,
 929                             t_blocka *block, char ***gn,
 930                             atom_id *nr, atom_id *index, char *gname)
 931 {
 932     static char   **names, *ostring;
 933     static gmx_bool bFirst = TRUE;
 934     int             j, n_names, sel_nr1;
 935     atom_id         i, nr1, *index1;
 936     char            c;
 937     gmx_bool        bRet, bCompl;
 938
 939     if (bFirst)
 940     {
 941         bFirst = FALSE;
 942         snew(names, MAXNAMES);
 943         for (i = 0; i < MAXNAMES; i++)
 944         {
 945             snew(names[i], NAME_LEN+1);
 946         }
 947     }
 948
 949     bRet    = FALSE;
 950     sel_nr1 = NOTSET;
 951
 952     while (*string[0] == ' ')
 953     {
 954         (*string)++;
 955     }
 956
 957     if ((*string)[0] == '!')
 958     {
 959         bCompl = TRUE;
 960         (*string)++;
 961         while (*string[0] == ' ')
 962         {
 963             (*string)++;
 964         }
 965     }
 966     else
 967     {
 968         bCompl = FALSE;
 969     }
 970
 971     ostring = *string;
 972
 973     if (parse_int(string, &sel_nr1) ||
 974         parse_string(string, &sel_nr1, block->nr, *gn))
 975     {
 976         if ((sel_nr1 >= 0) && (sel_nr1 < block->nr))
 977         {
 978             copy_group(sel_nr1, block, nr, index);
 979             strcpy(gname, (*gn)[sel_nr1]);
 980             printf("Copied index group %d '%s'\n", sel_nr1, (*gn)[sel_nr1]);
 981             bRet = TRUE;
 982         }
 983         else
 984         {
 985             printf("Group %d does not exist\n", sel_nr1);
 986         }
 987     }
 988     else if ((*string)[0] == 'a')
 989     {
 990         (*string)++;
 991         if (check_have_atoms(atoms, ostring))
 992         {
 993             if (parse_int(string, &sel_nr1))
 994             {
 995                 bRet = select_atomnumbers(string, atoms, sel_nr1, nr, index, gname);
 996             }
 997             else if (parse_names(string, &n_names, names))
 998             {
 999                 bRet = select_atomnames(atoms, n_names, names, nr, index, FALSE);
1000                 make_gname(n_names, names, gname);
1001             }
1002         }
1003     }
1004     else if ((*string)[0] == 't')
1005     {
1006         (*string)++;
1007         if (check_have_atoms(atoms, ostring) &&
1008             parse_names(string, &n_names, names))
1009         {
1010             if (atoms->atomtype == NULL)
1011             {
1012                 printf("Need a run input file to select atom types\n");
1013             }
1014             else
1015             {
1016                 bRet = select_atomnames(atoms, n_names, names, nr, index, TRUE);
1017                 make_gname(n_names, names, gname);
1018             }
1019         }
1020     }
1021     else if (strncmp(*string, "res", 3) == 0)
1022     {
1023         (*string) += 3;
1024         if (check_have_atoms(atoms, ostring) &&
1025             parse_int(string, &sel_nr1) &&
1026             (sel_nr1 >= 0) && (sel_nr1 < block->nr) )
1027         {
1028             bRet = atoms_from_residuenumbers(atoms,
1029                                              sel_nr1, block, nr, index, (*gn)[sel_nr1]);
1030             sprintf(gname, "atom_%s", (*gn)[sel_nr1]);
1031         }
1032     }
1033     else if (strncmp(*string, "ri", 2) == 0)
1034     {
1035         (*string) += 2;
1036         if (check_have_atoms(atoms, ostring) &&
1037             parse_int_char(string, &sel_nr1, &c))
1038         {
1039             bRet = select_residueindices(string, atoms, sel_nr1, c, nr, index, gname);
1040         }
1041     }
1042     else if ((*string)[0] == 'r')
1043     {
1044         (*string)++;
1045         if (check_have_atoms(atoms, ostring))
1046         {
1047             if (parse_int_char(string, &sel_nr1, &c))
1048             {
1049                 bRet = select_residuenumbers(string, atoms, sel_nr1, c, nr, index, gname);
1050             }
1051             else if (parse_names(string, &n_names, names))
1052             {
1053                 bRet = select_residuenames(atoms, n_names, names, nr, index);
1054                 make_gname(n_names, names, gname);
1055             }
1056         }
1057     }
1058     else if (strncmp(*string, "chain", 5) == 0)
1059     {
1060         (*string) += 5;
1061         if (check_have_atoms(atoms, ostring) &&
1062             parse_names(string, &n_names, names))
1063         {
1064             bRet = select_chainnames(atoms, n_names, names, nr, index);
1065             sprintf(gname, "ch%s", names[0]);
1066             for (i = 1; i < n_names; i++)
1067             {
1068                 strcat(gname, names[i]);
1069             }
1070         }
1071     }
1072     if (bRet && bCompl)
1073     {
1074         snew(index1, natoms-*nr);
1075         nr1 = 0;
1076         for (i = 0; i < natoms; i++)
1077         {
1078             j = 0;
1079             while ((j < *nr) && (index[j] != i))
1080             {
1081                 j++;
1082             }
1083             if (j == *nr)
1084             {
1085                 if (nr1 >= natoms-*nr)
1086                 {
1087                     printf("There are double atoms in your index group\n");
1088                     break;
1089                 }
1090                 index1[nr1] = i;
1091                 nr1++;
1092             }
1093         }
1094         *nr = nr1;
1095         for (i = 0; i < nr1; i++)
1096         {
1097             index[i] = index1[i];
1098         }
1099         sfree(index1);
1100
1101         for (i = strlen(gname)+1; i > 0; i--)
1102         {
1103             gname[i] = gname[i-1];
1104         }
1105         gname[0] = '!';
1106         printf("Complemented group: %d atoms\n", *nr);
1107     }
1108
1109     return bRet;
1110 }
1111
1112 static void list_residues(t_atoms *atoms)
1113 {
1114     int      i, j, start, end, prev_resind, resind;
1115     gmx_bool bDiff;
1116
1117     /* Print all the residues, assuming continuous resnr count */
1118     start       = atoms->atom[0].resind;
1119     prev_resind = start;
1120     for (i = 0; i < atoms->nr; i++)
1121     {
1122         resind = atoms->atom[i].resind;
1123         if ((resind != prev_resind) || (i == atoms->nr-1))
1124         {
1125             if ((bDiff = strcmp(*atoms->resinfo[resind].name,
1126                                 *atoms->resinfo[start].name)) ||
1127                 (i == atoms->nr-1))
1128             {
1129                 if (bDiff)
1130                 {
1131                     end = prev_resind;
1132                 }
1133                 else
1134                 {
1135                     end = resind;
1136                 }
1137                 if (end < start+3)
1138                 {
1139                     for (j = start; j <= end; j++)
1140                     {
1141                         printf("%4d %-5s",
1142                                j+1, *(atoms->resinfo[j].name));
1143                     }
1144                 }
1145                 else
1146                 {
1147                     printf(" %4d - %4d %-5s  ",
1148                            start+1, end+1, *(atoms->resinfo[start].name));
1149                 }
1150                 start = resind;
1151             }
1152         }
1153         prev_resind = resind;
1154     }
1155     printf("\n");
1156 }
1157
1158 static void edit_index(int natoms, t_atoms *atoms, rvec *x, t_blocka *block, char ***gn, gmx_bool bVerbose)
1159 {
1160     static char   **atnames, *ostring;
1161     static gmx_bool bFirst = TRUE;
1162     char            inp_string[STRLEN], *string;
1163     char            gname[STRLEN], gname1[STRLEN], gname2[STRLEN];
1164     int             i, i0, i1, sel_nr, sel_nr2, newgroup;
1165     atom_id         nr, nr1, nr2, *index, *index1, *index2;
1166     gmx_bool        bAnd, bOr, bPrintOnce;
1167
1168     if (bFirst)
1169     {
1170         bFirst = FALSE;
1171         snew(atnames, MAXNAMES);
1172         for (i = 0; i < MAXNAMES; i++)
1173         {
1174             snew(atnames[i], NAME_LEN+1);
1175         }
1176     }
1177
1178     string = NULL;
1179
1180     snew(index, natoms);
1181     snew(index1, natoms);
1182     snew(index2, natoms);
1183
1184     newgroup   = NOTSET;
1185     bPrintOnce = TRUE;
1186     do
1187     {
1188         gname1[0] = '\0';
1189         if (bVerbose || bPrintOnce || newgroup != NOTSET)
1190         {
1191             printf("\n");
1192             if (bVerbose || bPrintOnce || newgroup == NOTSET)
1193             {
1194                 i0 = 0;
1195                 i1 = block->nr;
1196             }
1197             else
1198             {
1199                 i0 = newgroup;
1200                 i1 = newgroup+1;
1201             }
1202             for (i = i0; i < i1; i++)
1203             {
1204                 printf("%3d %-20s: %5d atoms\n", i, (*gn)[i],
1205                        block->index[i+1]-block->index[i]);
1206             }
1207             newgroup = NOTSET;
1208         }
1209         if (bVerbose || bPrintOnce)
1210         {
1211             printf("\n");
1212             printf(" nr : group       !   'name' nr name   'splitch' nr    Enter: list groups\n");
1213             printf(" 'a': atom        &   'del' nr         'splitres' nr   'l': list residues\n");
1214             printf(" 't': atom type   |   'keep' nr        'splitat' nr    'h': help\n");
1215             printf(" 'r': residue         'res' nr         'chain' char\n");
1216             printf(" \"name\": group        'case': case %s         'q': save and quit\n",
1217                    bCase ? "insensitive" : "sensitive  ");
1218             printf(" 'ri': residue index\n");
1219             bPrintOnce = FALSE;
1220         }
1221         printf("\n");
1222         printf("> ");
1223         if (NULL == fgets(inp_string, STRLEN, stdin))
1224         {
1225             gmx_fatal(FARGS, "Error reading user input");
1226         }
1227         inp_string[strlen(inp_string)-1] = 0;
1228         printf("\n");
1229         string = inp_string;
1230         while (string[0] == ' ')
1231         {
1232             string++;
1233         }
1234
1235         ostring = string;
1236         nr      = 0;
1237         if (string[0] == 'h')
1238         {
1239             printf(" nr                : selects an index group by number or quoted string.\n");
1240             printf("                     The string is first matched against the whole group name,\n");
1241             printf("                     then against the beginning and finally against an\n");
1242             printf("                     arbitrary substring. A multiple match is an error.\n");
1243
1244             printf(" 'a' nr1 [nr2 ...] : selects atoms, atom numbering starts at 1.\n");
1245             printf(" 'a' nr1 - nr2     : selects atoms in the range from nr1 to nr2.\n");
1246             printf(" 'a' name1[*] [name2[*] ...] : selects atoms by name(s), '?' matches any char,\n");
1247             printf("                               wildcard '*' allowed at the end of a name.\n");
1248             printf(" 't' type1[*] [type2[*] ...] : as 'a', but for type, run input file required.\n");
1249             printf(" 'r' nr1[ic1] [nr2[ic2] ...] : selects residues by number and insertion code.\n");
1250             printf(" 'r' nr1 - nr2               : selects residues in the range from nr1 to nr2.\n");
1251             printf(" 'r' name1[*] [name2[*] ...] : as 'a', but for residue names.\n");
1252             printf(" 'ri' nr1 - nr2              : selects residue indices, 1-indexed, (as opposed to numbers) in the range from nr1 to nr2.\n");
1253             printf(" 'chain' ch1 [ch2 ...]       : selects atoms by chain identifier(s),\n");
1254             printf("                               not available with a .gro file as input.\n");
1255             printf(" !                 : takes the complement of a group with respect to all\n");
1256             printf("                     the atoms in the input file.\n");
1257             printf(" & |               : AND and OR, can be placed between any of the options\n");
1258             printf("                     above, the input is processed from left to right.\n");
1259             printf(" 'name' nr name    : rename group nr to name.\n");
1260             printf(" 'del' nr1 [- nr2] : deletes one group or groups in the range from nr1 to nr2.\n");
1261             printf(" 'keep' nr         : deletes all groups except nr.\n");
1262             printf(" 'case'            : make all name compares case (in)sensitive.\n");
1263             printf(" 'splitch' nr      : split group into chains using CA distances.\n");
1264             printf(" 'splitres' nr     : split group into residues.\n");
1265             printf(" 'splitat' nr      : split group into atoms.\n");
1266             printf(" 'res' nr          : interpret numbers in group as residue numbers\n");
1267             printf(" Enter             : list the currently defined groups and commands\n");
1268             printf(" 'l'               : list the residues.\n");
1269             printf(" 'h'               : show this help.\n");
1270             printf(" 'q'               : save and quit.\n");
1271             printf("\n");
1272             printf(" Examples:\n");
1273             printf(" > 2 | 4 & r 3-5\n");
1274             printf(" selects all atoms from group 2 and 4 that have residue numbers 3, 4 or 5\n");
1275             printf(" > a C* & !a C CA\n");
1276             printf(" selects all atoms starting with 'C' but not the atoms 'C' and 'CA'\n");
1277             printf(" > \"protein\" & ! \"backb\"\n");
1278             printf(" selects all atoms that are in group 'protein' and not in group 'backbone'\n");
1279             if (bVerbose)
1280             {
1281                 printf("\npress Enter ");
1282                 getchar();
1283             }
1284         }
1285         else if (strncmp(string, "del", 3) == 0)
1286         {
1287             string += 3;
1288             if (parse_int(&string, &sel_nr))
1289             {
1290                 while (string[0] == ' ')
1291                 {
1292                     string++;
1293                 }
1294                 if (string[0] == '-')
1295                 {
1296                     string++;
1297                     parse_int(&string, &sel_nr2);
1298                 }
1299                 else
1300                 {
1301                     sel_nr2 = NOTSET;
1302                 }
1303                 while (string[0] == ' ')
1304                 {
1305                     string++;
1306                 }
1307                 if (string[0] == '\0')
1308                 {
1309                     remove_group(sel_nr, sel_nr2, block, gn);
1310                 }
1311                 else
1312                 {
1313                     printf("\nSyntax error: \"%s\"\n", string);
1314                 }
1315             }
1316         }
1317         else if (strncmp(string, "keep", 4) == 0)
1318         {
1319             string += 4;
1320             if (parse_int(&string, &sel_nr))
1321             {
1322                 remove_group(sel_nr+1, block->nr-1, block, gn);
1323                 remove_group(0, sel_nr-1, block, gn);
1324             }
1325         }
1326         else if (strncmp(string, "name", 4) == 0)
1327         {
1328             string += 4;
1329             if (parse_int(&string, &sel_nr))
1330             {
1331                 if ((sel_nr >= 0) && (sel_nr < block->nr))
1332                 {
1333                     sscanf(string, "%s", gname);
1334                     sfree((*gn)[sel_nr]);
1335                     (*gn)[sel_nr] = gmx_strdup(gname);
1336                 }
1337             }
1338         }
1339         else if (strncmp(string, "case", 4) == 0)
1340         {
1341             bCase = !bCase;
1342             printf("Switched to case %s\n", bCase ? "sensitive" : "insensitive");
1343         }
1344         else if (string[0] == 'v')
1345         {
1346             bVerbose = !bVerbose;
1347             printf("Turned verbose %s\n", bVerbose ? "on" : "off");
1348         }
1349         else if (string[0] == 'l')
1350         {
1351             if (check_have_atoms(atoms, ostring) )
1352             {
1353                 list_residues(atoms);
1354             }
1355         }
1356         else if (strncmp(string, "splitch", 7) == 0)
1357         {
1358             string += 7;
1359             if (check_have_atoms(atoms, ostring) &&
1360                 parse_int(&string, &sel_nr) &&
1361                 (sel_nr >= 0) && (sel_nr < block->nr))
1362             {
1363                 split_chain(atoms, x, sel_nr, block, gn);
1364             }
1365         }
1366         else if (strncmp(string, "splitres", 8) == 0)
1367         {
1368             string += 8;
1369             if (check_have_atoms(atoms, ostring) &&
1370                 parse_int(&string, &sel_nr) &&
1371                 (sel_nr >= 0) && (sel_nr < block->nr))
1372             {
1373                 split_group(atoms, sel_nr, block, gn, FALSE);
1374             }
1375         }
1376         else if (strncmp(string, "splitat", 7) == 0)
1377         {
1378             string += 7;
1379             if (check_have_atoms(atoms, ostring) &&
1380                 parse_int(&string, &sel_nr) &&
1381                 (sel_nr >= 0) && (sel_nr < block->nr))
1382             {
1383                 split_group(atoms, sel_nr, block, gn, TRUE);
1384             }
1385         }
1386         else if (string[0] == '\0')
1387         {
1388             bPrintOnce = TRUE;
1389         }
1390         else if (string[0] != 'q')
1391         {
1392             nr1 = -1;
1393             nr2 = -1;
1394             if (parse_entry(&string, natoms, atoms, block, gn, &nr, index, gname))
1395             {
1396                 do
1397                 {
1398                     while (string[0] == ' ')
1399                     {
1400                         string++;
1401                     }
1402
1403                     bAnd = FALSE;
1404                     bOr  = FALSE;
1405                     if (string[0] == '&')
1406                     {
1407                         bAnd = TRUE;
1408                     }
1409                     else if (string[0] == '|')
1410                     {
1411                         bOr = TRUE;
1412                     }
1413
1414                     if (bAnd || bOr)
1415                     {
1416                         string++;
1417                         nr1 = nr;
1418                         for (i = 0; i < nr; i++)
1419                         {
1420                             index1[i] = index[i];
1421                         }
1422                         strcpy(gname1, gname);
1423                         if (parse_entry(&string, natoms, atoms, block, gn, &nr2, index2, gname2))
1424                         {
1425                             if (bOr)
1426                             {
1427                                 or_groups(nr1, index1, nr2, index2, &nr, index);
1428                                 sprintf(gname, "%s_%s", gname1, gname2);
1429                             }
1430                             else
1431                             {
1432                                 and_groups(nr1, index1, nr2, index2, &nr, index);
1433                                 sprintf(gname, "%s_&_%s", gname1, gname2);
1434                             }
1435                         }
1436                     }
1437                 }
1438                 while (bAnd || bOr);
1439             }
1440             while (string[0] == ' ')
1441             {
1442                 string++;
1443             }
1444             if (string[0])
1445             {
1446                 printf("\nSyntax error: \"%s\"\n", string);
1447             }
1448             else if (nr > 0)
1449             {
1450                 copy2block(nr, index, block);
1451                 srenew(*gn, block->nr);
1452                 newgroup        = block->nr-1;
1453                 (*gn)[newgroup] = gmx_strdup(gname);
1454             }
1455             else
1456             {
1457                 printf("Group is empty\n");
1458             }
1459         }
1460     }
1461     while (string[0] != 'q');
1462
1463     sfree(index);
1464     sfree(index1);
1465     sfree(index2);
1466 }
1467
1468 static int block2natoms(t_blocka *block)
1469 {
1470     int i, natoms;
1471
1472     natoms = 0;
1473     for (i = 0; i < block->nra; i++)
1474     {
1475         natoms = max(natoms, block->a[i]);
1476     }
1477     natoms++;
1478
1479     return natoms;
1480 }
1481
1482 void merge_blocks(t_blocka *dest, t_blocka *source)
1483 {
1484     int     i, nra0, i0;
1485
1486     /* count groups, srenew and fill */
1487     i0        = dest->nr;
1488     nra0      = dest->nra;
1489     dest->nr += source->nr;
1490     srenew(dest->index, dest->nr+1);
1491     for (i = 0; i < source->nr; i++)
1492     {
1493         dest->index[i0+i] = nra0 + source->index[i];
1494     }
1495     /* count atoms, srenew and fill */
1496     dest->nra += source->nra;
1497     srenew(dest->a, dest->nra);
1498     for (i = 0; i < source->nra; i++)
1499     {
1500         dest->a[nra0+i] = source->a[i];
1501     }
1502
1503     /* terminate list */
1504     dest->index[dest->nr] = dest->nra;
1505
1506 }
1507
1508 int gmx_make_ndx(int argc, char *argv[])
1509 {
1510     const char     *desc[] = {
1511         "Index groups are necessary for almost every GROMACS program.",
1512         "All these programs can generate default index groups. You ONLY",
1513         "have to use [THISMODULE] when you need SPECIAL index groups.",
1514         "There is a default index group for the whole system, 9 default",
1515         "index groups for proteins, and a default index group",
1516         "is generated for every other residue name.[PAR]",
1517         "When no index file is supplied, also [THISMODULE] will generate the",
1518         "default groups.",
1519         "With the index editor you can select on atom, residue and chain names",
1520         "and numbers.",
1521         "When a run input file is supplied you can also select on atom type.",
1522         "You can use NOT, AND and OR, you can split groups",
1523         "into chains, residues or atoms. You can delete and rename groups.[PAR]",
1524         "The atom numbering in the editor and the index file starts at 1.[PAR]",
1525         "The [TT]-twin[tt] switch duplicates all index groups with an offset of",
1526         "[TT]-natoms[tt], which is useful for Computational Electrophysiology",
1527         "double-layer membrane setups."
1528     };
1529
1530     static int      natoms     = 0;
1531     static gmx_bool bVerbose   = FALSE;
1532     static gmx_bool bDuplicate = FALSE;
1533     t_pargs         pa[]       = {
1534         { "-natoms",  FALSE, etINT, {&natoms},
1535           "set number of atoms (default: read from coordinate or index file)" },
1536         { "-twin",     FALSE, etBOOL, {&bDuplicate},
1537           "Duplicate all index groups with an offset of -natoms" },
1538         { "-verbose", FALSE, etBOOL, {&bVerbose},
1539           "HIDDENVerbose output" }
1540     };
1541 #define NPA asize(pa)
1542
1543     output_env_t oenv;
1544     char         title[STRLEN];
1545     int          nndxin;
1546     const char  *stxfile;
1547     char       **ndxinfiles;
1548     const char  *ndxoutfile;
1549     gmx_bool     bNatoms;
1550     int          i, j;
1551     t_atoms     *atoms;
1552     rvec        *x, *v;
1553     int          ePBC;
1554     matrix       box;
1555     t_blocka    *block, *block2;
1556     char       **gnames, **gnames2;
1557     t_filenm     fnm[] = {
1558         { efSTX, "-f", NULL,     ffOPTRD  },
1559         { efNDX, "-n", NULL,     ffOPTRDMULT },
1560         { efNDX, "-o", NULL,     ffWRITE }
1561     };
1562 #define NFILE asize(fnm)
1563
1564     if (!parse_common_args(&argc, argv, 0, NFILE, fnm, NPA, pa, asize(desc), desc,
1565                            0, NULL, &oenv))
1566     {
1567         return 0;
1568     }
1569
1570     stxfile = ftp2fn_null(efSTX, NFILE, fnm);
1571     if (opt2bSet("-n", NFILE, fnm))
1572     {
1573         nndxin = opt2fns(&ndxinfiles, "-n", NFILE, fnm);
1574     }
1575     else
1576     {
1577         nndxin = 0;
1578     }
1579     ndxoutfile = opt2fn("-o", NFILE, fnm);
1580     bNatoms    = opt2parg_bSet("-natoms", NPA, pa);
1581
1582     if (!stxfile && !nndxin)
1583     {
1584         gmx_fatal(FARGS, "No input files (structure or index)");
1585     }
1586
1587     if (stxfile)
1588     {
1589         snew(atoms, 1);
1590         get_stx_coordnum(stxfile, &(atoms->nr));
1591         init_t_atoms(atoms, atoms->nr, TRUE);
1592         snew(x, atoms->nr);
1593         snew(v, atoms->nr);
1594         fprintf(stderr, "\nReading structure file\n");
1595         read_stx_conf(stxfile, title, atoms, x, v, &ePBC, box);
1596         natoms  = atoms->nr;
1597         bNatoms = TRUE;
1598     }
1599     else
1600     {
1601         atoms = NULL;
1602         x     = NULL;
1603     }
1604
1605     /* read input file(s) */
1606     block  = new_blocka();
1607     gnames = NULL;
1608     printf("Going to read %d old index file(s)\n", nndxin);
1609     if (nndxin)
1610     {
1611         for (i = 0; i < nndxin; i++)
1612         {
1613             block2 = init_index(ndxinfiles[i], &gnames2);
1614             srenew(gnames, block->nr+block2->nr);
1615             for (j = 0; j < block2->nr; j++)
1616             {
1617                 gnames[block->nr+j] = gnames2[j];
1618             }
1619             sfree(gnames2);
1620             merge_blocks(block, block2);
1621             sfree(block2->a);
1622             sfree(block2->index);
1623 /*       done_block(block2); */
1624             sfree(block2);
1625         }
1626     }
1627     else
1628     {
1629         snew(gnames, 1);
1630         analyse(atoms, block, &gnames, FALSE, TRUE);
1631     }
1632
1633     if (!bNatoms)
1634     {
1635         natoms = block2natoms(block);
1636         printf("Counted atom numbers up to %d in index file\n", natoms);
1637     }
1638
1639     edit_index(natoms, atoms, x, block, &gnames, bVerbose);
1640
1641     write_index(ndxoutfile, block, gnames, bDuplicate, natoms);
1642
1643     return 0;
1644 }