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 * check out http://www.gromacs.org for more information.
7 * Copyright (c) 2012, by the GROMACS development team, led by
8 * David van der Spoel, Berk Hess, Erik Lindahl, and including many
9 * others, as listed in the AUTHORS file in the top-level source
10 * directory and at http://www.gromacs.org.
12 * GROMACS is free software; you can redistribute it and/or
13 * modify it under the terms of the GNU Lesser General Public License
14 * as published by the Free Software Foundation; either version 2.1
15 * of the License, or (at your option) any later version.
17 * GROMACS is distributed in the hope that it will be useful,
18 * but WITHOUT ANY WARRANTY; without even the implied warranty of
19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 * Lesser General Public License for more details.
22 * You should have received a copy of the GNU Lesser General Public
23 * License along with GROMACS; if not, see
24 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
25 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
27 * If you want to redistribute modifications to GROMACS, please
28 * consider that scientific software is very special. Version
29 * control is crucial - bugs must be traceable. We will be happy to
30 * consider code for inclusion in the official distribution, but
31 * derived work must not be called official GROMACS. Details are found
32 * in the README & COPYING files - if they are missing, get the
33 * official version at http://www.gromacs.org.
35 * To help us fund GROMACS development, we humbly ask that you cite
36 * the research papers on the package. Check out http://www.gromacs.org.
65 const char gmx_residuetype_undefined
[]="Other";
67 struct gmx_residuetype
76 static gmx_bool
gmx_ask_yesno(gmx_bool bASK
)
82 c
=toupper(fgetc(stdin
));
83 } while ((c
!= 'Y') && (c
!= 'N'));
91 t_blocka
*new_blocka(void)
101 void write_index(const char *outf
, t_blocka
*b
,char **gnames
)
106 out
=gmx_fio_fopen(outf
,"w");
107 /* fprintf(out,"%5d %5d\n",b->nr,b->nra); */
108 for(i
=0; (i
<b
->nr
); i
++) {
109 fprintf(out
,"[ %s ]\n",gnames
[i
]);
110 for(k
=0,j
=b
->index
[i
]; j
<b
->index
[i
+1]; j
++,k
++) {
111 fprintf(out
,"%4d ",b
->a
[j
]+1);
120 void add_grp(t_blocka
*b
,char ***gnames
,int nra
,atom_id a
[],const char *name
)
124 srenew(b
->index
,b
->nr
+2);
125 srenew(*gnames
,b
->nr
+1);
126 (*gnames
)[b
->nr
]=strdup(name
);
128 srenew(b
->a
,b
->nra
+nra
);
129 for(i
=0; (i
<nra
); i
++)
132 b
->index
[b
->nr
]=b
->nra
;
135 /* compare index in `a' with group in `b' at `index',
136 when `index'<0 it is relative to end of `b' */
137 static gmx_bool
grp_cmp(t_blocka
*b
, int nra
, atom_id a
[], int index
)
142 index
= b
->nr
-1+index
;
144 gmx_fatal(FARGS
,"no such index group %d in t_blocka (nr=%d)",index
,b
->nr
);
146 if ( nra
!= b
->index
[index
+1] - b
->index
[index
] )
149 if ( a
[i
] != b
->a
[b
->index
[index
]+i
] )
155 p_status(const char **restype
, int nres
, const char **typenames
, int ntypes
)
162 snew(counter
,ntypes
);
163 for(i
=0;i
<ntypes
;i
++)
170 for(j
=0;j
<ntypes
;j
++)
172 if(!gmx_strcasecmp(restype
[i
],typenames
[j
]))
179 for(i
=0; (i
<ntypes
); i
++)
183 printf("There are: %5d %10s residues\n",counter
[i
],typenames
[i
]);
192 mk_aid(t_atoms
*atoms
,const char ** restype
,const char * typestring
,int *nra
,gmx_bool bMatch
)
193 /* Make an array of atom_ids for all atoms with residuetypes matching typestring, or the opposite if bMatch is false */
201 for(i
=0; (i
<atoms
->nr
); i
++)
203 res
=!gmx_strcasecmp(restype
[atoms
->atom
[i
].resind
],typestring
);
223 static void analyse_other(const char ** restype
,t_atoms
*atoms
,
224 t_blocka
*gb
,char ***gn
,gmx_bool bASK
,gmx_bool bVerb
)
229 atom_id
*other_ndx
,*aid
,*aaid
;
230 int i
,j
,k
,l
,resind
,naid
,naaid
,natp
,nrestp
=0;
232 for(i
=0; (i
<atoms
->nres
); i
++)
234 if (gmx_strcasecmp(restype
[i
],"Protein") && gmx_strcasecmp(restype
[i
],"DNA") && gmx_strcasecmp(restype
[i
],"RNA") && gmx_strcasecmp(restype
[i
],"Water"))
239 if (i
< atoms
->nres
) {
242 printf("Analysing residues not classified as Protein/DNA/RNA/Water and splitting into groups...\n");
243 snew(other_ndx
,atoms
->nr
);
244 for(k
=0; (k
<atoms
->nr
); k
++) {
245 resind
= atoms
->atom
[k
].resind
;
246 rname
= *atoms
->resinfo
[resind
].name
;
247 if (gmx_strcasecmp(restype
[resind
],"Protein") && gmx_strcasecmp(restype
[resind
],"DNA") &&
248 gmx_strcasecmp(restype
[resind
],"RNA") && gmx_strcasecmp(restype
[resind
],"Water"))
251 for(l
=0; (l
<nrestp
); l
++)
252 if (strcmp(restp
[l
].rname
,rname
) == 0)
255 srenew(restp
,nrestp
+1);
256 restp
[nrestp
].rname
= strdup(rname
);
257 restp
[nrestp
].bNeg
= FALSE
;
258 restp
[nrestp
].gname
= strdup(rname
);
264 for(i
=0; (i
<nrestp
); i
++) {
267 for(j
=0; (j
<atoms
->nr
); j
++) {
268 rname
= *atoms
->resinfo
[atoms
->atom
[j
].resind
].name
;
269 if ((strcmp(restp
[i
].rname
,rname
) == 0 && !restp
[i
].bNeg
) ||
270 (strcmp(restp
[i
].rname
,rname
) != 0 && restp
[i
].bNeg
)) {
274 add_grp(gb
,gn
,naid
,aid
,restp
[i
].gname
);
276 printf("split %s into atoms (y/n) ? ",restp
[i
].gname
);
278 if (gmx_ask_yesno(bASK
)) {
280 for(k
=0; (k
<naid
); k
++) {
281 aname
=*atoms
->atomname
[aid
[k
]];
282 for(l
=0; (l
<natp
); l
++)
283 if (strcmp(aname
,attp
[l
]) == 0)
291 for(l
=0; (l
<natp
); l
++) {
294 for(k
=0; (k
<naid
); k
++) {
295 aname
=*atoms
->atomname
[aid
[k
]];
296 if (strcmp(aname
,attp
[l
])==0)
297 aaid
[naaid
++]=aid
[k
];
299 add_grp(gb
,gn
,naaid
,aaid
,attp
[l
]);
313 /*! /brief Instances of this struct contain the data necessary to
314 * construct a single (protein) index group in
316 typedef struct gmx_help_make_index_group
318 /* The set of atom names that will be used to form this index group */
319 const char **defining_atomnames
;
320 /* Size of the defining_atomnames array */
321 const int num_defining_atomnames
;
322 /* Name of this index group */
323 const char *group_name
;
324 /* Whether the above atom names name the atoms in the group, or
325 * those not in the group */
326 gmx_bool bTakeComplement
;
327 /* The index in wholename gives the first item in the arrays of
328 * atomnames that should be tested with 'gmx_strncasecmp' in stead of
329 * gmx_strcasecmp, or -1 if all items should be tested with strcasecmp
330 * This is comparable to using a '*' wildcard at the end of specific
331 * atom names, but that is more involved to implement...
334 /* Only create this index group if it differs from the one specified in compareto,
335 where -1 means to always create this group. */
337 } t_gmx_help_make_index_group
;
339 static void analyse_prot(const char ** restype
,t_atoms
*atoms
,
340 t_blocka
*gb
,char ***gn
,gmx_bool bASK
,gmx_bool bVerb
)
342 /* lists of atomnames to be used in constructing index groups: */
343 static const char *pnoh
[] = { "H", "HN" };
344 static const char *pnodum
[] = { "MN1", "MN2", "MCB1", "MCB2", "MCG1", "MCG2",
345 "MCD1", "MCD2", "MCE1", "MCE2", "MNZ1", "MNZ2" };
346 static const char *calpha
[] = { "CA" };
347 static const char *bb
[] = { "N","CA","C" };
348 static const char *mc
[] = { "N","CA","C","O","O1","O2","OC1","OC2","OT","OXT" };
349 static const char *mcb
[] = { "N","CA","CB","C","O","O1","O2","OC1","OC2","OT","OXT" };
350 static const char *mch
[] = { "N","CA","C","O","O1","O2","OC1","OC2","OT","OXT",
351 "H1","H2","H3","H","HN" };
353 static const t_gmx_help_make_index_group constructing_data
[] =
354 {{ NULL
, 0, "Protein", TRUE
, -1, -1},
355 { pnoh
, asize(pnoh
), "Protein-H", TRUE
, 0, -1},
356 { calpha
, asize(calpha
), "C-alpha", FALSE
, -1, -1},
357 { bb
, asize(bb
), "Backbone", FALSE
, -1, -1},
358 { mc
, asize(mc
), "MainChain", FALSE
, -1, -1},
359 { mcb
, asize(mcb
), "MainChain+Cb", FALSE
, -1, -1},
360 { mch
, asize(mch
), "MainChain+H", FALSE
, -1, -1},
361 { mch
, asize(mch
), "SideChain", TRUE
, -1, -1},
362 { mch
, asize(mch
), "SideChain-H", TRUE
, 11, -1},
363 { pnodum
, asize(pnodum
), "Prot-Masses", TRUE
, -1, 0},
365 const int num_index_groups
= asize(constructing_data
);
369 int nra
,nnpres
,npres
;
371 char ndx_name
[STRLEN
],*atnm
;
376 printf("Analysing Protein...\n");
380 /* calculate the number of protein residues */
382 for(i
=0; (i
<atoms
->nres
); i
++) {
383 if (0 == gmx_strcasecmp(restype
[i
],"Protein")) {
387 /* find matching or complement atoms */
388 for(i
=0; (i
<(int)num_index_groups
); i
++) {
390 for(n
=0; (n
<atoms
->nr
); n
++) {
391 if (0 == gmx_strcasecmp(restype
[atoms
->atom
[n
].resind
],"Protein")) {
393 for(j
=0; (j
<constructing_data
[i
].num_defining_atomnames
); j
++) {
394 /* skip digits at beginning of atomname, e.g. 1H */
395 atnm
=*atoms
->atomname
[n
];
396 while (isdigit(atnm
[0])) {
399 if ( (constructing_data
[i
].wholename
==-1) || (j
<constructing_data
[i
].wholename
) ) {
400 if (0 == gmx_strcasecmp(constructing_data
[i
].defining_atomnames
[j
],atnm
)) {
404 if (0 == gmx_strncasecmp(constructing_data
[i
].defining_atomnames
[j
],atnm
,strlen(constructing_data
[i
].defining_atomnames
[j
]))) {
409 if (constructing_data
[i
].bTakeComplement
!= match
) {
414 /* if we want to add this group always or it differs from previous
416 if ( -1 == constructing_data
[i
].compareto
|| !grp_cmp(gb
,nra
,aid
,constructing_data
[i
].compareto
-i
) ) {
417 add_grp(gb
,gn
,nra
,aid
,constructing_data
[i
].group_name
);
422 for(i
=0; (i
<(int)num_index_groups
); i
++) {
423 printf("Split %12s into %5d residues (y/n) ? ",constructing_data
[i
].group_name
,npres
);
424 if (gmx_ask_yesno(bASK
)) {
427 for(n
=0;((atoms
->atom
[n
].resind
< npres
) && (n
<atoms
->nr
));) {
428 resind
= atoms
->atom
[n
].resind
;
429 for(;((atoms
->atom
[n
].resind
==resind
) && (n
<atoms
->nr
));n
++) {
431 for(j
=0;(j
<constructing_data
[i
].num_defining_atomnames
); j
++) {
432 if (0 == gmx_strcasecmp(constructing_data
[i
].defining_atomnames
[j
],*atoms
->atomname
[n
])) {
436 if (constructing_data
[i
].bTakeComplement
!= match
) {
440 /* copy the residuename to the tail of the groupname */
443 ri
= &atoms
->resinfo
[resind
];
444 sprintf(ndx_name
,"%s_%s%d%c",
445 constructing_data
[i
].group_name
,*ri
->name
,ri
->nr
,ri
->ic
==' ' ? '\0' : ri
->ic
);
446 add_grp(gb
,gn
,nra
,aid
,ndx_name
);
452 printf("Make group with sidechain and C=O swapped (y/n) ? ");
453 if (gmx_ask_yesno(bASK
)) {
454 /* Make swap sidechain C=O index */
457 for(n
=0;((atoms
->atom
[n
].resind
< npres
) && (n
<atoms
->nr
));) {
458 resind
= atoms
->atom
[n
].resind
;
460 for(;((atoms
->atom
[n
].resind
==resind
) && (n
<atoms
->nr
));n
++)
461 if (strcmp("CA",*atoms
->atomname
[n
]) == 0) {
465 } else if (strcmp("C",*atoms
->atomname
[n
]) == 0) {
467 gmx_incons("Atom naming problem");
470 } else if (strcmp("O",*atoms
->atomname
[n
]) == 0) {
472 gmx_incons("Atom naming problem");
475 } else if (strcmp("O1",*atoms
->atomname
[n
]) == 0) {
477 gmx_incons("Atom naming problem");
483 /* copy the residuename to the tail of the groupname */
485 add_grp(gb
,gn
,nra
,aid
,"SwapSC-CO");
496 /* Return 0 if the name was found, otherwise -1.
497 * p_restype is set to a pointer to the type name, or 'Other' if we did not find it.
500 gmx_residuetype_get_type(gmx_residuetype_t rt
,const char * resname
, const char ** p_restype
)
505 for(i
=0;i
<rt
->n
&& rc
;i
++)
507 rc
=gmx_strcasecmp(rt
->resname
[i
],resname
);
510 *p_restype
= (rc
==0) ? rt
->restype
[i
-1] : gmx_residuetype_undefined
;
516 gmx_residuetype_add(gmx_residuetype_t rt
,const char *newresname
, const char *newrestype
)
520 const char * p_oldtype
;
522 found
= !gmx_residuetype_get_type(rt
,newresname
,&p_oldtype
);
524 if(found
&& gmx_strcasecmp(p_oldtype
,newrestype
))
526 fprintf(stderr
,"Warning: Residue '%s' already present with type '%s' in database, ignoring new type '%s'.",
527 newresname
,p_oldtype
,newrestype
);
532 srenew(rt
->resname
,rt
->n
+1);
533 srenew(rt
->restype
,rt
->n
+1);
534 rt
->resname
[rt
->n
]=strdup(newresname
);
535 rt
->restype
[rt
->n
]=strdup(newrestype
);
544 gmx_residuetype_init(gmx_residuetype_t
*prt
)
548 char resname
[STRLEN
],restype
[STRLEN
],dum
[STRLEN
];
551 struct gmx_residuetype
*rt
;
560 db
=libopen("residuetypes.dat");
562 while(get_a_line(db
,line
,STRLEN
))
568 if(sscanf(line
,"%s %s %s",resname
,restype
,dum
)!=2)
570 gmx_fatal(FARGS
,"Incorrect number of columns (2 expected) for line in residuetypes.dat");
572 gmx_residuetype_add(rt
,resname
,restype
);
584 gmx_residuetype_destroy(gmx_residuetype_t rt
)
590 sfree(rt
->resname
[i
]);
591 sfree(rt
->restype
[i
]);
601 gmx_residuetype_get_alltypes(gmx_residuetype_t rt
,
602 const char *** p_typenames
,
607 const char ** my_typename
;
617 for(j
=0;j
<n
&& !found
;j
++)
619 found
=!gmx_strcasecmp(p
,my_typename
[j
]);
624 srenew(my_typename
,n
+1);
630 *p_typenames
=my_typename
;
638 gmx_residuetype_is_protein(gmx_residuetype_t rt
, const char *resnm
)
643 if(gmx_residuetype_get_type(rt
,resnm
,&p_type
)==0 &&
644 gmx_strcasecmp(p_type
,"Protein")==0)
656 gmx_residuetype_is_dna(gmx_residuetype_t rt
, const char *resnm
)
661 if(gmx_residuetype_get_type(rt
,resnm
,&p_type
)==0 &&
662 gmx_strcasecmp(p_type
,"DNA")==0)
674 gmx_residuetype_is_rna(gmx_residuetype_t rt
, const char *resnm
)
679 if(gmx_residuetype_get_type(rt
,resnm
,&p_type
)==0 &&
680 gmx_strcasecmp(p_type
,"RNA")==0)
691 /* Return the size of the arrays */
693 gmx_residuetype_get_size(gmx_residuetype_t rt
)
698 /* Search for a residuetype with name resnm within the
699 * gmx_residuetype database. Return the index if found,
703 gmx_residuetype_get_index(gmx_residuetype_t rt
, const char *resnm
)
708 for(i
=0;i
<rt
->n
&& rc
;i
++)
710 rc
=gmx_strcasecmp(rt
->resname
[i
],resnm
);
713 return (0 == rc
) ? i
-1 : -1;
716 /* Return the name of the residuetype with the given index, or
717 * NULL if not found. */
719 gmx_residuetype_get_name(gmx_residuetype_t rt
, int index
)
721 if(index
>= 0 && index
< rt
->n
) {
722 return rt
->resname
[index
];
730 void analyse(t_atoms
*atoms
,t_blocka
*gb
,char ***gn
,gmx_bool bASK
,gmx_bool bVerb
)
732 gmx_residuetype_t rt
=NULL
;
735 const char ** restype
;
741 const char ** p_typename
;
748 printf("Analysing residue names:\n");
750 /* Create system group, every single atom */
752 for(i
=0;i
<atoms
->nr
;i
++)
756 add_grp(gb
,gn
,atoms
->nr
,aid
,"System");
759 /* For every residue, get a pointer to the residue type name */
760 gmx_residuetype_init(&rt
);
763 snew(restype
,atoms
->nres
);
766 for(i
=0;i
<atoms
->nres
;i
++)
768 resnm
= *atoms
->resinfo
[i
].name
;
769 gmx_residuetype_get_type(rt
,resnm
,&(restype
[i
]));
771 /* Note that this does not lead to a N*N loop, but N*K, where
772 * K is the number of residue _types_, which is small and independent of N.
775 for(k
=0;k
<ntypes
&& !found
;k
++)
777 found
= !strcmp(restype
[i
],p_typename
[k
]);
781 srenew(p_typename
,ntypes
+1);
782 p_typename
[ntypes
++] = strdup(restype
[i
]);
788 p_status(restype
,atoms
->nres
,p_typename
,ntypes
);
791 for(k
=0;k
<ntypes
;k
++)
793 aid
=mk_aid(atoms
,restype
,p_typename
[k
],&nra
,TRUE
);
795 /* Check for special types to do fancy stuff with */
797 if(!gmx_strcasecmp(p_typename
[k
],"Protein") && nra
>0)
801 analyse_prot(restype
,atoms
,gb
,gn
,bASK
,bVerb
);
803 /* Create a Non-Protein group */
804 aid
=mk_aid(atoms
,restype
,"Protein",&nra
,FALSE
);
805 if ((nra
> 0) && (nra
< atoms
->nr
))
807 add_grp(gb
,gn
,nra
,aid
,"non-Protein");
811 else if(!gmx_strcasecmp(p_typename
[k
],"Water") && nra
>0)
813 add_grp(gb
,gn
,nra
,aid
,p_typename
[k
]);
814 /* Add this group as 'SOL' too, for backward compatibility with older gromacs versions */
815 add_grp(gb
,gn
,nra
,aid
,"SOL");
819 /* Solvent, create a negated group too */
820 aid
=mk_aid(atoms
,restype
,"Water",&nra
,FALSE
);
821 if ((nra
> 0) && (nra
< atoms
->nr
))
823 add_grp(gb
,gn
,nra
,aid
,"non-Water");
830 add_grp(gb
,gn
,nra
,aid
,p_typename
[k
]);
832 analyse_other(restype
,atoms
,gb
,gn
,bASK
,bVerb
);
838 gmx_residuetype_destroy(rt
);
840 /* Create a merged water_and_ions group */
846 for(i
=0;i
<gb
->nr
;i
++)
848 if(!gmx_strcasecmp((*gn
)[i
],"Water"))
851 nwater
= gb
->index
[i
+1]-gb
->index
[i
];
853 else if(!gmx_strcasecmp((*gn
)[i
],"Ion"))
856 nion
= gb
->index
[i
+1]-gb
->index
[i
];
860 if(nwater
>0 && nion
>0)
862 srenew(gb
->index
,gb
->nr
+2);
863 srenew(*gn
,gb
->nr
+1);
864 (*gn
)[gb
->nr
] = strdup("Water_and_ions");
865 srenew(gb
->a
,gb
->nra
+nwater
+nion
);
868 for(i
=gb
->index
[iwater
];i
<gb
->index
[iwater
+1];i
++)
870 gb
->a
[gb
->nra
++] = gb
->a
[i
];
875 for(i
=gb
->index
[iion
];i
<gb
->index
[iion
+1];i
++)
877 gb
->a
[gb
->nra
++] = gb
->a
[i
];
881 gb
->index
[gb
->nr
]=gb
->nra
;
886 void check_index(char *gname
,int n
,atom_id index
[],char *traj
,int natoms
)
891 if (index
[i
] >= natoms
)
892 gmx_fatal(FARGS
,"%s atom number (index[%d]=%d) is larger than the number of atoms in %s (%d)",
893 gname
? gname
: "Index",i
+1, index
[i
]+1,
894 traj
? traj
: "the trajectory",natoms
);
895 else if (index
[i
] < 0)
896 gmx_fatal(FARGS
,"%s atom number (index[%d]=%d) is less than zero",
897 gname
? gname
: "Index",i
+1, index
[i
]+1);
900 t_blocka
*init_index(const char *gfile
, char ***grpname
)
906 char line
[STRLEN
],*pt
,str
[STRLEN
];
908 in
=gmx_fio_fopen(gfile
,"r");
910 get_a_line(in
,line
,STRLEN
);
911 if ( line
[0]=='[' ) {
920 if (get_header(line
,str
)) {
922 srenew(b
->index
,b
->nr
+1);
923 srenew(*grpname
,b
->nr
);
926 b
->index
[b
->nr
]=b
->index
[b
->nr
-1];
927 (*grpname
)[b
->nr
-1]=strdup(str
);
931 gmx_fatal(FARGS
,"The first header of your indexfile is invalid");
934 while (sscanf(pt
,"%s",str
) == 1) {
938 srenew(b
->a
,maxentries
);
940 b
->a
[i
]=strtol(str
, NULL
, 10)-1;
943 pt
=strstr(pt
,str
)+strlen(str
);
946 } while (get_a_line(in
,line
,STRLEN
));
950 sscanf(line
,"%d%d",&b
->nr
,&b
->nra
);
951 snew(b
->index
,b
->nr
+1);
952 snew(*grpname
,b
->nr
);
955 for (i
=0; (i
<b
->nr
); i
++) {
956 nread
=fscanf(in
,"%s%d",str
,&ng
);
957 (*grpname
)[i
]=strdup(str
);
958 b
->index
[i
+1]=b
->index
[i
]+ng
;
959 if (b
->index
[i
+1] > b
->nra
)
960 gmx_fatal(FARGS
,"Something wrong in your indexfile at group %s",str
);
961 for(j
=0; (j
<ng
); j
++) {
962 nread
=fscanf(in
,"%d",&a
);
963 b
->a
[b
->index
[i
]+j
]=a
;
969 for(i
=0; (i
<b
->nr
); i
++) {
970 for(j
=b
->index
[i
]; (j
<b
->index
[i
+1]); j
++) {
972 fprintf(stderr
,"\nWARNING: negative index %d in group %s\n\n",
973 b
->a
[j
],(*grpname
)[i
]);
980 static void minstring(char *str
)
984 for (i
=0; (i
< (int)strlen(str
)); i
++)
989 int find_group(char s
[], int ngrps
, char **grpname
)
998 /* first look for whole name match */
1000 for(i
=0; i
<ngrps
; i
++)
1001 if (gmx_strcasecmp_min(s
,grpname
[i
])==0) {
1006 /* second look for first string match */
1008 for(i
=0; i
<ngrps
; i
++)
1009 if (gmx_strncasecmp_min(s
,grpname
[i
],n
)==0) {
1014 /* last look for arbitrary substring match */
1018 for(i
=0; i
<ngrps
; i
++) {
1019 strcpy(string
, grpname
[i
]);
1022 if (strstr(string
,s
)!=NULL
) {
1030 printf("Error: Multiple groups '%s' selected\n", s
);
1036 static int qgroup(int *a
, int ngrps
, char **grpname
)
1044 fprintf(stderr
,"Select a group: ");
1046 if ( scanf("%s",s
)!=1 )
1047 gmx_fatal(FARGS
,"Cannot read from input");
1048 trim(s
); /* remove spaces */
1049 } while (strlen(s
)==0);
1050 aa
= strtol(s
, &end
, 10);
1051 if (aa
==0 && end
[0] != '\0') /* string entered */
1052 aa
= find_group(s
, ngrps
, grpname
);
1053 bInRange
= (aa
>= 0 && aa
< ngrps
);
1055 printf("Error: No such group '%s'\n", s
);
1056 } while (!bInRange
);
1057 printf("Selected %d: '%s'\n", aa
, grpname
[aa
]);
1062 static void rd_groups(t_blocka
*grps
,char **grpname
,char *gnames
[],
1063 int ngrps
,int isize
[],atom_id
*index
[],int grpnr
[])
1068 gmx_fatal(FARGS
,"Error: no groups in indexfile");
1069 for(i
=0; (i
<grps
->nr
); i
++)
1070 fprintf(stderr
,"Group %5d (%15s) has %5d elements\n",i
,grpname
[i
],
1071 grps
->index
[i
+1]-grps
->index
[i
]);
1072 for(i
=0; (i
<ngrps
); i
++) {
1075 gnr1
=qgroup(&grpnr
[i
], grps
->nr
, grpname
);
1076 if ((gnr1
<0) || (gnr1
>=grps
->nr
))
1077 fprintf(stderr
,"Select between %d and %d.\n",0,grps
->nr
-1);
1078 } while ((gnr1
<0) || (gnr1
>=grps
->nr
));
1080 fprintf(stderr
,"There is one group in the index\n");
1083 gnames
[i
]=strdup(grpname
[gnr1
]);
1084 isize
[i
]=grps
->index
[gnr1
+1]-grps
->index
[gnr1
];
1085 snew(index
[i
],isize
[i
]);
1086 for(j
=0; (j
<isize
[i
]); j
++)
1087 index
[i
][j
]=grps
->a
[grps
->index
[gnr1
]+j
];
1091 void rd_index(const char *statfile
,int ngrps
,int isize
[],
1092 atom_id
*index
[],char *grpnames
[])
1100 gmx_fatal(FARGS
,"No index file specified");
1101 grps
=init_index(statfile
,&gnames
);
1102 rd_groups(grps
,gnames
,grpnames
,ngrps
,isize
,index
,grpnr
);
1105 void rd_index_nrs(char *statfile
,int ngrps
,int isize
[],
1106 atom_id
*index
[],char *grpnames
[],int grpnr
[])
1112 gmx_fatal(FARGS
,"No index file specified");
1113 grps
=init_index(statfile
,&gnames
);
1115 rd_groups(grps
,gnames
,grpnames
,ngrps
,isize
,index
,grpnr
);
1118 void get_index(t_atoms
*atoms
, const char *fnm
, int ngrps
,
1119 int isize
[], atom_id
*index
[],char *grpnames
[])
1122 t_blocka
*grps
= NULL
;
1128 grps
=init_index(fnm
,gnames
);
1132 snew(grps
->index
,1);
1133 analyse(atoms
,grps
,gnames
,FALSE
,FALSE
);
1136 gmx_incons("You need to supply a valid atoms structure or a valid index file name");
1138 rd_groups(grps
,*gnames
,grpnames
,ngrps
,isize
,index
,grpnr
);
1141 t_cluster_ndx
*cluster_index(FILE *fplog
,const char *ndx
)
1147 c
->clust
= init_index(ndx
,&c
->grpname
);
1149 for(i
=0; (i
<c
->clust
->nra
); i
++)
1150 c
->maxframe
= max(c
->maxframe
,c
->clust
->a
[i
]);
1151 fprintf(fplog
? fplog
: stdout
,
1152 "There are %d clusters containing %d structures, highest framenr is %d\n",
1153 c
->clust
->nr
,c
->clust
->nra
,c
->maxframe
);
1155 pr_blocka(debug
,0,"clust",c
->clust
,TRUE
);
1156 for(i
=0; (i
<c
->clust
->nra
); i
++)
1157 if ((c
->clust
->a
[i
] < 0) || (c
->clust
->a
[i
] > c
->maxframe
))
1158 gmx_fatal(FARGS
,"Range check error for c->clust->a[%d] = %d\n"
1159 "should be within 0 and %d",i
,c
->clust
->a
[i
],c
->maxframe
+1);
1161 c
->inv_clust
=make_invblocka(c
->clust
,c
->maxframe
);