Project

General

Profile

index.c

Fixed index.c - Kyle Beauchamp, 08/22/2010 10:19 PM

 
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
 * GROningen Mixture of Alchemy and Childrens' Stories
34
 */
35
#ifdef HAVE_CONFIG_H
36
#include <config.h>
37
#endif
38

    
39
#include <ctype.h>
40
#include <string.h>
41
#include "sysstuff.h"
42
#include "strdb.h"
43
#include "futil.h"
44
#include "macros.h"
45
#include "names.h"
46
#include "string2.h"
47
#include "statutil.h"
48
#include "confio.h"
49
#include "main.h"
50
#include "copyrite.h"
51
#include "typedefs.h"
52
#include "smalloc.h"
53
#include "invblock.h"
54
#include "macros.h"
55
#include "index.h"
56
#include "txtdump.h"
57
#include "gmxfio.h"
58

    
59

    
60

    
61
const char gmx_residuetype_undefined[]="Other";
62

    
63
struct gmx_residuetype
64
{
65
    int      n; 
66
    char **  resname;
67
    char **  restype;
68
    
69
};
70

    
71

    
72
static bool gmx_ask_yesno(bool bASK)
73
{
74
  char c;
75

    
76
  if (bASK) {
77
    do {
78
      c=toupper(fgetc(stdin));
79
    } while ((c != 'Y') && (c != 'N'));
80

    
81
    return (c == 'Y');
82
  }
83
  else
84
    return FALSE;
85
}
86

    
87
t_blocka *new_blocka(void)
88
{
89
  t_blocka *block;
90

    
91
  snew(block,1);
92
  snew(block->index,1);
93

    
94
  return block;
95
}
96

    
97
void write_index(const char *outf, t_blocka *b,char **gnames)
98
{
99
  FILE *out;
100
  int  i,j,k;
101

    
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);
108
      if ((k % 15) == 14)
109
        fprintf(out,"\n");
110
    }
111
    fprintf(out,"\n");
112
  }
113
  gmx_fio_fclose(out);
114
}
115

    
116
void add_grp(t_blocka *b,char ***gnames,int nra,atom_id a[],const char *name)
117
{
118
  int i;
119

    
120
  srenew(b->index,b->nr+2);
121
  srenew(*gnames,b->nr+1);
122
  (*gnames)[b->nr]=strdup(name);
123

    
124
  srenew(b->a,b->nra+nra);
125
  for(i=0; (i<nra); i++)
126
    b->a[b->nra++]=a[i];
127
  b->nr++;
128
  b->index[b->nr]=b->nra;
129
}
130

    
131
/* compare index in `a' with group in `b' at `index', 
132
   when `index'<0 it is relative to end of `b' */
133
static bool grp_cmp(t_blocka *b, int nra, atom_id a[], int index)
134
{
135
  int i;
136
  
137
  if (index < 0)
138
    index = b->nr-1+index;
139
  if (index >= b->nr)
140
    gmx_fatal(FARGS,"no such index group %d in t_blocka (nr=%d)",index,b->nr);
141
  /* compare sizes */
142
  if ( nra != b->index[index+1] - b->index[index] )
143
    return FALSE;
144
  for(i=0; i<nra; i++)
145
    if ( a[i] != b->a[b->index[index]+i] )
146
      return FALSE;
147
  return TRUE;
148
}
149

    
150
static void 
151
p_status(const char **restype, int nres, const char **typenames, int ntypes, bool bVerb)
152
{
153
    int i,j;
154
    int found;
155
    
156
    int * counter;
157
    
158
    snew(counter,ntypes);
159
    for(i=0;i<ntypes;i++)
160
    {
161
        counter[i]=0;
162
    }
163
    for(i=0;i<nres;i++)
164
    {
165
        found=0;
166
        for(j=0;j<ntypes;j++)
167
        {
168
            if(!gmx_strcasecmp(restype[i],typenames[j]))
169
            {
170
                counter[j]++;
171
            }
172
        }
173
    }
174
    
175
    if (bVerb)
176
    {
177
        for(i=0; (i<ntypes); i++) 
178
        {
179
            if(counter[i]>0)
180
            {
181
                printf("There are: %5d %10s residues\n",counter[i],typenames[i]);
182
            }
183
        }
184
    }
185
}
186

    
187

    
188
atom_id *
189
mk_aid(t_atoms *atoms,const char ** restype,const char * typestring,int *nra,bool bMatch)
190
/* Make an array of atom_ids for all atoms with residuetypes matching typestring, or the opposite if bMatch is false */
191
{
192
    atom_id *a;
193
    int     i;
194
    int     res;
195
    
196
    snew(a,atoms->nr);
197
    *nra=0;
198
    for(i=0; (i<atoms->nr); i++) 
199
    {
200
        res=!gmx_strcasecmp(restype[atoms->atom[i].resind],typestring);
201
        if(bMatch==FALSE)
202
        {
203
            res=!res;
204
        }
205
        if(res)
206
        {
207
            a[(*nra)++]=i;
208
        }
209
    }
210
  
211
    return a;
212
}
213

    
214
typedef struct {
215
  char *rname;
216
  bool bNeg;
217
  char *gname;
218
} restp_t;
219

    
220
static void analyse_other(const char ** restype,t_atoms *atoms,
221
                          t_blocka *gb,char ***gn,bool bASK,bool bVerb)
