molecule.java

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

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

  protected Molecule (Molecule mol)
  {                             /* --- duplicate a molecule */
    int  i, k;                  /* loop variables */
    Atom a, src, dst;           /* to traverse the atoms */
    Bond b;                     /* to traverse the bonds */

    this.id      = mol.id;      /* copy the molecule identifier */
    this.group   = mol.group;   /* and the group index */
    this.map     = mol.map;     /* copy the label map */
    this.atomcnt = mol.atomcnt; /* and the molecule's atoms */
    this.atoms   = new Atom[mol.atomcnt];
    for (i = mol.atomcnt; --i >= 0; ) {
      a = mol.atoms[i];         /* traverse the atoms of the molecule */
      this.atoms[i] = new Atom(a.type, a.bondcnt);
      this.atoms[i].mark = a.mark;
    }                           /* create new atoms of the same types */
    this.bondcnt = mol.bondcnt; /* copy the bonds */
    this.bonds   = new Bond[mol.bondcnt];
    for (i = mol.bondcnt; --i >= 0; ) {
      b = mol.bonds[i];         /* traverse the bonds of the molecule */
      for (k = mol.atomcnt; --k >= 0; )
        if (mol.atoms[k] == b.src) break;
      src = this.atoms[k];      /* find the corresponding soure atom */
      for (k = mol.atomcnt; --k >= 0; )
        if (mol.atoms[k] == b.dst) break;
      dst = this.atoms[k];      /* find the corresponding dest. atom */
      this.bonds[i] = new Bond(src, dst, b.type);
      this.bonds[i].mark  = b.mark;
      this.bonds[i].flags = b.flags;
    }                           /* create new bonds of the same type */
    this.opt();                 /* optimize memory usage */
  }  /* Molecule() */

  /*--------------------------------------------------------------------
  The above function duplicates a molecule without changing anything
  in the original molecule. This leads to a somewhat awkward way of
  copying the bonds - using the atom markers would be more elegant.
  However, using the atom markers would render the function useless
  for debugging purposes, for which it is intended.
  --------------------------------------------------------------------*/

  private int bridges (Atom a, Bond in, int lvl)
  {                             /* --- recursively find bridges */
    int  i;                     /* loop variable */
    int  low, l;                /* lowest level of reachable atom */
    Bond out;                   /* to traverse the outgoing bonds */
    Atom d;                     /* destination of outgoing bond */

    low = a.mark = lvl++;       /* mark atom with current level */
    for (i = a.bondcnt; --i >= 0; ) {
      out = a.bonds[i];         /* traverse all outgoing bonds */
      if (out == in) continue;  /* (skip incoming bond) */
      d = (out.src != a) ? out.src : out.dst;
      l = d.mark;               /* get the destination atom marker */
      if (l < 0) {              /* if the dest. atom is not marked, */
        l = this.bridges(d, out, lvl);  /* find bridges recursively */
        out.markBridge(l >= lvl);       /* and mark bond as bridge  */
      }                         /* if no lower level can be reached */
      if (l < low) low = l;     /* update the lowest reachable level */
    }                           /* if necessary */
    return low;                 /* return lowest reachable level */
  }  /* bridges() */

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

  public int markBridges ()
  {                             /* --- mark bonds that are bridges */
    int i, n = 0;               /* loop variable, number of bridges */

    for (i = this.atomcnt; --i >= 0; )
      this.atoms[i].mark = -1;  /* unmark all atoms */
    for (i = this.atomcnt; --i >= 0; ) {
      if (this.atoms[i].mark < 0)
        this.bridges(this.atoms[i], null, 0);
    }                           /* find all bridges in all components */
    for (i = this.bondcnt; --i >= 0; )
      if (this.bonds[i].isBridge()) n++;
    return n;                   /* count the marked bridges */
  }  /* markBridges() */        /* and return their number */

  /*--------------------------------------------------------------------
  The bridge finding algorithm employed above is adapted from:
    R.E. Tarjan.
    Depth First Search and Linear Graph Algorithms.
    SIAM Journal of Computing 1(2):146--160.
    Society for Industrial and Applied Mathematics,
    Philadelphia, PA, USA 1972
  --------------------------------------------------------------------*/

  private void prune (Atom a)
  {                             /* --- prune a possible branch */
    int  i;                     /* loop variable */
    Bond b;                     /* to traverse the bonds */

    while (a.mark == 1) {       /* while at the end of a branch */
      a.mark = 0;               /* unmark (remove) the atom */
      for (i = a.bondcnt; --i >= 0; )
        if (a.bonds[i].mark > 0) break;
      b = a.bonds[i];           /* find the corresponding bond and */
      b.mark = -1;              /* unmark (remove) it from the mol. */
      a = (b.src != a) ? b.src : b.dst;
      a.mark--;                 /* get the atom on the other side */
    }                           /* and decrement the bond counter */
  }  /* prune() */

  /*--------------------------------------------------------------------
  The above function removes a branch from a molecule by unmarking the
  atoms and bonds contained in it. It starts from an atom a that has to
  be at the end of a branch (although it does no harm if it is not).
  --------------------------------------------------------------------*/

  private long rings (Atom atom, Bond bond, Atom term,
                      long rings, int min, int max)
  {                             /* --- recursive ring search */
    int  i;                     /* loop variable */
    long f;                     /* new ring flags */
    Bond b;                     /* to traverse the bonds */
    Atom d;                     /* to traverse the atoms */

    for (i = atom.bondcnt; --i >= 0; ) /* collect used (barred) */
      rings |= atom.bonds[i].flags;    /* ring flags on the path */
    --min; --max;               /* decrement the ring bond counters */
    if (atom == term) {         /* if a ring has been found */
      if (min > 0) return 0;    /* if min. size not reached, abort */
      for (f = 1; (f & Bond.RINGS) != 0; f <<= 1)
	if ((rings & f) == 0) break;
      if ((f & Bond.RINGS) == 0) return 0;
      // if ((f & Bond.RINGS) == 0) {
      //   System.out.println("warning: too many rings\n"); return 0; }
      bond.flags |= f;          /* find an unused ring flag, set it */
      return f;                 /* in the ring flags of the bond, */
    }                           /* and return the ring flag */
    if ((atom.mark < 0)         /* check for an already visited atom */
    ||  (max <= 0)) return 0;   /* and for the maximum search depth */
    atom.mark = -atom.mark;     /* mark the atom as visited */
    for (f = 0, i = atom.bondcnt; --i >= 0; ) {
      b = atom.bonds[i];        /* traverse the bonds of the atom */
      if ((b.mark <= 0) || (b == bond))
        continue;               /* skip unmarked bonds and back bond */
      d = (b.src != atom) ? b.src : b.dst;
      f |= this.rings(d, b, term, rings|f, min, max);
    }                           /* recursively mark rings */
    atom.mark = -atom.mark;     /* remove visited marker */
    bond.flags |= f;            /* add the rings flags to the bond */
    return f;                   /* return the new ring flags */
  }  /* rings() */

  /*--------------------------------------------------------------------
  The above function does a depth first search to find rings in a
  molecule and marks each ring it finds with a unique bit that does not
  conflict with a bit already used for another ring (coloring problem).
  --------------------------------------------------------------------*/

  public int markRings (int min, int max)
  { return this.markRings(min, max, Bond.RINGBOND); }

  public int markRings (int min, int max, int typeflag)
  {                             /* --- mark rings of size [min,max] */
    int  i, cnt = 0;            /* loop variable, counter */
    long r, f;                  /* buffer for ring bond flags */
    Bond b;                     /* to traverse the bonds */

    for (i = this.bondcnt; --i >= 0; ) {
      b = this.bonds[i];        /* note number of bonds for atoms, */
      b.mark   = 1;             /* mark the bonds (for pruning) and */
      b.flags &= ~Bond.RINGS;   /* clear the bond's ring flags */
      b.type  &= ~typeflag;     /* (individual ring flags as */
    }                           /* well as the general type flag) */
    if (max <= 0) return 0;     /* if nothing to mark, abort */
    for (i = this.atomcnt; --i >= 0; )
      this.atoms[i].mark = this.atoms[i].bondcnt;
    for (i = this.atomcnt; --i >= 0; )
      this.prune(this.atoms[i]);/* remove all initial branches */
    for (i = this.bondcnt; --i >= 0; ) {
      b = this.bonds[i];        /* traverse the remaining bonds and */
      if (b.mark <= 0) continue;/* find rings by depth first search */
      r = this.rings(b.dst, b, b.src, 0, min, max);
      for (f = 1; (r != 0) && (f != 0); f <<= 1)
        if ((r & f) != 0) { r &= ~f; cnt++; }
      b.mark = 0;               /* count the marked rings, */
      b.src.mark--;             /* unmark the processed bond, */
      b.dst.mark--;             /* and reduce the bond counters */
      this.prune(b.src);        /* prune possible new branches */
      this.prune(b.dst);        /* that have been created */
    }                           /* by removing the bond */
    if (typeflag == 0) return cnt;
    for (i = this.bondcnt; --i >= 0; ) {
      b = this.bonds[i];        /* traverse all bonds */
      if ((b.flags & Bond.RINGS) != 0) b.type |=  typeflag;
      else                             b.type &= ~typeflag;
    }                           /* set a ring bond flag in the type */
    return cnt;                 /* return the number of rings */
  }  /* markRings() */

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

  public int markPseudo (int max)
  {                             /* --- mark pseudo-rings up to max */
    int  i, k, cnt = 0;         /* loop variable, buffer, counter */
    long r, f;                  /* buffer for ring bond flags */
    Bond b;                     /* to traverse the bonds */
    Atom a;                     /* to traverse the atoms */

    if (max <= 0) return 0;     /* if nothing to mark, abort */
    for (i = this.bondcnt; --i >= 0; ) {
      b = this.bonds[i];        /* traverse the bonds */
      b.mark = ((b.type & Bond.RINGBOND) != 0) ? 1 : 0;
    }                           /* mark all ring bonds */
    for (i = this.atomcnt; --i >= 0; ) {
      a = this.atoms[i]; a.mark = 0;
      for (k = a.bondcnt; --k >= 0; )
        if (a.bonds[k].mark > 0) a.mark++;
    }                           /* count the marked incident bonds */
    for (i = this.bondcnt; --i >= 0; ) {
      b = this.bonds[i];        /* traverse the remaining bonds */
      if (b.mark <= 0) continue;/* and find and mark pseudo-rings */
      r = this.rings(b.dst, b, b.src, b.flags & Bond.RINGS, 0, max);
      for (f = 1; (r != 0) && (f != 0); f <<= 1)
        if ((r & f) != 0) { r &= ~f; cnt++; }
      b.mark = 0;               /* count the marked rings, */
      b.src.mark--;             /* unmark the processed bond, */
      b.dst.mark--;             /* and reduce the bond counters */
      this.prune(b.src);        /* prune possible new branches */
      this.prune(b.dst);        /* that have been created */
    }                           /* by removing the bond */
    for (i = this.bondcnt; --i >= 0; ) {
      b = this.bonds[i];        /* traverse all bonds */
      if ((b.flags & Bond.RINGS) != 0) b.type |=  Bond.RINGBOND;
      else                             b.type &= ~Bond.RINGBOND;
    }                           /* set a ring bond flag in the type */
    return cnt;                 /* return the number of rings */
  }  /* markPseudo() */

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

  protected int aromatize ()
  {                             /* --- convert Kekul\'e to aromatic */
    int    i, k, n, cnt = 0;    /* loop variables, counter */
    long   all, cur;            /* buffer for ring flags */
    Atom   a;                   /* to traverse the atoms */
    Bond   b, bb = null;        /* to traverse the bonds */
    Bond[] ring = new Bond[6];  /* bonds of an aromatic ring */

    this.markRings(6, 6);       /* mark rings of size 6 */
    for (i = this.bondcnt; --i >= 0; )
      this.bonds[i].mark = 1;   /* mark all bonds as unprocessed */
    for (i = this.bondcnt; --i >= 0; ) {
      b = this.bonds[i];        /* traverse the bonds that */
      if (b.mark <= 0) continue;/* have not been processed yet */
      k = b.type & Bond.TYPEMASK;
      if ((k != Bond.SINGLE) && (k != Bond.DOUBLE))
        continue;               /* process only single/double bonds */
      all     = b.getRings();   /* get the ring flags and */
      ring[0] = b;              /* note the first bond of the ring */
      for (cur = 1; all != 0; cur <<= 1) {
        if ((all & cur) == 0)   /* traverse the rings */
          continue;             /* the bond is part of */
        all &= ~cur;            /* remove processed ring flag */
        b = ring[n = 0]; a = b.dst;
        do {                    /* collect the ring bonds */
          for (k = a.bondcnt; --k >= 0; ) {