src/gromacs/gmxlib/index.c

   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 <assert.h>
  42 #include "sysstuff.h"
  43 #include "strdb.h"
  44 #include "futil.h"
  45 #include "macros.h"
  46 #include "names.h"
  47 #include "string2.h"
  48 #include "statutil.h"
  49 #include "confio.h"
  50 #include "main.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 gmx_bool gmx_ask_yesno(gmx_bool bASK)
  73 {
  74     char c;
  75
  76     if (bASK)
  77     {
  78         do
  79         {
  80             c = toupper(fgetc(stdin));
  81         }
  82         while ((c != 'Y') && (c != 'N'));
  83
  84         return (c == 'Y');
  85     }
  86     else
  87     {
  88         return FALSE;
  89     }
  90 }
  91
  92 t_blocka *new_blocka(void)
  93 {
  94     t_blocka *block;
  95
  96     snew(block, 1);
  97     snew(block->index, 1);
  98
  99     return block;
 100 }
 101
 102 void write_index(const char *outf, t_blocka *b, char **gnames)
 103 {
 104     FILE *out;
 105     int   i, j, k;
 106
 107     out = gmx_fio_fopen(outf, "w");
 108     /* fprintf(out,"%5d  %5d\n",b->nr,b->nra); */
 109     for (i = 0; (i < b->nr); i++)
 110     {
 111         fprintf(out, "[ %s ]\n", gnames[i]);
 112         for (k = 0, j = b->index[i]; j < b->index[i+1]; j++, k++)
 113         {
 114             fprintf(out, "%4d ", b->a[j]+1);
 115             if ((k % 15) == 14)
 116             {
 117                 fprintf(out, "\n");
 118             }
 119         }
 120         fprintf(out, "\n");
 121     }
 122     gmx_fio_fclose(out);
 123 }
 124
 125 void add_grp(t_blocka *b, char ***gnames, int nra, atom_id a[], const char *name)
 126 {
 127     int i;
 128
 129     srenew(b->index, b->nr+2);
 130     srenew(*gnames, b->nr+1);
 131     (*gnames)[b->nr] = strdup(name);
 132
 133     srenew(b->a, b->nra+nra);
 134     for (i = 0; (i < nra); i++)
 135     {
 136         b->a[b->nra++] = a[i];
 137     }
 138     b->nr++;
 139     b->index[b->nr] = b->nra;
 140 }
 141
 142 /* compare index in `a' with group in `b' at `index',
 143    when `index'<0 it is relative to end of `b' */
 144 static gmx_bool grp_cmp(t_blocka *b, int nra, atom_id a[], int index)
 145 {
 146     int i;
 147
 148     if (index < 0)
 149     {
 150         index = b->nr-1+index;
 151     }
 152     if (index >= b->nr)
 153     {
 154         gmx_fatal(FARGS, "no such index group %d in t_blocka (nr=%d)", index, b->nr);
 155     }
 156     /* compare sizes */
 157     if (nra != b->index[index+1] - b->index[index])
 158     {
 159         return FALSE;
 160     }
 161     for (i = 0; i < nra; i++)
 162     {
 163         if (a[i] != b->a[b->index[index]+i])
 164         {
 165             return FALSE;
 166         }
 167     }
 168     return TRUE;
 169 }
 170
 171 static void
 172 p_status(const char **restype, int nres, const char **typenames, int ntypes)
 173 {
 174     int   i, j;
 175     int   found;
 176
 177     int * counter;
 178
 179     snew(counter, ntypes);
 180     for (i = 0; i < ntypes; i++)
 181     {
 182         counter[i] = 0;
 183     }
 184     for (i = 0; i < nres; i++)
 185     {
 186         found = 0;
 187         for (j = 0; j < ntypes; j++)
 188         {
 189             if (!gmx_strcasecmp(restype[i], typenames[j]))
 190             {
 191                 counter[j]++;
 192             }
 193         }
 194     }
 195
 196     for (i = 0; (i < ntypes); i++)
 197     {
 198         if (counter[i] > 0)
 199         {
 200             printf("There are: %5d %10s residues\n", counter[i], typenames[i]);
 201         }
 202     }
 203
 204     sfree(counter);
 205 }
 206
 207
 208 atom_id *
 209 mk_aid(t_atoms *atoms, const char ** restype, const char * typestring, int *nra, gmx_bool bMatch)
 210 /* Make an array of atom_ids for all atoms with residuetypes matching typestring, or the opposite if bMatch is false */
 211 {
 212     atom_id *a;
 213     int      i;
 214     int      res;
 215
 216     snew(a, atoms->nr);
 217     *nra = 0;
 218     for (i = 0; (i < atoms->nr); i++)
 219     {
 220         res = !gmx_strcasecmp(restype[atoms->atom[i].resind], typestring);
 221         if (bMatch == FALSE)
 222         {
 223             res = !res;
 224         }
 225         if (res)
 226         {
 227             a[(*nra)++] = i;
 228         }
 229     }
 230
 231     return a;
 232 }
 233
 234 typedef struct {
 235     char    *rname;
 236     gmx_bool bNeg;
 237     char    *gname;
 238 } restp_t;
 239
 240 static void analyse_other(const char ** restype, t_atoms *atoms,
 241                           t_blocka *gb, char ***gn, gmx_bool bASK, gmx_bool bVerb)
 242 {
 243     restp_t *restp = NULL;
 244     char   **attp  = NULL;
 245     char    *rname, *aname;
 246     atom_id *other_ndx, *aid, *aaid;
 247     int      i, j, k, l, resind, naid, naaid, natp, nrestp = 0;
 248
 249     for (i = 0; (i < atoms->nres); i++)
 250     {
 251         if (gmx_strcasecmp(restype[i], "Protein") && gmx_strcasecmp(restype[i], "DNA") && gmx_strcasecmp(restype[i], "RNA") && gmx_strcasecmp(restype[i], "Water"))
 252         {
 253             break;
 254         }
 255     }
 256     if (i < atoms->nres)
 257     {
 258         /* we have others */
 259         if (bVerb)
 260         {
 261             printf("Analysing residues not classified as Protein/DNA/RNA/Water and splitting into groups...\n");
 262         }
 263         snew(other_ndx, atoms->nr);
 264         for (k = 0; (k < atoms->nr); k++)
 265         {
 266             resind = atoms->atom[k].resind;
 267             rname  = *atoms->resinfo[resind].name;
 268             if (gmx_strcasecmp(restype[resind], "Protein") && gmx_strcasecmp(restype[resind], "DNA") &&
 269                 gmx_strcasecmp(restype[resind], "RNA") && gmx_strcasecmp(restype[resind], "Water"))
 270             {
 271
 272                 for (l = 0; (l < nrestp); l++)
 273                 {
 274                     if (strcmp(restp[l].rname, rname) == 0)
 275                     {
 276                         break;
 277                     }
 278                 }
 279                 if (l == nrestp)
 280                 {
 281                     srenew(restp, nrestp+1);
 282                     restp[nrestp].rname = strdup(rname);
 283                     restp[nrestp].bNeg  = FALSE;
 284                     restp[nrestp].gname = strdup(rname);
 285                     nrestp++;
 286
 287                 }
 288             }
 289         }
 290         for (i = 0; (i < nrestp); i++)
 291         {
 292             snew(aid, atoms->nr);
 293             naid = 0;
 294             for (j = 0; (j < atoms->nr); j++)
 295             {
 296                 rname = *atoms->resinfo[atoms->atom[j].resind].name;
 297                 if ((strcmp(restp[i].rname, rname) == 0 && !restp[i].bNeg) ||
 298                     (strcmp(restp[i].rname, rname) != 0 &&  restp[i].bNeg))
 299                 {
 300                     aid[naid++] = j;
 301                 }
 302             }
 303             add_grp(gb, gn, naid, aid, restp[i].gname);
 304             if (bASK)
 305             {
 306                 printf("split %s into atoms (y/n) ? ", restp[i].gname);
 307                 fflush(stdout);
 308                 if (gmx_ask_yesno(bASK))
 309                 {
 310                     natp = 0;
 311                     for (k = 0; (k < naid); k++)
 312                     {
 313                         aname = *atoms->atomname[aid[k]];
 314                         for (l = 0; (l < natp); l++)
 315                         {
 316                             if (strcmp(aname, attp[l]) == 0)
 317                             {
 318                                 break;
 319                             }
 320                         }
 321                         if (l == natp)
 322                         {
 323                             srenew(attp, ++natp);
 324                             attp[natp-1] = aname;
 325                         }
 326                     }
 327                     if (natp > 1)
 328                     {
 329                         for (l = 0; (l < natp); l++)
 330                         {
 331                             snew(aaid, naid);
 332                             naaid = 0;
 333                             for (k = 0; (k < naid); k++)
 334                             {
 335                                 aname = *atoms->atomname[aid[k]];
 336                                 if (strcmp(aname, attp[l]) == 0)
 337                                 {
 338                                     aaid[naaid++] = aid[k];
 339                                 }
 340                             }
 341                             add_grp(gb, gn, naaid, aaid, attp[l]);
 342                             sfree(aaid);
 343                         }
 344                     }
 345                     sfree(attp);
 346                     attp = NULL;
 347                 }
 348                 sfree(aid);
 349             }
 350         }
 351         sfree(other_ndx);
 352     }
 353 }
 354
 355 /*! /brief Instances of this struct contain the data necessary to
 356  *         construct a single (protein) index group in
 357  *         analyse_prot(). */
 358 typedef struct gmx_help_make_index_group
 359 {
 360     /** The set of atom names that will be used to form this index group */
 361     const char **defining_atomnames;
 362     /** Size of the defining_atomnames array */
 363     const int    num_defining_atomnames;
 364     /** Name of this index group */
 365     const char  *group_name;
 366     /** Whether the above atom names name the atoms in the group, or
 367         those not in the group */
 368     gmx_bool bTakeComplement;
 369     /** The index in wholename gives the first item in the arrays of
 370        atomnames that should be tested with 'gmx_strncasecmp' in stead of
 371        gmx_strcasecmp, or -1 if all items should be tested with strcasecmp
 372        This is comparable to using a '*' wildcard at the end of specific
 373        atom names, but that is more involved to implement...
 374      */
 375     int wholename;
 376     /** Only create this index group if it differs from the one specified in compareto,
 377        where -1 means to always create this group. */
 378     int compareto;
 379 } t_gmx_help_make_index_group;
 380
 381 static void analyse_prot(const char ** restype, t_atoms *atoms,
 382                          t_blocka *gb, char ***gn, gmx_bool bASK, gmx_bool bVerb)
 383 {
 384     /* lists of atomnames to be used in constructing index groups: */
 385     static const char *pnoh[]    = { "H", "HN" };
 386     static const char *pnodum[]  = {
 387         "MN1",  "MN2",  "MCB1", "MCB2", "MCG1", "MCG2",
 388         "MCD1", "MCD2", "MCE1", "MCE2", "MNZ1", "MNZ2"
 389     };
 390     static const char *calpha[]  = { "CA" };
 391     static const char *bb[]      = { "N", "CA", "C" };
 392     static const char *mc[]      = { "N", "CA", "C", "O", "O1", "O2", "OC1", "OC2", "OT", "OXT" };
 393     static const char *mcb[]     = { "N", "CA", "CB", "C", "O", "O1", "O2", "OC1", "OC2", "OT", "OXT" };
 394     static const char *mch[]     = {
 395         "N", "CA", "C", "O", "O1", "O2", "OC1", "OC2", "OT", "OXT",
 396         "H1", "H2", "H3", "H", "HN"
 397     };
 398
 399     static const t_gmx_help_make_index_group constructing_data[] =
 400     {{ NULL,   0, "Protein",      TRUE,  -1, -1},
 401      { pnoh,   asize(pnoh),   "Protein-H",    TRUE,  0,  -1},
 402      { calpha, asize(calpha), "C-alpha",      FALSE, -1, -1},
 403      { bb,     asize(bb),     "Backbone",     FALSE, -1, -1},
 404      { mc,     asize(mc),     "MainChain",    FALSE, -1, -1},
 405      { mcb,    asize(mcb),    "MainChain+Cb", FALSE, -1, -1},
 406      { mch,    asize(mch),    "MainChain+H",  FALSE, -1, -1},
 407      { mch,    asize(mch),    "SideChain",    TRUE,  -1, -1},
 408      { mch,    asize(mch),    "SideChain-H",  TRUE,  11, -1},
 409      { pnodum, asize(pnodum), "Prot-Masses",  TRUE,  -1, 0}, };
 410     const int   num_index_groups = asize(constructing_data);
 411
 412     int         n, j;
 413     atom_id    *aid;
 414     int         nra, nnpres, npres;
 415     gmx_bool    match;
 416     char        ndx_name[STRLEN], *atnm;
 417     int         i;
 418
 419     if (bVerb)
 420     {
 421         printf("Analysing Protein...\n");
 422     }
 423     snew(aid, atoms->nr);
 424
 425     /* calculate the number of protein residues */
 426     npres = 0;
 427     for (i = 0; (i < atoms->nres); i++)
 428     {
 429         if (0 == gmx_strcasecmp(restype[i], "Protein"))
 430         {
 431             npres++;
 432         }
 433     }
 434     /* find matching or complement atoms */
 435     for (i = 0; (i < (int)num_index_groups); i++)
 436     {
 437         nra = 0;
 438         for (n = 0; (n < atoms->nr); n++)
 439         {
 440             if (0 == gmx_strcasecmp(restype[atoms->atom[n].resind], "Protein"))
 441             {
 442                 match = FALSE;
 443                 for (j = 0; (j < constructing_data[i].num_defining_atomnames); j++)
 444                 {
 445                     /* skip digits at beginning of atomname, e.g. 1H */
 446                     atnm = *atoms->atomname[n];
 447                     while (isdigit(atnm[0]))
 448                     {
 449                         atnm++;
 450                     }
 451                     if ( (constructing_data[i].wholename == -1) || (j < constructing_data[i].wholename) )
 452                     {
 453                         if (0 == gmx_strcasecmp(constructing_data[i].defining_atomnames[j], atnm))
 454                         {
 455                             match = TRUE;
 456                         }
 457                     }
 458                     else
 459                     {
 460                         if (0 == gmx_strncasecmp(constructing_data[i].defining_atomnames[j], atnm, strlen(constructing_data[i].defining_atomnames[j])))
 461                         {
 462                             match = TRUE;
 463                         }
 464                     }
 465                 }
 466                 if (constructing_data[i].bTakeComplement != match)
 467                 {
 468                     aid[nra++] = n;
 469                 }
 470             }
 471         }
 472         /* if we want to add this group always or it differs from previous
 473            group, add it: */
 474         if (-1 == constructing_data[i].compareto || !grp_cmp(gb, nra, aid, constructing_data[i].compareto-i) )
 475         {
 476             add_grp(gb, gn, nra, aid, constructing_data[i].group_name);
 477         }
 478     }
 479
 480     if (bASK)
 481     {
 482         for (i = 0; (i < (int)num_index_groups); i++)
 483         {
 484             printf("Split %12s into %5d residues (y/n) ? ", constructing_data[i].group_name, npres);
 485             if (gmx_ask_yesno(bASK))
 486             {
 487                 int resind;
 488                 nra = 0;
 489                 for (n = 0; ((atoms->atom[n].resind < npres) && (n < atoms->nr)); )
 490                 {
 491                     resind = atoms->atom[n].resind;
 492                     for (; ((atoms->atom[n].resind == resind) && (n < atoms->nr)); n++)
 493                     {
 494                         match = FALSE;
 495                         for (j = 0; (j < constructing_data[i].num_defining_atomnames); j++)
 496                         {
 497                             if (0 == gmx_strcasecmp(constructing_data[i].defining_atomnames[j], *atoms->atomname[n]))
 498                             {
 499                                 match = TRUE;
 500                             }
 501                         }
 502                         if (constructing_data[i].bTakeComplement != match)
 503                         {
 504                             aid[nra++] = n;
 505                         }
 506                     }
 507                     /* copy the residuename to the tail of the groupname */
 508                     if (nra > 0)
 509                     {
 510                         t_resinfo *ri;
 511                         ri = &atoms->resinfo[resind];
 512                         sprintf(ndx_name, "%s_%s%d%c",
 513                                 constructing_data[i].group_name, *ri->name, ri->nr, ri->ic == ' ' ? '\0' : ri->ic);
 514                         add_grp(gb, gn, nra, aid, ndx_name);
 515                         nra = 0;
 516                     }
 517                 }
 518             }
 519         }
 520         printf("Make group with sidechain and C=O swapped (y/n) ? ");
 521         if (gmx_ask_yesno(bASK))
 522         {
 523             /* Make swap sidechain C=O index */
 524             int resind, hold;
 525             nra = 0;
 526             for (n = 0; ((atoms->atom[n].resind < npres) && (n < atoms->nr)); )
 527             {
 528                 resind = atoms->atom[n].resind;
 529                 hold   = -1;
 530                 for (; ((atoms->atom[n].resind == resind) && (n < atoms->nr)); n++)
 531                 {
 532                     if (strcmp("CA", *atoms->atomname[n]) == 0)
 533                     {
 534                         aid[nra++] = n;
 535                         hold       = nra;
 536                         nra       += 2;
 537                     }
 538                     else if (strcmp("C", *atoms->atomname[n]) == 0)
 539                     {
 540                         if (hold == -1)
 541                         {
 542                             gmx_incons("Atom naming problem");
 543                         }
 544                         aid[hold] = n;
 545                     }
 546                     else if (strcmp("O", *atoms->atomname[n]) == 0)
 547                     {
 548                         if (hold == -1)
 549                         {
 550                             gmx_incons("Atom naming problem");
 551                         }
 552                         aid[hold+1] = n;
 553                     }
 554                     else if (strcmp("O1", *atoms->atomname[n]) == 0)
 555                     {
 556                         if (hold == -1)
 557                         {
 558                             gmx_incons("Atom naming problem");
 559                         }
 560                         aid[hold+1] = n;
 561                     }
 562                     else
 563                     {
 564                         aid[nra++] = n;
 565                     }
 566                 }
 567             }
 568             /* copy the residuename to the tail of the groupname */
 569             if (nra > 0)
 570             {
 571                 add_grp(gb, gn, nra, aid, "SwapSC-CO");
 572                 nra = 0;
 573             }
 574         }
 575     }
 576     sfree(aid);
 577 }
 578
 579
 580
 581
 582 /* Return 0 if the name was found, otherwise -1.
 583  * p_restype is set to a pointer to the type name, or 'Other' if we did not find it.
 584  */
 585 int
 586 gmx_residuetype_get_type(gmx_residuetype_t rt, const char * resname, const char ** p_restype)
 587 {
 588     int    i, rc;
 589
 590     rc = -1;
 591     for (i = 0; i < rt->n && rc; i++)
 592     {
 593         rc = gmx_strcasecmp(rt->resname[i], resname);
 594     }
 595
 596     *p_restype = (rc == 0) ? rt->restype[i-1] : gmx_residuetype_undefined;
 597
 598     return rc;
 599 }
 600
 601 int
 602 gmx_residuetype_add(gmx_residuetype_t rt, const char *newresname, const char *newrestype)
 603 {
 604     int           i;
 605     int           found;
 606     const char *  p_oldtype;
 607
 608     found = !gmx_residuetype_get_type(rt, newresname, &p_oldtype);
 609
 610     if (found && gmx_strcasecmp(p_oldtype, newrestype))
 611     {
 612         fprintf(stderr, "Warning: Residue '%s' already present with type '%s' in database, ignoring new type '%s'.",
 613                 newresname, p_oldtype, newrestype);
 614     }
 615
 616     if (found == 0)
 617     {
 618         srenew(rt->resname, rt->n+1);
 619         srenew(rt->restype, rt->n+1);
 620         rt->resname[rt->n] = strdup(newresname);
 621         rt->restype[rt->n] = strdup(newrestype);
 622         rt->n++;
 623     }
 624
 625     return 0;
 626 }
 627
 628
 629 int
 630 gmx_residuetype_init(gmx_residuetype_t *prt)
 631 {
 632     FILE                 *  db;
 633     char                    line[STRLEN];
 634     char                    resname[STRLEN], restype[STRLEN], dum[STRLEN];
 635     char                 *  p;
 636     int                     i;
 637     struct gmx_residuetype *rt;
 638
 639     snew(rt, 1);
 640     *prt = rt;
 641
 642     rt->n        = 0;
 643     rt->resname  = NULL;
 644     rt->restype  = NULL;
 645
 646     db = libopen("residuetypes.dat");
 647
 648     while (get_a_line(db, line, STRLEN))
 649     {
 650         strip_comment(line);
 651         trim(line);
 652         if (line[0] != '\0')
 653         {
 654             if (sscanf(line, "%s %s %s", resname, restype, dum) != 2)
 655             {
 656                 gmx_fatal(FARGS, "Incorrect number of columns (2 expected) for line in residuetypes.dat");
 657             }
 658             gmx_residuetype_add(rt, resname, restype);
 659         }
 660     }
 661
 662     fclose(db);
 663
 664     return 0;
 665 }
 666
 667
 668
 669 int
 670 gmx_residuetype_destroy(gmx_residuetype_t rt)
 671 {
 672     int i;
 673
 674     for (i = 0; i < rt->n; i++)
 675     {
 676         sfree(rt->resname[i]);
 677         sfree(rt->restype[i]);
 678     }
 679     sfree(rt->resname);
 680     sfree(rt->restype);
 681     sfree(rt);
 682
 683     return 0;
 684 }
 685
 686 int
 687 gmx_residuetype_get_alltypes(gmx_residuetype_t    rt,
 688                              const char ***       p_typenames,
 689                              int *                ntypes)
 690 {
 691     int            i, j, n;
 692     int            found;
 693     const char **  my_typename;
 694     char       *   p;
 695
 696     n = 0;
 697
 698     my_typename = NULL;
 699     for (i = 0; i < rt->n; i++)
 700     {
 701         p     = rt->restype[i];
 702         found = 0;
 703         for (j = 0; j < n && !found; j++)
 704         {
 705             found = !gmx_strcasecmp(p, my_typename[j]);
 706         }
 707
 708         if (!found)
 709         {
 710             srenew(my_typename, n+1);
 711             my_typename[n] = p;
 712             n++;
 713         }
 714     }
 715     *ntypes      = n;
 716     *p_typenames = my_typename;
 717
 718     return 0;
 719 }
 720
 721
 722
 723 gmx_bool
 724 gmx_residuetype_is_protein(gmx_residuetype_t rt, const char *resnm)
 725 {
 726     gmx_bool    rc;
 727     const char *p_type;
 728
 729     if (gmx_residuetype_get_type(rt, resnm, &p_type) == 0 &&
 730         gmx_strcasecmp(p_type, "Protein") == 0)
 731     {
 732         rc = TRUE;
 733     }
 734     else
 735     {
 736         rc = FALSE;
 737     }
 738     return rc;
 739 }
 740
 741 gmx_bool
 742 gmx_residuetype_is_dna(gmx_residuetype_t rt, const char *resnm)
 743 {
 744     gmx_bool    rc;
 745     const char *p_type;
 746
 747     if (gmx_residuetype_get_type(rt, resnm, &p_type) == 0 &&
 748         gmx_strcasecmp(p_type, "DNA") == 0)
 749     {
 750         rc = TRUE;
 751     }
 752     else
 753     {
 754         rc = FALSE;
 755     }
 756     return rc;
 757 }
 758
 759 gmx_bool
 760 gmx_residuetype_is_rna(gmx_residuetype_t rt, const char *resnm)
 761 {
 762     gmx_bool    rc;
 763     const char *p_type;
 764
 765     if (gmx_residuetype_get_type(rt, resnm, &p_type) == 0 &&
 766         gmx_strcasecmp(p_type, "RNA") == 0)
 767     {
 768         rc = TRUE;
 769     }
 770     else
 771     {
 772         rc = FALSE;
 773     }
 774     return rc;
 775 }
 776
 777 /* Return the size of the arrays */
 778 int
 779 gmx_residuetype_get_size(gmx_residuetype_t rt)
 780 {
 781     return rt->n;
 782 }
 783
 784 /* Search for a residuetype with name resnm within the
 785  * gmx_residuetype database. Return the index if found,
 786  * otherwise -1.
 787  */
 788 int
 789 gmx_residuetype_get_index(gmx_residuetype_t rt, const char *resnm)
 790 {
 791     int i, rc;
 792
 793     rc = -1;
 794     for (i = 0; i < rt->n && rc; i++)
 795     {
 796         rc = gmx_strcasecmp(rt->resname[i], resnm);
 797     }
 798
 799     return (0 == rc) ? i-1 : -1;
 800 }
 801
 802 /* Return the name of the residuetype with the given index, or
 803  * NULL if not found. */
 804 const char *
 805 gmx_residuetype_get_name(gmx_residuetype_t rt, int index)
 806 {
 807     if (index >= 0 && index < rt->n)
 808     {
 809         return rt->resname[index];
 810     }
 811     else
 812     {
 813         return NULL;
 814     }
 815 }
 816
 817
 818
 819 void analyse(t_atoms *atoms, t_blocka *gb, char ***gn, gmx_bool bASK, gmx_bool bVerb)
 820 {
 821     gmx_residuetype_t rt = NULL;
 822     char             *resnm;
 823     atom_id          *aid;
 824     const char    **  restype;
 825     int               nra;
 826     int               i, k;
 827     size_t            j;
 828     int               ntypes;
 829     char           *  p;
 830     const char     ** p_typename;
 831     int               iwater, iion;
 832     int               nwater, nion;
 833     int               found;
 834
 835     if (bVerb)
 836     {
 837         printf("Analysing residue names:\n");
 838     }
 839     /* Create system group, every single atom */
 840     snew(aid, atoms->nr);
 841     for (i = 0; i < atoms->nr; i++)
 842     {
 843         aid[i] = i;
 844     }
 845     add_grp(gb, gn, atoms->nr, aid, "System");
 846     sfree(aid);
 847
 848     /* For every residue, get a pointer to the residue type name */
 849     gmx_residuetype_init(&rt);
 850     assert(rt);
 851
 852     snew(restype, atoms->nres);
 853     ntypes     = 0;
 854     p_typename = NULL;
 855     for (i = 0; i < atoms->nres; i++)
 856     {
 857         resnm = *atoms->resinfo[i].name;
 858         gmx_residuetype_get_type(rt, resnm, &(restype[i]));
 859
 860         /* Note that this does not lead to a N*N loop, but N*K, where
 861          * K is the number of residue _types_, which is small and independent of N.
 862          */
 863         found = 0;
 864         for (k = 0; k < ntypes && !found; k++)
 865         {
 866             found = !strcmp(restype[i], p_typename[k]);
 867         }
 868         if (!found)
 869         {
 870             srenew(p_typename, ntypes+1);
 871             p_typename[ntypes++] = strdup(restype[i]);
 872         }
 873     }
 874
 875     if (bVerb)
 876     {
 877         p_status(restype, atoms->nres, p_typename, ntypes);
 878     }
 879
 880     for (k = 0; k < ntypes; k++)
 881     {
 882         aid = mk_aid(atoms, restype, p_typename[k], &nra, TRUE);
 883
 884         /* Check for special types to do fancy stuff with */
 885
 886         if (!gmx_strcasecmp(p_typename[k], "Protein") && nra > 0)
 887         {
 888             sfree(aid);
 889             /* PROTEIN */
 890             analyse_prot(restype, atoms, gb, gn, bASK, bVerb);
 891
 892             /* Create a Non-Protein group */
 893             aid = mk_aid(atoms, restype, "Protein", &nra, FALSE);
 894             if ((nra > 0) && (nra < atoms->nr))
 895             {
 896                 add_grp(gb, gn, nra, aid, "non-Protein");
 897             }
 898             sfree(aid);
 899         }
 900         else if (!gmx_strcasecmp(p_typename[k], "Water") && nra > 0)
 901         {
 902             add_grp(gb, gn, nra, aid, p_typename[k]);
 903             /* Add this group as 'SOL' too, for backward compatibility with older gromacs versions */
 904             add_grp(gb, gn, nra, aid, "SOL");
 905
 906             sfree(aid);
 907
 908             /* Solvent, create a negated group too */
 909             aid = mk_aid(atoms, restype, "Water", &nra, FALSE);
 910             if ((nra > 0) && (nra < atoms->nr))
 911             {
 912                 add_grp(gb, gn, nra, aid, "non-Water");
 913             }
 914             sfree(aid);
 915         }
 916         else if (nra > 0)
 917         {
 918             /* Other groups */
 919             add_grp(gb, gn, nra, aid, p_typename[k]);
 920             sfree(aid);
 921             analyse_other(restype, atoms, gb, gn, bASK, bVerb);
 922         }
 923     }
 924
 925     sfree(p_typename);
 926     sfree(restype);
 927     gmx_residuetype_destroy(rt);
 928
 929     /* Create a merged water_and_ions group */
 930     iwater = -1;
 931     iion   = -1;
 932     nwater = 0;
 933     nion   = 0;
 934
 935     for (i = 0; i < gb->nr; i++)
 936     {
 937         if (!gmx_strcasecmp((*gn)[i], "Water"))
 938         {
 939             iwater = i;
 940             nwater = gb->index[i+1]-gb->index[i];
 941         }
 942         else if (!gmx_strcasecmp((*gn)[i], "Ion"))
 943         {
 944             iion = i;
 945             nion = gb->index[i+1]-gb->index[i];
 946         }
 947     }
 948
 949     if (nwater > 0 && nion > 0)
 950     {
 951         srenew(gb->index, gb->nr+2);
 952         srenew(*gn, gb->nr+1);
 953         (*gn)[gb->nr] = strdup("Water_and_ions");
 954         srenew(gb->a, gb->nra+nwater+nion);
 955         if (nwater > 0)
 956         {
 957             for (i = gb->index[iwater]; i < gb->index[iwater+1]; i++)
 958             {
 959                 gb->a[gb->nra++] = gb->a[i];
 960             }
 961         }
 962         if (nion > 0)
 963         {
 964             for (i = gb->index[iion]; i < gb->index[iion+1]; i++)
 965             {
 966                 gb->a[gb->nra++] = gb->a[i];
 967             }
 968         }
 969         gb->nr++;
 970         gb->index[gb->nr] = gb->nra;
 971     }
 972 }
 973
 974
 975 void check_index(char *gname, int n, atom_id index[], char *traj, int natoms)
 976 {
 977     int i;
 978
 979     for (i = 0; i < n; i++)
 980     {
 981         if (index[i] >= natoms)
 982         {
 983             gmx_fatal(FARGS, "%s atom number (index[%d]=%d) is larger than the number of atoms in %s (%d)",
 984                       gname ? gname : "Index", i+1, index[i]+1,
 985                       traj ? traj : "the trajectory", natoms);
 986         }
 987         else if (index[i] < 0)
 988         {
 989             gmx_fatal(FARGS, "%s atom number (index[%d]=%d) is less than zero",
 990                       gname ? gname : "Index", i+1, index[i]+1);
 991         }
 992     }
 993 }
 994
 995 t_blocka *init_index(const char *gfile, char ***grpname)
 996 {
 997     FILE      *in;
 998     t_blocka  *b;
 999     int        a, maxentries;
1000     int        i, j, ng, nread;
1001     char       line[STRLEN], *pt, str[STRLEN];
1002
1003     in = gmx_fio_fopen(gfile, "r");
1004     snew(b, 1);
1005     get_a_line(in, line, STRLEN);
1006     if (line[0] == '[')
1007     {
1008         /* new format */
1009         b->nr      = 0;
1010         b->index   = NULL;
1011         b->nra     = 0;
1012         b->a       = NULL;
1013         *grpname   = NULL;
1014         maxentries = 0;
1015         do
1016         {
1017             if (get_header(line, str))
1018             {
1019                 b->nr++;
1020                 srenew(b->index, b->nr+1);
1021                 srenew(*grpname, b->nr);
1022                 if (b->nr == 1)
1023                 {
1024                     b->index[0] = 0;
1025                 }
1026                 b->index[b->nr]     = b->index[b->nr-1];
1027                 (*grpname)[b->nr-1] = strdup(str);
1028             }
1029             else
1030             {
1031                 if (b->nr == 0)
1032                 {
1033                     gmx_fatal(FARGS, "The first header of your indexfile is invalid");
1034                 }
1035                 pt = line;
1036                 while (sscanf(pt, "%s", str) == 1)
1037                 {
1038                     i = b->index[b->nr];
1039                     if (i >= maxentries)
1040                     {
1041                         maxentries += 1024;
1042                         srenew(b->a, maxentries);
1043                     }
1044                     b->a[i] = strtol(str, NULL, 10)-1;
1045                     b->index[b->nr]++;
1046                     (b->nra)++;
1047                     pt = strstr(pt, str)+strlen(str);
1048                 }
1049             }
1050         }
1051         while (get_a_line(in, line, STRLEN));
1052     }
1053     else
1054     {
1055         /* old format */
1056         sscanf(line, "%d%d", &b->nr, &b->nra);
1057         snew(b->index, b->nr+1);
1058         snew(*grpname, b->nr);
1059         b->index[0] = 0;
1060         snew(b->a, b->nra);
1061         for (i = 0; (i < b->nr); i++)
1062         {
1063             nread         = fscanf(in, "%s%d", str, &ng);
1064             (*grpname)[i] = strdup(str);
1065             b->index[i+1] = b->index[i]+ng;
1066             if (b->index[i+1] > b->nra)
1067             {
1068                 gmx_fatal(FARGS, "Something wrong in your indexfile at group %s", str);
1069             }
1070             for (j = 0; (j < ng); j++)
1071             {
1072                 nread               = fscanf(in, "%d", &a);
1073                 b->a[b->index[i]+j] = a;
1074             }
1075         }
1076     }
1077     gmx_fio_fclose(in);
1078
1079     for (i = 0; (i < b->nr); i++)
1080     {
1081         for (j = b->index[i]; (j < b->index[i+1]); j++)
1082         {
1083             if (b->a[j] < 0)
1084             {
1085                 fprintf(stderr, "\nWARNING: negative index %d in group %s\n\n",
1086                         b->a[j], (*grpname)[i]);
1087             }
1088         }
1089     }
1090
1091     return b;
1092 }
1093
1094 static void minstring(char *str)
1095 {
1096     int i;
1097
1098     for (i = 0; (i < (int)strlen(str)); i++)
1099     {
1100         if (str[i] == '-')
1101         {
1102             str[i] = '_';
1103         }
1104     }
1105 }
1106
1107 int find_group(char s[], int ngrps, char **grpname)
1108 {
1109     int      aa, i, n;
1110     char     string[STRLEN];
1111     gmx_bool bMultiple;
1112
1113     bMultiple = FALSE;
1114     n         = strlen(s);
1115     aa        = NOTSET;
1116     /* first look for whole name match */
1117     if (aa == NOTSET)
1118     {
1119         for (i = 0; i < ngrps; i++)
1120         {
1121             if (gmx_strcasecmp_min(s, grpname[i]) == 0)
1122             {
1123                 if (aa != NOTSET)
1124                 {
1125                     bMultiple = TRUE;
1126                 }
1127                 aa = i;
1128             }
1129         }
1130     }
1131     /* second look for first string match */
1132     if (aa == NOTSET)
1133     {
1134         for (i = 0; i < ngrps; i++)
1135         {
1136             if (gmx_strncasecmp_min(s, grpname[i], n) == 0)
1137             {
1138                 if (aa != NOTSET)
1139                 {
1140                     bMultiple = TRUE;
1141                 }
1142                 aa = i;
1143             }
1144         }
1145     }
1146     /* last look for arbitrary substring match */
1147     if (aa == NOTSET)
1148     {
1149         upstring(s);
1150         minstring(s);
1151         for (i = 0; i < ngrps; i++)
1152         {
1153             strcpy(string, grpname[i]);
1154             upstring(string);
1155             minstring(string);
1156             if (strstr(string, s) != NULL)
1157             {
1158                 if (aa != NOTSET)
1159                 {
1160                     bMultiple = TRUE;
1161                 }
1162                 aa = i;
1163             }
1164         }
1165     }
1166     if (bMultiple)
1167     {
1168         printf("Error: Multiple groups '%s' selected\n", s);
1169         aa = NOTSET;
1170     }
1171     return aa;
1172 }
1173
1174 static int qgroup(int *a, int ngrps, char **grpname)
1175 {
1176     char     s[STRLEN];
1177     int      aa;
1178     gmx_bool bInRange;
1179     char    *end;
1180
1181     do
1182     {
1183         fprintf(stderr, "Select a group: ");
1184         do
1185         {
1186             if (scanf("%s", s) != 1)
1187             {
1188                 gmx_fatal(FARGS, "Cannot read from input");
1189             }
1190             trim(s); /* remove spaces */
1191         }
1192         while (strlen(s) == 0);
1193         aa = strtol(s, &end, 10);
1194         if (aa == 0 && end[0] != '\0') /* string entered */
1195         {
1196             aa = find_group(s, ngrps, grpname);
1197         }
1198         bInRange = (aa >= 0 && aa < ngrps);
1199         if (!bInRange)
1200         {
1201             printf("Error: No such group '%s'\n", s);
1202         }
1203     }
1204     while (!bInRange);
1205     printf("Selected %d: '%s'\n", aa, grpname[aa]);
1206     *a = aa;
1207     return aa;
1208 }
1209
1210 static void rd_groups(t_blocka *grps, char **grpname, char *gnames[],
1211                       int ngrps, int isize[], atom_id *index[], int grpnr[])
1212 {
1213     int i, j, gnr1;
1214
1215     if (grps->nr == 0)
1216     {
1217         gmx_fatal(FARGS, "Error: no groups in indexfile");
1218     }
1219     for (i = 0; (i < grps->nr); i++)
1220     {
1221         fprintf(stderr, "Group %5d (%15s) has %5d elements\n", i, grpname[i],
1222                 grps->index[i+1]-grps->index[i]);
1223     }
1224     for (i = 0; (i < ngrps); i++)
1225     {
1226         if (grps->nr > 1)
1227         {
1228             do
1229             {
1230                 gnr1 = qgroup(&grpnr[i], grps->nr, grpname);
1231                 if ((gnr1 < 0) || (gnr1 >= grps->nr))
1232                 {
1233                     fprintf(stderr, "Select between %d and %d.\n", 0, grps->nr-1);
1234                 }
1235             }
1236             while ((gnr1 < 0) || (gnr1 >= grps->nr));
1237         }
1238         else
1239         {
1240             fprintf(stderr, "There is one group in the index\n");
1241             gnr1 = 0;
1242         }
1243         gnames[i] = strdup(grpname[gnr1]);
1244         isize[i]  = grps->index[gnr1+1]-grps->index[gnr1];
1245         snew(index[i], isize[i]);
1246         for (j = 0; (j < isize[i]); j++)
1247         {
1248             index[i][j] = grps->a[grps->index[gnr1]+j];
1249         }
1250     }
1251 }
1252
1253 void rd_index(const char *statfile, int ngrps, int isize[],
1254               atom_id *index[], char *grpnames[])
1255 {
1256     char    **gnames;
1257     t_blocka *grps;
1258     int      *grpnr;
1259
1260     snew(grpnr, ngrps);
1261     if (!statfile)
1262     {
1263         gmx_fatal(FARGS, "No index file specified");
1264     }
1265     grps = init_index(statfile, &gnames);
1266     rd_groups(grps, gnames, grpnames, ngrps, isize, index, grpnr);
1267 }
1268
1269 void rd_index_nrs(char *statfile, int ngrps, int isize[],
1270                   atom_id *index[], char *grpnames[], int grpnr[])
1271 {
1272     char    **gnames;
1273     t_blocka *grps;
1274
1275     if (!statfile)
1276     {
1277         gmx_fatal(FARGS, "No index file specified");
1278     }
1279     grps = init_index(statfile, &gnames);
1280
1281     rd_groups(grps, gnames, grpnames, ngrps, isize, index, grpnr);
1282 }
1283
1284 void get_index(t_atoms *atoms, const char *fnm, int ngrps,
1285                int isize[], atom_id *index[], char *grpnames[])
1286 {
1287     char    ***gnames;
1288     t_blocka  *grps = NULL;
1289     int       *grpnr;
1290
1291     snew(grpnr, ngrps);
1292     snew(gnames, 1);
1293     if (fnm != NULL)
1294     {
1295         grps = init_index(fnm, gnames);
1296     }
1297     else if (atoms)
1298     {
1299         snew(grps, 1);
1300         snew(grps->index, 1);
1301         analyse(atoms, grps, gnames, FALSE, FALSE);
1302     }
1303     else
1304     {
1305         gmx_incons("You need to supply a valid atoms structure or a valid index file name");
1306     }
1307
1308     rd_groups(grps, *gnames, grpnames, ngrps, isize, index, grpnr);
1309 }
1310
1311 t_cluster_ndx *cluster_index(FILE *fplog, const char *ndx)
1312 {
1313     t_cluster_ndx *c;
1314     int            i;
1315
1316     snew(c, 1);
1317     c->clust     = init_index(ndx, &c->grpname);
1318     c->maxframe  = -1;
1319     for (i = 0; (i < c->clust->nra); i++)
1320     {
1321         c->maxframe = max(c->maxframe, c->clust->a[i]);
1322     }
1323     fprintf(fplog ? fplog : stdout,
1324             "There are %d clusters containing %d structures, highest framenr is %d\n",
1325             c->clust->nr, c->clust->nra, c->maxframe);
1326     if (debug)
1327     {
1328         pr_blocka(debug, 0, "clust", c->clust, TRUE);
1329         for (i = 0; (i < c->clust->nra); i++)
1330         {
1331             if ((c->clust->a[i] < 0) || (c->clust->a[i] > c->maxframe))
1332             {
1333                 gmx_fatal(FARGS, "Range check error for c->clust->a[%d] = %d\n"
1334                           "should be within 0 and %d", i, c->clust->a[i], c->maxframe+1);
1335             }
1336         }
1337     }
1338     c->inv_clust = make_invblocka(c->clust, c->maxframe);
1339
1340     return c;
1341 }