?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
Parallel vs. sequentiual code: I get binary identical result (apart from comment...
[gromacs/qmmm-gamess-us.git] / src / tools / make_ndx.c
blob1d9b7c5ff8e935422fe45af6715e0f36bb6e8160
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.
14
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.
19 *
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.
26 *
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.
29 *
30 * For more info, check our website at http://www.gromacs.org
31 *
32 * And Hey:
33 * Green Red Orange Magenta Azure Cyan Skyblue
34 */
35 #ifdef HAVE_CONFIG_H
36 #include <config.h>
37 #endif
38
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"
53
54 #define MAXNAMES 30
55 #define NAME_LEN 30
56
57 bool bCase=FALSE;
58
59 static int or_groups(atom_id nr1,atom_id *at1,atom_id nr2,atom_id *at2,
60 atom_id *nr,atom_id *at)
61 {
62 atom_id i1,i2,max=0;
63 bool bNotIncr;
64
65 *nr=0;
66
67 bNotIncr=FALSE;
68 for(i1=0; i1<nr1; i1++) {
69 if ((i1>0) && (at1[i1] <= max))
70 bNotIncr=TRUE;
71 max=at1[i1];
72 }
73 for(i1=0; i1<nr2; i1++) {
74 if ((i1>0) && (at2[i1] <= max))
75 bNotIncr=TRUE;
76 max=at2[i1];
77 }
78
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++;
90 }
91 else {
92 if ((i2<nr2) && ((i1==nr1) || (at1[i1] > at2[i2]))) {
93 at[*nr]=at2[i2];
94 (*nr)++;
95 }
96 i2++;
97 }
98 }
99
100 printf("Merged two groups with OR: %u %u -> %u\n",nr1,nr2,*nr);
101 }
102
103 return *nr;
104 }
105
106 static int and_groups(atom_id nr1,atom_id *at1,atom_id nr2,atom_id *at2,
107 atom_id *nr,atom_id *at)
108 {
109 atom_id i1,i2;
110
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)++;
117 }
118 }
119 }
120
121 printf("Merged two groups with AND: %u %u -> %u\n",nr1,nr2,*nr);
122
123 return *nr;
124 }
125
126 static bool isalnum_star(char c)
127 {
128 return (isalnum(c) || (c=='_') || (c=='*') || (c=='?'));
129 }
130
131 static int parse_names(char **string,int *n_names,char **names)
132 {
133 int i;
134
135 *n_names=0;
136 while ((isalnum_star((*string)[0]) || ((*string)[0]==' '))) {
137 if (isalnum_star((*string)[0])) {
138 if (*n_names >= MAXNAMES)
139 gmx_fatal(FARGS,"To many names: %d\n",*n_names+1);
140 i=0;
141 while (isalnum_star((*string)[i])) {
142 names[*n_names][i]=(*string)[i];
143 i++;
144 if (i > NAME_LEN) {
145 printf("Name is too long, the maximum is %d characters\n",NAME_LEN);
146 return 0;
147 }
148 }
149 names[*n_names][i]='\0';
150 if (!bCase)
151 upstring(names[*n_names]);
152 *string += i;
153 (*n_names)++;
154 }
155 else
156 (*string)++;
157 }
158
159 return *n_names;
160 }
161
162 static bool parse_int_char(char **string,int *nr,char *c)
163 {
164 char *orig;
165 bool bRet;
166
167 orig = *string;
168
169 while ((*string)[0]==' ')
170 (*string)++;
171
172 bRet=FALSE;
173
174 *c = ' ';
175
176 if (isdigit((*string)[0])) {
177 *nr=(*string)[0]-'0';
178 (*string)++;
179 while (isdigit((*string)[0])) {
180 *nr = (*nr)*10+(*string)[0]-'0';
181 (*string)++;
182 }
183 if (isalpha((*string)[0])) {
184 *c = (*string)[0];
185 (*string)++;
186 }
187 /* Check if there is at most one non-digit character */
188 if (!isalnum((*string)[0])) {
189 bRet = TRUE;
190 } else {
191 *string = orig;
192 }
193 }
194 else
195 *nr=NOTSET;
196
197 return bRet;
198 }
199
200 static bool parse_int(char **string,int *nr)
201 {
202 char *orig,c;
203 bool bRet;
204
205 orig = *string;
206 bRet = parse_int_char(string,nr,&c);
207 if (bRet && c != ' ') {
208 *string = orig;
209 bRet = FALSE;
210 }
211
212 return bRet;
213 }
214
215 static bool isquote(char c)
216 {
217 char quotes[] = "'\"";
218
219 return strchr(quotes, c)!=NULL;
220 }
221
222 static bool parse_string(char **string,int *nr, int ngrps, char **grpname)
223 {
224 char *s, *sp;
225 char c;
226
227 while ((*string)[0]==' ')
228 (*string)++;
229
230 (*nr) = NOTSET;
231 if (isquote((*string)[0])) {
232 c=(*string)[0];
233 (*string)++;
234 s = strdup((*string));
235 sp = strchr(s, c);
236 if (sp!=NULL) {
237 (*string) += sp-s + 1;
238 sp[0]='\0';
239 (*nr) = find_group(s, ngrps, grpname);
240 }
241 }
242
243 return (*nr) != NOTSET;
244 }
245
246 static int select_atomnumbers(char **string,t_atoms *atoms,atom_id n1,
247 atom_id *nr,atom_id *index,char *gname)
248 {
249 char buf[STRLEN];
250 int i,up;
251
252 *nr=0;
253 while ((*string)[0]==' ')
254 (*string)++;
255 if ((*string)[0]=='-') {
256 (*string)++;
257 parse_int(string,&up);
258 if ((n1<1) || (n1>atoms->nr) || (up<1) || (up>atoms->nr))
259 printf("Invalid atom range\n");
260 else {
261 for(i=n1-1; i<=up-1; i++) {
262 index[*nr]=i;
263 (*nr)++;
264 }
265 printf("Found %u atom%s in range %u-%d\n",*nr,(*nr==1)?"":"s",n1,up);
266 if (n1==up)
267 sprintf(buf,"a_%u",n1);
268 else
269 sprintf(buf,"a_%u-%d",n1,up);
270 strcpy(gname,buf);
271 }
272 }
273 else {
274 i=n1;
275 sprintf(gname,"a");
276 do {
277 if ((i-1>=0) && (i-1<atoms->nr)) {
278 index[*nr] = i-1;
279 (*nr)++;
280 sprintf(buf,"_%d",i);
281 strcat(gname,buf);
282 } else {
283 printf("Invalid atom number %d\n",i);
284 *nr = 0;
285 }
286 } while ((*nr!=0) && (parse_int(string,&i)));
287 }
288
289 return *nr;
290 }
291
292 static int select_residuenumbers(char **string,t_atoms *atoms,
293 atom_id n1,char c,
294 atom_id *nr,atom_id *index,char *gname)
295 {
296 char buf[STRLEN];
297 int i,j,up;
298 t_resinfo *ri;
299
300 *nr=0;
301 while ((*string)[0]==' ')
302 (*string)++;
303 if ((*string)[0]=='-') {
304 /* Residue number range selection */
305 if (c != ' ') {
306 printf("Error: residue insertion codes can not be used with residue range selection\n");
307 return 0;
308 }
309 (*string)++;
310 parse_int(string,&up);
311
312 for(i=0; i<atoms->nr; i++) {
313 ri = &atoms->resinfo[atoms->atom[i].resind];
314 for(j=n1; (j<=up); j++) {
315 if (ri->nr == j && (c == ' ' || ri->ic == c)) {
316 index[*nr]=i;
317 (*nr)++;
318 }
319 }
320 }
321 printf("Found %u atom%s with res.nr. in range %u-%d\n",
322 *nr,(*nr==1)?"":"s",n1,up);
323 if (n1==up)
324 sprintf(buf,"r_%u",n1);
325 else
326 sprintf(buf,"r_%u-%d",n1,up);
327 strcpy(gname,buf);
328 }
329 else {
330 /* Individual residue number/insertion code selection */
331 j=n1;
332 sprintf(gname,"r");
333 do {
334 for(i=0; i<atoms->nr; i++) {
335 ri = &atoms->resinfo[atoms->atom[i].resind];
336 if (ri->nr == j && ri->ic == c) {
337 index[*nr]=i;
338 (*nr)++;
339 }
340 }
341 sprintf(buf,"_%d",j);
342 strcat(gname,buf);
343 } while (parse_int_char(string,&j,&c));
344 }
345
346 return *nr;
347 }
348
349 static int select_residueindices(char **string,t_atoms *atoms,
350 atom_id n1,char c,
351 atom_id *nr,atom_id *index,char *gname)
352 {
353 /*this should be similar to select_residuenumbers except select by index (sequential numbering in file)*/
354 /*resind+1 for 1-indexing*/
355 char buf[STRLEN];
356 int i,j,up;
357 t_resinfo *ri;
358
359 *nr=0;
360 while ((*string)[0]==' ')
361 (*string)++;
362 if ((*string)[0]=='-') {
363 /* Residue number range selection */
364 if (c != ' ') {
365 printf("Error: residue insertion codes can not be used with residue range selection\n");
366 return 0;
367 }
368 (*string)++;
369 parse_int(string,&up);
370
371 for(i=0; i<atoms->nr; i++) {
372 ri = &atoms->resinfo[atoms->atom[i].resind];
373 for(j=n1; (j<=up); j++) {
374 if (atoms->atom[i].resind+1 == j && (c == ' ' || ri->ic == c)) {
375 index[*nr]=i;
376 (*nr)++;
377 }
378 }
379 }
380 printf("Found %u atom%s with resind.+1 in range %u-%d\n",
381 *nr,(*nr==1)?"":"s",n1,up);
382 if (n1==up)
383 sprintf(buf,"r_%u",n1);
384 else
385 sprintf(buf,"r_%u-%d",n1,up);
386 strcpy(gname,buf);
387 }
388 else {
389 /* Individual residue number/insertion code selection */
390 j=n1;
391 sprintf(gname,"r");
392 do {
393 for(i=0; i<atoms->nr; i++) {
394 ri = &atoms->resinfo[atoms->atom[i].resind];
395 if (atoms->atom[i].resind+1 == j && ri->ic == c) {
396 index[*nr]=i;
397 (*nr)++;
398 }
399 }
400 sprintf(buf,"_%d",j);
401 strcat(gname,buf);
402 } while (parse_int_char(string,&j,&c));
403 }
404
405 return *nr;
406 }
407
408
409 static bool atoms_from_residuenumbers(t_atoms *atoms,int group,t_blocka *block,
410 atom_id *nr,atom_id *index,char *gname)
411 {
412 int i,j,j0,j1,resnr,nres;
413
414 j0=block->index[group];
415 j1=block->index[group+1];
416 nres = atoms->nres;
417 for(j=j0; j<j1; j++)
418 if (block->a[j]>=nres) {
419 printf("Index %s contains number>nres (%d>%d)\n",
420 gname,block->a[j]+1,nres);
421 return FALSE;
422 }
423 for(i=0; i<atoms->nr; i++) {
424 resnr = atoms->resinfo[atoms->atom[i].resind].nr;
425 for (j=j0; j<j1; j++)
426 if (block->a[j]+1 == resnr) {
427 index[*nr]=i;
428 (*nr)++;
429 break;
430 }
431 }
432 printf("Found %u atom%s in %d residues from group %s\n",
433 *nr,(*nr==1)?"":"s",j1-j0,gname);
434 return *nr;
435 }
436
437 static bool comp_name(char *name,char *search)
438 {
439 while (name[0] != '\0' && search[0] != '\0') {
440 switch (search[0]) {
441 case '?':
442 /* Always matches */
443 break;
444 case '*':
445 if (search[1] != '\0') {
446 printf("WARNING: Currently '*' is only supported at the end of an expression\n");
447 return FALSE;
448 } else {
449 return TRUE;
450 }
451 break;
452 default:
453 /* Compare a single character */
454 if (( bCase && strncmp(name,search,1)) ||
455 (!bCase && strncasecmp(name,search,1))) {
456 return FALSE;
457 }
458 }
459 name++;
460 search++;
461 }
462
463 return (name[0] == '\0' && (search[0] == '\0' || search[0] == '*'));
464 }
465
466 static int select_chainnames(t_atoms *atoms,int n_names,char **names,
467 atom_id *nr,atom_id *index)
468 {
469 char name[2];
470 int j;
471 atom_id i;
472
473 name[1]=0;
474 *nr=0;
475 for(i=0; i<atoms->nr; i++) {
476 name[0] = atoms->resinfo[atoms->atom[i].resind].chain;
477 j=0;
478 while (j<n_names && !comp_name(name,names[j]))
479 j++;
480 if (j<n_names) {
481 index[*nr]=i;
482 (*nr)++;
483 }
484 }
485 printf("Found %u atom%s with chain identifier%s",
486 *nr,(*nr==1)?"":"s",(n_names==1)?"":"s");
487 for(j=0; (j<n_names); j++)
488 printf(" %s",names[j]);
489 printf("\n");
490
491 return *nr;
492 }
493
494 static int select_atomnames(t_atoms *atoms,int n_names,char **names,
495 atom_id *nr,atom_id *index,bool bType)
496 {
497 char *name;
498 int j;
499 atom_id i;
500
501 *nr=0;
502 for(i=0; i<atoms->nr; i++) {
503 if (bType)
504 name=*(atoms->atomtype[i]);
505 else
506 name=*(atoms->atomname[i]);
507 j=0;
508 while (j<n_names && !comp_name(name,names[j]))
509 j++;
510 if (j<n_names) {
511 index[*nr]=i;
512 (*nr)++;
513 }
514 }
515 printf("Found %u atoms with %s%s",
516 *nr,bType ? "type" : "name",(n_names==1)?"":"s");
517 for(j=0; (j<n_names); j++)
518 printf(" %s",names[j]);
519 printf("\n");
520
521 return *nr;
522 }
523
524 static int select_residuenames(t_atoms *atoms,int n_names,char **names,
525 atom_id *nr,atom_id *index)
526 {
527 char *name;
528 int j;
529 atom_id i;
530
531 *nr=0;
532 for(i=0; i<atoms->nr; i++) {
533 name = *(atoms->resinfo[atoms->atom[i].resind].name);
534 j=0;
535 while (j<n_names && !comp_name(name,names[j]))
536 j++;
537 if (j<n_names) {
538 index[*nr]=i;
539 (*nr)++;
540 }
541 }
542 printf("Found %u atoms with residue name%s",*nr,(n_names==1)?"":"s");
543 for(j=0; (j<n_names); j++)
544 printf(" %s",names[j]);
545 printf("\n");
546
547 return *nr;
548 }
549
550 static void copy2block(int n,atom_id *index,t_blocka *block)
551 {
552 int i,n0;
553
554 block->nr++;
555 n0=block->nra;
556 block->nra=n0+n;
557 srenew(block->index,block->nr+1);
558 block->index[block->nr]=n0+n;
559 srenew(block->a,n0+n);
560 for(i=0; (i<n); i++)
561 block->a[n0+i]=index[i];
562 }
563
564 static void make_gname(int n,char **names,char *gname)
565 {
566 int i;
567
568 strcpy(gname,names[0]);
569 for (i=1; i<n; i++) {
570 strcat(gname,"_");
571 strcat(gname,names[i]);
572 }
573 }
574
575 static void copy_group(int g,t_blocka *block,atom_id *nr,atom_id *index)
576 {
577 int i,i0;
578
579 i0=block->index[g];
580 *nr=block->index[g+1]-i0;
581 for (i=0; i<*nr; i++)
582 index[i]=block->a[i0+i];
583 }
584
585 static void remove_group(int nr,int nr2,t_blocka *block,char ***gn)
586 {
587 int i,j,shift;
588 char *name;
589
590 if (nr2==NOTSET)
591 nr2=nr;
592
593 for(j=0; j<=nr2-nr; j++) {
594 if ((nr<0) || (nr>=block->nr))
595 printf("Group %d does not exist\n",nr+j);
596 else {
597 shift=block->index[nr+1]-block->index[nr];
598 for(i=block->index[nr+1]; i<block->nra; i++)
599 block->a[i-shift]=block->a[i];
600
601 for(i=nr; i<block->nr; i++)
602 block->index[i]=block->index[i+1]-shift;
603 name = strdup((*gn)[nr]);
604 sfree((*gn)[nr]);
605 for(i=nr; i<block->nr-1; i++) {
606 (*gn)[i]=(*gn)[i+1];
607 }
608 block->nr--;
609 block->nra=block->index[block->nr];
610 printf("Removed group %d '%s'\n",nr+j,name);
611 sfree(name);
612 }
613 }
614 }
615
616 static void split_group(t_atoms *atoms,int sel_nr,t_blocka *block,char ***gn,
617 bool bAtom)
618 {
619 char buf[STRLEN],*name;
620 int i,resind;
621 atom_id a,n0,n1;
622
623 printf("Splitting group %d '%s' into %s\n",sel_nr,(*gn)[sel_nr],
624 bAtom ? "atoms" : "residues");
625
626 n0=block->index[sel_nr];
627 n1=block->index[sel_nr+1];
628 srenew(block->a,block->nra+n1-n0);
629 for (i=n0; i<n1; i++) {
630 a=block->a[i];
631 resind = atoms->atom[a].resind;
632 name = *(atoms->resinfo[resind].name);
633 if (bAtom || (i==n0) || (atoms->atom[block->a[i-1]].resind!=resind)) {
634 if (i>n0)
635 block->index[block->nr]=block->nra;
636 block->nr++;
637 srenew(block->index,block->nr+1);
638 srenew(*gn,block->nr);
639 if (bAtom)
640 sprintf(buf,"%s_%s_%u",(*gn)[sel_nr],*atoms->atomname[a],a+1);
641 else
642 sprintf(buf,"%s_%s_%d",(*gn)[sel_nr],name,atoms->resinfo[resind].nr);
643 (*gn)[block->nr-1]=strdup(buf);
644 }
645 block->a[block->nra]=a;
646 block->nra++;
647 }
648 block->index[block->nr]=block->nra;
649 }
650
651 static int split_chain(t_atoms *atoms,rvec *x,
652 int sel_nr,t_blocka *block,char ***gn)
653 {
654 char buf[STRLEN];
655 int j,nchain;
656 atom_id i,a,natoms,*start=NULL,*end=NULL,ca_start,ca_end;
657 rvec vec;
658
659 natoms=atoms->nr;
660 nchain=0;
661 ca_start=0;
662
663 while (ca_start<natoms) {
664 while((ca_start<natoms) && strcmp(*atoms->atomname[ca_start],"CA"))
665 ca_start++;
666 if (ca_start<natoms) {
667 srenew(start,nchain+1);
668 srenew(end,nchain+1);
669 start[nchain]=ca_start;
670 while ((start[nchain]>0) &&
671 (atoms->atom[start[nchain]-1].resind ==
672 atoms->atom[ca_start].resind))
673 start[nchain]--;
674
675 i=ca_start;
676 do {
677 ca_end=i;
678 do {
679 i++;
680 } while ((i<natoms) && strcmp(*atoms->atomname[i],"CA"));
681 if (i<natoms)
682 rvec_sub(x[ca_end],x[i],vec);
683 } while ((i<natoms) && (norm(vec)<0.45));
684
685 end[nchain]=ca_end;
686 while ((end[nchain]+1<natoms) &&
687 (atoms->atom[end[nchain]+1].resind==atoms->atom[ca_end].resind))
688 end[nchain]++;
689 ca_start=end[nchain]+1;
690 nchain++;
691 }
692 }
693 if (nchain==1)
694 printf("Found 1 chain, will not split\n");
695 else
696 printf("Found %d chains\n",nchain);
697 for (j=0; j<nchain; j++)
698 printf("%d:%6u atoms (%u to %u)\n",
699 j+1,end[j]-start[j]+1,start[j]+1,end[j]+1);
700
701 if (nchain > 1) {
702 srenew(block->a,block->nra+block->index[sel_nr+1]-block->index[sel_nr]);
703 for (j=0; j<nchain; j++) {
704 block->nr++;
705 srenew(block->index,block->nr+1);
706 srenew(*gn,block->nr);
707 sprintf(buf,"%s_chain%d",(*gn)[sel_nr],j+1);
708 (*gn)[block->nr-1]=strdup(buf);
709 for (i=block->index[sel_nr]; i<block->index[sel_nr+1]; i++) {
710 a=block->a[i];
711 if ((a>=start[j]) && (a<=end[j])) {
712 block->a[block->nra]=a;
713 block->nra++;
714 }
715 }
716 block->index[block->nr]=block->nra;
717 if (block->index[block->nr-1]==block->index[block->nr])
718 remove_group(block->nr-1,NOTSET,block,gn);
719 }
720 }
721 sfree(start);
722 sfree(end);
723
724 return nchain;
725 }
726
727 static bool check_have_atoms(t_atoms *atoms, char *string)
728 {
729 if ( atoms==NULL ) {
730 printf("Can not process '%s' without atoms info\n", string);
731 return FALSE;
732 } else
733 return TRUE;
734 }
735
736 static bool parse_entry(char **string,int natoms,t_atoms *atoms,
737 t_blocka *block,char ***gn,
738 atom_id *nr,atom_id *index,char *gname)
739 {
740 static char **names, *ostring;
741 static bool bFirst=TRUE;
742 int j,n_names,sel_nr1;
743 atom_id i,nr1,*index1;
744 char c;
745 bool bRet,bCompl;
746
747 if (bFirst) {
748 bFirst=FALSE;
749 snew(names,MAXNAMES);
750 for (i=0; i<MAXNAMES; i++)
751 snew(names[i],NAME_LEN+1);
752 }
753
754 bRet=FALSE;
755 sel_nr1=NOTSET;
756
757 while(*string[0]==' ')
758 (*string)++;
759
760 if ((*string)[0]=='!') {
761 bCompl=TRUE;
762 (*string)++;
763 while(*string[0]==' ')
764 (*string)++;
765 } else
766 bCompl=FALSE;
767
768 ostring = *string;
769
770 if (parse_int(string,&sel_nr1) ||
771 parse_string(string,&sel_nr1,block->nr,*gn)) {
772 if ((sel_nr1>=0) && (sel_nr1<block->nr)) {
773 copy_group(sel_nr1,block,nr,index);
774 strcpy(gname,(*gn)[sel_nr1]);
775 printf("Copied index group %d '%s'\n",sel_nr1,(*gn)[sel_nr1]);
776 bRet=TRUE;
777 } else
778 printf("Group %d does not exist\n",sel_nr1);
779 }
780 else if ((*string)[0]=='a') {
781 (*string)++;
782 if (check_have_atoms(atoms, ostring)) {
783 if (parse_int(string,&sel_nr1)) {
784 bRet=select_atomnumbers(string,atoms,sel_nr1,nr,index,gname);
785 }
786 else if (parse_names(string,&n_names,names)) {
787 bRet=select_atomnames(atoms,n_names,names,nr,index,FALSE);
788 make_gname(n_names,names,gname);
789 }
790 }
791 }
792 else if ((*string)[0]=='t') {
793 (*string)++;
794 if (check_have_atoms(atoms, ostring) &&
795 parse_names(string,&n_names,names)) {
796 if (atoms->atomtype == NULL)
797 printf("Need a run input file to select atom types\n");
798 else {
799 bRet=select_atomnames(atoms,n_names,names,nr,index,TRUE);
800 make_gname(n_names,names,gname);
801 }
802 }
803 }
804 else if (strncmp(*string,"res",3)==0) {
805 (*string)+=3;
806 if ( check_have_atoms(atoms, ostring) &&
807 parse_int(string,&sel_nr1) &&
808 (sel_nr1>=0) && (sel_nr1<block->nr) ) {
809 bRet=atoms_from_residuenumbers(atoms,
810 sel_nr1,block,nr,index,(*gn)[sel_nr1]);
811 sprintf(gname,"atom_%s",(*gn)[sel_nr1]);
812 }
813 }
814 else if (strncmp(*string,"ri",2)==0) {
815 (*string)+=2;
816 if (check_have_atoms(atoms, ostring) &&
817 parse_int_char(string,&sel_nr1,&c)) {
818 bRet=select_residueindices(string,atoms,sel_nr1,c,nr,index,gname);
819 }
820 }
821 else if ((*string)[0]=='r') {
822 (*string)++;
823 if (check_have_atoms(atoms, ostring)) {
824 if (parse_int_char(string,&sel_nr1,&c)) {
825 bRet=select_residuenumbers(string,atoms,sel_nr1,c,nr,index,gname);
826 }
827 else if (parse_names(string,&n_names,names)) {
828 bRet=select_residuenames(atoms,n_names,names,nr,index);
829 make_gname(n_names,names,gname);
830 }
831 }
832 }
833 else if (strncmp(*string,"chain",5)==0) {
834 (*string)+=5;
835 if (check_have_atoms(atoms, ostring) &&
836 parse_names(string,&n_names,names)) {
837 bRet=select_chainnames(atoms,n_names,names,nr,index);
838 sprintf(gname,"ch%s",names[0]);
839 for (i=1; i<n_names; i++)
840 strcat(gname,names[i]);
841 }
842 }
843 if (bRet && bCompl) {
844 snew(index1,natoms-*nr);
845 nr1=0;
846 for(i=0; i<natoms; i++) {
847 j=0;
848 while ((j<*nr) && (index[j] != i))
849 j++;
850 if (j==*nr) {
851 if (nr1 >= natoms-*nr) {
852 printf("There are double atoms in your index group\n");
853 break;
854 }
855 index1[nr1]=i;
856 nr1++;
857 }
858 }
859 *nr=nr1;
860 for(i=0; i<nr1; i++)
861 index[i]=index1[i];
862 sfree(index1);
863
864 for (i=strlen(gname)+1; i>0; i--)
865 gname[i]=gname[i-1];
866 gname[0]='!';
867 printf("Complemented group: %u atoms\n",*nr);
868 }
869
870 return bRet;
871 }
872
873 static void list_residues(t_atoms *atoms)
874 {
875 int i,j,start,end,prev_resind,resind;
876 bool bDiff;
877
878 /* Print all the residues, assuming continuous resnr count */
879 start = atoms->atom[0].resind;
880 prev_resind = start;
881 for(i=0; i<atoms->nr; i++) {
882 resind = atoms->atom[i].resind;
883 if ((resind != prev_resind) || (i==atoms->nr-1)) {
884 if ((bDiff=strcmp(*atoms->resinfo[resind].name,
885 *atoms->resinfo[start].name)) ||
886 (i==atoms->nr-1)) {
887 if (bDiff)
888 end = prev_resind;
889 else
890 end = resind;
891 if (end < start+3)
892 for(j=start; j<=end; j++)
893 printf("%4d %-5s",
894 j+1,*(atoms->resinfo[j].name));
895 else
896 printf(" %4d - %4d %-5s ",
897 start+1,end+1,*(atoms->resinfo[start].name));
898 start = resind;
899 }
900 }
901 prev_resind = resind;
902 }
903 printf("\n");
904 }
905
906 static void edit_index(int natoms, t_atoms *atoms,rvec *x,t_blocka *block, char ***gn, bool bVerbose)
907 {
908 static char **atnames, *ostring;
909 static bool bFirst=TRUE;
910 char inp_string[STRLEN],*string;
911 char gname[STRLEN],gname1[STRLEN],gname2[STRLEN];
912 int i,i0,i1,sel_nr,sel_nr2,newgroup;
913 atom_id nr,nr1,nr2,*index,*index1,*index2;
914 bool bAnd,bOr,bPrintOnce;
915
916 if (bFirst) {
917 bFirst=FALSE;
918 snew(atnames,MAXNAMES);
919 for (i=0; i<MAXNAMES; i++)
920 snew(atnames[i],NAME_LEN+1);
921 }
922
923 string=NULL;
924
925 snew(index,natoms);
926 snew(index1,natoms);
927 snew(index2,natoms);
928
929 newgroup=NOTSET;
930 bPrintOnce = TRUE;
931 do {
932 gname1[0]='\0';
933 if (bVerbose || bPrintOnce || newgroup!=NOTSET) {
934 printf("\n");
935 if (bVerbose || bPrintOnce || newgroup==NOTSET) {
936 i0=0;
937 i1=block->nr;
938 } else {
939 i0=newgroup;
940 i1=newgroup+1;
941 }
942 for(i=i0; i<i1; i++)
943 printf("%3d %-20s: %5u atoms\n",i,(*gn)[i],
944 block->index[i+1]-block->index[i]);
945 newgroup=NOTSET;
946 }
947 if (bVerbose || bPrintOnce) {
948 printf("\n");
949 printf(" nr : group ! 'name' nr name 'splitch' nr Enter: list groups\n");
950 printf(" 'a': atom & 'del' nr 'splitres' nr 'l': list residues\n");
951 printf(" 't': atom type | 'keep' nr 'splitat' nr 'h': help\n");
952 printf(" 'r': residue 'res' nr 'chain' char\n");
953 printf(" \"name\": group 'case': case %s 'q': save and quit\n",
954 bCase ? "insensitive" : "sensitive ");
955 printf(" 'ri': residue index\n");
956 bPrintOnce = FALSE;
957 }
958 printf("\n");
959 printf("> ");
960 if(NULL==fgets(inp_string,STRLEN,stdin))
961 {
962 gmx_fatal(FARGS,"Error reading user input");
963 }
964 inp_string[strlen(inp_string)-1]=0;
965 printf("\n");
966 string=inp_string;
967 while (string[0]==' ')
968 string++;
969
970 ostring = string;
971 nr=0;
972 if (string[0] == 'h') {
973 printf(" nr : selects an index group by number or quoted string.\n");
974 printf(" The string is first matched against the whole group name,\n");
975 printf(" then against the beginning and finally against an\n");
976 printf(" arbitrary substring. A multiple match is an error.\n");
977
978 printf(" 'a' nr1 [nr2 ...] : selects atoms, atom numbering starts at 1.\n");
979 printf(" 'a' nr1 - nr2 : selects atoms in the range from nr1 to nr2.\n");
980 printf(" 'a' name1[*] [name2[*] ...] : selects atoms by name(s), '?' matches any char,\n");
981 printf(" wildcard '*' allowed at the end of a name.\n");
982 printf(" 't' type1[*] [type2[*] ...] : as 'a', but for type, run input file required.\n");
983 printf(" 'r' nr1[ic1] [nr2[ic2] ...] : selects residues by number and insertion code.\n");
984 printf(" 'r' nr1 - nr2 : selects residues in the range from nr1 to nr2.\n");
985 printf(" 'r' name1[*] [name2[*] ...] : as 'a', but for residue names.\n");
986 printf(" 'ri' nr1 - nr2 : selects residue indices, 1-indexed, (as opposed to numbers) in the range from nr1 to nr2.\n");
987 printf(" 'chain' ch1 [ch2 ...] : selects atoms by chain identifier(s),\n");
988 printf(" not available with a .gro file as input.\n");
989 printf(" ! : takes the complement of a group with respect to all\n");
990 printf(" the atoms in the input file.\n");
991 printf(" & | : AND and OR, can be placed between any of the options\n");
992 printf(" above, the input is processed from left to right.\n");
993 printf(" 'name' nr name : rename group nr to name.\n");
994 printf(" 'del' nr1 [- nr2] : deletes one group or groups in the range from nr1 to nr2.\n");
995 printf(" 'keep' nr : deletes all groups except nr.\n");
996 printf(" 'case' : make all name compares case (in)sensitive.\n");
997 printf(" 'splitch' nr : split group into chains using CA distances.\n");
998 printf(" 'splitres' nr : split group into residues.\n");
999 printf(" 'splitat' nr : split group into atoms.\n");
1000 printf(" 'res' nr : interpret numbers in group as residue numbers\n");
1001 printf(" Enter : list the currently defined groups and commands\n");
1002 printf(" 'l' : list the residues.\n");
1003 printf(" 'h' : show this help.\n");
1004 printf(" 'q' : save and quit.\n");
1005 printf("\n");
1006 printf(" Examples:\n");
1007 printf(" > 2 | 4 & r 3-5\n");
1008 printf(" selects all atoms from group 2 and 4 that have residue numbers 3, 4 or 5\n");
1009 printf(" > a C* & !a C CA\n");
1010 printf(" selects all atoms starting with 'C' but not the atoms 'C' and 'CA'\n");
1011 printf(" > \"protein\" & ! \"backb\"\n");
1012 printf(" selects all atoms that are in group 'protein' and not in group 'backbone'\n");
1013 if (bVerbose) {
1014 printf("\npress Enter ");
1015 getchar();
1016 }
1017 }
1018 else if (strncmp(string,"del",3)==0) {
1019 string+=3;
1020 if (parse_int(&string,&sel_nr)) {
1021 while(string[0]==' ')
1022 string++;
1023 if (string[0]=='-') {
1024 string++;
1025 parse_int(&string,&sel_nr2);
1026 } else
1027 sel_nr2=NOTSET;
1028 while(string[0]==' ')
1029 string++;
1030 if (string[0]=='\0')
1031 remove_group(sel_nr,sel_nr2,block,gn);
1032 else
1033 printf("\nSyntax error: \"%s\"\n",string);
1034 }
1035 }
1036 else if (strncmp(string,"keep",4)==0) {
1037 string+=4;
1038 if (parse_int(&string,&sel_nr)) {
1039 remove_group(sel_nr+1,block->nr-1,block,gn);
1040 remove_group(0,sel_nr-1,block,gn);
1041 }
1042 }
1043 else if (strncmp(string,"name",4)==0) {
1044 string+=4;
1045 if (parse_int(&string,&sel_nr)) {
1046 if ((sel_nr>=0) && (sel_nr<block->nr)) {
1047 sscanf(string,"%s",gname);
1048 sfree((*gn)[sel_nr]);
1049 (*gn)[sel_nr]=strdup(gname);
1050 }
1051 }
1052 }
1053 else if (strncmp(string,"case",4)==0) {
1054 bCase=!bCase;
1055 printf("Switched to case %s\n",bCase ? "sensitive" : "insensitive");
1056 }
1057 else if (string[0] == 'v' ) {
1058 bVerbose=!bVerbose;
1059 printf("Turned verbose %s\n",bVerbose ? "on" : "off");
1060 }
1061 else if (string[0] == 'l') {
1062 if ( check_have_atoms(atoms, ostring) )
1063 list_residues(atoms);
1064 }
1065 else if (strncmp(string,"splitch",7)==0) {
1066 string+=7;
1067 if ( check_have_atoms(atoms, ostring) &&
1068 parse_int(&string,&sel_nr) &&
1069 (sel_nr>=0) && (sel_nr<block->nr))
1070 split_chain(atoms,x,sel_nr,block,gn);
1071 }
1072 else if (strncmp(string,"splitres",8)==0 ) {
1073 string+=8;
1074 if ( check_have_atoms(atoms, ostring) &&
1075 parse_int(&string,&sel_nr) &&
1076 (sel_nr>=0) && (sel_nr<block->nr))
1077 split_group(atoms,sel_nr,block,gn,FALSE);
1078 }
1079 else if (strncmp(string,"splitat",7)==0 ) {
1080 string+=7;
1081 if ( check_have_atoms(atoms, ostring) &&
1082 parse_int(&string,&sel_nr) &&
1083 (sel_nr>=0) && (sel_nr<block->nr))
1084 split_group(atoms,sel_nr,block,gn,TRUE);
1085 }
1086 else if (string[0] == '\0') {
1087 bPrintOnce = TRUE;
1088 }
1089 else if (string[0] != 'q') {
1090 nr1=-1;
1091 nr2=-1;
1092 if (parse_entry(&string,natoms,atoms,block,gn,&nr,index,gname)) {
1093 do {
1094 while (string[0]==' ')
1095 string++;
1096
1097 bAnd=FALSE;
1098 bOr=FALSE;
1099 if (string[0]=='&')
1100 bAnd=TRUE;
1101 else if (string[0]=='|')
1102 bOr=TRUE;
1103
1104 if (bAnd || bOr) {
1105 string++;
1106 nr1=nr;
1107 for(i=0; i<nr; i++)
1108 index1[i]=index[i];
1109 strcpy(gname1,gname);
1110 if (parse_entry(&string,natoms,atoms,block,gn,&nr2,index2,gname2)) {
1111 if (bOr) {
1112 or_groups(nr1,index1,nr2,index2,&nr,index);
1113 sprintf(gname,"%s_%s",gname1,gname2);
1114 }
1115 else {
1116 and_groups(nr1,index1,nr2,index2,&nr,index);
1117 sprintf(gname,"%s_&_%s",gname1,gname2);
1118 }
1119 }
1120 }
1121 } while (bAnd || bOr);
1122 }
1123 while(string[0]==' ')
1124 string++;
1125 if (string[0])
1126 printf("\nSyntax error: \"%s\"\n",string);
1127 else if (nr>0) {
1128 copy2block(nr,index,block);
1129 srenew(*gn,block->nr);
1130 newgroup = block->nr-1;
1131 (*gn)[newgroup]=strdup(gname);
1132 }
1133 else
1134 printf("Group is empty\n");
1135 }
1136 } while (string[0]!='q');
1137
1138 sfree(index);
1139 sfree(index1);
1140 sfree(index2);
1141 }
1142
1143 static int block2natoms(t_blocka *block)
1144 {
1145 int i, natoms;
1146
1147 natoms = 0;
1148 for(i=0; i<block->nra; i++)
1149 natoms = max(natoms, block->a[i]);
1150 natoms++;
1151
1152 return natoms;
1153 }
1154
1155 void merge_blocks(t_blocka *dest, t_blocka *source)
1156 {
1157 int i,nra0,i0;
1158
1159 /* count groups, srenew and fill */
1160 i0 = dest->nr;
1161 nra0 = dest->nra;
1162 dest->nr+=source->nr;
1163 srenew(dest->index, dest->nr+1);
1164 for(i=0; i<source->nr; i++)
1165 dest->index[i0+i] = nra0 + source->index[i];
1166 /* count atoms, srenew and fill */
1167 dest->nra+=source->nra;
1168 srenew(dest->a, dest->nra);
1169 for(i=0; i<source->nra; i++)
1170 dest->a[nra0+i] = source->a[i];
1171
1172 /* terminate list */
1173 dest->index[dest->nr]=dest->nra;
1174
1175 }
1176
1177 int main(int argc,char *argv[])
1178 {
1179 const char *desc[] = {
1180 "Index groups are necessary for almost every gromacs program.",
1181 "All these programs can generate default index groups. You ONLY",
1182 "have to use make_ndx when you need SPECIAL index groups.",
1183 "There is a default index group for the whole system, 9 default",
1184 "index groups are generated for proteins, a default index group",
1185 "is generated for every other residue name.[PAR]"
1186 "When no index file is supplied, also make_ndx will generate the",
1187 "default groups.",
1188 "With the index editor you can select on atom, residue and chain names",
1189 "and numbers.",
1190 "When a run input file is supplied you can also select on atom type.",
1191 "You can use NOT, AND and OR, you can split groups",
1192 "into chains, residues or atoms. You can delete and rename groups.[PAR]",
1193 "The atom numbering in the editor and the index file starts at 1."
1194 };
1195
1196 static int natoms=0;
1197 static bool bVerbose=FALSE;
1198 t_pargs pa[] = {
1199 { "-natoms", FALSE, etINT, {&natoms},
1200 "set number of atoms (default: read from coordinate or index file)" },
1201 { "-verbose", FALSE, etBOOL, {&bVerbose},
1202 "HIDDENVerbose output" }
1203 };
1204 #define NPA asize(pa)
1205
1206 output_env_t oenv;
1207 char title[STRLEN];
1208 int nndxin;
1209 const char *stxfile;
1210 char **ndxinfiles;
1211 const char *ndxoutfile;
1212 bool bNatoms;
1213 int i,j;
1214 t_atoms *atoms;
1215 rvec *x,*v;
1216 int ePBC;
1217 matrix box;
1218 t_blocka *block,*block2;
1219 char **gnames,**gnames2;
1220 t_filenm fnm[] = {
1221 { efSTX, "-f", NULL, ffOPTRD },
1222 { efNDX, "-n", NULL, ffOPTRDMULT },
1223 { efNDX, "-o", NULL, ffWRITE }
1224 };
1225 #define NFILE asize(fnm)
1226
1227 CopyRight(stderr,argv[0]);
1228
1229 parse_common_args(&argc,argv,0,NFILE,fnm,NPA,pa,asize(desc),desc,
1230 0,NULL,&oenv);
1231
1232 stxfile = ftp2fn_null(efSTX,NFILE,fnm);
1233 if (opt2bSet("-n",NFILE,fnm)) {
1234 nndxin = opt2fns(&ndxinfiles,"-n",NFILE,fnm);
1235 } else {
1236 nndxin = 0;
1237 }
1238 ndxoutfile = opt2fn("-o",NFILE,fnm);
1239 bNatoms = opt2parg_bSet("-natoms",NPA,pa);
1240
1241 if (!stxfile && !nndxin)
1242 gmx_fatal(FARGS,"No input files (structure or index)");
1243
1244 if (stxfile) {
1245 snew(atoms,1);
1246 get_stx_coordnum(stxfile,&(atoms->nr));
1247 init_t_atoms(atoms,atoms->nr,TRUE);
1248 snew(x,atoms->nr);
1249 snew(v,atoms->nr);
1250 fprintf(stderr,"\nReading structure file\n");
1251 read_stx_conf(stxfile,title,atoms,x,v,&ePBC,box);
1252 natoms = atoms->nr;
1253 bNatoms=TRUE;
1254 } else {
1255 atoms = NULL;
1256 x = NULL;
1257 }
1258
1259 /* read input file(s) */
1260 block = new_blocka();
1261 gnames = NULL;
1262 printf("Going to read %d old index file(s)\n",nndxin);
1263 if (nndxin) {
1264 for(i=0; i<nndxin; i++) {
1265 block2 = init_index(ndxinfiles[i],&gnames2);
1266 srenew(gnames, block->nr+block2->nr);
1267 for(j=0; j<block2->nr; j++)
1268 gnames[block->nr+j]=gnames2[j];
1269 sfree(gnames2);
1270 merge_blocks(block, block2);
1271 sfree(block2->a);
1272 sfree(block2->index);
1273 /* done_block(block2); */
1274 sfree(block2);
1275 }
1276 }
1277 else {
1278 snew(gnames,1);
1279 analyse(atoms,block,&gnames,FALSE,TRUE);
1280 }
1281
1282 if (!bNatoms) {
1283 natoms = block2natoms(block);
1284 printf("Counted atom numbers up to %d in index file\n", natoms);
1285 }
1286
1287 edit_index(natoms,atoms,x,block,&gnames,bVerbose);
1288
1289 write_index(ndxoutfile,block,gnames);
1290
1291 thanx(stderr);
1292
1293 return 0;
1294 }
Atempt to write QMMM interface for Gamess-US