main.C

Go to the documentation of this file.

/* molecule, Copyright (c) 2001 Jamie Zawinski <jwz@jwz.org>
 * Draws molecules, based on coordinates from PDB (Protein Data Base) files.
 *
 * Permission to use, copy, modify, distribute, and sell this software and its
 * documentation for any purpose is hereby granted without fee, provided that
 * the above copyright notice appear in all copies and that both that
 * copyright notice and this permission notice appear in supporting
 * documentation.  No representations are made about the suitability of this
 * software for any purpose.  It is provided "as is" without express or 
 * implied warranty.
 */

#define countof(x) (sizeof((x))/sizeof((*x)))

#include <string.h>
#include <stdlib.h>
#include <stdio.h>
#include <ctype.h>
#include <time.h>
#include <sys/time.h>

#include <cfloat>

#include <archon/util/ref.H>
#include <archon/util/time.H>
#include <archon/util/thread.H>
#include <archon/util/color.H>
#include <archon/util/options.H>
#include <archon/util/text.H>

#include <archon/x3d/proxy/exception.H>
#include <archon/x3d/proxy/proxy.H>

using namespace Archon::Math;
using namespace Archon::Utilities;
using namespace Archon::X3D::Proxy;

namespace Archon
{
  const char *progname;

  struct AtomKind
  {
    string name;
    float size;
    string textureFile;
    // Vector3 color;
    // Vector3 text_color;
    Ref<ChildNode> model; // Filled in when needed first time
  };



  /* These are the traditional colors used to render these atoms,
     and their approximate size in angstroms.
  */
  static AtomKind allAtomKinds[] =
  {
    { "H", 1.17,  "h.png" /* Utilities::Color::white,           Utilities::Color::gray          */ },
    { "C", 1.75,  "c.png" /* Utilities::Color::gray,            Utilities::Color::white         */ },
    { "N", 1.55,  "n.png" /* Utilities::Color::lightsteelblue,  Utilities::Color::slateblue     */ },
    { "O", 1.40,  "o.png" /* Utilities::Color::red,             Utilities::Color::lightpink     */ },
    { "P", 1.28,  "p.png" /* Utilities::Color::mediumpurple,    Utilities::Color::palevioletred */ },
    { "S", 1.80,  "s.png" /* Utilities::Color::yellow,          Utilities::Color::goldenrod     */ },
    { "*", 1.40,  "x.png" /* Utilities::Color::green,           Utilities::Color::lightgreen    */ }
  };


  typedef struct
  {
    int id;             /* sequence number in the PDB file */
    const char *label;  /* The atom name */
    float x, y, z;      /* position in 3-space (angstroms) */
    AtomKind *kind;     /* computed: which style of atom this is */
  }
  molecule_atom;

  typedef struct
  {
    int from, to;               /* atom sequence numbers */
    int strength;               /* how many bonds are between these two atoms */
  }
  molecule_bond;

  typedef struct
  {
    const char *label;          /* description of this compound */
    int natoms, atoms_size;
    int nbonds, bonds_size;
    molecule_atom *atoms;
    molecule_bond *bonds;
  }
  molecule;

  static molecule *mol = 0;

  Ref<Application> application;
  Ref<Group> rootGroup;

  Rotation3 rotate_cylinder(Vector3 p1, Vector3 p2)
  {
    double p, angle;
    Vector3 v = p1;
    v -= p2;
    v.normalize();
    p = v[1];
    if(p < -1 + DBL_EPSILON*10) return Rotation3::zero(); 
    // Parallel same direction as negative y-axis
    if(p >  1 - DBL_EPSILON*10)         // Parallel opposite direction as negative y-axis
    {
      v = Vector3(1,0,0);
      angle = M_PI;
    }
    else
    {
      v *= Vector3(0,1,0);
      v.normalize();
      angle = acos(-p);
    }
    return Rotation3(v, angle);
  }

