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.
43 #include "gromacs/fileio/futil.h"
46 #include "gromacs/commandline/pargs.h"
47 #include "gromacs/fileio/confio.h"
57 gmx_bool bCase = FALSE;
59 static int or_groups(atom_id nr1, atom_id *at1, atom_id nr2, atom_id *at2,
60 atom_id *nr, atom_id *at)
62 atom_id i1, i2, max = 0;
68 for (i1 = 0; i1 < nr1; i1++)
70 if ((i1 > 0) && (at1[i1] <= max))
76 for (i1 = 0; i1 < nr2; i1++)
78 if ((i1 > 0) && (at2[i1] <= max))
87 printf("One of your groups is not ascending\n");
94 while ((i1 < nr1) || (i2 < nr2))
96 if ((i2 == nr2) || ((i1 < nr1) && (at1[i1] < at2[i2])))
104 if ((i2 < nr2) && ((i1 == nr1) || (at1[i1] > at2[i2])))
113 printf("Merged two groups with OR: %d %d -> %d\n", nr1, nr2, *nr);
119 static int and_groups(atom_id nr1, atom_id *at1, atom_id nr2, atom_id *at2,
120 atom_id *nr, atom_id *at)
125 for (i1 = 0; i1 < nr1; i1++)
127 for (i2 = 0; i2 < nr2; i2++)
129 if (at1[i1] == at2[i2])
137 printf("Merged two groups with AND: %d %d -> %d\n", nr1, nr2, *nr);
142 static gmx_bool is_name_char(char c)
144 /* This string should contain all characters that can not be
145 * the first letter of a name due to the make_ndx syntax.
147 const char *spec = " !&|";
149 return (c != '\0' && strchr(spec, c) == NULL);
152 static int parse_names(char **string, int *n_names, char **names)
157 while (is_name_char((*string)[0]) || (*string)[0] == ' ')
159 if (is_name_char((*string)[0]))
161 if (*n_names >= MAXNAMES)
163 gmx_fatal(FARGS, "To many names: %d\n", *n_names+1);
166 while (is_name_char((*string)[i]))
168 names[*n_names][i] = (*string)[i];
172 printf("Name is too long, the maximum is %d characters\n", NAME_LEN);
176 names[*n_names][i] = '\0';
179 upstring(names[*n_names]);
193 static gmx_bool parse_int_char(char **string, int *nr, char *c)
200 while ((*string)[0] == ' ')
209 if (isdigit((*string)[0]))
211 *nr = (*string)[0]-'0';
213 while (isdigit((*string)[0]))
215 *nr = (*nr)*10+(*string)[0]-'0';
218 if (isalpha((*string)[0]))
223 /* Check if there is at most one non-digit character */
224 if (!isalnum((*string)[0]))
241 static gmx_bool parse_int(char **string, int *nr)
247 bRet = parse_int_char(string, nr, &c);
248 if (bRet && c != ' ')
257 static gmx_bool isquote(char c)
262 static gmx_bool parse_string(char **string, int *nr, int ngrps, char **grpname)
267 while ((*string)[0] == ' ')
273 if (isquote((*string)[0]))
277 s = strdup((*string));
281 (*string) += sp-s + 1;
283 (*nr) = find_group(s, ngrps, grpname);
287 return (*nr) != NOTSET;
290 static int select_atomnumbers(char **string, t_atoms *atoms, atom_id n1,
291 atom_id *nr, atom_id *index, char *gname)
297 while ((*string)[0] == ' ')
301 if ((*string)[0] == '-')
304 parse_int(string, &up);
305 if ((n1 < 1) || (n1 > atoms->nr) || (up < 1) || (up > atoms->nr))
307 printf("Invalid atom range\n");
311 for (i = n1-1; i <= up-1; i++)
316 printf("Found %d atom%s in range %d-%d\n", *nr, (*nr == 1) ? "" : "s", n1, up);
319 sprintf(buf, "a_%d", n1);
323 sprintf(buf, "a_%d-%d", n1, up);
334 if ((i-1 >= 0) && (i-1 < atoms->nr))
338 sprintf(buf, "_%d", i);
343 printf("Invalid atom number %d\n", i);
347 while ((*nr != 0) && (parse_int(string, &i)));
353 static int select_residuenumbers(char **string, t_atoms *atoms,
355 atom_id *nr, atom_id *index, char *gname)
362 while ((*string)[0] == ' ')
366 if ((*string)[0] == '-')
368 /* Residue number range selection */
371 printf("Error: residue insertion codes can not be used with residue range selection\n");
375 parse_int(string, &up);
377 for (i = 0; i < atoms->nr; i++)
379 ri = &atoms->resinfo[atoms->atom[i].resind];
380 for (j = n1; (j <= up); j++)
382 if (ri->nr == j && (c == ' ' || ri->ic == c))
389 printf("Found %d atom%s with res.nr. in range %d-%d\n",
390 *nr, (*nr == 1) ? "" : "s", n1, up);
393 sprintf(buf, "r_%d", n1);
397 sprintf(buf, "r_%d-%d", n1, up);
403 /* Individual residue number/insertion code selection */
408 for (i = 0; i < atoms->nr; i++)
410 ri = &atoms->resinfo[atoms->atom[i].resind];
411 if (ri->nr == j && ri->ic == c)
417 sprintf(buf, "_%d", j);
420 while (parse_int_char(string, &j, &c));
426 static int select_residueindices(char **string, t_atoms *atoms,
428 atom_id *nr, atom_id *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",
465 *nr, (*nr == 1) ? "" : "s", n1, up);
468 sprintf(buf, "r_%d", n1);
472 sprintf(buf, "r_%d-%d", n1, up);
478 /* Individual residue number/insertion code selection */
483 for (i = 0; i < atoms->nr; i++)
485 ri = &atoms->resinfo[atoms->atom[i].resind];
486 if (atoms->atom[i].resind+1 == j && ri->ic == c)
492 sprintf(buf, "_%d", j);
495 while (parse_int_char(string, &j, &c));
502 static gmx_bool atoms_from_residuenumbers(t_atoms *atoms, int group, t_blocka *block,
503 atom_id *nr, atom_id *index, char *gname)
505 int i, j, j0, j1, resnr, nres;
507 j0 = block->index[group];
508 j1 = block->index[group+1];
510 for (j = j0; j < j1; j++)
512 if (block->a[j] >= nres)
514 printf("Index %s contains number>nres (%d>%d)\n",
515 gname, block->a[j]+1, nres);
519 for (i = 0; i < atoms->nr; i++)
521 resnr = atoms->resinfo[atoms->atom[i].resind].nr;
522 for (j = j0; j < j1; j++)
524 if (block->a[j]+1 == resnr)
532 printf("Found %d atom%s in %d residues from group %s\n",
533 *nr, (*nr == 1) ? "" : "s", j1-j0, gname);
537 static gmx_bool comp_name(char *name, char *search)
539 while (name[0] != '\0' && search[0] != '\0')
547 if (search[1] != '\0')
549 printf("WARNING: Currently '*' is only supported at the end of an expression\n");
558 /* Compare a single character */
559 if (( bCase && strncmp(name, search, 1)) ||
560 (!bCase && gmx_strncasecmp(name, search, 1)))
569 return (name[0] == '\0' && (search[0] == '\0' || search[0] == '*'));
572 static int select_chainnames(t_atoms *atoms, int n_names, char **names,
573 atom_id *nr, atom_id *index)
581 for (i = 0; i < atoms->nr; i++)
583 name[0] = atoms->resinfo[atoms->atom[i].resind].chainid;
585 while (j < n_names && !comp_name(name, names[j]))
595 printf("Found %d atom%s with chain identifier%s",
596 *nr, (*nr == 1) ? "" : "s", (n_names == 1) ? "" : "s");
597 for (j = 0; (j < n_names); j++)
599 printf(" %s", names[j]);
606 static int select_atomnames(t_atoms *atoms, int n_names, char **names,
607 atom_id *nr, atom_id *index, gmx_bool bType)
614 for (i = 0; i < atoms->nr; i++)
618 name = *(atoms->atomtype[i]);
622 name = *(atoms->atomname[i]);
625 while (j < n_names && !comp_name(name, names[j]))
635 printf("Found %d atoms with %s%s",
636 *nr, bType ? "type" : "name", (n_names == 1) ? "" : "s");
637 for (j = 0; (j < n_names); j++)
639 printf(" %s", names[j]);
646 static int select_residuenames(t_atoms *atoms, int n_names, char **names,
647 atom_id *nr, atom_id *index)
654 for (i = 0; i < atoms->nr; i++)
656 name = *(atoms->resinfo[atoms->atom[i].resind].name);
658 while (j < n_names && !comp_name(name, names[j]))
668 printf("Found %d atoms with residue name%s", *nr, (n_names == 1) ? "" : "s");
669 for (j = 0; (j < n_names); j++)
671 printf(" %s", names[j]);
678 static void copy2block(int n, atom_id *index, t_blocka *block)
685 srenew(block->index, block->nr+1);
686 block->index[block->nr] = n0+n;
687 srenew(block->a, n0+n);
688 for (i = 0; (i < n); i++)
690 block->a[n0+i] = index[i];
694 static void make_gname(int n, char **names, char *gname)
698 strcpy(gname, names[0]);
699 for (i = 1; i < n; i++)
702 strcat(gname, names[i]);
706 static void copy_group(int g, t_blocka *block, atom_id *nr, atom_id *index)
710 i0 = block->index[g];
711 *nr = block->index[g+1]-i0;
712 for (i = 0; i < *nr; i++)
714 index[i] = block->a[i0+i];
718 static void remove_group(int nr, int nr2, t_blocka *block, char ***gn)
728 for (j = 0; j <= nr2-nr; j++)
730 if ((nr < 0) || (nr >= block->nr))
732 printf("Group %d does not exist\n", nr+j);
736 shift = block->index[nr+1]-block->index[nr];
737 for (i = block->index[nr+1]; i < block->nra; i++)
739 block->a[i-shift] = block->a[i];
742 for (i = nr; i < block->nr; i++)
744 block->index[i] = block->index[i+1]-shift;
746 name = strdup((*gn)[nr]);
748 for (i = nr; i < block->nr-1; i++)
750 (*gn)[i] = (*gn)[i+1];
753 block->nra = block->index[block->nr];
754 printf("Removed group %d '%s'\n", nr+j, name);
760 static void split_group(t_atoms *atoms, int sel_nr, t_blocka *block, char ***gn,
763 char buf[STRLEN], *name;
767 printf("Splitting group %d '%s' into %s\n", sel_nr, (*gn)[sel_nr],
768 bAtom ? "atoms" : "residues");
770 n0 = block->index[sel_nr];
771 n1 = block->index[sel_nr+1];
772 srenew(block->a, block->nra+n1-n0);
773 for (i = n0; i < n1; i++)
776 resind = atoms->atom[a].resind;
777 name = *(atoms->resinfo[resind].name);
778 if (bAtom || (i == n0) || (atoms->atom[block->a[i-1]].resind != resind))
782 block->index[block->nr] = block->nra;
785 srenew(block->index, block->nr+1);
786 srenew(*gn, block->nr);
789 sprintf(buf, "%s_%s_%d", (*gn)[sel_nr], *atoms->atomname[a], a+1);
793 sprintf(buf, "%s_%s_%d", (*gn)[sel_nr], name, atoms->resinfo[resind].nr);
795 (*gn)[block->nr-1] = strdup(buf);
797 block->a[block->nra] = a;
800 block->index[block->nr] = block->nra;
803 static int split_chain(t_atoms *atoms, rvec *x,
804 int sel_nr, t_blocka *block, char ***gn)
808 atom_id i, a, natoms, *start = NULL, *end = NULL, ca_start, ca_end;
815 while (ca_start < natoms)
817 while ((ca_start < natoms) && strcmp(*atoms->atomname[ca_start], "CA"))
821 if (ca_start < natoms)
823 srenew(start, nchain+1);
824 srenew(end, nchain+1);
825 start[nchain] = ca_start;
826 while ((start[nchain] > 0) &&
827 (atoms->atom[start[nchain]-1].resind ==
828 atoms->atom[ca_start].resind))
841 while ((i < natoms) && strcmp(*atoms->atomname[i], "CA"));
844 rvec_sub(x[ca_end], x[i], vec);
847 while ((i < natoms) && (norm(vec) < 0.45));
849 end[nchain] = ca_end;
850 while ((end[nchain]+1 < natoms) &&
851 (atoms->atom[end[nchain]+1].resind == atoms->atom[ca_end].resind))
855 ca_start = end[nchain]+1;
861 printf("Found 1 chain, will not split\n");
865 printf("Found %d chains\n", nchain);
867 for (j = 0; j < nchain; j++)
869 printf("%d:%6d atoms (%d to %d)\n",
870 j+1, end[j]-start[j]+1, start[j]+1, end[j]+1);
875 srenew(block->a, block->nra+block->index[sel_nr+1]-block->index[sel_nr]);
876 for (j = 0; j < nchain; j++)
879 srenew(block->index, block->nr+1);
880 srenew(*gn, block->nr);
881 sprintf(buf, "%s_chain%d", (*gn)[sel_nr], j+1);
882 (*gn)[block->nr-1] = strdup(buf);
883 for (i = block->index[sel_nr]; i < block->index[sel_nr+1]; i++)
886 if ((a >= start[j]) && (a <= end[j]))
888 block->a[block->nra] = a;
892 block->index[block->nr] = block->nra;
893 if (block->index[block->nr-1] == block->index[block->nr])
895 remove_group(block->nr-1, NOTSET, block, gn);
905 static gmx_bool check_have_atoms(t_atoms *atoms, char *string)
909 printf("Can not process '%s' without atom info, use option -f\n", string);
918 static gmx_bool parse_entry(char **string, int natoms, t_atoms *atoms,
919 t_blocka *block, char ***gn,
920 atom_id *nr, atom_id *index, char *gname)
922 static char **names, *ostring;
923 static gmx_bool bFirst = TRUE;
924 int j, n_names, sel_nr1;
925 atom_id i, nr1, *index1;
927 gmx_bool bRet, bCompl;
932 snew(names, MAXNAMES);
933 for (i = 0; i < MAXNAMES; i++)
935 snew(names[i], NAME_LEN+1);
942 while (*string[0] == ' ')
947 if ((*string)[0] == '!')
951 while (*string[0] == ' ')
963 if (parse_int(string, &sel_nr1) ||
964 parse_string(string, &sel_nr1, block->nr, *gn))
966 if ((sel_nr1 >= 0) && (sel_nr1 < block->nr))
968 copy_group(sel_nr1, block, nr, index);
969 strcpy(gname, (*gn)[sel_nr1]);
970 printf("Copied index group %d '%s'\n", sel_nr1, (*gn)[sel_nr1]);
975 printf("Group %d does not exist\n", sel_nr1);
978 else if ((*string)[0] == 'a')
981 if (check_have_atoms(atoms, ostring))
983 if (parse_int(string, &sel_nr1))
985 bRet = select_atomnumbers(string, atoms, sel_nr1, nr, index, gname);
987 else if (parse_names(string, &n_names, names))
989 bRet = select_atomnames(atoms, n_names, names, nr, index, FALSE);
990 make_gname(n_names, names, gname);
994 else if ((*string)[0] == 't')
997 if (check_have_atoms(atoms, ostring) &&
998 parse_names(string, &n_names, names))
1000 if (atoms->atomtype == NULL)
1002 printf("Need a run input file to select atom types\n");
1006 bRet = select_atomnames(atoms, n_names, names, nr, index, TRUE);
1007 make_gname(n_names, names, gname);
1011 else if (strncmp(*string, "res", 3) == 0)
1014 if (check_have_atoms(atoms, ostring) &&
1015 parse_int(string, &sel_nr1) &&
1016 (sel_nr1 >= 0) && (sel_nr1 < block->nr) )
1018 bRet = atoms_from_residuenumbers(atoms,
1019 sel_nr1, block, nr, index, (*gn)[sel_nr1]);
1020 sprintf(gname, "atom_%s", (*gn)[sel_nr1]);
1023 else if (strncmp(*string, "ri", 2) == 0)
1026 if (check_have_atoms(atoms, ostring) &&
1027 parse_int_char(string, &sel_nr1, &c))
1029 bRet = select_residueindices(string, atoms, sel_nr1, c, nr, index, gname);
1032 else if ((*string)[0] == 'r')
1035 if (check_have_atoms(atoms, ostring))
1037 if (parse_int_char(string, &sel_nr1, &c))
1039 bRet = select_residuenumbers(string, atoms, sel_nr1, c, nr, index, gname);
1041 else if (parse_names(string, &n_names, names))
1043 bRet = select_residuenames(atoms, n_names, names, nr, index);
1044 make_gname(n_names, names, gname);
1048 else if (strncmp(*string, "chain", 5) == 0)
1051 if (check_have_atoms(atoms, ostring) &&
1052 parse_names(string, &n_names, names))
1054 bRet = select_chainnames(atoms, n_names, names, nr, index);
1055 sprintf(gname, "ch%s", names[0]);
1056 for (i = 1; i < n_names; i++)
1058 strcat(gname, names[i]);
1064 snew(index1, natoms-*nr);
1066 for (i = 0; i < natoms; i++)
1069 while ((j < *nr) && (index[j] != i))
1075 if (nr1 >= natoms-*nr)
1077 printf("There are double atoms in your index group\n");
1085 for (i = 0; i < nr1; i++)
1087 index[i] = index1[i];
1091 for (i = strlen(gname)+1; i > 0; i--)
1093 gname[i] = gname[i-1];
1096 printf("Complemented group: %d atoms\n", *nr);
1102 static void list_residues(t_atoms *atoms)
1104 int i, j, start, end, prev_resind, resind;
1107 /* Print all the residues, assuming continuous resnr count */
1108 start = atoms->atom[0].resind;
1109 prev_resind = start;
1110 for (i = 0; i < atoms->nr; i++)
1112 resind = atoms->atom[i].resind;
1113 if ((resind != prev_resind) || (i == atoms->nr-1))
1115 if ((bDiff = strcmp(*atoms->resinfo[resind].name,
1116 *atoms->resinfo[start].name)) ||
1129 for (j = start; j <= end; j++)
1132 j+1, *(atoms->resinfo[j].name));
1137 printf(" %4d - %4d %-5s ",
1138 start+1, end+1, *(atoms->resinfo[start].name));
1143 prev_resind = resind;
1148 static void edit_index(int natoms, t_atoms *atoms, rvec *x, t_blocka *block, char ***gn, gmx_bool bVerbose)
1150 static char **atnames, *ostring;
1151 static gmx_bool bFirst = TRUE;
1152 char inp_string[STRLEN], *string;
1153 char gname[STRLEN], gname1[STRLEN], gname2[STRLEN];
1154 int i, i0, i1, sel_nr, sel_nr2, newgroup;
1155 atom_id nr, nr1, nr2, *index, *index1, *index2;
1156 gmx_bool bAnd, bOr, bPrintOnce;
1161 snew(atnames, MAXNAMES);
1162 for (i = 0; i < MAXNAMES; i++)
1164 snew(atnames[i], NAME_LEN+1);
1170 snew(index, natoms);
1171 snew(index1, natoms);
1172 snew(index2, natoms);
1179 if (bVerbose || bPrintOnce || newgroup != NOTSET)
1182 if (bVerbose || bPrintOnce || newgroup == NOTSET)
1192 for (i = i0; i < i1; i++)
1194 printf("%3d %-20s: %5d atoms\n", i, (*gn)[i],
1195 block->index[i+1]-block->index[i]);
1199 if (bVerbose || bPrintOnce)
1202 printf(" nr : group ! 'name' nr name 'splitch' nr Enter: list groups\n");
1203 printf(" 'a': atom & 'del' nr 'splitres' nr 'l': list residues\n");
1204 printf(" 't': atom type | 'keep' nr 'splitat' nr 'h': help\n");
1205 printf(" 'r': residue 'res' nr 'chain' char\n");
1206 printf(" \"name\": group 'case': case %s 'q': save and quit\n",
1207 bCase ? "insensitive" : "sensitive ");
1208 printf(" 'ri': residue index\n");
1213 if (NULL == fgets(inp_string, STRLEN, stdin))
1215 gmx_fatal(FARGS, "Error reading user input");
1217 inp_string[strlen(inp_string)-1] = 0;
1219 string = inp_string;
1220 while (string[0] == ' ')
1227 if (string[0] == 'h')
1229 printf(" nr : selects an index group by number or quoted string.\n");
1230 printf(" The string is first matched against the whole group name,\n");
1231 printf(" then against the beginning and finally against an\n");
1232 printf(" arbitrary substring. A multiple match is an error.\n");
1234 printf(" 'a' nr1 [nr2 ...] : selects atoms, atom numbering starts at 1.\n");
1235 printf(" 'a' nr1 - nr2 : selects atoms in the range from nr1 to nr2.\n");
1236 printf(" 'a' name1[*] [name2[*] ...] : selects atoms by name(s), '?' matches any char,\n");
1237 printf(" wildcard '*' allowed at the end of a name.\n");
1238 printf(" 't' type1[*] [type2[*] ...] : as 'a', but for type, run input file required.\n");
1239 printf(" 'r' nr1[ic1] [nr2[ic2] ...] : selects residues by number and insertion code.\n");
1240 printf(" 'r' nr1 - nr2 : selects residues in the range from nr1 to nr2.\n");
1241 printf(" 'r' name1[*] [name2[*] ...] : as 'a', but for residue names.\n");
1242 printf(" 'ri' nr1 - nr2 : selects residue indices, 1-indexed, (as opposed to numbers) in the range from nr1 to nr2.\n");
1243 printf(" 'chain' ch1 [ch2 ...] : selects atoms by chain identifier(s),\n");
1244 printf(" not available with a .gro file as input.\n");
1245 printf(" ! : takes the complement of a group with respect to all\n");
1246 printf(" the atoms in the input file.\n");
1247 printf(" & | : AND and OR, can be placed between any of the options\n");
1248 printf(" above, the input is processed from left to right.\n");
1249 printf(" 'name' nr name : rename group nr to name.\n");
1250 printf(" 'del' nr1 [- nr2] : deletes one group or groups in the range from nr1 to nr2.\n");
1251 printf(" 'keep' nr : deletes all groups except nr.\n");
1252 printf(" 'case' : make all name compares case (in)sensitive.\n");
1253 printf(" 'splitch' nr : split group into chains using CA distances.\n");
1254 printf(" 'splitres' nr : split group into residues.\n");
1255 printf(" 'splitat' nr : split group into atoms.\n");
1256 printf(" 'res' nr : interpret numbers in group as residue numbers\n");
1257 printf(" Enter : list the currently defined groups and commands\n");
1258 printf(" 'l' : list the residues.\n");
1259 printf(" 'h' : show this help.\n");
1260 printf(" 'q' : save and quit.\n");
1262 printf(" Examples:\n");
1263 printf(" > 2 | 4 & r 3-5\n");
1264 printf(" selects all atoms from group 2 and 4 that have residue numbers 3, 4 or 5\n");
1265 printf(" > a C* & !a C CA\n");
1266 printf(" selects all atoms starting with 'C' but not the atoms 'C' and 'CA'\n");
1267 printf(" > \"protein\" & ! \"backb\"\n");
1268 printf(" selects all atoms that are in group 'protein' and not in group 'backbone'\n");
1271 printf("\npress Enter ");
1275 else if (strncmp(string, "del", 3) == 0)
1278 if (parse_int(&string, &sel_nr))
1280 while (string[0] == ' ')
1284 if (string[0] == '-')
1287 parse_int(&string, &sel_nr2);
1293 while (string[0] == ' ')
1297 if (string[0] == '\0')
1299 remove_group(sel_nr, sel_nr2, block, gn);
1303 printf("\nSyntax error: \"%s\"\n", string);
1307 else if (strncmp(string, "keep", 4) == 0)
1310 if (parse_int(&string, &sel_nr))
1312 remove_group(sel_nr+1, block->nr-1, block, gn);
1313 remove_group(0, sel_nr-1, block, gn);
1316 else if (strncmp(string, "name", 4) == 0)
1319 if (parse_int(&string, &sel_nr))
1321 if ((sel_nr >= 0) && (sel_nr < block->nr))
1323 sscanf(string, "%s", gname);
1324 sfree((*gn)[sel_nr]);
1325 (*gn)[sel_nr] = strdup(gname);
1329 else if (strncmp(string, "case", 4) == 0)
1332 printf("Switched to case %s\n", bCase ? "sensitive" : "insensitive");
1334 else if (string[0] == 'v')
1336 bVerbose = !bVerbose;
1337 printf("Turned verbose %s\n", bVerbose ? "on" : "off");
1339 else if (string[0] == 'l')
1341 if (check_have_atoms(atoms, ostring) )
1343 list_residues(atoms);
1346 else if (strncmp(string, "splitch", 7) == 0)
1349 if (check_have_atoms(atoms, ostring) &&
1350 parse_int(&string, &sel_nr) &&
1351 (sel_nr >= 0) && (sel_nr < block->nr))
1353 split_chain(atoms, x, sel_nr, block, gn);
1356 else if (strncmp(string, "splitres", 8) == 0)
1359 if (check_have_atoms(atoms, ostring) &&
1360 parse_int(&string, &sel_nr) &&
1361 (sel_nr >= 0) && (sel_nr < block->nr))
1363 split_group(atoms, sel_nr, block, gn, FALSE);
1366 else if (strncmp(string, "splitat", 7) == 0)
1369 if (check_have_atoms(atoms, ostring) &&
1370 parse_int(&string, &sel_nr) &&
1371 (sel_nr >= 0) && (sel_nr < block->nr))
1373 split_group(atoms, sel_nr, block, gn, TRUE);
1376 else if (string[0] == '\0')
1380 else if (string[0] != 'q')
1384 if (parse_entry(&string, natoms, atoms, block, gn, &nr, index, gname))
1388 while (string[0] == ' ')
1395 if (string[0] == '&')
1399 else if (string[0] == '|')
1408 for (i = 0; i < nr; i++)
1410 index1[i] = index[i];
1412 strcpy(gname1, gname);
1413 if (parse_entry(&string, natoms, atoms, block, gn, &nr2, index2, gname2))
1417 or_groups(nr1, index1, nr2, index2, &nr, index);
1418 sprintf(gname, "%s_%s", gname1, gname2);
1422 and_groups(nr1, index1, nr2, index2, &nr, index);
1423 sprintf(gname, "%s_&_%s", gname1, gname2);
1428 while (bAnd || bOr);
1430 while (string[0] == ' ')
1436 printf("\nSyntax error: \"%s\"\n", string);
1440 copy2block(nr, index, block);
1441 srenew(*gn, block->nr);
1442 newgroup = block->nr-1;
1443 (*gn)[newgroup] = strdup(gname);
1447 printf("Group is empty\n");
1451 while (string[0] != 'q');
1458 static int block2natoms(t_blocka *block)
1463 for (i = 0; i < block->nra; i++)
1465 natoms = max(natoms, block->a[i]);
1472 void merge_blocks(t_blocka *dest, t_blocka *source)
1476 /* count groups, srenew and fill */
1479 dest->nr += source->nr;
1480 srenew(dest->index, dest->nr+1);
1481 for (i = 0; i < source->nr; i++)
1483 dest->index[i0+i] = nra0 + source->index[i];
1485 /* count atoms, srenew and fill */
1486 dest->nra += source->nra;
1487 srenew(dest->a, dest->nra);
1488 for (i = 0; i < source->nra; i++)
1490 dest->a[nra0+i] = source->a[i];
1493 /* terminate list */
1494 dest->index[dest->nr] = dest->nra;
1498 int gmx_make_ndx(int argc, char *argv[])
1500 const char *desc[] = {
1501 "Index groups are necessary for almost every GROMACS program.",
1502 "All these programs can generate default index groups. You ONLY",
1503 "have to use [THISMODULE] when you need SPECIAL index groups.",
1504 "There is a default index group for the whole system, 9 default",
1505 "index groups for proteins, and a default index group",
1506 "is generated for every other residue name.[PAR]",
1507 "When no index file is supplied, also [THISMODULE] will generate the",
1509 "With the index editor you can select on atom, residue and chain names",
1511 "When a run input file is supplied you can also select on atom type.",
1512 "You can use NOT, AND and OR, you can split groups",
1513 "into chains, residues or atoms. You can delete and rename groups.[PAR]",
1514 "The atom numbering in the editor and the index file starts at 1.[PAR]",
1515 "The [TT]-twin[tt] switch duplicates all index groups with an offset of",
1516 "[TT]-natoms[tt], which is useful for Computational Electrophysiology",
1517 "double-layer membrane setups."
1520 static int natoms = 0;
1521 static gmx_bool bVerbose = FALSE;
1522 static gmx_bool bDuplicate = FALSE;
1524 { "-natoms", FALSE, etINT, {&natoms},
1525 "set number of atoms (default: read from coordinate or index file)" },
1526 { "-twin", FALSE, etBOOL, {&bDuplicate},
1527 "Duplicate all index groups with an offset of -natoms" },
1528 { "-verbose", FALSE, etBOOL, {&bVerbose},
1529 "HIDDENVerbose output" }
1531 #define NPA asize(pa)
1536 const char *stxfile;
1538 const char *ndxoutfile;
1545 t_blocka *block, *block2;
1546 char **gnames, **gnames2;
1548 { efSTX, "-f", NULL, ffOPTRD },
1549 { efNDX, "-n", NULL, ffOPTRDMULT },
1550 { efNDX, "-o", NULL, ffWRITE }
1552 #define NFILE asize(fnm)
1554 if (!parse_common_args(&argc, argv, 0, NFILE, fnm, NPA, pa, asize(desc), desc,
1560 stxfile = ftp2fn_null(efSTX, NFILE, fnm);
1561 if (opt2bSet("-n", NFILE, fnm))
1563 nndxin = opt2fns(&ndxinfiles, "-n", NFILE, fnm);
1569 ndxoutfile = opt2fn("-o", NFILE, fnm);
1570 bNatoms = opt2parg_bSet("-natoms", NPA, pa);
1572 if (!stxfile && !nndxin)
1574 gmx_fatal(FARGS, "No input files (structure or index)");
1580 get_stx_coordnum(stxfile, &(atoms->nr));
1581 init_t_atoms(atoms, atoms->nr, TRUE);
1584 fprintf(stderr, "\nReading structure file\n");
1585 read_stx_conf(stxfile, title, atoms, x, v, &ePBC, box);
1595 /* read input file(s) */
1596 block = new_blocka();
1598 printf("Going to read %d old index file(s)\n", nndxin);
1601 for (i = 0; i < nndxin; i++)
1603 block2 = init_index(ndxinfiles[i], &gnames2);
1604 srenew(gnames, block->nr+block2->nr);
1605 for (j = 0; j < block2->nr; j++)
1607 gnames[block->nr+j] = gnames2[j];
1610 merge_blocks(block, block2);
1612 sfree(block2->index);
1613 /* done_block(block2); */
1620 analyse(atoms, block, &gnames, FALSE, TRUE);
1625 natoms = block2natoms(block);
1626 printf("Counted atom numbers up to %d in index file\n", natoms);
1629 edit_index(natoms, atoms, x, block, &gnames, bVerbose);
1631 write_index(ndxoutfile, block, gnames, bDuplicate, natoms);