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, 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/fileio/futil.h"
47 #include "gromacs/commandline/pargs.h"
48 #include "gromacs/fileio/confio.h"
58 gmx_bool bCase = FALSE;
60 static int or_groups(atom_id nr1, atom_id *at1, atom_id nr2, atom_id *at2,
61 atom_id *nr, atom_id *at)
63 atom_id i1, i2, max = 0;
69 for (i1 = 0; i1 < nr1; i1++)
71 if ((i1 > 0) && (at1[i1] <= max))
77 for (i1 = 0; i1 < nr2; i1++)
79 if ((i1 > 0) && (at2[i1] <= max))
88 printf("One of your groups is not ascending\n");
95 while ((i1 < nr1) || (i2 < nr2))
97 if ((i2 == nr2) || ((i1 < nr1) && (at1[i1] < at2[i2])))
105 if ((i2 < nr2) && ((i1 == nr1) || (at1[i1] > at2[i2])))
114 printf("Merged two groups with OR: %d %d -> %d\n", nr1, nr2, *nr);
120 static int and_groups(atom_id nr1, atom_id *at1, atom_id nr2, atom_id *at2,
121 atom_id *nr, atom_id *at)
126 for (i1 = 0; i1 < nr1; i1++)
128 for (i2 = 0; i2 < nr2; i2++)
130 if (at1[i1] == at2[i2])
138 printf("Merged two groups with AND: %d %d -> %d\n", nr1, nr2, *nr);
143 static gmx_bool is_name_char(char c)
145 /* This string should contain all characters that can not be
146 * the first letter of a name due to the make_ndx syntax.
148 const char *spec = " !&|";
150 return (c != '\0' && strchr(spec, c) == NULL);
153 static int parse_names(char **string, int *n_names, char **names)
158 while (is_name_char((*string)[0]) || (*string)[0] == ' ')
160 if (is_name_char((*string)[0]))
162 if (*n_names >= MAXNAMES)
164 gmx_fatal(FARGS, "To many names: %d\n", *n_names+1);
167 while (is_name_char((*string)[i]))
169 names[*n_names][i] = (*string)[i];
173 printf("Name is too long, the maximum is %d characters\n", NAME_LEN);
177 names[*n_names][i] = '\0';
180 upstring(names[*n_names]);
194 static gmx_bool parse_int_char(char **string, int *nr, char *c)
201 while ((*string)[0] == ' ')
210 if (isdigit((*string)[0]))
212 *nr = (*string)[0]-'0';
214 while (isdigit((*string)[0]))
216 *nr = (*nr)*10+(*string)[0]-'0';
219 if (isalpha((*string)[0]))
224 /* Check if there is at most one non-digit character */
225 if (!isalnum((*string)[0]))
242 static gmx_bool parse_int(char **string, int *nr)
248 bRet = parse_int_char(string, nr, &c);
249 if (bRet && c != ' ')
258 static gmx_bool isquote(char c)
263 static gmx_bool parse_string(char **string, int *nr, int ngrps, char **grpname)
268 while ((*string)[0] == ' ')
274 if (isquote((*string)[0]))
278 s = strdup((*string));
282 (*string) += sp-s + 1;
284 (*nr) = find_group(s, ngrps, grpname);
288 return (*nr) != NOTSET;
291 static int select_atomnumbers(char **string, t_atoms *atoms, atom_id n1,
292 atom_id *nr, atom_id *index, char *gname)
298 while ((*string)[0] == ' ')
302 if ((*string)[0] == '-')
305 parse_int(string, &up);
306 if ((n1 < 1) || (n1 > atoms->nr) || (up < 1) || (up > atoms->nr))
308 printf("Invalid atom range\n");
312 for (i = n1-1; i <= up-1; i++)
317 printf("Found %d atom%s in range %d-%d\n", *nr, (*nr == 1) ? "" : "s", n1, up);
320 sprintf(buf, "a_%d", n1);
324 sprintf(buf, "a_%d-%d", n1, up);
335 if ((i-1 >= 0) && (i-1 < atoms->nr))
339 sprintf(buf, "_%d", i);
344 printf("Invalid atom number %d\n", i);
348 while ((*nr != 0) && (parse_int(string, &i)));
354 static int select_residuenumbers(char **string, t_atoms *atoms,
356 atom_id *nr, atom_id *index, char *gname)
363 while ((*string)[0] == ' ')
367 if ((*string)[0] == '-')
369 /* Residue number range selection */
372 printf("Error: residue insertion codes can not be used with residue range selection\n");
376 parse_int(string, &up);
378 for (i = 0; i < atoms->nr; i++)
380 ri = &atoms->resinfo[atoms->atom[i].resind];
381 for (j = n1; (j <= up); j++)
383 if (ri->nr == j && (c == ' ' || ri->ic == c))
390 printf("Found %d atom%s with res.nr. in range %d-%d\n",
391 *nr, (*nr == 1) ? "" : "s", n1, up);
394 sprintf(buf, "r_%d", n1);
398 sprintf(buf, "r_%d-%d", n1, up);
404 /* Individual residue number/insertion code selection */
409 for (i = 0; i < atoms->nr; i++)
411 ri = &atoms->resinfo[atoms->atom[i].resind];
412 if (ri->nr == j && ri->ic == c)
418 sprintf(buf, "_%d", j);
421 while (parse_int_char(string, &j, &c));
427 static int select_residueindices(char **string, t_atoms *atoms,
429 atom_id *nr, atom_id *index, char *gname)
431 /*this should be similar to select_residuenumbers except select by index (sequential numbering in file)*/
432 /*resind+1 for 1-indexing*/
438 while ((*string)[0] == ' ')
442 if ((*string)[0] == '-')
444 /* Residue number range selection */
447 printf("Error: residue insertion codes can not be used with residue range selection\n");
451 parse_int(string, &up);
453 for (i = 0; i < atoms->nr; i++)
455 ri = &atoms->resinfo[atoms->atom[i].resind];
456 for (j = n1; (j <= up); j++)
458 if (atoms->atom[i].resind+1 == j && (c == ' ' || ri->ic == c))
465 printf("Found %d atom%s with resind.+1 in range %d-%d\n",
466 *nr, (*nr == 1) ? "" : "s", n1, up);
469 sprintf(buf, "r_%d", n1);
473 sprintf(buf, "r_%d-%d", n1, up);
479 /* Individual residue number/insertion code selection */
484 for (i = 0; i < atoms->nr; i++)
486 ri = &atoms->resinfo[atoms->atom[i].resind];
487 if (atoms->atom[i].resind+1 == j && ri->ic == c)
493 sprintf(buf, "_%d", j);
496 while (parse_int_char(string, &j, &c));
503 static gmx_bool atoms_from_residuenumbers(t_atoms *atoms, int group, t_blocka *block,
504 atom_id *nr, atom_id *index, char *gname)
506 int i, j, j0, j1, resnr, nres;
508 j0 = block->index[group];
509 j1 = block->index[group+1];
511 for (j = j0; j < j1; j++)
513 if (block->a[j] >= nres)
515 printf("Index %s contains number>nres (%d>%d)\n",
516 gname, block->a[j]+1, nres);
520 for (i = 0; i < atoms->nr; i++)
522 resnr = atoms->resinfo[atoms->atom[i].resind].nr;
523 for (j = j0; j < j1; j++)
525 if (block->a[j]+1 == resnr)
533 printf("Found %d atom%s in %d residues from group %s\n",
534 *nr, (*nr == 1) ? "" : "s", j1-j0, gname);
538 static gmx_bool comp_name(char *name, char *search)
540 while (name[0] != '\0' && search[0] != '\0')
548 if (search[1] != '\0')
550 printf("WARNING: Currently '*' is only supported at the end of an expression\n");
559 /* Compare a single character */
560 if (( bCase && strncmp(name, search, 1)) ||
561 (!bCase && gmx_strncasecmp(name, search, 1)))
570 return (name[0] == '\0' && (search[0] == '\0' || search[0] == '*'));
573 static int select_chainnames(t_atoms *atoms, int n_names, char **names,
574 atom_id *nr, atom_id *index)
582 for (i = 0; i < atoms->nr; i++)
584 name[0] = atoms->resinfo[atoms->atom[i].resind].chainid;
586 while (j < n_names && !comp_name(name, names[j]))
596 printf("Found %d atom%s with chain identifier%s",
597 *nr, (*nr == 1) ? "" : "s", (n_names == 1) ? "" : "s");
598 for (j = 0; (j < n_names); j++)
600 printf(" %s", names[j]);
607 static int select_atomnames(t_atoms *atoms, int n_names, char **names,
608 atom_id *nr, atom_id *index, gmx_bool bType)
615 for (i = 0; i < atoms->nr; i++)
619 name = *(atoms->atomtype[i]);
623 name = *(atoms->atomname[i]);
626 while (j < n_names && !comp_name(name, names[j]))
636 printf("Found %d atoms with %s%s",
637 *nr, bType ? "type" : "name", (n_names == 1) ? "" : "s");
638 for (j = 0; (j < n_names); j++)
640 printf(" %s", names[j]);
647 static int select_residuenames(t_atoms *atoms, int n_names, char **names,
648 atom_id *nr, atom_id *index)
655 for (i = 0; i < atoms->nr; i++)
657 name = *(atoms->resinfo[atoms->atom[i].resind].name);
659 while (j < n_names && !comp_name(name, names[j]))
669 printf("Found %d atoms with residue name%s", *nr, (n_names == 1) ? "" : "s");
670 for (j = 0; (j < n_names); j++)
672 printf(" %s", names[j]);
679 static void copy2block(int n, atom_id *index, t_blocka *block)
686 srenew(block->index, block->nr+1);
687 block->index[block->nr] = n0+n;
688 srenew(block->a, n0+n);
689 for (i = 0; (i < n); i++)
691 block->a[n0+i] = index[i];
695 static void make_gname(int n, char **names, char *gname)
699 strcpy(gname, names[0]);
700 for (i = 1; i < n; i++)
703 strcat(gname, names[i]);
707 static void copy_group(int g, t_blocka *block, atom_id *nr, atom_id *index)
711 i0 = block->index[g];
712 *nr = block->index[g+1]-i0;
713 for (i = 0; i < *nr; i++)
715 index[i] = block->a[i0+i];
719 static void remove_group(int nr, int nr2, t_blocka *block, char ***gn)
729 for (j = 0; j <= nr2-nr; j++)
731 if ((nr < 0) || (nr >= block->nr))
733 printf("Group %d does not exist\n", nr+j);
737 shift = block->index[nr+1]-block->index[nr];
738 for (i = block->index[nr+1]; i < block->nra; i++)
740 block->a[i-shift] = block->a[i];
743 for (i = nr; i < block->nr; i++)
745 block->index[i] = block->index[i+1]-shift;
747 name = strdup((*gn)[nr]);
749 for (i = nr; i < block->nr-1; i++)
751 (*gn)[i] = (*gn)[i+1];
754 block->nra = block->index[block->nr];
755 printf("Removed group %d '%s'\n", nr+j, name);
761 static void split_group(t_atoms *atoms, int sel_nr, t_blocka *block, char ***gn,
764 char buf[STRLEN], *name;
768 printf("Splitting group %d '%s' into %s\n", sel_nr, (*gn)[sel_nr],
769 bAtom ? "atoms" : "residues");
771 n0 = block->index[sel_nr];
772 n1 = block->index[sel_nr+1];
773 srenew(block->a, block->nra+n1-n0);
774 for (i = n0; i < n1; i++)
777 resind = atoms->atom[a].resind;
778 name = *(atoms->resinfo[resind].name);
779 if (bAtom || (i == n0) || (atoms->atom[block->a[i-1]].resind != resind))
783 block->index[block->nr] = block->nra;
786 srenew(block->index, block->nr+1);
787 srenew(*gn, block->nr);
790 sprintf(buf, "%s_%s_%d", (*gn)[sel_nr], *atoms->atomname[a], a+1);
794 sprintf(buf, "%s_%s_%d", (*gn)[sel_nr], name, atoms->resinfo[resind].nr);
796 (*gn)[block->nr-1] = strdup(buf);
798 block->a[block->nra] = a;
801 block->index[block->nr] = block->nra;
804 static int split_chain(t_atoms *atoms, rvec *x,
805 int sel_nr, t_blocka *block, char ***gn)
809 atom_id i, a, natoms, *start = NULL, *end = NULL, ca_start, ca_end;
816 while (ca_start < natoms)
818 while ((ca_start < natoms) && strcmp(*atoms->atomname[ca_start], "CA"))
822 if (ca_start < natoms)
824 srenew(start, nchain+1);
825 srenew(end, nchain+1);
826 start[nchain] = ca_start;
827 while ((start[nchain] > 0) &&
828 (atoms->atom[start[nchain]-1].resind ==
829 atoms->atom[ca_start].resind))
842 while ((i < natoms) && strcmp(*atoms->atomname[i], "CA"));
845 rvec_sub(x[ca_end], x[i], vec);
848 while ((i < natoms) && (norm(vec) < 0.45));
850 end[nchain] = ca_end;
851 while ((end[nchain]+1 < natoms) &&
852 (atoms->atom[end[nchain]+1].resind == atoms->atom[ca_end].resind))
856 ca_start = end[nchain]+1;
862 printf("Found 1 chain, will not split\n");
866 printf("Found %d chains\n", nchain);
868 for (j = 0; j < nchain; j++)
870 printf("%d:%6d atoms (%d to %d)\n",
871 j+1, end[j]-start[j]+1, start[j]+1, end[j]+1);
876 srenew(block->a, block->nra+block->index[sel_nr+1]-block->index[sel_nr]);
877 for (j = 0; j < nchain; j++)
880 srenew(block->index, block->nr+1);
881 srenew(*gn, block->nr);
882 sprintf(buf, "%s_chain%d", (*gn)[sel_nr], j+1);
883 (*gn)[block->nr-1] = strdup(buf);
884 for (i = block->index[sel_nr]; i < block->index[sel_nr+1]; i++)
887 if ((a >= start[j]) && (a <= end[j]))
889 block->a[block->nra] = a;
893 block->index[block->nr] = block->nra;
894 if (block->index[block->nr-1] == block->index[block->nr])
896 remove_group(block->nr-1, NOTSET, block, gn);
906 static gmx_bool check_have_atoms(t_atoms *atoms, char *string)
910 printf("Can not process '%s' without atom info, use option -f\n", string);
919 static gmx_bool parse_entry(char **string, int natoms, t_atoms *atoms,
920 t_blocka *block, char ***gn,
921 atom_id *nr, atom_id *index, char *gname)
923 static char **names, *ostring;
924 static gmx_bool bFirst = TRUE;
925 int j, n_names, sel_nr1;
926 atom_id i, nr1, *index1;
928 gmx_bool bRet, bCompl;
933 snew(names, MAXNAMES);
934 for (i = 0; i < MAXNAMES; i++)
936 snew(names[i], NAME_LEN+1);
943 while (*string[0] == ' ')
948 if ((*string)[0] == '!')
952 while (*string[0] == ' ')
964 if (parse_int(string, &sel_nr1) ||
965 parse_string(string, &sel_nr1, block->nr, *gn))
967 if ((sel_nr1 >= 0) && (sel_nr1 < block->nr))
969 copy_group(sel_nr1, block, nr, index);
970 strcpy(gname, (*gn)[sel_nr1]);
971 printf("Copied index group %d '%s'\n", sel_nr1, (*gn)[sel_nr1]);
976 printf("Group %d does not exist\n", sel_nr1);
979 else if ((*string)[0] == 'a')
982 if (check_have_atoms(atoms, ostring))
984 if (parse_int(string, &sel_nr1))
986 bRet = select_atomnumbers(string, atoms, sel_nr1, nr, index, gname);
988 else if (parse_names(string, &n_names, names))
990 bRet = select_atomnames(atoms, n_names, names, nr, index, FALSE);
991 make_gname(n_names, names, gname);
995 else if ((*string)[0] == 't')
998 if (check_have_atoms(atoms, ostring) &&
999 parse_names(string, &n_names, names))
1001 if (atoms->atomtype == NULL)
1003 printf("Need a run input file to select atom types\n");
1007 bRet = select_atomnames(atoms, n_names, names, nr, index, TRUE);
1008 make_gname(n_names, names, gname);
1012 else if (strncmp(*string, "res", 3) == 0)
1015 if (check_have_atoms(atoms, ostring) &&
1016 parse_int(string, &sel_nr1) &&
1017 (sel_nr1 >= 0) && (sel_nr1 < block->nr) )
1019 bRet = atoms_from_residuenumbers(atoms,
1020 sel_nr1, block, nr, index, (*gn)[sel_nr1]);
1021 sprintf(gname, "atom_%s", (*gn)[sel_nr1]);
1024 else if (strncmp(*string, "ri", 2) == 0)
1027 if (check_have_atoms(atoms, ostring) &&
1028 parse_int_char(string, &sel_nr1, &c))
1030 bRet = select_residueindices(string, atoms, sel_nr1, c, nr, index, gname);
1033 else if ((*string)[0] == 'r')
1036 if (check_have_atoms(atoms, ostring))
1038 if (parse_int_char(string, &sel_nr1, &c))
1040 bRet = select_residuenumbers(string, atoms, sel_nr1, c, nr, index, gname);
1042 else if (parse_names(string, &n_names, names))
1044 bRet = select_residuenames(atoms, n_names, names, nr, index);
1045 make_gname(n_names, names, gname);
1049 else if (strncmp(*string, "chain", 5) == 0)
1052 if (check_have_atoms(atoms, ostring) &&
1053 parse_names(string, &n_names, names))
1055 bRet = select_chainnames(atoms, n_names, names, nr, index);
1056 sprintf(gname, "ch%s", names[0]);
1057 for (i = 1; i < n_names; i++)
1059 strcat(gname, names[i]);
1065 snew(index1, natoms-*nr);
1067 for (i = 0; i < natoms; i++)
1070 while ((j < *nr) && (index[j] != i))
1076 if (nr1 >= natoms-*nr)
1078 printf("There are double atoms in your index group\n");
1086 for (i = 0; i < nr1; i++)
1088 index[i] = index1[i];
1092 for (i = strlen(gname)+1; i > 0; i--)
1094 gname[i] = gname[i-1];
1097 printf("Complemented group: %d atoms\n", *nr);
1103 static void list_residues(t_atoms *atoms)
1105 int i, j, start, end, prev_resind, resind;
1108 /* Print all the residues, assuming continuous resnr count */
1109 start = atoms->atom[0].resind;
1110 prev_resind = start;
1111 for (i = 0; i < atoms->nr; i++)
1113 resind = atoms->atom[i].resind;
1114 if ((resind != prev_resind) || (i == atoms->nr-1))
1116 if ((bDiff = strcmp(*atoms->resinfo[resind].name,
1117 *atoms->resinfo[start].name)) ||
1130 for (j = start; j <= end; j++)
1133 j+1, *(atoms->resinfo[j].name));
1138 printf(" %4d - %4d %-5s ",
1139 start+1, end+1, *(atoms->resinfo[start].name));
1144 prev_resind = resind;
1149 static void edit_index(int natoms, t_atoms *atoms, rvec *x, t_blocka *block, char ***gn, gmx_bool bVerbose)
1151 static char **atnames, *ostring;
1152 static gmx_bool bFirst = TRUE;
1153 char inp_string[STRLEN], *string;
1154 char gname[STRLEN], gname1[STRLEN], gname2[STRLEN];
1155 int i, i0, i1, sel_nr, sel_nr2, newgroup;
1156 atom_id nr, nr1, nr2, *index, *index1, *index2;
1157 gmx_bool bAnd, bOr, bPrintOnce;
1162 snew(atnames, MAXNAMES);
1163 for (i = 0; i < MAXNAMES; i++)
1165 snew(atnames[i], NAME_LEN+1);
1171 snew(index, natoms);
1172 snew(index1, natoms);
1173 snew(index2, natoms);
1180 if (bVerbose || bPrintOnce || newgroup != NOTSET)
1183 if (bVerbose || bPrintOnce || newgroup == NOTSET)
1193 for (i = i0; i < i1; i++)
1195 printf("%3d %-20s: %5d atoms\n", i, (*gn)[i],
1196 block->index[i+1]-block->index[i]);
1200 if (bVerbose || bPrintOnce)
1203 printf(" nr : group ! 'name' nr name 'splitch' nr Enter: list groups\n");
1204 printf(" 'a': atom & 'del' nr 'splitres' nr 'l': list residues\n");
1205 printf(" 't': atom type | 'keep' nr 'splitat' nr 'h': help\n");
1206 printf(" 'r': residue 'res' nr 'chain' char\n");
1207 printf(" \"name\": group 'case': case %s 'q': save and quit\n",
1208 bCase ? "insensitive" : "sensitive ");
1209 printf(" 'ri': residue index\n");
1214 if (NULL == fgets(inp_string, STRLEN, stdin))
1216 gmx_fatal(FARGS, "Error reading user input");
1218 inp_string[strlen(inp_string)-1] = 0;
1220 string = inp_string;
1221 while (string[0] == ' ')
1228 if (string[0] == 'h')
1230 printf(" nr : selects an index group by number or quoted string.\n");
1231 printf(" The string is first matched against the whole group name,\n");
1232 printf(" then against the beginning and finally against an\n");
1233 printf(" arbitrary substring. A multiple match is an error.\n");
1235 printf(" 'a' nr1 [nr2 ...] : selects atoms, atom numbering starts at 1.\n");
1236 printf(" 'a' nr1 - nr2 : selects atoms in the range from nr1 to nr2.\n");
1237 printf(" 'a' name1[*] [name2[*] ...] : selects atoms by name(s), '?' matches any char,\n");
1238 printf(" wildcard '*' allowed at the end of a name.\n");
1239 printf(" 't' type1[*] [type2[*] ...] : as 'a', but for type, run input file required.\n");
1240 printf(" 'r' nr1[ic1] [nr2[ic2] ...] : selects residues by number and insertion code.\n");
1241 printf(" 'r' nr1 - nr2 : selects residues in the range from nr1 to nr2.\n");
1242 printf(" 'r' name1[*] [name2[*] ...] : as 'a', but for residue names.\n");
1243 printf(" 'ri' nr1 - nr2 : selects residue indices, 1-indexed, (as opposed to numbers) in the range from nr1 to nr2.\n");
1244 printf(" 'chain' ch1 [ch2 ...] : selects atoms by chain identifier(s),\n");
1245 printf(" not available with a .gro file as input.\n");
1246 printf(" ! : takes the complement of a group with respect to all\n");
1247 printf(" the atoms in the input file.\n");
1248 printf(" & | : AND and OR, can be placed between any of the options\n");
1249 printf(" above, the input is processed from left to right.\n");
1250 printf(" 'name' nr name : rename group nr to name.\n");
1251 printf(" 'del' nr1 [- nr2] : deletes one group or groups in the range from nr1 to nr2.\n");
1252 printf(" 'keep' nr : deletes all groups except nr.\n");
1253 printf(" 'case' : make all name compares case (in)sensitive.\n");
1254 printf(" 'splitch' nr : split group into chains using CA distances.\n");
1255 printf(" 'splitres' nr : split group into residues.\n");
1256 printf(" 'splitat' nr : split group into atoms.\n");
1257 printf(" 'res' nr : interpret numbers in group as residue numbers\n");
1258 printf(" Enter : list the currently defined groups and commands\n");
1259 printf(" 'l' : list the residues.\n");
1260 printf(" 'h' : show this help.\n");
1261 printf(" 'q' : save and quit.\n");
1263 printf(" Examples:\n");
1264 printf(" > 2 | 4 & r 3-5\n");
1265 printf(" selects all atoms from group 2 and 4 that have residue numbers 3, 4 or 5\n");
1266 printf(" > a C* & !a C CA\n");
1267 printf(" selects all atoms starting with 'C' but not the atoms 'C' and 'CA'\n");
1268 printf(" > \"protein\" & ! \"backb\"\n");
1269 printf(" selects all atoms that are in group 'protein' and not in group 'backbone'\n");
1272 printf("\npress Enter ");
1276 else if (strncmp(string, "del", 3) == 0)
1279 if (parse_int(&string, &sel_nr))
1281 while (string[0] == ' ')
1285 if (string[0] == '-')
1288 parse_int(&string, &sel_nr2);
1294 while (string[0] == ' ')
1298 if (string[0] == '\0')
1300 remove_group(sel_nr, sel_nr2, block, gn);
1304 printf("\nSyntax error: \"%s\"\n", string);
1308 else if (strncmp(string, "keep", 4) == 0)
1311 if (parse_int(&string, &sel_nr))
1313 remove_group(sel_nr+1, block->nr-1, block, gn);
1314 remove_group(0, sel_nr-1, block, gn);
1317 else if (strncmp(string, "name", 4) == 0)
1320 if (parse_int(&string, &sel_nr))
1322 if ((sel_nr >= 0) && (sel_nr < block->nr))
1324 sscanf(string, "%s", gname);
1325 sfree((*gn)[sel_nr]);
1326 (*gn)[sel_nr] = strdup(gname);
1330 else if (strncmp(string, "case", 4) == 0)
1333 printf("Switched to case %s\n", bCase ? "sensitive" : "insensitive");
1335 else if (string[0] == 'v')
1337 bVerbose = !bVerbose;
1338 printf("Turned verbose %s\n", bVerbose ? "on" : "off");
1340 else if (string[0] == 'l')
1342 if (check_have_atoms(atoms, ostring) )
1344 list_residues(atoms);
1347 else if (strncmp(string, "splitch", 7) == 0)
1350 if (check_have_atoms(atoms, ostring) &&
1351 parse_int(&string, &sel_nr) &&
1352 (sel_nr >= 0) && (sel_nr < block->nr))
1354 split_chain(atoms, x, sel_nr, block, gn);
1357 else if (strncmp(string, "splitres", 8) == 0)
1360 if (check_have_atoms(atoms, ostring) &&
1361 parse_int(&string, &sel_nr) &&
1362 (sel_nr >= 0) && (sel_nr < block->nr))
1364 split_group(atoms, sel_nr, block, gn, FALSE);
1367 else if (strncmp(string, "splitat", 7) == 0)
1370 if (check_have_atoms(atoms, ostring) &&
1371 parse_int(&string, &sel_nr) &&
1372 (sel_nr >= 0) && (sel_nr < block->nr))
1374 split_group(atoms, sel_nr, block, gn, TRUE);
1377 else if (string[0] == '\0')
1381 else if (string[0] != 'q')
1385 if (parse_entry(&string, natoms, atoms, block, gn, &nr, index, gname))
1389 while (string[0] == ' ')
1396 if (string[0] == '&')
1400 else if (string[0] == '|')
1409 for (i = 0; i < nr; i++)
1411 index1[i] = index[i];
1413 strcpy(gname1, gname);
1414 if (parse_entry(&string, natoms, atoms, block, gn, &nr2, index2, gname2))
1418 or_groups(nr1, index1, nr2, index2, &nr, index);
1419 sprintf(gname, "%s_%s", gname1, gname2);
1423 and_groups(nr1, index1, nr2, index2, &nr, index);
1424 sprintf(gname, "%s_&_%s", gname1, gname2);
1429 while (bAnd || bOr);
1431 while (string[0] == ' ')
1437 printf("\nSyntax error: \"%s\"\n", string);
1441 copy2block(nr, index, block);
1442 srenew(*gn, block->nr);
1443 newgroup = block->nr-1;
1444 (*gn)[newgroup] = strdup(gname);
1448 printf("Group is empty\n");
1452 while (string[0] != 'q');
1459 static int block2natoms(t_blocka *block)
1464 for (i = 0; i < block->nra; i++)
1466 natoms = max(natoms, block->a[i]);
1473 void merge_blocks(t_blocka *dest, t_blocka *source)
1477 /* count groups, srenew and fill */
1480 dest->nr += source->nr;
1481 srenew(dest->index, dest->nr+1);
1482 for (i = 0; i < source->nr; i++)
1484 dest->index[i0+i] = nra0 + source->index[i];
1486 /* count atoms, srenew and fill */
1487 dest->nra += source->nra;
1488 srenew(dest->a, dest->nra);
1489 for (i = 0; i < source->nra; i++)
1491 dest->a[nra0+i] = source->a[i];
1494 /* terminate list */
1495 dest->index[dest->nr] = dest->nra;
1499 int gmx_make_ndx(int argc, char *argv[])
1501 const char *desc[] = {
1502 "Index groups are necessary for almost every GROMACS program.",
1503 "All these programs can generate default index groups. You ONLY",
1504 "have to use [THISMODULE] when you need SPECIAL index groups.",
1505 "There is a default index group for the whole system, 9 default",
1506 "index groups for proteins, and a default index group",
1507 "is generated for every other residue name.[PAR]",
1508 "When no index file is supplied, also [THISMODULE] will generate the",
1510 "With the index editor you can select on atom, residue and chain names",
1512 "When a run input file is supplied you can also select on atom type.",
1513 "You can use NOT, AND and OR, you can split groups",
1514 "into chains, residues or atoms. You can delete and rename groups.[PAR]",
1515 "The atom numbering in the editor and the index file starts at 1."
1518 static int natoms = 0;
1519 static gmx_bool bVerbose = FALSE;
1521 { "-natoms", FALSE, etINT, {&natoms},
1522 "set number of atoms (default: read from coordinate or index file)" },
1523 { "-verbose", FALSE, etBOOL, {&bVerbose},
1524 "HIDDENVerbose output" }
1526 #define NPA asize(pa)
1531 const char *stxfile;
1533 const char *ndxoutfile;
1540 t_blocka *block, *block2;
1541 char **gnames, **gnames2;
1543 { efSTX, "-f", NULL, ffOPTRD },
1544 { efNDX, "-n", NULL, ffOPTRDMULT },
1545 { efNDX, "-o", NULL, ffWRITE }
1547 #define NFILE asize(fnm)
1549 if (!parse_common_args(&argc, argv, 0, NFILE, fnm, NPA, pa, asize(desc), desc,
1555 stxfile = ftp2fn_null(efSTX, NFILE, fnm);
1556 if (opt2bSet("-n", NFILE, fnm))
1558 nndxin = opt2fns(&ndxinfiles, "-n", NFILE, fnm);
1564 ndxoutfile = opt2fn("-o", NFILE, fnm);
1565 bNatoms = opt2parg_bSet("-natoms", NPA, pa);
1567 if (!stxfile && !nndxin)
1569 gmx_fatal(FARGS, "No input files (structure or index)");
1575 get_stx_coordnum(stxfile, &(atoms->nr));
1576 init_t_atoms(atoms, atoms->nr, TRUE);
1579 fprintf(stderr, "\nReading structure file\n");
1580 read_stx_conf(stxfile, title, atoms, x, v, &ePBC, box);
1590 /* read input file(s) */
1591 block = new_blocka();
1593 printf("Going to read %d old index file(s)\n", nndxin);
1596 for (i = 0; i < nndxin; i++)
1598 block2 = init_index(ndxinfiles[i], &gnames2);
1599 srenew(gnames, block->nr+block2->nr);
1600 for (j = 0; j < block2->nr; j++)
1602 gnames[block->nr+j] = gnames2[j];
1605 merge_blocks(block, block2);
1607 sfree(block2->index);
1608 /* done_block(block2); */
1615 analyse(atoms, block, &gnames, FALSE, TRUE);
1620 natoms = block2natoms(block);
1621 printf("Counted atom numbers up to %d in index file\n", natoms);
1624 edit_index(natoms, atoms, x, block, &gnames, bVerbose);
1626 write_index(ndxoutfile, block, gnames);