Enforced rotation: fixed torque calculation for FLEX potential when using mass-weighting
[gromacs/adressmacs.git] / src / tools / make_ndx.c
blobee4311d6c98b2a9549dc049b6818a254ec787f75
1 /*
2 *
3 * This source code is part of
4 *
5 * G R O M A C S
6 *
7 * GROningen MAchine for Chemical Simulations
8 *
9 * VERSION 3.2.0
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
32 * And Hey:
33 * Green Red Orange Magenta Azure Cyan Skyblue
35 #ifdef HAVE_CONFIG_H
36 #include <config.h>
37 #endif
39 #include <ctype.h>
40 #include "sysstuff.h"
41 #include "strdb.h"
42 #include "futil.h"
43 #include "macros.h"
44 #include "string2.h"
45 #include "statutil.h"
46 #include "confio.h"
47 #include "copyrite.h"
48 #include "typedefs.h"
49 #include "index.h"
50 #include "smalloc.h"
51 #include "vec.h"
52 #include "index.h"
54 #define MAXNAMES 30
55 #define NAME_LEN 30
57 gmx_bool bCase=FALSE;
59 static int or_groups(atom_id nr1,atom_id *at1,atom_id nr2,atom_id *at2,
60 atom_id *nr,atom_id *at)
62 atom_id i1,i2,max=0;
63 gmx_bool bNotIncr;
65 *nr=0;
67 bNotIncr=FALSE;
68 for(i1=0; i1<nr1; i1++) {
69 if ((i1>0) && (at1[i1] <= max))
70 bNotIncr=TRUE;
71 max=at1[i1];
73 for(i1=0; i1<nr2; i1++) {
74 if ((i1>0) && (at2[i1] <= max))
75 bNotIncr=TRUE;
76 max=at2[i1];
79 if (bNotIncr)
80 printf("One of your groups is not ascending\n");
81 else {
82 i1=0;
83 i2=0;
84 *nr=0;
85 while ((i1 < nr1) || (i2 < nr2)) {
86 if ((i2 == nr2) || ((i1<nr1) && (at1[i1] < at2[i2]))) {
87 at[*nr]=at1[i1];
88 (*nr)++;
89 i1++;
91 else {
92 if ((i2<nr2) && ((i1==nr1) || (at1[i1] > at2[i2]))) {
93 at[*nr]=at2[i2];
94 (*nr)++;
96 i2++;
100 printf("Merged two groups with OR: %u %u -> %u\n",nr1,nr2,*nr);
103 return *nr;
106 static int and_groups(atom_id nr1,atom_id *at1,atom_id nr2,atom_id *at2,
107 atom_id *nr,atom_id *at)
109 atom_id i1,i2;
111 *nr=0;
112 for (i1=0; i1<nr1; i1++) {
113 for (i2=0; i2<nr2; i2++) {
114 if (at1[i1]==at2[i2]) {
115 at[*nr]=at1[i1];
116 (*nr)++;
121 printf("Merged two groups with AND: %u %u -> %u\n",nr1,nr2,*nr);
123 return *nr;
126 static gmx_bool is_name_char(char c)
128 /* This string should contain all characters that can not be
129 * the first letter of a name due to the make_ndx syntax.
131 const char *spec=" !&|";
133 return (c != '\0' && strchr(spec,c) == NULL);
136 static int parse_names(char **string,int *n_names,char **names)
138 int i;
140 *n_names=0;
141 while (is_name_char((*string)[0]) || (*string)[0] == ' ') {
142 if (is_name_char((*string)[0])) {
143 if (*n_names >= MAXNAMES)
144 gmx_fatal(FARGS,"To many names: %d\n",*n_names+1);
145 i=0;
146 while (is_name_char((*string)[i])) {
147 names[*n_names][i]=(*string)[i];
148 i++;
149 if (i > NAME_LEN) {
150 printf("Name is too long, the maximum is %d characters\n",NAME_LEN);
151 return 0;
154 names[*n_names][i]='\0';
155 if (!bCase)
156 upstring(names[*n_names]);
157 *string += i;
158 (*n_names)++;
160 else
161 (*string)++;
164 return *n_names;
167 static gmx_bool parse_int_char(char **string,int *nr,char *c)
169 char *orig;
170 gmx_bool bRet;
172 orig = *string;
174 while ((*string)[0]==' ')
175 (*string)++;
177 bRet=FALSE;
179 *c = ' ';
181 if (isdigit((*string)[0])) {
182 *nr=(*string)[0]-'0';
183 (*string)++;
184 while (isdigit((*string)[0])) {
185 *nr = (*nr)*10+(*string)[0]-'0';
186 (*string)++;
188 if (isalpha((*string)[0])) {
189 *c = (*string)[0];
190 (*string)++;
192 /* Check if there is at most one non-digit character */
193 if (!isalnum((*string)[0])) {
194 bRet = TRUE;
195 } else {
196 *string = orig;
199 else
200 *nr=NOTSET;
202 return bRet;
205 static gmx_bool parse_int(char **string,int *nr)
207 char *orig,c;
208 gmx_bool bRet;
210 orig = *string;
211 bRet = parse_int_char(string,nr,&c);
212 if (bRet && c != ' ') {
213 *string = orig;
214 bRet = FALSE;
217 return bRet;
220 static gmx_bool isquote(char c)
222 return (c == '\"');
225 static gmx_bool parse_string(char **string,int *nr, int ngrps, char **grpname)
227 char *s, *sp;
228 char c;
230 while ((*string)[0]==' ')
231 (*string)++;
233 (*nr) = NOTSET;
234 if (isquote((*string)[0])) {
235 c=(*string)[0];
236 (*string)++;
237 s = strdup((*string));
238 sp = strchr(s, c);
239 if (sp!=NULL) {
240 (*string) += sp-s + 1;
241 sp[0]='\0';
242 (*nr) = find_group(s, ngrps, grpname);
246 return (*nr) != NOTSET;
249 static int select_atomnumbers(char **string,t_atoms *atoms,atom_id n1,
250 atom_id *nr,atom_id *index,char *gname)
252 char buf[STRLEN];
253 int i,up;
255 *nr=0;
256 while ((*string)[0]==' ')
257 (*string)++;
258 if ((*string)[0]=='-') {
259 (*string)++;
260 parse_int(string,&up);
261 if ((n1<1) || (n1>atoms->nr) || (up<1) || (up>atoms->nr))
262 printf("Invalid atom range\n");
263 else {
264 for(i=n1-1; i<=up-1; i++) {
265 index[*nr]=i;
266 (*nr)++;
268 printf("Found %u atom%s in range %u-%d\n",*nr,(*nr==1)?"":"s",n1,up);
269 if (n1==up)
270 sprintf(buf,"a_%u",n1);
271 else
272 sprintf(buf,"a_%u-%d",n1,up);
273 strcpy(gname,buf);
276 else {
277 i=n1;
278 sprintf(gname,"a");
279 do {
280 if ((i-1>=0) && (i-1<atoms->nr)) {
281 index[*nr] = i-1;
282 (*nr)++;
283 sprintf(buf,"_%d",i);
284 strcat(gname,buf);
285 } else {
286 printf("Invalid atom number %d\n",i);
287 *nr = 0;
289 } while ((*nr!=0) && (parse_int(string,&i)));
292 return *nr;
295 static int select_residuenumbers(char **string,t_atoms *atoms,
296 atom_id n1,char c,
297 atom_id *nr,atom_id *index,char *gname)
299 char buf[STRLEN];
300 int i,j,up;
301 t_resinfo *ri;
303 *nr=0;
304 while ((*string)[0]==' ')
305 (*string)++;
306 if ((*string)[0]=='-') {
307 /* Residue number range selection */
308 if (c != ' ') {
309 printf("Error: residue insertion codes can not be used with residue range selection\n");
310 return 0;
312 (*string)++;
313 parse_int(string,&up);
315 for(i=0; i<atoms->nr; i++) {
316 ri = &atoms->resinfo[atoms->atom[i].resind];
317 for(j=n1; (j<=up); j++) {
318 if (ri->nr == j && (c == ' ' || ri->ic == c)) {
319 index[*nr]=i;
320 (*nr)++;
324 printf("Found %u atom%s with res.nr. in range %u-%d\n",
325 *nr,(*nr==1)?"":"s",n1,up);
326 if (n1==up)
327 sprintf(buf,"r_%u",n1);
328 else
329 sprintf(buf,"r_%u-%d",n1,up);
330 strcpy(gname,buf);
332 else {
333 /* Individual residue number/insertion code selection */
334 j=n1;
335 sprintf(gname,"r");
336 do {
337 for(i=0; i<atoms->nr; i++) {
338 ri = &atoms->resinfo[atoms->atom[i].resind];
339 if (ri->nr == j && ri->ic == c) {
340 index[*nr]=i;
341 (*nr)++;
344 sprintf(buf,"_%d",j);
345 strcat(gname,buf);
346 } while (parse_int_char(string,&j,&c));
349 return *nr;
352 static int select_residueindices(char **string,t_atoms *atoms,
353 atom_id n1,char c,
354 atom_id *nr,atom_id *index,char *gname)
356 /*this should be similar to select_residuenumbers except select by index (sequential numbering in file)*/
357 /*resind+1 for 1-indexing*/
358 char buf[STRLEN];
359 int i,j,up;
360 t_resinfo *ri;
362 *nr=0;
363 while ((*string)[0]==' ')
364 (*string)++;
365 if ((*string)[0]=='-') {
366 /* Residue number range selection */
367 if (c != ' ') {
368 printf("Error: residue insertion codes can not be used with residue range selection\n");
369 return 0;
371 (*string)++;
372 parse_int(string,&up);
374 for(i=0; i<atoms->nr; i++) {
375 ri = &atoms->resinfo[atoms->atom[i].resind];
376 for(j=n1; (j<=up); j++) {
377 if (atoms->atom[i].resind+1 == j && (c == ' ' || ri->ic == c)) {
378 index[*nr]=i;
379 (*nr)++;
383 printf("Found %u atom%s with resind.+1 in range %u-%d\n",
384 *nr,(*nr==1)?"":"s",n1,up);
385 if (n1==up)
386 sprintf(buf,"r_%u",n1);
387 else
388 sprintf(buf,"r_%u-%d",n1,up);
389 strcpy(gname,buf);
391 else {
392 /* Individual residue number/insertion code selection */
393 j=n1;
394 sprintf(gname,"r");
395 do {
396 for(i=0; i<atoms->nr; i++) {
397 ri = &atoms->resinfo[atoms->atom[i].resind];
398 if (atoms->atom[i].resind+1 == j && ri->ic == c) {
399 index[*nr]=i;
400 (*nr)++;
403 sprintf(buf,"_%d",j);
404 strcat(gname,buf);
405 } while (parse_int_char(string,&j,&c));
408 return *nr;
412 static gmx_bool atoms_from_residuenumbers(t_atoms *atoms,int group,t_blocka *block,
413 atom_id *nr,atom_id *index,char *gname)
415 int i,j,j0,j1,resnr,nres;
417 j0=block->index[group];
418 j1=block->index[group+1];
419 nres = atoms->nres;
420 for(j=j0; j<j1; j++)
421 if (block->a[j]>=nres) {
422 printf("Index %s contains number>nres (%d>%d)\n",
423 gname,block->a[j]+1,nres);
424 return FALSE;
426 for(i=0; i<atoms->nr; i++) {
427 resnr = atoms->resinfo[atoms->atom[i].resind].nr;
428 for (j=j0; j<j1; j++)
429 if (block->a[j]+1 == resnr) {
430 index[*nr]=i;
431 (*nr)++;
432 break;
435 printf("Found %u atom%s in %d residues from group %s\n",
436 *nr,(*nr==1)?"":"s",j1-j0,gname);
437 return *nr;
440 static gmx_bool comp_name(char *name,char *search)
442 while (name[0] != '\0' && search[0] != '\0') {
443 switch (search[0]) {
444 case '?':
445 /* Always matches */
446 break;
447 case '*':
448 if (search[1] != '\0') {
449 printf("WARNING: Currently '*' is only supported at the end of an expression\n");
450 return FALSE;
451 } else {
452 return TRUE;
454 break;
455 default:
456 /* Compare a single character */
457 if (( bCase && strncmp(name,search,1)) ||
458 (!bCase && gmx_strncasecmp(name,search,1))) {
459 return FALSE;
462 name++;
463 search++;
466 return (name[0] == '\0' && (search[0] == '\0' || search[0] == '*'));
469 static int select_chainnames(t_atoms *atoms,int n_names,char **names,
470 atom_id *nr,atom_id *index)
472 char name[2];
473 int j;
474 atom_id i;
476 name[1]=0;
477 *nr=0;
478 for(i=0; i<atoms->nr; i++) {
479 name[0] = atoms->resinfo[atoms->atom[i].resind].chainid;
480 j=0;
481 while (j<n_names && !comp_name(name,names[j]))
482 j++;
483 if (j<n_names) {
484 index[*nr]=i;
485 (*nr)++;
488 printf("Found %u atom%s with chain identifier%s",
489 *nr,(*nr==1)?"":"s",(n_names==1)?"":"s");
490 for(j=0; (j<n_names); j++)
491 printf(" %s",names[j]);
492 printf("\n");
494 return *nr;
497 static int select_atomnames(t_atoms *atoms,int n_names,char **names,
498 atom_id *nr,atom_id *index,gmx_bool bType)
500 char *name;
501 int j;
502 atom_id i;
504 *nr=0;
505 for(i=0; i<atoms->nr; i++) {
506 if (bType)
507 name=*(atoms->atomtype[i]);
508 else
509 name=*(atoms->atomname[i]);
510 j=0;
511 while (j<n_names && !comp_name(name,names[j]))
512 j++;
513 if (j<n_names) {
514 index[*nr]=i;
515 (*nr)++;
518 printf("Found %u atoms with %s%s",
519 *nr,bType ? "type" : "name",(n_names==1)?"":"s");
520 for(j=0; (j<n_names); j++)
521 printf(" %s",names[j]);
522 printf("\n");
524 return *nr;
527 static int select_residuenames(t_atoms *atoms,int n_names,char **names,
528 atom_id *nr,atom_id *index)
530 char *name;
531 int j;
532 atom_id i;
534 *nr=0;
535 for(i=0; i<atoms->nr; i++) {
536 name = *(atoms->resinfo[atoms->atom[i].resind].name);
537 j=0;
538 while (j<n_names && !comp_name(name,names[j]))
539 j++;
540 if (j<n_names) {
541 index[*nr]=i;
542 (*nr)++;
545 printf("Found %u atoms with residue name%s",*nr,(n_names==1)?"":"s");
546 for(j=0; (j<n_names); j++)
547 printf(" %s",names[j]);
548 printf("\n");
550 return *nr;
553 static void copy2block(int n,atom_id *index,t_blocka *block)
555 int i,n0;
557 block->nr++;
558 n0=block->nra;
559 block->nra=n0+n;
560 srenew(block->index,block->nr+1);
561 block->index[block->nr]=n0+n;
562 srenew(block->a,n0+n);
563 for(i=0; (i<n); i++)
564 block->a[n0+i]=index[i];
567 static void make_gname(int n,char **names,char *gname)
569 int i;
571 strcpy(gname,names[0]);
572 for (i=1; i<n; i++) {
573 strcat(gname,"_");
574 strcat(gname,names[i]);
578 static void copy_group(int g,t_blocka *block,atom_id *nr,atom_id *index)
580 int i,i0;
582 i0=block->index[g];
583 *nr=block->index[g+1]-i0;
584 for (i=0; i<*nr; i++)
585 index[i]=block->a[i0+i];
588 static void remove_group(int nr,int nr2,t_blocka *block,char ***gn)
590 int i,j,shift;
591 char *name;
593 if (nr2==NOTSET)
594 nr2=nr;
596 for(j=0; j<=nr2-nr; j++) {
597 if ((nr<0) || (nr>=block->nr))
598 printf("Group %d does not exist\n",nr+j);
599 else {
600 shift=block->index[nr+1]-block->index[nr];
601 for(i=block->index[nr+1]; i<block->nra; i++)
602 block->a[i-shift]=block->a[i];
604 for(i=nr; i<block->nr; i++)
605 block->index[i]=block->index[i+1]-shift;
606 name = strdup((*gn)[nr]);
607 sfree((*gn)[nr]);
608 for(i=nr; i<block->nr-1; i++) {
609 (*gn)[i]=(*gn)[i+1];
611 block->nr--;
612 block->nra=block->index[block->nr];
613 printf("Removed group %d '%s'\n",nr+j,name);
614 sfree(name);
619 static void split_group(t_atoms *atoms,int sel_nr,t_blocka *block,char ***gn,
620 gmx_bool bAtom)
622 char buf[STRLEN],*name;
623 int i,resind;
624 atom_id a,n0,n1;
626 printf("Splitting group %d '%s' into %s\n",sel_nr,(*gn)[sel_nr],
627 bAtom ? "atoms" : "residues");
629 n0=block->index[sel_nr];
630 n1=block->index[sel_nr+1];
631 srenew(block->a,block->nra+n1-n0);
632 for (i=n0; i<n1; i++) {
633 a=block->a[i];
634 resind = atoms->atom[a].resind;
635 name = *(atoms->resinfo[resind].name);
636 if (bAtom || (i==n0) || (atoms->atom[block->a[i-1]].resind!=resind)) {
637 if (i>n0)
638 block->index[block->nr]=block->nra;
639 block->nr++;
640 srenew(block->index,block->nr+1);
641 srenew(*gn,block->nr);
642 if (bAtom)
643 sprintf(buf,"%s_%s_%u",(*gn)[sel_nr],*atoms->atomname[a],a+1);
644 else
645 sprintf(buf,"%s_%s_%d",(*gn)[sel_nr],name,atoms->resinfo[resind].nr);
646 (*gn)[block->nr-1]=strdup(buf);
648 block->a[block->nra]=a;
649 block->nra++;
651 block->index[block->nr]=block->nra;
654 static int split_chain(t_atoms *atoms,rvec *x,
655 int sel_nr,t_blocka *block,char ***gn)
657 char buf[STRLEN];
658 int j,nchain;
659 atom_id i,a,natoms,*start=NULL,*end=NULL,ca_start,ca_end;
660 rvec vec;
662 natoms=atoms->nr;
663 nchain=0;
664 ca_start=0;
666 while (ca_start<natoms) {
667 while((ca_start<natoms) && strcmp(*atoms->atomname[ca_start],"CA"))
668 ca_start++;
669 if (ca_start<natoms) {
670 srenew(start,nchain+1);
671 srenew(end,nchain+1);
672 start[nchain]=ca_start;
673 while ((start[nchain]>0) &&
674 (atoms->atom[start[nchain]-1].resind ==
675 atoms->atom[ca_start].resind))
676 start[nchain]--;
678 i=ca_start;
679 do {
680 ca_end=i;
681 do {
682 i++;
683 } while ((i<natoms) && strcmp(*atoms->atomname[i],"CA"));
684 if (i<natoms)
685 rvec_sub(x[ca_end],x[i],vec);
686 } while ((i<natoms) && (norm(vec)<0.45));
688 end[nchain]=ca_end;
689 while ((end[nchain]+1<natoms) &&
690 (atoms->atom[end[nchain]+1].resind==atoms->atom[ca_end].resind))
691 end[nchain]++;
692 ca_start=end[nchain]+1;
693 nchain++;
696 if (nchain==1)
697 printf("Found 1 chain, will not split\n");
698 else
699 printf("Found %d chains\n",nchain);
700 for (j=0; j<nchain; j++)
701 printf("%d:%6u atoms (%u to %u)\n",
702 j+1,end[j]-start[j]+1,start[j]+1,end[j]+1);
704 if (nchain > 1) {
705 srenew(block->a,block->nra+block->index[sel_nr+1]-block->index[sel_nr]);
706 for (j=0; j<nchain; j++) {
707 block->nr++;
708 srenew(block->index,block->nr+1);
709 srenew(*gn,block->nr);
710 sprintf(buf,"%s_chain%d",(*gn)[sel_nr],j+1);
711 (*gn)[block->nr-1]=strdup(buf);
712 for (i=block->index[sel_nr]; i<block->index[sel_nr+1]; i++) {
713 a=block->a[i];
714 if ((a>=start[j]) && (a<=end[j])) {
715 block->a[block->nra]=a;
716 block->nra++;
719 block->index[block->nr]=block->nra;
720 if (block->index[block->nr-1]==block->index[block->nr])
721 remove_group(block->nr-1,NOTSET,block,gn);
724 sfree(start);
725 sfree(end);
727 return nchain;
730 static gmx_bool check_have_atoms(t_atoms *atoms, char *string)
732 if ( atoms==NULL ) {
733 printf("Can not process '%s' without atom info, use option -f\n", string);
734 return FALSE;
735 } else
736 return TRUE;
739 static gmx_bool parse_entry(char **string,int natoms,t_atoms *atoms,
740 t_blocka *block,char ***gn,
741 atom_id *nr,atom_id *index,char *gname)
743 static char **names, *ostring;
744 static gmx_bool bFirst=TRUE;
745 int j,n_names,sel_nr1;
746 atom_id i,nr1,*index1;
747 char c;
748 gmx_bool bRet,bCompl;
750 if (bFirst) {
751 bFirst=FALSE;
752 snew(names,MAXNAMES);
753 for (i=0; i<MAXNAMES; i++)
754 snew(names[i],NAME_LEN+1);
757 bRet=FALSE;
758 sel_nr1=NOTSET;
760 while(*string[0]==' ')
761 (*string)++;
763 if ((*string)[0]=='!') {
764 bCompl=TRUE;
765 (*string)++;
766 while(*string[0]==' ')
767 (*string)++;
768 } else
769 bCompl=FALSE;
771 ostring = *string;
773 if (parse_int(string,&sel_nr1) ||
774 parse_string(string,&sel_nr1,block->nr,*gn)) {
775 if ((sel_nr1>=0) && (sel_nr1<block->nr)) {
776 copy_group(sel_nr1,block,nr,index);
777 strcpy(gname,(*gn)[sel_nr1]);
778 printf("Copied index group %d '%s'\n",sel_nr1,(*gn)[sel_nr1]);
779 bRet=TRUE;
780 } else
781 printf("Group %d does not exist\n",sel_nr1);
783 else if ((*string)[0]=='a') {
784 (*string)++;
785 if (check_have_atoms(atoms, ostring)) {
786 if (parse_int(string,&sel_nr1)) {
787 bRet=select_atomnumbers(string,atoms,sel_nr1,nr,index,gname);
789 else if (parse_names(string,&n_names,names)) {
790 bRet=select_atomnames(atoms,n_names,names,nr,index,FALSE);
791 make_gname(n_names,names,gname);
795 else if ((*string)[0]=='t') {
796 (*string)++;
797 if (check_have_atoms(atoms, ostring) &&
798 parse_names(string,&n_names,names)) {
799 if (atoms->atomtype == NULL)
800 printf("Need a run input file to select atom types\n");
801 else {
802 bRet=select_atomnames(atoms,n_names,names,nr,index,TRUE);
803 make_gname(n_names,names,gname);
807 else if (strncmp(*string,"res",3)==0) {
808 (*string)+=3;
809 if ( check_have_atoms(atoms, ostring) &&
810 parse_int(string,&sel_nr1) &&
811 (sel_nr1>=0) && (sel_nr1<block->nr) ) {
812 bRet=atoms_from_residuenumbers(atoms,
813 sel_nr1,block,nr,index,(*gn)[sel_nr1]);
814 sprintf(gname,"atom_%s",(*gn)[sel_nr1]);
817 else if (strncmp(*string,"ri",2)==0) {
818 (*string)+=2;
819 if (check_have_atoms(atoms, ostring) &&
820 parse_int_char(string,&sel_nr1,&c)) {
821 bRet=select_residueindices(string,atoms,sel_nr1,c,nr,index,gname);
824 else if ((*string)[0]=='r') {
825 (*string)++;
826 if (check_have_atoms(atoms, ostring)) {
827 if (parse_int_char(string,&sel_nr1,&c)) {
828 bRet=select_residuenumbers(string,atoms,sel_nr1,c,nr,index,gname);
830 else if (parse_names(string,&n_names,names)) {
831 bRet=select_residuenames(atoms,n_names,names,nr,index);
832 make_gname(n_names,names,gname);
836 else if (strncmp(*string,"chain",5)==0) {
837 (*string)+=5;
838 if (check_have_atoms(atoms, ostring) &&
839 parse_names(string,&n_names,names)) {
840 bRet=select_chainnames(atoms,n_names,names,nr,index);
841 sprintf(gname,"ch%s",names[0]);
842 for (i=1; i<n_names; i++)
843 strcat(gname,names[i]);
846 if (bRet && bCompl) {
847 snew(index1,natoms-*nr);
848 nr1=0;
849 for(i=0; i<natoms; i++) {
850 j=0;
851 while ((j<*nr) && (index[j] != i))
852 j++;
853 if (j==*nr) {
854 if (nr1 >= natoms-*nr) {
855 printf("There are double atoms in your index group\n");
856 break;
858 index1[nr1]=i;
859 nr1++;
862 *nr=nr1;
863 for(i=0; i<nr1; i++)
864 index[i]=index1[i];
865 sfree(index1);
867 for (i=strlen(gname)+1; i>0; i--)
868 gname[i]=gname[i-1];
869 gname[0]='!';
870 printf("Complemented group: %u atoms\n",*nr);
873 return bRet;
876 static void list_residues(t_atoms *atoms)
878 int i,j,start,end,prev_resind,resind;
879 gmx_bool bDiff;
881 /* Print all the residues, assuming continuous resnr count */
882 start = atoms->atom[0].resind;
883 prev_resind = start;
884 for(i=0; i<atoms->nr; i++) {
885 resind = atoms->atom[i].resind;
886 if ((resind != prev_resind) || (i==atoms->nr-1)) {
887 if ((bDiff=strcmp(*atoms->resinfo[resind].name,
888 *atoms->resinfo[start].name)) ||
889 (i==atoms->nr-1)) {
890 if (bDiff)
891 end = prev_resind;
892 else
893 end = resind;
894 if (end < start+3)
895 for(j=start; j<=end; j++)
896 printf("%4d %-5s",
897 j+1,*(atoms->resinfo[j].name));
898 else
899 printf(" %4d - %4d %-5s ",
900 start+1,end+1,*(atoms->resinfo[start].name));
901 start = resind;
904 prev_resind = resind;
906 printf("\n");
909 static void edit_index(int natoms, t_atoms *atoms,rvec *x,t_blocka *block, char ***gn, gmx_bool bVerbose)
911 static char **atnames, *ostring;
912 static gmx_bool bFirst=TRUE;
913 char inp_string[STRLEN],*string;
914 char gname[STRLEN],gname1[STRLEN],gname2[STRLEN];
915 int i,i0,i1,sel_nr,sel_nr2,newgroup;
916 atom_id nr,nr1,nr2,*index,*index1,*index2;
917 gmx_bool bAnd,bOr,bPrintOnce;
919 if (bFirst) {
920 bFirst=FALSE;
921 snew(atnames,MAXNAMES);
922 for (i=0; i<MAXNAMES; i++)
923 snew(atnames[i],NAME_LEN+1);
926 string=NULL;
928 snew(index,natoms);
929 snew(index1,natoms);
930 snew(index2,natoms);
932 newgroup=NOTSET;
933 bPrintOnce = TRUE;
934 do {
935 gname1[0]='\0';
936 if (bVerbose || bPrintOnce || newgroup!=NOTSET) {
937 printf("\n");
938 if (bVerbose || bPrintOnce || newgroup==NOTSET) {
939 i0=0;
940 i1=block->nr;
941 } else {
942 i0=newgroup;
943 i1=newgroup+1;
945 for(i=i0; i<i1; i++)
946 printf("%3d %-20s: %5u atoms\n",i,(*gn)[i],
947 block->index[i+1]-block->index[i]);
948 newgroup=NOTSET;
950 if (bVerbose || bPrintOnce) {
951 printf("\n");
952 printf(" nr : group ! 'name' nr name 'splitch' nr Enter: list groups\n");
953 printf(" 'a': atom & 'del' nr 'splitres' nr 'l': list residues\n");
954 printf(" 't': atom type | 'keep' nr 'splitat' nr 'h': help\n");
955 printf(" 'r': residue 'res' nr 'chain' char\n");
956 printf(" \"name\": group 'case': case %s 'q': save and quit\n",
957 bCase ? "insensitive" : "sensitive ");
958 printf(" 'ri': residue index\n");
959 bPrintOnce = FALSE;
961 printf("\n");
962 printf("> ");
963 if(NULL==fgets(inp_string,STRLEN,stdin))
965 gmx_fatal(FARGS,"Error reading user input");
967 inp_string[strlen(inp_string)-1]=0;
968 printf("\n");
969 string=inp_string;
970 while (string[0]==' ')
971 string++;
973 ostring = string;
974 nr=0;
975 if (string[0] == 'h') {
976 printf(" nr : selects an index group by number or quoted string.\n");
977 printf(" The string is first matched against the whole group name,\n");
978 printf(" then against the beginning and finally against an\n");
979 printf(" arbitrary substring. A multiple match is an error.\n");
981 printf(" 'a' nr1 [nr2 ...] : selects atoms, atom numbering starts at 1.\n");
982 printf(" 'a' nr1 - nr2 : selects atoms in the range from nr1 to nr2.\n");
983 printf(" 'a' name1[*] [name2[*] ...] : selects atoms by name(s), '?' matches any char,\n");
984 printf(" wildcard '*' allowed at the end of a name.\n");
985 printf(" 't' type1[*] [type2[*] ...] : as 'a', but for type, run input file required.\n");
986 printf(" 'r' nr1[ic1] [nr2[ic2] ...] : selects residues by number and insertion code.\n");
987 printf(" 'r' nr1 - nr2 : selects residues in the range from nr1 to nr2.\n");
988 printf(" 'r' name1[*] [name2[*] ...] : as 'a', but for residue names.\n");
989 printf(" 'ri' nr1 - nr2 : selects residue indices, 1-indexed, (as opposed to numbers) in the range from nr1 to nr2.\n");
990 printf(" 'chain' ch1 [ch2 ...] : selects atoms by chain identifier(s),\n");
991 printf(" not available with a .gro file as input.\n");
992 printf(" ! : takes the complement of a group with respect to all\n");
993 printf(" the atoms in the input file.\n");
994 printf(" & | : AND and OR, can be placed between any of the options\n");
995 printf(" above, the input is processed from left to right.\n");
996 printf(" 'name' nr name : rename group nr to name.\n");
997 printf(" 'del' nr1 [- nr2] : deletes one group or groups in the range from nr1 to nr2.\n");
998 printf(" 'keep' nr : deletes all groups except nr.\n");
999 printf(" 'case' : make all name compares case (in)sensitive.\n");
1000 printf(" 'splitch' nr : split group into chains using CA distances.\n");
1001 printf(" 'splitres' nr : split group into residues.\n");
1002 printf(" 'splitat' nr : split group into atoms.\n");
1003 printf(" 'res' nr : interpret numbers in group as residue numbers\n");
1004 printf(" Enter : list the currently defined groups and commands\n");
1005 printf(" 'l' : list the residues.\n");
1006 printf(" 'h' : show this help.\n");
1007 printf(" 'q' : save and quit.\n");
1008 printf("\n");
1009 printf(" Examples:\n");
1010 printf(" > 2 | 4 & r 3-5\n");
1011 printf(" selects all atoms from group 2 and 4 that have residue numbers 3, 4 or 5\n");
1012 printf(" > a C* & !a C CA\n");
1013 printf(" selects all atoms starting with 'C' but not the atoms 'C' and 'CA'\n");
1014 printf(" > \"protein\" & ! \"backb\"\n");
1015 printf(" selects all atoms that are in group 'protein' and not in group 'backbone'\n");
1016 if (bVerbose) {
1017 printf("\npress Enter ");
1018 getchar();
1021 else if (strncmp(string,"del",3)==0) {
1022 string+=3;
1023 if (parse_int(&string,&sel_nr)) {
1024 while(string[0]==' ')
1025 string++;
1026 if (string[0]=='-') {
1027 string++;
1028 parse_int(&string,&sel_nr2);
1029 } else
1030 sel_nr2=NOTSET;
1031 while(string[0]==' ')
1032 string++;
1033 if (string[0]=='\0')
1034 remove_group(sel_nr,sel_nr2,block,gn);
1035 else
1036 printf("\nSyntax error: \"%s\"\n",string);
1039 else if (strncmp(string,"keep",4)==0) {
1040 string+=4;
1041 if (parse_int(&string,&sel_nr)) {
1042 remove_group(sel_nr+1,block->nr-1,block,gn);
1043 remove_group(0,sel_nr-1,block,gn);
1046 else if (strncmp(string,"name",4)==0) {
1047 string+=4;
1048 if (parse_int(&string,&sel_nr)) {
1049 if ((sel_nr>=0) && (sel_nr<block->nr)) {
1050 sscanf(string,"%s",gname);
1051 sfree((*gn)[sel_nr]);
1052 (*gn)[sel_nr]=strdup(gname);
1056 else if (strncmp(string,"case",4)==0) {
1057 bCase=!bCase;
1058 printf("Switched to case %s\n",bCase ? "sensitive" : "insensitive");
1060 else if (string[0] == 'v' ) {
1061 bVerbose=!bVerbose;
1062 printf("Turned verbose %s\n",bVerbose ? "on" : "off");
1064 else if (string[0] == 'l') {
1065 if ( check_have_atoms(atoms, ostring) )
1066 list_residues(atoms);
1068 else if (strncmp(string,"splitch",7)==0) {
1069 string+=7;
1070 if ( check_have_atoms(atoms, ostring) &&
1071 parse_int(&string,&sel_nr) &&
1072 (sel_nr>=0) && (sel_nr<block->nr))
1073 split_chain(atoms,x,sel_nr,block,gn);
1075 else if (strncmp(string,"splitres",8)==0 ) {
1076 string+=8;
1077 if ( check_have_atoms(atoms, ostring) &&
1078 parse_int(&string,&sel_nr) &&
1079 (sel_nr>=0) && (sel_nr<block->nr))
1080 split_group(atoms,sel_nr,block,gn,FALSE);
1082 else if (strncmp(string,"splitat",7)==0 ) {
1083 string+=7;
1084 if ( check_have_atoms(atoms, ostring) &&
1085 parse_int(&string,&sel_nr) &&
1086 (sel_nr>=0) && (sel_nr<block->nr))
1087 split_group(atoms,sel_nr,block,gn,TRUE);
1089 else if (string[0] == '\0') {
1090 bPrintOnce = TRUE;
1092 else if (string[0] != 'q') {
1093 nr1=-1;
1094 nr2=-1;
1095 if (parse_entry(&string,natoms,atoms,block,gn,&nr,index,gname)) {
1096 do {
1097 while (string[0]==' ')
1098 string++;
1100 bAnd=FALSE;
1101 bOr=FALSE;
1102 if (string[0]=='&')
1103 bAnd=TRUE;
1104 else if (string[0]=='|')
1105 bOr=TRUE;
1107 if (bAnd || bOr) {
1108 string++;
1109 nr1=nr;
1110 for(i=0; i<nr; i++)
1111 index1[i]=index[i];
1112 strcpy(gname1,gname);
1113 if (parse_entry(&string,natoms,atoms,block,gn,&nr2,index2,gname2)) {
1114 if (bOr) {
1115 or_groups(nr1,index1,nr2,index2,&nr,index);
1116 sprintf(gname,"%s_%s",gname1,gname2);
1118 else {
1119 and_groups(nr1,index1,nr2,index2,&nr,index);
1120 sprintf(gname,"%s_&_%s",gname1,gname2);
1124 } while (bAnd || bOr);
1126 while(string[0]==' ')
1127 string++;
1128 if (string[0])
1129 printf("\nSyntax error: \"%s\"\n",string);
1130 else if (nr>0) {
1131 copy2block(nr,index,block);
1132 srenew(*gn,block->nr);
1133 newgroup = block->nr-1;
1134 (*gn)[newgroup]=strdup(gname);
1136 else
1137 printf("Group is empty\n");
1139 } while (string[0]!='q');
1141 sfree(index);
1142 sfree(index1);
1143 sfree(index2);
1146 static int block2natoms(t_blocka *block)
1148 int i, natoms;
1150 natoms = 0;
1151 for(i=0; i<block->nra; i++)
1152 natoms = max(natoms, block->a[i]);
1153 natoms++;
1155 return natoms;
1158 void merge_blocks(t_blocka *dest, t_blocka *source)
1160 int i,nra0,i0;
1162 /* count groups, srenew and fill */
1163 i0 = dest->nr;
1164 nra0 = dest->nra;
1165 dest->nr+=source->nr;
1166 srenew(dest->index, dest->nr+1);
1167 for(i=0; i<source->nr; i++)
1168 dest->index[i0+i] = nra0 + source->index[i];
1169 /* count atoms, srenew and fill */
1170 dest->nra+=source->nra;
1171 srenew(dest->a, dest->nra);
1172 for(i=0; i<source->nra; i++)
1173 dest->a[nra0+i] = source->a[i];
1175 /* terminate list */
1176 dest->index[dest->nr]=dest->nra;
1180 int main(int argc,char *argv[])
1182 const char *desc[] = {
1183 "Index groups are necessary for almost every gromacs program.",
1184 "All these programs can generate default index groups. You ONLY",
1185 "have to use make_ndx when you need SPECIAL index groups.",
1186 "There is a default index group for the whole system, 9 default",
1187 "index groups are generated for proteins, a default index group",
1188 "is generated for every other residue name.[PAR]",
1189 "When no index file is supplied, also make_ndx will generate the",
1190 "default groups.",
1191 "With the index editor you can select on atom, residue and chain names",
1192 "and numbers.",
1193 "When a run input file is supplied you can also select on atom type.",
1194 "You can use NOT, AND and OR, you can split groups",
1195 "into chains, residues or atoms. You can delete and rename groups.[PAR]",
1196 "The atom numbering in the editor and the index file starts at 1."
1199 static int natoms=0;
1200 static gmx_bool bVerbose=FALSE;
1201 t_pargs pa[] = {
1202 { "-natoms", FALSE, etINT, {&natoms},
1203 "set number of atoms (default: read from coordinate or index file)" },
1204 { "-verbose", FALSE, etBOOL, {&bVerbose},
1205 "HIDDENVerbose output" }
1207 #define NPA asize(pa)
1209 output_env_t oenv;
1210 char title[STRLEN];
1211 int nndxin;
1212 const char *stxfile;
1213 char **ndxinfiles;
1214 const char *ndxoutfile;
1215 gmx_bool bNatoms;
1216 int i,j;
1217 t_atoms *atoms;
1218 rvec *x,*v;
1219 int ePBC;
1220 matrix box;
1221 t_blocka *block,*block2;
1222 char **gnames,**gnames2;
1223 t_filenm fnm[] = {
1224 { efSTX, "-f", NULL, ffOPTRD },
1225 { efNDX, "-n", NULL, ffOPTRDMULT },
1226 { efNDX, "-o", NULL, ffWRITE }
1228 #define NFILE asize(fnm)
1230 CopyRight(stderr,argv[0]);
1232 parse_common_args(&argc,argv,0,NFILE,fnm,NPA,pa,asize(desc),desc,
1233 0,NULL,&oenv);
1235 stxfile = ftp2fn_null(efSTX,NFILE,fnm);
1236 if (opt2bSet("-n",NFILE,fnm)) {
1237 nndxin = opt2fns(&ndxinfiles,"-n",NFILE,fnm);
1238 } else {
1239 nndxin = 0;
1241 ndxoutfile = opt2fn("-o",NFILE,fnm);
1242 bNatoms = opt2parg_bSet("-natoms",NPA,pa);
1244 if (!stxfile && !nndxin)
1245 gmx_fatal(FARGS,"No input files (structure or index)");
1247 if (stxfile) {
1248 snew(atoms,1);
1249 get_stx_coordnum(stxfile,&(atoms->nr));
1250 init_t_atoms(atoms,atoms->nr,TRUE);
1251 snew(x,atoms->nr);
1252 snew(v,atoms->nr);
1253 fprintf(stderr,"\nReading structure file\n");
1254 read_stx_conf(stxfile,title,atoms,x,v,&ePBC,box);
1255 natoms = atoms->nr;
1256 bNatoms=TRUE;
1257 } else {
1258 atoms = NULL;
1259 x = NULL;
1262 /* read input file(s) */
1263 block = new_blocka();
1264 gnames = NULL;
1265 printf("Going to read %d old index file(s)\n",nndxin);
1266 if (nndxin) {
1267 for(i=0; i<nndxin; i++) {
1268 block2 = init_index(ndxinfiles[i],&gnames2);
1269 srenew(gnames, block->nr+block2->nr);
1270 for(j=0; j<block2->nr; j++)
1271 gnames[block->nr+j]=gnames2[j];
1272 sfree(gnames2);
1273 merge_blocks(block, block2);
1274 sfree(block2->a);
1275 sfree(block2->index);
1276 /* done_block(block2); */
1277 sfree(block2);
1280 else {
1281 snew(gnames,1);
1282 analyse(atoms,block,&gnames,FALSE,TRUE);
1285 if (!bNatoms) {
1286 natoms = block2natoms(block);
1287 printf("Counted atom numbers up to %d in index file\n", natoms);
1290 edit_index(natoms,atoms,x,block,&gnames,bVerbose);
1292 write_index(ndxoutfile,block,gnames);
1294 thanx(stderr);
1296 return 0;