222
{
223
  restp_t *restp=NULL;
224
  char **attp=NULL;
225
  char *rname,*aname;
226
  atom_id *other_ndx,*aid,*aaid;
227
  int  i,j,k,l,resind,naid,naaid,natp,nrestp=0;
228
  
229
    for(i=0; (i<atoms->nres); i++)
230
    {
231
        if (gmx_strcasecmp(restype[i],"Protein") && gmx_strcasecmp(restype[i],"DNA") && gmx_strcasecmp(restype[i],"RNA") && gmx_strcasecmp(restype[i],"Water"))
232
        {
233
            break;
234
        }
235
    }
236
  if (i < atoms->nres) {
237
    /* we have others */
238
    if (bVerb)
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")) 
246
        {
247

    
248
        for(l=0; (l<nrestp); l++)
249
          if (strcmp(restp[l].rname,rname) == 0)
250
            break;
251
        if (l==nrestp) {
252
          srenew(restp,nrestp+1);
253
          restp[nrestp].rname = strdup(rname);
254
          restp[nrestp].bNeg  = FALSE;
255
          restp[nrestp].gname = strdup(rname);
256
          nrestp++;
257
        
258
    }
259
      }
260
    }
261
    for(i=0; (i<nrestp); i++) {
262
      snew(aid,atoms->nr);
263
      naid=0;
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)) {
268
          aid[naid++] = j;
269
        }
270
      }
271
      add_grp(gb,gn,naid,aid,restp[i].gname);
272
      if (bASK) {
273
        printf("split %s into atoms (y/n) ? ",restp[i].gname);
274
        fflush(stdout);
275
        if (gmx_ask_yesno(bASK)) {
276
          natp=0;
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)
281
                break;
282
            if (l == natp) {
283
              srenew(attp,++natp);
284
              attp[natp-1]=aname;
285
            }
286
          }
287
          if (natp > 1) {
288
            for(l=0; (l<natp); l++) {
289
              snew(aaid,naid);
290
              naaid=0;
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];
295
              }
296
              add_grp(gb,gn,naaid,aaid,attp[l]);
297
              sfree(aaid);
298
            }
299
          }
300
          sfree(attp);
301
          attp=NULL;
302
        }
303
        sfree(aid);
304
      }
305
    }
306
    sfree(other_ndx);
307
  }
308
}
309

    
310
static void analyse_prot(const char ** restype,t_atoms *atoms,
311
                         t_blocka *gb,char ***gn,bool bASK,bool bVerb)
312
{
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)
330
  };
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", 
335
    "Prot-Masses"
336
  };
337
  /* construct index group containing (TRUE) or excluding (FALSE)
338
     given atom names */
339
  const bool complement[NCH] = { 
340
    TRUE, TRUE, FALSE, FALSE, FALSE, FALSE, FALSE, TRUE, TRUE, TRUE
341
  };
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 'strncasecmp' in stead of
345
   * 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...
348
   */
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 };
352

    
353
  int     n,j;
354
  atom_id *aid;
355
  int     nra,nnpres,npres;
356
  bool    match;
357
  char    ndx_name[STRLEN],*atnm;
358
  int i;
359

    
360
  if (bVerb)
361
  {
362
    printf("Analysing Protein...\n");
363
  }
