---

AtomSel.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: AtomSel.C,v $
 *      $Author: johns $        $Locker:  $                $State: Exp $
 *      $Revision: 1.157 $      $Date: 2008/03/27 19:36:35 $
 *
 ***************************************************************************
 * DESCRIPTION:
 * 
 * Parse and maintain the data for selecting atoms.
 *
 ***************************************************************************/

#include <math.h>
#include <stdlib.h>
#include <stdio.h>
#include <ctype.h>
#include <string.h>

#if defined(ARCH_AIX4)
#include <strings.h>
#endif

#include "Atom.h"
#include "AtomSel.h"
#include "DrawMolecule.h"
#include "MoleculeList.h"
#include "VMDApp.h"
#include "Inform.h"
#include "SymbolTable.h"
#include "ParseTree.h"
#include "JRegex.h"
#include "VolumetricData.h"
#include "PeriodicTable.h"

// 'all'
static int atomsel_all(void * ,int, int *) {
  return 1;
}

// 'none'
static int atomsel_none(void *, int num, int *flgs) {
  for (int i=num-1; i>=0; i--) {
    flgs[i] = FALSE;
  }
  return 1;
}

#define generic_atom_data(fctnname, datatype, field)                   \
static int fctnname(void *v, int num, datatype *data, int *flgs) {     \
  DrawMolecule *atom_sel_mol = ((atomsel_ctxt *)v)->atom_sel_mol;      \
  for (int i=0; i<num; i++) {                                          \
    if (flgs[i]) {                                                     \
      data[i] = atom_sel_mol->atom(i)->field;                          \
    }                                                                  \
  }                                                                    \
  return 1;                                                            \
}


#define generic_set_atom_data(fctnname, datatype, fieldtype, field)    \
static int fctnname(void *v, int num, datatype *data, int *flgs) {     \
  DrawMolecule *atom_sel_mol = ((atomsel_ctxt *)v)->atom_sel_mol;      \
  int i;                                                               \
  for (i=0; i<num; i++) {                                              \
    if (flgs[i]) {                                                     \
      atom_sel_mol->atom(i)->field = (fieldtype) data[i];              \
    }                                                                  \
  }                                                                    \
  return 1;                                                            \
}


// 'name'
static int atomsel_name(void *v, int num, const char **data, int *flgs) {
  DrawMolecule *atom_sel_mol = ((atomsel_ctxt *)v)->atom_sel_mol;
  for (int i=0; i<num; i++) {
    if (flgs[i])
      data[i] = (atom_sel_mol->atomNames).name(
                  atom_sel_mol->atom(i)->nameindex);
  }
  return 1;
}


// 'type'
static int atomsel_type(void *v, int num, const char **data, int *flgs) {
  DrawMolecule *atom_sel_mol = ((atomsel_ctxt *)v)->atom_sel_mol;
  for (int i=0; i<num; i++) {
    if (flgs[i])
      data[i] = (atom_sel_mol->atomTypes).name(
                  atom_sel_mol->atom(i)->typeindex);
  }
  return 1;
}

// XXX probably only use this for internal testing
// 'backbonetype'
static int atomsel_backbonetype(void *v, int num, const char **data, int *flgs) {
  DrawMolecule *atom_sel_mol = ((atomsel_ctxt *)v)->atom_sel_mol;
  for (int i=0; i<num; i++) {
    if (flgs[i]) {
      switch (atom_sel_mol->atom(i)->atomType) {
        case ATOMNORMAL:
          data[i] = "normal";
          break;

        case ATOMPROTEINBACK:
          data[i] = "proteinback";
          break;

        case ATOMNUCLEICBACK:
          data[i] = "nucleicback";
          break;

        case ATOMHYDROGEN:
          data[i] = "hydrogen";
          break;

        default:
          data[i] = "unknown";
          break;
      }
    }
  }
  return 1;
}


// XXX probably only use this for internal testing
// 'residuetype'
static int atomsel_residuetype(void *v, int num, const char **data, int *flgs) {
  DrawMolecule *atom_sel_mol = ((atomsel_ctxt *)v)->atom_sel_mol;
  for (int i=0; i<num; i++) {
    if (flgs[i]) {
      switch (atom_sel_mol->atom(i)->residueType) {
        case RESNOTHING:
          data[i] = "nothing";
          break;

        case RESPROTEIN:
          data[i] = "protein";
          break;

        case RESNUCLEIC:
          data[i] = "nucleic";
          break;

        case RESWATERS:
          data[i] = "waters";
          break;

        default:
          data[i] = "unknown";
          break;
      }
    }
  }
  return 1;
}


// 'atomicnumber'
generic_atom_data(atomsel_atomicnumber, int, atomicnumber)
generic_set_atom_data(atomsel_set_atomicnumber, int, int, atomicnumber)


// 'element'
static int atomsel_element(void *v, int num, const char **data, int *flgs) {
  DrawMolecule *atom_sel_mol = ((atomsel_ctxt *)v)->atom_sel_mol;
  for (int i=0; i<num; i++) {
    if (flgs[i])
      data[i] = get_pte_label(atom_sel_mol->atom(i)->atomicnumber);
  }
  return 1;
}
static int atomsel_set_element(void *v, int num, const char **data, int *flgs) {
  DrawMolecule *atom_sel_mol = ((atomsel_ctxt *)v)->atom_sel_mol;
  for (int i=0; i<num; i++)
    if (flgs[i]) {
      int idx = get_pte_idx_from_string((const char *)(data[i]));
      atom_sel_mol->atom(i)->atomicnumber = idx;
    }
  return 1;
}


