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.
45 #include "gromacs/fileio/futil.h"
47 #include "gromacs/utility/cstringutil.h"
48 #include "gromacs/commandline/pargs.h"
49 #include "gromacs/fileio/confio.h"
52 #include "gromacs/utility/smalloc.h"
56 #include "gromacs/legacyheaders/gmx_fatal.h"
58 /* It's not nice to have size limits, but we should not spend more time
59 * on this ancient tool, but instead use the new selection library.
64 gmx_bool bCase = FALSE;
66 static int or_groups(atom_id nr1, atom_id *at1, atom_id nr2, atom_id *at2,
67 atom_id *nr, atom_id *at)
69 atom_id i1, i2, max = 0;
75 for (i1 = 0; i1 < nr1; i1++)
77 if ((i1 > 0) && (at1[i1] <= max))
83 for (i1 = 0; i1 < nr2; i1++)
85 if ((i1 > 0) && (at2[i1] <= max))
94 printf("One of your groups is not ascending\n");
101 while ((i1 < nr1) || (i2 < nr2))
103 if ((i2 == nr2) || ((i1 < nr1) && (at1[i1] < at2[i2])))
111 if ((i2 < nr2) && ((i1 == nr1) || (at1[i1] > at2[i2])))
120 printf("Merged two groups with OR: %d %d -> %d\n", nr1, nr2, *nr);
126 static int and_groups(atom_id nr1, atom_id *at1, atom_id nr2, atom_id *at2,
127 atom_id *nr, atom_id *at)
132 for (i1 = 0; i1 < nr1; i1++)
134 for (i2 = 0; i2 < nr2; i2++)
136 if (at1[i1] == at2[i2])
144 printf("Merged two groups with AND: %d %d -> %d\n", nr1, nr2, *nr);
149 static gmx_bool is_name_char(char c)
151 /* This string should contain all characters that can not be
152 * the first letter of a name due to the make_ndx syntax.
154 const char *spec = " !&|";
156 return (c != '\0' && strchr(spec, c) == NULL);
159 static int parse_names(char **string, int *n_names, char **names)
164 while (is_name_char((*string)[0]) || (*string)[0] == ' ')
166 if (is_name_char((*string)[0]))
168 if (*n_names >= MAXNAMES)
170 gmx_fatal(FARGS, "To many names: %d\n", *n_names+1);
173 while (is_name_char((*string)[i]))
175 names[*n_names][i] = (*string)[i];
179 printf("Name is too long, the maximum is %d characters\n", NAME_LEN);
183 names[*n_names][i] = '\0';
186 upstring(names[*n_names]);
200 static gmx_bool parse_int_char(char **string, int *nr, char *c)
207 while ((*string)[0] == ' ')
216 if (isdigit((*string)[0]))
218 *nr = (*string)[0]-'0';
220 while (isdigit((*string)[0]))
222 *nr = (*nr)*10+(*string)[0]-'0';
225 if (isalpha((*string)[0]))
230 /* Check if there is at most one non-digit character */
231 if (!isalnum((*string)[0]))
248 static gmx_bool parse_int(char **string, int *nr)
254 bRet = parse_int_char(string, nr, &c);
255 if (bRet && c != ' ')
264 static gmx_bool isquote(char c)
269 static gmx_bool parse_string(char **string, int *nr, int ngrps, char **grpname)
274 while ((*string)[0] == ' ')
280 if (isquote((*string)[0]))
284 s = strdup((*string));
288 (*string) += sp-s + 1;
290 (*nr) = find_group(s, ngrps, grpname);
294 return (*nr) != NOTSET;
297 static int select_atomnumbers(char **string, t_atoms *atoms, atom_id n1,
298 atom_id *nr, atom_id *index, char *gname)
304 while ((*string)[0] == ' ')
308 if ((*string)[0] == '-')
311 parse_int(string, &up);
312 if ((n1 < 1) || (n1 > atoms->nr) || (up < 1) || (up > atoms->nr))
314 printf("Invalid atom range\n");
318 for (i = n1-1; i <= up-1; i++)
323 printf("Found %d atom%s in range %d-%d\n", *nr, (*nr == 1) ? "" : "s", n1, up);
326 sprintf(buf, "a_%d", n1);
330 sprintf(buf, "a_%d-%d", n1, up);
341 if ((i-1 >= 0) && (i-1 < atoms->nr))
345 sprintf(buf, "_%d", i);
350 printf("Invalid atom number %d\n", i);
354 while ((*nr != 0) && (parse_int(string, &i)));
360 static int select_residuenumbers(char **string, t_atoms *atoms,
362 atom_id *nr, atom_id *index, char *gname)
369 while ((*string)[0] == ' ')
373 if ((*string)[0] == '-')
375 /* Residue number range selection */
378 printf("Error: residue insertion codes can not be used with residue range selection\n");
382 parse_int(string, &up);
384 for (i = 0; i < atoms->nr; i++)
386 ri = &atoms->resinfo[atoms->atom[i].resind];
387 for (j = n1; (j <= up); j++)
389 if (ri->nr == j && (c == ' ' || ri->ic == c))
396 printf("Found %d atom%s with res.nr. in range %d-%d\n",
397 *nr, (*nr == 1) ? "" : "s", n1, up);
400 sprintf(buf, "r_%d", n1);
404 sprintf(buf, "r_%d-%d", n1, up);
410 /* Individual residue number/insertion code selection */
415 for (i = 0; i < atoms->nr; i++)
417 ri = &atoms->resinfo[atoms->atom[i].resind];
418 if (ri->nr == j && ri->ic == c)
424 sprintf(buf, "_%d", j);
427 while (parse_int_char(string, &j, &c));
433 static int select_residueindices(char **string, t_atoms *atoms,
435 atom_id *nr, atom_id *index, char *gname)
437 /*this should be similar to select_residuenumbers except select by index (sequential numbering in file)*/
438 /*resind+1 for 1-indexing*/
444 while ((*string)[0] == ' ')
448 if ((*string)[0] == '-')
450 /* Residue number range selection */
453 printf("Error: residue insertion codes can not be used with residue range selection\n");
457 parse_int(string, &up);
459 for (i = 0; i < atoms->nr; i++)
461 ri = &atoms->resinfo[atoms->atom[i].resind];
462 for (j = n1; (j <= up); j++)
464 if (atoms->atom[i].resind+1 == j && (c == ' ' || ri->ic == c))
471 printf("Found %d atom%s with resind.+1 in range %d-%d\n",
472 *nr, (*nr == 1) ? "" : "s", n1, up);
475 sprintf(buf, "r_%d", n1);
479 sprintf(buf, "r_%d-%d", n1, up);
485 /* Individual residue number/insertion code selection */
490 for (i = 0; i < atoms->nr; i++)
492 ri = &atoms->resinfo[atoms->atom[i].resind];
493 if (atoms->atom[i].resind+1 == j && ri->ic == c)
499 sprintf(buf, "_%d", j);
502 while (parse_int_char(string, &j, &c));
509 static gmx_bool atoms_from_residuenumbers(t_atoms *atoms, int group, t_blocka *block,
510 atom_id *nr, atom_id *index, char *gname)
512 int i, j, j0, j1, resnr, nres;
514 j0 = block->index[group];
515 j1 = block->index[group+1];
517 for (j = j0; j < j1; j++)
519 if (block->a[j] >= nres)
521 printf("Index %s contains number>nres (%d>%d)\n",
522 gname, block->a[j]+1, nres);
526 for (i = 0; i < atoms->nr; i++)
528 resnr = atoms->resinfo[atoms->atom[i].resind].nr;
529 for (j = j0; j < j1; j++)
531 if (block->a[j]+1 == resnr)
539 printf("Found %d atom%s in %d residues from group %s\n",
540 *nr, (*nr == 1) ? "" : "s", j1-j0, gname);
544 static gmx_bool comp_name(char *name, char *search)
546 while (name[0] != '\0' && search[0] != '\0')
554 if (search[1] != '\0')
556 printf("WARNING: Currently '*' is only supported at the end of an expression\n");
565 /* Compare a single character */
566 if (( bCase && strncmp(name, search, 1)) ||
567 (!bCase && gmx_strncasecmp(name, search, 1)))
576 return (name[0] == '\0' && (search[0] == '\0' || search[0] == '*'));
579 static int select_chainnames(t_atoms *atoms, int n_names, char **names,
580 atom_id *nr, atom_id *index)
588 for (i = 0; i < atoms->nr; i++)
590 name[0] = atoms->resinfo[atoms->atom[i].resind].chainid;
592 while (j < n_names && !comp_name(name, names[j]))
602 printf("Found %d atom%s with chain identifier%s",
603 *nr, (*nr == 1) ? "" : "s", (n_names == 1) ? "" : "s");
604 for (j = 0; (j < n_names); j++)
606 printf(" %s", names[j]);
613 static int select_atomnames(t_atoms *atoms, int n_names, char **names,
614 atom_id *nr, atom_id *index, gmx_bool bType)
621 for (i = 0; i < atoms->nr; i++)
625 name = *(atoms->atomtype[i]);
629 name = *(atoms->atomname[i]);
632 while (j < n_names && !comp_name(name, names[j]))
642 printf("Found %d atoms with %s%s",
643 *nr, bType ? "type" : "name", (n_names == 1) ? "" : "s");
644 for (j = 0; (j < n_names); j++)
646 printf(" %s", names[j]);
653 static int select_residuenames(t_atoms *atoms, int n_names, char **names,
654 atom_id *nr, atom_id *index)
661 for (i = 0; i < atoms->nr; i++)
663 name = *(atoms->resinfo[atoms->atom[i].resind].name);
665 while (j < n_names && !comp_name(name, names[j]))
675 printf("Found %d atoms with residue name%s", *nr, (n_names == 1) ? "" : "s");
676 for (j = 0; (j < n_names); j++)
678 printf(" %s", names[j]);
685 static void copy2block(int n, atom_id *index, t_blocka *block)
692 srenew(block->index, block->nr+1);
693 block->index[block->nr] = n0+n;
694 srenew(block->a, n0+n);
695 for (i = 0; (i < n); i++)
697 block->a[n0+i] = index[i];
701 static void make_gname(int n, char **names, char *gname)
705 strcpy(gname, names[0]);
706 for (i = 1; i < n; i++)
709 strcat(gname, names[i]);
713 static void copy_group(int g, t_blocka *block, atom_id *nr, atom_id *index)
717 i0 = block->index[g];
718 *nr = block->index[g+1]-i0;
719 for (i = 0; i < *nr; i++)
721 index[i] = block->a[i0+i];
725 static void remove_group(int nr, int nr2, t_blocka *block, char ***gn)
735 for (j = 0; j <= nr2-nr; j++)
737 if ((nr < 0) || (nr >= block->nr))
739 printf("Group %d does not exist\n", nr+j);
743 shift = block->index[nr+1]-block->index[nr];
744 for (i = block->index[nr+1]; i < block->nra; i++)
746 block->a[i-shift] = block->a[i];
749 for (i = nr; i < block->nr; i++)
751 block->index[i] = block->index[i+1]-shift;
753 name = strdup((*gn)[nr]);
755 for (i = nr; i < block->nr-1; i++)
757 (*gn)[i] = (*gn)[i+1];
760 block->nra = block->index[block->nr];
761 printf("Removed group %d '%s'\n", nr+j, name);
767 static void split_group(t_atoms *atoms, int sel_nr, t_blocka *block, char ***gn,
770 char buf[STRLEN], *name;
774 printf("Splitting group %d '%s' into %s\n", sel_nr, (*gn)[sel_nr],
775 bAtom ? "atoms" : "residues");
777 n0 = block->index[sel_nr];
778 n1 = block->index[sel_nr+1];
779 srenew(block->a, block->nra+n1-n0);
780 for (i = n0; i < n1; i++)
783 resind = atoms->atom[a].resind;
784 name = *(atoms->resinfo[resind].name);
785 if (bAtom || (i == n0) || (atoms->atom[block->a[i-1]].resind != resind))
789 block->index[block->nr] = block->nra;
792 srenew(block->index, block->nr+1);
793 srenew(*gn, block->nr);
796 sprintf(buf, "%s_%s_%d", (*gn)[sel_nr], *atoms->atomname[a], a+1);
800 sprintf(buf, "%s_%s_%d", (*gn)[sel_nr], name, atoms->resinfo[resind].nr);
802 (*gn)[block->nr-1] = strdup(buf);
804 block->a[block->nra] = a;
807 block->index[block->nr] = block->nra;
810 static int split_chain(t_atoms *atoms, rvec *x,
811 int sel_nr, t_blocka *block, char ***gn)
815 atom_id i, a, natoms, *start = NULL, *end = NULL, ca_start, ca_end;
822 while (ca_start < natoms)
824 while ((ca_start < natoms) && strcmp(*atoms->atomname[ca_start], "CA"))
828 if (ca_start < natoms)
830 srenew(start, nchain+1);
831 srenew(end, nchain+1);
832 start[nchain] = ca_start;
833 while ((start[nchain] > 0) &&
834 (atoms->atom[start[nchain]-1].resind ==
835 atoms->atom[ca_start].resind))
848 while ((i < natoms) && strcmp(*atoms->atomname[i], "CA"));
851 rvec_sub(x[ca_end], x[i], vec);
854 while ((i < natoms) && (norm(vec) < 0.45));
856 end[nchain] = ca_end;
857 while ((end[nchain]+1 < natoms) &&
858 (atoms->atom[end[nchain]+1].resind == atoms->atom[ca_end].resind))
862 ca_start = end[nchain]+1;
868 printf("Found 1 chain, will not split\n");
872 printf("Found %d chains\n", nchain);
874 for (j = 0; j < nchain; j++)
876 printf("%d:%6d atoms (%d to %d)\n",
877 j+1, end[j]-start[j]+1, start[j]+1, end[j]+1);
882 srenew(block->a, block->nra+block->index[sel_nr+1]-block->index[sel_nr]);
883 for (j = 0; j < nchain; j++)
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] = strdup(buf);
890 for (i = block->index[sel_nr]; i < block->index[sel_nr+1]; i++)
893 if ((a >= start[j]) && (a <= end[j]))
895 block->a[block->nra] = a;
899 block->index[block->nr] = block->nra;
900 if (block->index[block->nr-1] == block->index[block->nr])
902 remove_group(block->nr-1, NOTSET, block, gn);
912 static gmx_bool check_have_atoms(t_atoms *atoms, char *string)
916 printf("Can not process '%s' without atom info, use option -f\n", string);
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)
929 static char **names, *ostring;
930 static gmx_bool bFirst = TRUE;
931 int j, n_names, sel_nr1;
932 atom_id i, nr1, *index1;
934 gmx_bool bRet, bCompl;
939 snew(names, MAXNAMES);
940 for (i = 0; i < MAXNAMES; i++)
942 snew(names[i], NAME_LEN+1);
949 while (*string[0] == ' ')
954 if ((*string)[0] == '!')
958 while (*string[0] == ' ')
970 if (parse_int(string, &sel_nr1) ||
971 parse_string(string, &sel_nr1, block->nr, *gn))
973 if ((sel_nr1 >= 0) && (sel_nr1 < block->nr))
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]);
982 printf("Group %d does not exist\n", sel_nr1);
985 else if ((*string)[0] == 'a')
988 if (check_have_atoms(atoms, ostring))
990 if (parse_int(string, &sel_nr1))
992 bRet = select_atomnumbers(string, atoms, sel_nr1, nr, index, gname);
994 else if (parse_names(string, &n_names, names))
996 bRet = select_atomnames(atoms, n_names, names, nr, index, FALSE);
997 make_gname(n_names, names, gname);
1001 else if ((*string)[0] == 't')
1004 if (check_have_atoms(atoms, ostring) &&
1005 parse_names(string, &n_names, names))
1007 if (atoms->atomtype == NULL)
1009 printf("Need a run input file to select atom types\n");
1013 bRet = select_atomnames(atoms, n_names, names, nr, index, TRUE);
1014 make_gname(n_names, names, gname);
1018 else if (strncmp(*string, "res", 3) == 0)
1021 if (check_have_atoms(atoms, ostring) &&
1022 parse_int(string, &sel_nr1) &&
1023 (sel_nr1 >= 0) && (sel_nr1 < block->nr) )
1025 bRet = atoms_from_residuenumbers(atoms,
1026 sel_nr1, block, nr, index, (*gn)[sel_nr1]);
1027 sprintf(gname, "atom_%s", (*gn)[sel_nr1]);
1030 else if (strncmp(*string, "ri", 2) == 0)
1033 if (check_have_atoms(atoms, ostring) &&
1034 parse_int_char(string, &sel_nr1, &c))
1036 bRet = select_residueindices(string, atoms, sel_nr1, c, nr, index, gname);
1039 else if ((*string)[0] == 'r')
1042 if (check_have_atoms(atoms, ostring))
1044 if (parse_int_char(string, &sel_nr1, &c))
1046 bRet = select_residuenumbers(string, atoms, sel_nr1, c, nr, index, gname);
1048 else if (parse_names(string, &n_names, names))
1050 bRet = select_residuenames(atoms, n_names, names, nr, index);
1051 make_gname(n_names, names, gname);
1055 else if (strncmp(*string, "chain", 5) == 0)
1058 if (check_have_atoms(atoms, ostring) &&
1059 parse_names(string, &n_names, names))
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++)
1065 strcat(gname, names[i]);
1071 snew(index1, natoms-*nr);
1073 for (i = 0; i < natoms; i++)
1076 while ((j < *nr) && (index[j] != i))
1082 if (nr1 >= natoms-*nr)
1084 printf("There are double atoms in your index group\n");
1092 for (i = 0; i < nr1; i++)
1094 index[i] = index1[i];
1098 for (i = strlen(gname)+1; i > 0; i--)
1100 gname[i] = gname[i-1];
1103 printf("Complemented group: %d atoms\n", *nr);
1109 static void list_residues(t_atoms *atoms)
1111 int i, j, start, end, prev_resind, resind;
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++)
1119 resind = atoms->atom[i].resind;
1120 if ((resind != prev_resind) || (i == atoms->nr-1))
1122 if ((bDiff = strcmp(*atoms->resinfo[resind].name,
1123 *atoms->resinfo[start].name)) ||
1136 for (j = start; j <= end; j++)
1139 j+1, *(atoms->resinfo[j].name));
1144 printf(" %4d - %4d %-5s ",
1145 start+1, end+1, *(atoms->resinfo[start].name));
1150 prev_resind = resind;
1155 static void edit_index(int natoms, t_atoms *atoms, rvec *x, t_blocka *block, char ***gn, gmx_bool bVerbose)
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;
1168 snew(atnames, MAXNAMES);
1169 for (i = 0; i < MAXNAMES; i++)
1171 snew(atnames[i], NAME_LEN+1);
1177 snew(index, natoms);
1178 snew(index1, natoms);
1179 snew(index2, natoms);
1186 if (bVerbose || bPrintOnce || newgroup != NOTSET)
1189 if (bVerbose || bPrintOnce || newgroup == NOTSET)
1199 for (i = i0; i < i1; i++)
1201 printf("%3d %-20s: %5d atoms\n", i, (*gn)[i],
1202 block->index[i+1]-block->index[i]);
1206 if (bVerbose || bPrintOnce)
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");
1220 if (NULL == fgets(inp_string, STRLEN, stdin))
1222 gmx_fatal(FARGS, "Error reading user input");
1224 inp_string[strlen(inp_string)-1] = 0;
1226 string = inp_string;
1227 while (string[0] == ' ')
1234 if (string[0] == 'h')
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");
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");
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");
1278 printf("\npress Enter ");
1282 else if (strncmp(string, "del", 3) == 0)
1285 if (parse_int(&string, &sel_nr))
1287 while (string[0] == ' ')
1291 if (string[0] == '-')
1294 parse_int(&string, &sel_nr2);
1300 while (string[0] == ' ')
1304 if (string[0] == '\0')
1306 remove_group(sel_nr, sel_nr2, block, gn);
1310 printf("\nSyntax error: \"%s\"\n", string);
1314 else if (strncmp(string, "keep", 4) == 0)
1317 if (parse_int(&string, &sel_nr))
1319 remove_group(sel_nr+1, block->nr-1, block, gn);
1320 remove_group(0, sel_nr-1, block, gn);
1323 else if (strncmp(string, "name", 4) == 0)
1326 if (parse_int(&string, &sel_nr))
1328 if ((sel_nr >= 0) && (sel_nr < block->nr))
1330 sscanf(string, "%s", gname);
1331 sfree((*gn)[sel_nr]);
1332 (*gn)[sel_nr] = strdup(gname);
1336 else if (strncmp(string, "case", 4) == 0)
1339 printf("Switched to case %s\n", bCase ? "sensitive" : "insensitive");
1341 else if (string[0] == 'v')
1343 bVerbose = !bVerbose;
1344 printf("Turned verbose %s\n", bVerbose ? "on" : "off");
1346 else if (string[0] == 'l')
1348 if (check_have_atoms(atoms, ostring) )
1350 list_residues(atoms);
1353 else if (strncmp(string, "splitch", 7) == 0)
1356 if (check_have_atoms(atoms, ostring) &&
1357 parse_int(&string, &sel_nr) &&
1358 (sel_nr >= 0) && (sel_nr < block->nr))
1360 split_chain(atoms, x, sel_nr, block, gn);
1363 else if (strncmp(string, "splitres", 8) == 0)
1366 if (check_have_atoms(atoms, ostring) &&
1367 parse_int(&string, &sel_nr) &&
1368 (sel_nr >= 0) && (sel_nr < block->nr))
1370 split_group(atoms, sel_nr, block, gn, FALSE);
1373 else if (strncmp(string, "splitat", 7) == 0)
1376 if (check_have_atoms(atoms, ostring) &&
1377 parse_int(&string, &sel_nr) &&
1378 (sel_nr >= 0) && (sel_nr < block->nr))
1380 split_group(atoms, sel_nr, block, gn, TRUE);
1383 else if (string[0] == '\0')
1387 else if (string[0] != 'q')
1391 if (parse_entry(&string, natoms, atoms, block, gn, &nr, index, gname))
1395 while (string[0] == ' ')
1402 if (string[0] == '&')
1406 else if (string[0] == '|')
1415 for (i = 0; i < nr; i++)
1417 index1[i] = index[i];
1419 strcpy(gname1, gname);
1420 if (parse_entry(&string, natoms, atoms, block, gn, &nr2, index2, gname2))
1424 or_groups(nr1, index1, nr2, index2, &nr, index);
1425 sprintf(gname, "%s_%s", gname1, gname2);
1429 and_groups(nr1, index1, nr2, index2, &nr, index);
1430 sprintf(gname, "%s_&_%s", gname1, gname2);
1435 while (bAnd || bOr);
1437 while (string[0] == ' ')
1443 printf("\nSyntax error: \"%s\"\n", string);
1447 copy2block(nr, index, block);
1448 srenew(*gn, block->nr);
1449 newgroup = block->nr-1;
1450 (*gn)[newgroup] = strdup(gname);
1454 printf("Group is empty\n");
1458 while (string[0] != 'q');
1465 static int block2natoms(t_blocka *block)
1470 for (i = 0; i < block->nra; i++)
1472 natoms = max(natoms, block->a[i]);
1479 void merge_blocks(t_blocka *dest, t_blocka *source)
1483 /* count groups, srenew and fill */
1486 dest->nr += source->nr;
1487 srenew(dest->index, dest->nr+1);
1488 for (i = 0; i < source->nr; i++)
1490 dest->index[i0+i] = nra0 + source->index[i];
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++)
1497 dest->a[nra0+i] = source->a[i];
1500 /* terminate list */
1501 dest->index[dest->nr] = dest->nra;
1505 int gmx_make_ndx(int argc, char *argv[])
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",
1516 "With the index editor you can select on atom, residue and chain names",
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."
1527 static int natoms = 0;
1528 static gmx_bool bVerbose = FALSE;
1529 static gmx_bool bDuplicate = FALSE;
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" }
1538 #define NPA asize(pa)
1543 const char *stxfile;
1545 const char *ndxoutfile;
1552 t_blocka *block, *block2;
1553 char **gnames, **gnames2;
1555 { efSTX, "-f", NULL, ffOPTRD },
1556 { efNDX, "-n", NULL, ffOPTRDMULT },
1557 { efNDX, "-o", NULL, ffWRITE }
1559 #define NFILE asize(fnm)
1561 if (!parse_common_args(&argc, argv, 0, NFILE, fnm, NPA, pa, asize(desc), desc,
1567 stxfile = ftp2fn_null(efSTX, NFILE, fnm);
1568 if (opt2bSet("-n", NFILE, fnm))
1570 nndxin = opt2fns(&ndxinfiles, "-n", NFILE, fnm);
1576 ndxoutfile = opt2fn("-o", NFILE, fnm);
1577 bNatoms = opt2parg_bSet("-natoms", NPA, pa);
1579 if (!stxfile && !nndxin)
1581 gmx_fatal(FARGS, "No input files (structure or index)");
1587 get_stx_coordnum(stxfile, &(atoms->nr));
1588 init_t_atoms(atoms, atoms->nr, TRUE);
1591 fprintf(stderr, "\nReading structure file\n");
1592 read_stx_conf(stxfile, title, atoms, x, v, &ePBC, box);
1602 /* read input file(s) */
1603 block = new_blocka();
1605 printf("Going to read %d old index file(s)\n", nndxin);
1608 for (i = 0; i < nndxin; i++)
1610 block2 = init_index(ndxinfiles[i], &gnames2);
1611 srenew(gnames, block->nr+block2->nr);
1612 for (j = 0; j < block2->nr; j++)
1614 gnames[block->nr+j] = gnames2[j];
1617 merge_blocks(block, block2);
1619 sfree(block2->index);
1620 /* done_block(block2); */
1627 analyse(atoms, block, &gnames, FALSE, TRUE);
1632 natoms = block2natoms(block);
1633 printf("Counted atom numbers up to %d in index file\n", natoms);
1636 edit_index(natoms, atoms, x, block, &gnames, bVerbose);
1638 write_index(ndxoutfile, block, gnames, bDuplicate, natoms);