364
  snew(aid,atoms->nr);
365

    
366
  /* calculate the number of protein residues */
367
  npres=0;
368
  for(i=0; (i<atoms->nres); i++)
369
  if (!gmx_strcasecmp(restype[i],"Protein"))
370
  {
371
      npres++;
372
  }
373
  /* find matching or complement atoms */
374
  for(i=0; (i<(int)NCH); i++) {
375
    nra=0;
376
    for(n=0; (n<atoms->nr); n++) {
377
        if (!gmx_strcasecmp(restype[atoms->atom[n].resind],"Protein")) {
378

    
379
        match=FALSE;
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]))
384
            atnm++;
385
          if ( (wholename[i]==-1) || (j<wholename[i]) ) {
386
            if (gmx_strcasecmp(chains[i][j],atnm) == 0)
387
              match=TRUE;
388
          } else {
389
            if (gmx_strncasecmp(chains[i][j],atnm,strlen(chains[i][j])) == 0)
390
              match=TRUE;
391
          }
392
        }
393
        if (match != complement[i])
394
          aid[nra++]=n;
395
      }
396
    }
397
    /* if we want to add this group always or it differs from previous 
398
       group, add it: */
399
    if ( compareto[i] == -1 || !grp_cmp(gb,nra,aid,compareto[i]-i) )
400
      add_grp(gb,gn,nra,aid,ch_name[i]); 
401
  }
402
  
403
  if (bASK) {
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)) {
407
        int resind;
408
        nra = 0;
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++) {
412
            match=FALSE;
413
            for(j=0;(j<sizes[i]); j++) 
414
              if (gmx_strcasecmp(chains[i][j],*atoms->atomname[n]) == 0)
415
                match=TRUE;
416
            if (match != complement[i])
417
              aid[nra++]=n;
418
          }
419
          /* copy the residuename to the tail of the groupname */
420
          if (nra > 0) {
421
            t_resinfo *ri;
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);
426
            nra = 0;
427
          }
428
        }
429
      } 
430
    }
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 */
434
      int resind,hold;
435
      nra = 0;
436
      for(n=0;((atoms->atom[n].resind < npres) && (n<atoms->nr));) {
437
        resind = atoms->atom[n].resind;
438
        hold  = -1;
439
        for(;((atoms->atom[n].resind==resind) && (n<atoms->nr));n++)
440
          if (strcmp("CA",*atoms->atomname[n]) == 0) {
441
            aid[nra++]=n;
442
            hold=nra;
443
            nra+=2;
444
          } else if (strcmp("C",*atoms->atomname[n]) == 0) {
445
            if (hold == -1)
446
              gmx_incons("Atom naming problem");
447
            aid[hold]=n;
448
          } else if (strcmp("O",*atoms->atomname[n]) == 0) {
449
            if (hold == -1)
450
              gmx_incons("Atom naming problem");
451
            aid[hold+1]=n;
452
          } else if (strcmp("O1",*atoms->atomname[n]) == 0) {
453
            if (hold == -1)
454
              gmx_incons("Atom naming problem");
455
            aid[hold+1]=n;
456
          } else 
457
            aid[nra++]=n;
458
      }
459
      /* copy the residuename to the tail of the groupname */
460
      if (nra > 0) {
461
        add_grp(gb,gn,nra,aid,"SwapSC-CO");
462
        nra = 0;
463
      } 
464
    }
465
  }
466
  sfree(aid);
467
}
468

    
469

    
470

    
471

    
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.
474
 */
475
int
476
gmx_residuetype_get_type(gmx_residuetype_t rt,const char * resname, const char ** p_restype)
477
{
478
    int    i,rc;
479
    
480
    rc=-1;
481
    for(i=0;i<rt->n && rc;i++)
482
    {
483
        rc=gmx_strcasecmp(rt->resname[i],resname);
484
    }
485
    
486
    *p_restype = (rc==0) ? rt->restype[i-1] : gmx_residuetype_undefined;
487
    
488
    return rc;
489
}
490

    
491
int
492
gmx_residuetype_add(gmx_residuetype_t rt,const char *newresname, const char *newrestype)
493
{
494
    int     i;
495
    int     found;
496
    const char *  p_oldtype;
497
    
498
    found = !gmx_residuetype_get_type(rt,newresname,&p_oldtype);
499
    
500
    if(found && gmx_strcasecmp(p_oldtype,newrestype))
501
    {
502
        fprintf(stderr,"Warning: Residue '%s' already present with type '%s' in database, ignoring new type '%s'.",
503
                newresname,p_oldtype,newrestype);
504
    }
505
    
506
    if(found==0)
507
    {
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);
512
        rt->n++;
513
    }