  void tube(double fx, double fy, double fz,
            double tx, double ty, double tz,
            double thickness, Vector3 &color)
  {
    Vector3 a(fx, fy, fz);
    Vector3 b(tx, ty, tz);

    Ref<Transform> t_tube = Transform::create(application);
    t_tube->rotation = rotate_cylinder(a, b);
    t_tube->translation = (a+b)/2;

    Ref<Cylinder> cyl = Cylinder::create(application);
    cyl->radius = thickness;
    cyl->height = dist(a,b);

    Ref<Material> mat = Material::create(application);
    mat->diffuseColor = color;

    Ref<Appearance> app = Appearance::create(application);
    app->material = mat;

    Ref<Shape> s = Shape::create(application);
    s->geometry = cyl;
    s->appearance = app;

    t_tube->children.add(s);
    rootGroup->children.add(t_tube);
  }


static AtomKind *get_atom_data(const char *atom_name)
{
  AtomKind *d = 0;
  char *n = strdup(atom_name);
  char *n2 = n;
  int L;

  while(!isalpha(*n)) n++;
  L = strlen(n);
  while(L > 0 && !isalpha(n[L-1])) n[--L] = 0;

  for(unsigned i = 0; i < countof(allAtomKinds); i++)
  {
    d = &allAtomKinds[i];
    if(!strcmp(n, allAtomKinds[i].name.c_str())) break;
  }

  free(n2);
  return d;
}


struct ParseException: Utilities::Exception
{
  ParseException(string l, string m): Utilities::Exception(l, m) {}
};

static molecule_atom *get_atom(molecule_atom *atoms, int natoms, int id)
{
  for(int i = 0; i < natoms; i++) if(id == atoms[i].id) return &atoms[i];

  ARCHON_THROW1(ParseException, "Undefined atom id '" + Utilities::Text::toString(id) + "'");
}

static Ref<ChildNode> buildAtomModel(molecule_atom *atom)
{
  /*
   * let the molecules have the same relative sizes, but scale
   * them to a smaller range, so that the bond-tubes are
   * actually visible...
   */
  double bot = 0.4;
  double top = 0.6;
  double min = 1.17;
  double max = 1.80;
  double ratio = (atom->kind->size - min) / (max - min);
  double size = bot + (ratio * (top - bot));

  Ref<Sphere> sphere = Sphere::create(application);
  sphere->radius = size;

  Ref<Material> material = Material::create(application);
  material->diffuseColor = Vector3(1, 1, 1);

  Ref<ImageTexture> texture = ImageTexture::create(application);
  //texture->repeatS = false;
  //texture->repeatT = false;
  texture->url.add(atom->kind->textureFile);

  Ref<TextureTransform> textureTransform = TextureTransform::create(application);
  Vector2 scale(4, 2);
  scale /= 4.0/3;
  textureTransform->scale = scale;
  textureTransform->translation = Vector2(0.5/scale[0]-0.5, 0.5/scale[1]-0.5);

  Ref<Appearance> appearance = Appearance::create(application);
  appearance->material = material;
  appearance->texture = texture;
  appearance->textureTransform = textureTransform;

  Ref<Shape> shape = Shape::create(application);
  shape->geometry = sphere;
  shape->appearance = appearance;

  Ref<Billboard> billboard = Billboard::create(application);
  billboard->axisOfRotation = Vector3::zero();
  billboard->children.add(shape);

  return billboard;
}

// Construct the graphical model of the molecule
static void build_molecule()
{
  molecule *m = mol;

  for(int i = 0; i < m->nbonds; i++)
    {
      Vector3 color = Utilities::Color::gray;

      molecule_bond *b = &m->bonds[i];
      molecule_atom *from = get_atom(m->atoms, m->natoms, b->from);
      molecule_atom *to   = get_atom(m->atoms, m->natoms, b->to);
      
      float thickness = 0.07 * b->strength;
      if(thickness > 0.3) thickness = 0.3;

      tube(from->x, from->y, from->z, to->x, to->y, to->z, thickness, color);
    }

  for(int i = 0; i < m->natoms; i++)
    {
      molecule_atom *atom = &m->atoms[i];

      if(!atom->kind->model) atom->kind->model = buildAtomModel(atom);

      Ref<Transform> t = Transform::create(application);
      t->translation = Vector3(atom->x, atom->y, atom->z);
      t->children.add(atom->kind->model);

      rootGroup->children.add(t);
    }
}

/* loading */

static void push_atom(molecule *m, int id, const char *label, float x, float y, float z)
{
  m->natoms++;
  if(m->atoms_size < m->natoms)
    {
      m->atoms_size += 20;
      m->atoms = (molecule_atom *) realloc (m->atoms,
                                            m->atoms_size * sizeof(*m->atoms));
    }
  m->atoms[m->natoms-1].id = id;
  m->atoms[m->natoms-1].label = label;
  m->atoms[m->natoms-1].x = x;
  m->atoms[m->natoms-1].y = y;
  m->atoms[m->natoms-1].z = z;
  m->atoms[m->natoms-1].kind = get_atom_data (label);
}


static void push_bond(molecule *m, int from, int to)
{
  int i;

  for(i = 0; i < m->nbonds; i++)
    if ((m->bonds[i].from == from && m->bonds[i].to   == to) ||
        (m->bonds[i].to   == from && m->bonds[i].from == to))
      {
        m->bonds[i].strength++;
        return;
      }

  m->nbonds++;
  if (m->bonds_size < m->nbonds)
    {
      m->bonds_size += 20;
      m->bonds = (molecule_bond *) realloc (m->bonds,
                                            m->bonds_size * sizeof(*m->bonds));
    }
  m->bonds[m->nbonds-1].from = from;
  m->bonds[m->nbonds-1].to = to;
  m->bonds[m->nbonds-1].strength = 1;
}



/* This function is crap.
 */
static void parse_pdb_data(molecule *m, const char *data, string file, int line)
{
  const char *s = data;
  char *ss;
  while (*s)
  {
    if((!m->label || !*m->label) &&
        (!strncmp (s, "HEADER", 6) || !strncmp (s, "COMPND", 6)))
    {
      char *name = new char[100];
      char *n2 = name;
      int L = strlen(s);
      if (L > 99) L = 99;

      strncpy (n2, s, L);
      n2 += 7;
      while (isspace(*n2)) n2++;

      ss = strchr (n2, '\n');
      if (ss) *ss = 0;
      ss = strchr (n2, '\r');
      if (ss) *ss = 0;

      ss = n2+strlen(n2)-1;
      while (isspace(*ss) && ss > n2) *ss-- = 0;

      if (strlen (n2) > 4 && !strcmp (n2 + strlen(n2) - 4, ".pdb"))
        n2[strlen(n2)-4] = 0;

      if (m->label) free ((char *) m->label);
      m->label = strdup (n2);
      free (name);
    }
    else if(!strncmp (s, "TITLE ", 6) ||
            !strncmp (s, "HEADER", 6) ||
            !strncmp (s, "COMPND", 6) ||
            !strncmp (s, "AUTHOR", 6) ||
            !strncmp (s, "REVDAT", 6) ||
            !strncmp (s, "SOURCE", 6) ||
            !strncmp (s, "EXPDTA", 6) ||
            !strncmp (s, "JRNL  ", 6) ||
            !strncmp (s, "REMARK", 6) ||
            !strncmp (s, "SEQRES", 6) ||
            !strncmp (s, "HET   ", 6) ||
            !strncmp (s, "FORMUL", 6) ||
            !strncmp (s, "CRYST1", 6) ||
            !strncmp (s, "ORIGX1", 6) ||
            !strncmp (s, "ORIGX2", 6) ||
            !strncmp (s, "ORIGX3", 6) ||
            !strncmp (s, "SCALE1", 6) ||
            !strncmp (s, "SCALE2", 6) ||
            !strncmp (s, "SCALE3", 6) ||
            !strncmp (s, "MASTER", 6) ||
            !strncmp (s, "KEYWDS", 6) ||
            !strncmp (s, "DBREF ", 6) ||
            !strncmp (s, "HETNAM", 6) ||
            !strncmp (s, "HETSYN", 6) ||
            !strncmp (s, "HELIX ", 6) ||
            !strncmp (s, "LINK  ", 6) ||
            !strncmp (s, "MTRIX1", 6) ||
            !strncmp (s, "MTRIX2", 6) ||
            !strncmp (s, "MTRIX3", 6) ||
            !strncmp (s, "SHEET ", 6) ||
            !strncmp (s, "CISPEP", 6) ||
            !strncmp (s, "GENERATED BY", 12) ||
            !strncmp (s, "TER ", 4) ||
            !strncmp (s, "END ", 4) ||
            !strncmp (s, "TER\n", 4) ||
            !strncmp (s, "END\n", 4) ||
            !strncmp (s, "\n", 1))
      /* ignored. */;
    else if (!strncmp (s, "ATOM   ", 7))
    {
      int id;
      char *name = (char *) calloc (1, 4);
      float x = -999, y = -999, z = -999;

      sscanf (s+7, " %d ", &id);

      strncpy (name, s+12, 3);
      while (isspace(*name)) name++;
      ss = name + strlen(name)-1;
      while (isspace(*ss) && ss > name) *ss-- = 0;
      sscanf (s + 32, " %f %f %f ", &x, &y, &z);

      push_atom (m, id, name, x, y, z);
    }
    else if(!strncmp (s, "HETATM ", 7))
    {
      int id;
      char *name = (char *) calloc (1, 4);
      float x = -999, y = -999, z = -999;

      sscanf (s+7, " %d ", &id);

      strncpy (name, s+12, 3);
      while (isspace(*name)) name++;
      ss = name + strlen(name)-1;
      while (isspace(*ss) && ss > name) *ss-- = 0;
      sscanf (s + 30, " %f %f %f ", &x, &y, &z);

      push_atom (m, id, name, x, y, z);
    }
    else if (!strncmp (s, "CONECT ", 7))
    {
      int atoms[11];
      int i = sscanf(s + 8, " %d %d %d %d %d %d %d %d %d %d %d %d ",
                     &atoms[0], &atoms[1], &atoms[2], &atoms[3], 
                     &atoms[4], &atoms[5], &atoms[6], &atoms[7], 
                     &atoms[8], &atoms[9], &atoms[10], &atoms[11]);
      int j;
      for(j = 1; j < i; j++)
        if(atoms[j] > 0) push_bond (m, atoms[0], atoms[j]);
    }
    else
    {
      char *s1 = strdup (s);
      for (ss = s1; *ss && *ss != '\n'; ss++)  /* WHAT??? */;
      *ss = 0;
      fprintf(stderr, "%s: %s: %d: unrecognised line: %s\n",
              progname, file.c_str(), line, s1);
    }

    while (*s && *s != '\n') s++;
    if(*s == '\n') s++;
    line++;
  }
}


static void parse_pdb_file(molecule *m, string file)
{
  FILE *in;
  int buf_size = 40960;
  char *buf;
  int line = 1;

  in = fopen(file.c_str(), "r");
  if(!in)
  {
    char *buf = (char *) malloc(1024 + file.size());
    sprintf(buf, "%s: error reading \"%s\"", progname, file.c_str());
    perror(buf);
    exit (1);
  }

  buf = (char *) malloc (buf_size);

  while (fgets (buf, buf_size-1, in))
  {
    char *s;
    for (s = buf; *s; s++)
      if (*s == '\r') *s = '\n';
    parse_pdb_data (m, buf, file, line++);
  }

  free (buf);
  fclose (in);

  if(!m->natoms)
  {
    fprintf(stderr, "%s: file %s contains no atomic coordinates!\n",
            progname, file.c_str());
    exit (1);
  }

  if(!m->nbonds)
  {
    fprintf(stderr, "%s: warning: file %s contains no atomic bond info.\n",
            progname, file.c_str());
  }
}


typedef struct { char *atom; int count; } atom_and_count;

/* When listing the components of a molecule, the convention is to put the
   carbon atoms first, the hydrogen atoms second, and the other atom types
   sorted alphabetically after that (although for some molecules, the usual
   order is different: we special-case a few of those.)
 */
static int cmp_atoms(const void *aa, const void *bb)
{
  const atom_and_count *a = (atom_and_count *) aa;
  const atom_and_count *b = (atom_and_count *) bb;
  if (!a->atom) return  1;
  if (!b->atom) return -1;
  if (!strcmp(a->atom, "C")) return -1;
  if (!strcmp(b->atom, "C")) return  1;
  if (!strcmp(a->atom, "H")) return -1;
  if (!strcmp(b->atom, "H")) return  1;
  return strcmp (a->atom, b->atom);
}

/* thanks to Rene Uittenbogaard <ruittenb@wish.nl> */
static void special_case_formula (char *f)
{
  if      (!strcmp(f, "H(2)Be"))   strcpy(f, "BeH(2)");
  else if (!strcmp(f, "H(3)B"))    strcpy(f, "BH(3)");
  else if (!strcmp(f, "H(3)N"))    strcpy(f, "NH(3)");
  else if (!strcmp(f, "CHN"))      strcpy(f, "HCN");
  else if (!strcmp(f, "CKN"))      strcpy(f, "KCN");
  else if (!strcmp(f, "H(4)N(2)")) strcpy(f, "N(2)H(4)");
  else if (!strcmp(f, "Cl(3)P"))   strcpy(f, "PCl(3)");
  else if (!strcmp(f, "Cl(5)P"))   strcpy(f, "PCl(5)");
}

static void generate_molecule_formula(molecule *m)
{
  char *buf = (char *) malloc (m->natoms * 10);
  char *s = buf;
  int i;
  atom_and_count counts[200];
  memset (counts, 0, sizeof(counts));
  *s = 0;
  for(i = 0; i < m->natoms; i++)
  {
    int j = 0;
    char *a = (char *) m->atoms[i].label;
    char *e;
    while(!isalpha(*a)) a++;
    a = strdup (a);
    for(e = a; isalpha(*e); e++);
    *e = 0;
    while (counts[j].atom && !!strcmp(a, counts[j].atom)) j++;
    if (counts[j].atom) free (a);
    else counts[j].atom = a;
    counts[j].count++;
  }

  i = 0;
  while (counts[i].atom) i++;
  qsort (counts, i, sizeof(*counts), cmp_atoms);

  i = 0;
  while(counts[i].atom)
  {
    strcat (s, counts[i].atom);
    free (counts[i].atom);
    s += strlen (s);
    if(counts[i].count > 1) sprintf (s, "(%d)", counts[i].count);
    s += strlen (s);
    i++;
  }

  special_case_formula (buf);

  if (!m->label) m->label = strdup("");
  s = (char *) malloc (strlen (m->label) + strlen (buf) + 2);
  strcpy (s, m->label);
  strcat (s, "\n");
  strcat (s, buf);
  free ((char *) m->label);
  free (buf);
  m->label = s;
}



static void insert_vertical_whitespace (char *string)
{
  while (*string)
  {
    if((string[0] == ',' || string[0] == ';' || string[0] == ':') && string[1] == ' ')
      string[0] = ' ', string[1] = '\n';
    string++;
  }
}


void init_molecule(string file)
{
  mol = (molecule *)calloc(sizeof (molecule), 1);
  parse_pdb_file(mol, file);
  generate_molecule_formula(mol);
  insert_vertical_whitespace((char *)mol->label);
}




