src/tools/make_ndx.c

   1 /*
   2  *
   3  *                This source code is part of
   4  *
   5  *                 G   R   O   M   A   C   S
   6  *
   7  *          GROningen MAchine for Chemical Simulations
   8  *
   9  *                        VERSION 3.2.0
  10  * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
  11  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  12  * Copyright (c) 2001-2004, The GROMACS development team,
  13  * check out http://www.gromacs.org for more information.
  14
  15  * This program is free software; you can redistribute it and/or
  16  * modify it under the terms of the GNU General Public License
  17  * as published by the Free Software Foundation; either version 2
  18  * of the License, or (at your option) any later version.
  19  *
  20  * If you want to redistribute modifications, please consider that
  21  * scientific software is very special. Version control is crucial -
  22  * bugs must be traceable. We will be happy to consider code for
  23  * inclusion in the official distribution, but derived work must not
  24  * be called official GROMACS. Details are found in the README & COPYING
  25  * files - if they are missing, get the official version at www.gromacs.org.
  26  *
  27  * To help us fund GROMACS development, we humbly ask that you cite
  28  * the papers on the package - you can find them in the top README file.
  29  *
  30  * For more info, check our website at http://www.gromacs.org
  31  *
  32  * And Hey:
  33  * Green Red Orange Magenta Azure Cyan Skyblue
  34  */
  35 #ifdef HAVE_CONFIG_H
  36 #include <config.h>
  37 #endif
  38
  39 #include <ctype.h>
  40 #include "sysstuff.h"
  41 #include "strdb.h"
  42 #include "futil.h"
  43 #include "macros.h"
  44 #include "string2.h"
  45 #include "statutil.h"
  46 #include "confio.h"
  47 #include "copyrite.h"
  48 #include "typedefs.h"
  49 #include "index.h"
  50 #include "smalloc.h"
  51 #include "vec.h"
  52 #include "index.h"
  53
  54 #define MAXNAMES 30
  55 #define NAME_LEN 30
  56
  57 bool bCase=FALSE;
  58
  59 static int or_groups(atom_id nr1,atom_id *at1,atom_id nr2,atom_id *at2,
  60                      atom_id *nr,atom_id *at)
  61 {
  62   atom_id i1,i2,max=0;
  63   bool bNotIncr;
  64
  65   *nr=0;
  66
  67   bNotIncr=FALSE;
  68   for(i1=0; i1<nr1; i1++) {
  69     if ((i1>0) && (at1[i1] <= max))
  70       bNotIncr=TRUE;
  71     max=at1[i1];
  72   }
  73   for(i1=0; i1<nr2; i1++) {
  74     if ((i1>0) && (at2[i1] <= max))
  75       bNotIncr=TRUE;
  76     max=at2[i1];
  77   }
  78
  79   if (bNotIncr)
  80     printf("One of your groups is not ascending\n");
  81   else {
  82     i1=0;
  83     i2=0;
  84     *nr=0;
  85     while ((i1 < nr1) || (i2 < nr2)) {
  86       if ((i2 == nr2) || ((i1<nr1) && (at1[i1] < at2[i2]))) {
  87         at[*nr]=at1[i1];
  88         (*nr)++;
  89         i1++;
  90       }
  91       else {
  92         if ((i2<nr2) && ((i1==nr1) || (at1[i1] > at2[i2]))) {
  93           at[*nr]=at2[i2];
  94           (*nr)++;
  95         }
  96         i2++;
  97       }
  98     }
  99
 100     printf("Merged two groups with OR: %u %u -> %u\n",nr1,nr2,*nr);
 101   }
 102
 103   return *nr;
 104 }
 105
 106 static int and_groups(atom_id nr1,atom_id *at1,atom_id nr2,atom_id *at2,
 107                       atom_id *nr,atom_id *at)
 108 {
 109   atom_id i1,i2;
 110
 111   *nr=0;
 112   for (i1=0; i1<nr1; i1++) {
 113     for (i2=0; i2<nr2; i2++) {
 114       if (at1[i1]==at2[i2]) {
 115         at[*nr]=at1[i1];
 116         (*nr)++;
 117       }
 118     }
 119   }
 120
 121   printf("Merged two groups with AND: %u %u -> %u\n",nr1,nr2,*nr);
 122
 123   return *nr;
 124 }
 125
 126 static bool is_name_char(char c)
 127 {
 128   /* This string should contain all characters that can not be
 129    * the first letter of a name due to the make_ndx syntax.
 130    */
 131   const char *spec=" !&|";
 132
 133   return (c != '\0' && strchr(spec,c) == NULL);
 134 }
 135
 136 static int parse_names(char **string,int *n_names,char **names)
 137 {
 138   int i;
 139
 140   *n_names=0;
 141   while (is_name_char((*string)[0]) || (*string)[0] == ' ') {
 142     if (is_name_char((*string)[0])) {
 143       if (*n_names >= MAXNAMES)
 144         gmx_fatal(FARGS,"To many names: %d\n",*n_names+1);
 145       i=0;
 146       while (is_name_char((*string)[i])) {
 147         names[*n_names][i]=(*string)[i];
 148         i++;
 149         if (i > NAME_LEN) {
 150           printf("Name is too long, the maximum is %d characters\n",NAME_LEN);
 151           return 0;
 152         }
 153       }
 154       names[*n_names][i]='\0';
 155       if (!bCase)
 156         upstring(names[*n_names]);
 157       *string += i;
 158       (*n_names)++;
 159     }
 160     else
 161       (*string)++;
 162   }
 163
 164   return *n_names;
 165 }
 166
 167 static bool parse_int_char(char **string,int *nr,char *c)
 168 {
 169   char *orig;
 170   bool bRet;
 171
 172   orig = *string;
 173
 174   while ((*string)[0]==' ')
 175     (*string)++;
 176
 177   bRet=FALSE;
 178
 179   *c = ' ';
 180
 181   if (isdigit((*string)[0])) {
 182     *nr=(*string)[0]-'0';
 183     (*string)++;
 184     while (isdigit((*string)[0])) {
 185       *nr = (*nr)*10+(*string)[0]-'0';
 186       (*string)++;
 187     }
 188     if (isalpha((*string)[0])) {
 189       *c = (*string)[0];
 190       (*string)++;
 191     }
 192     /* Check if there is at most one non-digit character */
 193     if (!isalnum((*string)[0])) {
 194       bRet = TRUE;
 195     } else {
 196       *string = orig;
 197     }
 198   }
 199   else
 200     *nr=NOTSET;
 201
 202   return bRet;
 203 }
 204
 205 static bool parse_int(char **string,int *nr)
 206 {
 207   char *orig,c;
 208   bool bRet;
 209
 210   orig = *string;
 211   bRet = parse_int_char(string,nr,&c);
 212   if (bRet && c != ' ') {
 213     *string = orig;
 214     bRet = FALSE;
 215   }
 216
 217   return bRet;
 218 }
 219
 220 static bool isquote(char c)
 221 {
 222   return (c == '\"');
 223 }
 224
 225 static bool parse_string(char **string,int *nr, int ngrps, char **grpname)
 226 {
 227   char *s, *sp;
 228   char c;
 229
 230   while ((*string)[0]==' ')
 231     (*string)++;
 232
 233   (*nr) = NOTSET;
 234   if (isquote((*string)[0])) {
 235     c=(*string)[0];
 236     (*string)++;
 237     s = strdup((*string));
 238     sp = strchr(s, c);
 239     if (sp!=NULL) {
 240       (*string) += sp-s + 1;
 241       sp[0]='\0';
 242       (*nr) = find_group(s, ngrps, grpname);
 243     }
 244   }
 245
 246   return (*nr) != NOTSET;
 247 }
 248
 249 static int select_atomnumbers(char **string,t_atoms *atoms,atom_id n1,
 250                               atom_id *nr,atom_id *index,char *gname)
 251 {
 252   char    buf[STRLEN];
 253   int     i,up;
 254
 255   *nr=0;
 256   while ((*string)[0]==' ')
 257     (*string)++;
 258   if ((*string)[0]=='-') {
 259     (*string)++;
 260     parse_int(string,&up);
 261     if ((n1<1) || (n1>atoms->nr) || (up<1) || (up>atoms->nr))
 262       printf("Invalid atom range\n");
 263     else {
 264       for(i=n1-1; i<=up-1; i++) {
 265         index[*nr]=i;
 266         (*nr)++;
 267       }
 268       printf("Found %u atom%s in range %u-%d\n",*nr,(*nr==1)?"":"s",n1,up);
 269       if (n1==up)
 270         sprintf(buf,"a_%u",n1);
 271       else
 272         sprintf(buf,"a_%u-%d",n1,up);
 273       strcpy(gname,buf);
 274     }
 275   }
 276   else {
 277     i=n1;
 278     sprintf(gname,"a");
 279     do {
 280       if ((i-1>=0) && (i-1<atoms->nr)) {
 281         index[*nr] = i-1;
 282         (*nr)++;
 283         sprintf(buf,"_%d",i);
 284         strcat(gname,buf);
 285       } else {
 286         printf("Invalid atom number %d\n",i);
 287         *nr = 0;
 288       }
 289     } while ((*nr!=0) && (parse_int(string,&i)));
 290   }
 291
 292   return *nr;
 293 }
 294
 295 static int select_residuenumbers(char **string,t_atoms *atoms,
 296                                  atom_id n1,char c,
 297                                  atom_id *nr,atom_id *index,char *gname)
 298 {
 299   char    buf[STRLEN];
 300   int     i,j,up;
 301   t_resinfo *ri;
 302
 303   *nr=0;
 304   while ((*string)[0]==' ')
 305     (*string)++;
 306   if ((*string)[0]=='-') {
 307     /* Residue number range selection */
 308     if (c != ' ') {
 309       printf("Error: residue insertion codes can not be used with residue range selection\n");
 310       return 0;
 311     }
 312     (*string)++;
 313     parse_int(string,&up);
 314
 315     for(i=0; i<atoms->nr; i++) {
 316       ri = &atoms->resinfo[atoms->atom[i].resind];
 317       for(j=n1; (j<=up); j++) {
 318         if (ri->nr == j && (c == ' ' || ri->ic == c)) {
 319           index[*nr]=i;
 320           (*nr)++;
 321         }
 322       }
 323     }
 324     printf("Found %u atom%s with res.nr. in range %u-%d\n",
 325            *nr,(*nr==1)?"":"s",n1,up);
 326     if (n1==up)
 327       sprintf(buf,"r_%u",n1);
 328     else
 329       sprintf(buf,"r_%u-%d",n1,up);
 330     strcpy(gname,buf);
 331   }
 332   else {
 333     /* Individual residue number/insertion code selection */
 334     j=n1;
 335     sprintf(gname,"r");
 336     do {
 337       for(i=0; i<atoms->nr; i++) {
 338         ri = &atoms->resinfo[atoms->atom[i].resind];
 339         if (ri->nr == j && ri->ic == c) {
 340         index[*nr]=i;
 341         (*nr)++;
 342         }
 343       }
 344       sprintf(buf,"_%d",j);
 345       strcat(gname,buf);
 346     } while (parse_int_char(string,&j,&c));
 347   }
 348
 349   return *nr;
 350 }
 351
 352 static int select_residueindices(char **string,t_atoms *atoms,
 353                                  atom_id n1,char c,
 354                                  atom_id *nr,atom_id *index,char *gname)
 355 {
 356   /*this should be similar to select_residuenumbers except select by index (sequential numbering in file)*/
 357   /*resind+1 for 1-indexing*/
 358   char    buf[STRLEN];
 359   int     i,j,up;
 360   t_resinfo *ri;
 361
 362   *nr=0;
 363   while ((*string)[0]==' ')
 364     (*string)++;
 365   if ((*string)[0]=='-') {
 366     /* Residue number range selection */
 367     if (c != ' ') {
 368       printf("Error: residue insertion codes can not be used with residue range selection\n");
 369       return 0;
 370     }
 371     (*string)++;
 372     parse_int(string,&up);
 373
 374     for(i=0; i<atoms->nr; i++) {
 375       ri = &atoms->resinfo[atoms->atom[i].resind];
 376       for(j=n1; (j<=up); j++) {
 377         if (atoms->atom[i].resind+1 == j && (c == ' ' || ri->ic == c)) {
 378           index[*nr]=i;
 379           (*nr)++;
 380         }
 381       }
 382     }
 383     printf("Found %u atom%s with resind.+1 in range %u-%d\n",
 384            *nr,(*nr==1)?"":"s",n1,up);
 385     if (n1==up)
 386       sprintf(buf,"r_%u",n1);
 387     else
 388       sprintf(buf,"r_%u-%d",n1,up);
 389     strcpy(gname,buf);
 390   }
 391   else {
 392     /* Individual residue number/insertion code selection */
 393     j=n1;
 394     sprintf(gname,"r");
 395     do {
 396       for(i=0; i<atoms->nr; i++) {
 397         ri = &atoms->resinfo[atoms->atom[i].resind];
 398         if (atoms->atom[i].resind+1 == j && ri->ic == c) {
 399         index[*nr]=i;
 400         (*nr)++;
 401         }
 402       }
 403       sprintf(buf,"_%d",j);
 404       strcat(gname,buf);
 405     } while (parse_int_char(string,&j,&c));
 406   }
 407
 408   return *nr;
 409 }
 410
 411
 412 static bool atoms_from_residuenumbers(t_atoms *atoms,int group,t_blocka *block,
 413                                       atom_id *nr,atom_id *index,char *gname)
 414 {
 415   int i,j,j0,j1,resnr,nres;
 416
 417   j0=block->index[group];
 418   j1=block->index[group+1];
 419   nres = atoms->nres;
 420   for(j=j0; j<j1; j++)
 421     if (block->a[j]>=nres) {
 422       printf("Index %s contains number>nres (%d>%d)\n",
 423              gname,block->a[j]+1,nres);
 424       return FALSE;
 425     }
 426   for(i=0; i<atoms->nr; i++) {
 427     resnr = atoms->resinfo[atoms->atom[i].resind].nr;
 428     for (j=j0; j<j1; j++)
 429       if (block->a[j]+1 == resnr) {
 430         index[*nr]=i;
 431         (*nr)++;
 432         break;
 433       }
 434   }
 435   printf("Found %u atom%s in %d residues from group %s\n",
 436          *nr,(*nr==1)?"":"s",j1-j0,gname);
 437   return *nr;
 438 }
 439
 440 static bool comp_name(char *name,char *search)
 441 {
 442   while (name[0] != '\0' && search[0] != '\0') {
 443     switch (search[0]) {
 444     case '?':
 445       /* Always matches */
 446       break;
 447     case '*':
 448       if (search[1] != '\0') {
 449         printf("WARNING: Currently '*' is only supported at the end of an expression\n");
 450         return FALSE;
 451       } else {
 452         return TRUE;
 453       }
 454       break;
 455     default:
 456       /* Compare a single character */
 457       if (( bCase && strncmp(name,search,1)) ||
 458           (!bCase && gmx_strncasecmp(name,search,1))) {
 459         return FALSE;
 460       }
 461     }
 462     name++;
 463     search++;
 464   }
 465
 466   return (name[0] == '\0' && (search[0] == '\0' || search[0] == '*'));
 467 }
 468
 469 static int select_chainnames(t_atoms *atoms,int n_names,char **names,
 470                              atom_id *nr,atom_id *index)
 471 {
 472   char name[2];
 473   int j;
 474   atom_id i;
 475
 476   name[1]=0;
 477   *nr=0;
 478   for(i=0; i<atoms->nr; i++) {
 479     name[0] = atoms->resinfo[atoms->atom[i].resind].chainid;
 480     j=0;
 481     while (j<n_names && !comp_name(name,names[j]))
 482       j++;
 483     if (j<n_names) {
 484       index[*nr]=i;
 485       (*nr)++;
 486     }
 487   }
 488   printf("Found %u atom%s with chain identifier%s",
 489          *nr,(*nr==1)?"":"s",(n_names==1)?"":"s");
 490   for(j=0; (j<n_names); j++)
 491     printf(" %s",names[j]);
 492   printf("\n");
 493
 494   return *nr;
 495 }
 496
 497 static int select_atomnames(t_atoms *atoms,int n_names,char **names,
 498                             atom_id *nr,atom_id *index,bool bType)
 499 {
 500   char *name;
 501   int j;
 502   atom_id i;
 503
 504   *nr=0;
 505   for(i=0; i<atoms->nr; i++) {
 506     if (bType)
 507       name=*(atoms->atomtype[i]);
 508     else
 509       name=*(atoms->atomname[i]);
 510     j=0;
 511     while (j<n_names && !comp_name(name,names[j]))
 512       j++;
 513     if (j<n_names) {
 514       index[*nr]=i;
 515       (*nr)++;
 516     }
 517   }
 518   printf("Found %u atoms with %s%s",
 519          *nr,bType ? "type" : "name",(n_names==1)?"":"s");
 520   for(j=0; (j<n_names); j++)
 521     printf(" %s",names[j]);
 522   printf("\n");
 523
 524   return *nr;
 525 }
 526
 527 static int select_residuenames(t_atoms *atoms,int n_names,char **names,
 528                                atom_id *nr,atom_id *index)
 529 {
 530   char *name;
 531   int j;
 532   atom_id i;
 533
 534   *nr=0;
 535   for(i=0; i<atoms->nr; i++) {
 536     name = *(atoms->resinfo[atoms->atom[i].resind].name);
 537     j=0;
 538     while (j<n_names && !comp_name(name,names[j]))
 539       j++;
 540     if (j<n_names) {
 541       index[*nr]=i;
 542       (*nr)++;
 543     }
 544   }
 545   printf("Found %u atoms with residue name%s",*nr,(n_names==1)?"":"s");
 546   for(j=0; (j<n_names); j++)
 547     printf(" %s",names[j]);
 548   printf("\n");
 549
 550   return *nr;
 551 }
 552
 553 static void copy2block(int n,atom_id *index,t_blocka *block)
 554 {
 555   int i,n0;
 556
 557   block->nr++;
 558   n0=block->nra;
 559   block->nra=n0+n;
 560   srenew(block->index,block->nr+1);
 561   block->index[block->nr]=n0+n;
 562   srenew(block->a,n0+n);
 563   for(i=0; (i<n); i++)
 564     block->a[n0+i]=index[i];
 565 }
 566
 567 static void make_gname(int n,char **names,char *gname)
 568 {
 569   int i;
 570
 571   strcpy(gname,names[0]);
 572   for (i=1; i<n; i++) {
 573     strcat(gname,"_");
 574     strcat(gname,names[i]);
 575   }
 576 }
 577
 578 static void copy_group(int g,t_blocka *block,atom_id *nr,atom_id *index)
 579 {
 580   int i,i0;
 581
 582   i0=block->index[g];
 583   *nr=block->index[g+1]-i0;
 584   for (i=0; i<*nr; i++)
 585     index[i]=block->a[i0+i];
 586 }
 587
 588 static void remove_group(int nr,int nr2,t_blocka *block,char ***gn)
 589 {
 590   int i,j,shift;
 591   char *name;
 592
 593   if (nr2==NOTSET)
 594     nr2=nr;
 595
 596   for(j=0; j<=nr2-nr; j++) {
 597     if ((nr<0) || (nr>=block->nr))
 598       printf("Group %d does not exist\n",nr+j);
 599     else {
 600       shift=block->index[nr+1]-block->index[nr];
 601       for(i=block->index[nr+1]; i<block->nra; i++)
 602         block->a[i-shift]=block->a[i];
 603
 604       for(i=nr; i<block->nr; i++)
 605         block->index[i]=block->index[i+1]-shift;
 606       name = strdup((*gn)[nr]);
 607       sfree((*gn)[nr]);
 608       for(i=nr; i<block->nr-1; i++) {
 609         (*gn)[i]=(*gn)[i+1];
 610       }
 611       block->nr--;
 612       block->nra=block->index[block->nr];
 613       printf("Removed group %d '%s'\n",nr+j,name);
 614       sfree(name);
 615     }
 616   }
 617 }
 618
 619 static void split_group(t_atoms *atoms,int sel_nr,t_blocka *block,char ***gn,
 620                         bool bAtom)
 621 {
 622   char buf[STRLEN],*name;
 623   int i,resind;
 624   atom_id a,n0,n1;
 625
 626   printf("Splitting group %d '%s' into %s\n",sel_nr,(*gn)[sel_nr],
 627          bAtom ? "atoms" : "residues");
 628
 629   n0=block->index[sel_nr];
 630   n1=block->index[sel_nr+1];
 631   srenew(block->a,block->nra+n1-n0);
 632   for (i=n0; i<n1; i++) {
 633     a=block->a[i];
 634     resind = atoms->atom[a].resind;
 635     name = *(atoms->resinfo[resind].name);
 636     if (bAtom || (i==n0) || (atoms->atom[block->a[i-1]].resind!=resind)) {
 637       if (i>n0)
 638         block->index[block->nr]=block->nra;
 639       block->nr++;
 640       srenew(block->index,block->nr+1);
 641       srenew(*gn,block->nr);
 642       if (bAtom)
 643         sprintf(buf,"%s_%s_%u",(*gn)[sel_nr],*atoms->atomname[a],a+1);
 644       else
 645         sprintf(buf,"%s_%s_%d",(*gn)[sel_nr],name,atoms->resinfo[resind].nr);
 646       (*gn)[block->nr-1]=strdup(buf);
 647     }
 648     block->a[block->nra]=a;
 649     block->nra++;
 650   }
 651   block->index[block->nr]=block->nra;
 652 }
 653
 654 static int split_chain(t_atoms *atoms,rvec *x,
 655                         int sel_nr,t_blocka *block,char ***gn)
 656 {
 657   char    buf[STRLEN];
 658   int     j,nchain;
 659   atom_id i,a,natoms,*start=NULL,*end=NULL,ca_start,ca_end;
 660   rvec    vec;
 661
 662   natoms=atoms->nr;
 663   nchain=0;
 664   ca_start=0;
 665
 666   while (ca_start<natoms) {
 667     while((ca_start<natoms) && strcmp(*atoms->atomname[ca_start],"CA"))
 668       ca_start++;
 669     if (ca_start<natoms) {
 670       srenew(start,nchain+1);
 671       srenew(end,nchain+1);
 672       start[nchain]=ca_start;
 673       while ((start[nchain]>0) &&
 674              (atoms->atom[start[nchain]-1].resind ==
 675               atoms->atom[ca_start].resind))
 676         start[nchain]--;
 677
 678       i=ca_start;
 679       do {
 680         ca_end=i;
 681         do {
 682           i++;
 683         } while ((i<natoms) && strcmp(*atoms->atomname[i],"CA"));
 684         if (i<natoms)
 685           rvec_sub(x[ca_end],x[i],vec);
 686       } while ((i<natoms) && (norm(vec)<0.45));
 687
 688       end[nchain]=ca_end;
 689       while ((end[nchain]+1<natoms) &&
 690             (atoms->atom[end[nchain]+1].resind==atoms->atom[ca_end].resind))
 691         end[nchain]++;
 692       ca_start=end[nchain]+1;
 693       nchain++;
 694     }
 695   }
 696   if (nchain==1)
 697     printf("Found 1 chain, will not split\n");
 698   else
 699     printf("Found %d chains\n",nchain);
 700   for (j=0; j<nchain; j++)
 701     printf("%d:%6u atoms (%u to %u)\n",
 702            j+1,end[j]-start[j]+1,start[j]+1,end[j]+1);
 703
 704   if (nchain > 1) {
 705     srenew(block->a,block->nra+block->index[sel_nr+1]-block->index[sel_nr]);
 706     for (j=0; j<nchain; j++) {
 707       block->nr++;
 708       srenew(block->index,block->nr+1);
 709       srenew(*gn,block->nr);
 710       sprintf(buf,"%s_chain%d",(*gn)[sel_nr],j+1);
 711       (*gn)[block->nr-1]=strdup(buf);
 712       for (i=block->index[sel_nr]; i<block->index[sel_nr+1]; i++) {
 713         a=block->a[i];
 714         if ((a>=start[j]) && (a<=end[j])) {
 715           block->a[block->nra]=a;
 716           block->nra++;
 717         }
 718       }
 719       block->index[block->nr]=block->nra;
 720       if (block->index[block->nr-1]==block->index[block->nr])
 721         remove_group(block->nr-1,NOTSET,block,gn);
 722     }
 723   }
 724   sfree(start);
 725   sfree(end);
 726
 727   return nchain;
 728 }
 729
 730 static bool check_have_atoms(t_atoms *atoms, char *string)
 731 {
 732   if ( atoms==NULL ) {
 733     printf("Can not process '%s' without atoms info\n", string);
 734     return FALSE;
 735   } else
 736     return TRUE;
 737 }
 738
 739 static bool parse_entry(char **string,int natoms,t_atoms *atoms,
 740                         t_blocka *block,char ***gn,
 741                         atom_id *nr,atom_id *index,char *gname)
 742 {
 743   static char **names, *ostring;
 744   static bool bFirst=TRUE;
 745   int         j,n_names,sel_nr1;
 746   atom_id     i,nr1,*index1;
 747   char        c;
 748   bool        bRet,bCompl;
 749
 750   if (bFirst) {
 751     bFirst=FALSE;
 752     snew(names,MAXNAMES);
 753     for (i=0; i<MAXNAMES; i++)
 754       snew(names[i],NAME_LEN+1);
 755   }
 756
 757   bRet=FALSE;
 758   sel_nr1=NOTSET;
 759
 760   while(*string[0]==' ')
 761     (*string)++;
 762
 763   if ((*string)[0]=='!') {
 764     bCompl=TRUE;
 765     (*string)++;
 766      while(*string[0]==' ')
 767        (*string)++;
 768   } else
 769     bCompl=FALSE;
 770
 771   ostring = *string;
 772
 773   if (parse_int(string,&sel_nr1) ||
 774       parse_string(string,&sel_nr1,block->nr,*gn)) {
 775     if ((sel_nr1>=0) && (sel_nr1<block->nr)) {
 776       copy_group(sel_nr1,block,nr,index);
 777       strcpy(gname,(*gn)[sel_nr1]);
 778       printf("Copied index group %d '%s'\n",sel_nr1,(*gn)[sel_nr1]);
 779       bRet=TRUE;
 780     } else
 781       printf("Group %d does not exist\n",sel_nr1);
 782   }
 783   else if ((*string)[0]=='a') {
 784     (*string)++;
 785     if (check_have_atoms(atoms, ostring)) {
 786     if (parse_int(string,&sel_nr1)) {
 787       bRet=select_atomnumbers(string,atoms,sel_nr1,nr,index,gname);
 788     }
 789     else if (parse_names(string,&n_names,names)) {
 790       bRet=select_atomnames(atoms,n_names,names,nr,index,FALSE);
 791       make_gname(n_names,names,gname);
 792       }
 793     }
 794   }
 795   else if ((*string)[0]=='t') {
 796     (*string)++;
 797     if (check_have_atoms(atoms, ostring) &&
 798         parse_names(string,&n_names,names)) {
 799       if (atoms->atomtype == NULL)
 800         printf("Need a run input file to select atom types\n");
 801       else {
 802         bRet=select_atomnames(atoms,n_names,names,nr,index,TRUE);
 803         make_gname(n_names,names,gname);
 804       }
 805     }
 806   }
 807   else if (strncmp(*string,"res",3)==0) {
 808     (*string)+=3;
 809     if ( check_have_atoms(atoms, ostring) &&
 810          parse_int(string,&sel_nr1) &&
 811          (sel_nr1>=0) && (sel_nr1<block->nr) ) {
 812       bRet=atoms_from_residuenumbers(atoms,
 813                                      sel_nr1,block,nr,index,(*gn)[sel_nr1]);
 814       sprintf(gname,"atom_%s",(*gn)[sel_nr1]);
 815     }
 816   }
 817   else if (strncmp(*string,"ri",2)==0) {
 818     (*string)+=2;
 819     if (check_have_atoms(atoms, ostring) &&
 820                 parse_int_char(string,&sel_nr1,&c)) {
 821                 bRet=select_residueindices(string,atoms,sel_nr1,c,nr,index,gname);
 822     }
 823   }
 824   else if ((*string)[0]=='r') {
 825     (*string)++;
 826     if (check_have_atoms(atoms, ostring)) {
 827       if (parse_int_char(string,&sel_nr1,&c)) {
 828         bRet=select_residuenumbers(string,atoms,sel_nr1,c,nr,index,gname);
 829     }
 830     else if (parse_names(string,&n_names,names)) {
 831       bRet=select_residuenames(atoms,n_names,names,nr,index);
 832       make_gname(n_names,names,gname);
 833     }
 834   }
 835   }
 836   else if (strncmp(*string,"chain",5)==0) {
 837     (*string)+=5;
 838     if (check_have_atoms(atoms, ostring) &&
 839         parse_names(string,&n_names,names)) {
 840       bRet=select_chainnames(atoms,n_names,names,nr,index);
 841       sprintf(gname,"ch%s",names[0]);
 842       for (i=1; i<n_names; i++)
 843         strcat(gname,names[i]);
 844     }
 845   }
 846   if (bRet && bCompl) {
 847     snew(index1,natoms-*nr);
 848     nr1=0;
 849     for(i=0; i<natoms; i++) {
 850       j=0;
 851       while ((j<*nr) && (index[j] != i))
 852         j++;
 853       if (j==*nr) {
 854         if (nr1 >= natoms-*nr) {
 855           printf("There are double atoms in your index group\n");
 856           break;
 857         }
 858         index1[nr1]=i;
 859         nr1++;
 860       }
 861     }
 862     *nr=nr1;
 863     for(i=0; i<nr1; i++)
 864       index[i]=index1[i];
 865     sfree(index1);
 866
 867     for (i=strlen(gname)+1; i>0; i--)
 868       gname[i]=gname[i-1];
 869     gname[0]='!';
 870     printf("Complemented group: %u atoms\n",*nr);
 871   }
 872
 873   return bRet;
 874 }
 875
 876 static void list_residues(t_atoms *atoms)
 877 {
 878   int i,j,start,end,prev_resind,resind;
 879   bool bDiff;
 880
 881   /* Print all the residues, assuming continuous resnr count */
 882   start = atoms->atom[0].resind;
 883   prev_resind = start;
 884   for(i=0; i<atoms->nr; i++) {
 885     resind = atoms->atom[i].resind;
 886     if ((resind != prev_resind) || (i==atoms->nr-1)) {
 887       if ((bDiff=strcmp(*atoms->resinfo[resind].name,
 888                         *atoms->resinfo[start].name)) ||
 889           (i==atoms->nr-1)) {
 890         if (bDiff)
 891           end = prev_resind;
 892         else
 893           end = resind;
 894         if (end < start+3)
 895           for(j=start; j<=end; j++)
 896             printf("%4d %-5s",
 897                    j+1,*(atoms->resinfo[j].name));
 898         else
 899           printf(" %4d - %4d %-5s  ",
 900                  start+1,end+1,*(atoms->resinfo[start].name));
 901         start = resind;
 902       }
 903     }
 904     prev_resind = resind;
 905   }
 906   printf("\n");
 907 }
 908
 909 static void edit_index(int natoms, t_atoms *atoms,rvec *x,t_blocka *block, char ***gn, bool bVerbose)
 910 {
 911   static char **atnames, *ostring;
 912   static bool bFirst=TRUE;
 913   char inp_string[STRLEN],*string;
 914   char gname[STRLEN],gname1[STRLEN],gname2[STRLEN];
 915   int  i,i0,i1,sel_nr,sel_nr2,newgroup;
 916   atom_id nr,nr1,nr2,*index,*index1,*index2;
 917   bool bAnd,bOr,bPrintOnce;
 918
 919   if (bFirst) {
 920     bFirst=FALSE;
 921     snew(atnames,MAXNAMES);
 922     for (i=0; i<MAXNAMES; i++)
 923       snew(atnames[i],NAME_LEN+1);
 924   }
 925
 926   string=NULL;
 927
 928   snew(index,natoms);
 929   snew(index1,natoms);
 930   snew(index2,natoms);
 931
 932   newgroup=NOTSET;
 933   bPrintOnce = TRUE;
 934   do {
 935     gname1[0]='\0';
 936     if (bVerbose || bPrintOnce || newgroup!=NOTSET) {
 937       printf("\n");
 938       if (bVerbose || bPrintOnce || newgroup==NOTSET) {
 939         i0=0;
 940         i1=block->nr;
 941       } else {
 942         i0=newgroup;
 943         i1=newgroup+1;
 944       }
 945       for(i=i0; i<i1; i++)
 946       printf("%3d %-20s: %5u atoms\n",i,(*gn)[i],
 947              block->index[i+1]-block->index[i]);
 948       newgroup=NOTSET;
 949     }
 950     if (bVerbose || bPrintOnce) {
 951       printf("\n");
 952       printf(" nr : group       !   'name' nr name   'splitch' nr    Enter: list groups\n");
 953       printf(" 'a': atom        &   'del' nr         'splitres' nr   'l': list residues\n");
 954       printf(" 't': atom type   |   'keep' nr        'splitat' nr    'h': help\n");
 955       printf(" 'r': residue         'res' nr         'chain' char\n");
 956       printf(" \"name\": group        'case': case %s         'q': save and quit\n",
 957              bCase ? "insensitive" : "sensitive  ");
 958           printf(" 'ri': residue index\n");
 959       bPrintOnce = FALSE;
 960     }
 961     printf("\n");
 962     printf("> ");
 963     if(NULL==fgets(inp_string,STRLEN,stdin))
 964     {
 965         gmx_fatal(FARGS,"Error reading user input");
 966     }
 967     inp_string[strlen(inp_string)-1]=0;
 968     printf("\n");
 969     string=inp_string;
 970     while (string[0]==' ')
 971       string++;
 972
 973     ostring = string;
 974     nr=0;
 975     if (string[0] == 'h') {
 976       printf(" nr                : selects an index group by number or quoted string.\n");
 977       printf("                     The string is first matched against the whole group name,\n");
 978       printf("                     then against the beginning and finally against an\n");
 979       printf("                     arbitrary substring. A multiple match is an error.\n");
 980
 981       printf(" 'a' nr1 [nr2 ...] : selects atoms, atom numbering starts at 1.\n");
 982       printf(" 'a' nr1 - nr2     : selects atoms in the range from nr1 to nr2.\n");
 983       printf(" 'a' name1[*] [name2[*] ...] : selects atoms by name(s), '?' matches any char,\n");
 984       printf("                               wildcard '*' allowed at the end of a name.\n");
 985       printf(" 't' type1[*] [type2[*] ...] : as 'a', but for type, run input file required.\n");
 986       printf(" 'r' nr1[ic1] [nr2[ic2] ...] : selects residues by number and insertion code.\n");
 987       printf(" 'r' nr1 - nr2               : selects residues in the range from nr1 to nr2.\n");
 988       printf(" 'r' name1[*] [name2[*] ...] : as 'a', but for residue names.\n");
 989       printf(" 'ri' nr1 - nr2              : selects residue indices, 1-indexed, (as opposed to numbers) in the range from nr1 to nr2.\n");
 990       printf(" 'chain' ch1 [ch2 ...]       : selects atoms by chain identifier(s),\n");
 991       printf("                               not available with a .gro file as input.\n");
 992       printf(" !                 : takes the complement of a group with respect to all\n");
 993       printf("                     the atoms in the input file.\n");
 994       printf(" & |               : AND and OR, can be placed between any of the options\n");
 995       printf("                     above, the input is processed from left to right.\n");
 996       printf(" 'name' nr name    : rename group nr to name.\n");
 997       printf(" 'del' nr1 [- nr2] : deletes one group or groups in the range from nr1 to nr2.\n");
 998       printf(" 'keep' nr         : deletes all groups except nr.\n");
 999       printf(" 'case'            : make all name compares case (in)sensitive.\n");
1000       printf(" 'splitch' nr      : split group into chains using CA distances.\n");
1001       printf(" 'splitres' nr     : split group into residues.\n");
1002       printf(" 'splitat' nr      : split group into atoms.\n");
1003       printf(" 'res' nr          : interpret numbers in group as residue numbers\n");
1004       printf(" Enter             : list the currently defined groups and commands\n");
1005       printf(" 'l'               : list the residues.\n");
1006       printf(" 'h'               : show this help.\n");
1007       printf(" 'q'               : save and quit.\n");
1008       printf("\n");
1009       printf(" Examples:\n");
1010       printf(" > 2 | 4 & r 3-5\n");
1011       printf(" selects all atoms from group 2 and 4 that have residue numbers 3, 4 or 5\n");
1012       printf(" > a C* & !a C CA\n");
1013       printf(" selects all atoms starting with 'C' but not the atoms 'C' and 'CA'\n");
1014       printf(" > \"protein\" & ! \"backb\"\n");
1015       printf(" selects all atoms that are in group 'protein' and not in group 'backbone'\n");
1016       if (bVerbose) {
1017         printf("\npress Enter ");
1018         getchar();
1019       }
1020     }
1021     else if (strncmp(string,"del",3)==0) {
1022       string+=3;
1023       if (parse_int(&string,&sel_nr)) {
1024         while(string[0]==' ')
1025           string++;
1026         if (string[0]=='-') {
1027           string++;
1028           parse_int(&string,&sel_nr2);
1029         } else
1030           sel_nr2=NOTSET;
1031         while(string[0]==' ')
1032           string++;
1033         if (string[0]=='\0')
1034           remove_group(sel_nr,sel_nr2,block,gn);
1035         else
1036           printf("\nSyntax error: \"%s\"\n",string);
1037       }
1038     }
1039     else if (strncmp(string,"keep",4)==0) {
1040       string+=4;
1041       if (parse_int(&string,&sel_nr)) {
1042         remove_group(sel_nr+1,block->nr-1,block,gn);
1043         remove_group(0,sel_nr-1,block,gn);
1044       }
1045     }
1046     else if (strncmp(string,"name",4)==0) {
1047       string+=4;
1048       if (parse_int(&string,&sel_nr)) {
1049         if ((sel_nr>=0) && (sel_nr<block->nr)) {
1050           sscanf(string,"%s",gname);
1051           sfree((*gn)[sel_nr]);
1052           (*gn)[sel_nr]=strdup(gname);
1053         }
1054       }
1055     }
1056     else if (strncmp(string,"case",4)==0) {
1057       bCase=!bCase;
1058       printf("Switched to case %s\n",bCase ? "sensitive" : "insensitive");
1059     }
1060     else if (string[0] == 'v' ) {
1061       bVerbose=!bVerbose;
1062       printf("Turned verbose %s\n",bVerbose ? "on" : "off");
1063     }
1064     else if (string[0] == 'l') {
1065       if ( check_have_atoms(atoms, ostring) )
1066       list_residues(atoms);
1067     }
1068     else if (strncmp(string,"splitch",7)==0) {
1069       string+=7;
1070       if ( check_have_atoms(atoms, ostring) &&
1071            parse_int(&string,&sel_nr) &&
1072            (sel_nr>=0) && (sel_nr<block->nr))
1073           split_chain(atoms,x,sel_nr,block,gn);
1074     }
1075     else if (strncmp(string,"splitres",8)==0 ) {
1076       string+=8;
1077       if ( check_have_atoms(atoms, ostring) &&
1078            parse_int(&string,&sel_nr) &&
1079            (sel_nr>=0) && (sel_nr<block->nr))
1080           split_group(atoms,sel_nr,block,gn,FALSE);
1081     }
1082     else if (strncmp(string,"splitat",7)==0 ) {
1083      string+=7;
1084       if ( check_have_atoms(atoms, ostring) &&
1085            parse_int(&string,&sel_nr) &&
1086            (sel_nr>=0) && (sel_nr<block->nr))
1087           split_group(atoms,sel_nr,block,gn,TRUE);
1088     }
1089     else if (string[0] == '\0') {
1090       bPrintOnce = TRUE;
1091     }
1092     else if (string[0] != 'q') {
1093       nr1=-1;
1094       nr2=-1;
1095       if (parse_entry(&string,natoms,atoms,block,gn,&nr,index,gname)) {
1096         do {
1097           while (string[0]==' ')
1098             string++;
1099
1100           bAnd=FALSE;
1101           bOr=FALSE;
1102           if (string[0]=='&')
1103             bAnd=TRUE;
1104           else if (string[0]=='|')
1105             bOr=TRUE;
1106
1107           if (bAnd || bOr) {
1108             string++;
1109             nr1=nr;
1110             for(i=0; i<nr; i++)
1111               index1[i]=index[i];
1112             strcpy(gname1,gname);
1113             if (parse_entry(&string,natoms,atoms,block,gn,&nr2,index2,gname2)) {
1114               if (bOr) {
1115                 or_groups(nr1,index1,nr2,index2,&nr,index);
1116                 sprintf(gname,"%s_%s",gname1,gname2);
1117               }
1118               else {
1119                 and_groups(nr1,index1,nr2,index2,&nr,index);
1120                 sprintf(gname,"%s_&_%s",gname1,gname2);
1121               }
1122             }
1123           }
1124         } while (bAnd || bOr);
1125       }
1126       while(string[0]==' ')
1127         string++;
1128       if (string[0])
1129         printf("\nSyntax error: \"%s\"\n",string);
1130       else if (nr>0) {
1131         copy2block(nr,index,block);
1132         srenew(*gn,block->nr);
1133         newgroup = block->nr-1;
1134         (*gn)[newgroup]=strdup(gname);
1135       }
1136       else
1137         printf("Group is empty\n");
1138     }
1139   } while (string[0]!='q');
1140
1141   sfree(index);
1142   sfree(index1);
1143   sfree(index2);
1144 }
1145
1146 static int block2natoms(t_blocka *block)
1147 {
1148   int i, natoms;
1149
1150   natoms = 0;
1151   for(i=0; i<block->nra; i++)
1152     natoms = max(natoms, block->a[i]);
1153   natoms++;
1154
1155   return natoms;
1156 }
1157
1158 void merge_blocks(t_blocka *dest, t_blocka *source)
1159 {
1160   int     i,nra0,i0;
1161
1162   /* count groups, srenew and fill */
1163   i0 = dest->nr;
1164   nra0 = dest->nra;
1165   dest->nr+=source->nr;
1166   srenew(dest->index, dest->nr+1);
1167   for(i=0; i<source->nr; i++)
1168     dest->index[i0+i] = nra0 + source->index[i];
1169   /* count atoms, srenew and fill */
1170   dest->nra+=source->nra;
1171   srenew(dest->a, dest->nra);
1172   for(i=0; i<source->nra; i++)
1173     dest->a[nra0+i] = source->a[i];
1174
1175   /* terminate list */
1176   dest->index[dest->nr]=dest->nra;
1177
1178 }
1179
1180 int main(int argc,char *argv[])
1181 {
1182   const char *desc[] = {
1183     "Index groups are necessary for almost every gromacs program.",
1184     "All these programs can generate default index groups. You ONLY",
1185     "have to use make_ndx when you need SPECIAL index groups.",
1186     "There is a default index group for the whole system, 9 default",
1187     "index groups are generated for proteins, a default index group",
1188     "is generated for every other residue name.[PAR]",
1189     "When no index file is supplied, also make_ndx will generate the",
1190     "default groups.",
1191     "With the index editor you can select on atom, residue and chain names",
1192     "and numbers.",
1193     "When a run input file is supplied you can also select on atom type.",
1194     "You can use NOT, AND and OR, you can split groups",
1195     "into chains, residues or atoms. You can delete and rename groups.[PAR]",
1196     "The atom numbering in the editor and the index file starts at 1."
1197   };
1198
1199   static int natoms=0;
1200   static bool bVerbose=FALSE;
1201   t_pargs pa[] = {
1202     { "-natoms",  FALSE, etINT, {&natoms},
1203       "set number of atoms (default: read from coordinate or index file)" },
1204     { "-verbose", FALSE, etBOOL, {&bVerbose},
1205       "HIDDENVerbose output" }
1206   };
1207 #define NPA asize(pa)
1208
1209   output_env_t oenv;
1210   char     title[STRLEN];
1211   int      nndxin;
1212   const char *stxfile;
1213   char     **ndxinfiles;
1214   const char *ndxoutfile;
1215   bool     bNatoms;
1216   int      i,j;
1217   t_atoms  *atoms;
1218   rvec     *x,*v;
1219   int      ePBC;
1220   matrix   box;
1221   t_blocka *block,*block2;
1222   char     **gnames,**gnames2;
1223   t_filenm fnm[] = {
1224     { efSTX, "-f", NULL,     ffOPTRD  },
1225     { efNDX, "-n", NULL,     ffOPTRDMULT },
1226     { efNDX, "-o", NULL,     ffWRITE }
1227   };
1228 #define NFILE asize(fnm)
1229
1230   CopyRight(stderr,argv[0]);
1231
1232   parse_common_args(&argc,argv,0,NFILE,fnm,NPA,pa,asize(desc),desc,
1233                     0,NULL,&oenv);
1234
1235   stxfile = ftp2fn_null(efSTX,NFILE,fnm);
1236   if (opt2bSet("-n",NFILE,fnm)) {
1237     nndxin = opt2fns(&ndxinfiles,"-n",NFILE,fnm);
1238   } else {
1239     nndxin = 0;
1240   }
1241   ndxoutfile = opt2fn("-o",NFILE,fnm);
1242   bNatoms = opt2parg_bSet("-natoms",NPA,pa);
1243
1244   if (!stxfile && !nndxin)
1245     gmx_fatal(FARGS,"No input files (structure or index)");
1246
1247   if (stxfile) {
1248     snew(atoms,1);
1249     get_stx_coordnum(stxfile,&(atoms->nr));
1250     init_t_atoms(atoms,atoms->nr,TRUE);
1251     snew(x,atoms->nr);
1252     snew(v,atoms->nr);
1253     fprintf(stderr,"\nReading structure file\n");
1254     read_stx_conf(stxfile,title,atoms,x,v,&ePBC,box);
1255     natoms = atoms->nr;
1256     bNatoms=TRUE;
1257   } else {
1258     atoms = NULL;
1259     x = NULL;
1260   }
1261
1262   /* read input file(s) */
1263   block = new_blocka();
1264   gnames = NULL;
1265   printf("Going to read %d old index file(s)\n",nndxin);
1266   if (nndxin) {
1267     for(i=0; i<nndxin; i++) {
1268       block2 = init_index(ndxinfiles[i],&gnames2);
1269       srenew(gnames, block->nr+block2->nr);
1270       for(j=0; j<block2->nr; j++)
1271         gnames[block->nr+j]=gnames2[j];
1272       sfree(gnames2);
1273       merge_blocks(block, block2);
1274       sfree(block2->a);
1275       sfree(block2->index);
1276 /*       done_block(block2); */
1277       sfree(block2);
1278     }
1279   }
1280   else {
1281     snew(gnames,1);
1282     analyse(atoms,block,&gnames,FALSE,TRUE);
1283   }
1284
1285   if (!bNatoms) {
1286     natoms = block2natoms(block);
1287     printf("Counted atom numbers up to %d in index file\n", natoms);
1288   }
1289
1290   edit_index(natoms,atoms,x,block,&gnames,bVerbose);
1291
1292   write_index(ndxoutfile,block,gnames);
1293
1294   thanx(stderr);
1295
1296   return 0;
1297 }