514
  
515
    return 0;
516
}
517

    
518

    
519
int
520
gmx_residuetype_init(gmx_residuetype_t *prt)
521
{
522
    FILE *  db;
523
    char    line[STRLEN];
524
    char    resname[STRLEN],restype[STRLEN],dum[STRLEN];
525
    char *  p;
526
    int     i;
527
    struct gmx_residuetype *rt;
528
    
529
    snew(rt,1);
530
    *prt=rt;
531
    
532
    rt->n        = 0;
533
    rt->resname  = NULL;
534
    rt->restype = NULL;
535
    
536
    db=libopen("residuetypes.dat");
537
    
538
    while(get_a_line(db,line,STRLEN)) 
539
    {
540
        strip_comment(line);
541
        trim(line);
542
        if(line[0]!='\0')
543
        {
544
            if(sscanf(line,"%s %s %s",resname,restype,dum)!=2)
545
            {
546
                gmx_fatal(FARGS,"Incorrect number of columns (2 expected) for line in residuetypes.dat");
547
            }
548
            gmx_residuetype_add(rt,resname,restype);
549
        }
550
    }
551
    
552
    fclose(db);
553
    
554
    return 0;
555
}
556

    
557

    
558

    
559
int
560
gmx_residuetype_destroy(gmx_residuetype_t rt)
561
{
562
    int i;
563
    
564
    for(i=0;i<rt->n;i++)
565
    {
566
        free(rt->resname[i]);
567
        free(rt->restype[i]);
568
    }
569
    free(rt);
570
    
571
    return 0;
572
}
573

    
574
int
575
gmx_residuetype_get_alltypes(gmx_residuetype_t    rt,
576
                             const char ***       p_typenames,
577
                             int *                ntypes)
578
{
579
    int      i,j,n;
580
    int      found;
581
    const char **  my_typename;
582
    char *   p;
583
    
584
    n=0;
585
    
586
    my_typename=NULL;
587
    for(i=0;i<rt->n;i++)
588
    {
589
        p=rt->restype[i];
590
        found=0;
591
        for(j=0;j<n && !found;j++)
592
        {
593
            found=!gmx_strcasecmp(p,my_typename[j]);
594
        }
595
        
596
        if(!found)
597
        {
598
            srenew(my_typename,n+1);
599
            my_typename[n]=p;
600
            n++;
601
        }
602
    }
603
    *ntypes=n;
604
    *p_typenames=my_typename; 
605
    
606
    return 0;
607
}
608
    
