---

BaseMolecule.C

Go to the documentation of this file.

/***************************************************************************
 *cr                                                                       
 *cr            (C) Copyright 1995-2008 The Board of Trustees of the           
 *cr                        University of Illinois                       
 *cr                         All Rights Reserved                        
 *cr                                                                   
 ***************************************************************************/

/***************************************************************************
 * RCS INFORMATION:
 *
 *      $RCSfile: BaseMolecule.C,v $
 *      $Author: johns $        $Locker:  $             $State: Exp $
 *      $Revision: 1.209 $      $Date: 2008/04/28 22:08:46 $
 *
 ***************************************************************************
 * DESCRIPTION:
 *
 * Base class for all molecules, without display-specific information.  This
 * portion of a molecule contains the structural data, and all routines to
 * find the structure (backbone, residues, etc).  It does NOT contain the
 * animation list; that is maintained by Molecule (which is derived from
 * this class).
 *
 ***************************************************************************/

#include <ctype.h>
#include <stdlib.h>
#include <math.h>
#include "Inform.h"
#include "utilities.h"
#include "intstack.h"

#include "BaseMolecule.h"
#include "VolumetricData.h"

#define MAXBONDERRORS 25




BaseMolecule::BaseMolecule(int myID) : residueList(10) , fragList(1),
#if defined(VMDWITHCARBS)
     pfragList(1), nfragList(1), smallringList(), smallringLinkages(), currentMaxRingSize(-1), currentMaxPathLength(-1), ID(myID) {
#else
     pfragList(1), nfragList(1), ID(myID) {
#endif

  // initialize all variables
  nAtoms = 0;
  cur_atom = 0;
  atomList = NULL;
  moleculename = NULL;
  lastbonderratomid=-1;
  bonderrorcount=0;
}



BaseMolecule::~BaseMolecule(void) {
  int i;

  // delete structural data
  delete [] atomList;
  for (i=0; i<residueList.num(); i++) {
    delete residueList[i];
  }
  for (i=0; i<nfragList.num(); i++) {
    delete nfragList[i];
  }
#if defined(VMDFASTRIBBONS)
  for (i=0; i<nfragCPList.num(); i++) {
    delete [] nfragCPList[i];
  }
#endif
  for (i=0; i<pfragList.num(); i++) {
    delete pfragList[i];
  }
#if defined(VMDFASTRIBBONS)
  for (i=0; i<pfragCPList.num(); i++) {
    delete [] pfragCPList[i];
  }
#endif
  for (i=0; i<fragList.num(); i++) {
    delete fragList[i];
  }
  for (i=0; i<volumeList.num(); i++) {
    delete volumeList[i];
  }
  for (i=0; i<extra.num(); i++) {
    delete [] extra.data(i);
  }
}



// initialize the atom list ... should be called before adding any atoms
int BaseMolecule::init_atoms(int n) {

  if(n <= 0) {
    msgErr << "BaseMolecule: init_atoms called with invalid number of atoms: "
           << n << sendmsg;
    return FALSE;
  }
  if (cur_atom != 0 && nAtoms != n) {
    msgErr << "BaseMolecule: attempt to init atoms while structure building in progress!" << sendmsg;
    return FALSE;
  }

  if (!atomList) {
    // first call to init_atoms
    nAtoms = n; // only place where nAtoms is set!
    atomList = new MolAtom[nAtoms];
    memset(atomList, 0, nAtoms*sizeof(MolAtom));

    // initialize NULL extra data field, which is returned when
    // querying a non-existent field with extra.data("fielddoesntexist")
    extra.add_name("NULL", NULL);

    // initialize "default" extra data fields.
    extra.add_name("beta", new float[nAtoms]);
    extra.add_name("occupancy", new float[nAtoms]);
    extra.add_name("charge", new float[nAtoms]);
    extra.add_name("mass", new float[nAtoms]);
    extra.add_name("radius", new float[nAtoms]);
    for (int i=0; i<extra.num(); i++) {
      void *data = extra.data(i);
      if (data != NULL) 
        memset(data, 0, nAtoms*sizeof(float));
    }
    return TRUE;
  }
  if (n != nAtoms) {
    msgErr << "The number of atoms in this molecule has already been assigned."
           << sendmsg;
    return FALSE;
  }
  return TRUE;
}

// add a new atom; return it's index, or (-1) if error.
int BaseMolecule::add_atom(char *name, char *atomtype, int atomicnumber, 
    char *resname, int resid, const char *chain, const char *segname, 
    char *insertion, const char *altloc) {
  int nameindex, typeindex;
  int resnameindex, segnameindex, altlocindex, chainindex;

  if (!atomList || cur_atom >= nAtoms) {
    msgErr << "BaseMolecule: Cannot add new atom; currently " << nAtoms;
    msgErr << " atoms in structure." << sendmsg;
    return (-1);
  }

  // create atom  
  MolAtom *newatom = atom(cur_atom);
  newatom->init(cur_atom, resid, insertion);

  // add names to namelist, and put indices in MolAtom object
  nameindex = atomNames.add_name(name,atomNames.num());
  typeindex = atomTypes.add_name(atomtype, atomTypes.num());
  resnameindex = resNames.add_name(resname, resNames.num());
  segnameindex = segNames.add_name(segname, segNames.num());
  altlocindex = altlocNames.add_name(altloc, altlocNames.num());

  // use default of 'X' for chain if not given
  if(!chain || ! (*chain) || *chain == ' ')
    chainindex = chainNames.add_name("X", chainNames.num());
  else
    chainindex = chainNames.add_name(chain, chainNames.num());

  // set atom member variables
  newatom->nameindex = nameindex;
  newatom->typeindex = typeindex;
  newatom->atomicnumber = atomicnumber;
  newatom->resnameindex = resnameindex;
  newatom->segnameindex = segnameindex;
  newatom->altlocindex = altlocindex;
  newatom->chainindex = chainindex;

  // check for integer overflow/wraparound condition, which can occur
  // if an evil plugin defines 100,000 unique atom names, for example
  if (newatom->nameindex != nameindex ||
      newatom->typeindex != typeindex ||
      newatom->atomicnumber != atomicnumber ||
      newatom->resnameindex != resnameindex ||
      newatom->segnameindex != segnameindex ||
      newatom->altlocindex != altlocindex ||
      newatom->chainindex != chainindex) {
    msgErr << "BaseMolecule: Cannot add atom; namelist index value too large." << sendmsg;
    msgErr << "Recompile VMD with larger index types." << sendmsg;
    msgErr << "Atom namelist index values at time of overflow:" << sendmsg;
    msgErr << "  nameindex: " << nameindex << sendmsg;;
    msgErr << "  typeindex: " << typeindex << sendmsg;;
    msgErr << "  resnameindex: " << resnameindex << sendmsg;;
    msgErr << "  segnameindex: " << segnameindex << sendmsg;;
    msgErr << "  altlocindex: " << altlocindex << sendmsg;;
    msgErr << "  chainindex: " << chainindex << sendmsg;
    return -1;
  }

  return cur_atom++;
}


// add a new bond; return 0 on success, or (-1) if error.
int BaseMolecule::add_bond(int a, int b, float bondorder, int backbonetype) {
  if (!nAtoms || a >= nAtoms || b >= nAtoms) {
    msgErr << "BaseMolecule: Atoms must be added before bonds." << sendmsg;
    return (-1);
  } 

  if (a == b) {
    msgErr << "BaseMolecule: Cannot bond atom " <<a<< " to itself." << sendmsg;
    return (-1);
  }

  // put the bond in the atom list
  if (atom(a)->add_bond(b, backbonetype)) {
    if (bonderrorcount < MAXBONDERRORS) {
      if (lastbonderratomid != a) {
        msgErr << "MolAtom " << a << ": Exceeded maximum number of bonds ("
               << atom(a)->bonds << ")." << sendmsg;
        lastbonderratomid=a;
        bonderrorcount++;
      }
    } else if (bonderrorcount == MAXBONDERRORS) {
      msgErr << "BaseMolecule: Excessive bonding errors encountered, perhaps atom coordinates are in the wrong units?" << sendmsg;
      msgErr << "BaseMolecule: Silencing bonding error messages." << sendmsg;
      bonderrorcount++;
    }
    return (-1);
  }

  if (atom(b)->add_bond(a, backbonetype)) {
    if (bonderrorcount < MAXBONDERRORS) {
      if (lastbonderratomid != b) {
        msgErr << "MolAtom " << b << ": Exceeded maximum number of bonds ("
               << atom(b)->bonds << ")." << sendmsg;
        lastbonderratomid=b;
        bonderrorcount++;
      }
    } else if (bonderrorcount == MAXBONDERRORS) {
      msgErr << "BaseMolecule: Excessive bonding errors encountered, perhaps atom coordinates are in the wrong units?" << sendmsg;
      msgErr << "BaseMolecule: Silencing bonding error messages." << sendmsg;
      bonderrorcount++;
    }
    return (-1);
  }

  // store bond orders
  setbondorder(a, atom(a)->bonds-1, bondorder);
  setbondorder(b, atom(b)->bonds-1, bondorder);

  return 0;
}

// Add a bond to a structure but check to make sure there isn't a 
// duplicate, as may be the case when merging bondlists from a file 
// and from a distance-based bond search
int BaseMolecule::add_bond_dupcheck(int a, int b, float bondorder) {
  int i;

  if (!nAtoms || a >= nAtoms || b >= nAtoms) {
    msgErr << "BaseMolecule: Atoms must be added before bonds." << sendmsg;
    return (-1);
  }

  MolAtom *atm = atom(a);
  int nbonds = atm->bonds;
  const int *bonds = &atm->bondTo[0];
  for (i=0; i<nbonds; i++) {
    if (bonds[i] == b) {
      return 0; // skip bond that already exists
    }
  }
  add_bond(a, b, bondorder); // add it if it doesn't already exist

  return 0;
}


void BaseMolecule::setbondorder(int atom, int bond, float order) {
  float *bondOrders = extra.data("bondorders");

  // if not already there, add it
  if (bondOrders == NULL) {
    if (order != 1) {
      int i;
      extra.add_name("bondorders", new float[nAtoms*MAXATOMBONDS]);
      bondOrders = extra.data("bondorders");
      for (i=0; i<nAtoms*MAXATOMBONDS; i++)
        bondOrders[i] = 1.0f;    

      bondOrders[atom * MAXATOMBONDS + bond] = order;
    } 
    return;
  }

  bondOrders[atom * MAXATOMBONDS + bond] = order;
}

float BaseMolecule::getbondorder(int atom, int bond) {
  float *bondOrders = extra.data("bondorders");

  // if not already there, add it
  if (bondOrders == NULL) { 
    return 1;
  }
   
  return bondOrders[atom * MAXATOMBONDS + bond];
}


// return the Nth residue
Residue *BaseMolecule::residue(int n) {
  return residueList[n];
}


// return the Nth fragment
Fragment *BaseMolecule::fragment(int n) {
  return fragList[n];
}


// given an atom index, return the residue object for the residue it
// is in.  If it is not in a residue, return NULL.
Residue *BaseMolecule::atom_residue(int n) {
  MolAtom *atm = atom(n);
  if(atm->uniq_resid < 0)
    return NULL;
  else
    return residue(atm->uniq_resid);
}


// given an atom index, return the fragment object for the fragment it
// is in.  If it is not in a fragment, return NULL.
Fragment *BaseMolecule::atom_fragment(int n) {
  MolAtom *atm = atom(n);
  int frag = residue(atm->uniq_resid)->fragment;
  if(frag < 0)
    return NULL;
  else
    return fragment(frag);
}

// return a 'default' value for a given atom name
float BaseMolecule::default_radius(char *nm) {
  float val = 1.5;
  // some names start with a number
  while (*nm && isdigit(*nm))
    nm++;
  if(nm) {
    switch(toupper(nm[0])) {
      // These are similar to the values used by X-PLOR with sigma=0.8
      // see page 50 of the X-PLOR 3.1 manual
      case 'H' : val = 1.00f; break;
      case 'C' : val = 1.50f; break;
      case 'N' : val = 1.40f; break;
      case 'O' : val = 1.30f; break;
      case 'F' : val = 1.20f; break;
      case 'S' : val = 1.90f; break;
    }
  }

  return val;
}


// return a 'default' value for a given atom name
float BaseMolecule::default_mass(char *nm) {
  float val = 12.0;

  // some names start with a number
  while (*nm && isdigit(*nm)) nm++;
  if(nm) {
    switch(toupper(nm[0])) {
      case 'H' : val = 1.00800f; break;
      case 'C' : val = 12.01100f; break;
      case 'N' : val = 14.00700f; break;
      case 'O' : val = 15.99900f; break;
      case 'F' : val = 55.84700f; break;
      case 'P' : val = 30.97376f; break;
      case 'S' : val = 32.06000f; break;
    }
    if      (toupper(nm[0] == 'C') && toupper(nm[1] == 'L')) val = 35.453f;
    else if (toupper(nm[0] == 'N') && toupper(nm[1] == 'A')) val = 22.989770f;
    else if (toupper(nm[0] == 'M') && toupper(nm[1] == 'G')) val = 24.3050f;
  }

  return val;
}


float BaseMolecule::default_charge(char *) {
  // return 0 for everything; later, when we put in a more reliable
  // system for determining the charge that's user-configurable,
  // we can start assigning more realistic charges.
  return 0.0f;
}


// count the number of unique bonds in the structure
int BaseMolecule::count_bonds(void) {
  int i, j;
  int count=0;

  for (i=0; i<nAtoms; i++) {
    int nbonds = atomList[i].bonds;
    const int *bonds = &atomList[i].bondTo[0];

    for (j=0; j<nbonds; j++) {
      if (bonds[j] > i)
        count++;
    }
  }

  return count;
}

void BaseMolecule::clear_bonds(void) {
  int i;
  for (i=0; i<nAtoms; i++)
    atomList[i].bonds = 0;
}


// analyze the molecule for more than just the atom/bond information
// This is here since it is called _after_ the molecule is added to
// the MoleculeList.  Thus, there is a Tcl callback to allow the
// user to update the bond information (or other fields?) before
// the actual search.
void BaseMolecule::analyze(void) {
  need_find_bonds = 0; // at this point it's too late

  // I have to let 0 atoms in because I want to be able to read things
  // like electron density maps, which have no atoms.
  // It is kinda wierd, then to make BaseMolecule be at the top of the
  // heirarchy.  Oh well.
  if(nAtoms < 1)
    return;

  // call routines to find different characteristics of the molecule
  msgInfo << "Analyzing structure ..." << sendmsg;
  msgInfo << "   Atoms: " << nAtoms << sendmsg;

  // count unique bonds
  msgInfo << "   Bonds: " << count_bonds() << sendmsg;

  // restore residue and fragment lists to pristine state
  residueList.clear();
  fragList.clear();
  pfragList.clear();   
  pfragCyclic.clear(); 
#if defined(VMDFASTRIBBONS)
  pfragCPList.clear(); 
#endif
  nfragList.clear();   
  nfragCyclic.clear(); 
#if defined(VMDFASTRIBBONS)
  nfragCPList.clear(); 
#endif

  // assign per-atom backbone types
  find_backbone();

  // find all the atoms in a resid connected to DNA/RNA/PROTEIN/WATER
  // also, assign a unique resid (uniq_resid) to each atom
  nResidues = find_residues();
  msgInfo << "   Residues: " << nResidues << sendmsg;

  nWaters = find_waters();
  msgInfo << "   Waters: " << nWaters << sendmsg;
  
  // determine which residues are connected to each other
  bonderrorcount=0; // reset error count before residue connectivity search
  find_connected_residues(nResidues); 

 
  nSegments = find_segments(); 
  msgInfo << "   Segments: " << nSegments << sendmsg;

  nFragments = find_fragments();
  msgInfo << "   Fragments: " << nFragments;

  nProteinFragments = pfragList.num();
  msgInfo << "   Protein: " << nProteinFragments;

  nNucleicFragments = nfragList.num();
  msgInfo << "   Nucleic: " << nNucleicFragments << sendmsg;
  
  // NOTE: The current procedure incorrectly identifies some lipid 
  // atoms as "ATOMNUCLEICBACK" (but not as "nucleic") as well
  // as some single water oxygens as "backbone". Here, we 
  // correct this by setting all atoms of non-polymeric residue types
  // to be "ATOMNORMAL" (i.e.: not backbone).
  MolAtom *a;
  int i;
  for (i=0; i<nAtoms; i++) {
    a = atom(i); 
    if ((a->residueType != RESNUCLEIC) && (a->residueType != RESPROTEIN)) 
      a->atomType = ATOMNORMAL;
  }

  // Search for hydrogens
  // XXX Must be done after the rest of the structure finding routines,
  // because those routines assume that anything that isn't NORMAL is
  // a backbone atom.
  // We use the name-based definition used in the IS_HYDROGEN macro
  for (i=0; i<nAtoms; i++) {
    MolAtom *a = atom(i);
    const char *aname = atomNames.name(a->nameindex);
    if (IS_HYDROGEN(aname))
      a->atomType = ATOMHYDROGEN;
  }

#if defined(VMDFASTRIBBONS)
  calculate_ribbon_controlpoints();
#endif
  
}


int BaseMolecule::find_backbone(void) {
  int i, j, k;

  // Search for the protein backbone
  int protypes[4];
  protypes[0] = atomNames.typecode((char *) "CA");
  protypes[1] = atomNames.typecode((char *) "C");
  protypes[2] = atomNames.typecode((char *) "O");
  protypes[3] = atomNames.typecode((char *) "N");

  // special case for terminal oxygens that miss the search for O
  // by looking for ones connected to a C
  int termtypes[4];
  termtypes[0] = atomNames.typecode((char *) "OT1"); // standard names
  termtypes[1] = atomNames.typecode((char *) "OT2");
  termtypes[2] = atomNames.typecode((char *) "O1");  // gromacs force field 
  termtypes[3] = atomNames.typecode((char *) "O2");  // atom names

  // search for the DNA/RNA backbone;  the atom names are:
  // for the phosphate:  P, O1P, O2P, OP1, OP2
  // for the rest: O3', C3', C4', C5', O5'
  // (or O3*, C3*, C4*, C5*, O5*)
  int nuctypes[15];
  nuctypes[ 0] = atomNames.typecode((char *) "P");
  nuctypes[ 1] = atomNames.typecode((char *) "O1P"); // old PDB files
  nuctypes[ 2] = atomNames.typecode((char *) "O2P"); // old PDB files
  nuctypes[ 3] = atomNames.typecode((char *) "OP1"); // new PDB files
  nuctypes[ 4] = atomNames.typecode((char *) "OP2"); // new PDB files
  nuctypes[ 5] = atomNames.typecode((char *) "C3*");
  nuctypes[ 6] = atomNames.typecode((char *) "C3'");
  nuctypes[ 7] = atomNames.typecode((char *) "O3*");
  nuctypes[ 8] = atomNames.typecode((char *) "O3'");
  nuctypes[ 9] = atomNames.typecode((char *) "C4*");
  nuctypes[10] = atomNames.typecode((char *) "C4'");
  nuctypes[11] = atomNames.typecode((char *) "C5*");
  nuctypes[12] = atomNames.typecode((char *) "C5'");
  nuctypes[13] = atomNames.typecode((char *) "O5*");
  nuctypes[14] = atomNames.typecode((char *) "O5'");

#if 0
  // non-backbone nucleic acid atom names
  nuctypes[ 3] = atomNames.typecode((char *) "C1*");
  nuctypes[ 4] = atomNames.typecode((char *) "C1'");
  nuctypes[ 5] = atomNames.typecode((char *) "C2*");
  nuctypes[ 6] = atomNames.typecode((char *) "C2'");
  nuctypes[ 7] = atomNames.typecode((char *) "O2*");
  nuctypes[ 8] = atomNames.typecode((char *) "O2'");
  nuctypes[15] = atomNames.typecode((char *) "O4*");
  nuctypes[16] = atomNames.typecode((char *) "O4'");
#endif

  // special case for terminal nucleic residues
  int nuctermtypes[2];
  nuctermtypes[0] = atomNames.typecode((char *) "H5T"); // standard names
  nuctermtypes[1] = atomNames.typecode((char *) "H3T");


  // loop over all atoms assigning atom backbone type flags
  for (i=0; i<nAtoms; i++) {
    MolAtom *a = atom(i);
 
    // initialize atom type to non-backbone
    a->atomType = ATOMNORMAL;

    // check for protein backbone atom names
    for (j=0; j < 4; j++) {
      if (a->nameindex == protypes[j]) {
        a->atomType = ATOMPROTEINBACK;
        break;
      }
    }

    // check terminal residue names as well
    for (j=0; j < 4; j++) {
      if (a->nameindex == termtypes[j]) { // check if OT1, OT2
        for (k=0; k < a->bonds; k++) {
          if (atom(a->bondTo[k])->atomType == ATOMPROTEINBACK) {
            a->atomType = ATOMPROTEINBACK;
            break;
          }
        }
      }
    }
  
    // check if in nucleic backbone, if not already set
    if(!(a->atomType)) {
      for (j=0; j < 13; j++) {
        if (a->nameindex == nuctypes[j]) {
          a->atomType = ATOMNUCLEICBACK;
          break;
        }
      }
    }

    // check if nucleic terminal atom names
    for (j=0; j < 2; j++) {
      if (a->nameindex == nuctermtypes[j]) {
        for (k=0; k < a->bonds; k++) {
          if (atom(a->bondTo[k])->atomType == ATOMNUCLEICBACK) {
            a->atomType = ATOMNUCLEICBACK;
            break;
          }
        }
      }
    }
  }

  return 0; 
}



// find water molecules based on the residue name
// from the documentation for molscript, these are possible
// waters:
// type H2O HH0 OHH HOH OH2 SOL WAT
// as well, I add TIP, TIP2, TIP3, and TIP4
// The count is the number of sets of connected RESWATERS
int BaseMolecule::find_waters(void) {
  int i, j;
  MolAtom *a;

  int watertypes[12];
  watertypes[0] = resNames.typecode((char *) "H2O");
  watertypes[1] = resNames.typecode((char *) "HH0");
  watertypes[2] = resNames.typecode((char *) "OHH");
  watertypes[3] = resNames.typecode((char *) "HOH");
  watertypes[4] = resNames.typecode((char *) "OH2");
  watertypes[5] = resNames.typecode((char *) "SOL");
  watertypes[6] = resNames.typecode((char *) "WAT");
  watertypes[7] = resNames.typecode((char *) "TIP");
  watertypes[8] = resNames.typecode((char *) "TIP2");
  watertypes[9] = resNames.typecode((char *) "TIP3");
  watertypes[10] = resNames.typecode((char *) "TIP4");
  // this conflicts with a PDB compound:
  //   http://minerva.roca.csic.es/hicup/SPC/spc_pdb.txt
  //
  // XXX we should add a check to make sure that there are only 
  //     three atoms in the residue when all is said and done. If there
  //     are more, then we should undo the assignment as a water and mark it
  //     as protein etc as appropriate.  This is tricky since at this stage
  //     we haven't done any connectivity tests etc, and we're only working
  //     with individual atoms.  Perhaps its time to re-think this logic.
  watertypes[11] = resNames.typecode((char *) "SPC");

  for (i=0; i<nAtoms; i++) {
    a = atom(i);
    if (a->residueType == RESNOTHING) {  // make sure it isn't named yet
      for (j=0; j<12; j++) {
        if (watertypes[j] == a->resnameindex) {
          a->residueType = RESWATERS;
          break;
        }
      }
    }
  }
 
  int count = find_connected_waters2();

  return count;   
}


// if this is a RESWATERS with index idx, mark it and find if
// any of its neighbors are RESWATERS
// this does a depth-first search with RECURSION.
void BaseMolecule::find_connected_waters(int i, char *tmp) {
  MolAtom *a = atom(i);
  int j;
  if (a->residueType == RESWATERS && !tmp[i]) {
    tmp[i] = TRUE;
    for (j=0; j<a->bonds; j++) {
      find_connected_waters(a->bondTo[j], tmp);
    }
  }
}


// if this is a RESWATERS with index idx, mark it and find if
// any of its neighbors are RESWATERS
int BaseMolecule::find_connected_waters2(void) {
  MolAtom *a;
  int count, i;
  IntStackHandle s;

  char *tmp = new char[nAtoms];
  memset(tmp, 0, nAtoms * sizeof(char));

  s = intstack_create(nAtoms);

  for (count=0, i=0; i<nAtoms; i++) {
    if (atom(i)->residueType == RESWATERS && !tmp[i]) {
      int nextatom;

      count++;
      intstack_push(s, i);
    
      // find and mark all connected waters 
      while (!intstack_pop(s, &nextatom)) { 
        int j;

        a = atom(nextatom);
        tmp[nextatom] = TRUE;

        for (j=a->bonds - 1; j>=0; j--) {
          int bi = a->bondTo[j];
          MolAtom *b = atom(bi);
          if (b->residueType == RESWATERS && !tmp[bi])
            intstack_push(s, bi);
        }
      }
    }
  }

  intstack_destroy(s);
  delete [] tmp;

  return count;
}


// find n backbone atoms connected together with the given residueid
// return the total count
// this assumes that the given atom (atomidx) is correct
int BaseMolecule::find_connected_backbone(IntStackHandle s, int backbone,
                         int atomidx, int residueid, int tmpid, int *flgs) {
  if (flgs[atomidx] != 0)
    return 0; // already done

  MolAtom *x = atom(atomidx);
  if (x->atomType != backbone || x->resid != residueid)
    return 0; // not a backbone atom, or resid doesn't match

  intstack_popall(s); // just in case
  intstack_push(s, atomidx);
  int nextatom;
  int count = 0;
   
  // find and mark connected backbone atoms
  while (!intstack_pop(s, &nextatom)) {
    MolAtom *a = atom(nextatom);
    flgs[nextatom] = tmpid;
    count++;

    int j;
    for (j=a->bonds - 1; j>=0; j--) {
      int bi = a->bondTo[j];
      if (flgs[bi] == 0) {
        MolAtom *b = atom(bi);

        // skip connections to atoms on different chains/segnames
        if (a->chainindex != b->chainindex || 
            a->segnameindex != b->segnameindex)
          continue;

        if (b->atomType == backbone && b->resid == residueid)
          intstack_push(s, bi);
      }
    }
  }

  return count;
}


// the find_connected_backbone left terms of flgs which need to be cleaned up
void BaseMolecule::clean_up_connection(IntStackHandle s, int i, int tmpid, int *flgs) {
  if (flgs[i] != tmpid)  // been here before
    return;

  intstack_popall(s); // just in case
  intstack_push(s, i);
  int nextatom;
 
  // find and null out non-matching atom flags
  while (!intstack_pop(s, &nextatom)) {
    flgs[nextatom] = 0;
    MolAtom *a = atom(nextatom);
    int j;
    for (j=a->bonds - 1; j>=0; j--) {
      int bi = a->bondTo[j];
      if (flgs[bi] == tmpid) {
        intstack_push(s, bi);
      }
    }
  }
}



// now that I know this resid is okay, mark it so
void BaseMolecule::find_connected_atoms_in_resid(IntStackHandle s,
    int restype, int i, int residueid, int tmpid, int *flgs)
{
  if (flgs[i] != 0 || atom(i)->resid != residueid)
    return;

  intstack_popall(s); // just in case
  intstack_push(s, i);
  int nextatom;

  // find and mark all connected residues in the same chain/segname
  while (!intstack_pop(s, &nextatom)) {
    flgs[nextatom] = tmpid;
    MolAtom *a = atom(nextatom);
    a->residueType = restype;

    int j;
    for (j=a->bonds - 1; j>=0; j--) {
      int bi = a->bondTo[j];
      MolAtom *b = atom(bi);
      if (flgs[bi] == 0 &&
          a->chainindex == b->chainindex &&
          a->segnameindex == b->segnameindex &&
          b->resid == residueid) {
        intstack_push(s, bi);
      }
    }
  }
}



// Find connected backbone atoms with the same resid
// if found, find all the atoms with the same resid
// which are connected to those backbone atoms only through
// atoms of the same resid
void BaseMolecule::find_and_mark(int n, int backbone,
  int restype, int *tmpid, int *flgs)
{
  int i;
  MolAtom *a;
  int residueid; // the real resid
  IntStackHandle s = intstack_create(nAtoms);

  for (i=0; i<nAtoms; i++) {
    a = atom(i);   // look for a new backbone atom
    if (a->atomType == backbone && flgs[i] == 0) {
      residueid = a->resid;
      if (find_connected_backbone(s, backbone, i, residueid, *tmpid, flgs) >= n) {
        // if find was successful, start all over again
        clean_up_connection(s, i, *tmpid, flgs);
        // but mark all the Atoms connected to here
        find_connected_atoms_in_resid(s, restype, i, residueid, *tmpid, flgs);
        // and one more was made
        (*tmpid)++;
      } else {
        // clean things up so I won't have problems later
        clean_up_connection(s, i, *tmpid, flgs);
      }
    }
  }

  intstack_destroy(s);
}



// assign a uniq resid (uniq_resid) to each set of connected atoms
// with the same residue id.  There could be many residues with the
// same id, but not connected (the SSN problem - SSNs are not unique
// so don't use them as the primary key)
int BaseMolecule::make_uniq_resids(int *flgs) {
  int i;
  int num_residues = 0;
  IntStackHandle s = intstack_create(nAtoms);

  for (i=0; i<nAtoms; i++) {
    if (!flgs[i]) {  // not been numbered
      // find connected atoms to i with the same resid and label
      // it with the uniq_resid
      MolAtom *a = atom(i);
      int resid = a->resid;
      char *insertion = a->insertionstr;

      intstack_push(s, i);
      int nextatom;

      // Loop over all atoms we're bonded to in the same chain/segname
      while (!intstack_pop(s, &nextatom)) {
        MolAtom *a = atom(nextatom);
        a->uniq_resid = num_residues;  // give it the new resid number
        flgs[nextatom] = TRUE;         // mark this atom done
  
        int j;
        for (j=a->bonds - 1; j>=0; j--) {
          int bi = a->bondTo[j];
          if (flgs[bi] == 0) {
            MolAtom *b = atom(bi);
            if (a->chainindex == b->chainindex && 
                a->segnameindex == b->segnameindex &&
                b->resid == resid && !strcmp(b->insertionstr, insertion))
              intstack_push(s, bi);
          }
        }
      }

      num_residues++;
    }
  }

  intstack_destroy(s);

  return num_residues;
}



int BaseMolecule::find_residues(void) {
  int *flgs = new int[nAtoms]; // flags used for connected atom searches
  memset(flgs, 0, nAtoms * sizeof(int)); // clear flags array
  
  // assign a uniq resid (uniq_resid) to each set of connected atoms
  // with the same residue id.  There could be many residues with the
  // same id, but not connected (the SSN problem - SSNs are not unique
  // so don't use them as the primary key)
  int num_residues = make_uniq_resids(flgs);
   
  int back_res_count = 1; // tmp count of number of residues on some backbone
  memset(flgs, 0, nAtoms * sizeof(int)); // clear flags array
  
  //  hunt for the proteins
  // there must be 4 PROTEINBACK atoms connected and with the same resid
  // then all connected atoms will be marked as PROTEIN atoms
  // this gets everything except the terminals
  find_and_mark(4, ATOMPROTEINBACK, RESPROTEIN, &back_res_count, flgs);
  
  // do the same for nucleic acids
  // XXX we might not want to check for the phosphate (P and 2 O's).  Here's
  // the quick way I can almost do that.  Unfortionately, that
  // messes up nfragList, since it needs a P to detect an end
  find_and_mark(4, ATOMNUCLEICBACK, RESNUCLEIC, &back_res_count, flgs);
  
  delete [] flgs;
  return num_residues;
}

int BaseMolecule::find_atom_in_residue(const char *name, int residue) {
  int nametype = atomNames.typecode(name);
  if (nametype < 0)
    return -2;

  return find_atom_in_residue(nametype, residue);
}


// find which residues are connected to which
// remember, I already have the uniq_id for each atom
void BaseMolecule::find_connected_residues(int num_residues) {
  int i, j;
  for (i=0; i<num_residues; i++) {   // init the list to NULLs
    residueList.append(NULL);
  }
  
  for (i=nAtoms-1; i>=0; i--) {      // go through all the atoms
    j = atom(i)->uniq_resid;
    if (residueList[j] == NULL) {    // then init the residue
      residueList[j] = new Residue(atom(i)->resid, atom(i)->residueType);
    }
    // Tell the residue that this atom is in it
    residueList[j]->add_atom(i);
  }

  // double check that everything was created
  for (i=0; i<num_residues; i++) {
    if (residueList[i] == NULL) { // no atom was found for this residue
      msgErr << "Mysterious residue creation in ";
      msgErr << "BaseMolecule::find_connected_residues." << sendmsg;
      residueList[i] = new Residue((int) -1, RESNOTHING);
    }
  }

  // finally, check for unusual connections between residues, e.g. between
  // protein and water.
  for (i=0; i<num_residues; i++) {
    Residue *res = residueList[i];
    int bondfromtype = res->residueType;
    int numatoms = res->atoms.num();
    for (j=0; j<numatoms; j++) {
      MolAtom *a = atom(res->atoms[j]);

      // find off-residue bonds to residues of the same chain/segname
      int k;
      for (k=0; k<a->bonds; k++) {
        MolAtom *b = atom(a->bondTo[k]);

        // skip connections to atoms on different chains/segnames
        if (a->chainindex != b->chainindex || 
            a->segnameindex != b->segnameindex)
          continue;
         
        if (b->uniq_resid != i) {
          int bondtotype = residueList[b->uniq_resid]->residueType;

          if (bondfromtype != bondtotype) {
            if (i < b->uniq_resid) { // so that we only warn once
              msgWarn << "Unusual bond between residues:  ";
              msgWarn << residueList[i]->resid;
              switch (bondfromtype) {
                case RESPROTEIN: msgWarn << " (protein)"; break;
                case RESNUCLEIC: msgWarn << " (nucleic)"; break;
                case RESWATERS:  msgWarn << " (waters)"; break;
                default:
                case RESNOTHING: msgWarn << " (none)"; break;
              }
              msgWarn << " and ";
              msgWarn << residueList[b->uniq_resid]->resid;
              switch (bondtotype) {
                case RESPROTEIN: msgWarn << " (protein)"; break;
                case RESNUCLEIC: msgWarn << " (nucleic)"; break;
                case RESWATERS:  msgWarn << " (waters)"; break;
                default:
                case RESNOTHING: msgWarn << " (none)"; break;
              }
              msgWarn << sendmsg;
            }
          }
        }
      }
    }
  }
}


// find all the residues connected to a specific residue
int BaseMolecule::find_connected_fragments(void) {
  int i;
  int count = 0;
  char *flgs = new char[residueList.num()]; // set up temp space
  memset(flgs, 0, residueList.num() * sizeof(char)); // clear flags
  IntStackHandle s = intstack_create(residueList.num());

  int atomsg = atomNames.typecode((char *) "SG"); // to find disulfide bonds

  int nextres;
  for (i=0; i<residueList.num(); i++) { // find unmarked fragment
    if (!flgs[i]) {
      fragList.append(new Fragment);
      intstack_push(s, i);

      // find and mark all connected residues with the same chain/segname
      while (!intstack_pop(s, &nextres)) {
        fragList[count]->append(nextres);
        Residue *res = residueList[nextres];
        res->fragment = count; // store residue's fragment

        int numatoms = res->atoms.num();
        int j;
        for (j=0; j<numatoms; j++) {
          MolAtom *a = atom(res->atoms[j]);

          // find all bonds to residues of the same chain/segname 
          int k;
          for (k=0; k<a->bonds; k++) {
            MolAtom *b = atom(a->bondTo[k]);
            int ri = b->uniq_resid;

            // skip connections to residues with different chains/segnames,
            // and don't follow disulfide bonds, as we want the order of
            // residue traversal to be correct so we can use it to build
            // subfragment lists later on
            if ((ri != i) &&
                (flgs[ri] == 0) &&
                (a->chainindex == b->chainindex) &&
                (a->segnameindex == b->segnameindex) &&
                ((a->nameindex != atomsg) || (b->nameindex != atomsg))) {
              flgs[ri] = TRUE;
              intstack_push(s, ri);
            }
          }
        }
      }

      count++;
    }
  }

  intstack_destroy(s);
  delete [] flgs;

  return count;
}


// find each collection of connected fragments
int BaseMolecule::find_fragments(void) {
  int count = find_connected_fragments();  // find and mark its neighbors

#if 1
  // find the protein subfragments
  find_subfragments(atomNames.typecode((char *) "N"), 
     -1,
     -1,
     atomNames.typecode((char *) "C"), 
     -1,
     -1,
     -1,
     RESPROTEIN, &pfragList);

#if 0
  // find the nucleic acid subfragments
  find_subfragments(atomNames.typecode((char *) "P"), 
     atomNames.typecode((char *) "H5T"),
     -1,
     atomNames.typecode((char *) "O3'"),
     atomNames.typecode((char *) "O3*"),
     atomNames.typecode((char *) "H3T"),
     -1,
     RESNUCLEIC, &nfragList);
#else
  // find the nucleic acid subfragments
  find_subfragments_topologically(
     RESNUCLEIC, &nfragList,
     atomNames.typecode((char *) "O3'"),
     atomNames.typecode((char *) "O3*"),
     atomNames.typecode((char *) "H3T"),
     -1);
#endif
#else
  find_subfragments_cyclic(&pfragList, RESPROTEIN);
  find_subfragments_cyclic(&nfragList, RESNUCLEIC);
#endif

  // determine whether fragments are cyclic or not
  find_cyclic_subfragments(&pfragList, &pfragCyclic);
  find_cyclic_subfragments(&nfragList, &nfragCyclic);

  return count;
}


void BaseMolecule::find_subfragments_cyclic(ResizeArray<Fragment *> *subfragList, int restype) {
  int numfrags = fragList.num();
  int i, frag;

  // test each fragment to see if it's a candidate for the subfraglist
  for (frag=0; frag<numfrags; frag++) {
    int numres = fragList[frag]->num();       // residues in this frag
    int match=1; // start true, and falsify

    // check each residue to see they are all the right restype
    for (i=0; i<numres; i++) {
      int resid = (*fragList[frag])[i];
      if (residueList[resid]->residueType != restype) {
        match=0;
        break;
      }
    }

    // if we found a matching fragment, add it to the subfraglist
    if (match) {
      Fragment *frg = new Fragment;

      // add all of the residues for this fragment to the subfraglist
      for (i=0; i<numres; i++) {
        int resid = (*fragList[frag])[i];
        frg->append(resid);
      } 

      subfragList->append(frg);
    }    
  }
}



void BaseMolecule::find_cyclic_subfragments(ResizeArray<Fragment *> *subfragList, ResizeArray<int> *subfragCyclic) {
  int i, j, frag;
  int numfrags = subfragList->num();

  // check each fragment for cycles
  for (frag=0; frag<numfrags; frag++) {
    int numres   = (*subfragList)[frag]->num();       // residues in this frag

    // skip testing fragments containing zero residues
    if (numres < 1)
      continue;

    int startres = (*(*subfragList)[frag])[0];        // first residue
    int endres   = (*(*subfragList)[frag])[numres-1]; // last residue
    int cyclic   = 0;

    // check for bonds between startres and endres
    int numatoms = residueList[endres]->atoms.num();
    int done = 0;
    for (i=0; (i < numatoms) && (!done); i++) {
      MolAtom *a = atom(residueList[endres]->atoms[i]);
      int nbonds = a->bonds;
      for (j=0; j < nbonds; j++) {
        MolAtom *b = atom(a->bondTo[j]);

        if (b->uniq_resid == startres) {
          cyclic=1;
          done=1;
          break;
        }
      }  
    }

    // record whether this fragment is cyclic or not
    subfragCyclic->append(cyclic);
  }
}


// this adds the current residue type to the *subfragList,
// this finds the residue connected to the endatom atom type
// and calls this function recursively on that residue
// this will NOT work across NORMAL bonds
void BaseMolecule::find_connected_subfragment(int resnum, int fragnum, 
         char *flgs, int endatom,  int altendatom, 
         int alt2endatom, int alt3endatom,
         int restype, 
         ResizeArray<Fragment *> *subfragList)
{
  if (flgs[resnum] || residueList[resnum]->residueType != restype) 
      return;  // been here before, or this is no good
  (*subfragList)[fragnum]->append(resnum);    // add to the list
  flgs[resnum] = TRUE;                        // and prevent repeats

  // find the atom in this residue with the right type
  int i, j, nextres;
  MolAtom *a;
  for (i=residueList[resnum]->atoms.num() - 1; i>=0; i--) {
    a = atom(residueList[resnum]->atoms[i]);
    if (a->nameindex == endatom ||
        a->nameindex == altendatom ||
        a->nameindex == alt2endatom ||
        a->nameindex == alt3endatom) {   // found the end atom
      for (j=a->bonds-1; j>=0; j--) {    // look at the bonds
        // I can't look at if the residue "bond" is a PRO-PRO or NUC-NUC, since
        // that won't tell me if the atom involved is the endatom atom
        // This is important because I need to avoid things like S-S bonds
        // (note that I never checked that the end was bonded to a start on
        //  the next residue! - c'est la vie, or something like that
        if ((!(a->atomType == ATOMNORMAL && atom(a->bondTo[j])->atomType == ATOMNORMAL)) && // not backbone 
            (nextres = atom(a->bondTo[j])->uniq_resid) != resnum &&
            !flgs[nextres] ) { // found next residue, and unvisited
          find_connected_subfragment(nextres, fragnum, flgs, endatom,
              altendatom, alt2endatom, alt3endatom, restype, subfragList);
          return; // only find one; assume no branching
        }
      } // end of for
    } // end of found correct endtype
  } // searching atoms
} // end of finding connected subfragment


// find a class of fragments, and add them to the subfragment list
void BaseMolecule::find_subfragments(int startatom, 
          int altstartatom, int alt2startatom,
          int endatom, int altendatom, int alt2endatom, int alt3endatom,
          int restype, ResizeArray<Fragment *> *subfragList)
{
  int i, j, k;
  MolAtom *a;
  char *flgs = new char[residueList.num()];
  memset(flgs, 0, residueList.num() * sizeof(char));  // clear flags

  // Loop over all residues looking for candidate residues that start
  // a fragment.  A fragment starting residue must be an unvisited 
  // residue which has an startatom with no off residue bond to 
  // the same restype
  for (i=residueList.num()-1; i>=0; i--) {
    // test for previous visit, and whether it's the restype we want
    if (!flgs[i] && residueList[i]->residueType == restype) {
      // does this residue have a matching startatom
      for (j=residueList[i]->atoms.num()-1; j>=0; j--) { 
        int satom = (a=atom(residueList[i]->atoms[j]))->nameindex;
        if (satom == startatom || 
            satom == altstartatom || 
            satom == alt2startatom){
          for (k=a->bonds-1; k>=0; k--) {
            MolAtom *bondto = atom(a->bondTo[k]);
            // are there any off-residue bonds to the same restype
            if (bondto->uniq_resid != i && bondto->residueType == restype) {
              break; // if so then stop, so that k>=0
            }
          }

          // if we found a valid fragment start atom, find residues downchain
          if (k<0) { 
            subfragList->append(new Fragment);
            find_connected_subfragment(i, subfragList->num()-1, flgs, 
                  endatom, altendatom, alt2endatom, alt3endatom,
                  restype, subfragList);
          } // found starting residue
        } // found startatom
      } // going through atoms
    } // found restype
  } // going through residues

  // found 'em all
  delete [] flgs;
} 


// find a class of fragments, and add them to the subfragment list
void BaseMolecule::find_subfragments_topologically(int restype, 
  ResizeArray<Fragment *> *subfragList, 
  int endatom, int altendatom, int alt2endatom, int alt3endatom) {
  int i; 
  char *flgs = new char[residueList.num()];
  memset(flgs, 0, residueList.num() * sizeof(char));  // clear flags
  int numres = residueList.num();

  // Loop over all residues looking for candidate residues that start
  // a fragment.  A fragment starting residue must be an unvisited
  // residue which has an startatom with no off residue bond to
  // the same restype
  for (i=0; i<numres; i++) {
    Residue *res = residueList[i];

    // test for previous visit, and whether it's the restype we want
    if (!flgs[i] && res->residueType == restype) {
      // if this residue only has 1 bond to a residue of the same restype
      // it must be a terminal residue
      int offresbondcount = 0;
      int j, k;
      int numatoms = res->atoms.num();
      for (j=0; j<numatoms; j++) {
        MolAtom *a = atom(res->atoms[j]);

        // find off-residue bonds
        for (k=0; k<a->bonds; k++) {
          MolAtom *b = atom(a->bondTo[k]);
          if (b->uniq_resid != i && 
              residueList[b->uniq_resid]->residueType == restype) {
            offresbondcount++;
          }
        }
      }

      // if we found a valid fragment start atom, find residues downchain
      if (offresbondcount == 1) {
        subfragList->append(new Fragment);
        find_connected_subfragment(i, subfragList->num()-1, flgs,
              endatom, altendatom, alt2endatom, alt3endatom,
              restype, subfragList);
      }
    } // found restype
  } // going through residues

  // found 'em all
  delete [] flgs;
}


#if defined(VMDFASTRIBBONS)
// routine to pre-calculate lists of control point indices to
// be used by the ribbon/cartoon representations
void BaseMolecule::calculate_ribbon_controlpoints() {
  int onum, canum, frag, num, res;

  vmd_timerhandle tm = vmd_timer_create();

  msgInfo << "XXX Building pre-computed control point lists" << sendmsg;

  vmd_timer_start(tm);

  // Lookup atom typecodes ahead of time so we can use the most
  // efficient variant of find_atom_in_residue() in the main loop.
  // If we can't find the atom types we need, bail out immediately
  int CAtypecode  = atomNames.typecode("CA");
  int Otypecode   = atomNames.typecode("O");
  int OT1typecode = atomNames.typecode("OT1");
   
  // We can't draw a ribbon without CA and O atoms for guidance!
  if (CAtypecode < 0) {
    msgErr << "This structure has no identifiable CA atoms,"
           << "ribbon and cartoon representations will not be usable." 
           << sendmsg;
    return;
  }
  if ((Otypecode < 0) && (OT1typecode < 0)) {
    msgErr << "This structure has no identifiable Oxygen atoms,"
           << "ribbon and cartoon representations will not be usable." 
           << sendmsg;
    return;
  }

  //
  // calculate control point lists for the protein fragments
  //

  msgInfo << "XXX Building protein control point lists" << sendmsg;

  // go through each protein and find the CA and O atoms which are
  // eventually used to construct control points and perpendicular vectors
  ResizeArray<int> cpind;
  for (frag=0; frag<pfragList.num(); frag++) {
    num = pfragList[frag]->num();  // number of residues in this fragment
    if (num < 2) {
      // with less than two residues, we can't do anything useful, so skip
      nfragCPList.append(NULL);
      continue;
    }

    // check that we have a valid structure before continuing
      res = (*pfragList[fra