// 'index'
static int atomsel_index(void *v, int num, int *data, int *flgs) { 
  for (int i=0; i<num; i++) {
    if (flgs[i]) {
      data[i] = i; // zero-based
    }
  }
  return 1;
}


// 'serial' (one-based atom index)
static int atomsel_serial(void *v, int num, int *data, int *flgs) { 
  for (int i=0; i<num; i++) {
    if (flgs[i]) {
      data[i] = i + 1; // one-based
    }
  }
  return 1;
}


// 'fragment'
static int atomsel_fragment(void *v, int num, int *data, int *flgs) {
  DrawMolecule *atom_sel_mol = ((atomsel_ctxt *)v)->atom_sel_mol;
  for (int i=0; i<num; i++) {
    if (flgs[i]) {
      int residue = atom_sel_mol->atom(i)->uniq_resid;
      data[i] = (atom_sel_mol->residue(residue))->fragment;
    }
  }
  return 1;
}


// 'numbonds'
generic_atom_data(atomsel_numbonds, int, bonds)


// 'residue'
generic_atom_data(atomsel_residue, int, uniq_resid)


// 'resname'
static int atomsel_resname(void *v, int num, const char **data, int *flgs) {
  DrawMolecule *atom_sel_mol = ((atomsel_ctxt *)v)->atom_sel_mol;
  for (int i=0; i<num; i++) {
    if (flgs[i])
      data[i] = (atom_sel_mol->resNames).name(
                 atom_sel_mol->atom(i)->resnameindex);
  }
  return 1;
}


// 'altloc'
static int atomsel_altloc(void *v, int num, const char **data, int *flgs) {
  DrawMolecule *atom_sel_mol = ((atomsel_ctxt *)v)->atom_sel_mol;                 
  for (int i=0; i<num; i++) {
    if (flgs[i])
      data[i] = (atom_sel_mol->altlocNames).name(
        atom_sel_mol->atom(i)->altlocindex);
  }
  return 1;
}
static int atomsel_set_altloc(void *v, int num, const char **data, int *flgs) {
  DrawMolecule *atom_sel_mol = ((atomsel_ctxt *)v)->atom_sel_mol;                 
  NameList<int> &arr = atom_sel_mol->altlocNames;
  for (int i=0; i<num; i++)
    if (flgs[i]) {
      int ind = arr.add_name((const char *)(data[i]), arr.num());
      atom_sel_mol->atom(i)->altlocindex = ind;
    }
  return 1;
}


// 'insertion'
generic_atom_data(atomsel_insertion, const char *, insertionstr)


// 'chain'
static int atomsel_chain(void *v, int num, const char **data, int *flgs) {
  DrawMolecule *atom_sel_mol = ((atomsel_ctxt *)v)->atom_sel_mol;                 
  for (int i=0; i<num; i++) {
    if (flgs[i])
      data[i] = (atom_sel_mol->chainNames).name(
                  atom_sel_mol->atom(i)->chainindex);
  }
  return 1;
}


// 'segname'
static int atomsel_segname(void *v, int num, const char **data, int *flgs) {
  DrawMolecule *atom_sel_mol = ((atomsel_ctxt *)v)->atom_sel_mol;
  for (int i=0; i<num; i++) {
    if (flgs[i])
      data[i] = (atom_sel_mol->segNames).name(
                 atom_sel_mol->atom(i)->segnameindex);
  }
  return 1;
}


// The next few set functions affect the namelists kept in BaseMolecule
static int atomsel_set_name(void *v, int num, const char **data, int *flgs) {
  DrawMolecule *atom_sel_mol = ((atomsel_ctxt *)v)->atom_sel_mol;
  NameList<int> &arr = atom_sel_mol->atomNames;
  for (int i=0; i<num; i++) {
    if (flgs[i]) {
      int ind = arr.add_name((const char *)(data[i]), arr.num());
      atom_sel_mol->atom(i)->nameindex = ind;
    }
  }
  return 1;
}

static int atomsel_set_type(void *v, int num, const char **data, int *flgs) {
  DrawMolecule *atom_sel_mol = ((atomsel_ctxt *)v)->atom_sel_mol;
  NameList<int> &arr = atom_sel_mol->atomTypes;
  for (int i=0; i<num; i++) {
    if (flgs[i]) {
      int ind = arr.add_name((const char *)(data[i]), arr.num());
      atom_sel_mol->atom(i)->typeindex = ind;
    }
  }
  return 1;
}

static int atomsel_set_resname(void *v, int num, const char **data, int *flgs) {
  DrawMolecule *atom_sel_mol = ((atomsel_ctxt *)v)->atom_sel_mol;
  NameList<int> &arr = atom_sel_mol->resNames;
  for (int i=0; i<num; i++) {
    if (flgs[i]) {
      int ind = arr.add_name((const char *)(data[i]), arr.num());
      atom_sel_mol->atom(i)->resnameindex = ind;
    }
  }
  return 1;
}

static int atomsel_set_chain(void *v, int num, const char **data, int *flgs) {
  DrawMolecule *atom_sel_mol = ((atomsel_ctxt *)v)->atom_sel_mol;
  NameList<int> &arr = atom_sel_mol->chainNames;
  for (int i=0; i<num; i++) {
    if (flgs[i]) {
      int ind = arr.add_name((const char *)(data[i]), arr.num());
      atom_sel_mol->atom(i)->chainindex = ind;
    }
  }
  return 1;
}

static int atomsel_set_segname(void *v, int num, const char **data, int *flgs) {
  DrawMolecule *atom_sel_mol = ((atomsel_ctxt *)v)->atom_sel_mol;
  NameList<int> &arr = atom_sel_mol->segNames;
  for (int i=0; i<num; i++) {
    if (flgs[i]) {
      int ind = arr.add_name((const char *)(data[i]), arr.num());
      atom_sel_mol->atom(i)->segnameindex = ind;
    }
  }
  return 1;
}