609

    
610

    
611
bool 
612
gmx_residuetype_is_protein(gmx_residuetype_t rt, const char *resnm)
613
{
614
    bool rc;
615
    const char *p_type;
616
    
617
    if(gmx_residuetype_get_type(rt,resnm,&p_type)==0 &&
618
       gmx_strcasecmp(p_type,"Protein")==0)
619
    {
620
        rc=TRUE;
621
    }
622
    else
623
    {
624
        rc=FALSE;
625
    }
626
    return rc;
627
}
628

    
629
bool 
630
gmx_residuetype_is_dna(gmx_residuetype_t rt, const char *resnm)
631
{
632
    bool rc;
633
    const char *p_type;
634

    
635
    if(gmx_residuetype_get_type(rt,resnm,&p_type)==0 &&
636
       gmx_strcasecmp(p_type,"DNA")==0)
637
    {
638
        rc=TRUE;
639
    }
640
    else
641
    {
642
        rc=FALSE;
643
    }
644
    return rc;
645
}
646

    
647
bool 
648
gmx_residuetype_is_rna(gmx_residuetype_t rt, const char *resnm)
649
{
650
    bool rc;
651
    const char *p_type;
652

    
653
    if(gmx_residuetype_get_type(rt,resnm,&p_type)==0 &&
654
       gmx_strcasecmp(p_type,"RNA")==0)
655
    {
656
        rc=TRUE;
657
    }
658
    else
659
    {
660
        rc=FALSE;
661
    }
662
    return rc;
663
}
664

    
665

    
666

    
667

    
668
void analyse(t_atoms *atoms,t_blocka *gb,char ***gn,bool bASK,bool bVerb)
669
{
670
    gmx_residuetype_t rt;
671
    char    *resnm;
672
    atom_id *aid;
673
    const char **  restype;
674
    int     nra;
675
    int     i,k;
676
    size_t  j;
677
    int     ntypes;
678
    char *  p;
679
    const char ** p_typename;
680
    int     iwater,iion;
681
    int     nwater,nion;
682
    int     found;
683
    
684
    if (bVerb)
685
    {
686
        printf("Analysing residue names:\n");
687
    }
688
    /* Create system group, every single atom */
689
    snew(aid,atoms->nr);
690
    for(i=0;i<atoms->nr;i++)
691
    {
692
        aid[i]=i;
693
    }
694
    add_grp(gb,gn,atoms->nr,aid,"System"); 
695
    sfree(aid);
696

    
697
    /* For every residue, get a pointer to the residue type name */
698
    gmx_residuetype_init(&rt);
699

    
700
    snew(restype,atoms->nres);
701
    ntypes = 0;
702
    p_typename = NULL;
703
    for(i=0;i<atoms->nres;i++)
704
    {
705
        resnm = *atoms->resinfo[i].name;
706
        gmx_residuetype_get_type(rt,resnm,&(restype[i]));
707

    
708
        if(i==0)
709
        {
710
            snew(p_typename,1);
711
            p_typename[ntypes++] = strdup(restype[i]);
712
        }
713
        else
714
        {
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.
717
             */
718
            found = 0;
719
            for(k=0;k<i && !found;k++)
720
            {
721
                found = !strcmp(restype[i],restype[k]);
722
            }
723
            if(!found)
724
            {
725
                srenew(p_typename,ntypes+1);
726
                p_typename[ntypes++] = strdup(restype[i]);
727
            }
728
        }
729
    }    
730
    
731
    p_status(restype,atoms->nres,p_typename,ntypes,bVerb);
732

    
733
    for(k=0;k<ntypes;k++)
734
    {              
735
        aid=mk_aid(atoms,restype,p_typename[k],&nra,TRUE);
736

    
737
        /* Check for special types to do fancy stuff with */
738
        
739
        if(!gmx_strcasecmp(p_typename[k],"Protein") && nra>0)
740
        {
741
            sfree(aid);
742
            /* PROTEIN */
743
            analyse_prot(restype,atoms,gb,gn,bASK,bVerb);
744
            
745
            /* Create a Non-Protein group */
746
            aid=mk_aid(atoms,restype,"Protein",&nra,FALSE);
747
            if ((nra > 0) && (nra < atoms->nr))
748
            {
749
                add_grp(gb,gn,nra,aid,"non-Protein"); 
750
            }
751
            sfree(aid);
752
        }
753
        else if(!gmx_strcasecmp(p_typename[k],"Water") && nra>0)
754
        {
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"); 
758

    
759
            sfree(aid);
760

    
761
            /* Solvent, create a negated group too */
762
            aid=mk_aid(atoms,restype,"Water",&nra,FALSE);
763
            if ((nra > 0) && (nra < atoms->nr))
764
            {
765
                add_grp(gb,gn,nra,aid,"non-Water"); 
766
            }
767
            sfree(aid);
768
        }
769
        else if(nra>0)
770
        {
771
            /* Other groups */
772
            add_grp(gb,gn,nra,aid,p_typename[k]); 
773
            sfree(aid);
774
            analyse_other(restype,atoms,gb,gn,bASK,bVerb);
775
        }
776
    }
777
    
778
    sfree(p_typename);
779
    sfree(restype);
780
    gmx_residuetype_destroy(rt);      
781
    
782
    /* Create a merged water_and_ions group */
783
    iwater = -1;
784
    iion   = -1;
785
    nwater = 0;
786
    nion   = 0;
787
        
788
    for(i=0;i<gb->nr;i++)
789
    {        
790
        if(!gmx_strcasecmp((*gn)[i],"Water"))
791
        {
792
            iwater = i;
793
            nwater = gb->index[i+1]-gb->index[i];
794
        }
795
        else if(!gmx_strcasecmp((*gn)[i],"Ion"))
796
        {
797
            iion = i;
798
            nion = gb->index[i+1]-gb->index[i];
799
        }
800
    }
801
    
802
    if(nwater>0 && nion>0)
803
    {
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);
808
        if(nwater>0)
809
        {
810
            for(i=gb->index[iwater];i<gb->index[iwater+1];i++)
811
            {
812
                gb->a[gb->nra++] = gb->a[i];
813
            }
814
        }
815
        if(nion>0)
816
        {
817
            for(i=gb->index[iion];i<gb->index[iion+1];i++)
818
            {
819
                gb->a[gb->nra++] = gb->a[i];
820
            }
821
        }
822
        gb->nr++;
823
        gb->index[gb->nr]=gb->nra;
824
    }
