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