// 'radius'
static int atomsel_radius(void *v, int num, double *data, int *flgs) {
  DrawMolecule *atom_sel_mol = ((atomsel_ctxt *)v)->atom_sel_mol;
  float *radius = atom_sel_mol->radius();
  for (int i=0; i<num; i++)
    if (flgs[i]) data[i] = radius[i];
  return 1;
}
static int atomsel_set_radius(void *v, int num, double *data, int *flgs) {
  DrawMolecule *atom_sel_mol = ((atomsel_ctxt *)v)->atom_sel_mol;
  float *radius = atom_sel_mol->radius();
  for (int i=0; i<num; i++)
    if (flgs[i]) radius[i] = (float) data[i];
  return 1;
}


// 'mass'
static int atomsel_mass(void *v, int num, double *data, int *flgs) {
  DrawMolecule *atom_sel_mol = ((atomsel_ctxt *)v)->atom_sel_mol;
  float *mass = atom_sel_mol->mass();
  for (int i=0; i<num; i++)
    if (flgs[i]) data[i] = mass[i];
  return 1;
}
static int atomsel_set_mass(void *v, int num, double *data, int *flgs) {
  DrawMolecule *atom_sel_mol = ((atomsel_ctxt *)v)->atom_sel_mol;
  float *mass = atom_sel_mol->mass();
  for (int i=0; i<num; i++)
    if (flgs[i]) mass[i] = (float) data[i];
  return 1;
}


// 'charge'
static int atomsel_charge(void *v, int num, double *data, int *flgs) {
  DrawMolecule *atom_sel_mol = ((atomsel_ctxt *)v)->atom_sel_mol;
  float *charge = atom_sel_mol->charge();
  for (int i=0; i<num; i++)
    if (flgs[i]) data[i] = charge[i];
  return 1;
}
static int atomsel_set_charge(void *v, int num, double *data, int *flgs) {
  DrawMolecule *atom_sel_mol = ((atomsel_ctxt *)v)->atom_sel_mol;
  float *charge = atom_sel_mol->charge();
  for (int i=0; i<num; i++)
    if (flgs[i]) charge[i] = (float) data[i];
  return 1;
}


// 'beta'
static int atomsel_beta(void *v, int num, double *data, int *flgs) {
  DrawMolecule *atom_sel_mol = ((atomsel_ctxt *)v)->atom_sel_mol;
  float *beta = atom_sel_mol->beta();
  for (int i=0; i<num; i++)
    if (flgs[i]) data[i] = beta[i];
  return 1;
}
static int atomsel_set_beta(void *v, int num, double *data, int *flgs) {
  DrawMolecule *atom_sel_mol = ((atomsel_ctxt *)v)->atom_sel_mol;
  float *beta = atom_sel_mol->beta();
  for (int i=0; i<num; i++)
    if (flgs[i]) beta[i] = (float) data[i];
  return 1;
}


// 'occupancy?'
static int atomsel_occupancy(void *v, int num, double *data, int *flgs) {
  DrawMolecule *atom_sel_mol = ((atomsel_ctxt *)v)->atom_sel_mol;
  float *occupancy = atom_sel_mol->occupancy();
  for (int i=0; i<num; i++)
    if (flgs[i]) data[i] = occupancy[i];
  return 1;
}
static int atomsel_set_occupancy(void *v, int num, double *data, int *flgs) {
  DrawMolecule *atom_sel_mol = ((atomsel_ctxt *)v)->atom_sel_mol;
  float *occupancy = atom_sel_mol->occupancy();
  for (int i=0; i<num; i++)
    if (flgs[i]) occupancy[i] = (float) data[i];
  return 1;
}


// 'resid'
static int atomsel_resid(void *v, int num, int *data, int *flgs) {
  DrawMolecule *atom_sel_mol = ((atomsel_ctxt *)v)->atom_sel_mol;
  for (int i=0; i<num; i++) {
    if (flgs[i]) {
      data[i] = atom_sel_mol->atom(i)->resid;
    }
  }
  return 1;
}
static int atomsel_set_resid(void *v, int num, int *data, int *flgs) {
  DrawMolecule *atom_sel_mol = ((atomsel_ctxt *)v)->atom_sel_mol;
  for (int i=0; i<num; i++) {
    if (flgs[i])
      atom_sel_mol->atom(i)->resid = data[i];
  }
  return 1;
}



#define generic_atom_boolean(fctnname, comparison)                    \
static int fctnname(void *v, int num, int *flgs) {                    \
  DrawMolecule *atom_sel_mol = ((atomsel_ctxt *)v)->atom_sel_mol;     \
  for (int i=0; i<num; i++) {                                         \
    if (flgs[i]) {                                                    \
      flgs[i] = atom_sel_mol->atom(i)->comparison;                    \
    }                                                                 \
  }                                                                   \
  return 1;                                                           \
}

// 'backbone'
static int atomsel_backbone(void *v, int num, int *flgs) {
  DrawMolecule *atom_sel_mol = ((atomsel_ctxt *)v)->atom_sel_mol;
  for (int i=0; i<num; i++) {
    if (flgs[i]) {
      const MolAtom *a = atom_sel_mol->atom(i);
      flgs[i] = (a->atomType == ATOMPROTEINBACK ||
                 a->atomType == ATOMNUCLEICBACK);
    }
  }
  return 1;
}

// 'h' (hydrogen)
generic_atom_boolean(atomsel_hydrogen, atomType == ATOMHYDROGEN)


// 'protein'
generic_atom_boolean(atomsel_protein, residueType == RESPROTEIN)


