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, 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 ]\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
++)
269 if (strcmp(restp
[l
].rname
, rname
) == 0)
276 srenew(restp
, nrestp
+1);
277 restp
[nrestp
].rname
= gmx_strdup(rname
);
278 restp
[nrestp
].bNeg
= FALSE
;
279 restp
[nrestp
].gname
= gmx_strdup(rname
);
284 for (i
= 0; (i
< nrestp
); i
++)
286 snew(aid
, atoms
->nr
);
288 for (j
= 0; (j
< atoms
->nr
); j
++)
290 rname
= *atoms
->resinfo
[atoms
->atom
[j
].resind
].name
;
291 if ((strcmp(restp
[i
].rname
, rname
) == 0 && !restp
[i
].bNeg
) ||
292 (strcmp(restp
[i
].rname
, rname
) != 0 && restp
[i
].bNeg
))
297 add_grp(gb
, gn
, naid
, aid
, restp
[i
].gname
);
300 printf("split %s into atoms (y/n) ? ", restp
[i
].gname
);
302 if (gmx_ask_yesno(bASK
))
305 for (k
= 0; (k
< naid
); k
++)
307 aname
= *atoms
->atomname
[aid
[k
]];
308 for (l
= 0; (l
< natp
); l
++)
310 if (strcmp(aname
, attp
[l
]) == 0)
317 srenew(attp
, ++natp
);
318 attp
[natp
-1] = aname
;
323 for (l
= 0; (l
< natp
); l
++)
327 for (k
= 0; (k
< naid
); k
++)
329 aname
= *atoms
->atomname
[aid
[k
]];
330 if (strcmp(aname
, attp
[l
]) == 0)
332 aaid
[naaid
++] = aid
[k
];
335 add_grp(gb
, gn
, naaid
, aaid
, attp
[l
]);
344 sfree(restp
[i
].rname
);
345 sfree(restp
[i
].gname
);
352 * Cata necessary to construct a single (protein) index group in
355 typedef struct gmx_help_make_index_group
357 /** The set of atom names that will be used to form this index group */
358 const char **defining_atomnames
;
359 /** Size of the defining_atomnames array */
360 int num_defining_atomnames
;
361 /** Name of this index group */
362 const char *group_name
;
363 /** Whether the above atom names name the atoms in the group, or
364 those not in the group */
365 gmx_bool bTakeComplement
;
366 /** The index in wholename gives the first item in the arrays of
367 atomnames that should be tested with 'gmx_strncasecmp' in stead of
368 gmx_strcasecmp, or -1 if all items should be tested with strcasecmp
369 This is comparable to using a '*' wildcard at the end of specific
370 atom names, but that is more involved to implement...
373 /** Only create this index group if it differs from the one specified in compareto,
374 where -1 means to always create this group. */
376 } t_gmx_help_make_index_group
;
378 static void analyse_prot(const char ** restype
, t_atoms
*atoms
,
379 t_blocka
*gb
, char ***gn
, gmx_bool bASK
, gmx_bool bVerb
)
381 /* lists of atomnames to be used in constructing index groups: */
382 static const char *pnoh
[] = { "H", "HN" };
383 static const char *pnodum
[] = {
384 "MN1", "MN2", "MCB1", "MCB2", "MCG1", "MCG2",
385 "MCD1", "MCD2", "MCE1", "MCE2", "MNZ1", "MNZ2"
387 static const char *calpha
[] = { "CA" };
388 static const char *bb
[] = { "N", "CA", "C" };
389 static const char *mc
[] = { "N", "CA", "C", "O", "O1", "O2", "OC1", "OC2", "OT", "OXT" };
390 static const char *mcb
[] = { "N", "CA", "CB", "C", "O", "O1", "O2", "OC1", "OC2", "OT", "OXT" };
391 static const char *mch
[] = {
392 "N", "CA", "C", "O", "O1", "O2", "OC1", "OC2", "OT", "OXT",
393 "H1", "H2", "H3", "H", "HN"
396 static const t_gmx_help_make_index_group constructing_data
[] =
397 {{ NULL
, 0, "Protein", TRUE
, -1, -1},
398 { pnoh
, asize(pnoh
), "Protein-H", TRUE
, 0, -1},
399 { calpha
, asize(calpha
), "C-alpha", FALSE
, -1, -1},
400 { bb
, asize(bb
), "Backbone", FALSE
, -1, -1},
401 { mc
, asize(mc
), "MainChain", FALSE
, -1, -1},
402 { mcb
, asize(mcb
), "MainChain+Cb", FALSE
, -1, -1},
403 { mch
, asize(mch
), "MainChain+H", FALSE
, -1, -1},
404 { mch
, asize(mch
), "SideChain", TRUE
, -1, -1},
405 { mch
, asize(mch
), "SideChain-H", TRUE
, 11, -1},
406 { pnodum
, asize(pnodum
), "Prot-Masses", TRUE
, -1, 0}, };
407 const int num_index_groups
= asize(constructing_data
);
413 char ndx_name
[STRLEN
], *atnm
;
418 printf("Analysing Protein...\n");
420 snew(aid
, atoms
->nr
);
422 /* calculate the number of protein residues */
424 for (i
= 0; (i
< atoms
->nres
); i
++)
426 if (0 == gmx_strcasecmp(restype
[i
], "Protein"))
431 /* find matching or complement atoms */
432 for (i
= 0; (i
< (int)num_index_groups
); i
++)
435 for (n
= 0; (n
< atoms
->nr
); n
++)
437 if (0 == gmx_strcasecmp(restype
[atoms
->atom
[n
].resind
], "Protein"))
440 for (j
= 0; (j
< constructing_data
[i
].num_defining_atomnames
); j
++)
442 /* skip digits at beginning of atomname, e.g. 1H */
443 atnm
= *atoms
->atomname
[n
];
444 while (isdigit(atnm
[0]))
448 if ( (constructing_data
[i
].wholename
== -1) || (j
< constructing_data
[i
].wholename
) )
450 if (0 == gmx_strcasecmp(constructing_data
[i
].defining_atomnames
[j
], atnm
))
457 if (0 == gmx_strncasecmp(constructing_data
[i
].defining_atomnames
[j
], atnm
, strlen(constructing_data
[i
].defining_atomnames
[j
])))
463 if (constructing_data
[i
].bTakeComplement
!= match
)
469 /* if we want to add this group always or it differs from previous
471 if (-1 == constructing_data
[i
].compareto
|| !grp_cmp(gb
, nra
, aid
, constructing_data
[i
].compareto
-i
) )
473 add_grp(gb
, gn
, nra
, aid
, constructing_data
[i
].group_name
);
479 for (i
= 0; (i
< (int)num_index_groups
); i
++)
481 printf("Split %12s into %5d residues (y/n) ? ", constructing_data
[i
].group_name
, npres
);
482 if (gmx_ask_yesno(bASK
))
486 for (n
= 0; ((atoms
->atom
[n
].resind
< npres
) && (n
< atoms
->nr
)); )
488 resind
= atoms
->atom
[n
].resind
;
489 for (; ((atoms
->atom
[n
].resind
== resind
) && (n
< atoms
->nr
)); n
++)
492 for (j
= 0; (j
< constructing_data
[i
].num_defining_atomnames
); j
++)
494 if (0 == gmx_strcasecmp(constructing_data
[i
].defining_atomnames
[j
], *atoms
->atomname
[n
]))
499 if (constructing_data
[i
].bTakeComplement
!= match
)
504 /* copy the residuename to the tail of the groupname */
508 ri
= &atoms
->resinfo
[resind
];
509 sprintf(ndx_name
, "%s_%s%d%c",
510 constructing_data
[i
].group_name
, *ri
->name
, ri
->nr
, ri
->ic
== ' ' ? '\0' : ri
->ic
);
511 add_grp(gb
, gn
, nra
, aid
, ndx_name
);
517 printf("Make group with sidechain and C=O swapped (y/n) ? ");
518 if (gmx_ask_yesno(bASK
))
520 /* Make swap sidechain C=O index */
523 for (n
= 0; ((atoms
->atom
[n
].resind
< npres
) && (n
< atoms
->nr
)); )
525 resind
= atoms
->atom
[n
].resind
;
527 for (; ((atoms
->atom
[n
].resind
== resind
) && (n
< atoms
->nr
)); n
++)
529 if (strcmp("CA", *atoms
->atomname
[n
]) == 0)
535 else if (strcmp("C", *atoms
->atomname
[n
]) == 0)
539 gmx_incons("Atom naming problem");
543 else if (strcmp("O", *atoms
->atomname
[n
]) == 0)
547 gmx_incons("Atom naming problem");
551 else if (strcmp("O1", *atoms
->atomname
[n
]) == 0)
555 gmx_incons("Atom naming problem");
565 /* copy the residuename to the tail of the groupname */
568 add_grp(gb
, gn
, nra
, aid
, "SwapSC-CO");
576 void analyse(t_atoms
*atoms
, t_blocka
*gb
, char ***gn
, gmx_bool bASK
, gmx_bool bVerb
)
578 gmx_residuetype_t
*rt
= NULL
;
581 const char ** restype
;
592 printf("Analysing residue names:\n");
594 /* Create system group, every single atom */
595 snew(aid
, atoms
->nr
);
596 for (i
= 0; i
< atoms
->nr
; i
++)
600 add_grp(gb
, gn
, atoms
->nr
, aid
, "System");
603 /* For every residue, get a pointer to the residue type name */
604 gmx_residuetype_init(&rt
);
607 snew(restype
, atoms
->nres
);
614 resnm
= *atoms
->resinfo
[i
].name
;
615 gmx_residuetype_get_type(rt
, resnm
, &(restype
[i
]));
616 snew(p_typename
, ntypes
+1);
617 p_typename
[ntypes
] = gmx_strdup(restype
[i
]);
620 for (i
= 1; i
< atoms
->nres
; i
++)
622 resnm
= *atoms
->resinfo
[i
].name
;
623 gmx_residuetype_get_type(rt
, resnm
, &(restype
[i
]));
625 /* Note that this does not lead to a N*N loop, but N*K, where
626 * K is the number of residue _types_, which is small and independent of N.
629 for (k
= 0; k
< ntypes
&& !found
; k
++)
631 found
= !strcmp(restype
[i
], p_typename
[k
]);
635 srenew(p_typename
, ntypes
+1);
636 p_typename
[ntypes
] = gmx_strdup(restype
[i
]);
644 p_status(restype
, atoms
->nres
, p_typename
, ntypes
);
647 for (k
= 0; k
< ntypes
; k
++)
649 aid
= mk_aid(atoms
, restype
, p_typename
[k
], &nra
, TRUE
);
651 /* Check for special types to do fancy stuff with */
653 if (!gmx_strcasecmp(p_typename
[k
], "Protein") && nra
> 0)
657 analyse_prot(restype
, atoms
, gb
, gn
, bASK
, bVerb
);
659 /* Create a Non-Protein group */
660 aid
= mk_aid(atoms
, restype
, "Protein", &nra
, FALSE
);
661 if ((nra
> 0) && (nra
< atoms
->nr
))
663 add_grp(gb
, gn
, nra
, aid
, "non-Protein");
667 else if (!gmx_strcasecmp(p_typename
[k
], "Water") && nra
> 0)
669 add_grp(gb
, gn
, nra
, aid
, p_typename
[k
]);
670 /* Add this group as 'SOL' too, for backward compatibility with older gromacs versions */
671 add_grp(gb
, gn
, nra
, aid
, "SOL");
675 /* Solvent, create a negated group too */
676 aid
= mk_aid(atoms
, restype
, "Water", &nra
, FALSE
);
677 if ((nra
> 0) && (nra
< atoms
->nr
))
679 add_grp(gb
, gn
, nra
, aid
, "non-Water");
686 add_grp(gb
, gn
, nra
, aid
, p_typename
[k
]);
688 analyse_other(restype
, atoms
, gb
, gn
, bASK
, bVerb
);
690 sfree(p_typename
[k
]);
695 gmx_residuetype_destroy(rt
);
697 /* Create a merged water_and_ions group */
703 for (i
= 0; i
< gb
->nr
; i
++)
705 if (!gmx_strcasecmp((*gn
)[i
], "Water"))
708 nwater
= gb
->index
[i
+1]-gb
->index
[i
];
710 else if (!gmx_strcasecmp((*gn
)[i
], "Ion"))
713 nion
= gb
->index
[i
+1]-gb
->index
[i
];
717 if (nwater
> 0 && nion
> 0)
719 srenew(gb
->index
, gb
->nr
+2);
720 srenew(*gn
, gb
->nr
+1);
721 (*gn
)[gb
->nr
] = gmx_strdup("Water_and_ions");
722 srenew(gb
->a
, gb
->nra
+nwater
+nion
);
725 for (i
= gb
->index
[iwater
]; i
< gb
->index
[iwater
+1]; i
++)
727 gb
->a
[gb
->nra
++] = gb
->a
[i
];
732 for (i
= gb
->index
[iion
]; i
< gb
->index
[iion
+1]; i
++)
734 gb
->a
[gb
->nra
++] = gb
->a
[i
];
738 gb
->index
[gb
->nr
] = gb
->nra
;
743 void check_index(char *gname
, int n
, int index
[], char *traj
, int natoms
)
747 for (i
= 0; i
< n
; i
++)
749 if (index
[i
] >= natoms
)
751 gmx_fatal(FARGS
, "%s atom number (index[%d]=%d) is larger than the number of atoms in %s (%d)",
752 gname
? gname
: "Index", i
+1, index
[i
]+1,
753 traj
? traj
: "the trajectory", natoms
);
755 else if (index
[i
] < 0)
757 gmx_fatal(FARGS
, "%s atom number (index[%d]=%d) is less than zero",
758 gname
? gname
: "Index", i
+1, index
[i
]+1);
763 t_blocka
*init_index(const char *gfile
, char ***grpname
)
769 char line
[STRLEN
], *pt
, str
[STRLEN
];
771 in
= gmx_ffopen(gfile
, "r");
779 while (get_a_line(in
, line
, STRLEN
))
781 if (get_header(line
, str
))
784 srenew(b
->index
, b
->nr
+1);
785 srenew(*grpname
, b
->nr
);
790 b
->index
[b
->nr
] = b
->index
[b
->nr
-1];
791 (*grpname
)[b
->nr
-1] = gmx_strdup(str
);
797 gmx_fatal(FARGS
, "The first header of your indexfile is invalid");
800 while (sscanf(pt
, "%s", str
) == 1)
806 srenew(b
->a
, maxentries
);
808 assert(b
->a
!= NULL
); // for clang analyzer
809 b
->a
[i
] = strtol(str
, NULL
, 10)-1;
812 pt
= strstr(pt
, str
)+strlen(str
);
818 for (i
= 0; (i
< b
->nr
); i
++)
820 assert(b
->a
!= NULL
); // for clang analyzer
821 for (j
= b
->index
[i
]; (j
< b
->index
[i
+1]); j
++)
825 fprintf(stderr
, "\nWARNING: negative index %d in group %s\n\n",
826 b
->a
[j
], (*grpname
)[i
]);
834 static void minstring(char *str
)
838 for (i
= 0; (i
< (int)strlen(str
)); i
++)
847 int find_group(const char *s
, int ngrps
, char **grpname
)
855 /* first look for whole name match */
858 for (i
= 0; i
< ngrps
; i
++)
860 if (gmx_strcasecmp_min(s
, grpname
[i
]) == 0)
870 /* second look for first string match */
873 for (i
= 0; i
< ngrps
; i
++)
875 if (gmx_strncasecmp_min(s
, grpname
[i
], n
) == 0)
885 /* last look for arbitrary substring match */
889 strncpy(key
, s
, sizeof(key
)-1);
890 key
[STRLEN
-1] = '\0';
893 for (i
= 0; i
< ngrps
; i
++)
895 strcpy(string
, grpname
[i
]);
898 if (strstr(string
, key
) != NULL
)
910 printf("Error: Multiple groups '%s' selected\n", s
);
916 static int qgroup(int *a
, int ngrps
, char **grpname
)
925 fprintf(stderr
, "Select a group: ");
928 if (scanf("%s", s
) != 1)
930 gmx_fatal(FARGS
, "Cannot read from input");
932 trim(s
); /* remove spaces */
934 while (strlen(s
) == 0);
935 aa
= strtol(s
, &end
, 10);
936 if (aa
== 0 && end
[0] != '\0') /* string entered */
938 aa
= find_group(s
, ngrps
, grpname
);
940 bInRange
= (aa
>= 0 && aa
< ngrps
);
943 printf("Error: No such group '%s'\n", s
);
947 printf("Selected %d: '%s'\n", aa
, grpname
[aa
]);
952 static void rd_groups(t_blocka
*grps
, char **grpname
, char *gnames
[],
953 int ngrps
, int isize
[], int *index
[], int grpnr
[])
959 gmx_fatal(FARGS
, "Error: no groups in indexfile");
961 for (i
= 0; (i
< grps
->nr
); i
++)
963 fprintf(stderr
, "Group %5d (%15s) has %5d elements\n", i
, grpname
[i
],
964 grps
->index
[i
+1]-grps
->index
[i
]);
966 for (i
= 0; (i
< ngrps
); i
++)
972 gnr1
= qgroup(&grpnr
[i
], grps
->nr
, grpname
);
973 if ((gnr1
< 0) || (gnr1
>= grps
->nr
))
975 fprintf(stderr
, "Select between %d and %d.\n", 0, grps
->nr
-1);
978 while ((gnr1
< 0) || (gnr1
>= grps
->nr
));
982 fprintf(stderr
, "There is one group in the index\n");
985 gnames
[i
] = gmx_strdup(grpname
[gnr1
]);
986 isize
[i
] = grps
->index
[gnr1
+1]-grps
->index
[gnr1
];
987 snew(index
[i
], isize
[i
]);
988 for (j
= 0; (j
< isize
[i
]); j
++)
990 index
[i
][j
] = grps
->a
[grps
->index
[gnr1
]+j
];
995 void rd_index(const char *statfile
, int ngrps
, int isize
[],
996 int *index
[], char *grpnames
[])
1005 gmx_fatal(FARGS
, "No index file specified");
1007 grps
= init_index(statfile
, &gnames
);
1008 rd_groups(grps
, gnames
, grpnames
, ngrps
, isize
, index
, grpnr
);
1011 void rd_index_nrs(char *statfile
, int ngrps
, int isize
[],
1012 int *index
[], char *grpnames
[], int grpnr
[])
1019 gmx_fatal(FARGS
, "No index file specified");
1021 grps
= init_index(statfile
, &gnames
);
1023 rd_groups(grps
, gnames
, grpnames
, ngrps
, isize
, index
, grpnr
);
1026 void get_index(t_atoms
*atoms
, const char *fnm
, int ngrps
,
1027 int isize
[], int *index
[], char *grpnames
[])
1030 t_blocka
*grps
= NULL
;
1037 grps
= init_index(fnm
, gnames
);
1042 snew(grps
->index
, 1);
1043 analyse(atoms
, grps
, gnames
, FALSE
, FALSE
);
1047 gmx_incons("You need to supply a valid atoms structure or a valid index file name");
1050 rd_groups(grps
, *gnames
, grpnames
, ngrps
, isize
, index
, grpnr
);
1053 t_cluster_ndx
*cluster_index(FILE *fplog
, const char *ndx
)
1059 c
->clust
= init_index(ndx
, &c
->grpname
);
1061 for (i
= 0; (i
< c
->clust
->nra
); i
++)
1063 c
->maxframe
= std::max(c
->maxframe
, c
->clust
->a
[i
]);
1065 fprintf(fplog
? fplog
: stdout
,
1066 "There are %d clusters containing %d structures, highest framenr is %d\n",
1067 c
->clust
->nr
, c
->clust
->nra
, c
->maxframe
);
1070 pr_blocka(debug
, 0, "clust", c
->clust
, TRUE
);
1071 for (i
= 0; (i
< c
->clust
->nra
); i
++)
1073 if ((c
->clust
->a
[i
] < 0) || (c
->clust
->a
[i
] > c
->maxframe
))
1075 gmx_fatal(FARGS
, "Range check error for c->clust->a[%d] = %d\n"
1076 "should be within 0 and %d", i
, c
->clust
->a
[i
], c
->maxframe
+1);
1080 c
->inv_clust
= make_invblocka(c
->clust
, c
->maxframe
);