825
}
826

    
827

    
828
void check_index(char *gname,int n,atom_id index[],char *traj,int natoms)
829
{
830
  int i;
831
  
832
  for(i=0; i<n; i++)
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);
840
}
841

    
842
t_blocka *init_index(const char *gfile, char ***grpname)
843
{
844
  FILE     *in;
845
  t_blocka  *b;
846
  int      a,maxentries;
847
  int      i,j,ng,nread;
848
  char     line[STRLEN],*pt,str[STRLEN];
849

    
850
  in=gmx_fio_fopen(gfile,"r");
851
  snew(b,1);
852
  get_a_line(in,line,STRLEN);
853
  if ( line[0]=='[' ) {
854
    /* new format */
855
    b->nr=0;
856
    b->index=NULL;
857
    b->nra=0;
858
    b->a=NULL;
859
    *grpname=NULL;
860
    maxentries=0;
861
    do {
862
      if (get_header(line,str)) {
863
        b->nr++;
864
        srenew(b->index,b->nr+1);
865
        srenew(*grpname,b->nr);
866
        if (b->nr==1)
867
          b->index[0]=0;
868
        b->index[b->nr]=b->index[b->nr-1];
869
        (*grpname)[b->nr-1]=strdup(str);
870
      } else {
871
        pt=line;
872
        while ((i=sscanf(pt,"%s",str)) == 1) {
873
          i=b->index[b->nr];
874
          if (i>=maxentries) {
875
            maxentries+=1024;
876
            srenew(b->a,maxentries);
877
          }
878
          b->a[i]=strtol(str, NULL, 10)-1;
879
          b->index[b->nr]++;
880
          (b->nra)++;
881
          pt=strstr(pt,str)+strlen(str);
882
        }
883
      }
884
    } while (get_a_line(in,line,STRLEN));
885
  } 
886
  else {
887
    /* old format */
888
    sscanf(line,"%d%d",&b->nr,&b->nra);
889
    snew(b->index,b->nr+1);
890
    snew(*grpname,b->nr);
891
    b->index[0]=0;
892
    snew(b->a,b->nra);
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;
902
      }
903
    }
904
  }
905
  gmx_fio_fclose(in);
906

    
907
  for(i=0; (i<b->nr); i++) {
908
    for(j=b->index[i]; (j<b->index[i+1]); j++) {
909
      if (b->a[j] < 0) 
910
        fprintf(stderr,"\nWARNING: negative index %d in group %s\n\n",
911
                b->a[j],(*grpname)[i]);
912
    }
913
  }
914
  
915
  return b;