// 'nucleic'
generic_atom_boolean(atomsel_nucleic, residueType == RESNUCLEIC)


// 'water'
generic_atom_boolean(atomsel_water, residueType == RESWATERS)

#define generic_sstruct_boolean(fctnname, comparison)                    \
static int fctnname(void *v, int num, int *flgs) {                       \
  DrawMolecule *atom_sel_mol = ((atomsel_ctxt *)v)->atom_sel_mol;        \
  atom_sel_mol->need_secondary_structure(1);                             \
  for (int i=0; i<num; i++) {                                            \
    int s;                                                               \
    if (flgs[i]) {                                                       \
      s = atom_sel_mol->residue(                                         \
                                  atom_sel_mol->atom(i)->uniq_resid      \
                                  )->sstruct;                            \
      if (!comparison) {                                                 \
        flgs[i] = 0;                                                     \
      }                                                                  \
    }                                                                    \
  }                                                                      \
  return 1;                                                              \
}

// once I define a structure, I don't need to recompute; hence the 
//   need_secondary_structure(0);
#define generic_set_sstruct_boolean(fctnname, newval)                    \
static int fctnname(void *v, int num, int *data, int *flgs) {            \
  DrawMolecule *atom_sel_mol = ((atomsel_ctxt *)v)->atom_sel_mol;        \
  atom_sel_mol->need_secondary_structure(0);                             \
  for (int i=0; i<num; i++) {                                            \
    if (flgs[i] && data[i]) {                                            \
        atom_sel_mol->residue(atom_sel_mol->atom(i)->uniq_resid)         \
          ->sstruct = newval;                                            \
    }                                                                    \
  }                                                                      \
  return 1;                                                              \
}


// XXX recursive routines should be replaced by an iterative version with
// it's own stack, so that huge molecules don't overflow the stack
static void recursive_find_sidechain_atoms(BaseMolecule *mol, int *sidechain,
                                           int atom_index) {
  // Have I been here before?
  if (sidechain[atom_index] == 2) 
    return;

  // Is this a backbone atom
  MolAtom *atom = mol->atom(atom_index);
  if (atom->atomType == ATOMPROTEINBACK ||
      atom->atomType == ATOMNUCLEICBACK) return;
  
  // mark this atom
  sidechain[atom_index] = 2;

  // try the atoms to which this is bonded
  for (int i=0; i<atom->bonds; i++) {
    recursive_find_sidechain_atoms(mol, sidechain, atom->bondTo[i]);
  }
}

// give an array where 1 indicates an atom on the sidechain, find the
// connected atoms which are also on the sidechain
static void find_sidechain_atoms(BaseMolecule *mol, int *sidechain) {
  for (int i=0; i<mol->nAtoms; i++) {
    if (sidechain[i] == 1) {
      recursive_find_sidechain_atoms(mol, sidechain, i);
    }
  }
}


// 'sidechain' is tricky.  I start from a protein CA, pick a bond which
// isn't along the backbone (and discard it once if it is a hydrogen),
// and follow the atoms until I stop at another backbone atom or run out
// of atoms
static int atomsel_sidechain(void *v, int num, int *flgs) {
  DrawMolecule *atom_sel_mol = ((atomsel_ctxt *)v)->atom_sel_mol;               
  const float *mass = atom_sel_mol->mass();
  int *seed = new int[num];
  int i;

  // generate a list of the "CB" atoms (or whatever they are)
  for (i=0; i<num; i++) {
    seed[i] = 0;
  }

  // get the CA and HA2 name index
  int CA = atom_sel_mol->atomNames.typecode((char *) "CA");

  // for each protein
  for (int pfrag=atom_sel_mol->pfragList.num()-1; pfrag>=0; pfrag--) {
    // get a residue
    Fragment &frag = *(atom_sel_mol->pfragList[pfrag]);
    for (int res = frag.num()-1; res >=0; res--) {
      // find the CA
      int atomid = atom_sel_mol->find_atom_in_residue(CA, frag[res]);
      if (atomid < 0) {
        msgErr << "atomselection: sidechain: cannot find CA" << sendmsg;
        continue;
      }
      // find at most two neighbors which are not on the backbone
      MolAtom *atom = atom_sel_mol->atom(atomid);
      int b1 = -1, b2 = -1;
      for (i=0; i<atom->bonds; i++) {
        int bondtype = atom_sel_mol->atom(atom->bondTo[i])->atomType;
        if (bondtype == ATOMNORMAL || bondtype == ATOMHYDROGEN) {
          if (b1 == -1) {
            b1 = atom->bondTo[i];
          } else {
            if (b2 == -1) {
              b2 = atom->bondTo[i];
            } else {
              msgErr << "atomselection: sidechain: protein residue index " 
                     << res << ", C-alpha index " << i << " has more than "
                     << "two non-backbone bonds; ignoring the others"
                     << sendmsg;
            }
          }
        }
      }
      if (b1 == -1) 
        continue;

      if (b2 != -1) {  // find the right one
        // first check the number of atoms and see if I have a lone H
        int c1 = atom_sel_mol->atom(b1)->bonds;
        int c2 = atom_sel_mol->atom(b2)->bonds;
        if (c1 == 1 && c2 > 1) {
          b1 = b2;
        } else if (c2 == 1 && c1 > 1) {
          b1 = b1;
        } else if (c1 ==1 && c2 == 1) {
          // check the masses
          float m1 = mass[b1];
          float m2 = mass[b2];
          if (m1 > 2.3 && m2 <= 2.3) {
            b1 = b2;
          } else if (m2 > 2.3 && m1 <= 2.3) {
            b1 = b1;
          } else if (m1 <= 2.0 && m2 <= 2.3) {
            // should have two H's, find the "first" of these
            if (strcmp(
              (atom_sel_mol->atomNames).name(atom_sel_mol->atom(b1)->nameindex),
              (atom_sel_mol->atomNames).name(atom_sel_mol->atom(b2)->nameindex)
              ) > 0) {
              b1 = b2;
            } // else b1 = b1
          } else {
            msgErr << "atomselect: sidechain:  protein residue index " 
                   << res << ", C-alpha index " << i << " has sidechain-like "
                   << "atom (indices " << b1 << " and " << b2 << "), and "
                   << "I cannot determine which to call a sidechain -- "
                   << "I'll guess" << sendmsg;
            if (strcmp(
              (atom_sel_mol->atomNames).name(atom_sel_mol->atom(b1)->nameindex),
              (atom_sel_mol->atomNames).name(atom_sel_mol->atom(b2)->nameindex)
              ) > 0) {
              b1 = b2;
            }
          } // punted
        } // checked connections and masses
      } // found the right one; it is b1.
      seed[b1] = 1;

    } // loop over residues
  } // loop over protein fragments

  // do the search for all the sidechain atoms (returned in seed)
  find_sidechain_atoms(atom_sel_mol, seed);

  // set the return values
  for (i=0; i<num; i++) {
    if (flgs[i]) 
      flgs[i] = (seed[i] != 0);
  }

  delete [] seed;

  return 1;
}


