2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5 * Copyright (c) 2001-2004, The GROMACS development team.
6 * Copyright (c) 2012,2013,2014, by the GROMACS development team, led by
7 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
8 * and including many others, as listed in the AUTHORS file in the
9 * top-level source directory and at http://www.gromacs.org.
11 * GROMACS is free software; you can redistribute it and/or
12 * modify it under the terms of the GNU Lesser General Public License
13 * as published by the Free Software Foundation; either version 2.1
14 * of the License, or (at your option) any later version.
16 * GROMACS is distributed in the hope that it will be useful,
17 * but WITHOUT ANY WARRANTY; without even the implied warranty of
18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19 * Lesser General Public License for more details.
21 * You should have received a copy of the GNU Lesser General Public
22 * License along with GROMACS; if not, see
23 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
24 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
26 * If you want to redistribute modifications to GROMACS, please
27 * consider that scientific software is very special. Version
28 * control is crucial - bugs must be traceable. We will be happy to
29 * consider code for inclusion in the official distribution, but
30 * derived work must not be called official GROMACS. Details are found
31 * in the README & COPYING files - if they are missing, get the
32 * official version at http://www.gromacs.org.
34 * To help us fund GROMACS development, we humbly ask that you cite
35 * the research papers on the package. Check out http://www.gromacs.org.
45 #include "gromacs/fileio/futil.h"
47 #include "gromacs/utility/cstringutil.h"
48 #include "gromacs/commandline/pargs.h"
49 #include "gromacs/fileio/confio.h"
52 #include "gromacs/utility/smalloc.h"
56 #include "gromacs/legacyheaders/gmx_fatal.h"
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 = 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 = 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] = 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);
851 while ((i < natoms) && (norm(vec) < 0.45));
853 end[nchain] = ca_end;
854 while ((end[nchain]+1 < natoms) &&
855 (atoms->atom[end[nchain]+1].resind == atoms->atom[ca_end].resind))
859 ca_start = end[nchain]+1;
865 printf("Found 1 chain, will not split\n");
869 printf("Found %d chains\n", nchain);
871 for (j = 0; j < nchain; j++)
873 printf("%d:%6d atoms (%d to %d)\n",
874 j+1, end[j]-start[j]+1, start[j]+1, end[j]+1);
879 srenew(block->a, block->nra+block->index[sel_nr+1]-block->index[sel_nr]);
880 for (j = 0; j < nchain; j++)
883 srenew(block->index, block->nr+1);
884 srenew(*gn, block->nr);
885 sprintf(buf, "%s_chain%d", (*gn)[sel_nr], j+1);
886 (*gn)[block->nr-1] = strdup(buf);
887 for (i = block->index[sel_nr]; i < block->index[sel_nr+1]; i++)
890 if ((a >= start[j]) && (a <= end[j]))
892 block->a[block->nra] = a;
896 block->index[block->nr] = block->nra;
897 if (block->index[block->nr-1] == block->index[block->nr])
899 remove_group(block->nr-1, NOTSET, block, gn);
909 static gmx_bool check_have_atoms(t_atoms *atoms, char *string)
913 printf("Can not process '%s' without atom info, use option -f\n", string);
922 static gmx_bool parse_entry(char **string, int natoms, t_atoms *atoms,
923 t_blocka *block, char ***gn,
924 atom_id *nr, atom_id *index, char *gname)
926 static char **names, *ostring;
927 static gmx_bool bFirst = TRUE;
928 int j, n_names, sel_nr1;
929 atom_id i, nr1, *index1;
931 gmx_bool bRet, bCompl;
936 snew(names, MAXNAMES);
937 for (i = 0; i < MAXNAMES; i++)
939 snew(names[i], NAME_LEN+1);
946 while (*string[0] == ' ')
951 if ((*string)[0] == '!')
955 while (*string[0] == ' ')
967 if (parse_int(string, &sel_nr1) ||
968 parse_string(string, &sel_nr1, block->nr, *gn))
970 if ((sel_nr1 >= 0) && (sel_nr1 < block->nr))
972 copy_group(sel_nr1, block, nr, index);
973 strcpy(gname, (*gn)[sel_nr1]);
974 printf("Copied index group %d '%s'\n", sel_nr1, (*gn)[sel_nr1]);
979 printf("Group %d does not exist\n", sel_nr1);
982 else if ((*string)[0] == 'a')
985 if (check_have_atoms(atoms, ostring))
987 if (parse_int(string, &sel_nr1))
989 bRet = select_atomnumbers(string, atoms, sel_nr1, nr, index, gname);
991 else if (parse_names(string, &n_names, names))
993 bRet = select_atomnames(atoms, n_names, names, nr, index, FALSE);
994 make_gname(n_names, names, gname);
998 else if ((*string)[0] == 't')
1001 if (check_have_atoms(atoms, ostring) &&
1002 parse_names(string, &n_names, names))
1004 if (atoms->atomtype == NULL)
1006 printf("Need a run input file to select atom types\n");
1010 bRet = select_atomnames(atoms, n_names, names, nr, index, TRUE);
1011 make_gname(n_names, names, gname);
1015 else if (strncmp(*string, "res", 3) == 0)
1018 if (check_have_atoms(atoms, ostring) &&
1019 parse_int(string, &sel_nr1) &&
1020 (sel_nr1 >= 0) && (sel_nr1 < block->nr) )
1022 bRet = atoms_from_residuenumbers(atoms,
1023 sel_nr1, block, nr, index, (*gn)[sel_nr1]);
1024 sprintf(gname, "atom_%s", (*gn)[sel_nr1]);
1027 else if (strncmp(*string, "ri", 2) == 0)
1030 if (check_have_atoms(atoms, ostring) &&
1031 parse_int_char(string, &sel_nr1, &c))
1033 bRet = select_residueindices(string, atoms, sel_nr1, c, nr, index, gname);
1036 else if ((*string)[0] == 'r')
1039 if (check_have_atoms(atoms, ostring))
1041 if (parse_int_char(string, &sel_nr1, &c))
1043 bRet = select_residuenumbers(string, atoms, sel_nr1, c, nr, index, gname);
1045 else if (parse_names(string, &n_names, names))
1047 bRet = select_residuenames(atoms, n_names, names, nr, index);
1048 make_gname(n_names, names, gname);
1052 else if (strncmp(*string, "chain", 5) == 0)
1055 if (check_have_atoms(atoms, ostring) &&
1056 parse_names(string, &n_names, names))
1058 bRet = select_chainnames(atoms, n_names, names, nr, index);
1059 sprintf(gname, "ch%s", names[0]);
1060 for (i = 1; i < n_names; i++)
1062 strcat(gname, names[i]);
1068 snew(index1, natoms-*nr);
1070 for (i = 0; i < natoms; i++)
1073 while ((j < *nr) && (index[j] != i))
1079 if (nr1 >= natoms-*nr)
1081 printf("There are double atoms in your index group\n");
1089 for (i = 0; i < nr1; i++)
1091 index[i] = index1[i];
1095 for (i = strlen(gname)+1; i > 0; i--)
1097 gname[i] = gname[i-1];
1100 printf("Complemented group: %d atoms\n", *nr);
1106 static void list_residues(t_atoms *atoms)
1108 int i, j, start, end, prev_resind, resind;
1111 /* Print all the residues, assuming continuous resnr count */
1112 start = atoms->atom[0].resind;
1113 prev_resind = start;
1114 for (i = 0; i < atoms->nr; i++)
1116 resind = atoms->atom[i].resind;
1117 if ((resind != prev_resind) || (i == atoms->nr-1))
1119 if ((bDiff = strcmp(*atoms->resinfo[resind].name,
1120 *atoms->resinfo[start].name)) ||
1133 for (j = start; j <= end; j++)
1136 j+1, *(atoms->resinfo[j].name));
1141 printf(" %4d - %4d %-5s ",
1142 start+1, end+1, *(atoms->resinfo[start].name));
1147 prev_resind = resind;
1152 static void edit_index(int natoms, t_atoms *atoms, rvec *x, t_blocka *block, char ***gn, gmx_bool bVerbose)
1154 static char **atnames, *ostring;
1155 static gmx_bool bFirst = TRUE;
1156 char inp_string[STRLEN], *string;
1157 char gname[STRLEN], gname1[STRLEN], gname2[STRLEN];
1158 int i, i0, i1, sel_nr, sel_nr2, newgroup;
1159 atom_id nr, nr1, nr2, *index, *index1, *index2;
1160 gmx_bool bAnd, bOr, bPrintOnce;
1165 snew(atnames, MAXNAMES);
1166 for (i = 0; i < MAXNAMES; i++)
1168 snew(atnames[i], NAME_LEN+1);
1174 snew(index, natoms);
1175 snew(index1, natoms);
1176 snew(index2, natoms);
1183 if (bVerbose || bPrintOnce || newgroup != NOTSET)
1186 if (bVerbose || bPrintOnce || newgroup == NOTSET)
1196 for (i = i0; i < i1; i++)
1198 printf("%3d %-20s: %5d atoms\n", i, (*gn)[i],
1199 block->index[i+1]-block->index[i]);
1203 if (bVerbose || bPrintOnce)
1206 printf(" nr : group ! 'name' nr name 'splitch' nr Enter: list groups\n");
1207 printf(" 'a': atom & 'del' nr 'splitres' nr 'l': list residues\n");
1208 printf(" 't': atom type | 'keep' nr 'splitat' nr 'h': help\n");
1209 printf(" 'r': residue 'res' nr 'chain' char\n");
1210 printf(" \"name\": group 'case': case %s 'q': save and quit\n",
1211 bCase ? "insensitive" : "sensitive ");
1212 printf(" 'ri': residue index\n");
1217 if (NULL == fgets(inp_string, STRLEN, stdin))
1219 gmx_fatal(FARGS, "Error reading user input");
1221 inp_string[strlen(inp_string)-1] = 0;
1223 string = inp_string;
1224 while (string[0] == ' ')
1231 if (string[0] == 'h')
1233 printf(" nr : selects an index group by number or quoted string.\n");
1234 printf(" The string is first matched against the whole group name,\n");
1235 printf(" then against the beginning and finally against an\n");
1236 printf(" arbitrary substring. A multiple match is an error.\n");
1238 printf(" 'a' nr1 [nr2 ...] : selects atoms, atom numbering starts at 1.\n");
1239 printf(" 'a' nr1 - nr2 : selects atoms in the range from nr1 to nr2.\n");
1240 printf(" 'a' name1[*] [name2[*] ...] : selects atoms by name(s), '?' matches any char,\n");
1241 printf(" wildcard '*' allowed at the end of a name.\n");
1242 printf(" 't' type1[*] [type2[*] ...] : as 'a', but for type, run input file required.\n");
1243 printf(" 'r' nr1[ic1] [nr2[ic2] ...] : selects residues by number and insertion code.\n");
1244 printf(" 'r' nr1 - nr2 : selects residues in the range from nr1 to nr2.\n");
1245 printf(" 'r' name1[*] [name2[*] ...] : as 'a', but for residue names.\n");
1246 printf(" 'ri' nr1 - nr2 : selects residue indices, 1-indexed, (as opposed to numbers) in the range from nr1 to nr2.\n");
1247 printf(" 'chain' ch1 [ch2 ...] : selects atoms by chain identifier(s),\n");
1248 printf(" not available with a .gro file as input.\n");
1249 printf(" ! : takes the complement of a group with respect to all\n");
1250 printf(" the atoms in the input file.\n");
1251 printf(" & | : AND and OR, can be placed between any of the options\n");
1252 printf(" above, the input is processed from left to right.\n");
1253 printf(" 'name' nr name : rename group nr to name.\n");
1254 printf(" 'del' nr1 [- nr2] : deletes one group or groups in the range from nr1 to nr2.\n");
1255 printf(" 'keep' nr : deletes all groups except nr.\n");
1256 printf(" 'case' : make all name compares case (in)sensitive.\n");
1257 printf(" 'splitch' nr : split group into chains using CA distances.\n");
1258 printf(" 'splitres' nr : split group into residues.\n");
1259 printf(" 'splitat' nr : split group into atoms.\n");
1260 printf(" 'res' nr : interpret numbers in group as residue numbers\n");
1261 printf(" Enter : list the currently defined groups and commands\n");
1262 printf(" 'l' : list the residues.\n");
1263 printf(" 'h' : show this help.\n");
1264 printf(" 'q' : save and quit.\n");
1266 printf(" Examples:\n");
1267 printf(" > 2 | 4 & r 3-5\n");
1268 printf(" selects all atoms from group 2 and 4 that have residue numbers 3, 4 or 5\n");
1269 printf(" > a C* & !a C CA\n");
1270 printf(" selects all atoms starting with 'C' but not the atoms 'C' and 'CA'\n");
1271 printf(" > \"protein\" & ! \"backb\"\n");
1272 printf(" selects all atoms that are in group 'protein' and not in group 'backbone'\n");
1275 printf("\npress Enter ");
1279 else if (strncmp(string, "del", 3) == 0)
1282 if (parse_int(&string, &sel_nr))
1284 while (string[0] == ' ')
1288 if (string[0] == '-')
1291 parse_int(&string, &sel_nr2);
1297 while (string[0] == ' ')
1301 if (string[0] == '\0')
1303 remove_group(sel_nr, sel_nr2, block, gn);
1307 printf("\nSyntax error: \"%s\"\n", string);
1311 else if (strncmp(string, "keep", 4) == 0)
1314 if (parse_int(&string, &sel_nr))
1316 remove_group(sel_nr+1, block->nr-1, block, gn);
1317 remove_group(0, sel_nr-1, block, gn);
1320 else if (strncmp(string, "name", 4) == 0)
1323 if (parse_int(&string, &sel_nr))
1325 if ((sel_nr >= 0) && (sel_nr < block->nr))
1327 sscanf(string, "%s", gname);
1328 sfree((*gn)[sel_nr]);
1329 (*gn)[sel_nr] = strdup(gname);
1333 else if (strncmp(string, "case", 4) == 0)
1336 printf("Switched to case %s\n", bCase ? "sensitive" : "insensitive");
1338 else if (string[0] == 'v')
1340 bVerbose = !bVerbose;
1341 printf("Turned verbose %s\n", bVerbose ? "on" : "off");
1343 else if (string[0] == 'l')
1345 if (check_have_atoms(atoms, ostring) )
1347 list_residues(atoms);
1350 else if (strncmp(string, "splitch", 7) == 0)
1353 if (check_have_atoms(atoms, ostring) &&
1354 parse_int(&string, &sel_nr) &&
1355 (sel_nr >= 0) && (sel_nr < block->nr))
1357 split_chain(atoms, x, sel_nr, block, gn);
1360 else if (strncmp(string, "splitres", 8) == 0)
1363 if (check_have_atoms(atoms, ostring) &&
1364 parse_int(&string, &sel_nr) &&
1365 (sel_nr >= 0) && (sel_nr < block->nr))
1367 split_group(atoms, sel_nr, block, gn, FALSE);
1370 else if (strncmp(string, "splitat", 7) == 0)
1373 if (check_have_atoms(atoms, ostring) &&
1374 parse_int(&string, &sel_nr) &&
1375 (sel_nr >= 0) && (sel_nr < block->nr))
1377 split_group(atoms, sel_nr, block, gn, TRUE);
1380 else if (string[0] == '\0')
1384 else if (string[0] != 'q')
1388 if (parse_entry(&string, natoms, atoms, block, gn, &nr, index, gname))
1392 while (string[0] == ' ')
1399 if (string[0] == '&')
1403 else if (string[0] == '|')
1412 for (i = 0; i < nr; i++)
1414 index1[i] = index[i];
1416 strcpy(gname1, gname);
1417 if (parse_entry(&string, natoms, atoms, block, gn, &nr2, index2, gname2))
1421 or_groups(nr1, index1, nr2, index2, &nr, index);
1422 sprintf(gname, "%s_%s", gname1, gname2);
1426 and_groups(nr1, index1, nr2, index2, &nr, index);
1427 sprintf(gname, "%s_&_%s", gname1, gname2);
1432 while (bAnd || bOr);
1434 while (string[0] == ' ')
1440 printf("\nSyntax error: \"%s\"\n", string);
1444 copy2block(nr, index, block);
1445 srenew(*gn, block->nr);
1446 newgroup = block->nr-1;
1447 (*gn)[newgroup] = strdup(gname);
1451 printf("Group is empty\n");
1455 while (string[0] != 'q');
1462 static int block2natoms(t_blocka *block)
1467 for (i = 0; i < block->nra; i++)
1469 natoms = max(natoms, block->a[i]);
1476 void merge_blocks(t_blocka *dest, t_blocka *source)
1480 /* count groups, srenew and fill */
1483 dest->nr += source->nr;
1484 srenew(dest->index, dest->nr+1);
1485 for (i = 0; i < source->nr; i++)
1487 dest->index[i0+i] = nra0 + source->index[i];
1489 /* count atoms, srenew and fill */
1490 dest->nra += source->nra;
1491 srenew(dest->a, dest->nra);
1492 for (i = 0; i < source->nra; i++)
1494 dest->a[nra0+i] = source->a[i];
1497 /* terminate list */
1498 dest->index[dest->nr] = dest->nra;
1502 int gmx_make_ndx(int argc, char *argv[])
1504 const char *desc[] = {
1505 "Index groups are necessary for almost every GROMACS program.",
1506 "All these programs can generate default index groups. You ONLY",
1507 "have to use [THISMODULE] when you need SPECIAL index groups.",
1508 "There is a default index group for the whole system, 9 default",
1509 "index groups for proteins, and a default index group",
1510 "is generated for every other residue name.[PAR]",
1511 "When no index file is supplied, also [THISMODULE] will generate the",
1513 "With the index editor you can select on atom, residue and chain names",
1515 "When a run input file is supplied you can also select on atom type.",
1516 "You can use NOT, AND and OR, you can split groups",
1517 "into chains, residues or atoms. You can delete and rename groups.[PAR]",
1518 "The atom numbering in the editor and the index file starts at 1.[PAR]",
1519 "The [TT]-twin[tt] switch duplicates all index groups with an offset of",
1520 "[TT]-natoms[tt], which is useful for Computational Electrophysiology",
1521 "double-layer membrane setups."
1524 static int natoms = 0;
1525 static gmx_bool bVerbose = FALSE;
1526 static gmx_bool bDuplicate = FALSE;
1528 { "-natoms", FALSE, etINT, {&natoms},
1529 "set number of atoms (default: read from coordinate or index file)" },
1530 { "-twin", FALSE, etBOOL, {&bDuplicate},
1531 "Duplicate all index groups with an offset of -natoms" },
1532 { "-verbose", FALSE, etBOOL, {&bVerbose},
1533 "HIDDENVerbose output" }
1535 #define NPA asize(pa)
1540 const char *stxfile;
1542 const char *ndxoutfile;
1549 t_blocka *block, *block2;
1550 char **gnames, **gnames2;
1552 { efSTX, "-f", NULL, ffOPTRD },
1553 { efNDX, "-n", NULL, ffOPTRDMULT },
1554 { efNDX, "-o", NULL, ffWRITE }
1556 #define NFILE asize(fnm)
1558 if (!parse_common_args(&argc, argv, 0, NFILE, fnm, NPA, pa, asize(desc), desc,
1564 stxfile = ftp2fn_null(efSTX, NFILE, fnm);
1565 if (opt2bSet("-n", NFILE, fnm))
1567 nndxin = opt2fns(&ndxinfiles, "-n", NFILE, fnm);
1573 ndxoutfile = opt2fn("-o", NFILE, fnm);
1574 bNatoms = opt2parg_bSet("-natoms", NPA, pa);
1576 if (!stxfile && !nndxin)
1578 gmx_fatal(FARGS, "No input files (structure or index)");
1584 get_stx_coordnum(stxfile, &(atoms->nr));
1585 init_t_atoms(atoms, atoms->nr, TRUE);
1588 fprintf(stderr, "\nReading structure file\n");
1589 read_stx_conf(stxfile, title, atoms, x, v, &ePBC, box);
1599 /* read input file(s) */
1600 block = new_blocka();
1602 printf("Going to read %d old index file(s)\n", nndxin);
1605 for (i = 0; i < nndxin; i++)
1607 block2 = init_index(ndxinfiles[i], &gnames2);
1608 srenew(gnames, block->nr+block2->nr);
1609 for (j = 0; j < block2->nr; j++)
1611 gnames[block->nr+j] = gnames2[j];
1614 merge_blocks(block, block2);
1616 sfree(block2->index);
1617 /* done_block(block2); */
1624 analyse(atoms, block, &gnames, FALSE, TRUE);
1629 natoms = block2natoms(block);
1630 printf("Counted atom numbers up to %d in index file\n", natoms);
1633 edit_index(natoms, atoms, x, block, &gnames, bVerbose);
1635 write_index(ndxoutfile, block, gnames, bDuplicate, natoms);