3 * This source code is part of
7 * GROningen MAchine for Chemical Simulations
10 * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
11 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
12 * Copyright (c) 2001-2004, The GROMACS development team,
13 * check out http://www.gromacs.org for more information.
15 * This program is free software; you can redistribute it and/or
16 * modify it under the terms of the GNU General Public License
17 * as published by the Free Software Foundation; either version 2
18 * of the License, or (at your option) any later version.
20 * If you want to redistribute modifications, please consider that
21 * scientific software is very special. Version control is crucial -
22 * bugs must be traceable. We will be happy to consider code for
23 * inclusion in the official distribution, but derived work must not
24 * be called official GROMACS. Details are found in the README & COPYING
25 * files - if they are missing, get the official version at www.gromacs.org.
27 * To help us fund GROMACS development, we humbly ask that you cite
28 * the papers on the package - you can find them in the top README file.
30 * For more info, check our website at http://www.gromacs.org
33 * GROningen Mixture of Alchemy and Childrens' Stories
61 const char gmx_residuetype_undefined
[]="Other";
63 struct gmx_residuetype
72 static gmx_bool
gmx_ask_yesno(gmx_bool bASK
)
78 c
=toupper(fgetc(stdin
));
79 } while ((c
!= 'Y') && (c
!= 'N'));
87 t_blocka
*new_blocka(void)
97 void write_index(const char *outf
, t_blocka
*b
,char **gnames
)
102 out
=gmx_fio_fopen(outf
,"w");
103 /* fprintf(out,"%5d %5d\n",b->nr,b->nra); */
104 for(i
=0; (i
<b
->nr
); i
++) {
105 fprintf(out
,"[ %s ]\n",gnames
[i
]);
106 for(k
=0,j
=b
->index
[i
]; j
<b
->index
[i
+1]; j
++,k
++) {
107 fprintf(out
,"%4d ",b
->a
[j
]+1);
116 void add_grp(t_blocka
*b
,char ***gnames
,int nra
,atom_id a
[],const char *name
)
120 srenew(b
->index
,b
->nr
+2);
121 srenew(*gnames
,b
->nr
+1);
122 (*gnames
)[b
->nr
]=strdup(name
);
124 srenew(b
->a
,b
->nra
+nra
);
125 for(i
=0; (i
<nra
); i
++)
128 b
->index
[b
->nr
]=b
->nra
;
131 /* compare index in `a' with group in `b' at `index',
132 when `index'<0 it is relative to end of `b' */
133 static gmx_bool
grp_cmp(t_blocka
*b
, int nra
, atom_id a
[], int index
)
138 index
= b
->nr
-1+index
;
140 gmx_fatal(FARGS
,"no such index group %d in t_blocka (nr=%d)",index
,b
->nr
);
142 if ( nra
!= b
->index
[index
+1] - b
->index
[index
] )
145 if ( a
[i
] != b
->a
[b
->index
[index
]+i
] )
151 p_status(const char **restype
, int nres
, const char **typenames
, int ntypes
, gmx_bool bVerb
)
158 snew(counter
,ntypes
);
159 for(i
=0;i
<ntypes
;i
++)
166 for(j
=0;j
<ntypes
;j
++)
168 if(!gmx_strcasecmp(restype
[i
],typenames
[j
]))
177 for(i
=0; (i
<ntypes
); i
++)
181 printf("There are: %5d %10s residues\n",counter
[i
],typenames
[i
]);
189 mk_aid(t_atoms
*atoms
,const char ** restype
,const char * typestring
,int *nra
,gmx_bool bMatch
)
190 /* Make an array of atom_ids for all atoms with residuetypes matching typestring, or the opposite if bMatch is false */
198 for(i
=0; (i
<atoms
->nr
); i
++)
200 res
=!gmx_strcasecmp(restype
[atoms
->atom
[i
].resind
],typestring
);
220 static void analyse_other(const char ** restype
,t_atoms
*atoms
,
221 t_blocka
*gb
,char ***gn
,gmx_bool bASK
,gmx_bool bVerb
)
226 atom_id
*other_ndx
,*aid
,*aaid
;
227 int i
,j
,k
,l
,resind
,naid
,naaid
,natp
,nrestp
=0;
229 for(i
=0; (i
<atoms
->nres
); i
++)
231 if (gmx_strcasecmp(restype
[i
],"Protein") && gmx_strcasecmp(restype
[i
],"DNA") && gmx_strcasecmp(restype
[i
],"RNA") && gmx_strcasecmp(restype
[i
],"Water"))
236 if (i
< atoms
->nres
) {
239 printf("Analysing residues not classified as Protein/DNA/RNA/Water and splitting into groups...\n");
240 snew(other_ndx
,atoms
->nr
);
241 for(k
=0; (k
<atoms
->nr
); k
++) {
242 resind
= atoms
->atom
[k
].resind
;
243 rname
= *atoms
->resinfo
[resind
].name
;
244 if (gmx_strcasecmp(restype
[resind
],"Protein") && gmx_strcasecmp(restype
[resind
],"DNA") &&
245 gmx_strcasecmp(restype
[resind
],"RNA") && gmx_strcasecmp(restype
[resind
],"Water"))
248 for(l
=0; (l
<nrestp
); l
++)
249 if (strcmp(restp
[l
].rname
,rname
) == 0)
252 srenew(restp
,nrestp
+1);
253 restp
[nrestp
].rname
= strdup(rname
);
254 restp
[nrestp
].bNeg
= FALSE
;
255 restp
[nrestp
].gname
= strdup(rname
);
261 for(i
=0; (i
<nrestp
); i
++) {
264 for(j
=0; (j
<atoms
->nr
); j
++) {
265 rname
= *atoms
->resinfo
[atoms
->atom
[j
].resind
].name
;
266 if ((strcmp(restp
[i
].rname
,rname
) == 0 && !restp
[i
].bNeg
) ||
267 (strcmp(restp
[i
].rname
,rname
) != 0 && restp
[i
].bNeg
)) {
271 add_grp(gb
,gn
,naid
,aid
,restp
[i
].gname
);
273 printf("split %s into atoms (y/n) ? ",restp
[i
].gname
);
275 if (gmx_ask_yesno(bASK
)) {
277 for(k
=0; (k
<naid
); k
++) {
278 aname
=*atoms
->atomname
[aid
[k
]];
279 for(l
=0; (l
<natp
); l
++)
280 if (strcmp(aname
,attp
[l
]) == 0)
288 for(l
=0; (l
<natp
); l
++) {
291 for(k
=0; (k
<naid
); k
++) {
292 aname
=*atoms
->atomname
[aid
[k
]];
293 if (strcmp(aname
,attp
[l
])==0)
294 aaid
[naaid
++]=aid
[k
];
296 add_grp(gb
,gn
,naaid
,aaid
,attp
[l
]);
310 static void analyse_prot(const char ** restype
,t_atoms
*atoms
,
311 t_blocka
*gb
,char ***gn
,gmx_bool bASK
,gmx_bool bVerb
)
313 /* atomnames to be used in constructing index groups: */
314 static const char *pnoh
[] = { "H" };
315 static const char *pnodum
[] = { "MN1", "MN2", "MCB1", "MCB2", "MCG1", "MCG2",
316 "MCD1", "MCD2", "MCE1", "MCE2", "MNZ1", "MNZ2" };
317 static const char *calpha
[] = { "CA" };
318 static const char *bb
[] = { "N","CA","C" };
319 static const char *mc
[] = { "N","CA","C","O","O1","O2","OC1","OC2","OT","OXT" };
320 static const char *mcb
[] = { "N","CA","CB","C","O","O1","O2","OC1","OC2","OT","OXT" };
321 static const char *mch
[] = { "N","CA","C","O","O1","O2","OC1","OC2","OT","OXT",
322 "H1","H2","H3","H" };
323 /* array of arrays of atomnames: */
324 static const char **chains
[] = { NULL
,pnoh
,calpha
,bb
,mc
,mcb
,mch
,mch
,mch
,pnodum
};
325 #define NCH asize(chains)
326 /* array of sizes of arrays of atomnames: */
327 const int sizes
[NCH
] = {
328 0, asize(pnoh
), asize(calpha
), asize(bb
),
329 asize(mc
), asize(mcb
), asize(mch
), asize(mch
), asize(mch
), asize(pnodum
)
331 /* descriptive names of index groups */
332 const char *ch_name
[NCH
] = {
333 "Protein", "Protein-H", "C-alpha", "Backbone",
334 "MainChain", "MainChain+Cb", "MainChain+H", "SideChain", "SideChain-H",
337 /* construct index group containing (TRUE) or excluding (FALSE)
339 const gmx_bool complement
[NCH
] = {
340 TRUE
, TRUE
, FALSE
, FALSE
, FALSE
, FALSE
, FALSE
, TRUE
, TRUE
, TRUE
342 const int wholename
[NCH
] = { -1, 0,-1,-1,-1,-1,-1,-1, 11,-1 };
343 /* the index in wholename gives the first item in the arrays of
344 * atomtypes that should be tested with 'gmx_strncasecmp' in stead of
345 * gmx_strcasecmp, or -1 if all items should be tested with strcasecmp
346 * This is comparable to using a '*' wildcard at the end of specific
347 * atom names, but that is more involved to implement...
349 /* only add index group if it differs from the specified one,
350 specify -1 to always add group */
351 const int compareto
[NCH
] = { -1,-1,-1,-1,-1,-1,-1,-1,-1, 0 };
355 int nra
,nnpres
,npres
;
357 char ndx_name
[STRLEN
],*atnm
;
362 printf("Analysing Protein...\n");
366 /* calculate the number of protein residues */
368 for(i
=0; (i
<atoms
->nres
); i
++)
369 if (!gmx_strcasecmp(restype
[i
],"Protein"))
373 /* find matching or complement atoms */
374 for(i
=0; (i
<(int)NCH
); i
++) {
376 for(n
=0; (n
<atoms
->nr
); n
++) {
377 if (!gmx_strcasecmp(restype
[atoms
->atom
[n
].resind
],"Protein")) {
380 for(j
=0; (j
<sizes
[i
]); j
++) {
381 /* skip digits at beginning of atomname, e.g. 1H */
382 atnm
=*atoms
->atomname
[n
];
383 while (isdigit(atnm
[0]))
385 if ( (wholename
[i
]==-1) || (j
<wholename
[i
]) ) {
386 if (gmx_strcasecmp(chains
[i
][j
],atnm
) == 0)
389 if (gmx_strncasecmp(chains
[i
][j
],atnm
,strlen(chains
[i
][j
])) == 0)
393 if (match
!= complement
[i
])
397 /* if we want to add this group always or it differs from previous
399 if ( compareto
[i
] == -1 || !grp_cmp(gb
,nra
,aid
,compareto
[i
]-i
) )
400 add_grp(gb
,gn
,nra
,aid
,ch_name
[i
]);
404 for(i
=0; (i
<(int)NCH
); i
++) {
405 printf("Split %12s into %5d residues (y/n) ? ",ch_name
[i
],npres
);
406 if (gmx_ask_yesno(bASK
)) {
409 for(n
=0;((atoms
->atom
[n
].resind
< npres
) && (n
<atoms
->nr
));) {
410 resind
= atoms
->atom
[n
].resind
;
411 for(;((atoms
->atom
[n
].resind
==resind
) && (n
<atoms
->nr
));n
++) {
413 for(j
=0;(j
<sizes
[i
]); j
++)
414 if (gmx_strcasecmp(chains
[i
][j
],*atoms
->atomname
[n
]) == 0)
416 if (match
!= complement
[i
])
419 /* copy the residuename to the tail of the groupname */
422 ri
= &atoms
->resinfo
[resind
];
423 sprintf(ndx_name
,"%s_%s%d%c",
424 ch_name
[i
],*ri
->name
,ri
->nr
,ri
->ic
==' ' ? '\0' : ri
->ic
);
425 add_grp(gb
,gn
,nra
,aid
,ndx_name
);
431 printf("Make group with sidechain and C=O swapped (y/n) ? ");
432 if (gmx_ask_yesno(bASK
)) {
433 /* Make swap sidechain C=O index */
436 for(n
=0;((atoms
->atom
[n
].resind
< npres
) && (n
<atoms
->nr
));) {
437 resind
= atoms
->atom
[n
].resind
;
439 for(;((atoms
->atom
[n
].resind
==resind
) && (n
<atoms
->nr
));n
++)
440 if (strcmp("CA",*atoms
->atomname
[n
]) == 0) {
444 } else if (strcmp("C",*atoms
->atomname
[n
]) == 0) {
446 gmx_incons("Atom naming problem");
448 } else if (strcmp("O",*atoms
->atomname
[n
]) == 0) {
450 gmx_incons("Atom naming problem");
452 } else if (strcmp("O1",*atoms
->atomname
[n
]) == 0) {
454 gmx_incons("Atom naming problem");
459 /* copy the residuename to the tail of the groupname */
461 add_grp(gb
,gn
,nra
,aid
,"SwapSC-CO");
472 /* Return 0 if the name was found, otherwise -1.
473 * p_restype is set to a pointer to the type name, or 'Other' if we did not find it.
476 gmx_residuetype_get_type(gmx_residuetype_t rt
,const char * resname
, const char ** p_restype
)
481 for(i
=0;i
<rt
->n
&& rc
;i
++)
483 rc
=gmx_strcasecmp(rt
->resname
[i
],resname
);
486 *p_restype
= (rc
==0) ? rt
->restype
[i
-1] : gmx_residuetype_undefined
;
492 gmx_residuetype_add(gmx_residuetype_t rt
,const char *newresname
, const char *newrestype
)
496 const char * p_oldtype
;
498 found
= !gmx_residuetype_get_type(rt
,newresname
,&p_oldtype
);
500 if(found
&& gmx_strcasecmp(p_oldtype
,newrestype
))
502 fprintf(stderr
,"Warning: Residue '%s' already present with type '%s' in database, ignoring new type '%s'.",
503 newresname
,p_oldtype
,newrestype
);
508 srenew(rt
->resname
,rt
->n
+1);
509 srenew(rt
->restype
,rt
->n
+1);
510 rt
->resname
[rt
->n
]=strdup(newresname
);
511 rt
->restype
[rt
->n
]=strdup(newrestype
);
520 gmx_residuetype_init(gmx_residuetype_t
*prt
)
524 char resname
[STRLEN
],restype
[STRLEN
],dum
[STRLEN
];
527 struct gmx_residuetype
*rt
;
536 db
=libopen("residuetypes.dat");
538 while(get_a_line(db
,line
,STRLEN
))
544 if(sscanf(line
,"%s %s %s",resname
,restype
,dum
)!=2)
546 gmx_fatal(FARGS
,"Incorrect number of columns (2 expected) for line in residuetypes.dat");
548 gmx_residuetype_add(rt
,resname
,restype
);
560 gmx_residuetype_destroy(gmx_residuetype_t rt
)
566 free(rt
->resname
[i
]);
567 free(rt
->restype
[i
]);
575 gmx_residuetype_get_alltypes(gmx_residuetype_t rt
,
576 const char *** p_typenames
,
581 const char ** my_typename
;
591 for(j
=0;j
<n
&& !found
;j
++)
593 found
=!gmx_strcasecmp(p
,my_typename
[j
]);
598 srenew(my_typename
,n
+1);
604 *p_typenames
=my_typename
;
612 gmx_residuetype_is_protein(gmx_residuetype_t rt
, const char *resnm
)
617 if(gmx_residuetype_get_type(rt
,resnm
,&p_type
)==0 &&
618 gmx_strcasecmp(p_type
,"Protein")==0)
630 gmx_residuetype_is_dna(gmx_residuetype_t rt
, const char *resnm
)
635 if(gmx_residuetype_get_type(rt
,resnm
,&p_type
)==0 &&
636 gmx_strcasecmp(p_type
,"DNA")==0)
648 gmx_residuetype_is_rna(gmx_residuetype_t rt
, const char *resnm
)
653 if(gmx_residuetype_get_type(rt
,resnm
,&p_type
)==0 &&
654 gmx_strcasecmp(p_type
,"RNA")==0)
668 void analyse(t_atoms
*atoms
,t_blocka
*gb
,char ***gn
,gmx_bool bASK
,gmx_bool bVerb
)
670 gmx_residuetype_t rt
;
673 const char ** restype
;
679 const char ** p_typename
;
686 printf("Analysing residue names:\n");
688 /* Create system group, every single atom */
690 for(i
=0;i
<atoms
->nr
;i
++)
694 add_grp(gb
,gn
,atoms
->nr
,aid
,"System");
697 /* For every residue, get a pointer to the residue type name */
698 gmx_residuetype_init(&rt
);
700 snew(restype
,atoms
->nres
);
703 for(i
=0;i
<atoms
->nres
;i
++)
705 resnm
= *atoms
->resinfo
[i
].name
;
706 gmx_residuetype_get_type(rt
,resnm
,&(restype
[i
]));
711 p_typename
[ntypes
++] = strdup(restype
[i
]);
715 /* Note that this does not lead to a N*N loop, but N*K, where
716 * K is the number of residue _types_, which is small and independent of N.
719 for(k
=0;k
<i
&& !found
;k
++)
721 found
= !strcmp(restype
[i
],restype
[k
]);
725 srenew(p_typename
,ntypes
+1);
726 p_typename
[ntypes
++] = strdup(restype
[i
]);
731 p_status(restype
,atoms
->nres
,p_typename
,ntypes
,bVerb
);
733 for(k
=0;k
<ntypes
;k
++)
735 aid
=mk_aid(atoms
,restype
,p_typename
[k
],&nra
,TRUE
);
737 /* Check for special types to do fancy stuff with */
739 if(!gmx_strcasecmp(p_typename
[k
],"Protein") && nra
>0)
743 analyse_prot(restype
,atoms
,gb
,gn
,bASK
,bVerb
);
745 /* Create a Non-Protein group */
746 aid
=mk_aid(atoms
,restype
,"Protein",&nra
,FALSE
);
747 if ((nra
> 0) && (nra
< atoms
->nr
))
749 add_grp(gb
,gn
,nra
,aid
,"non-Protein");
753 else if(!gmx_strcasecmp(p_typename
[k
],"Water") && nra
>0)
755 add_grp(gb
,gn
,nra
,aid
,p_typename
[k
]);
756 /* Add this group as 'SOL' too, for backward compatibility with older gromacs versions */
757 add_grp(gb
,gn
,nra
,aid
,"SOL");
761 /* Solvent, create a negated group too */
762 aid
=mk_aid(atoms
,restype
,"Water",&nra
,FALSE
);
763 if ((nra
> 0) && (nra
< atoms
->nr
))
765 add_grp(gb
,gn
,nra
,aid
,"non-Water");
772 add_grp(gb
,gn
,nra
,aid
,p_typename
[k
]);
774 analyse_other(restype
,atoms
,gb
,gn
,bASK
,bVerb
);
780 gmx_residuetype_destroy(rt
);
782 /* Create a merged water_and_ions group */
788 for(i
=0;i
<gb
->nr
;i
++)
790 if(!gmx_strcasecmp((*gn
)[i
],"Water"))
793 nwater
= gb
->index
[i
+1]-gb
->index
[i
];
795 else if(!gmx_strcasecmp((*gn
)[i
],"Ion"))
798 nion
= gb
->index
[i
+1]-gb
->index
[i
];
802 if(nwater
>0 && nion
>0)
804 srenew(gb
->index
,gb
->nr
+2);
805 srenew(*gn
,gb
->nr
+1);
806 (*gn
)[gb
->nr
] = strdup("Water_and_ions");
807 srenew(gb
->a
,gb
->nra
+nwater
+nion
);
810 for(i
=gb
->index
[iwater
];i
<gb
->index
[iwater
+1];i
++)
812 gb
->a
[gb
->nra
++] = gb
->a
[i
];
817 for(i
=gb
->index
[iion
];i
<gb
->index
[iion
+1];i
++)
819 gb
->a
[gb
->nra
++] = gb
->a
[i
];
823 gb
->index
[gb
->nr
]=gb
->nra
;
828 void check_index(char *gname
,int n
,atom_id index
[],char *traj
,int natoms
)
833 if (index
[i
] >= natoms
)
834 gmx_fatal(FARGS
,"%s atom number (index[%d]=%d) is larger than the number of atoms in %s (%d)",
835 gname
? gname
: "Index",i
+1, index
[i
]+1,
836 traj
? traj
: "the trajectory",natoms
);
837 else if (index
[i
] < 0)
838 gmx_fatal(FARGS
,"%s atom number (index[%d]=%d) is less than zero",
839 gname
? gname
: "Index",i
+1, index
[i
]+1);
842 t_blocka
*init_index(const char *gfile
, char ***grpname
)
848 char line
[STRLEN
],*pt
,str
[STRLEN
];
850 in
=gmx_fio_fopen(gfile
,"r");
852 get_a_line(in
,line
,STRLEN
);
853 if ( line
[0]=='[' ) {
862 if (get_header(line
,str
)) {
864 srenew(b
->index
,b
->nr
+1);
865 srenew(*grpname
,b
->nr
);
868 b
->index
[b
->nr
]=b
->index
[b
->nr
-1];
869 (*grpname
)[b
->nr
-1]=strdup(str
);
872 while ((i
=sscanf(pt
,"%s",str
)) == 1) {
876 srenew(b
->a
,maxentries
);
878 b
->a
[i
]=strtol(str
, NULL
, 10)-1;
881 pt
=strstr(pt
,str
)+strlen(str
);
884 } while (get_a_line(in
,line
,STRLEN
));
888 sscanf(line
,"%d%d",&b
->nr
,&b
->nra
);
889 snew(b
->index
,b
->nr
+1);
890 snew(*grpname
,b
->nr
);
893 for (i
=0; (i
<b
->nr
); i
++) {
894 nread
=fscanf(in
,"%s%d",str
,&ng
);
895 (*grpname
)[i
]=strdup(str
);
896 b
->index
[i
+1]=b
->index
[i
]+ng
;
897 if (b
->index
[i
+1] > b
->nra
)
898 gmx_fatal(FARGS
,"Something wrong in your indexfile at group %s",str
);
899 for(j
=0; (j
<ng
); j
++) {
900 nread
=fscanf(in
,"%d",&a
);
901 b
->a
[b
->index
[i
]+j
]=a
;
907 for(i
=0; (i
<b
->nr
); i
++) {
908 for(j
=b
->index
[i
]; (j
<b
->index
[i
+1]); j
++) {
910 fprintf(stderr
,"\nWARNING: negative index %d in group %s\n\n",
911 b
->a
[j
],(*grpname
)[i
]);
918 static void minstring(char *str
)
922 for (i
=0; (i
< (int)strlen(str
)); i
++)
927 int find_group(char s
[], int ngrps
, char **grpname
)
936 /* first look for whole name match */
938 for(i
=0; i
<ngrps
; i
++)
939 if (gmx_strcasecmp_min(s
,grpname
[i
])==0) {
944 /* second look for first string match */
946 for(i
=0; i
<ngrps
; i
++)
947 if (gmx_strncasecmp_min(s
,grpname
[i
],n
)==0) {
952 /* last look for arbitrary substring match */
956 for(i
=0; i
<ngrps
; i
++) {
957 strcpy(string
, grpname
[i
]);
960 if (strstr(string
,s
)!=NULL
) {
968 printf("Error: Multiple groups '%s' selected\n", s
);
974 static int qgroup(int *a
, int ngrps
, char **grpname
)
982 fprintf(stderr
,"Select a group: ");
984 if ( scanf("%s",s
)!=1 )
985 gmx_fatal(FARGS
,"Cannot read from input");
986 trim(s
); /* remove spaces */
987 } while (strlen(s
)==0);
988 aa
= strtol(s
, &end
, 10);
989 if (aa
==0 && end
[0] != '\0') /* string entered */
990 aa
= find_group(s
, ngrps
, grpname
);
991 bInRange
= (aa
>= 0 && aa
< ngrps
);
993 printf("Error: No such group '%s'\n", s
);
995 printf("Selected %d: '%s'\n", aa
, grpname
[aa
]);
1000 static void rd_groups(t_blocka
*grps
,char **grpname
,char *gnames
[],
1001 int ngrps
,int isize
[],atom_id
*index
[],int grpnr
[])
1006 gmx_fatal(FARGS
,"Error: no groups in indexfile");
1007 for(i
=0; (i
<grps
->nr
); i
++)
1008 fprintf(stderr
,"Group %5d (%15s) has %5d elements\n",i
,grpname
[i
],
1009 grps
->index
[i
+1]-grps
->index
[i
]);
1010 for(i
=0; (i
<ngrps
); i
++) {
1013 gnr1
=qgroup(&grpnr
[i
], grps
->nr
, grpname
);
1014 if ((gnr1
<0) || (gnr1
>=grps
->nr
))
1015 fprintf(stderr
,"Select between %d and %d.\n",0,grps
->nr
-1);
1016 } while ((gnr1
<0) || (gnr1
>=grps
->nr
));
1018 fprintf(stderr
,"There is one group in the index\n");
1021 gnames
[i
]=strdup(grpname
[gnr1
]);
1022 isize
[i
]=grps
->index
[gnr1
+1]-grps
->index
[gnr1
];
1023 snew(index
[i
],isize
[i
]);
1024 for(j
=0; (j
<isize
[i
]); j
++)
1025 index
[i
][j
]=grps
->a
[grps
->index
[gnr1
]+j
];
1029 void rd_index(const char *statfile
,int ngrps
,int isize
[],
1030 atom_id
*index
[],char *grpnames
[])
1038 gmx_fatal(FARGS
,"No index file specified");
1039 grps
=init_index(statfile
,&gnames
);
1040 rd_groups(grps
,gnames
,grpnames
,ngrps
,isize
,index
,grpnr
);
1043 void rd_index_nrs(char *statfile
,int ngrps
,int isize
[],
1044 atom_id
*index
[],char *grpnames
[],int grpnr
[])
1050 gmx_fatal(FARGS
,"No index file specified");
1051 grps
=init_index(statfile
,&gnames
);
1053 rd_groups(grps
,gnames
,grpnames
,ngrps
,isize
,index
,grpnr
);
1056 void get_index(t_atoms
*atoms
, const char *fnm
, int ngrps
,
1057 int isize
[], atom_id
*index
[],char *grpnames
[])
1060 t_blocka
*grps
= NULL
;
1066 grps
=init_index(fnm
,gnames
);
1070 snew(grps
->index
,1);
1071 analyse(atoms
,grps
,gnames
,FALSE
,FALSE
);
1074 gmx_incons("You need to supply a valid atoms structure or a valid index file name");
1076 rd_groups(grps
,*gnames
,grpnames
,ngrps
,isize
,index
,grpnr
);
1079 t_cluster_ndx
*cluster_index(FILE *fplog
,const char *ndx
)
1085 c
->clust
= init_index(ndx
,&c
->grpname
);
1087 for(i
=0; (i
<c
->clust
->nra
); i
++)
1088 c
->maxframe
= max(c
->maxframe
,c
->clust
->a
[i
]);
1089 fprintf(fplog
? fplog
: stdout
,
1090 "There are %d clusters containing %d structures, highest framenr is %d\n",
1091 c
->clust
->nr
,c
->clust
->nra
,c
->maxframe
);
1093 pr_blocka(debug
,0,"clust",c
->clust
,TRUE
);
1094 for(i
=0; (i
<c
->clust
->nra
); i
++)
1095 if ((c
->clust
->a
[i
] < 0) || (c
->clust
->a
[i
] > c
->maxframe
))
1096 gmx_fatal(FARGS
,"Range check error for c->clust->a[%d] = %d\n"
1097 "should be within 0 and %d",i
,c
->clust
->a
[i
],c
->maxframe
+1);
1099 c
->inv_clust
=make_invblocka(c
->clust
,c
->maxframe
);