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,2015,2016,2017,2018,2019, 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/commandline/pargs.h"
45 #include "gromacs/fileio/confio.h"
46 #include "gromacs/gmxana/gmx_ana.h"
47 #include "gromacs/math/vec.h"
48 #include "gromacs/topology/block.h"
49 #include "gromacs/topology/index.h"
50 #include "gromacs/topology/topology.h"
51 #include "gromacs/utility/arraysize.h"
52 #include "gromacs/utility/cstringutil.h"
53 #include "gromacs/utility/fatalerror.h"
54 #include "gromacs/utility/futil.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.
62 static const int NOTSET = -92637;
64 static gmx_bool bCase = FALSE;
66 static int or_groups(int nr1, const int *at1, int nr2, const int *at2,
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(int nr1, const int *at1, int nr2, const int *at2,
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' && std::strchr(spec, c) == nullptr);
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 (std::isdigit((*string)[0]))
218 *nr = (*string)[0]-'0';
220 while (std::isdigit((*string)[0]))
222 *nr = (*nr)*10+(*string)[0]-'0';
225 if (std::isalpha((*string)[0]))
230 /* Check if there is at most one non-digit character */
231 if (!std::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 = gmx_strdup((*string));
285 sp = std::strchr(s, c);
288 (*string) += sp-s + 1;
290 (*nr) = find_group(s, ngrps, grpname);
294 return (*nr) != NOTSET;
297 static int select_atomnumbers(char **string, const t_atoms *atoms, int n1,
298 int *nr, int *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);
332 std::strcpy(gname, buf);
341 if ((i-1 >= 0) && (i-1 < atoms->nr))
345 sprintf(buf, "_%d", i);
346 std::strcat(gname, buf);
350 printf("Invalid atom number %d\n", i);
354 while ((*nr != 0) && (parse_int(string, &i)));
360 static int select_residuenumbers(char **string, const t_atoms *atoms,
362 int *nr, int *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);
406 std::strcpy(gname, buf);
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);
425 std::strcat(gname, buf);
427 while (parse_int_char(string, &j, &c));
433 static int select_residueindices(char **string, const t_atoms *atoms,
435 int *nr, int *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);
481 std::strcpy(gname, buf);
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);
500 std::strcat(gname, buf);
502 while (parse_int_char(string, &j, &c));
509 static gmx_bool atoms_from_residuenumbers(const t_atoms *atoms, int group, t_blocka *block,
510 int *nr, int *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(const char *name, const char *search)
546 gmx_bool matches = TRUE;
548 // Loop while name and search are not end-of-string and matches is true
549 for (; *name && *search && matches; name++, search++)
553 // still matching, continue to next character
556 else if (*search == '*')
560 printf("WARNING: Currently '*' is only supported at the end of an expression\n");
562 // if * is the last char in search string, we have a match,
563 // otherwise we just failed. Return in either case, we're done.
564 return (*(search+1) == '\0');
567 // Compare one character
570 matches = (*name == *search);
574 matches = (std::toupper(*name) == std::toupper(*search));
578 matches = matches && (*name == '\0' && (*search == '\0' || *search == '*'));
583 static int select_chainnames(const t_atoms *atoms, int n_names, char **names,
592 for (i = 0; i < atoms->nr; i++)
594 name[0] = atoms->resinfo[atoms->atom[i].resind].chainid;
596 while (j < n_names && !comp_name(name, names[j]))
606 printf("Found %d atom%s with chain identifier%s",
607 *nr, (*nr == 1) ? "" : "s", (n_names == 1) ? "" : "s");
608 for (j = 0; (j < n_names); j++)
610 printf(" %s", names[j]);
617 static int select_atomnames(const t_atoms *atoms, int n_names, char **names,
618 int *nr, int *index, gmx_bool bType)
625 for (i = 0; i < atoms->nr; i++)
629 name = *(atoms->atomtype[i]);
633 name = *(atoms->atomname[i]);
636 while (j < n_names && !comp_name(name, names[j]))
646 printf("Found %d atoms with %s%s",
647 *nr, bType ? "type" : "name", (n_names == 1) ? "" : "s");
648 for (j = 0; (j < n_names); j++)
650 printf(" %s", names[j]);
657 static int select_residuenames(const t_atoms *atoms, int n_names, char **names,
665 for (i = 0; i < atoms->nr; i++)
667 name = *(atoms->resinfo[atoms->atom[i].resind].name);
669 while (j < n_names && !comp_name(name, names[j]))
679 printf("Found %d atoms with residue name%s", *nr, (n_names == 1) ? "" : "s");
680 for (j = 0; (j < n_names); j++)
682 printf(" %s", names[j]);
689 static void copy2block(int n, const int *index, t_blocka *block)
696 srenew(block->index, block->nr+1);
697 block->index[block->nr] = n0+n;
698 srenew(block->a, n0+n);
699 for (i = 0; (i < n); i++)
701 block->a[n0+i] = index[i];
705 static void make_gname(int n, char **names, char *gname)
709 std::strcpy(gname, names[0]);
710 for (i = 1; i < n; i++)
712 std::strcat(gname, "_");
713 std::strcat(gname, names[i]);
717 static void copy_group(int g, t_blocka *block, int *nr, int *index)
721 i0 = block->index[g];
722 *nr = block->index[g+1]-i0;
723 for (i = 0; i < *nr; i++)
725 index[i] = block->a[i0+i];
729 static void remove_group(int nr, int nr2, t_blocka *block, char ***gn)
739 for (j = 0; j <= nr2-nr; j++)
741 if ((nr < 0) || (nr >= block->nr))
743 printf("Group %d does not exist\n", nr+j);
747 shift = block->index[nr+1]-block->index[nr];
748 for (i = block->index[nr+1]; i < block->nra; i++)
750 block->a[i-shift] = block->a[i];
753 for (i = nr; i < block->nr; i++)
755 block->index[i] = block->index[i+1]-shift;
757 name = gmx_strdup((*gn)[nr]);
759 for (i = nr; i < block->nr-1; i++)
761 (*gn)[i] = (*gn)[i+1];
764 block->nra = block->index[block->nr];
765 printf("Removed group %d '%s'\n", nr+j, name);
771 static void split_group(const t_atoms *atoms, int sel_nr, t_blocka *block, char ***gn,
774 char buf[STRLEN], *name;
778 printf("Splitting group %d '%s' into %s\n", sel_nr, (*gn)[sel_nr],
779 bAtom ? "atoms" : "residues");
781 n0 = block->index[sel_nr];
782 n1 = block->index[sel_nr+1];
783 srenew(block->a, block->nra+n1-n0);
784 for (i = n0; i < n1; i++)
787 resind = atoms->atom[a].resind;
788 name = *(atoms->resinfo[resind].name);
789 if (bAtom || (i == n0) || (atoms->atom[block->a[i-1]].resind != resind))
793 block->index[block->nr] = block->nra;
796 srenew(block->index, block->nr+1);
797 srenew(*gn, block->nr);
800 sprintf(buf, "%s_%s_%d", (*gn)[sel_nr], *atoms->atomname[a], a+1);
804 sprintf(buf, "%s_%s_%d", (*gn)[sel_nr], name, atoms->resinfo[resind].nr);
806 (*gn)[block->nr-1] = gmx_strdup(buf);
808 block->a[block->nra] = a;
811 block->index[block->nr] = block->nra;
814 static int split_chain(const t_atoms *atoms, const rvec *x,
815 int sel_nr, t_blocka *block, char ***gn)
819 int i, a, natoms, *start = nullptr, *end = nullptr, ca_start, ca_end;
826 while (ca_start < natoms)
828 while ((ca_start < natoms) && std::strcmp(*atoms->atomname[ca_start], "CA") != 0)
832 if (ca_start < natoms)
834 srenew(start, nchain+1);
835 srenew(end, nchain+1);
836 start[nchain] = ca_start;
837 while ((start[nchain] > 0) &&
838 (atoms->atom[start[nchain]-1].resind ==
839 atoms->atom[ca_start].resind))
852 while ((i < natoms) && std::strcmp(*atoms->atomname[i], "CA") != 0);
855 rvec_sub(x[ca_end], x[i], vec);
862 while (norm(vec) < 0.45);
864 end[nchain] = ca_end;
865 while ((end[nchain]+1 < natoms) &&
866 (atoms->atom[end[nchain]+1].resind == atoms->atom[ca_end].resind))
870 ca_start = end[nchain]+1;
876 printf("Found 1 chain, will not split\n");
880 printf("Found %d chains\n", nchain);
882 for (j = 0; j < nchain; j++)
884 printf("%d:%6d atoms (%d to %d)\n",
885 j+1, end[j]-start[j]+1, start[j]+1, end[j]+1);
890 srenew(block->a, block->nra+block->index[sel_nr+1]-block->index[sel_nr]);
891 for (j = 0; j < nchain; j++)
894 srenew(block->index, block->nr+1);
895 srenew(*gn, block->nr);
896 sprintf(buf, "%s_chain%d", (*gn)[sel_nr], j+1);
897 (*gn)[block->nr-1] = gmx_strdup(buf);
898 for (i = block->index[sel_nr]; i < block->index[sel_nr+1]; i++)
901 if ((a >= start[j]) && (a <= end[j]))
903 block->a[block->nra] = a;
907 block->index[block->nr] = block->nra;
908 if (block->index[block->nr-1] == block->index[block->nr])
910 remove_group(block->nr-1, NOTSET, block, gn);
920 static gmx_bool check_have_atoms(const t_atoms *atoms, char *string)
922 if (atoms == nullptr)
924 printf("Can not process '%s' without atom info, use option -f\n", string);
933 static gmx_bool parse_entry(char **string, int natoms, const t_atoms *atoms,
934 t_blocka *block, char ***gn,
935 int *nr, int *index, char *gname)
937 static char **names, *ostring;
938 static gmx_bool bFirst = TRUE;
939 int j, n_names, sel_nr1;
942 gmx_bool bRet, bCompl;
947 snew(names, MAXNAMES);
948 for (i = 0; i < MAXNAMES; i++)
950 snew(names[i], NAME_LEN+1);
957 while (*string[0] == ' ')
962 if ((*string)[0] == '!')
966 while (*string[0] == ' ')
978 if (parse_int(string, &sel_nr1) ||
979 parse_string(string, &sel_nr1, block->nr, *gn))
981 if ((sel_nr1 >= 0) && (sel_nr1 < block->nr))
983 copy_group(sel_nr1, block, nr, index);
984 std::strcpy(gname, (*gn)[sel_nr1]);
985 printf("Copied index group %d '%s'\n", sel_nr1, (*gn)[sel_nr1]);
990 printf("Group %d does not exist\n", sel_nr1);
993 else if ((*string)[0] == 'a')
996 if (check_have_atoms(atoms, ostring))
998 if (parse_int(string, &sel_nr1))
1000 bRet = (select_atomnumbers(string, atoms, sel_nr1, nr, index, gname) != 0);
1002 else if (parse_names(string, &n_names, names))
1004 bRet = (select_atomnames(atoms, n_names, names, nr, index, FALSE) != 0);
1005 make_gname(n_names, names, gname);
1009 else if ((*string)[0] == 't')
1012 if (check_have_atoms(atoms, ostring) &&
1013 parse_names(string, &n_names, names))
1015 if (!(atoms->haveType))
1017 printf("Need a run input file to select atom types\n");
1021 bRet = (select_atomnames(atoms, n_names, names, nr, index, TRUE) != 0);
1022 make_gname(n_names, names, gname);
1026 else if (std::strncmp(*string, "res", 3) == 0)
1029 if (check_have_atoms(atoms, ostring) &&
1030 parse_int(string, &sel_nr1) &&
1031 (sel_nr1 >= 0) && (sel_nr1 < block->nr) )
1033 bRet = atoms_from_residuenumbers(atoms,
1034 sel_nr1, block, nr, index, (*gn)[sel_nr1]);
1035 sprintf(gname, "atom_%s", (*gn)[sel_nr1]);
1038 else if (std::strncmp(*string, "ri", 2) == 0)
1041 if (check_have_atoms(atoms, ostring) &&
1042 parse_int_char(string, &sel_nr1, &c))
1044 bRet = (select_residueindices(string, atoms, sel_nr1, c, nr, index, gname) != 0);
1047 else if ((*string)[0] == 'r')
1050 if (check_have_atoms(atoms, ostring))
1052 if (parse_int_char(string, &sel_nr1, &c))
1054 bRet = (select_residuenumbers(string, atoms, sel_nr1, c, nr, index, gname) != 0);
1056 else if (parse_names(string, &n_names, names))
1058 bRet = (select_residuenames(atoms, n_names, names, nr, index) != 0);
1059 make_gname(n_names, names, gname);
1063 else if (std::strncmp(*string, "chain", 5) == 0)
1066 if (check_have_atoms(atoms, ostring) &&
1067 parse_names(string, &n_names, names))
1069 bRet = (select_chainnames(atoms, n_names, names, nr, index) != 0);
1070 sprintf(gname, "ch%s", names[0]);
1071 for (i = 1; i < n_names; i++)
1073 std::strcat(gname, names[i]);
1079 snew(index1, natoms-*nr);
1081 for (i = 0; i < natoms; i++)
1084 while ((j < *nr) && (index[j] != i))
1090 if (nr1 >= natoms-*nr)
1092 printf("There are double atoms in your index group\n");
1100 for (i = 0; i < nr1; i++)
1102 index[i] = index1[i];
1106 for (i = std::strlen(gname)+1; i > 0; i--)
1108 gname[i] = gname[i-1];
1111 printf("Complemented group: %d atoms\n", *nr);
1117 static void list_residues(const t_atoms *atoms)
1119 int i, j, start, end, prev_resind, resind;
1122 /* Print all the residues, assuming continuous resnr count */
1123 start = atoms->atom[0].resind;
1124 prev_resind = start;
1125 for (i = 0; i < atoms->nr; i++)
1127 resind = atoms->atom[i].resind;
1128 if ((resind != prev_resind) || (i == atoms->nr-1))
1130 if ((bDiff = (std::strcmp(*atoms->resinfo[resind].name,
1131 *atoms->resinfo[start].name) != 0)) ||
1144 for (j = start; j <= end; j++)
1147 j+1, *(atoms->resinfo[j].name));
1152 printf(" %4d - %4d %-5s ",
1153 start+1, end+1, *(atoms->resinfo[start].name));
1158 prev_resind = resind;
1163 static void edit_index(int natoms, const t_atoms *atoms, const rvec *x, t_blocka *block, char ***gn, gmx_bool bVerbose)
1165 static char **atnames, *ostring;
1166 static gmx_bool bFirst = TRUE;
1167 char inp_string[STRLEN], *string;
1168 char gname[STRLEN*3], gname1[STRLEN], gname2[STRLEN];
1169 int i, i0, i1, sel_nr, sel_nr2, newgroup;
1170 int nr, nr1, nr2, *index, *index1, *index2;
1171 gmx_bool bAnd, bOr, bPrintOnce;
1176 snew(atnames, MAXNAMES);
1177 for (i = 0; i < MAXNAMES; i++)
1179 snew(atnames[i], NAME_LEN+1);
1185 snew(index, natoms);
1186 snew(index1, natoms);
1187 snew(index2, natoms);
1194 if (bVerbose || bPrintOnce || newgroup != NOTSET)
1197 if (bVerbose || bPrintOnce || newgroup == NOTSET)
1207 for (i = i0; i < i1; i++)
1209 printf("%3d %-20s: %5d atoms\n", i, (*gn)[i],
1210 block->index[i+1]-block->index[i]);
1214 if (bVerbose || bPrintOnce)
1217 printf(" nr : group '!': not 'name' nr name 'splitch' nr Enter: list groups\n");
1218 printf(" 'a': atom '&': and 'del' nr 'splitres' nr 'l': list residues\n");
1219 printf(" 't': atom type '|': or 'keep' nr 'splitat' nr 'h': help\n");
1220 printf(" 'r': residue 'res' nr 'chain' char\n");
1221 printf(" \"name\": group 'case': case %s 'q': save and quit\n",
1222 bCase ? "insensitive" : "sensitive ");
1223 printf(" 'ri': residue index\n");
1228 if (nullptr == fgets(inp_string, STRLEN, stdin))
1230 gmx_fatal(FARGS, "Error reading user input");
1232 inp_string[std::strlen(inp_string)-1] = 0;
1234 string = inp_string;
1235 while (string[0] == ' ')
1242 if (string[0] == 'h')
1244 printf(" nr : selects an index group by number or quoted string.\n");
1245 printf(" The string is first matched against the whole group name,\n");
1246 printf(" then against the beginning and finally against an\n");
1247 printf(" arbitrary substring. A multiple match is an error.\n");
1249 printf(" 'a' nr1 [nr2 ...] : selects atoms, atom numbering starts at 1.\n");
1250 printf(" 'a' nr1 - nr2 : selects atoms in the range from nr1 to nr2.\n");
1251 printf(" 'a' name1[*] [name2[*] ...] : selects atoms by name(s), '?' matches any char,\n");
1252 printf(" wildcard '*' allowed at the end of a name.\n");
1253 printf(" 't' type1[*] [type2[*] ...] : as 'a', but for type, run input file required.\n");
1254 printf(" 'r' nr1[ic1] [nr2[ic2] ...] : selects residues by number and insertion code.\n");
1255 printf(" 'r' nr1 - nr2 : selects residues in the range from nr1 to nr2.\n");
1256 printf(" 'r' name1[*] [name2[*] ...] : as 'a', but for residue names.\n");
1257 printf(" 'ri' nr1 - nr2 : selects residue indices, 1-indexed, (as opposed to numbers) in the range from nr1 to nr2.\n");
1258 printf(" 'chain' ch1 [ch2 ...] : selects atoms by chain identifier(s),\n");
1259 printf(" not available with a .gro file as input.\n");
1260 printf(" ! : takes the complement of a group with respect to all\n");
1261 printf(" the atoms in the input file.\n");
1262 printf(" & | : AND and OR, can be placed between any of the options\n");
1263 printf(" above, the input is processed from left to right.\n");
1264 printf(" 'name' nr name : rename group nr to name.\n");
1265 printf(" 'del' nr1 [- nr2] : deletes one group or groups in the range from nr1 to nr2.\n");
1266 printf(" 'keep' nr : deletes all groups except nr.\n");
1267 printf(" 'case' : make all name compares case (in)sensitive.\n");
1268 printf(" 'splitch' nr : split group into chains using CA distances.\n");
1269 printf(" 'splitres' nr : split group into residues.\n");
1270 printf(" 'splitat' nr : split group into atoms.\n");
1271 printf(" 'res' nr : interpret numbers in group as residue numbers\n");
1272 printf(" Enter : list the currently defined groups and commands\n");
1273 printf(" 'l' : list the residues.\n");
1274 printf(" 'h' : show this help.\n");
1275 printf(" 'q' : save and quit.\n");
1277 printf(" Examples:\n");
1278 printf(" > 2 | 4 & r 3-5\n");
1279 printf(" selects all atoms from group 2 and 4 that have residue numbers 3, 4 or 5\n");
1280 printf(" > a C* & !a C CA\n");
1281 printf(" selects all atoms starting with 'C' but not the atoms 'C' and 'CA'\n");
1282 printf(" > \"protein\" & ! \"backb\"\n");
1283 printf(" selects all atoms that are in group 'protein' and not in group 'backbone'\n");
1286 printf("\npress Enter ");
1290 else if (std::strncmp(string, "del", 3) == 0)
1293 if (parse_int(&string, &sel_nr))
1295 while (string[0] == ' ')
1299 if (string[0] == '-')
1302 parse_int(&string, &sel_nr2);
1308 while (string[0] == ' ')
1312 if (string[0] == '\0')
1314 remove_group(sel_nr, sel_nr2, block, gn);
1318 printf("\nSyntax error: \"%s\"\n", string);
1322 else if (std::strncmp(string, "keep", 4) == 0)
1325 if (parse_int(&string, &sel_nr))
1327 remove_group(sel_nr+1, block->nr-1, block, gn);
1328 remove_group(0, sel_nr-1, block, gn);
1331 else if (std::strncmp(string, "name", 4) == 0)
1334 if (parse_int(&string, &sel_nr))
1336 if ((sel_nr >= 0) && (sel_nr < block->nr))
1338 sscanf(string, "%s", gname);
1339 sfree((*gn)[sel_nr]);
1340 (*gn)[sel_nr] = gmx_strdup(gname);
1344 else if (std::strncmp(string, "case", 4) == 0)
1347 printf("Switched to case %s\n", bCase ? "sensitive" : "insensitive");
1349 else if (string[0] == 'v')
1351 bVerbose = !bVerbose;
1352 printf("Turned verbose %s\n", bVerbose ? "on" : "off");
1354 else if (string[0] == 'l')
1356 if (check_have_atoms(atoms, ostring) )
1358 list_residues(atoms);
1361 else if (std::strncmp(string, "splitch", 7) == 0)
1364 if (check_have_atoms(atoms, ostring) &&
1365 parse_int(&string, &sel_nr) &&
1366 (sel_nr >= 0) && (sel_nr < block->nr))
1368 split_chain(atoms, x, sel_nr, block, gn);
1371 else if (std::strncmp(string, "splitres", 8) == 0)
1374 if (check_have_atoms(atoms, ostring) &&
1375 parse_int(&string, &sel_nr) &&
1376 (sel_nr >= 0) && (sel_nr < block->nr))
1378 split_group(atoms, sel_nr, block, gn, FALSE);
1381 else if (std::strncmp(string, "splitat", 7) == 0)
1384 if (check_have_atoms(atoms, ostring) &&
1385 parse_int(&string, &sel_nr) &&
1386 (sel_nr >= 0) && (sel_nr < block->nr))
1388 split_group(atoms, sel_nr, block, gn, TRUE);
1391 else if (string[0] == '\0')
1395 else if (string[0] != 'q')
1398 if (parse_entry(&string, natoms, atoms, block, gn, &nr, index, gname))
1402 while (string[0] == ' ')
1409 if (string[0] == '&')
1413 else if (string[0] == '|')
1422 for (i = 0; i < nr; i++)
1424 index1[i] = index[i];
1426 std::strcpy(gname1, gname);
1427 if (parse_entry(&string, natoms, atoms, block, gn, &nr2, index2, gname2))
1431 or_groups(nr1, index1, nr2, index2, &nr, index);
1432 sprintf(gname, "%s_%s", gname1, gname2);
1436 and_groups(nr1, index1, nr2, index2, &nr, index);
1437 sprintf(gname, "%s_&_%s", gname1, gname2);
1442 while (bAnd || bOr);
1444 while (string[0] == ' ')
1450 printf("\nSyntax error: \"%s\"\n", string);
1454 copy2block(nr, index, block);
1455 srenew(*gn, block->nr);
1456 newgroup = block->nr-1;
1457 (*gn)[newgroup] = gmx_strdup(gname);
1461 printf("Group is empty\n");
1465 while (string[0] != 'q');
1472 static int block2natoms(t_blocka *block)
1477 for (i = 0; i < block->nra; i++)
1479 natoms = std::max(natoms, block->a[i]);
1486 static void merge_blocks(t_blocka *dest, t_blocka *source)
1490 /* count groups, srenew and fill */
1493 dest->nr += source->nr;
1494 srenew(dest->index, dest->nr+1);
1495 for (i = 0; i < source->nr; i++)
1497 dest->index[i0+i] = nra0 + source->index[i];
1499 /* count atoms, srenew and fill */
1500 dest->nra += source->nra;
1501 srenew(dest->a, dest->nra);
1502 for (i = 0; i < source->nra; i++)
1504 dest->a[nra0+i] = source->a[i];
1507 /* terminate list */
1508 dest->index[dest->nr] = dest->nra;
1512 int gmx_make_ndx(int argc, char *argv[])
1514 const char *desc[] = {
1515 "Index groups are necessary for almost every GROMACS program.",
1516 "All these programs can generate default index groups. You ONLY",
1517 "have to use [THISMODULE] when you need SPECIAL index groups.",
1518 "There is a default index group for the whole system, 9 default",
1519 "index groups for proteins, and a default index group",
1520 "is generated for every other residue name.",
1522 "When no index file is supplied, also [THISMODULE] will generate the",
1524 "With the index editor you can select on atom, residue and chain names",
1526 "When a run input file is supplied you can also select on atom type.",
1527 "You can use boolean operations, you can split groups",
1528 "into chains, residues or atoms. You can delete and rename groups.",
1529 "Type 'h' in the editor for more details.",
1531 "The atom numbering in the editor and the index file starts at 1.",
1533 "The [TT]-twin[tt] switch duplicates all index groups with an offset of",
1534 "[TT]-natoms[tt], which is useful for Computational Electrophysiology",
1535 "double-layer membrane setups.",
1537 "See also [gmx-select] [TT]-on[tt], which provides an alternative way",
1538 "for constructing index groups. It covers nearly all of [THISMODULE]",
1539 "functionality, and in many cases much more."
1542 static int natoms = 0;
1543 static gmx_bool bVerbose = FALSE;
1544 static gmx_bool bDuplicate = FALSE;
1546 { "-natoms", FALSE, etINT, {&natoms},
1547 "set number of atoms (default: read from coordinate or index file)" },
1548 { "-twin", FALSE, etBOOL, {&bDuplicate},
1549 "Duplicate all index groups with an offset of -natoms" },
1550 { "-verbose", FALSE, etBOOL, {&bVerbose},
1551 "HIDDENVerbose output" }
1553 #define NPA asize(pa)
1555 gmx_output_env_t *oenv;
1556 const char *stxfile;
1557 const char *ndxoutfile;
1564 t_blocka *block, *block2;
1565 char **gnames, **gnames2;
1567 { efSTX, "-f", nullptr, ffOPTRD },
1568 { efNDX, "-n", nullptr, ffOPTRDMULT },
1569 { efNDX, "-o", nullptr, ffWRITE }
1571 #define NFILE asize(fnm)
1573 if (!parse_common_args(&argc, argv, 0, NFILE, fnm, NPA, pa, asize(desc), desc,
1579 stxfile = ftp2fn_null(efSTX, NFILE, fnm);
1580 gmx::ArrayRef<const std::string> ndxInFiles = opt2fnsIfOptionSet("-n", NFILE, fnm);
1581 ndxoutfile = opt2fn("-o", NFILE, fnm);
1582 bNatoms = opt2parg_bSet("-natoms", NPA, pa);
1584 if (!stxfile && ndxInFiles.empty())
1586 gmx_fatal(FARGS, "No input files (structure or index)");
1593 fprintf(stderr, "\nReading structure file\n");
1594 read_tps_conf(stxfile, top, &ePBC, &x, &v, box, FALSE);
1595 atoms = &top->atoms;
1596 if (atoms->pdbinfo == nullptr)
1598 snew(atoms->pdbinfo, atoms->nr);
1609 /* read input file(s) */
1610 block = new_blocka();
1612 printf("Going to read %td old index file(s)\n", ndxInFiles.ssize());
1613 if (!ndxInFiles.empty())
1615 for (const std::string &ndxInFile : ndxInFiles)
1617 block2 = init_index(ndxInFile.c_str(), &gnames2);
1618 srenew(gnames, block->nr+block2->nr);
1619 for (j = 0; j < block2->nr; j++)
1621 gnames[block->nr+j] = gnames2[j];
1624 merge_blocks(block, block2);
1626 sfree(block2->index);
1627 /* done_block(block2); */
1634 analyse(atoms, block, &gnames, FALSE, TRUE);
1639 natoms = block2natoms(block);
1640 printf("Counted atom numbers up to %d in index file\n", natoms);
1643 edit_index(natoms, atoms, x, block, &gnames, bVerbose);
1645 write_index(ndxoutfile, block, gnames, bDuplicate, natoms);