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 by the GROMACS development team.
7 * Copyright (c) 2017,2018,2019,2020, by the GROMACS development team, led by
8 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
9 * and including many others, as listed in the AUTHORS file in the
10 * top-level source directory and at http://www.gromacs.org.
12 * GROMACS is free software; you can redistribute it and/or
13 * modify it under the terms of the GNU Lesser General Public License
14 * as published by the Free Software Foundation; either version 2.1
15 * of the License, or (at your option) any later version.
17 * GROMACS is distributed in the hope that it will be useful,
18 * but WITHOUT ANY WARRANTY; without even the implied warranty of
19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 * Lesser General Public License for more details.
22 * You should have received a copy of the GNU Lesser General Public
23 * License along with GROMACS; if not, see
24 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
25 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
27 * If you want to redistribute modifications to GROMACS, please
28 * consider that scientific software is very special. Version
29 * control is crucial - bugs must be traceable. We will be happy to
30 * consider code for inclusion in the official distribution, but
31 * derived work must not be called official GROMACS. Details are found
32 * in the README & COPYING files - if they are missing, get the
33 * official version at http://www.gromacs.org.
35 * To help us fund GROMACS development, we humbly ask that you cite
36 * the research papers on the package. Check out http://www.gromacs.org.
47 #include "gromacs/commandline/pargs.h"
48 #include "gromacs/fileio/confio.h"
49 #include "gromacs/math/vec.h"
50 #include "gromacs/topology/block.h"
51 #include "gromacs/topology/index.h"
52 #include "gromacs/topology/mtop_util.h"
53 #include "gromacs/topology/topology.h"
54 #include "gromacs/utility/arraysize.h"
55 #include "gromacs/utility/cstringutil.h"
56 #include "gromacs/utility/fatalerror.h"
57 #include "gromacs/utility/futil.h"
58 #include "gromacs/utility/smalloc.h"
60 /* It's not nice to have size limits, but we should not spend more time
61 * on this ancient tool, but instead use the new selection library.
65 static const int NOTSET = -92637;
67 static gmx_bool bCase = FALSE;
69 static int or_groups(int nr1, const int* at1, int nr2, const int* at2, int* nr, int* at)
77 for (i1 = 0; i1 < nr1; i1++)
79 if ((i1 > 0) && (at1[i1] <= max))
85 for (i1 = 0; i1 < nr2; i1++)
87 if ((i1 > 0) && (at2[i1] <= max))
96 printf("One of your groups is not ascending\n");
103 while ((i1 < nr1) || (i2 < nr2))
105 if ((i2 == nr2) || ((i1 < nr1) && (at1[i1] < at2[i2])))
113 if ((i2 < nr2) && ((i1 == nr1) || (at1[i1] > at2[i2])))
122 printf("Merged two groups with OR: %d %d -> %d\n", nr1, nr2, *nr);
128 static int and_groups(int nr1, const int* at1, int nr2, const int* at2, int* nr, int* at)
133 for (i1 = 0; i1 < nr1; i1++)
135 for (i2 = 0; i2 < nr2; i2++)
137 if (at1[i1] == at2[i2])
145 printf("Merged two groups with AND: %d %d -> %d\n", nr1, nr2, *nr);
150 static gmx_bool is_name_char(char c)
152 /* This string should contain all characters that can not be
153 * the first letter of a name due to the make_ndx syntax.
155 const char* spec = " !&|";
157 return (c != '\0' && std::strchr(spec, c) == nullptr);
160 static int parse_names(char** string, int* n_names, char** names)
165 while (is_name_char((*string)[0]) || (*string)[0] == ' ')
167 if (is_name_char((*string)[0]))
169 if (*n_names >= MAXNAMES)
171 gmx_fatal(FARGS, "To many names: %d\n", *n_names + 1);
174 while (is_name_char((*string)[i]))
176 names[*n_names][i] = (*string)[i];
180 printf("Name is too long, the maximum is %d characters\n", NAME_LEN);
184 names[*n_names][i] = '\0';
187 upstring(names[*n_names]);
201 static gmx_bool parse_int_char(char** string, int* nr, char* c)
208 while ((*string)[0] == ' ')
217 if (std::isdigit((*string)[0]))
219 *nr = (*string)[0] - '0';
221 while (std::isdigit((*string)[0]))
223 *nr = (*nr) * 10 + (*string)[0] - '0';
226 if (std::isalpha((*string)[0]))
231 /* Check if there is at most one non-digit character */
232 if (!std::isalnum((*string)[0]))
249 static gmx_bool parse_int(char** string, int* nr)
255 bRet = parse_int_char(string, nr, &c);
256 if (bRet && c != ' ')
265 static gmx_bool isquote(char c)
270 static gmx_bool parse_string(char** string, int* nr, int ngrps, char** grpname)
275 while ((*string)[0] == ' ')
281 if (isquote((*string)[0]))
285 s = gmx_strdup((*string));
286 sp = std::strchr(s, c);
289 (*string) += sp - s + 1;
291 (*nr) = find_group(s, ngrps, grpname);
295 return (*nr) != NOTSET;
298 static int select_atomnumbers(char** string, const t_atoms* atoms, int n1, 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);
353 } while ((*nr != 0) && (parse_int(string, &i)));
359 static int select_residuenumbers(char** string, const t_atoms* atoms, int n1, char c, int* nr, int* 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", *nr, (*nr == 1) ? "" : "s", n1, up);
396 sprintf(buf, "r_%d", n1);
400 sprintf(buf, "r_%d-%d", n1, up);
402 std::strcpy(gname, buf);
406 /* Individual residue number/insertion code selection */
411 for (i = 0; i < atoms->nr; i++)
413 ri = &atoms->resinfo[atoms->atom[i].resind];
414 if (ri->nr == j && ri->ic == c)
420 sprintf(buf, "_%d", j);
421 std::strcat(gname, buf);
422 } while (parse_int_char(string, &j, &c));
428 static int select_residueindices(char** string, const t_atoms* atoms, int n1, char c, int* nr, int* index, char* gname)
430 /*this should be similar to select_residuenumbers except select by index (sequential numbering in file)*/
431 /*resind+1 for 1-indexing*/
437 while ((*string)[0] == ' ')
441 if ((*string)[0] == '-')
443 /* Residue number range selection */
446 printf("Error: residue insertion codes can not be used with residue range selection\n");
450 parse_int(string, &up);
452 for (i = 0; i < atoms->nr; i++)
454 ri = &atoms->resinfo[atoms->atom[i].resind];
455 for (j = n1; (j <= up); j++)
457 if (atoms->atom[i].resind + 1 == j && (c == ' ' || ri->ic == c))
464 printf("Found %d atom%s with resind.+1 in range %d-%d\n", *nr, (*nr == 1) ? "" : "s", n1, up);
467 sprintf(buf, "r_%d", n1);
471 sprintf(buf, "r_%d-%d", n1, up);
473 std::strcpy(gname, buf);
477 /* Individual residue number/insertion code selection */
482 for (i = 0; i < atoms->nr; i++)
484 ri = &atoms->resinfo[atoms->atom[i].resind];
485 if (atoms->atom[i].resind + 1 == j && ri->ic == c)
491 sprintf(buf, "_%d", j);
492 std::strcat(gname, buf);
493 } while (parse_int_char(string, &j, &c));
501 atoms_from_residuenumbers(const t_atoms* atoms, int group, t_blocka* block, int* nr, int* index, char* gname)
503 int i, j, j0, j1, resnr, nres;
505 j0 = block->index[group];
506 j1 = block->index[group + 1];
508 for (j = j0; j < j1; j++)
510 if (block->a[j] >= nres)
512 printf("Index %s contains number>nres (%d>%d)\n", gname, block->a[j] + 1, nres);
516 for (i = 0; i < atoms->nr; i++)
518 resnr = atoms->resinfo[atoms->atom[i].resind].nr;
519 for (j = j0; j < j1; j++)
521 if (block->a[j] + 1 == resnr)
529 printf("Found %d atom%s in %d residues from group %s\n", *nr, (*nr == 1) ? "" : "s", j1 - j0, gname);
533 static gmx_bool comp_name(const char* name, const char* search)
535 gmx_bool matches = TRUE;
537 // Loop while name and search are not end-of-string and matches is true
538 for (; *name && *search && matches; name++, search++)
542 // still matching, continue to next character
545 else if (*search == '*')
549 printf("WARNING: Currently '*' is only supported at the end of an expression\n");
551 // if * is the last char in search string, we have a match,
552 // otherwise we just failed. Return in either case, we're done.
553 return (*(search + 1) == '\0');
556 // Compare one character
559 matches = (*name == *search);
563 matches = (std::toupper(*name) == std::toupper(*search));
567 matches = matches && (*name == '\0' && (*search == '\0' || *search == '*'));
572 static int select_chainnames(const t_atoms* atoms, int n_names, char** names, int* nr, int* index)
580 for (i = 0; i < atoms->nr; i++)
582 name[0] = atoms->resinfo[atoms->atom[i].resind].chainid;
584 while (j < n_names && !comp_name(name, names[j]))
594 printf("Found %d atom%s with chain identifier%s", *nr, (*nr == 1) ? "" : "s", (n_names == 1) ? "" : "s");
595 for (j = 0; (j < n_names); j++)
597 printf(" %s", names[j]);
604 static int select_atomnames(const t_atoms* atoms, int n_names, char** names, int* nr, int* index, gmx_bool bType)
611 for (i = 0; i < atoms->nr; i++)
615 name = *(atoms->atomtype[i]);
619 name = *(atoms->atomname[i]);
622 while (j < n_names && !comp_name(name, names[j]))
632 printf("Found %d atoms with %s%s", *nr, bType ? "type" : "name", (n_names == 1) ? "" : "s");
633 for (j = 0; (j < n_names); j++)
635 printf(" %s", names[j]);
642 static int select_residuenames(const t_atoms* atoms, int n_names, char** names, int* nr, int* index)
649 for (i = 0; i < atoms->nr; i++)
651 name = *(atoms->resinfo[atoms->atom[i].resind].name);
653 while (j < n_names && !comp_name(name, names[j]))
663 printf("Found %d atoms with residue name%s", *nr, (n_names == 1) ? "" : "s");
664 for (j = 0; (j < n_names); j++)
666 printf(" %s", names[j]);
673 static void copy2block(int n, const int* index, t_blocka* block)
680 srenew(block->index, block->nr + 1);
681 block->index[block->nr] = n0 + n;
682 srenew(block->a, n0 + n);
683 for (i = 0; (i < n); i++)
685 block->a[n0 + i] = index[i];
689 static void make_gname(int n, char** names, char* gname)
693 std::strcpy(gname, names[0]);
694 for (i = 1; i < n; i++)
696 std::strcat(gname, "_");
697 std::strcat(gname, names[i]);
701 static void copy_group(int g, t_blocka* block, int* nr, int* index)
705 i0 = block->index[g];
706 *nr = block->index[g + 1] - i0;
707 for (i = 0; i < *nr; i++)
709 index[i] = block->a[i0 + i];
713 static void remove_group(int nr, int nr2, t_blocka* block, char*** gn)
723 for (j = 0; j <= nr2 - nr; j++)
725 if ((nr < 0) || (nr >= block->nr))
727 printf("Group %d does not exist\n", nr + j);
731 shift = block->index[nr + 1] - block->index[nr];
732 for (i = block->index[nr + 1]; i < block->nra; i++)
734 block->a[i - shift] = block->a[i];
737 for (i = nr; i < block->nr; i++)
739 block->index[i] = block->index[i + 1] - shift;
741 name = gmx_strdup((*gn)[nr]);
743 for (i = nr; i < block->nr - 1; i++)
745 (*gn)[i] = (*gn)[i + 1];
748 block->nra = block->index[block->nr];
749 printf("Removed group %d '%s'\n", nr + j, name);
755 static void split_group(const t_atoms* atoms, int sel_nr, t_blocka* block, char*** gn, gmx_bool bAtom)
757 char buf[STRLEN], *name;
761 printf("Splitting group %d '%s' into %s\n", sel_nr, (*gn)[sel_nr], bAtom ? "atoms" : "residues");
763 n0 = block->index[sel_nr];
764 n1 = block->index[sel_nr + 1];
765 srenew(block->a, block->nra + n1 - n0);
766 for (i = n0; i < n1; i++)
769 resind = atoms->atom[a].resind;
770 name = *(atoms->resinfo[resind].name);
771 if (bAtom || (i == n0) || (atoms->atom[block->a[i - 1]].resind != resind))
775 block->index[block->nr] = block->nra;
778 srenew(block->index, block->nr + 1);
779 srenew(*gn, block->nr);
782 sprintf(buf, "%s_%s_%d", (*gn)[sel_nr], *atoms->atomname[a], a + 1);
786 sprintf(buf, "%s_%s_%d", (*gn)[sel_nr], name, atoms->resinfo[resind].nr);
788 (*gn)[block->nr - 1] = gmx_strdup(buf);
790 block->a[block->nra] = a;
793 block->index[block->nr] = block->nra;
796 static int split_chain(const t_atoms* atoms, const rvec* x, int sel_nr, t_blocka* block, char*** gn)
800 int i, a, natoms, *start = nullptr, *end = nullptr, ca_start, ca_end;
807 while (ca_start < natoms)
809 while ((ca_start < natoms) && std::strcmp(*atoms->atomname[ca_start], "CA") != 0)
813 if (ca_start < natoms)
815 srenew(start, nchain + 1);
816 srenew(end, nchain + 1);
817 start[nchain] = ca_start;
818 while ((start[nchain] > 0)
819 && (atoms->atom[start[nchain] - 1].resind == atoms->atom[ca_start].resind))
831 } while ((i < natoms) && std::strcmp(*atoms->atomname[i], "CA") != 0);
834 rvec_sub(x[ca_end], x[i], vec);
840 } while (norm(vec) < 0.45);
842 end[nchain] = ca_end;
843 while ((end[nchain] + 1 < natoms)
844 && (atoms->atom[end[nchain] + 1].resind == atoms->atom[ca_end].resind))
848 ca_start = end[nchain] + 1;
854 printf("Found 1 chain, will not split\n");
858 printf("Found %d chains\n", nchain);
860 for (j = 0; j < nchain; j++)
862 printf("%d:%6d atoms (%d to %d)\n", j + 1, end[j] - start[j] + 1, start[j] + 1, end[j] + 1);
867 srenew(block->a, block->nra + block->index[sel_nr + 1] - block->index[sel_nr]);
868 for (j = 0; j < nchain; j++)
871 srenew(block->index, block->nr + 1);
872 srenew(*gn, block->nr);
873 sprintf(buf, "%s_chain%d", (*gn)[sel_nr], j + 1);
874 (*gn)[block->nr - 1] = gmx_strdup(buf);
875 for (i = block->index[sel_nr]; i < block->index[sel_nr + 1]; i++)
878 if ((a >= start[j]) && (a <= end[j]))
880 block->a[block->nra] = a;
884 block->index[block->nr] = block->nra;
885 if (block->index[block->nr - 1] == block->index[block->nr])
887 remove_group(block->nr - 1, NOTSET, block, gn);
897 static gmx_bool check_have_atoms(const t_atoms* atoms, char* string)
899 if (atoms == nullptr)
901 printf("Can not process '%s' without atom info, use option -f\n", string);
910 static gmx_bool parse_entry(char** string,
912 const t_atoms* atoms,
919 static char ** names, *ostring;
920 static gmx_bool bFirst = TRUE;
921 int j, n_names, sel_nr1;
924 gmx_bool bRet, bCompl;
929 snew(names, MAXNAMES);
930 for (i = 0; i < MAXNAMES; i++)
932 snew(names[i], NAME_LEN + 1);
939 while (*string[0] == ' ')
944 if ((*string)[0] == '!')
948 while (*string[0] == ' ')
960 if (parse_int(string, &sel_nr1) || parse_string(string, &sel_nr1, block->nr, *gn))
962 if ((sel_nr1 >= 0) && (sel_nr1 < block->nr))
964 copy_group(sel_nr1, block, nr, index);
965 std::strcpy(gname, (*gn)[sel_nr1]);
966 printf("Copied index group %d '%s'\n", sel_nr1, (*gn)[sel_nr1]);
971 printf("Group %d does not exist\n", sel_nr1);
974 else if ((*string)[0] == 'a')
977 if (check_have_atoms(atoms, ostring))
979 if (parse_int(string, &sel_nr1))
981 bRet = (select_atomnumbers(string, atoms, sel_nr1, nr, index, gname) != 0);
983 else if (parse_names(string, &n_names, names))
985 bRet = (select_atomnames(atoms, n_names, names, nr, index, FALSE) != 0);
986 make_gname(n_names, names, gname);
990 else if ((*string)[0] == 't')
993 if (check_have_atoms(atoms, ostring) && parse_names(string, &n_names, names))
995 if (!(atoms->haveType))
997 printf("Need a run input file to select atom types\n");
1001 bRet = (select_atomnames(atoms, n_names, names, nr, index, TRUE) != 0);
1002 make_gname(n_names, names, gname);
1006 else if (std::strncmp(*string, "res", 3) == 0)
1009 if (check_have_atoms(atoms, ostring) && parse_int(string, &sel_nr1) && (sel_nr1 >= 0)
1010 && (sel_nr1 < block->nr))
1012 bRet = atoms_from_residuenumbers(atoms, sel_nr1, block, nr, index, (*gn)[sel_nr1]);
1013 sprintf(gname, "atom_%s", (*gn)[sel_nr1]);
1016 else if (std::strncmp(*string, "ri", 2) == 0)
1019 if (check_have_atoms(atoms, ostring) && parse_int_char(string, &sel_nr1, &c))
1021 bRet = (select_residueindices(string, atoms, sel_nr1, c, nr, index, gname) != 0);
1024 else if ((*string)[0] == 'r')
1027 if (check_have_atoms(atoms, ostring))
1029 if (parse_int_char(string, &sel_nr1, &c))
1031 bRet = (select_residuenumbers(string, atoms, sel_nr1, c, nr, index, gname) != 0);
1033 else if (parse_names(string, &n_names, names))
1035 bRet = (select_residuenames(atoms, n_names, names, nr, index) != 0);
1036 make_gname(n_names, names, gname);
1040 else if (std::strncmp(*string, "chain", 5) == 0)
1043 if (check_have_atoms(atoms, ostring) && parse_names(string, &n_names, names))
1045 bRet = (select_chainnames(atoms, n_names, names, nr, index) != 0);
1046 sprintf(gname, "ch%s", names[0]);
1047 for (i = 1; i < n_names; i++)
1049 std::strcat(gname, names[i]);
1055 snew(index1, natoms - *nr);
1057 for (i = 0; i < natoms; i++)
1060 while ((j < *nr) && (index[j] != i))
1066 if (nr1 >= natoms - *nr)
1068 printf("There are double atoms in your index group\n");
1076 for (i = 0; i < nr1; i++)
1078 index[i] = index1[i];
1082 for (i = std::strlen(gname) + 1; i > 0; i--)
1084 gname[i] = gname[i - 1];
1087 printf("Complemented group: %d atoms\n", *nr);
1093 static void list_residues(const t_atoms* atoms)
1095 int i, j, start, end, prev_resind, resind;
1098 /* Print all the residues, assuming continuous resnr count */
1099 start = atoms->atom[0].resind;
1100 prev_resind = start;
1101 for (i = 0; i < atoms->nr; i++)
1103 resind = atoms->atom[i].resind;
1104 if ((resind != prev_resind) || (i == atoms->nr - 1))
1106 if ((bDiff = (std::strcmp(*atoms->resinfo[resind].name, *atoms->resinfo[start].name) != 0))
1107 || (i == atoms->nr - 1))
1117 if (end < start + 3)
1119 for (j = start; j <= end; j++)
1121 printf("%4d %-5s", j + 1, *(atoms->resinfo[j].name));
1126 printf(" %4d - %4d %-5s ", start + 1, end + 1, *(atoms->resinfo[start].name));
1131 prev_resind = resind;
1136 static void edit_index(int natoms, const t_atoms* atoms, const rvec* x, t_blocka* block, char*** gn, gmx_bool bVerbose)
1138 static char ** atnames, *ostring;
1139 static gmx_bool bFirst = TRUE;
1140 char inp_string[STRLEN], *string;
1141 char gname[STRLEN * 3], gname1[STRLEN], gname2[STRLEN];
1142 int i, i0, i1, sel_nr, sel_nr2, newgroup;
1143 int nr, nr1, nr2, *index, *index1, *index2;
1144 gmx_bool bAnd, bOr, bPrintOnce;
1149 snew(atnames, MAXNAMES);
1150 for (i = 0; i < MAXNAMES; i++)
1152 snew(atnames[i], NAME_LEN + 1);
1158 snew(index, natoms);
1159 snew(index1, natoms);
1160 snew(index2, natoms);
1167 if (bVerbose || bPrintOnce || newgroup != NOTSET)
1170 if (bVerbose || bPrintOnce || newgroup == NOTSET)
1180 for (i = i0; i < i1; i++)
1182 printf("%3d %-20s: %5d atoms\n", i, (*gn)[i], block->index[i + 1] - block->index[i]);
1186 if (bVerbose || bPrintOnce)
1189 printf(" nr : group '!': not 'name' nr name 'splitch' nr Enter: list "
1191 printf(" 'a': atom '&': and 'del' nr 'splitres' nr 'l': list "
1193 printf(" 't': atom type '|': or 'keep' nr 'splitat' nr 'h': help\n");
1194 printf(" 'r': residue 'res' nr 'chain' char\n");
1195 printf(" \"name\": group 'case': case %s 'q': save and quit\n",
1196 bCase ? "insensitive" : "sensitive ");
1197 printf(" 'ri': residue index\n");
1202 if (nullptr == fgets(inp_string, STRLEN, stdin))
1204 gmx_fatal(FARGS, "Error reading user input");
1206 inp_string[std::strlen(inp_string) - 1] = 0;
1208 string = inp_string;
1209 while (string[0] == ' ')
1216 if (string[0] == 'h')
1218 printf(" nr : selects an index group by number or quoted string.\n");
1219 printf(" The string is first matched against the whole group "
1221 printf(" then against the beginning and finally against an\n");
1222 printf(" arbitrary substring. A multiple match is an error.\n");
1224 printf(" 'a' nr1 [nr2 ...] : selects atoms, atom numbering starts at 1.\n");
1225 printf(" 'a' nr1 - nr2 : selects atoms in the range from nr1 to nr2.\n");
1226 printf(" 'a' name1[*] [name2[*] ...] : selects atoms by name(s), '?' matches any "
1228 printf(" wildcard '*' allowed at the end of a name.\n");
1229 printf(" 't' type1[*] [type2[*] ...] : as 'a', but for type, run input file "
1231 printf(" 'r' nr1[ic1] [nr2[ic2] ...] : selects residues by number and insertion "
1233 printf(" 'r' nr1 - nr2 : selects residues in the range from nr1 to "
1235 printf(" 'r' name1[*] [name2[*] ...] : as 'a', but for residue names.\n");
1236 printf(" 'ri' nr1 - nr2 : selects residue indices, 1-indexed, (as opposed "
1237 "to numbers) in the range from nr1 to nr2.\n");
1238 printf(" 'chain' ch1 [ch2 ...] : selects atoms by chain identifier(s),\n");
1239 printf(" not available with a .gro file as input.\n");
1240 printf(" ! : takes the complement of a group with respect to all\n");
1241 printf(" the atoms in the input file.\n");
1242 printf(" & | : AND and OR, can be placed between any of the options\n");
1243 printf(" above, the input is processed from left to right.\n");
1244 printf(" 'name' nr name : rename group nr to name.\n");
1245 printf(" 'del' nr1 [- nr2] : deletes one group or groups in the range from nr1 to "
1247 printf(" 'keep' nr : deletes all groups except nr.\n");
1248 printf(" 'case' : make all name compares case (in)sensitive.\n");
1249 printf(" 'splitch' nr : split group into chains using CA distances.\n");
1250 printf(" 'splitres' nr : split group into residues.\n");
1251 printf(" 'splitat' nr : split group into atoms.\n");
1252 printf(" 'res' nr : interpret numbers in group as residue numbers\n");
1253 printf(" Enter : list the currently defined groups and commands\n");
1254 printf(" 'l' : list the residues.\n");
1255 printf(" 'h' : show this help.\n");
1256 printf(" 'q' : save and quit.\n");
1258 printf(" Examples:\n");
1259 printf(" > 2 | 4 & r 3-5\n");
1260 printf(" selects all atoms from group 2 and 4 that have residue numbers 3, 4 or 5\n");
1261 printf(" > a C* & !a C CA\n");
1262 printf(" selects all atoms starting with 'C' but not the atoms 'C' and 'CA'\n");
1263 printf(" > \"protein\" & ! \"backb\"\n");
1264 printf(" selects all atoms that are in group 'protein' and not in group 'backbone'\n");
1267 printf("\npress Enter ");
1271 else if (std::strncmp(string, "del", 3) == 0)
1274 if (parse_int(&string, &sel_nr))
1276 while (string[0] == ' ')
1280 if (string[0] == '-')
1283 parse_int(&string, &sel_nr2);
1289 while (string[0] == ' ')
1293 if (string[0] == '\0')
1295 remove_group(sel_nr, sel_nr2, block, gn);
1299 printf("\nSyntax error: \"%s\"\n", string);
1303 else if (std::strncmp(string, "keep", 4) == 0)
1306 if (parse_int(&string, &sel_nr))
1308 remove_group(sel_nr + 1, block->nr - 1, block, gn);
1309 remove_group(0, sel_nr - 1, block, gn);
1312 else if (std::strncmp(string, "name", 4) == 0)
1315 if (parse_int(&string, &sel_nr))
1317 if ((sel_nr >= 0) && (sel_nr < block->nr))
1319 sscanf(string, "%s", gname);
1320 sfree((*gn)[sel_nr]);
1321 (*gn)[sel_nr] = gmx_strdup(gname);
1325 else if (std::strncmp(string, "case", 4) == 0)
1328 printf("Switched to case %s\n", bCase ? "sensitive" : "insensitive");
1330 else if (string[0] == 'v')
1332 bVerbose = !bVerbose;
1333 printf("Turned verbose %s\n", bVerbose ? "on" : "off");
1335 else if (string[0] == 'l')
1337 if (check_have_atoms(atoms, ostring))
1339 list_residues(atoms);
1342 else if (std::strncmp(string, "splitch", 7) == 0)
1345 if (check_have_atoms(atoms, ostring) && parse_int(&string, &sel_nr) && (sel_nr >= 0)
1346 && (sel_nr < block->nr))
1348 split_chain(atoms, x, sel_nr, block, gn);
1351 else if (std::strncmp(string, "splitres", 8) == 0)
1354 if (check_have_atoms(atoms, ostring) && parse_int(&string, &sel_nr) && (sel_nr >= 0)
1355 && (sel_nr < block->nr))
1357 split_group(atoms, sel_nr, block, gn, FALSE);
1360 else if (std::strncmp(string, "splitat", 7) == 0)
1363 if (check_have_atoms(atoms, ostring) && parse_int(&string, &sel_nr) && (sel_nr >= 0)
1364 && (sel_nr < block->nr))
1366 split_group(atoms, sel_nr, block, gn, TRUE);
1369 else if (string[0] == '\0')
1373 else if (string[0] != 'q')
1376 if (parse_entry(&string, natoms, atoms, block, gn, &nr, index, gname))
1380 while (string[0] == ' ')
1387 if (string[0] == '&')
1391 else if (string[0] == '|')
1400 for (i = 0; i < nr; i++)
1402 index1[i] = index[i];
1404 std::strcpy(gname1, gname);
1405 if (parse_entry(&string, natoms, atoms, block, gn, &nr2, index2, gname2))
1409 or_groups(nr1, index1, nr2, index2, &nr, index);
1410 sprintf(gname, "%s_%s", gname1, gname2);
1414 and_groups(nr1, index1, nr2, index2, &nr, index);
1415 sprintf(gname, "%s_&_%s", gname1, gname2);
1419 } while (bAnd || bOr);
1421 while (string[0] == ' ')
1427 printf("\nSyntax error: \"%s\"\n", string);
1431 copy2block(nr, index, block);
1432 srenew(*gn, block->nr);
1433 newgroup = block->nr - 1;
1434 (*gn)[newgroup] = gmx_strdup(gname);
1438 printf("Group is empty\n");
1441 } while (string[0] != 'q');
1448 static int block2natoms(t_blocka* block)
1453 for (i = 0; i < block->nra; i++)
1455 natoms = std::max(natoms, block->a[i]);
1462 static void merge_blocks(t_blocka* dest, t_blocka* source)
1466 /* count groups, srenew and fill */
1469 dest->nr += source->nr;
1470 srenew(dest->index, dest->nr + 1);
1471 for (i = 0; i < source->nr; i++)
1473 dest->index[i0 + i] = nra0 + source->index[i];
1475 /* count atoms, srenew and fill */
1476 dest->nra += source->nra;
1477 srenew(dest->a, dest->nra);
1478 for (i = 0; i < source->nra; i++)
1480 dest->a[nra0 + i] = source->a[i];
1483 /* terminate list */
1484 dest->index[dest->nr] = dest->nra;
1487 int gmx_make_ndx(int argc, char* argv[])
1489 const char* desc[] = { "Index groups are necessary for almost every GROMACS program.",
1490 "All these programs can generate default index groups. You ONLY",
1491 "have to use [THISMODULE] when you need SPECIAL index groups.",
1492 "There is a default index group for the whole system, 9 default",
1493 "index groups for proteins, and a default index group",
1494 "is generated for every other residue name.",
1496 "When no index file is supplied, also [THISMODULE] will generate the",
1498 "With the index editor you can select on atom, residue and chain names",
1500 "When a run input file is supplied you can also select on atom type.",
1501 "You can use boolean operations, you can split groups",
1502 "into chains, residues or atoms. You can delete and rename groups.",
1503 "Type 'h' in the editor for more details.",
1505 "The atom numbering in the editor and the index file starts at 1.",
1507 "The [TT]-twin[tt] switch duplicates all index groups with an offset of",
1508 "[TT]-natoms[tt], which is useful for Computational Electrophysiology",
1509 "double-layer membrane setups.",
1511 "See also [gmx-select] [TT]-on[tt], which provides an alternative way",
1512 "for constructing index groups. It covers nearly all of [THISMODULE]",
1513 "functionality, and in many cases much more." };
1515 static int natoms = 0;
1516 static gmx_bool bVerbose = FALSE;
1517 static gmx_bool bDuplicate = FALSE;
1518 t_pargs pa[] = { { "-natoms",
1522 "set number of atoms (default: read from coordinate or index file)" },
1527 "Duplicate all index groups with an offset of -natoms" },
1528 { "-verbose", FALSE, etBOOL, { &bVerbose }, "HIDDENVerbose output" } };
1529 #define NPA asize(pa)
1531 gmx_output_env_t* oenv;
1532 const char* stxfile;
1533 const char* ndxoutfile;
1540 t_blocka * block, *block2;
1541 char ** gnames, **gnames2;
1542 t_filenm fnm[] = { { efSTX, "-f", nullptr, ffOPTRD },
1543 { efNDX, "-n", nullptr, ffOPTRDMULT },
1544 { efNDX, "-o", nullptr, ffWRITE } };
1545 #define NFILE asize(fnm)
1547 if (!parse_common_args(&argc, argv, 0, NFILE, fnm, NPA, pa, asize(desc), desc, 0, nullptr, &oenv))
1552 stxfile = ftp2fn_null(efSTX, NFILE, fnm);
1553 gmx::ArrayRef<const std::string> ndxInFiles = opt2fnsIfOptionSet("-n", NFILE, fnm);
1554 ndxoutfile = opt2fn("-o", NFILE, fnm);
1555 bNatoms = opt2parg_bSet("-natoms", NPA, pa);
1557 if (!stxfile && ndxInFiles.empty())
1559 gmx_fatal(FARGS, "No input files (structure or index)");
1565 bool haveFullTopology = false;
1566 fprintf(stderr, "\nReading structure file\n");
1567 readConfAndTopology(stxfile, &haveFullTopology, &mtop, &pbcType, &x, &v, box);
1568 atoms = gmx_mtop_global_atoms(&mtop);
1569 if (atoms.pdbinfo == nullptr)
1571 snew(atoms.pdbinfo, atoms.nr);
1581 /* read input file(s) */
1582 block = new_blocka();
1584 printf("Going to read %td old index file(s)\n", ndxInFiles.ssize());
1585 if (!ndxInFiles.empty())
1587 for (const std::string& ndxInFile : ndxInFiles)
1589 block2 = init_index(ndxInFile.c_str(), &gnames2);
1590 srenew(gnames, block->nr + block2->nr);
1591 for (j = 0; j < block2->nr; j++)
1593 gnames[block->nr + j] = gnames2[j];
1596 merge_blocks(block, block2);
1598 sfree(block2->index);
1599 /* done_block(block2); */
1606 analyse(&atoms, block, &gnames, FALSE, TRUE);
1611 natoms = block2natoms(block);
1612 printf("Counted atom numbers up to %d in index file\n", natoms);
1615 edit_index(natoms, &atoms, x, block, &gnames, bVerbose);
1617 write_index(ndxoutfile, block, gnames, bDuplicate, natoms);