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) 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.
48 #include "gromacs/topology/atoms.h"
49 #include "gromacs/topology/block.h"
50 #include "gromacs/topology/invblock.h"
51 #include "gromacs/topology/residuetypes.h"
52 #include "gromacs/utility/arraysize.h"
53 #include "gromacs/utility/cstringutil.h"
54 #include "gromacs/utility/fatalerror.h"
55 #include "gromacs/utility/futil.h"
56 #include "gromacs/utility/smalloc.h"
57 #include "gromacs/utility/strdb.h"
59 static gmx_bool gmx_ask_yesno(gmx_bool bASK)
67 c = toupper(fgetc(stdin));
69 while ((c != 'Y') && (c != 'N'));
79 void write_index(const char *outf, t_blocka *b, char **gnames, gmx_bool bDuplicate, int natoms)
84 out = gmx_ffopen(outf, "w");
85 /* fprintf(out,"%5d %5d\n",b->nr,b->nra); */
86 for (i = 0; (i < b->nr); i++)
88 fprintf(out, "[ %s ]", gnames[i]);
89 for (k = 0, j = b->index[i]; j < b->index[i+1]; j++, k++)
91 const char sep = (k % 15 == 0 ? '\n' : ' ');
92 fprintf(out, "%c%4d", sep, b->a[j]+1);
97 /* Duplicate copy, useful for computational electrophysiology double-layer setups */
100 fprintf(stderr, "Duplicating the whole system with an atom offset of %d atoms.\n", natoms);
101 for (i = 0; (i < b->nr); i++)
103 fprintf(out, "[ %s_copy ]", gnames[i]);
104 for (k = 0, j = b->index[i]; j < b->index[i+1]; j++, k++)
106 const char sep = (k % 15 == 0 ? '\n' : ' ');
107 fprintf(out, "%c%4d", sep, b->a[j]+1 + natoms);
116 void add_grp(t_blocka *b, char ***gnames, gmx::ArrayRef<const int> a, const std::string &name)
118 srenew(b->index, b->nr+2);
119 srenew(*gnames, b->nr+1);
120 (*gnames)[b->nr] = gmx_strdup(name.c_str());
122 srenew(b->a, b->nra+a.size());
123 for (int i = 0; (i < a.size()); i++)
125 b->a[b->nra++] = a[i];
128 b->index[b->nr] = b->nra;
132 * Compare index in groups.
134 * Checks the index in \p a to the one in \p b at \p index.
135 * If \p index is < 0, it is taken as relative to the end of \p b.
137 * \param[in] b Block with groups.
138 * \param[in] a New group to compare to.
139 * \param[in] index The index to check.
140 * \returns True if groups are the same.
142 static bool grp_cmp(t_blocka *b, gmx::ArrayRef<const int> a, int index)
146 index = b->nr-1+index;
150 gmx_fatal(FARGS, "no such index group %d in t_blocka (nr=%d)", index, b->nr);
153 if (a.size() != b->index[index+1] - b->index[index])
157 for (int i = 0; i < a.size(); i++)
159 if (a[i] != b->a[b->index[index]+i])
166 //! Print out how many residues of a certain type are present.
167 static void p_status(gmx::ArrayRef<const std::string> restype,
168 gmx::ArrayRef<const std::string> typenames)
170 std::vector<int> counter(typenames.size(), 0);
171 std::transform(typenames.begin(), typenames.end(), counter.begin(),
172 [&restype](const std::string &typenm) {
173 return std::count_if(restype.begin(), restype.end(),
174 [&typenm](const std::string &res) {
175 return gmx_strcasecmp(res.c_str(), typenm.c_str()) == 0;
181 for (int i = 0; (i < typenames.size()); i++)
185 printf("There are: %5d %10s residues\n", counter[i], typenames[i].c_str());
191 * Return a new group of atoms with \p restype that match \p typestring,
192 * or not matching if \p bMatch.
194 * \param[in] atoms Atoms data to use.
195 * \param[in] restype Residuetypes to match.
196 * \param[in] typestring Which type the residue should have.
197 * \param[in] bMatch whether to return matching atoms or those that don't.
198 * \returns Vector of atoms that match.
200 static std::vector<int> mk_aid(const t_atoms *atoms,
201 gmx::ArrayRef<const std::string> restype,
202 const std::string &typestring,
206 for (int i = 0; (i < atoms->nr); i++)
208 bool res = gmx_strcasecmp(restype[atoms->atom[i].resind].c_str(), typestring.c_str()) == 0;
228 static void analyse_other(gmx::ArrayRef<std::string> restype, const t_atoms *atoms,
229 t_blocka *gb, char ***gn, gmx_bool bASK, gmx_bool bVerb)
231 restp_t *restp = nullptr;
232 char **attp = nullptr;
234 int i, resind, natp, nrestp = 0;
236 for (i = 0; (i < atoms->nres); i++)
238 if (gmx_strcasecmp(restype[i].c_str(), "Protein") && gmx_strcasecmp(restype[i].c_str(), "DNA") && gmx_strcasecmp(restype[i].c_str(), "RNA") && gmx_strcasecmp(restype[i].c_str(), "Water"))
248 printf("Analysing residues not classified as Protein/DNA/RNA/Water and splitting into groups...\n");
250 for (int k = 0; (k < atoms->nr); k++)
252 resind = atoms->atom[k].resind;
253 rname = *atoms->resinfo[resind].name;
254 if (gmx_strcasecmp(restype[resind].c_str(), "Protein") && gmx_strcasecmp(restype[resind].c_str(), "DNA") &&
255 gmx_strcasecmp(restype[resind].c_str(), "RNA") && gmx_strcasecmp(restype[resind].c_str(), "Water"))
258 for (l = 0; (l < nrestp); l++)
261 if (strcmp(restp[l].rname, rname) == 0)
268 srenew(restp, nrestp+1);
269 restp[nrestp].rname = gmx_strdup(rname);
270 restp[nrestp].bNeg = FALSE;
271 restp[nrestp].gname = gmx_strdup(rname);
276 for (int i = 0; (i < nrestp); i++)
278 std::vector<int> aid;
279 for (int j = 0; (j < atoms->nr); j++)
281 rname = *atoms->resinfo[atoms->atom[j].resind].name;
282 if ((strcmp(restp[i].rname, rname) == 0 && !restp[i].bNeg) ||
283 (strcmp(restp[i].rname, rname) != 0 && restp[i].bNeg))
288 add_grp(gb, gn, aid, restp[i].gname);
291 printf("split %s into atoms (y/n) ? ", restp[i].gname);
293 if (gmx_ask_yesno(bASK))
296 for (size_t k = 0; (k < aid.size()); k++)
298 aname = *atoms->atomname[aid[k]];
300 for (l = 0; (l < natp); l++)
302 if (strcmp(aname, attp[l]) == 0)
309 srenew(attp, ++natp);
310 attp[natp-1] = aname;
315 for (int l = 0; (l < natp); l++)
317 std::vector<int> aaid;
318 for (size_t k = 0; (k < aid.size()); k++)
320 aname = *atoms->atomname[aid[k]];
321 if (strcmp(aname, attp[l]) == 0)
323 aaid.push_back(aid[k]);
326 add_grp(gb, gn, aaid, attp[l]);
333 sfree(restp[i].rname);
334 sfree(restp[i].gname);
341 * Data necessary to construct a single (protein) index group in
344 typedef struct gmx_help_make_index_group // NOLINT(clang-analyzer-optin.performance.Padding)
346 /** The set of atom names that will be used to form this index group */
347 const char **defining_atomnames;
348 /** Size of the defining_atomnames array */
349 int num_defining_atomnames;
350 /** Name of this index group */
351 const char *group_name;
352 /** Whether the above atom names name the atoms in the group, or
353 those not in the group */
354 gmx_bool bTakeComplement;
355 /** The index in wholename gives the first item in the arrays of
356 atomnames that should be tested with 'gmx_strncasecmp' in stead of
357 gmx_strcasecmp, or -1 if all items should be tested with strcasecmp
358 This is comparable to using a '*' wildcard at the end of specific
359 atom names, but that is more involved to implement...
362 /** Only create this index group if it differs from the one specified in compareto,
363 where -1 means to always create this group. */
365 } t_gmx_help_make_index_group;
367 static void analyse_prot(gmx::ArrayRef<const std::string> restype, const t_atoms *atoms,
368 t_blocka *gb, char ***gn, gmx_bool bASK, gmx_bool bVerb)
370 /* lists of atomnames to be used in constructing index groups: */
371 static const char *pnoh[] = { "H", "HN" };
372 static const char *pnodum[] = {
373 "MN1", "MN2", "MCB1", "MCB2", "MCG1", "MCG2",
374 "MCD1", "MCD2", "MCE1", "MCE2", "MNZ1", "MNZ2"
376 static const char *calpha[] = { "CA" };
377 static const char *bb[] = { "N", "CA", "C" };
378 static const char *mc[] = { "N", "CA", "C", "O", "O1", "O2", "OC1", "OC2", "OT", "OXT" };
379 static const char *mcb[] = { "N", "CA", "CB", "C", "O", "O1", "O2", "OC1", "OC2", "OT", "OXT" };
380 static const char *mch[] = {
381 "N", "CA", "C", "O", "O1", "O2", "OC1", "OC2", "OT", "OXT",
382 "H1", "H2", "H3", "H", "HN"
385 static const t_gmx_help_make_index_group constructing_data[] =
386 {{ nullptr, 0, "Protein", TRUE, -1, -1},
387 { pnoh, asize(pnoh), "Protein-H", TRUE, 0, -1},
388 { calpha, asize(calpha), "C-alpha", FALSE, -1, -1},
389 { bb, asize(bb), "Backbone", FALSE, -1, -1},
390 { mc, asize(mc), "MainChain", FALSE, -1, -1},
391 { mcb, asize(mcb), "MainChain+Cb", FALSE, -1, -1},
392 { mch, asize(mch), "MainChain+H", FALSE, -1, -1},
393 { mch, asize(mch), "SideChain", TRUE, -1, -1},
394 { mch, asize(mch), "SideChain-H", TRUE, 11, -1},
395 { pnodum, asize(pnodum), "Prot-Masses", TRUE, -1, 0}, };
396 const int num_index_groups = asize(constructing_data);
401 char ndx_name[STRLEN], *atnm;
406 printf("Analysing Protein...\n");
408 std::vector<int> aid;
410 /* calculate the number of protein residues */
412 for (i = 0; (i < atoms->nres); i++)
414 if (0 == gmx_strcasecmp(restype[i].c_str(), "Protein"))
419 /* find matching or complement atoms */
420 for (i = 0; (i < num_index_groups); i++)
422 for (n = 0; (n < atoms->nr); n++)
424 if (0 == gmx_strcasecmp(restype[atoms->atom[n].resind].c_str(), "Protein"))
427 for (j = 0; (j < constructing_data[i].num_defining_atomnames); j++)
429 /* skip digits at beginning of atomname, e.g. 1H */
430 atnm = *atoms->atomname[n];
431 while (isdigit(atnm[0]))
435 if ( (constructing_data[i].wholename == -1) || (j < constructing_data[i].wholename) )
437 if (0 == gmx_strcasecmp(constructing_data[i].defining_atomnames[j], atnm))
444 if (0 == gmx_strncasecmp(constructing_data[i].defining_atomnames[j], atnm, strlen(constructing_data[i].defining_atomnames[j])))
450 if (constructing_data[i].bTakeComplement != match)
456 /* if we want to add this group always or it differs from previous
458 if (-1 == constructing_data[i].compareto || !grp_cmp(gb, aid, constructing_data[i].compareto-i) )
460 add_grp(gb, gn, aid, constructing_data[i].group_name);
467 for (i = 0; (i < num_index_groups); i++)
469 printf("Split %12s into %5d residues (y/n) ? ", constructing_data[i].group_name, npres);
470 if (gmx_ask_yesno(bASK))
474 for (n = 0; ((atoms->atom[n].resind < npres) && (n < atoms->nr)); )
476 resind = atoms->atom[n].resind;
477 for (; ((atoms->atom[n].resind == resind) && (n < atoms->nr)); n++)
480 for (j = 0; (j < constructing_data[i].num_defining_atomnames); j++)
482 if (0 == gmx_strcasecmp(constructing_data[i].defining_atomnames[j], *atoms->atomname[n]))
487 if (constructing_data[i].bTakeComplement != match)
492 /* copy the residuename to the tail of the groupname */
496 ri = &atoms->resinfo[resind];
497 sprintf(ndx_name, "%s_%s%d%c",
498 constructing_data[i].group_name, *ri->name, ri->nr, ri->ic == ' ' ? '\0' : ri->ic);
499 add_grp(gb, gn, aid, ndx_name);
505 printf("Make group with sidechain and C=O swapped (y/n) ? ");
506 if (gmx_ask_yesno(bASK))
508 /* Make swap sidechain C=O index */
510 for (int n = 0; ((atoms->atom[n].resind < npres) && (n < atoms->nr)); )
512 int resind = atoms->atom[n].resind;
514 for (; ((atoms->atom[n].resind == resind) && (n < atoms->nr)); n++)
516 if (strcmp("CA", *atoms->atomname[n]) == 0)
520 aid.resize(aid.size() + 3);
522 else if (strcmp("C", *atoms->atomname[n]) == 0)
526 gmx_incons("Atom naming problem");
530 else if (strcmp("O", *atoms->atomname[n]) == 0)
534 gmx_incons("Atom naming problem");
538 else if (strcmp("O1", *atoms->atomname[n]) == 0)
542 gmx_incons("Atom naming problem");
552 /* copy the residuename to the tail of the groupname */
555 add_grp(gb, gn, aid, "SwapSC-CO");
562 void analyse(const t_atoms *atoms, t_blocka *gb, char ***gn, gmx_bool bASK, gmx_bool bVerb)
571 printf("Analysing residue names:\n");
573 /* Create system group, every single atom */
574 std::vector<int> aid(atoms->nr);
575 for (i = 0; i < atoms->nr; i++)
579 add_grp(gb, gn, aid, "System");
581 /* For every residue, get a pointer to the residue type name */
584 std::vector<std::string> restype;
585 std::vector<std::string> previousTypename;
590 resnm = *atoms->resinfo[i].name;
591 restype.emplace_back(rt.typeNameForIndexedResidue(resnm));
592 previousTypename.push_back(restype[i]);
594 for (i = 1; i < atoms->nres; i++)
596 resnm = *atoms->resinfo[i].name;
597 restype.emplace_back(rt.typeNameForIndexedResidue(resnm));
599 /* Note that this does not lead to a N*N loop, but N*K, where
600 * K is the number of residue _types_, which is small and independent of N.
603 for (size_t k = 0; k < previousTypename.size() && !found; k++)
605 found = strcmp(restype[i].c_str(), previousTypename[k].c_str()) == 0;
609 previousTypename.push_back(restype[i]);
616 p_status(restype, previousTypename);
619 for (int k = 0; k < gmx::index(previousTypename.size()); k++)
621 aid = mk_aid(atoms, restype, previousTypename[k], TRUE);
623 /* Check for special types to do fancy stuff with */
625 if (!gmx_strcasecmp(previousTypename[k].c_str(), "Protein") && !aid.empty())
628 analyse_prot(restype, atoms, gb, gn, bASK, bVerb);
630 /* Create a Non-Protein group */
631 aid = mk_aid(atoms, restype, "Protein", FALSE);
632 if ((!aid.empty()) && (static_cast<int>(aid.size()) < atoms->nr))
634 add_grp(gb, gn, aid, "non-Protein");
637 else if (!gmx_strcasecmp(previousTypename[k].c_str(), "Water") && !aid.empty())
639 add_grp(gb, gn, aid, previousTypename[k]);
640 /* Add this group as 'SOL' too, for backward compatibility with older gromacs versions */
641 add_grp(gb, gn, aid, "SOL");
644 /* Solvent, create a negated group too */
645 aid = mk_aid(atoms, restype, "Water", FALSE);
646 if ((!aid.empty()) && (static_cast<int>(aid.size()) < atoms->nr))
648 add_grp(gb, gn, aid, "non-Water");
651 else if (!aid.empty())
654 add_grp(gb, gn, aid, previousTypename[k]);
655 analyse_other(restype, atoms, gb, gn, bASK, bVerb);
660 /* Create a merged water_and_ions group */
666 for (i = 0; i < gb->nr; i++)
668 if (!gmx_strcasecmp((*gn)[i], "Water"))
671 nwater = gb->index[i+1]-gb->index[i];
673 else if (!gmx_strcasecmp((*gn)[i], "Ion"))
676 nion = gb->index[i+1]-gb->index[i];
680 if (nwater > 0 && nion > 0)
682 srenew(gb->index, gb->nr+2);
683 srenew(*gn, gb->nr+1);
684 (*gn)[gb->nr] = gmx_strdup("Water_and_ions");
685 srenew(gb->a, gb->nra+nwater+nion);
688 for (i = gb->index[iwater]; i < gb->index[iwater+1]; i++)
690 gb->a[gb->nra++] = gb->a[i];
695 for (i = gb->index[iion]; i < gb->index[iion+1]; i++)
697 gb->a[gb->nra++] = gb->a[i];
701 gb->index[gb->nr] = gb->nra;
706 void check_index(const char *gname, int n, int index[], const char *traj, int natoms)
710 for (i = 0; i < n; i++)
712 if (index[i] >= natoms)
714 gmx_fatal(FARGS, "%s atom number (index[%d]=%d) is larger than the number of atoms in %s (%d)",
715 gname ? gname : "Index", i+1, index[i]+1,
716 traj ? traj : "the trajectory", natoms);
718 else if (index[i] < 0)
720 gmx_fatal(FARGS, "%s atom number (index[%d]=%d) is less than zero",
721 gname ? gname : "Index", i+1, index[i]+1);
726 t_blocka *init_index(const char *gfile, char ***grpname)
732 char line[STRLEN], *pt, str[STRLEN];
734 in = gmx_ffopen(gfile, "r");
742 while (get_a_line(in, line, STRLEN))
744 if (get_header(line, str))
747 srenew(b->index, b->nr+1);
748 srenew(*grpname, b->nr);
753 b->index[b->nr] = b->index[b->nr-1];
754 (*grpname)[b->nr-1] = gmx_strdup(str);
760 gmx_fatal(FARGS, "The first header of your indexfile is invalid");
763 while (sscanf(pt, "%s", str) == 1)
769 srenew(b->a, maxentries);
771 assert(b->a != nullptr); // for clang analyzer
772 b->a[i] = strtol(str, nullptr, 10)-1;
775 pt = strstr(pt, str)+strlen(str);
781 for (i = 0; (i < b->nr); i++)
783 assert(b->a != nullptr); // for clang analyzer
784 for (j = b->index[i]; (j < b->index[i+1]); j++)
788 fprintf(stderr, "\nWARNING: negative index %d in group %s\n\n",
789 b->a[j], (*grpname)[i]);
797 static void minstring(char *str)
801 for (i = 0; (i < static_cast<int>(strlen(str))); i++)
810 int find_group(const char *s, int ngrps, char **grpname)
818 /* first look for whole name match */
820 for (i = 0; i < ngrps; i++)
822 if (gmx_strcasecmp_min(s, grpname[i]) == 0)
832 /* second look for first string match */
835 for (i = 0; i < ngrps; i++)
837 if (gmx_strncasecmp_min(s, grpname[i], n) == 0)
847 /* last look for arbitrary substring match */
851 strncpy(key, s, sizeof(key)-1);
852 key[STRLEN-1] = '\0';
855 for (i = 0; i < ngrps; i++)
857 strcpy(string, grpname[i]);
860 if (strstr(string, key) != nullptr)
872 printf("Error: Multiple groups '%s' selected\n", s);
878 static int qgroup(int *a, int ngrps, char **grpname)
887 fprintf(stderr, "Select a group: ");
890 if (scanf("%s", s) != 1)
892 gmx_fatal(FARGS, "Cannot read from input");
894 trim(s); /* remove spaces */
896 while (strlen(s) == 0);
897 aa = strtol(s, &end, 10);
898 if (aa == 0 && end[0] != '\0') /* string entered */
900 aa = find_group(s, ngrps, grpname);
902 bInRange = (aa >= 0 && aa < ngrps);
905 printf("Error: No such group '%s'\n", s);
909 printf("Selected %d: '%s'\n", aa, grpname[aa]);
914 static void rd_groups(t_blocka *grps, char **grpname, char *gnames[],
915 int ngrps, int isize[], int *index[], int grpnr[])
921 gmx_fatal(FARGS, "Error: no groups in indexfile");
923 for (i = 0; (i < grps->nr); i++)
925 fprintf(stderr, "Group %5d (%15s) has %5d elements\n", i, grpname[i],
926 grps->index[i+1]-grps->index[i]);
928 for (i = 0; (i < ngrps); i++)
934 gnr1 = qgroup(&grpnr[i], grps->nr, grpname);
935 if ((gnr1 < 0) || (gnr1 >= grps->nr))
937 fprintf(stderr, "Select between %d and %d.\n", 0, grps->nr-1);
940 while ((gnr1 < 0) || (gnr1 >= grps->nr));
944 fprintf(stderr, "There is one group in the index\n");
947 gnames[i] = gmx_strdup(grpname[gnr1]);
948 isize[i] = grps->index[gnr1+1]-grps->index[gnr1];
949 snew(index[i], isize[i]);
950 for (j = 0; (j < isize[i]); j++)
952 index[i][j] = grps->a[grps->index[gnr1]+j];
957 void rd_index(const char *statfile, int ngrps, int isize[],
958 int *index[], char *grpnames[])
967 gmx_fatal(FARGS, "No index file specified");
969 grps = init_index(statfile, &gnames);
970 rd_groups(grps, gnames, grpnames, ngrps, isize, index, grpnr);
971 for (int i = 0; i < grps->nr; i++)
981 void get_index(const t_atoms *atoms, const char *fnm, int ngrps,
982 int isize[], int *index[], char *grpnames[])
985 t_blocka *grps = nullptr;
992 grps = init_index(fnm, gnames);
997 snew(grps->index, 1);
998 analyse(atoms, grps, gnames, FALSE, FALSE);
1002 gmx_incons("You need to supply a valid atoms structure or a valid index file name");
1005 rd_groups(grps, *gnames, grpnames, ngrps, isize, index, grpnr);
1006 for (int i = 0; i < grps->nr; ++i)
1008 sfree((*gnames)[i]);
1017 t_cluster_ndx *cluster_index(FILE *fplog, const char *ndx)
1023 c->clust = init_index(ndx, &c->grpname);
1025 for (i = 0; (i < c->clust->nra); i++)
1027 c->maxframe = std::max(c->maxframe, c->clust->a[i]);
1029 fprintf(fplog ? fplog : stdout,
1030 "There are %d clusters containing %d structures, highest framenr is %d\n",
1031 c->clust->nr, c->clust->nra, c->maxframe);
1034 pr_blocka(debug, 0, "clust", c->clust, TRUE);
1035 for (i = 0; (i < c->clust->nra); i++)
1037 if ((c->clust->a[i] < 0) || (c->clust->a[i] > c->maxframe))
1039 gmx_fatal(FARGS, "Range check error for c->clust->a[%d] = %d\n"
1040 "should be within 0 and %d", i, c->clust->a[i], c->maxframe+1);
1044 c->inv_clust = make_invblocka(c->clust, c->maxframe);