// which of these are helices?
generic_sstruct_boolean(atomsel_helix, (s == SS_HELIX_ALPHA ||
                                        s == SS_HELIX_3_10  ||
                                        s == SS_HELIX_PI))
generic_sstruct_boolean(atomsel_alpha_helix, (s == SS_HELIX_ALPHA))
generic_sstruct_boolean(atomsel_3_10_helix, (s == SS_HELIX_3_10))
generic_sstruct_boolean(atomsel_pi_helix, (s == SS_HELIX_PI))

// Makes the residue into a HELIX_ALPHA
generic_set_sstruct_boolean(atomsel_set_helix, SS_HELIX_ALPHA)
// Makes the residue into a HELIX_3_10
generic_set_sstruct_boolean(atomsel_set_3_10_helix, SS_HELIX_3_10)
// Makes the residue into a HELIX_PI
generic_set_sstruct_boolean(atomsel_set_pi_helix, SS_HELIX_PI)

// which of these are beta sheets?
generic_sstruct_boolean(atomsel_sheet, (s == SS_BETA ||
                                        s == SS_BRIDGE ))
generic_sstruct_boolean(atomsel_extended_sheet, (s == SS_BETA))
generic_sstruct_boolean(atomsel_bridge_sheet, (s == SS_BRIDGE ))

// Makes the residue into a BETA
generic_set_sstruct_boolean(atomsel_set_sheet, SS_BETA)
generic_set_sstruct_boolean(atomsel_set_bridge_sheet, SS_BRIDGE)

// which of these are coils?
generic_sstruct_boolean(atomsel_coil, (s == SS_COIL))

// Makes the residue into COIL
generic_set_sstruct_boolean(atomsel_set_coil, SS_COIL)

// which of these are TURNS?
generic_sstruct_boolean(atomsel_turn, (s == SS_TURN))

// Makes the residue into TURN
generic_set_sstruct_boolean(atomsel_set_turn, SS_TURN)

// return the sstruct information as a 1 character string
static int atomsel_sstruct(void *v, int num, const char **data, int *flgs) {
  DrawMolecule *atom_sel_mol = ((atomsel_ctxt *)v)->atom_sel_mol;
  atom_sel_mol->need_secondary_structure(1);
  for (int i=0; i<num; i++) {
    if (flgs[i]) {
      switch(atom_sel_mol->residue(atom_sel_mol->atom(i)->uniq_resid)->sstruct) {
        case SS_HELIX_ALPHA: data[i] = "H"; break;
        case SS_HELIX_3_10 : data[i] = "G"; break;
        case SS_HELIX_PI   : data[i] = "I"; break;
        case SS_BETA       : data[i] = "E"; break;
        case SS_BRIDGE     : data[i] = "B"; break;
        case SS_TURN       : data[i] = "T"; break;
        default:
        case SS_COIL       : data[i] = "C"; break;
      }
    }
  }
  return 1;
}


