2 * This file is part of the GROMACS molecular simulation package.
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.
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.
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.
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.
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.
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.
44 #include "gromacs/utility/futil.h"
45 #include "gromacs/legacyheaders/macros.h"
46 #include "gromacs/fileio/confio.h"
47 #include "gromacs/legacyheaders/typedefs.h"
48 #include "gromacs/topology/index.h"
49 #include "gromacs/math/vec.h"
51 #include "gromacs/commandline/pargs.h"
52 #include "gromacs/topology/block.h"
53 #include "gromacs/utility/cstringutil.h"
54 #include "gromacs/utility/fatalerror.h"
55 #include "gromacs/utility/smalloc.h"
57 /* It's not nice to have size limits, but we should not spend more time
58 * on this ancient tool, but instead use the new selection library.
63 gmx_bool bCase = FALSE;
65 static int or_groups(atom_id nr1, atom_id *at1, atom_id nr2, atom_id *at2,
66 atom_id *nr, atom_id *at)
68 atom_id i1, i2, max = 0;
74 for (i1 = 0; i1 < nr1; i1++)
76 if ((i1 > 0) && (at1[i1] <= max))
82 for (i1 = 0; i1 < nr2; i1++)
84 if ((i1 > 0) && (at2[i1] <= max))
93 printf("One of your groups is not ascending\n");
100 while ((i1 < nr1) || (i2 < nr2))
102 if ((i2 == nr2) || ((i1 < nr1) && (at1[i1] < at2[i2])))
110 if ((i2 < nr2) && ((i1 == nr1) || (at1[i1] > at2[i2])))
119 printf("Merged two groups with OR: %d %d -> %d\n", nr1, nr2, *nr);
125 static int and_groups(atom_id nr1, atom_id *at1, atom_id nr2, atom_id *at2,
126 atom_id *nr, atom_id *at)
131 for (i1 = 0; i1 < nr1; i1++)
133 for (i2 = 0; i2 < nr2; i2++)
135 if (at1[i1] == at2[i2])
143 printf("Merged two groups with AND: %d %d -> %d\n", nr1, nr2, *nr);
148 static gmx_bool is_name_char(char c)
150 /* This string should contain all characters that can not be
151 * the first letter of a name due to the make_ndx syntax.
153 const char *spec = " !&|";
155 return (c != '\0' && strchr(spec, c) == NULL);
158 static int parse_names(char **string, int *n_names, char **names)
163 while (is_name_char((*string)[0]) || (*string)[0] == ' ')
165 if (is_name_char((*string)[0]))
167 if (*n_names >= MAXNAMES)
169 gmx_fatal(FARGS, "To many names: %d\n", *n_names+1);
172 while (is_name_char((*string)[i]))
174 names[*n_names][i] = (*string)[i];
178 printf("Name is too long, the maximum is %d characters\n", NAME_LEN);
182 names[*n_names][i] = '\0';
185 upstring(names[*n_names]);
199 static gmx_bool parse_int_char(char **string, int *nr, char *c)
206 while ((*string)[0] == ' ')
215 if (isdigit((*string)[0]))
217 *nr = (*string)[0]-'0';
219 while (isdigit((*string)[0]))
221 *nr = (*nr)*10+(*string)[0]-'0';
224 if (isalpha((*string)[0]))
229 /* Check if there is at most one non-digit character */
230 if (!isalnum((*string)[0]))
247 static gmx_bool parse_int(char **string, int *nr)
253 bRet = parse_int_char(string, nr, &c);
254 if (bRet && c != ' ')
263 static gmx_bool isquote(char c)
268 static gmx_bool parse_string(char **string, int *nr, int ngrps, char **grpname)
273 while ((*string)[0] == ' ')
279 if (isquote((*string)[0]))
283 s = gmx_strdup((*string));
287 (*string) += sp-s + 1;
289 (*nr) = find_group(s, ngrps, grpname);
293 return (*nr) != NOTSET;
296 static int select_atomnumbers(char **string, t_atoms *atoms, atom_id n1,
297 atom_id *nr, atom_id *index, char *gname)
303 while ((*string)[0] == ' ')
307 if ((*string)[0] == '-')
310 parse_int(string, &up);
311 if ((n1 < 1) || (n1 > atoms->nr) || (up < 1) || (up > atoms->nr))
313 printf("Invalid atom range\n");
317 for (i = n1-1; i <= up-1; i++)
322 printf("Found %d atom%s in range %d-%d\n", *nr, (*nr == 1) ? "" : "s", n1, up);
325 sprintf(buf, "a_%d", n1);
329 sprintf(buf, "a_%d-%d", n1, up);
340 if ((i-1 >= 0) && (i-1 < atoms->nr))
344 sprintf(buf, "_%d", i);
349 printf("Invalid atom number %d\n", i);
353 while ((*nr != 0) && (parse_int(string, &i)));
359 static int select_residuenumbers(char **string, t_atoms *atoms,
361 atom_id *nr, atom_id *index, char *gname)
368 while ((*string)[0] == ' ')
372 if ((*string)[0] == '-')
374 /* Residue number range selection */
377 printf("Error: residue insertion codes can not be used with residue range selection\n");
381 parse_int(string, &up);
383 for (i = 0; i < atoms->nr; i++)
385 ri = &atoms->resinfo[atoms->atom[i].resind];
386 for (j = n1; (j <= up); j++)
388 if (ri->nr == j && (c == ' ' || ri->ic == c))
395 printf("Found %d atom%s with res.nr. in range %d-%d\n",
396 *nr, (*nr == 1) ? "" : "s", n1, up);
399 sprintf(buf, "r_%d", n1);
403 sprintf(buf, "r_%d-%d", n1, up);
409 /* Individual residue number/insertion code selection */
414 for (i = 0; i < atoms->nr; i++)
416 ri = &atoms->resinfo[atoms->atom[i].resind];
417 if (ri->nr == j && ri->ic == c)
423 sprintf(buf, "_%d", j);
426 while (parse_int_char(string, &j, &c));
432 static int select_residueindices(char **string, t_atoms *atoms,
434 atom_id *nr, atom_id *index, char *gname)
436 /*this should be similar to select_residuenumbers except select by index (sequential numbering in file)*/
437 /*resind+1 for 1-indexing*/
443 while ((*string)[0] == ' ')
447 if ((*string)[0] == '-')
449 /* Residue number range selection */
452 printf("Error: residue insertion codes can not be used with residue range selection\n");
456 parse_int(string, &up);
458 for (i = 0; i < atoms->nr; i++)
460 ri = &atoms->resinfo[atoms->atom[i].resind];
461 for (j = n1; (j <= up); j++)
463 if (atoms->atom[i].resind+1 == j && (c == ' ' || ri->ic == c))
470 printf("Found %d atom%s with resind.+1 in range %d-%d\n",
471 *nr, (*nr == 1) ? "" : "s", n1, up);
474 sprintf(buf, "r_%d", n1);
478 sprintf(buf, "r_%d-%d", n1, up);
484 /* Individual residue number/insertion code selection */
489 for (i = 0; i < atoms->nr; i++)
491 ri = &atoms->resinfo[atoms->atom[i].resind];
492 if (atoms->atom[i].resind+1 == j && ri->ic == c)
498 sprintf(buf, "_%d", j);
501 while (parse_int_char(string, &j, &c));
508 static gmx_bool atoms_from_residuenumbers(t_atoms *atoms, int group, t_blocka *block,
509 atom_id *nr, atom_id *index, char *gname)
511 int i, j, j0, j1, resnr, nres;
513 j0 = block->index[group];
514 j1 = block->index[group+1];
516 for (j = j0; j < j1; j++)
518 if (block->a[j] >= nres)
520 printf("Index %s contains number>nres (%d>%d)\n",
521 gname, block->a[j]+1, nres);
525 for (i = 0; i < atoms->nr; i++)
527 resnr = atoms->resinfo[atoms->atom[i].resind].nr;
528 for (j = j0; j < j1; j++)
530 if (block->a[j]+1 == resnr)
538 printf("Found %d atom%s in %d residues from group %s\n",
539 *nr, (*nr == 1) ? "" : "s", j1-j0, gname);
543 static gmx_bool comp_name(char *name, char *search)
545 while (name[0] != '\0' && search[0] != '\0')
553 if (search[1] != '\0')
555 printf("WARNING: Currently '*' is only supported at the end of an expression\n");
564 /* Compare a single character */
565 if (( bCase && strncmp(name, search, 1)) ||
566 (!bCase && gmx_strncasecmp(name, search, 1)))
575 return (name[0] == '\0' && (search[0] == '\0' || search[0] == '*'));
578 static int select_chainnames(t_atoms *atoms, int n_names, char **names,
579 atom_id *nr, atom_id *index)
587 for (i = 0; i < atoms->nr; i++)
589 name[0] = atoms->resinfo[atoms->atom[i].resind].chainid;
591 while (j < n_names && !comp_name(name, names[j]))
601 printf("Found %d atom%s with chain identifier%s",
602 *nr, (*nr == 1) ? "" : "s", (n_names == 1) ? "" : "s");
603 for (j = 0; (j < n_names); j++)
605 printf(" %s", names[j]);
612 static int select_atomnames(t_atoms *atoms, int n_names, char **names,
613 atom_id *nr, atom_id *index, gmx_bool bType)
620 for (i = 0; i < atoms->nr; i++)
624 name = *(atoms->atomtype[i]);
628 name = *(atoms->atomname[i]);
631 while (j < n_names && !comp_name(name, names[j]))
641 printf("Found %d atoms with %s%s",
642 *nr, bType ? "type" : "name", (n_names == 1) ? "" : "s");
643 for (j = 0; (j < n_names); j++)
645 printf(" %s", names[j]);
652 static int select_residuenames(t_atoms *atoms, int n_names, char **names,
653 atom_id *nr, atom_id *index)
660 for (i = 0; i < atoms->nr; i++)
662 name = *(atoms->resinfo[atoms->atom[i].resind].name);
664 while (j < n_names && !comp_name(name, names[j]))
674 printf("Found %d atoms with residue name%s", *nr, (n_names == 1) ? "" : "s");
675 for (j = 0; (j < n_names); j++)
677 printf(" %s", names[j]);
684 static void copy2block(int n, atom_id *index, t_blocka *block)
691 srenew(block->index, block->nr+1);
692 block->index[block->nr] = n0+n;
693 srenew(block->a, n0+n);
694 for (i = 0; (i < n); i++)
696 block->a[n0+i] = index[i];
700 static void make_gname(int n, char **names, char *gname)
704 strcpy(gname, names[0]);
705 for (i = 1; i < n; i++)
708 strcat(gname, names[i]);
712 static void copy_group(int g, t_blocka *block, atom_id *nr, atom_id *index)
716 i0 = block->index[g];
717 *nr = block->index[g+1]-i0;
718 for (i = 0; i < *nr; i++)
720 index[i] = block->a[i0+i];
724 static void remove_group(int nr, int nr2, t_blocka *block, char ***gn)
734 for (j = 0; j <= nr2-nr; j++)
736 if ((nr < 0) || (nr >= block->nr))
738 printf("Group %d does not exist\n", nr+j);
742 shift = block->index[nr+1]-block->index[nr];
743 for (i = block->index[nr+1]; i < block->nra; i++)
745 block->a[i-shift] = block->a[i];
748 for (i = nr; i < block->nr; i++)
750 block->index[i] = block->index[i+1]-shift;
752 name = gmx_strdup((*gn)[nr]);
754 for (i = nr; i < block->nr-1; i++)
756 (*gn)[i] = (*gn)[i+1];
759 block->nra = block->index[block->nr];
760 printf("Removed group %d '%s'\n", nr+j, name);
766 static void split_group(t_atoms *atoms, int sel_nr, t_blocka *block, char ***gn,
769 char buf[STRLEN], *name;
773 printf("Splitting group %d '%s' into %s\n", sel_nr, (*gn)[sel_nr],
774 bAtom ? "atoms" : "residues");
776 n0 = block->index[sel_nr];
777 n1 = block->index[sel_nr+1];
778 srenew(block->a, block->nra+n1-n0);
779 for (i = n0; i < n1; i++)
782 resind = atoms->atom[a].resind;
783 name = *(atoms->resinfo[resind].name);
784 if (bAtom || (i == n0) || (atoms->atom[block->a[i-1]].resind != resind))
788 block->index[block->nr] = block->nra;
791 srenew(block->index, block->nr+1);
792 srenew(*gn, block->nr);
795 sprintf(buf, "%s_%s_%d", (*gn)[sel_nr], *atoms->atomname[a], a+1);
799 sprintf(buf, "%s_%s_%d", (*gn)[sel_nr], name, atoms->resinfo[resind].nr);
801 (*gn)[block->nr-1] = gmx_strdup(buf);
803 block->a[block->nra] = a;
806 block->index[block->nr] = block->nra;
809 static int split_chain(t_atoms *atoms, rvec *x,
810 int sel_nr, t_blocka *block, char ***gn)
814 atom_id i, a, natoms, *start = NULL, *end = NULL, ca_start, ca_end;
821 while (ca_start < natoms)
823 while ((ca_start < natoms) && strcmp(*atoms->atomname[ca_start], "CA"))
827 if (ca_start < natoms)
829 srenew(start, nchain+1);
830 srenew(end, nchain+1);
831 start[nchain] = ca_start;
832 while ((start[nchain] > 0) &&
833 (atoms->atom[start[nchain]-1].resind ==
834 atoms->atom[ca_start].resind))
847 while ((i < natoms) && strcmp(*atoms->atomname[i], "CA"));
850 rvec_sub(x[ca_end], x[i], vec);
857 while (norm(vec) < 0.45);
859 end[nchain] = ca_end;
860 while ((end[nchain]+1 < natoms) &&
861 (atoms->atom[end[nchain]+1].resind == atoms->atom[ca_end].resind))
865 ca_start = end[nchain]+1;
871 printf("Found 1 chain, will not split\n");
875 printf("Found %d chains\n", nchain);
877 for (j = 0; j < nchain; j++)
879 printf("%d:%6d atoms (%d to %d)\n",
880 j+1, end[j]-start[j]+1, start[j]+1, end[j]+1);
885 srenew(block->a, block->nra+block->index[sel_nr+1]-block->index[sel_nr]);
886 for (j = 0; j < nchain; j++)
889 srenew(block->index, block->nr+1);
890 srenew(*gn, block->nr);
891 sprintf(buf, "%s_chain%d", (*gn)[sel_nr], j+1);
892 (*gn)[block->nr-1] = gmx_strdup(buf);
893 for (i = block->index[sel_nr]; i < block->index[sel_nr+1]; i++)
896 if ((a >= start[j]) && (a <= end[j]))
898 block->a[block->nra] = a;
902 block->index[block->nr] = block->nra;
903 if (block->index[block->nr-1] == block->index[block->nr])
905 remove_group(block->nr-1, NOTSET, block, gn);
915 static gmx_bool check_have_atoms(t_atoms *atoms, char *string)
919 printf("Can not process '%s' without atom info, use option -f\n", string);
928 static gmx_bool parse_entry(char **string, int natoms, t_atoms *atoms,
929 t_blocka *block, char ***gn,
930 atom_id *nr, atom_id *index, char *gname)
932 static char **names, *ostring;
933 static gmx_bool bFirst = TRUE;
934 int j, n_names, sel_nr1;
935 atom_id i, nr1, *index1;
937 gmx_bool bRet, bCompl;
942 snew(names, MAXNAMES);
943 for (i = 0; i < MAXNAMES; i++)
945 snew(names[i], NAME_LEN+1);
952 while (*string[0] == ' ')
957 if ((*string)[0] == '!')
961 while (*string[0] == ' ')
973 if (parse_int(string, &sel_nr1) ||
974 parse_string(string, &sel_nr1, block->nr, *gn))
976 if ((sel_nr1 >= 0) && (sel_nr1 < block->nr))
978 copy_group(sel_nr1, block, nr, index);
979 strcpy(gname, (*gn)[sel_nr1]);
980 printf("Copied index group %d '%s'\n", sel_nr1, (*gn)[sel_nr1]);
985 printf("Group %d does not exist\n", sel_nr1);
988 else if ((*string)[0] == 'a')
991 if (check_have_atoms(atoms, ostring))
993 if (parse_int(string, &sel_nr1))
995 bRet = select_atomnumbers(string, atoms, sel_nr1, nr, index, gname);
997 else if (parse_names(string, &n_names, names))
999 bRet = select_atomnames(atoms, n_names, names, nr, index, FALSE);
1000 make_gname(n_names, names, gname);
1004 else if ((*string)[0] == 't')
1007 if (check_have_atoms(atoms, ostring) &&
1008 parse_names(string, &n_names, names))
1010 if (atoms->atomtype == NULL)
1012 printf("Need a run input file to select atom types\n");
1016 bRet = select_atomnames(atoms, n_names, names, nr, index, TRUE);
1017 make_gname(n_names, names, gname);
1021 else if (strncmp(*string, "res", 3) == 0)
1024 if (check_have_atoms(atoms, ostring) &&
1025 parse_int(string, &sel_nr1) &&
1026 (sel_nr1 >= 0) && (sel_nr1 < block->nr) )
1028 bRet = atoms_from_residuenumbers(atoms,
1029 sel_nr1, block, nr, index, (*gn)[sel_nr1]);
1030 sprintf(gname, "atom_%s", (*gn)[sel_nr1]);
1033 else if (strncmp(*string, "ri", 2) == 0)
1036 if (check_have_atoms(atoms, ostring) &&
1037 parse_int_char(string, &sel_nr1, &c))
1039 bRet = select_residueindices(string, atoms, sel_nr1, c, nr, index, gname);
1042 else if ((*string)[0] == 'r')
1045 if (check_have_atoms(atoms, ostring))
1047 if (parse_int_char(string, &sel_nr1, &c))
1049 bRet = select_residuenumbers(string, atoms, sel_nr1, c, nr, index, gname);
1051 else if (parse_names(string, &n_names, names))
1053 bRet = select_residuenames(atoms, n_names, names, nr, index);
1054 make_gname(n_names, names, gname);
1058 else if (strncmp(*string, "chain", 5) == 0)
1061 if (check_have_atoms(atoms, ostring) &&
1062 parse_names(string, &n_names, names))
1064 bRet = select_chainnames(atoms, n_names, names, nr, index);
1065 sprintf(gname, "ch%s", names[0]);
1066 for (i = 1; i < n_names; i++)
1068 strcat(gname, names[i]);
1074 snew(index1, natoms-*nr);
1076 for (i = 0; i < natoms; i++)
1079 while ((j < *nr) && (index[j] != i))
1085 if (nr1 >= natoms-*nr)
1087 printf("There are double atoms in your index group\n");
1095 for (i = 0; i < nr1; i++)
1097 index[i] = index1[i];
1101 for (i = strlen(gname)+1; i > 0; i--)
1103 gname[i] = gname[i-1];
1106 printf("Complemented group: %d atoms\n", *nr);
1112 static void list_residues(t_atoms *atoms)
1114 int i, j, start, end, prev_resind, resind;
1117 /* Print all the residues, assuming continuous resnr count */
1118 start = atoms->atom[0].resind;
1119 prev_resind = start;
1120 for (i = 0; i < atoms->nr; i++)
1122 resind = atoms->atom[i].resind;
1123 if ((resind != prev_resind) || (i == atoms->nr-1))
1125 if ((bDiff = strcmp(*atoms->resinfo[resind].name,
1126 *atoms->resinfo[start].name)) ||
1139 for (j = start; j <= end; j++)
1142 j+1, *(atoms->resinfo[j].name));
1147 printf(" %4d - %4d %-5s ",
1148 start+1, end+1, *(atoms->resinfo[start].name));
1153 prev_resind = resind;
1158 static void edit_index(int natoms, t_atoms *atoms, rvec *x, t_blocka *block, char ***gn, gmx_bool bVerbose)
1160 static char **atnames, *ostring;
1161 static gmx_bool bFirst = TRUE;
1162 char inp_string[STRLEN], *string;
1163 char gname[STRLEN], gname1[STRLEN], gname2[STRLEN];
1164 int i, i0, i1, sel_nr, sel_nr2, newgroup;
1165 atom_id nr, nr1, nr2, *index, *index1, *index2;
1166 gmx_bool bAnd, bOr, bPrintOnce;
1171 snew(atnames, MAXNAMES);
1172 for (i = 0; i < MAXNAMES; i++)
1174 snew(atnames[i], NAME_LEN+1);
1180 snew(index, natoms);
1181 snew(index1, natoms);
1182 snew(index2, natoms);
1189 if (bVerbose || bPrintOnce || newgroup != NOTSET)
1192 if (bVerbose || bPrintOnce || newgroup == NOTSET)
1202 for (i = i0; i < i1; i++)
1204 printf("%3d %-20s: %5d atoms\n", i, (*gn)[i],
1205 block->index[i+1]-block->index[i]);
1209 if (bVerbose || bPrintOnce)
1212 printf(" nr : group ! 'name' nr name 'splitch' nr Enter: list groups\n");
1213 printf(" 'a': atom & 'del' nr 'splitres' nr 'l': list residues\n");
1214 printf(" 't': atom type | 'keep' nr 'splitat' nr 'h': help\n");
1215 printf(" 'r': residue 'res' nr 'chain' char\n");
1216 printf(" \"name\": group 'case': case %s 'q': save and quit\n",
1217 bCase ? "insensitive" : "sensitive ");
1218 printf(" 'ri': residue index\n");
1223 if (NULL == fgets(inp_string, STRLEN, stdin))
1225 gmx_fatal(FARGS, "Error reading user input");
1227 inp_string[strlen(inp_string)-1] = 0;
1229 string = inp_string;
1230 while (string[0] == ' ')
1237 if (string[0] == 'h')
1239 printf(" nr : selects an index group by number or quoted string.\n");
1240 printf(" The string is first matched against the whole group name,\n");
1241 printf(" then against the beginning and finally against an\n");
1242 printf(" arbitrary substring. A multiple match is an error.\n");
1244 printf(" 'a' nr1 [nr2 ...] : selects atoms, atom numbering starts at 1.\n");
1245 printf(" 'a' nr1 - nr2 : selects atoms in the range from nr1 to nr2.\n");
1246 printf(" 'a' name1[*] [name2[*] ...] : selects atoms by name(s), '?' matches any char,\n");
1247 printf(" wildcard '*' allowed at the end of a name.\n");
1248 printf(" 't' type1[*] [type2[*] ...] : as 'a', but for type, run input file required.\n");
1249 printf(" 'r' nr1[ic1] [nr2[ic2] ...] : selects residues by number and insertion code.\n");
1250 printf(" 'r' nr1 - nr2 : selects residues in the range from nr1 to nr2.\n");
1251 printf(" 'r' name1[*] [name2[*] ...] : as 'a', but for residue names.\n");
1252 printf(" 'ri' nr1 - nr2 : selects residue indices, 1-indexed, (as opposed to numbers) in the range from nr1 to nr2.\n");
1253 printf(" 'chain' ch1 [ch2 ...] : selects atoms by chain identifier(s),\n");
1254 printf(" not available with a .gro file as input.\n");
1255 printf(" ! : takes the complement of a group with respect to all\n");
1256 printf(" the atoms in the input file.\n");
1257 printf(" & | : AND and OR, can be placed between any of the options\n");
1258 printf(" above, the input is processed from left to right.\n");
1259 printf(" 'name' nr name : rename group nr to name.\n");
1260 printf(" 'del' nr1 [- nr2] : deletes one group or groups in the range from nr1 to nr2.\n");
1261 printf(" 'keep' nr : deletes all groups except nr.\n");
1262 printf(" 'case' : make all name compares case (in)sensitive.\n");
1263 printf(" 'splitch' nr : split group into chains using CA distances.\n");
1264 printf(" 'splitres' nr : split group into residues.\n");
1265 printf(" 'splitat' nr : split group into atoms.\n");
1266 printf(" 'res' nr : interpret numbers in group as residue numbers\n");
1267 printf(" Enter : list the currently defined groups and commands\n");
1268 printf(" 'l' : list the residues.\n");
1269 printf(" 'h' : show this help.\n");
1270 printf(" 'q' : save and quit.\n");
1272 printf(" Examples:\n");
1273 printf(" > 2 | 4 & r 3-5\n");
1274 printf(" selects all atoms from group 2 and 4 that have residue numbers 3, 4 or 5\n");
1275 printf(" > a C* & !a C CA\n");
1276 printf(" selects all atoms starting with 'C' but not the atoms 'C' and 'CA'\n");
1277 printf(" > \"protein\" & ! \"backb\"\n");
1278 printf(" selects all atoms that are in group 'protein' and not in group 'backbone'\n");
1281 printf("\npress Enter ");
1285 else if (strncmp(string, "del", 3) == 0)
1288 if (parse_int(&string, &sel_nr))
1290 while (string[0] == ' ')
1294 if (string[0] == '-')
1297 parse_int(&string, &sel_nr2);
1303 while (string[0] == ' ')
1307 if (string[0] == '\0')
1309 remove_group(sel_nr, sel_nr2, block, gn);
1313 printf("\nSyntax error: \"%s\"\n", string);
1317 else if (strncmp(string, "keep", 4) == 0)
1320 if (parse_int(&string, &sel_nr))
1322 remove_group(sel_nr+1, block->nr-1, block, gn);
1323 remove_group(0, sel_nr-1, block, gn);
1326 else if (strncmp(string, "name", 4) == 0)
1329 if (parse_int(&string, &sel_nr))
1331 if ((sel_nr >= 0) && (sel_nr < block->nr))
1333 sscanf(string, "%s", gname);
1334 sfree((*gn)[sel_nr]);
1335 (*gn)[sel_nr] = gmx_strdup(gname);
1339 else if (strncmp(string, "case", 4) == 0)
1342 printf("Switched to case %s\n", bCase ? "sensitive" : "insensitive");
1344 else if (string[0] == 'v')
1346 bVerbose = !bVerbose;
1347 printf("Turned verbose %s\n", bVerbose ? "on" : "off");
1349 else if (string[0] == 'l')
1351 if (check_have_atoms(atoms, ostring) )
1353 list_residues(atoms);
1356 else if (strncmp(string, "splitch", 7) == 0)
1359 if (check_have_atoms(atoms, ostring) &&
1360 parse_int(&string, &sel_nr) &&
1361 (sel_nr >= 0) && (sel_nr < block->nr))
1363 split_chain(atoms, x, sel_nr, block, gn);
1366 else if (strncmp(string, "splitres", 8) == 0)
1369 if (check_have_atoms(atoms, ostring) &&
1370 parse_int(&string, &sel_nr) &&
1371 (sel_nr >= 0) && (sel_nr < block->nr))
1373 split_group(atoms, sel_nr, block, gn, FALSE);
1376 else if (strncmp(string, "splitat", 7) == 0)
1379 if (check_have_atoms(atoms, ostring) &&
1380 parse_int(&string, &sel_nr) &&
1381 (sel_nr >= 0) && (sel_nr < block->nr))
1383 split_group(atoms, sel_nr, block, gn, TRUE);
1386 else if (string[0] == '\0')
1390 else if (string[0] != 'q')
1394 if (parse_entry(&string, natoms, atoms, block, gn, &nr, index, gname))
1398 while (string[0] == ' ')
1405 if (string[0] == '&')
1409 else if (string[0] == '|')
1418 for (i = 0; i < nr; i++)
1420 index1[i] = index[i];
1422 strcpy(gname1, gname);
1423 if (parse_entry(&string, natoms, atoms, block, gn, &nr2, index2, gname2))
1427 or_groups(nr1, index1, nr2, index2, &nr, index);
1428 sprintf(gname, "%s_%s", gname1, gname2);
1432 and_groups(nr1, index1, nr2, index2, &nr, index);
1433 sprintf(gname, "%s_&_%s", gname1, gname2);
1438 while (bAnd || bOr);
1440 while (string[0] == ' ')
1446 printf("\nSyntax error: \"%s\"\n", string);
1450 copy2block(nr, index, block);
1451 srenew(*gn, block->nr);
1452 newgroup = block->nr-1;
1453 (*gn)[newgroup] = gmx_strdup(gname);
1457 printf("Group is empty\n");
1461 while (string[0] != 'q');
1468 static int block2natoms(t_blocka *block)
1473 for (i = 0; i < block->nra; i++)
1475 natoms = max(natoms, block->a[i]);
1482 void merge_blocks(t_blocka *dest, t_blocka *source)
1486 /* count groups, srenew and fill */
1489 dest->nr += source->nr;
1490 srenew(dest->index, dest->nr+1);
1491 for (i = 0; i < source->nr; i++)
1493 dest->index[i0+i] = nra0 + source->index[i];
1495 /* count atoms, srenew and fill */
1496 dest->nra += source->nra;
1497 srenew(dest->a, dest->nra);
1498 for (i = 0; i < source->nra; i++)
1500 dest->a[nra0+i] = source->a[i];
1503 /* terminate list */
1504 dest->index[dest->nr] = dest->nra;
1508 int gmx_make_ndx(int argc, char *argv[])
1510 const char *desc[] = {
1511 "Index groups are necessary for almost every GROMACS program.",
1512 "All these programs can generate default index groups. You ONLY",
1513 "have to use [THISMODULE] when you need SPECIAL index groups.",
1514 "There is a default index group for the whole system, 9 default",
1515 "index groups for proteins, and a default index group",
1516 "is generated for every other residue name.[PAR]",
1517 "When no index file is supplied, also [THISMODULE] will generate the",
1519 "With the index editor you can select on atom, residue and chain names",
1521 "When a run input file is supplied you can also select on atom type.",
1522 "You can use NOT, AND and OR, you can split groups",
1523 "into chains, residues or atoms. You can delete and rename groups.[PAR]",
1524 "The atom numbering in the editor and the index file starts at 1.[PAR]",
1525 "The [TT]-twin[tt] switch duplicates all index groups with an offset of",
1526 "[TT]-natoms[tt], which is useful for Computational Electrophysiology",
1527 "double-layer membrane setups."
1530 static int natoms = 0;
1531 static gmx_bool bVerbose = FALSE;
1532 static gmx_bool bDuplicate = FALSE;
1534 { "-natoms", FALSE, etINT, {&natoms},
1535 "set number of atoms (default: read from coordinate or index file)" },
1536 { "-twin", FALSE, etBOOL, {&bDuplicate},
1537 "Duplicate all index groups with an offset of -natoms" },
1538 { "-verbose", FALSE, etBOOL, {&bVerbose},
1539 "HIDDENVerbose output" }
1541 #define NPA asize(pa)
1546 const char *stxfile;
1548 const char *ndxoutfile;
1555 t_blocka *block, *block2;
1556 char **gnames, **gnames2;
1558 { efSTX, "-f", NULL, ffOPTRD },
1559 { efNDX, "-n", NULL, ffOPTRDMULT },
1560 { efNDX, "-o", NULL, ffWRITE }
1562 #define NFILE asize(fnm)
1564 if (!parse_common_args(&argc, argv, 0, NFILE, fnm, NPA, pa, asize(desc), desc,
1570 stxfile = ftp2fn_null(efSTX, NFILE, fnm);
1571 if (opt2bSet("-n", NFILE, fnm))
1573 nndxin = opt2fns(&ndxinfiles, "-n", NFILE, fnm);
1579 ndxoutfile = opt2fn("-o", NFILE, fnm);
1580 bNatoms = opt2parg_bSet("-natoms", NPA, pa);
1582 if (!stxfile && !nndxin)
1584 gmx_fatal(FARGS, "No input files (structure or index)");
1590 get_stx_coordnum(stxfile, &(atoms->nr));
1591 init_t_atoms(atoms, atoms->nr, TRUE);
1594 fprintf(stderr, "\nReading structure file\n");
1595 read_stx_conf(stxfile, title, atoms, x, v, &ePBC, box);
1605 /* read input file(s) */
1606 block = new_blocka();
1608 printf("Going to read %d old index file(s)\n", nndxin);
1611 for (i = 0; i < nndxin; i++)
1613 block2 = init_index(ndxinfiles[i], &gnames2);
1614 srenew(gnames, block->nr+block2->nr);
1615 for (j = 0; j < block2->nr; j++)
1617 gnames[block->nr+j] = gnames2[j];
1620 merge_blocks(block, block2);
1622 sfree(block2->index);
1623 /* done_block(block2); */
1630 analyse(atoms, block, &gnames, FALSE, TRUE);
1635 natoms = block2natoms(block);
1636 printf("Counted atom numbers up to %d in index file\n", natoms);
1639 edit_index(natoms, atoms, x, block, &gnames, bVerbose);
1641 write_index(ndxoutfile, block, gnames, bDuplicate, natoms);