  namespace SaiTestApps
  {
    namespace Molecule
    {
      static bool quit = false;

      int main(int argc, const char *argv[])
      {
        bool   opt_help              = false;
        int    opt_corbaEndpointPort = -1; // Auto assign
        string opt_serverName        = "archon";
        
        Options o;

        o.addSwitch("h", "help", opt_help, true,
                    "Describe the parameters");
        o.addConfig("p", "corba-endpoint-port", opt_corbaEndpointPort, opt_corbaEndpointPort,
                    "The TCP/IP port on which the CORBA Orb listens (-1 for auto assign)",
                    Options::wantArg_always);
        o.addConfig("n", "server-name", opt_serverName, opt_serverName,
                    "The name of the server to connect to",
                    Options::wantArg_always);
        
        if(o.processCommandLine(argc, argv))
        {
          cerr << "Try --help\n";
          return 1;
        }

        if(opt_help)
        {
          cout <<
            "Test Application for 3D-Console by Brian Kristiansen & Kristian Spangsege\n"
            "\n"
            "Synopsis: " << argv[0] << " FILE\n"
            "\n"
            "Available options:\n";
          cout << o.list();
          return 0;
        }

        if(argc < 2 || argc > 2)
        {
          cerr << "Wrong number of aguments\n";
          cerr << "Try --help\n";
          return 1;
        }

        Ref<Session> session = Session::create("archon", -1);
  
        application = Application::create(session, Uri());
        rootGroup = Group::create(application);
  
        progname = argv[0];

        init_molecule(argv[1]);
        build_molecule();

        application->setRootGroup(rootGroup);
        application->launch();

        Archon::Utilities::Time delay(0, 100000000l);
        while(!quit)
        {
          session->beginUpdate();
          session->endUpdate();

          Archon::Utilities::Thread::sleep(delay);
        }

        return 0;
      }
    }
  }
}

int main(int argc, const char *argv[]) throw()
{
  using namespace std;
  using namespace Archon::X3D;

  set_unexpected(Exception::terminal<Proxy::exceptionCatchInfo>);
  set_terminate (Exception::terminal<Proxy::exceptionCatchInfo>);

  return Archon::SaiTestApps::Molecule::main(argc, argv);
}

---