// set the secondary structure based on a string value
static int atomsel_set_sstruct(void *v, int num, const char **data, int *flgs) {
  DrawMolecule *atom_sel_mol = ((atomsel_ctxt *)v)->atom_sel_mol;
  char c;
  // Defining it here so remind myself that I don't need STRIDE (or
  // whoever) to do it automatically
  atom_sel_mol->need_secondary_structure(0);
  for (int i=0; i<num; i++) {
    if (flgs[i]) {
      if (strlen(data[i]) == 0) {
        msgErr << "cannot set a 0 length secondary structure string"
               << sendmsg;
      } else {
        c = ((const char *) data[i])[0];
        if (strlen(data[i]) > 1) {
          while (1) {
 if (!strcasecmp((const char *) data[i], "helix")) { c = 'H'; break;}
 if (!strcasecmp((const char *) data[i], "alpha")) { c = 'H'; break;}
 if (!strcasecmp((const char *) data[i], "alpha_helix")) { c = 'H'; break;}
 if (!strcasecmp((const char *) data[i], "alphahelix"))  { c = 'H'; break;}
 if (!strcasecmp((const char *) data[i], "alpha helix")) { c = 'H'; break;}

 if (!strcasecmp((const char *) data[i], "pi"))        { c = 'I'; break;}
 if (!strcasecmp((const char *) data[i], "pi_helix"))  { c = 'I'; break;}
 if (!strcasecmp((const char *) data[i], "pihelix"))   { c = 'I'; break;}
 if (!strcasecmp((const char *) data[i], "pi helix"))  { c = 'I'; break;}

 if (!strcasecmp((const char *) data[i], "310"))       { c = 'G'; break;}
 if (!strcasecmp((const char *) data[i], "310_helix")) { c = 'G'; break;}
 if (!strcasecmp((const char *) data[i], "3_10"))      { c = 'G'; break;}
 if (!strcasecmp((const char *) data[i], "3"))         { c = 'G'; break;}
 if (!strcasecmp((const char *) data[i], "310 helix")) { c = 'G'; break;}
 if (!strcasecmp((const char *) data[i], "3_10_helix")){ c = 'G'; break;}
 if (!strcasecmp((const char *) data[i], "3_10 helix")){ c = 'G'; break;}
 if (!strcasecmp((const char *) data[i], "3 10 helix")){ c = 'G'; break;}
 if (!strcasecmp((const char *) data[i], "helix_3_10")){ c = 'G'; break;}
 if (!strcasecmp((const char *) data[i], "helix 3 10")){ c = 'G'; break;}
 if (!strcasecmp((const char *) data[i], "helix3_10")) { c = 'G'; break;}
 if (!strcasecmp((const char *) data[i], "helix310"))  { c = 'G'; break;}

 if (!strcasecmp((const char *) data[i], "beta"))      { c = 'E'; break;}
 if (!strcasecmp((const char *) data[i], "betasheet")) { c = 'E'; break;}
 if (!strcasecmp((const char *) data[i], "beta_sheet")){ c = 'E'; break;}
 if (!strcasecmp((const char *) data[i], "beta sheet")){ c = 'E'; break;}
 if (!strcasecmp((const char *) data[i], "sheet"))     { c = 'E'; break;}
 if (!strcasecmp((const char *) data[i], "strand"))    { c = 'E'; break;}
 if (!strcasecmp((const char *) data[i], "beta_strand"))  { c = 'E'; break;}
 if (!strcasecmp((const char *) data[i], "beta strand"))  { c = 'E'; break;}

 if (!strcasecmp((const char *) data[i], "turn"))  { c = 'T'; break;}

 if (!strcasecmp((const char *) data[i], "coil"))     { c = 'C'; break;}
 if (!strcasecmp((const char *) data[i], "unknown"))  { c = 'C'; break;}
 c = 'X';
 break;
          } // while (1)
        }
        // and set the value
        int s = SS_COIL;
        switch ( c ) {
        case 'H':
        case 'h': s = SS_HELIX_ALPHA; break;
        case '3':
        case 'G':
        case 'g': s = SS_HELIX_3_10; break;
        case 'P':  // so you can say 'pi'
        case 'p':
        case 'I':
        case 'i': s = SS_HELIX_PI; break;
        case 'S':  // for "sheet"
        case 's':
        case 'E':
        case 'e': s = SS_BETA; break;
        case 'B':
        case 'b': s = SS_BRIDGE; break;
        case 'T':
        case 't': s = SS_TURN; break;
        case 'L': // L is from PHD and may be an alternate form
        case 'l':
        case 'C':
        case 'c': s = SS_COIL; break;
        default: {
          msgErr << "Unknown sstruct assignment: '"
            << (const char *) data[i] << "'" << sendmsg;
          s = SS_COIL; break;
        }
        }
        atom_sel_mol->residue(atom_sel_mol->atom(i)->uniq_resid)->sstruct = s;
      }
    }
  }
  return 1;
}


//      and leave the rest alone.
// It is slower this way, but easier to understand
static void mark_atoms_given_residue(DrawMolecule *mol, int residue, int *on)
{
  ResizeArray<int> *atoms = &(mol->residueList[residue]->atoms);
  for (int i= atoms->num()-1; i>=0; i--) {
     on[(*atoms)[i]] = TRUE;
  }
}


// macro for either protein or nucleic fragments
#define fragment_data(fctn, fragtype)                                         \
static int fctn(void *v, int num, int *data, int *)                           \
{                                                                             \
  DrawMolecule *atom_sel_mol = ((atomsel_ctxt *)v)->atom_sel_mol;                     \
   int *tmp = new int[num];                                                   \
   int i;                                                                     \
   for (i=num-1; i>=0; i--) {  /* clear the arrays */                 \
      tmp[i] = 0;                                                             \
      data[i] = -1;  /* default is -1 for 'not a [np]frag' */                 \
   }                                                                          \
   /* for each fragment */                                                    \
   for ( i=atom_sel_mol->fragtype.num()-1; i>=0; i--) {                       \
      /* for each residues of the fragment */                                 \
      int j;                                                                  \
      for (j=atom_sel_mol->fragtype[i]->num()-1; j>=0; j--) {                 \
         /* mark the atoms in the fragment */                                 \
         mark_atoms_given_residue(atom_sel_mol,(*atom_sel_mol->fragtype[i])[j], tmp);      \
      }                                                                       \
      /* and label them with the appropriate number */                        \
      for (j=num-1; j>=0; j--) {                                              \
         if (tmp[j]) {                                                        \
            data[j] = i;                                                      \
            tmp[j] = 0;                                                       \
         }                                                                    \
      }                                                                       \
   }                                                                          \
   delete [] tmp;                                                             \
   return 1;                                                                  \
}

fragment_data(atomsel_pfrag, pfragList)
fragment_data(atomsel_nfrag, nfragList)

