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, 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/fileio/gmxfio.h"
49 #include "gromacs/topology/atoms.h"
50 #include "gromacs/topology/block.h"
51 #include "gromacs/topology/invblock.h"
52 #include "gromacs/topology/residuetypes.h"
53 #include "gromacs/utility/arraysize.h"
54 #include "gromacs/utility/cstringutil.h"
55 #include "gromacs/utility/fatalerror.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_fio_fopen(outf
, "w");
85 /* fprintf(out,"%5d %5d\n",b->nr,b->nra); */
86 for (i
= 0; (i
< b
->nr
); i
++)
88 fprintf(out
, "[ %s ]\n", gnames
[i
]);
89 for (k
= 0, j
= b
->index
[i
]; j
< b
->index
[i
+1]; j
++, k
++)
91 fprintf(out
, "%4d ", b
->a
[j
]+1);
100 /* Duplicate copy, useful for computational electrophysiology double-layer setups */
103 fprintf(stderr
, "Duplicating the whole system with an atom offset of %d atoms.\n", natoms
);
104 for (i
= 0; (i
< b
->nr
); i
++)
106 fprintf(out
, "[ %s_copy ]\n", gnames
[i
]);
107 for (k
= 0, j
= b
->index
[i
]; j
< b
->index
[i
+1]; j
++, k
++)
109 fprintf(out
, "%4d ", b
->a
[j
]+1 + natoms
);
122 void add_grp(t_blocka
*b
, char ***gnames
, int nra
, int a
[], const char *name
)
126 srenew(b
->index
, b
->nr
+2);
127 srenew(*gnames
, b
->nr
+1);
128 (*gnames
)[b
->nr
] = gmx_strdup(name
);
130 srenew(b
->a
, b
->nra
+nra
);
131 for (i
= 0; (i
< nra
); i
++)
133 b
->a
[b
->nra
++] = a
[i
];
136 b
->index
[b
->nr
] = b
->nra
;
139 /* compare index in `a' with group in `b' at `index',
140 when `index'<0 it is relative to end of `b' */
141 static gmx_bool
grp_cmp(t_blocka
*b
, int nra
, int a
[], int index
)
147 index
= b
->nr
-1+index
;
151 gmx_fatal(FARGS
, "no such index group %d in t_blocka (nr=%d)", index
, b
->nr
);
154 if (nra
!= b
->index
[index
+1] - b
->index
[index
])
158 for (i
= 0; i
< nra
; i
++)
160 if (a
[i
] != b
->a
[b
->index
[index
]+i
])
169 p_status(const char *const *restype
, int nres
,
170 const char *const *typenames
, int ntypes
)
175 snew(counter
, ntypes
);
176 for (i
= 0; i
< ntypes
; i
++)
180 for (i
= 0; i
< nres
; i
++)
182 for (j
= 0; j
< ntypes
; j
++)
184 if (!gmx_strcasecmp(restype
[i
], typenames
[j
]))
191 for (i
= 0; (i
< ntypes
); i
++)
195 printf("There are: %5d %10s residues\n", counter
[i
], typenames
[i
]);
204 mk_aid(t_atoms
*atoms
, const char ** restype
, const char * typestring
, int *nra
, gmx_bool bMatch
)
205 /* Make an array of ints for all atoms with residuetypes matching typestring, or the opposite if bMatch is false */
213 for (i
= 0; (i
< atoms
->nr
); i
++)
215 res
= !gmx_strcasecmp(restype
[atoms
->atom
[i
].resind
], typestring
);
235 static void analyse_other(const char ** restype
, t_atoms
*atoms
,
236 t_blocka
*gb
, char ***gn
, gmx_bool bASK
, gmx_bool bVerb
)
238 restp_t
*restp
= NULL
;
242 int i
, j
, k
, l
, resind
, naid
, naaid
, natp
, nrestp
= 0;
244 for (i
= 0; (i
< atoms
->nres
); i
++)
246 if (gmx_strcasecmp(restype
[i
], "Protein") && gmx_strcasecmp(restype
[i
], "DNA") && gmx_strcasecmp(restype
[i
], "RNA") && gmx_strcasecmp(restype
[i
], "Water"))
256 printf("Analysing residues not classified as Protein/DNA/RNA/Water and splitting into groups...\n");
258 for (k
= 0; (k
< atoms
->nr
); k
++)
260 resind
= atoms
->atom
[k
].resind
;
261 rname
= *atoms
->resinfo
[resind
].name
;
262 if (gmx_strcasecmp(restype
[resind
], "Protein") && gmx_strcasecmp(restype
[resind
], "DNA") &&
263 gmx_strcasecmp(restype
[resind
], "RNA") && gmx_strcasecmp(restype
[resind
], "Water"))
266 for (l
= 0; (l
< nrestp
); l
++)
268 if (strcmp(restp
[l
].rname
, rname
) == 0)
275 srenew(restp
, nrestp
+1);
276 restp
[nrestp
].rname
= gmx_strdup(rname
);
277 restp
[nrestp
].bNeg
= FALSE
;
278 restp
[nrestp
].gname
= gmx_strdup(rname
);
283 for (i
= 0; (i
< nrestp
); i
++)
285 snew(aid
, atoms
->nr
);
287 for (j
= 0; (j
< atoms
->nr
); j
++)
289 rname
= *atoms
->resinfo
[atoms
->atom
[j
].resind
].name
;
290 if ((strcmp(restp
[i
].rname
, rname
) == 0 && !restp
[i
].bNeg
) ||
291 (strcmp(restp
[i
].rname
, rname
) != 0 && restp
[i
].bNeg
))
296 add_grp(gb
, gn
, naid
, aid
, restp
[i
].gname
);
299 printf("split %s into atoms (y/n) ? ", restp
[i
].gname
);
301 if (gmx_ask_yesno(bASK
))
304 for (k
= 0; (k
< naid
); k
++)
306 aname
= *atoms
->atomname
[aid
[k
]];
307 for (l
= 0; (l
< natp
); l
++)
309 if (strcmp(aname
, attp
[l
]) == 0)
316 srenew(attp
, ++natp
);
317 attp
[natp
-1] = aname
;
322 for (l
= 0; (l
< natp
); l
++)
326 for (k
= 0; (k
< naid
); k
++)
328 aname
= *atoms
->atomname
[aid
[k
]];
329 if (strcmp(aname
, attp
[l
]) == 0)
331 aaid
[naaid
++] = aid
[k
];
334 add_grp(gb
, gn
, naaid
, aaid
, attp
[l
]);
343 sfree(restp
[i
].rname
);
344 sfree(restp
[i
].gname
);
351 * Cata necessary to construct a single (protein) index group in
354 typedef struct gmx_help_make_index_group
356 /** The set of atom names that will be used to form this index group */
357 const char **defining_atomnames
;
358 /** Size of the defining_atomnames array */
359 int num_defining_atomnames
;
360 /** Name of this index group */
361 const char *group_name
;
362 /** Whether the above atom names name the atoms in the group, or
363 those not in the group */
364 gmx_bool bTakeComplement
;
365 /** The index in wholename gives the first item in the arrays of
366 atomnames that should be tested with 'gmx_strncasecmp' in stead of
367 gmx_strcasecmp, or -1 if all items should be tested with strcasecmp
368 This is comparable to using a '*' wildcard at the end of specific
369 atom names, but that is more involved to implement...
372 /** Only create this index group if it differs from the one specified in compareto,
373 where -1 means to always create this group. */
375 } t_gmx_help_make_index_group
;
377 static void analyse_prot(const char ** restype
, t_atoms
*atoms
,
378 t_blocka
*gb
, char ***gn
, gmx_bool bASK
, gmx_bool bVerb
)
380 /* lists of atomnames to be used in constructing index groups: */
381 static const char *pnoh
[] = { "H", "HN" };
382 static const char *pnodum
[] = {
383 "MN1", "MN2", "MCB1", "MCB2", "MCG1", "MCG2",
384 "MCD1", "MCD2", "MCE1", "MCE2", "MNZ1", "MNZ2"
386 static const char *calpha
[] = { "CA" };
387 static const char *bb
[] = { "N", "CA", "C" };
388 static const char *mc
[] = { "N", "CA", "C", "O", "O1", "O2", "OC1", "OC2", "OT", "OXT" };
389 static const char *mcb
[] = { "N", "CA", "CB", "C", "O", "O1", "O2", "OC1", "OC2", "OT", "OXT" };
390 static const char *mch
[] = {
391 "N", "CA", "C", "O", "O1", "O2", "OC1", "OC2", "OT", "OXT",
392 "H1", "H2", "H3", "H", "HN"
395 static const t_gmx_help_make_index_group constructing_data
[] =
396 {{ NULL
, 0, "Protein", TRUE
, -1, -1},
397 { pnoh
, asize(pnoh
), "Protein-H", TRUE
, 0, -1},
398 { calpha
, asize(calpha
), "C-alpha", FALSE
, -1, -1},
399 { bb
, asize(bb
), "Backbone", FALSE
, -1, -1},
400 { mc
, asize(mc
), "MainChain", FALSE
, -1, -1},
401 { mcb
, asize(mcb
), "MainChain+Cb", FALSE
, -1, -1},
402 { mch
, asize(mch
), "MainChain+H", FALSE
, -1, -1},
403 { mch
, asize(mch
), "SideChain", TRUE
, -1, -1},
404 { mch
, asize(mch
), "SideChain-H", TRUE
, 11, -1},
405 { pnodum
, asize(pnodum
), "Prot-Masses", TRUE
, -1, 0}, };
406 const int num_index_groups
= asize(constructing_data
);
412 char ndx_name
[STRLEN
], *atnm
;
417 printf("Analysing Protein...\n");
419 snew(aid
, atoms
->nr
);
421 /* calculate the number of protein residues */
423 for (i
= 0; (i
< atoms
->nres
); i
++)
425 if (0 == gmx_strcasecmp(restype
[i
], "Protein"))
430 /* find matching or complement atoms */
431 for (i
= 0; (i
< (int)num_index_groups
); i
++)
434 for (n
= 0; (n
< atoms
->nr
); n
++)
436 if (0 == gmx_strcasecmp(restype
[atoms
->atom
[n
].resind
], "Protein"))
439 for (j
= 0; (j
< constructing_data
[i
].num_defining_atomnames
); j
++)
441 /* skip digits at beginning of atomname, e.g. 1H */
442 atnm
= *atoms
->atomname
[n
];
443 while (isdigit(atnm
[0]))
447 if ( (constructing_data
[i
].wholename
== -1) || (j
< constructing_data
[i
].wholename
) )
449 if (0 == gmx_strcasecmp(constructing_data
[i
].defining_atomnames
[j
], atnm
))
456 if (0 == gmx_strncasecmp(constructing_data
[i
].defining_atomnames
[j
], atnm
, strlen(constructing_data
[i
].defining_atomnames
[j
])))
462 if (constructing_data
[i
].bTakeComplement
!= match
)
468 /* if we want to add this group always or it differs from previous
470 if (-1 == constructing_data
[i
].compareto
|| !grp_cmp(gb
, nra
, aid
, constructing_data
[i
].compareto
-i
) )
472 add_grp(gb
, gn
, nra
, aid
, constructing_data
[i
].group_name
);
478 for (i
= 0; (i
< (int)num_index_groups
); i
++)
480 printf("Split %12s into %5d residues (y/n) ? ", constructing_data
[i
].group_name
, npres
);
481 if (gmx_ask_yesno(bASK
))
485 for (n
= 0; ((atoms
->atom
[n
].resind
< npres
) && (n
< atoms
->nr
)); )
487 resind
= atoms
->atom
[n
].resind
;
488 for (; ((atoms
->atom
[n
].resind
== resind
) && (n
< atoms
->nr
)); n
++)
491 for (j
= 0; (j
< constructing_data
[i
].num_defining_atomnames
); j
++)
493 if (0 == gmx_strcasecmp(constructing_data
[i
].defining_atomnames
[j
], *atoms
->atomname
[n
]))
498 if (constructing_data
[i
].bTakeComplement
!= match
)
503 /* copy the residuename to the tail of the groupname */
507 ri
= &atoms
->resinfo
[resind
];
508 sprintf(ndx_name
, "%s_%s%d%c",
509 constructing_data
[i
].group_name
, *ri
->name
, ri
->nr
, ri
->ic
== ' ' ? '\0' : ri
->ic
);
510 add_grp(gb
, gn
, nra
, aid
, ndx_name
);
516 printf("Make group with sidechain and C=O swapped (y/n) ? ");
517 if (gmx_ask_yesno(bASK
))
519 /* Make swap sidechain C=O index */
522 for (n
= 0; ((atoms
->atom
[n
].resind
< npres
) && (n
< atoms
->nr
)); )
524 resind
= atoms
->atom
[n
].resind
;
526 for (; ((atoms
->atom
[n
].resind
== resind
) && (n
< atoms
->nr
)); n
++)
528 if (strcmp("CA", *atoms
->atomname
[n
]) == 0)
534 else if (strcmp("C", *atoms
->atomname
[n
]) == 0)
538 gmx_incons("Atom naming problem");
542 else if (strcmp("O", *atoms
->atomname
[n
]) == 0)
546 gmx_incons("Atom naming problem");
550 else if (strcmp("O1", *atoms
->atomname
[n
]) == 0)
554 gmx_incons("Atom naming problem");
564 /* copy the residuename to the tail of the groupname */
567 add_grp(gb
, gn
, nra
, aid
, "SwapSC-CO");
575 void analyse(t_atoms
*atoms
, t_blocka
*gb
, char ***gn
, gmx_bool bASK
, gmx_bool bVerb
)
577 gmx_residuetype_t
*rt
= NULL
;
580 const char ** restype
;
591 printf("Analysing residue names:\n");
593 /* Create system group, every single atom */
594 snew(aid
, atoms
->nr
);
595 for (i
= 0; i
< atoms
->nr
; i
++)
599 add_grp(gb
, gn
, atoms
->nr
, aid
, "System");
602 /* For every residue, get a pointer to the residue type name */
603 gmx_residuetype_init(&rt
);
606 snew(restype
, atoms
->nres
);
609 for (i
= 0; i
< atoms
->nres
; i
++)
611 resnm
= *atoms
->resinfo
[i
].name
;
612 gmx_residuetype_get_type(rt
, resnm
, &(restype
[i
]));
614 /* Note that this does not lead to a N*N loop, but N*K, where
615 * K is the number of residue _types_, which is small and independent of N.
618 for (k
= 0; k
< ntypes
&& !found
; k
++)
620 assert(p_typename
!= NULL
);
621 found
= !strcmp(restype
[i
], p_typename
[k
]);
625 srenew(p_typename
, ntypes
+1);
626 p_typename
[ntypes
++] = gmx_strdup(restype
[i
]);
632 p_status(restype
, atoms
->nres
, p_typename
, ntypes
);
635 for (k
= 0; k
< ntypes
; k
++)
637 aid
= mk_aid(atoms
, restype
, p_typename
[k
], &nra
, TRUE
);
639 /* Check for special types to do fancy stuff with */
641 if (!gmx_strcasecmp(p_typename
[k
], "Protein") && nra
> 0)
645 analyse_prot(restype
, atoms
, gb
, gn
, bASK
, bVerb
);
647 /* Create a Non-Protein group */
648 aid
= mk_aid(atoms
, restype
, "Protein", &nra
, FALSE
);
649 if ((nra
> 0) && (nra
< atoms
->nr
))
651 add_grp(gb
, gn
, nra
, aid
, "non-Protein");
655 else if (!gmx_strcasecmp(p_typename
[k
], "Water") && nra
> 0)
657 add_grp(gb
, gn
, nra
, aid
, p_typename
[k
]);
658 /* Add this group as 'SOL' too, for backward compatibility with older gromacs versions */
659 add_grp(gb
, gn
, nra
, aid
, "SOL");
663 /* Solvent, create a negated group too */
664 aid
= mk_aid(atoms
, restype
, "Water", &nra
, FALSE
);
665 if ((nra
> 0) && (nra
< atoms
->nr
))
667 add_grp(gb
, gn
, nra
, aid
, "non-Water");
674 add_grp(gb
, gn
, nra
, aid
, p_typename
[k
]);
676 analyse_other(restype
, atoms
, gb
, gn
, bASK
, bVerb
);
678 sfree(p_typename
[k
]);
683 gmx_residuetype_destroy(rt
);
685 /* Create a merged water_and_ions group */
691 for (i
= 0; i
< gb
->nr
; i
++)
693 if (!gmx_strcasecmp((*gn
)[i
], "Water"))
696 nwater
= gb
->index
[i
+1]-gb
->index
[i
];
698 else if (!gmx_strcasecmp((*gn
)[i
], "Ion"))
701 nion
= gb
->index
[i
+1]-gb
->index
[i
];
705 if (nwater
> 0 && nion
> 0)
707 srenew(gb
->index
, gb
->nr
+2);
708 srenew(*gn
, gb
->nr
+1);
709 (*gn
)[gb
->nr
] = gmx_strdup("Water_and_ions");
710 srenew(gb
->a
, gb
->nra
+nwater
+nion
);
713 for (i
= gb
->index
[iwater
]; i
< gb
->index
[iwater
+1]; i
++)
715 gb
->a
[gb
->nra
++] = gb
->a
[i
];
720 for (i
= gb
->index
[iion
]; i
< gb
->index
[iion
+1]; i
++)
722 gb
->a
[gb
->nra
++] = gb
->a
[i
];
726 gb
->index
[gb
->nr
] = gb
->nra
;
731 void check_index(char *gname
, int n
, int index
[], char *traj
, int natoms
)
735 for (i
= 0; i
< n
; i
++)
737 if (index
[i
] >= natoms
)
739 gmx_fatal(FARGS
, "%s atom number (index[%d]=%d) is larger than the number of atoms in %s (%d)",
740 gname
? gname
: "Index", i
+1, index
[i
]+1,
741 traj
? traj
: "the trajectory", natoms
);
743 else if (index
[i
] < 0)
745 gmx_fatal(FARGS
, "%s atom number (index[%d]=%d) is less than zero",
746 gname
? gname
: "Index", i
+1, index
[i
]+1);
751 t_blocka
*init_index(const char *gfile
, char ***grpname
)
757 char line
[STRLEN
], *pt
, str
[STRLEN
];
759 in
= gmx_fio_fopen(gfile
, "r");
767 while (get_a_line(in
, line
, STRLEN
))
769 if (get_header(line
, str
))
772 srenew(b
->index
, b
->nr
+1);
773 srenew(*grpname
, b
->nr
);
778 b
->index
[b
->nr
] = b
->index
[b
->nr
-1];
779 (*grpname
)[b
->nr
-1] = gmx_strdup(str
);
785 gmx_fatal(FARGS
, "The first header of your indexfile is invalid");
788 while (sscanf(pt
, "%s", str
) == 1)
794 srenew(b
->a
, maxentries
);
796 assert(b
->a
!= NULL
); // for clang analyzer
797 b
->a
[i
] = strtol(str
, NULL
, 10)-1;
800 pt
= strstr(pt
, str
)+strlen(str
);
806 for (i
= 0; (i
< b
->nr
); i
++)
808 assert(b
->a
!= NULL
); // for clang analyzer
809 for (j
= b
->index
[i
]; (j
< b
->index
[i
+1]); j
++)
813 fprintf(stderr
, "\nWARNING: negative index %d in group %s\n\n",
814 b
->a
[j
], (*grpname
)[i
]);
822 static void minstring(char *str
)
826 for (i
= 0; (i
< (int)strlen(str
)); i
++)
835 int find_group(const char *s
, int ngrps
, char **grpname
)
843 /* first look for whole name match */
846 for (i
= 0; i
< ngrps
; i
++)
848 if (gmx_strcasecmp_min(s
, grpname
[i
]) == 0)
858 /* second look for first string match */
861 for (i
= 0; i
< ngrps
; i
++)
863 if (gmx_strncasecmp_min(s
, grpname
[i
], n
) == 0)
873 /* last look for arbitrary substring match */
877 strncpy(key
, s
, sizeof(key
)-1);
878 key
[STRLEN
-1] = '\0';
881 for (i
= 0; i
< ngrps
; i
++)
883 strcpy(string
, grpname
[i
]);
886 if (strstr(string
, key
) != NULL
)
898 printf("Error: Multiple groups '%s' selected\n", s
);
904 static int qgroup(int *a
, int ngrps
, char **grpname
)
913 fprintf(stderr
, "Select a group: ");
916 if (scanf("%s", s
) != 1)
918 gmx_fatal(FARGS
, "Cannot read from input");
920 trim(s
); /* remove spaces */
922 while (strlen(s
) == 0);
923 aa
= strtol(s
, &end
, 10);
924 if (aa
== 0 && end
[0] != '\0') /* string entered */
926 aa
= find_group(s
, ngrps
, grpname
);
928 bInRange
= (aa
>= 0 && aa
< ngrps
);
931 printf("Error: No such group '%s'\n", s
);
935 printf("Selected %d: '%s'\n", aa
, grpname
[aa
]);
940 static void rd_groups(t_blocka
*grps
, char **grpname
, char *gnames
[],
941 int ngrps
, int isize
[], int *index
[], int grpnr
[])
947 gmx_fatal(FARGS
, "Error: no groups in indexfile");
949 for (i
= 0; (i
< grps
->nr
); i
++)
951 fprintf(stderr
, "Group %5d (%15s) has %5d elements\n", i
, grpname
[i
],
952 grps
->index
[i
+1]-grps
->index
[i
]);
954 for (i
= 0; (i
< ngrps
); i
++)
960 gnr1
= qgroup(&grpnr
[i
], grps
->nr
, grpname
);
961 if ((gnr1
< 0) || (gnr1
>= grps
->nr
))
963 fprintf(stderr
, "Select between %d and %d.\n", 0, grps
->nr
-1);
966 while ((gnr1
< 0) || (gnr1
>= grps
->nr
));
970 fprintf(stderr
, "There is one group in the index\n");
973 gnames
[i
] = gmx_strdup(grpname
[gnr1
]);
974 isize
[i
] = grps
->index
[gnr1
+1]-grps
->index
[gnr1
];
975 snew(index
[i
], isize
[i
]);
976 for (j
= 0; (j
< isize
[i
]); j
++)
978 index
[i
][j
] = grps
->a
[grps
->index
[gnr1
]+j
];
983 void rd_index(const char *statfile
, int ngrps
, int isize
[],
984 int *index
[], char *grpnames
[])
993 gmx_fatal(FARGS
, "No index file specified");
995 grps
= init_index(statfile
, &gnames
);
996 rd_groups(grps
, gnames
, grpnames
, ngrps
, isize
, index
, grpnr
);
999 void rd_index_nrs(char *statfile
, int ngrps
, int isize
[],
1000 int *index
[], char *grpnames
[], int grpnr
[])
1007 gmx_fatal(FARGS
, "No index file specified");
1009 grps
= init_index(statfile
, &gnames
);
1011 rd_groups(grps
, gnames
, grpnames
, ngrps
, isize
, index
, grpnr
);
1014 void get_index(t_atoms
*atoms
, const char *fnm
, int ngrps
,
1015 int isize
[], int *index
[], char *grpnames
[])
1018 t_blocka
*grps
= NULL
;
1025 grps
= init_index(fnm
, gnames
);
1030 snew(grps
->index
, 1);
1031 analyse(atoms
, grps
, gnames
, FALSE
, FALSE
);
1035 gmx_incons("You need to supply a valid atoms structure or a valid index file name");
1038 rd_groups(grps
, *gnames
, grpnames
, ngrps
, isize
, index
, grpnr
);
1041 t_cluster_ndx
*cluster_index(FILE *fplog
, const char *ndx
)
1047 c
->clust
= init_index(ndx
, &c
->grpname
);
1049 for (i
= 0; (i
< c
->clust
->nra
); i
++)
1051 c
->maxframe
= std::max(c
->maxframe
, c
->clust
->a
[i
]);
1053 fprintf(fplog
? fplog
: stdout
,
1054 "There are %d clusters containing %d structures, highest framenr is %d\n",
1055 c
->clust
->nr
, c
->clust
->nra
, c
->maxframe
);
1058 pr_blocka(debug
, 0, "clust", c
->clust
, TRUE
);
1059 for (i
= 0; (i
< c
->clust
->nra
); i
++)
1061 if ((c
->clust
->a
[i
] < 0) || (c
->clust
->a
[i
] > c
->maxframe
))
1063 gmx_fatal(FARGS
, "Range check error for c->clust->a[%d] = %d\n"
1064 "should be within 0 and %d", i
, c
->clust
->a
[i
], c
->maxframe
+1);
1068 c
->inv_clust
= make_invblocka(c
->clust
, c
->maxframe
);