molecule.java

A program to find frequent molecular substructures and discriminative fragments in a database of mol

            bb = a.bonds[k];    /* traverse the bonds of the atom */
            if (((bb.flags & cur) != 0) && (bb != b)) break;
          }                     /* find the next ring bond */
          if (++n < 6) ring[n] = bb;
          b = bb; b.mark = 0;   /* record up to six bonds of a ring */
          a = (b.src != a) ? b.src : b.dst;
        } while (a != ring[0].src);
        if (++n != 6) continue; /* process only rings with 6 bonds */
        k = ring[1].type & Bond.TYPEMASK;
        if (((k != Bond.SINGLE) && (k != Bond.DOUBLE))
        ||  (ring[0].type == ring[1].type)
        ||  (ring[0].type != ring[2].type)
        ||  (ring[2].type != ring[4].type)
        ||  (ring[1].type != ring[3].type)
        ||  (ring[3].type != ring[5].type))
          continue;             /* check for a Kekul\'e represent. */
        cnt++;                  /* count the converted ring */
        for (k = 6; --k >= 0; ) {
          b = ring[k];          /* traverse the ring found */
          b.type      = Bond.AROMATIC;
          b.src.type |= Atom.AROMATIC;
          b.dst.type |= Atom.AROMATIC;
        }                       /* remove the old bond types and */
      }                         /* set all bond types to aromatic, */
    }                           /* also mark all atoms as aromatic */
    this.markRings(0, 0);       /* unmark the rings again */
    return cnt;                 /* return the number of rings */
  }  /* aromatize() */

  /*--------------------------------------------------------------------
  The above function turns all Kekul\'e representations of aromatic
  rings into true aromatic rings (that is, it replaces the alternating
  single and double bonds with aromatic bonds).
  --------------------------------------------------------------------*/

  public boolean hasOpenRings (int min, int max)
  {                             /* --- check for open rings */
    int  i;                     /* loop variable */
    Bond b;                     /* to traverse the bonds */

    this.markRings(min, max, 0);/* mark rings, but keep type flag */
    for (i = this.bondcnt; --i >= 0; ) {
      b = this.bonds[i];        /* traverse all bonds */
      if (((b.flags & Bond.RINGS)    == 0)
      &&  ((b.type  & Bond.RINGBOND) != 0))
        return true;            /* check for a bond with ring flag, */
    }                           /* but without individual rings */
    return false;               /* if there is, there is an open ring */
  }  /* hasOpenRings() */       /* otherwise all rings are closed */

  /*------------------------------------------------------------------*/

  protected void maskTypes (int masks[])
  {                             /* --- mask atom and bond types */
    int  i, k;                  /* loop variable, index */
    Atom a;                     /* to traverse the atoms */
    Bond b;                     /* to traverse the bonds */

    for (i = this.bondcnt; --i >= 0; ) {
      b = this.bonds[i];        /* traverse the bonds and */
      k = b.isInRing() ? 2 : 0; /* get the index depending on */
      b.type     &= masks[k+1]; /* whether the bond is in a ring, */
      b.src.type &= masks[k];   /* then mask the bond type and */
      b.dst.type &= masks[k];   /* the types of the incident atoms */
    }
    for (i = this.atomcnt; --i >= 0; ) {
      a = this.atoms[i];        /* traverse the atoms */
      if (a.bondcnt <= 0) a.type &= masks[0];
    }                           /* mask atoms without incident bonds */
  }  /* maskTypes() */

  /*------------------------------------------------------------------*/

  public void encode (TypeMap map)
  {                             /* --- encode atom types */
    this.map = map;             /* store the label map */
    for (int i = this.atomcnt; --i >= 0; )
      this.atoms[i].type = this.map.encode(this.atoms[i].type);
  }  /* encode() */             /* encode the atoms */

  /*------------------------------------------------------------------*/

  public void decode ()
  {                             /* --- decode atom types */
    if (this.map == null)       /* if the molecule is not encoded, */
      return;                   /* abort the function */
    for (int i = this.atomcnt; --i >= 0; )
      this.atoms[i].type = this.map.decode(this.atoms[i].type);
    this.map = null;            /* decode the atoms and */
  }  /* decode() */             /* delete the label map */

  /*------------------------------------------------------------------*/

  protected boolean trim (boolean remove)
  {                             /* --- trim molecules */
    int  i, k, n = 0;           /* loop variables */
    Atom a;                     /* to traverse the atoms */
    Bond b;                     /* to traverse the bonds */

    if (this.map != null) {     /* if there may be excluded atoms */
      for (i = this.atomcnt; --i >= 0; ) {
        a = this.atoms[i];      /* traverse the atoms */
        if (!this.map.isExcluded(a.type))
          continue;             /* skip not excluded atoms */
        a.mark = -2; n++;       /* mark excluded atoms for deletion */
        for (k = a.bondcnt; --k >= 0; )
          a.bonds[k].mark = -2; /* also mark all incident bonds */
      }                         /* (actual deletion is done later) */
    }
    if (n <= 0)  return false;  /* if no excluded atoms found, abort */
    if (!remove) return true;   /* if not to remove atoms, abort */
    for (i = this.atomcnt; --i >= 0; ) {
      a = this.atoms[i];        /* traverse the atoms again, */
      if (a.mark < -1) continue;/* but skip marked (deleted) atoms */
      for (k = n = 0; k < a.bondcnt; k++) {
        b = a.bonds[k];         /* traverse the bonds of the atom */
        if (b.mark >= -1) a.bonds[n++] = b;
      }                         /* delete all edges that */
      a.bondcnt = n;            /* lead to deleted atoms and */
    }                           /* set the new bond counter */
    for (i = n = 0; i < this.atomcnt; i++) {
      a = this.atoms[i];        /* traverse the atoms again */
      if (a.mark >= -1) this.atoms[n++] = a;
    }                           /* remove the deleted atoms */
    this.atomcnt = n;           /* set the new atom counter */
    for (i = n = 0; i < this.bondcnt; i++) {
      b = this.bonds[i];        /* traverse the bonds */
      if (b.mark >= -1) this.bonds[n++] = b;
    }                           /* remove the deleted bonds */
    this.bondcnt = n;           /* set the new bond counter */
    return true;                /* return 'atoms removed' */
  }  /* trim() */  

  /*------------------------------------------------------------------*/

  protected boolean sort () { return this.sort(false); }

  protected boolean sort (boolean seed)
  {                             /* --- sort bonds and atoms */
    int  i, k, na, nb;          /* loop variables, counter */
    Atom a, d, x;               /* to traverse the atoms */
    Bond b, v[];                /* to traverse the bonds, buffer */

    for (i = this.atomcnt; --i >= 0; )
      this.atoms[i].sortBonds();/* sort the bonds of all atoms */
    if (!seed) return true;     /* if not to sort atoms, abort */
    this.markId();              /* mark atoms and bonds */
    for (x = this.atoms[i = 0]; ++i < this.atomcnt; ) {
      a = this.atoms[i];        /* traverse the atoms and */
      if (a.type < x.type) x = a;
    }                           /* find the "minimal" atom */
    a = this.atoms[0];          /* and swap it to front */
    this.atoms[a.mark = x.mark] = a;
    this.atoms[0] = x; x.mark = -1;
    v = new Bond[this.bondcnt]; /* create a bond vector buffer */
    for (i = na = nb = 0; i < this.atomcnt; i++) {
      if (i > na) {             /* check for a connected molecule */
        this.unmark(); return false; }
      a = this.atoms[i];        /* traverse the atoms breadth-first */
      for (k = 0; k < a.bondcnt; k++) {
        b = a.bonds[k];         /* traverse the bonds of each atom */
        if (b.mark < 0) continue;
        b.mark  = -1;           /* mark the bond as processed */
        v[nb++] =  b;           /* and store it (sort bonds) */
        d = (b.src != a) ? b.src : b.dst;
        if (d.mark < 0) continue;
        x = this.atoms[++na];   /* swap dest. atom to next position */
        this.atoms[x.mark = d.mark] = x;
        this.atoms[na] = d; d.mark = -1;
      }                         /* bring atoms into depth-first */
    }                           /* search order for embedding */
    this.bonds = v;             /* set the new bond vector */
    return true;                /* return "connected molecule" */
  }  /* sort() */

  /*--------------------------------------------------------------------
  The bonds of the atoms are sorted, because their order is exploited
  to restrict the search for common substructures (search tree pruning).
  If the parameter "seed" is true, the atoms are also brought into a
  depth first search order for easier seed/substructure embedding.
  --------------------------------------------------------------------*/

  protected Embedding embed (int type)
  {                             /* --- find embeddings of an atom */
    int       i;                /* loop variable */
    Embedding emb, list = null; /* created (list of) embedding(s) */

    for (i = this.atomcnt; --i >= 0; ) {
      if (this.atoms[i].type != type)
        continue;               /* traverse atoms of the given type */
      emb = new Embedding(this, i);
      emb.succ = list;          /* create a new embedding */
      list     = emb;           /* and add it at the head */
    }                           /* of the embedding list */
    return list;                /* return the created embeddings */
  }  /* embed() */

  /*------------------------------------------------------------------*/

  private Embedding embed (Molecule seed, int bid, int aid,
                           Atom[] avec, Embedding list)
  {                             /* --- embed a seed recursively */
    int       i, m;             /* loop variable, buffer */
    Atom      rd, a, d;         /* to access/traverse the atoms */
    Bond      rb, b;            /* to access/traverse the bonds */
    Embedding emb;              /* created embedding */

    if (bid >= seed.bondcnt) {  /* if all bonds have been matched */
      if (list == Molecule.check)  /* if only to check containment, */
        return Molecule.found;  /* return a special dummy embedding */
      emb = new Embedding(this, avec, bid);
      emb.succ = list;          /* create a new embedding */
      return list = emb;        /* and add it at the head */
    }                           /* of the embedding list */
    rb = seed.bonds[bid];       /* get the next bond to embed */
    rd = rb.dst; a = rb.src;    /* and its incident atoms */
    if ((a.mark < 0) || ((rd.mark >= 0) && (a.mark > rd.mark))) {
      rd = a; a = rb.dst; }     /* identify source and dest. atom */
    m = rd.mark;                /* note the dest. atom marker */
    a = avec[a.mark];           /* get corresponding source atom */
    for (i = 0; i < a.bondcnt; i++) {
      b = a.bonds[i];           /* traverse the unmarked bonds */
      if (b.mark >= 0)       continue;     /* compare */
      if (b.type <  rb.type) continue;     /* bond type */
      if (b.type >  rb.type) return list;  /* and flags */
      d = (b.src != a) ? b.src : b.dst;
      if (d.type <  rd.type) continue;     /* compare dest. */
      if (d.type >  rd.type) return list;  /* atom type */
      if (d.mark != rd.mark) continue;     /* and index */
      b.mark = bid;             /* mark the bond (and the atom) */
      if (m < 0) { avec[d.mark = rd.mark = aid++] = d; }
      list = this.embed(seed, bid+1, aid, avec, list);
      if (m < 0) { d.mark = rd.mark = -1;  aid--; }
      b.mark = -1;              /* find embeddings recursively */
      if (list == found) break; /* if only to check containment, */
    }                           /* check the recursion result */
    return list;                /* return the list of embeddings */
  }  /* embed() */

  /*------------------------------------------------------------------*/

  private Embedding embed (Molecule seed, Embedding list)
  {                             /* --- find embeddings of a seed */
    int  i, t;                  /* loop variable, type buffer */
    Atom a;                     /* to traverse the atoms */
    Atom avec[];                /* atoms of the embedding */

    a = seed.atoms[0];          /* get the first atom of the seed */
    t = a.type;                 /* and note its type for comparisons */
    if (seed.atomcnt <= 1) {    /* if there is only one seed atom */
      if (list == null)         /* if to embed only a single atom, */
        return this.embed(t);   /* call specialized function */
      for (i = this.atomcnt; --i >= 0; )
        if (this.atoms[i].type == t)
          return found;         /* search for the single seed atom */
      return null;              /* among the molecule's atoms */
    }                           /* and report the result */
    if ((seed.atomcnt > this.atomcnt)
    ||  (seed.bondcnt > this.bondcnt))
      return null;              /* the seed must not be larger */
    a.mark = 0;                 /* mark the first seed atom */
    avec   = new Atom[seed.atomcnt];
    for (i = this.atomcnt; --i >= 0; ) {