static Timestep *selframe(DrawMolecule *atom_sel_mol, int which_frame) {
  Timestep *r;
  switch (which_frame) {
   case AtomSel::TS_LAST: r = atom_sel_mol->get_last_frame(); break;

   case AtomSel::TS_NOW : r = atom_sel_mol->current(); break;
   default: {
     if (!atom_sel_mol->get_frame(which_frame)) {
       r = atom_sel_mol->get_last_frame();

     } else {
       r = atom_sel_mol->get_frame(which_frame);
     }
   }
  }
  return r;
}

static int atomsel_user(void *v, int num, double *data, int *flgs) {
  DrawMolecule *atom_sel_mol = ((atomsel_ctxt *)v)->atom_sel_mol;                     
  int which_frame = ((atomsel_ctxt *)v)->which_frame;
  Timestep *ts = selframe(atom_sel_mol, which_frame);
  if (!ts || !ts->user) {
    memset(data, 0, num*sizeof(double));
    return 1;                                                           
  }
  for (int i=0; i<num; i++) {                                         
    if (flgs[i]) {                                                            
      data[i] = ts->user[i];                                          
    }                                                                         
  }                                                                           
  return 1;                                                                  
}
static int atomsel_set_user(void *v, int num, double *data, int *flgs) {
  DrawMolecule *atom_sel_mol = ((atomsel_ctxt *)v)->atom_sel_mol;                     
  int which_frame = ((atomsel_ctxt *)v)->which_frame;
  Timestep *ts = selframe(atom_sel_mol, which_frame);
  if (!ts) return 0;
  if (!ts->user) {
    ts->user= new float[num];
    memset(ts->user, 0, num*sizeof(float));
  }
  for (int i=0; i<num; i++) {                                         
    if (flgs[i]) {                                                            
      ts->user[i] = (float)data[i];                                           
    }                                                                         
  }                                                                           
  return 1;                                                                  
}

static int atomsel_user2(void *v, int num, double *data, int *flgs) {
  DrawMolecule *atom_sel_mol = ((atomsel_ctxt *)v)->atom_sel_mol;                     
  int which_frame = ((atomsel_ctxt *)v)->which_frame;
  Timestep *ts = selframe(atom_sel_mol, which_frame);
  if (!ts || !ts->user2) {
    memset(data, 0, num*sizeof(double));
    return 1;                                                           
  }
  for (int i=0; i<num; i++) {                                         
    if (flgs[i]) {                                                            
      data[i] = ts->user2[i];                                         
    }                                                                         
  }                                                                           
  return 1;                                                                  
}
static int atomsel_set_user2(void *v, int num, double *data, int *flgs) {
  DrawMolecule *atom_sel_mol = ((atomsel_ctxt *)v)->atom_sel_mol;                     
  int which_frame = ((atomsel_ctxt *)v)->which_frame;
  Timestep *ts = selframe(atom_sel_mol, which_frame);
  if (!ts) return 0;
  if (!ts->user2) {
    ts->user2= new float[num];
    memset(ts->user2, 0, num*sizeof(float));
  }
  for (int i=0; i<num; i++) {                                         
    if (flgs[i]) {                                                            
      ts->user2[i] = (float)data[i];                                          
    }                                                                         
  }                                                                           
  return 1;                                                                  
}

static int atomsel_user3(void *v, int num, double *data, int *flgs) {
  DrawMolecule *atom_sel_mol = ((atomsel_ctxt *)v)->atom_sel_mol;                     
  int which_frame = ((atomsel_ctxt *)v)->which_frame;
  Timestep *ts = selframe(atom_sel_mol, which_frame);
  if (!ts || !ts->user3) {
    memset(data, 0, num*sizeof(double));
    return 1;                                                           
  }
  for (int i=0; i<num; i++) {                                         
    if (flgs[i]) {                                                            
      data[i] = ts->user3[i];                                         
    }                                                                         
  }                                                                           
  return 1;                                                                  
}
static int atomsel_set_user3(void *v, int num, double *data, int *flgs) {
  DrawMolecule *atom_sel_mol = ((atomsel_ctxt *)v)->atom_sel_mol;                     
  int which_frame = ((atomsel_ctxt *)v)->which_frame;
  Timestep *ts = selframe(atom_sel_mol, which_frame);
  if (!ts) return 0;
  if (!ts->user3) {
    ts->user3= new float[num];
    memset(ts->user3, 0, num*sizeof(float));
  }
  for (int i=0; i<num; i++) {                                         
    if (flgs[i]) {                                                            
      ts->user3[i] = (float)data[i];                                          
    }                                                                         
  }                                                                           
  return 1;                                                                  
}

static int atomsel_user4(void *v, int num, double *data, int *flgs) {
  DrawMolecule *atom_sel_mol = ((atomsel_ctxt *)v)->atom_sel_mol;                     
  int which_frame = ((atomsel_ctxt *)v)->which_frame;
  Timestep *ts = selframe(atom_sel_mol, which_frame);
  if (!ts || !ts->user4) {
    memset(data, 0, num*sizeof(double));
    return 1;                                                           
  }
  for (int i=0; i<num; i++) {                                         
    if (flgs[i]) {                                                            
      data[i] = ts->user4[i];                                         
    }                                                                         
  }                                                                           
  return 1;                                                                  
}
static int atomsel_set_user4(void *v, int num, double *data, int *flgs) {
  DrawMolecule *atom_sel_mol = ((atomsel_ctxt *)v)->atom_sel_mol;                     
  int which_frame = ((atomsel_ctxt *)v)->which_frame;
  Timestep *ts = selframe(atom_sel_mol, which_frame);
  if (!ts) return 0;
  if (!ts->user4) {
    ts->user4= new float[num];
    memset(ts->user4, 0, num*sizeof(float));
  }
  for (int i=0; i<num; i++) {                                         
    if (flgs[i]) {                                                            
      ts->user4[i] = (float)data[i];                                          
    }                                                                         
  }                                                                           
  return 1;                                                                  
}



