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.
42 #include "gromacs/utility/futil.h"
43 #include "gromacs/legacyheaders/macros.h"
44 #include "gromacs/fileio/confio.h"
45 #include "gromacs/legacyheaders/typedefs.h"
46 #include "gromacs/topology/index.h"
47 #include "gromacs/math/vec.h"
49 #include "gromacs/commandline/pargs.h"
50 #include "gromacs/topology/block.h"
51 #include "gromacs/utility/cstringutil.h"
52 #include "gromacs/utility/fatalerror.h"
53 #include "gromacs/utility/smalloc.h"
55 /* It's not nice to have size limits, but we should not spend more time
56 * on this ancient tool, but instead use the new selection library.
61 gmx_bool bCase = FALSE;
63 static int or_groups(atom_id nr1, atom_id *at1, atom_id nr2, atom_id *at2,
64 atom_id *nr, atom_id *at)
66 atom_id i1, i2, max = 0;
72 for (i1 = 0; i1 < nr1; i1++)
74 if ((i1 > 0) && (at1[i1] <= max))
80 for (i1 = 0; i1 < nr2; i1++)
82 if ((i1 > 0) && (at2[i1] <= max))
91 printf("One of your groups is not ascending\n");
98 while ((i1 < nr1) || (i2 < nr2))
100 if ((i2 == nr2) || ((i1 < nr1) && (at1[i1] < at2[i2])))
108 if ((i2 < nr2) && ((i1 == nr1) || (at1[i1] > at2[i2])))
117 printf("Merged two groups with OR: %d %d -> %d\n", nr1, nr2, *nr);
123 static int and_groups(atom_id nr1, atom_id *at1, atom_id nr2, atom_id *at2,
124 atom_id *nr, atom_id *at)
129 for (i1 = 0; i1 < nr1; i1++)
131 for (i2 = 0; i2 < nr2; i2++)
133 if (at1[i1] == at2[i2])
141 printf("Merged two groups with AND: %d %d -> %d\n", nr1, nr2, *nr);
146 static gmx_bool is_name_char(char c)
148 /* This string should contain all characters that can not be
149 * the first letter of a name due to the make_ndx syntax.
151 const char *spec = " !&|";
153 return (c != '\0' && strchr(spec, c) == NULL);
156 static int parse_names(char **string, int *n_names, char **names)
161 while (is_name_char((*string)[0]) || (*string)[0] == ' ')
163 if (is_name_char((*string)[0]))
165 if (*n_names >= MAXNAMES)
167 gmx_fatal(FARGS, "To many names: %d\n", *n_names+1);
170 while (is_name_char((*string)[i]))
172 names[*n_names][i] = (*string)[i];
176 printf("Name is too long, the maximum is %d characters\n", NAME_LEN);
180 names[*n_names][i] = '\0';
183 upstring(names[*n_names]);
197 static gmx_bool parse_int_char(char **string, int *nr, char *c)
204 while ((*string)[0] == ' ')
213 if (isdigit((*string)[0]))
215 *nr = (*string)[0]-'0';
217 while (isdigit((*string)[0]))
219 *nr = (*nr)*10+(*string)[0]-'0';
222 if (isalpha((*string)[0]))
227 /* Check if there is at most one non-digit character */
228 if (!isalnum((*string)[0]))
245 static gmx_bool parse_int(char **string, int *nr)
251 bRet = parse_int_char(string, nr, &c);
252 if (bRet && c != ' ')
261 static gmx_bool isquote(char c)
266 static gmx_bool parse_string(char **string, int *nr, int ngrps, char **grpname)
271 while ((*string)[0] == ' ')
277 if (isquote((*string)[0]))
281 s = gmx_strdup((*string));
285 (*string) += sp-s + 1;
287 (*nr) = find_group(s, ngrps, grpname);
291 return (*nr) != NOTSET;
294 static int select_atomnumbers(char **string, t_atoms *atoms, atom_id n1,
295 atom_id *nr, atom_id *index, char *gname)
301 while ((*string)[0] == ' ')
305 if ((*string)[0] == '-')
308 parse_int(string, &up);
309 if ((n1 < 1) || (n1 > atoms->nr) || (up < 1) || (up > atoms->nr))
311 printf("Invalid atom range\n");
315 for (i = n1-1; i <= up-1; i++)
320 printf("Found %d atom%s in range %d-%d\n", *nr, (*nr == 1) ? "" : "s", n1, up);
323 sprintf(buf, "a_%d", n1);
327 sprintf(buf, "a_%d-%d", n1, up);
338 if ((i-1 >= 0) && (i-1 < atoms->nr))
342 sprintf(buf, "_%d", i);
347 printf("Invalid atom number %d\n", i);
351 while ((*nr != 0) && (parse_int(string, &i)));
357 static int select_residuenumbers(char **string, t_atoms *atoms,
359 atom_id *nr, atom_id *index, char *gname)
366 while ((*string)[0] == ' ')
370 if ((*string)[0] == '-')
372 /* Residue number range selection */
375 printf("Error: residue insertion codes can not be used with residue range selection\n");
379 parse_int(string, &up);
381 for (i = 0; i < atoms->nr; i++)
383 ri = &atoms->resinfo[atoms->atom[i].resind];
384 for (j = n1; (j <= up); j++)
386 if (ri->nr == j && (c == ' ' || ri->ic == c))
393 printf("Found %d atom%s with res.nr. in range %d-%d\n",
394 *nr, (*nr == 1) ? "" : "s", n1, up);
397 sprintf(buf, "r_%d", n1);
401 sprintf(buf, "r_%d-%d", n1, up);
407 /* Individual residue number/insertion code selection */
412 for (i = 0; i < atoms->nr; i++)
414 ri = &atoms->resinfo[atoms->atom[i].resind];
415 if (ri->nr == j && ri->ic == c)
421 sprintf(buf, "_%d", j);
424 while (parse_int_char(string, &j, &c));
430 static int select_residueindices(char **string, t_atoms *atoms,
432 atom_id *nr, atom_id *index, char *gname)
434 /*this should be similar to select_residuenumbers except select by index (sequential numbering in file)*/
435 /*resind+1 for 1-indexing*/
441 while ((*string)[0] == ' ')
445 if ((*string)[0] == '-')
447 /* Residue number range selection */
450 printf("Error: residue insertion codes can not be used with residue range selection\n");
454 parse_int(string, &up);
456 for (i = 0; i < atoms->nr; i++)
458 ri = &atoms->resinfo[atoms->atom[i].resind];
459 for (j = n1; (j <= up); j++)
461 if (atoms->atom[i].resind+1 == j && (c == ' ' || ri->ic == c))
468 printf("Found %d atom%s with resind.+1 in range %d-%d\n",
469 *nr, (*nr == 1) ? "" : "s", n1, up);
472 sprintf(buf, "r_%d", n1);
476 sprintf(buf, "r_%d-%d", n1, up);
482 /* Individual residue number/insertion code selection */
487 for (i = 0; i < atoms->nr; i++)
489 ri = &atoms->resinfo[atoms->atom[i].resind];
490 if (atoms->atom[i].resind+1 == j && ri->ic == c)
496 sprintf(buf, "_%d", j);
499 while (parse_int_char(string, &j, &c));
506 static gmx_bool atoms_from_residuenumbers(t_atoms *atoms, int group, t_blocka *block,
507 atom_id *nr, atom_id *index, char *gname)
509 int i, j, j0, j1, resnr, nres;
511 j0 = block->index[group];
512 j1 = block->index[group+1];
514 for (j = j0; j < j1; j++)
516 if (block->a[j] >= nres)
518 printf("Index %s contains number>nres (%d>%d)\n",
519 gname, block->a[j]+1, nres);
523 for (i = 0; i < atoms->nr; i++)
525 resnr = atoms->resinfo[atoms->atom[i].resind].nr;
526 for (j = j0; j < j1; j++)
528 if (block->a[j]+1 == resnr)
536 printf("Found %d atom%s in %d residues from group %s\n",
537 *nr, (*nr == 1) ? "" : "s", j1-j0, gname);
541 static gmx_bool comp_name(char *name, char *search)
543 while (name[0] != '\0' && search[0] != '\0')
551 if (search[1] != '\0')
553 printf("WARNING: Currently '*' is only supported at the end of an expression\n");
562 /* Compare a single character */
563 if (( bCase && strncmp(name, search, 1)) ||
564 (!bCase && gmx_strncasecmp(name, search, 1)))
573 return (name[0] == '\0' && (search[0] == '\0' || search[0] == '*'));
576 static int select_chainnames(t_atoms *atoms, int n_names, char **names,
577 atom_id *nr, atom_id *index)
585 for (i = 0; i < atoms->nr; i++)
587 name[0] = atoms->resinfo[atoms->atom[i].resind].chainid;
589 while (j < n_names && !comp_name(name, names[j]))
599 printf("Found %d atom%s with chain identifier%s",
600 *nr, (*nr == 1) ? "" : "s", (n_names == 1) ? "" : "s");
601 for (j = 0; (j < n_names); j++)
603 printf(" %s", names[j]);
610 static int select_atomnames(t_atoms *atoms, int n_names, char **names,
611 atom_id *nr, atom_id *index, gmx_bool bType)
618 for (i = 0; i < atoms->nr; i++)
622 name = *(atoms->atomtype[i]);
626 name = *(atoms->atomname[i]);
629 while (j < n_names && !comp_name(name, names[j]))
639 printf("Found %d atoms with %s%s",
640 *nr, bType ? "type" : "name", (n_names == 1) ? "" : "s");
641 for (j = 0; (j < n_names); j++)
643 printf(" %s", names[j]);
650 static int select_residuenames(t_atoms *atoms, int n_names, char **names,
651 atom_id *nr, atom_id *index)
658 for (i = 0; i < atoms->nr; i++)
660 name = *(atoms->resinfo[atoms->atom[i].resind].name);
662 while (j < n_names && !comp_name(name, names[j]))
672 printf("Found %d atoms with residue name%s", *nr, (n_names == 1) ? "" : "s");
673 for (j = 0; (j < n_names); j++)
675 printf(" %s", names[j]);
682 static void copy2block(int n, atom_id *index, t_blocka *block)
689 srenew(block->index, block->nr+1);
690 block->index[block->nr] = n0+n;
691 srenew(block->a, n0+n);
692 for (i = 0; (i < n); i++)
694 block->a[n0+i] = index[i];
698 static void make_gname(int n, char **names, char *gname)
702 strcpy(gname, names[0]);
703 for (i = 1; i < n; i++)
706 strcat(gname, names[i]);
710 static void copy_group(int g, t_blocka *block, atom_id *nr, atom_id *index)
714 i0 = block->index[g];
715 *nr = block->index[g+1]-i0;
716 for (i = 0; i < *nr; i++)
718 index[i] = block->a[i0+i];
722 static void remove_group(int nr, int nr2, t_blocka *block, char ***gn)
732 for (j = 0; j <= nr2-nr; j++)
734 if ((nr < 0) || (nr >= block->nr))
736 printf("Group %d does not exist\n", nr+j);
740 shift = block->index[nr+1]-block->index[nr];
741 for (i = block->index[nr+1]; i < block->nra; i++)
743 block->a[i-shift] = block->a[i];
746 for (i = nr; i < block->nr; i++)
748 block->index[i] = block->index[i+1]-shift;
750 name = gmx_strdup((*gn)[nr]);
752 for (i = nr; i < block->nr-1; i++)
754 (*gn)[i] = (*gn)[i+1];
757 block->nra = block->index[block->nr];
758 printf("Removed group %d '%s'\n", nr+j, name);
764 static void split_group(t_atoms *atoms, int sel_nr, t_blocka *block, char ***gn,
767 char buf[STRLEN], *name;
771 printf("Splitting group %d '%s' into %s\n", sel_nr, (*gn)[sel_nr],
772 bAtom ? "atoms" : "residues");
774 n0 = block->index[sel_nr];
775 n1 = block->index[sel_nr+1];
776 srenew(block->a, block->nra+n1-n0);
777 for (i = n0; i < n1; i++)
780 resind = atoms->atom[a].resind;
781 name = *(atoms->resinfo[resind].name);
782 if (bAtom || (i == n0) || (atoms->atom[block->a[i-1]].resind != resind))
786 block->index[block->nr] = block->nra;
789 srenew(block->index, block->nr+1);
790 srenew(*gn, block->nr);
793 sprintf(buf, "%s_%s_%d", (*gn)[sel_nr], *atoms->atomname[a], a+1);
797 sprintf(buf, "%s_%s_%d", (*gn)[sel_nr], name, atoms->resinfo[resind].nr);
799 (*gn)[block->nr-1] = gmx_strdup(buf);
801 block->a[block->nra] = a;
804 block->index[block->nr] = block->nra;
807 static int split_chain(t_atoms *atoms, rvec *x,
808 int sel_nr, t_blocka *block, char ***gn)
812 atom_id i, a, natoms, *start = NULL, *end = NULL, ca_start, ca_end;
819 while (ca_start < natoms)
821 while ((ca_start < natoms) && strcmp(*atoms->atomname[ca_start], "CA"))
825 if (ca_start < natoms)
827 srenew(start, nchain+1);
828 srenew(end, nchain+1);
829 start[nchain] = ca_start;
830 while ((start[nchain] > 0) &&
831 (atoms->atom[start[nchain]-1].resind ==
832 atoms->atom[ca_start].resind))
845 while ((i < natoms) && strcmp(*atoms->atomname[i], "CA"));
848 rvec_sub(x[ca_end], x[i], vec);
855 while (norm(vec) < 0.45);
857 end[nchain] = ca_end;
858 while ((end[nchain]+1 < natoms) &&
859 (atoms->atom[end[nchain]+1].resind == atoms->atom[ca_end].resind))
863 ca_start = end[nchain]+1;
869 printf("Found 1 chain, will not split\n");
873 printf("Found %d chains\n", nchain);
875 for (j = 0; j < nchain; j++)
877 printf("%d:%6d atoms (%d to %d)\n",
878 j+1, end[j]-start[j]+1, start[j]+1, end[j]+1);
883 srenew(block->a, block->nra+block->index[sel_nr+1]-block->index[sel_nr]);
884 for (j = 0; j < nchain; j++)
887 srenew(block->index, block->nr+1);
888 srenew(*gn, block->nr);
889 sprintf(buf, "%s_chain%d", (*gn)[sel_nr], j+1);
890 (*gn)[block->nr-1] = gmx_strdup(buf);
891 for (i = block->index[sel_nr]; i < block->index[sel_nr+1]; i++)
894 if ((a >= start[j]) && (a <= end[j]))
896 block->a[block->nra] = a;
900 block->index[block->nr] = block->nra;
901 if (block->index[block->nr-1] == block->index[block->nr])
903 remove_group(block->nr-1, NOTSET, block, gn);
913 static gmx_bool check_have_atoms(t_atoms *atoms, char *string)
917 printf("Can not process '%s' without atom info, use option -f\n", string);
926 static gmx_bool parse_entry(char **string, int natoms, t_atoms *atoms,
927 t_blocka *block, char ***gn,
928 atom_id *nr, atom_id *index, char *gname)
930 static char **names, *ostring;
931 static gmx_bool bFirst = TRUE;
932 int j, n_names, sel_nr1;
933 atom_id i, nr1, *index1;
935 gmx_bool bRet, bCompl;
940 snew(names, MAXNAMES);
941 for (i = 0; i < MAXNAMES; i++)
943 snew(names[i], NAME_LEN+1);
950 while (*string[0] == ' ')
955 if ((*string)[0] == '!')
959 while (*string[0] == ' ')
971 if (parse_int(string, &sel_nr1) ||
972 parse_string(string, &sel_nr1, block->nr, *gn))
974 if ((sel_nr1 >= 0) && (sel_nr1 < block->nr))
976 copy_group(sel_nr1, block, nr, index);
977 strcpy(gname, (*gn)[sel_nr1]);
978 printf("Copied index group %d '%s'\n", sel_nr1, (*gn)[sel_nr1]);
983 printf("Group %d does not exist\n", sel_nr1);
986 else if ((*string)[0] == 'a')
989 if (check_have_atoms(atoms, ostring))
991 if (parse_int(string, &sel_nr1))
993 bRet = select_atomnumbers(string, atoms, sel_nr1, nr, index, gname);
995 else if (parse_names(string, &n_names, names))
997 bRet = select_atomnames(atoms, n_names, names, nr, index, FALSE);
998 make_gname(n_names, names, gname);
1002 else if ((*string)[0] == 't')
1005 if (check_have_atoms(atoms, ostring) &&
1006 parse_names(string, &n_names, names))
1008 if (atoms->atomtype == NULL)
1010 printf("Need a run input file to select atom types\n");
1014 bRet = select_atomnames(atoms, n_names, names, nr, index, TRUE);
1015 make_gname(n_names, names, gname);
1019 else if (strncmp(*string, "res", 3) == 0)
1022 if (check_have_atoms(atoms, ostring) &&
1023 parse_int(string, &sel_nr1) &&
1024 (sel_nr1 >= 0) && (sel_nr1 < block->nr) )
1026 bRet = atoms_from_residuenumbers(atoms,
1027 sel_nr1, block, nr, index, (*gn)[sel_nr1]);
1028 sprintf(gname, "atom_%s", (*gn)[sel_nr1]);
1031 else if (strncmp(*string, "ri", 2) == 0)
1034 if (check_have_atoms(atoms, ostring) &&
1035 parse_int_char(string, &sel_nr1, &c))
1037 bRet = select_residueindices(string, atoms, sel_nr1, c, nr, index, gname);
1040 else if ((*string)[0] == 'r')
1043 if (check_have_atoms(atoms, ostring))
1045 if (parse_int_char(string, &sel_nr1, &c))
1047 bRet = select_residuenumbers(string, atoms, sel_nr1, c, nr, index, gname);
1049 else if (parse_names(string, &n_names, names))
1051 bRet = select_residuenames(atoms, n_names, names, nr, index);
1052 make_gname(n_names, names, gname);
1056 else if (strncmp(*string, "chain", 5) == 0)
1059 if (check_have_atoms(atoms, ostring) &&
1060 parse_names(string, &n_names, names))
1062 bRet = select_chainnames(atoms, n_names, names, nr, index);
1063 sprintf(gname, "ch%s", names[0]);
1064 for (i = 1; i < n_names; i++)
1066 strcat(gname, names[i]);
1072 snew(index1, natoms-*nr);
1074 for (i = 0; i < natoms; i++)
1077 while ((j < *nr) && (index[j] != i))
1083 if (nr1 >= natoms-*nr)
1085 printf("There are double atoms in your index group\n");
1093 for (i = 0; i < nr1; i++)
1095 index[i] = index1[i];
1099 for (i = strlen(gname)+1; i > 0; i--)
1101 gname[i] = gname[i-1];
1104 printf("Complemented group: %d atoms\n", *nr);
1110 static void list_residues(t_atoms *atoms)
1112 int i, j, start, end, prev_resind, resind;
1115 /* Print all the residues, assuming continuous resnr count */
1116 start = atoms->atom[0].resind;
1117 prev_resind = start;
1118 for (i = 0; i < atoms->nr; i++)
1120 resind = atoms->atom[i].resind;
1121 if ((resind != prev_resind) || (i == atoms->nr-1))
1123 if ((bDiff = strcmp(*atoms->resinfo[resind].name,
1124 *atoms->resinfo[start].name)) ||
1137 for (j = start; j <= end; j++)
1140 j+1, *(atoms->resinfo[j].name));
1145 printf(" %4d - %4d %-5s ",
1146 start+1, end+1, *(atoms->resinfo[start].name));
1151 prev_resind = resind;
1156 static void edit_index(int natoms, t_atoms *atoms, rvec *x, t_blocka *block, char ***gn, gmx_bool bVerbose)
1158 static char **atnames, *ostring;
1159 static gmx_bool bFirst = TRUE;
1160 char inp_string[STRLEN], *string;
1161 char gname[STRLEN], gname1[STRLEN], gname2[STRLEN];
1162 int i, i0, i1, sel_nr, sel_nr2, newgroup;
1163 atom_id nr, nr1, nr2, *index, *index1, *index2;
1164 gmx_bool bAnd, bOr, bPrintOnce;
1169 snew(atnames, MAXNAMES);
1170 for (i = 0; i < MAXNAMES; i++)
1172 snew(atnames[i], NAME_LEN+1);
1178 snew(index, natoms);
1179 snew(index1, natoms);
1180 snew(index2, natoms);
1187 if (bVerbose || bPrintOnce || newgroup != NOTSET)
1190 if (bVerbose || bPrintOnce || newgroup == NOTSET)
1200 for (i = i0; i < i1; i++)
1202 printf("%3d %-20s: %5d atoms\n", i, (*gn)[i],
1203 block->index[i+1]-block->index[i]);
1207 if (bVerbose || bPrintOnce)
1210 printf(" nr : group ! 'name' nr name 'splitch' nr Enter: list groups\n");
1211 printf(" 'a': atom & 'del' nr 'splitres' nr 'l': list residues\n");
1212 printf(" 't': atom type | 'keep' nr 'splitat' nr 'h': help\n");
1213 printf(" 'r': residue 'res' nr 'chain' char\n");
1214 printf(" \"name\": group 'case': case %s 'q': save and quit\n",
1215 bCase ? "insensitive" : "sensitive ");
1216 printf(" 'ri': residue index\n");
1221 if (NULL == fgets(inp_string, STRLEN, stdin))
1223 gmx_fatal(FARGS, "Error reading user input");
1225 inp_string[strlen(inp_string)-1] = 0;
1227 string = inp_string;
1228 while (string[0] == ' ')
1235 if (string[0] == 'h')
1237 printf(" nr : selects an index group by number or quoted string.\n");
1238 printf(" The string is first matched against the whole group name,\n");
1239 printf(" then against the beginning and finally against an\n");
1240 printf(" arbitrary substring. A multiple match is an error.\n");
1242 printf(" 'a' nr1 [nr2 ...] : selects atoms, atom numbering starts at 1.\n");
1243 printf(" 'a' nr1 - nr2 : selects atoms in the range from nr1 to nr2.\n");
1244 printf(" 'a' name1[*] [name2[*] ...] : selects atoms by name(s), '?' matches any char,\n");
1245 printf(" wildcard '*' allowed at the end of a name.\n");
1246 printf(" 't' type1[*] [type2[*] ...] : as 'a', but for type, run input file required.\n");
1247 printf(" 'r' nr1[ic1] [nr2[ic2] ...] : selects residues by number and insertion code.\n");
1248 printf(" 'r' nr1 - nr2 : selects residues in the range from nr1 to nr2.\n");
1249 printf(" 'r' name1[*] [name2[*] ...] : as 'a', but for residue names.\n");
1250 printf(" 'ri' nr1 - nr2 : selects residue indices, 1-indexed, (as opposed to numbers) in the range from nr1 to nr2.\n");
1251 printf(" 'chain' ch1 [ch2 ...] : selects atoms by chain identifier(s),\n");
1252 printf(" not available with a .gro file as input.\n");
1253 printf(" ! : takes the complement of a group with respect to all\n");
1254 printf(" the atoms in the input file.\n");
1255 printf(" & | : AND and OR, can be placed between any of the options\n");
1256 printf(" above, the input is processed from left to right.\n");
1257 printf(" 'name' nr name : rename group nr to name.\n");
1258 printf(" 'del' nr1 [- nr2] : deletes one group or groups in the range from nr1 to nr2.\n");
1259 printf(" 'keep' nr : deletes all groups except nr.\n");
1260 printf(" 'case' : make all name compares case (in)sensitive.\n");
1261 printf(" 'splitch' nr : split group into chains using CA distances.\n");
1262 printf(" 'splitres' nr : split group into residues.\n");
1263 printf(" 'splitat' nr : split group into atoms.\n");
1264 printf(" 'res' nr : interpret numbers in group as residue numbers\n");
1265 printf(" Enter : list the currently defined groups and commands\n");
1266 printf(" 'l' : list the residues.\n");
1267 printf(" 'h' : show this help.\n");
1268 printf(" 'q' : save and quit.\n");
1270 printf(" Examples:\n");
1271 printf(" > 2 | 4 & r 3-5\n");
1272 printf(" selects all atoms from group 2 and 4 that have residue numbers 3, 4 or 5\n");
1273 printf(" > a C* & !a C CA\n");
1274 printf(" selects all atoms starting with 'C' but not the atoms 'C' and 'CA'\n");
1275 printf(" > \"protein\" & ! \"backb\"\n");
1276 printf(" selects all atoms that are in group 'protein' and not in group 'backbone'\n");
1279 printf("\npress Enter ");
1283 else if (strncmp(string, "del", 3) == 0)
1286 if (parse_int(&string, &sel_nr))
1288 while (string[0] == ' ')
1292 if (string[0] == '-')
1295 parse_int(&string, &sel_nr2);
1301 while (string[0] == ' ')
1305 if (string[0] == '\0')
1307 remove_group(sel_nr, sel_nr2, block, gn);
1311 printf("\nSyntax error: \"%s\"\n", string);
1315 else if (strncmp(string, "keep", 4) == 0)
1318 if (parse_int(&string, &sel_nr))
1320 remove_group(sel_nr+1, block->nr-1, block, gn);
1321 remove_group(0, sel_nr-1, block, gn);
1324 else if (strncmp(string, "name", 4) == 0)
1327 if (parse_int(&string, &sel_nr))
1329 if ((sel_nr >= 0) && (sel_nr < block->nr))
1331 sscanf(string, "%s", gname);
1332 sfree((*gn)[sel_nr]);
1333 (*gn)[sel_nr] = gmx_strdup(gname);
1337 else if (strncmp(string, "case", 4) == 0)
1340 printf("Switched to case %s\n", bCase ? "sensitive" : "insensitive");
1342 else if (string[0] == 'v')
1344 bVerbose = !bVerbose;
1345 printf("Turned verbose %s\n", bVerbose ? "on" : "off");
1347 else if (string[0] == 'l')
1349 if (check_have_atoms(atoms, ostring) )
1351 list_residues(atoms);
1354 else if (strncmp(string, "splitch", 7) == 0)
1357 if (check_have_atoms(atoms, ostring) &&
1358 parse_int(&string, &sel_nr) &&
1359 (sel_nr >= 0) && (sel_nr < block->nr))
1361 split_chain(atoms, x, sel_nr, block, gn);
1364 else if (strncmp(string, "splitres", 8) == 0)
1367 if (check_have_atoms(atoms, ostring) &&
1368 parse_int(&string, &sel_nr) &&
1369 (sel_nr >= 0) && (sel_nr < block->nr))
1371 split_group(atoms, sel_nr, block, gn, FALSE);
1374 else if (strncmp(string, "splitat", 7) == 0)
1377 if (check_have_atoms(atoms, ostring) &&
1378 parse_int(&string, &sel_nr) &&
1379 (sel_nr >= 0) && (sel_nr < block->nr))
1381 split_group(atoms, sel_nr, block, gn, TRUE);
1384 else if (string[0] == '\0')
1388 else if (string[0] != 'q')
1392 if (parse_entry(&string, natoms, atoms, block, gn, &nr, index, gname))
1396 while (string[0] == ' ')
1403 if (string[0] == '&')
1407 else if (string[0] == '|')
1416 for (i = 0; i < nr; i++)
1418 index1[i] = index[i];
1420 strcpy(gname1, gname);
1421 if (parse_entry(&string, natoms, atoms, block, gn, &nr2, index2, gname2))
1425 or_groups(nr1, index1, nr2, index2, &nr, index);
1426 sprintf(gname, "%s_%s", gname1, gname2);
1430 and_groups(nr1, index1, nr2, index2, &nr, index);
1431 sprintf(gname, "%s_&_%s", gname1, gname2);
1436 while (bAnd || bOr);
1438 while (string[0] == ' ')
1444 printf("\nSyntax error: \"%s\"\n", string);
1448 copy2block(nr, index, block);
1449 srenew(*gn, block->nr);
1450 newgroup = block->nr-1;
1451 (*gn)[newgroup] = gmx_strdup(gname);
1455 printf("Group is empty\n");
1459 while (string[0] != 'q');
1466 static int block2natoms(t_blocka *block)
1471 for (i = 0; i < block->nra; i++)
1473 natoms = max(natoms, block->a[i]);
1480 void merge_blocks(t_blocka *dest, t_blocka *source)
1484 /* count groups, srenew and fill */
1487 dest->nr += source->nr;
1488 srenew(dest->index, dest->nr+1);
1489 for (i = 0; i < source->nr; i++)
1491 dest->index[i0+i] = nra0 + source->index[i];
1493 /* count atoms, srenew and fill */
1494 dest->nra += source->nra;
1495 srenew(dest->a, dest->nra);
1496 for (i = 0; i < source->nra; i++)
1498 dest->a[nra0+i] = source->a[i];
1501 /* terminate list */
1502 dest->index[dest->nr] = dest->nra;
1506 int gmx_make_ndx(int argc, char *argv[])
1508 const char *desc[] = {
1509 "Index groups are necessary for almost every GROMACS program.",
1510 "All these programs can generate default index groups. You ONLY",
1511 "have to use [THISMODULE] when you need SPECIAL index groups.",
1512 "There is a default index group for the whole system, 9 default",
1513 "index groups for proteins, and a default index group",
1514 "is generated for every other residue name.[PAR]",
1515 "When no index file is supplied, also [THISMODULE] will generate the",
1517 "With the index editor you can select on atom, residue and chain names",
1519 "When a run input file is supplied you can also select on atom type.",
1520 "You can use NOT, AND and OR, you can split groups",
1521 "into chains, residues or atoms. You can delete and rename groups.[PAR]",
1522 "The atom numbering in the editor and the index file starts at 1.[PAR]",
1523 "The [TT]-twin[tt] switch duplicates all index groups with an offset of",
1524 "[TT]-natoms[tt], which is useful for Computational Electrophysiology",
1525 "double-layer membrane setups."
1528 static int natoms = 0;
1529 static gmx_bool bVerbose = FALSE;
1530 static gmx_bool bDuplicate = FALSE;
1532 { "-natoms", FALSE, etINT, {&natoms},
1533 "set number of atoms (default: read from coordinate or index file)" },
1534 { "-twin", FALSE, etBOOL, {&bDuplicate},
1535 "Duplicate all index groups with an offset of -natoms" },
1536 { "-verbose", FALSE, etBOOL, {&bVerbose},
1537 "HIDDENVerbose output" }
1539 #define NPA asize(pa)
1544 const char *stxfile;
1546 const char *ndxoutfile;
1553 t_blocka *block, *block2;
1554 char **gnames, **gnames2;
1556 { efSTX, "-f", NULL, ffOPTRD },
1557 { efNDX, "-n", NULL, ffOPTRDMULT },
1558 { efNDX, "-o", NULL, ffWRITE }
1560 #define NFILE asize(fnm)
1562 if (!parse_common_args(&argc, argv, 0, NFILE, fnm, NPA, pa, asize(desc), desc,
1568 stxfile = ftp2fn_null(efSTX, NFILE, fnm);
1569 if (opt2bSet("-n", NFILE, fnm))
1571 nndxin = opt2fns(&ndxinfiles, "-n", NFILE, fnm);
1577 ndxoutfile = opt2fn("-o", NFILE, fnm);
1578 bNatoms = opt2parg_bSet("-natoms", NPA, pa);
1580 if (!stxfile && !nndxin)
1582 gmx_fatal(FARGS, "No input files (structure or index)");
1588 get_stx_coordnum(stxfile, &(atoms->nr));
1589 init_t_atoms(atoms, atoms->nr, TRUE);
1592 fprintf(stderr, "\nReading structure file\n");
1593 read_stx_conf(stxfile, title, atoms, x, v, &ePBC, box);
1603 /* read input file(s) */
1604 block = new_blocka();
1606 printf("Going to read %d old index file(s)\n", nndxin);
1609 for (i = 0; i < nndxin; i++)
1611 block2 = init_index(ndxinfiles[i], &gnames2);
1612 srenew(gnames, block->nr+block2->nr);
1613 for (j = 0; j < block2->nr; j++)
1615 gnames[block->nr+j] = gnames2[j];
1618 merge_blocks(block, block2);
1620 sfree(block2->index);
1621 /* done_block(block2); */
1628 analyse(atoms, block, &gnames, FALSE, TRUE);
1633 natoms = block2natoms(block);
1634 printf("Counted atom numbers up to %d in index file\n", natoms);
1637 edit_index(natoms, atoms, x, block, &gnames, bVerbose);
1639 write_index(ndxoutfile, block, gnames, bDuplicate, natoms);