916
}
917

    
918
static void minstring(char *str)
919
{
920
  int i;
921

    
922
  for (i=0; (i < (int)strlen(str)); i++) 
923
    if (str[i]=='-')
924
      str[i]='_';
925
}
926

    
927
int find_group(char s[], int ngrps, char **grpname)
928
{
929
  int aa, i, n;
930
  char string[STRLEN];
931
  bool bMultiple;
932
  
933
  bMultiple = FALSE;
934
  n = strlen(s);
935
  aa=NOTSET;
936
  /* first look for whole name match */
937
  if (aa==NOTSET)
938
    for(i=0; i<ngrps; i++)
939
      if (strcasecmp_min(s,grpname[i])==0) {
940
        if(aa!=NOTSET)
941
          bMultiple = TRUE;
942
        aa=i;
943
      }
944
  /* second look for first string match */
945
  if (aa==NOTSET)
946
    for(i=0; i<ngrps; i++)
947
      if (strncasecmp_min(s,grpname[i],n)==0) {
948
        if(aa!=NOTSET)
949
          bMultiple = TRUE;
950
        aa=i;
951
      }
952
  /* last look for arbitrary substring match */
953
  if (aa==NOTSET) {
954
    upstring(s);
955
    minstring(s);
956
    for(i=0; i<ngrps; i++) {
957
      strcpy(string, grpname[i]);
958
      upstring(string);
959
      minstring(string);
960
      if (strstr(string,s)!=NULL) {
961
        if(aa!=NOTSET)
962
          bMultiple = TRUE;
963
        aa=i;
964
      }
965
    }
966
  }
967
  if (bMultiple) {
968
    printf("Error: Multiple groups '%s' selected\n", s);
969
    aa=NOTSET;
970
  }
971
  return aa;
972
}
973

    
974
static int qgroup(int *a, int ngrps, char **grpname)
975
{
976
    char s[STRLEN];
977
    int  aa;
978
    bool bInRange;
979
    char *end;
980

    
981
    do {
982
        fprintf(stderr,"Select a group: ");
983
        do {
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);
992
        if (!bInRange)
993
            printf("Error: No such group '%s'\n", s);
994
    } while (!bInRange);
995
    printf("Selected %d: '%s'\n", aa, grpname[aa]);
996
    *a = aa;
997
    return aa;
998
}
999

    
1000
static void rd_groups(t_blocka *grps,char **grpname,char *gnames[],
1001
                      int ngrps,int isize[],atom_id *index[],int grpnr[])
1002
{
1003
  int i,j,gnr1;
1004

    
1005
  if (grps->nr==0)
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++) {
1011
    if (grps->nr > 1)
1012
      do {
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));
1017
    else {
1018
      fprintf(stderr,"There is one group in the index\n");
1019
      gnr1=0;
1020
    }
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];
1026
  }
1027
}
1028

    
1029
void rd_index(const char *statfile,int ngrps,int isize[],
1030
              atom_id *index[],char *grpnames[])
1031
{
1032
  char    **gnames;
1033
  t_blocka *grps;
1034
  int     *grpnr;
1035
  
1036
  snew(grpnr,ngrps);
1037
  if (!statfile)
1038
    gmx_fatal(FARGS,"No index file specified");
1039
  grps=init_index(statfile,&gnames);
1040
  rd_groups(grps,gnames,grpnames,ngrps,isize,index,grpnr);
1041
}
1042

    
1043
void rd_index_nrs(char *statfile,int ngrps,int isize[],
1044
                  atom_id *index[],char *grpnames[],int grpnr[])
1045
{
1046
  char    **gnames;
1047
  t_blocka *grps;
1048
  
1049
  if (!statfile)
1050
    gmx_fatal(FARGS,"No index file specified");
1051
  grps=init_index(statfile,&gnames);
1052
  
1053
  rd_groups(grps,gnames,grpnames,ngrps,isize,index,grpnr);
1054
}
1055

    
1056
void get_index(t_atoms *atoms, const char *fnm, int ngrps,
1057
               int isize[], atom_id *index[],char *grpnames[])
1058
{
1059
  char    ***gnames;
1060
  t_blocka *grps = NULL; 
1061
  int     *grpnr;
1062
  
1063
  snew(grpnr,ngrps);
1064
  snew(gnames,1);
1065
  if (fnm != NULL) {
1066
    grps=init_index(fnm,gnames);
1067
  }
1068
  else if (atoms) {
1069
    snew(grps,1);
1070
    snew(grps->index,1);
1071
    analyse(atoms,grps,gnames,FALSE,FALSE);
1072
  }
1073
  else 
1074
    gmx_incons("You need to supply a valid atoms structure or a valid index file name");
1075
  
1076
  rd_groups(grps,*gnames,grpnames,ngrps,isize,index,grpnr);
1077
}
1078

    
1079
t_cluster_ndx *cluster_index(FILE *fplog,const char *ndx)
1080
{
1081
  t_cluster_ndx *c;
1082
  int i;
1083
  
1084
  snew(c,1);
1085
  c->clust     = init_index(ndx,&c->grpname);
1086
  c->maxframe = -1;
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);
1092
  if (debug) {
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);
1098
  }
1099
  c->inv_clust=make_invblocka(c->clust,c->maxframe);
1100
          
1101
  return c;
1102
}
1103