static int atomsel_xpos(void *v, int num, double *data, int *flgs) {
  DrawMolecule *atom_sel_mol = ((atomsel_ctxt *)v)->atom_sel_mol;                     
  int which_frame = ((atomsel_ctxt *)v)->which_frame;
  Timestep *ts = selframe(atom_sel_mol, which_frame);
  int i;
  if (!ts) {
    for (i=0; i<num; i++) 
      if (flgs[i]) data[i] = 0.0;
    return 0;                                                           
  }
  for (i=0; i<num; i++) {                                             
    if (flgs[i]) {                                                            
      data[i] = ts->pos[3*i + 0];                                             
    }                                                                         
  }                                                                           
  return 1;                                                                  
}

static int atomsel_ypos(void *v, int num, double *data, int *flgs) {
  DrawMolecule *atom_sel_mol = ((atomsel_ctxt *)v)->atom_sel_mol;                     
  int which_frame = ((atomsel_ctxt *)v)->which_frame;
  Timestep *ts = selframe(atom_sel_mol, which_frame);
  int i;
  if (!ts) {
    for (i=0; i<num; i++) 
      if (flgs[i]) data[i] = 0.0;
    return 0;                                                           
  }
  for (i=0; i<num; i++) {                                             
    if (flgs[i]) {                                                            
      data[i] = ts->pos[3*i + 1];                                             
    }                                                                         
  }                                                                           
  return 1;                                                                  
}

static int atomsel_zpos(void *v, int num, double *data, int *flgs) {
  DrawMolecule *atom_sel_mol = ((atomsel_ctxt *)v)->atom_sel_mol;                     
  int which_frame = ((atomsel_ctxt *)v)->which_frame;
  Timestep *ts = selframe(atom_sel_mol, which_frame);
  int i;
  if (!ts) {
    for (i=0; i<num; i++) 
      if (flgs[i]) data[i] = 0.0;
    return 0;                                                           
  }
  for (i=0; i<num; i++) {                                             
    if (flgs[i]) {                                                            
      data[i] = ts->pos[3*i + 2];                                             
    }                                                                         
  }                                                                           
  return 1;                                                                  
}

static int atomsel_ufx(void *v, int num, double *data, int *flgs) {
  DrawMolecule *atom_sel_mol = ((atomsel_ctxt *)v)->atom_sel_mol;                     
  int which_frame = ((atomsel_ctxt *)v)->which_frame;
  Timestep *ts = selframe(atom_sel_mol, which_frame);
  int i;
  if (!ts || !ts->force) {
    for (i=0; i<num; i++) 
      if (flgs[i]) data[i] = 0.0;
    return 0;                                                           
  }
  for (i=0; i<num; i++) {                                             
    if (flgs[i]) {                                                            
      data[i] = ts->force[3*i + 0];                                           
    }                                                                         
  }                                                                           
  return 1;                                                                  
}

static int atomsel_ufy(void *v, int num, double *data, int *flgs) {
  DrawMolecule *atom_sel_mol = ((atomsel_ctxt *)v)->atom_sel_mol;                     
  int which_frame = ((atomsel_ctxt *)v)->which_frame;
  Timestep *ts = selframe(atom_sel_mol, which_frame);
  int i;
  if (!ts || !ts->force) {
    for (i=0; i<num; i++) 
      if (flgs[i]) data[i] = 0.0;
    return 0;                                                           
  }
  for (i=0; i<num; i++) {                                             
    if (flgs[i]) {                                                            
      data[i] = ts->force[3*i + 1];                                           
    }                                                                         
  }                                                                           
  return 1;                                                                  
}

static int atomsel_ufz(void *v, int num, double *data, int *flgs) {
  DrawMolecule *atom_sel_mol = ((atomsel_ctxt *)v)->atom_sel_mol;                     
  int which_frame = ((atomsel_ctxt *)v)->which_frame;
  Timestep *ts = selframe(atom_sel_mol, which_frame);
  int i;
  if (!ts || !ts->force) {
    for (i=0; i<num; i++) 
      if (flgs[i]) data[i] = 0.0;
    return 0;                                                           
  }
  for (i=0; i<num; i++) {                                             
    if (flgs[i]) {                                                            
      data[i] = ts->force[3*i + 2];                                           
    }                                                                         
  }                                                                           
  return 1;                                                                  
}

#define atomsel_get_vel(name, offset) \
static int atomsel_##name(void *v, int num, double *data, int *flgs) {  \
  DrawMolecule *atom_sel_mol = ((atomsel_ctxt *)v)->atom_sel_mol;       \
  int which_frame = ((atomsel_ctxt *)v)->which_frame;  \
  Timestep *ts = selframe(atom_sel_mol, which_frame);  \
  int i;                                               \
  if (!ts || !ts->vel) {                               \
    for (i=0; i<num; i++)                              \
      if (flgs[i]) data[i] = 0.0;                      \
    return 0;                                          \
  }                                                    \
  for (i=0; i<num; i++) {                              \
    if (flgs[i]) {                                     \
      data[i] = ts->vel[3*i + (offset)];               \
    }                                                  \
  }                                                    \
  return 1;                                            \
}                                                      

atomsel_get_vel(vx, 0)
atomsel_get_vel(vy, 1)
atomsel_get_vel(vz, 2)

static int atomsel_phi(void *v, int num, double *data, int *flgs) {
  DrawMolecule *atom_sel_mol = ((atomsel_ctxt *)v)->atom_sel_mol;                     
  int which_frame = ((atomsel_ctxt *)v)->which_frame;
  Timestep *ts =