Logo Search packages:      
Sourcecode: pymol version File versions  Download package

ObjectMolecule2.c

/* 
   A* -------------------------------------------------------------------
   B* This file contains source code for the PyMOL computer program
   C* copyright 1998-2000 by Warren Lyford Delano of DeLano Scientific. 
   D* -------------------------------------------------------------------
   E* It is unlawful to modify or remove this copyright notice.
   F* -------------------------------------------------------------------
   G* Please see the accompanying LICENSE file for further information. 
   H* -------------------------------------------------------------------
   I* Additional authors of this source file include:
   -* 
   -* 
   -*
   Z* -------------------------------------------------------------------
*/

#include"os_predef.h"
#include"os_std.h"
#include"os_gl.h"

#include"Base.h"
#include"Debug.h"
#include"Parse.h"
#include"OOMac.h"
#include"Vector.h"
#include"PConv.h"
#include"ObjectMolecule.h"
#include"Feedback.h"
#include"Util.h"
#include"AtomInfo.h"
#include"Selector.h"
#include"ObjectDist.h"
#include"Executive.h"
#include"P.h"
#include"ObjectCGO.h"
#include"Scene.h"

#define ntrim ParseNTrim
#define nextline ParseNextLine
#define ncopy ParseNCopy
#define nskip ParseNSkip

#define cResvMask 0x7FFF

#define cMaxOther 6

00048 typedef struct {
  int n_cyclic_arom, cyclic_arom[cMaxOther];
  int n_arom, arom[cMaxOther];
  int n_high_val, high_val[cMaxOther];
  int n_cyclic, cyclic[cMaxOther];
  int n_planer, planer[cMaxOther];
  int n_rest, rest[cMaxOther];
  int score;
} OtherRec;

static int populate_other(OtherRec * other, int at, AtomInfoType * ai, BondType * bd,
                          int *neighbor)
{
  int five_cycle = false;
  int six_cycle = false;

  {
    int mem[9], nbr[7];
    const int ESCAPE_MAX = 500;
    register int escape_count;

    escape_count = ESCAPE_MAX;  /* don't get bogged down with structures 
                                   that have unreasonable connectivity */
    mem[0] = bd->index[0];
    mem[1] = bd->index[1];
    nbr[1] = neighbor[mem[1]] + 1;
    while(((mem[2] = neighbor[nbr[1]]) >= 0)) {
      if(mem[2] != mem[0]) {
        nbr[2] = neighbor[mem[2]] + 1;
        while(((mem[3] = neighbor[nbr[2]]) >= 0)) {
          if(mem[3] != mem[1]) {
            nbr[3] = neighbor[mem[3]] + 1;
            while(((mem[4] = neighbor[nbr[3]]) >= 0)) {
              if((mem[4] != mem[2]) && (mem[4] != mem[1]) && (mem[4] != mem[0])) {
                nbr[4] = neighbor[mem[4]] + 1;
                while(((mem[5] = neighbor[nbr[4]]) >= 0)) {
                  if(!(escape_count--))
                    goto escape;
                  if((mem[5] != mem[3]) && (mem[5] != mem[2]) && (mem[5] != mem[1])) {
                    if(mem[5] == mem[0]) {      /* five-cycle */
                      five_cycle = true;
                    }
                    nbr[5] = neighbor[mem[5]] + 1;
                    while(((mem[6] = neighbor[nbr[5]]) >= 0)) {
                      if((mem[6] != mem[4]) && (mem[6] != mem[3]) && (mem[6] != mem[2])
                         && (mem[6] != mem[1])) {
                        if(mem[6] == mem[0]) {  /* six-cycle */
                          six_cycle = true;
                        }
                      }
                      nbr[5] += 2;
                    }
                  }
                  nbr[4] += 2;
                }
              }
              nbr[3] += 2;
            }
          }
          nbr[2] += 2;
        }
      }
      nbr[1] += 2;
    }
  }
escape:

  if(bd->order == 4) {          /* aromatic */
    if((five_cycle || six_cycle) && (other->n_cyclic_arom < cMaxOther)) {
      other->cyclic_arom[other->n_cyclic_arom++] = at;
      if(five_cycle && six_cycle)
        other->score += 34;
      else if(five_cycle)
        other->score += 33;
      else
        other->score += 32;
      return 1;
    } else if(other->n_arom < cMaxOther) {
      other->arom[other->n_arom++] = at;
      other->score += 64;
      return 1;
    }
  }
  if(bd->order > 1) {
    if(other->n_high_val < cMaxOther) {
      other->high_val[other->n_high_val++] = at;
      other->score += 16;
      return 1;
    }
  }
  if(five_cycle || six_cycle) {
    if(other->n_cyclic < cMaxOther) {
      other->cyclic[other->n_cyclic++] = at;
      other->score += 8;
      return 1;
    }
  }
  if(ai->geom == cAtomInfoPlanar) {
    if(other->n_planer < cMaxOther) {
      other->planer[other->n_planer++] = at;
      other->score += 4;
      return 1;
    }
  }
  if(other->n_rest < cMaxOther) {
    other->rest[other->n_rest++] = at;
    other->score += 1;
    return 1;
  }
  return 0;
}

static int append_index(int *result, int offset, int a1, int a2, int score, int ar_count)
{
  int c;
  c = result[a1];
  while(c < offset) {
    if(result[c] == a2) {       /* already entered */
      if(result[c + 1] < score) {
        result[c + 1] = score;
        result[c + 2] = ar_count;
      }
      return offset;
    }
    c += 3;
  }
  result[offset++] = a2;
  result[offset++] = score;
  result[offset++] = ar_count;
  return offset;
}

int ObjectMoleculeAddPseudoatom(ObjectMolecule * I, int sele_index, char *name,
                                char *resn, char *resi, char *chain,
                                char *segi, char *elem, float vdw,
                                int hetatm, float b, float q, char *label,
                                float *pos, int color, int state, int mode, int quiet)
{
  PyMOLGlobals *G = I->Obj.G;
  int start_state = 0, stop_state = 0;
  int nAtom = 1;
  int extant_only = false;
  int ai_merged = false;
  float pos_array[3] = { 0.0F, 0.0F, 0.0F };
  int ok = true;

  AtomInfoType *atInfo = VLACalloc(AtomInfoType, 1);

  if(state >= 0) {              /* specific state */
    start_state = state;
    stop_state = state + 1;
  } else if(state == -1) {      /* current state */
    start_state = ObjectGetCurrentState(&I->Obj, true);
    stop_state = start_state + 1;
  } else {                      /* all states */
    if(sele_index >= 0) {
      start_state = 0;
      stop_state = SelectorCountStates(G, sele_index);
      if(state == -3)
        extant_only = true;
    } else {
      start_state = 0;
      stop_state = ExecutiveCountStates(G, cKeywordAll);
      if(stop_state < 1)
        stop_state = 1;
    }
  }
  {
    /* match existing properties of the old atom */
    int auto_show_lines = (int) SettingGet(G, cSetting_auto_show_lines);
    int auto_show_spheres = (int) SettingGet(G, cSetting_auto_show_spheres);
    int auto_show_nonbonded = (int) SettingGet(G, cSetting_auto_show_nonbonded);
    AtomInfoType *ai = atInfo;
    ai->resv = AtomResvFromResi(resi);
    ai->hetatm = hetatm;
    ai->geom = cAtomInfoNone;
    ai->q = q;
    ai->b = b;
    strcpy(ai->chain, chain);
    strcpy(ai->resi, resi);
    strcpy(ai->segi, segi);
    strcpy(ai->resn, resn);
    strcpy(ai->name, name);
    strcpy(ai->elem, elem);
    ai->visRep[cRepLine] = auto_show_lines;
    ai->visRep[cRepNonbonded] = auto_show_nonbonded;
    ai->visRep[cRepSphere] = auto_show_spheres;
    ai->id = -1;
    ai->rank = -1;
    if(vdw >= 0.0F)
      ai->vdw = vdw;
    else
      ai->vdw = 1.0F;
    if(label[0]) {
      OVreturn_word ret = OVLexicon_GetFromCString(G->Lexicon, label);
      if(OVreturn_IS_OK(ret)) {
        ai->label = ret.word;
        ai->visRep[cRepLabel] = true;
        ai->visRep[cRepLine] = false;
        ai->visRep[cRepNonbonded] = false;
        ai->visRep[cRepSphere] = false;
      }
    }
    if(color < 0) {
      AtomInfoAssignColors(I->Obj.G, ai);
      if((ai->elem[0] == 'C') && (ai->elem[1] == 0))
        /* carbons are always colored according to the object color */
        ai->color = I->Obj.Color;
    } else {
      ai->color = color;
    }
    AtomInfoAssignParameters(I->Obj.G, ai);
    AtomInfoUniquefyNames(I->Obj.G, I->AtomInfo, I->NAtom, ai, NULL, 1);
    if(!quiet) {
      PRINTFB(G, FB_ObjectMolecule, FB_Actions)
        " ObjMol: created %s/%s/%s/%s`%s/%s\n",
        I->Obj.Name, ai->segi, ai->chain, ai->resn, ai->resi, ai->name ENDFB(G);
    }
  }

  for(state = start_state; state < stop_state; state++) {

    CoordSet *cset = NULL;

    if((extant_only && (state < I->NCSet) && I->CSet[state]) || !extant_only) {

      if(sele_index >= 0) {
        ObjectMoleculeOpRec op;
        ObjectMoleculeOpRecInit(&op);
        op.code = OMOP_CSetSumVertices;
        op.cs1 = state;

        ExecutiveObjMolSeleOp(I->Obj.G, sele_index, &op);

        if(op.i1) {
          float factor = 1.0F / op.i1;
          scale3f(op.v1, factor, pos_array);
          pos = pos_array;

          if(vdw < 0.0F) {
            switch (mode) {
            case 1:
              ObjectMoleculeOpRecInit(&op);
              op.code = OMOP_CSetMaxDistToPt;
              copy3f(pos_array, op.v1);
              op.cs1 = state;
              ExecutiveObjMolSeleOp(I->Obj.G, sele_index, &op);
              vdw = op.f1;
              break;
            case 2:
              ObjectMoleculeOpRecInit(&op);
              op.code = OMOP_CSetSumSqDistToPt;
              copy3f(pos_array, op.v1);
              op.cs1 = state;
              ExecutiveObjMolSeleOp(I->Obj.G, sele_index, &op);
              vdw = sqrt1f(op.d1 / op.i1);
              break;
            case 0:
            default:
              vdw = 0.5F;
              break;
            }
            if(vdw >= 0.0F)
              atInfo->vdw = vdw;        /* NOTE: only uses vdw from first state selection... */
          }
        } else {
          pos = NULL;           /* skip this state */
        }
      } else if(!pos) {
        pos = pos_array;
        SceneGetPos(I->Obj.G, pos);
      }

      if(pos) {                 /* only add coordinate to state if we have position for it */

        float *coord = VLAlloc(float, 3 * nAtom);

        copy3f(pos, coord);

        cset = CoordSetNew(I->Obj.G);
        cset->NIndex = nAtom;
        cset->Coord = coord;
        cset->TmpBond = NULL;
        cset->NTmpBond = 0;

        cset->Obj = I;
        if(cset->fEnumIndices)
          cset->fEnumIndices(cset);
        if(!ai_merged) {
          ObjectMoleculeMerge(I, atInfo, cset, false, cAIC_AllMask, true);      /* NOTE: will release atInfo */
          ObjectMoleculeExtendIndices(I, -1);
          ObjectMoleculeUpdateNeighbors(I);
          ai_merged = true;
        }
        if(state >= I->NCSet) {
          VLACheck(I->CSet, CoordSet *, state);
          I->NCSet = state + 1;
        }
        if(!I->CSet[state]) {
          /* new coordinate set */
          I->CSet[state] = cset;
          cset = NULL;
        } else {
          /* merge coordinate set */
          CoordSetMerge(I->CSet[state], cset);
          if(cset->fFree) {
            cset->fFree(cset);
            cset = NULL;
          }
        }
      }
    }
  }
  if(ai_merged) {
    ObjectMoleculeSort(I);
    ObjectMoleculeUpdateIDNumbers(I);
    ObjectMoleculeUpdateNonbonded(I);
    ObjectMoleculeInvalidate(I, cRepAll, cRepInvAtoms, -1);
  } else {
    VLAFreeP(atInfo);
  }
  return (ok);
}

int *ObjectMoleculeGetPrioritizedOtherIndexList(ObjectMolecule * I, CoordSet * cs)
{
  int a, b;
  int b1, b2, a1, a2, a3;
  OtherRec *o;
  OtherRec *other = Calloc(OtherRec, cs->NIndex);
  int *result = NULL;
  int offset;
  int n_alloc = 0;
  BondType *bd;

  ObjectMoleculeUpdateNeighbors(I);
  bd = I->Bond;
  for(a = 0; a < I->NBond; a++) {
    b1 = bd->index[0];
    b2 = bd->index[1];
    if(I->DiscreteFlag) {
      if((cs == I->DiscreteCSet[b1]) && (cs == I->DiscreteCSet[b2])) {
        a1 = I->DiscreteAtmToIdx[b1];
        a2 = I->DiscreteAtmToIdx[b2];
      } else {
        a1 = -1;
        a2 = -1;
      }
    } else {
      a1 = cs->AtmToIdx[b1];
      a2 = cs->AtmToIdx[b2];
    }
    if((a1 >= 0) && (a2 >= 0)) {
      n_alloc += populate_other(other + a1, a2, I->AtomInfo + b2, bd, I->Neighbor);
      n_alloc += populate_other(other + a2, a1, I->AtomInfo + b1, bd, I->Neighbor);
    }
    bd++;
  }

  n_alloc = 3 * (n_alloc + cs->NIndex);
  o = other;
  result = Alloc(int, n_alloc);
  for(a = 0; a < cs->NIndex; a++) {
    result[a] = -1;
  }
  offset = cs->NIndex;
  bd = I->Bond;
  for(a = 0; a < I->NBond; a++) {
    b1 = bd->index[0];
    b2 = bd->index[1];
    if(I->DiscreteFlag) {
      if((cs == I->DiscreteCSet[b1]) && (cs == I->DiscreteCSet[b2])) {
        a1 = I->DiscreteAtmToIdx[b1];
        a2 = I->DiscreteAtmToIdx[b2];
      } else {
        a1 = -1;
        a2 = -1;
      }
    } else {
      a1 = cs->AtmToIdx[b1];
      a2 = cs->AtmToIdx[b2];
    }
    if((a1 >= 0) && (a2 >= 0)) {
      if(result[a1] < 0) {
        o = other + a1;
        result[a1] = offset;
        for(b = 0; b < o->n_cyclic_arom; b++) {
          a3 = o->cyclic_arom[b];
          offset = append_index(result, offset, a1, a3, 128 + other[a3].score, 1);
        }
        for(b = 0; b < o->n_arom; b++) {
          a3 = o->arom[b];
          offset = append_index(result, offset, a1, a3, 64 + other[a3].score, 1);
        }
        for(b = 0; b < o->n_high_val; b++) {
          a3 = o->high_val[b];
          offset = append_index(result, offset, a1, a3, 16 + other[a3].score, 0);
        }
        for(b = 0; b < o->n_cyclic; b++) {
          a3 = o->cyclic[b];
          offset = append_index(result, offset, a1, a3, 8 + other[a3].score, 0);
        }
        for(b = 0; b < o->n_planer; b++) {
          a3 = o->planer[b];
          offset = append_index(result, offset, a1, a3, 2 + other[a3].score, 0);
        }
        for(b = 0; b < o->n_rest; b++) {
          a3 = o->rest[b];
          offset = append_index(result, offset, a1, a3, 1 + other[a3].score, 0);
        }
        result[offset++] = -1;
      }

      if(result[a2] < 0) {
        o = other + a2;
        result[a2] = offset;
        for(b = 0; b < o->n_cyclic_arom; b++) {
          a3 = o->cyclic_arom[b];
          offset = append_index(result, offset, a2, a3, 128 + other[a3].score, 1);
        }
        for(b = 0; b < o->n_arom; b++) {
          a3 = o->arom[b];
          offset = append_index(result, offset, a2, a3, 64 + other[a3].score, 1);
        }
        for(b = 0; b < o->n_high_val; b++) {
          a3 = o->high_val[b];
          offset = append_index(result, offset, a2, a3, 16 + other[a3].score, 0);
        }
        for(b = 0; b < o->n_cyclic; b++) {
          a3 = o->cyclic[b];
          offset = append_index(result, offset, a2, a3, 8 + other[a3].score, 0);
        }
        for(b = 0; b < o->n_planer; b++) {
          a3 = o->planer[b];
          offset = append_index(result, offset, a2, a3, 2 + other[a3].score, 0);
        }
        for(b = 0; b < o->n_rest; b++) {
          a3 = o->rest[b];
          offset = append_index(result, offset, a2, a3, 1 + other[a3].score, 0);
        }
        result[offset++] = -1;
      }

    }
    bd++;
  }
  FreeP(other);
  return result;
}

int ObjectMoleculeGetNearestBlendedColor(ObjectMolecule * I, float *point,
                                         float cutoff, int state, float *dist,
                                         float *color, int sub_vdw)
{
  int result = -1;
  float tot_weight = 0.0F;
  float cutoff2 = cutoff * cutoff;
  register float nearest = -1.0F;

  color[0] = 0.0F;
  color[1] = 0.0F;
  color[2] = 0.0F;

  if(state < 0)
    state = ObjectGetCurrentState(&I->Obj, true);

  if((state >= 0) && (state < I->NCSet)) {
    CoordSet *cs = I->CSet[state];
    if(cs) {
      MapType *map;
      CoordSetUpdateCoord2IdxMap(cs, cutoff);
      if(sub_vdw) {
        cutoff -= MAX_VDW;
        cutoff2 = cutoff * cutoff;
      }
      nearest = cutoff2;
      if((map = cs->Coord2Idx)) {
        int a, b, c, d, e, f, j;
        register float test;
        register float *v;
        MapLocus(map, point, &a, &b, &c);
        for(d = a - 1; d <= a + 1; d++)
          for(e = b - 1; e <= b + 1; e++)
            for(f = c - 1; f <= c + 1; f++) {
              j = *(MapFirst(map, d, e, f));
              while(j >= 0) {
                v = cs->Coord + (3 * j);
                test = diffsq3f(v, point);
                if(sub_vdw) {
                  test = sqrt1f(test);
                  test -= I->AtomInfo[cs->IdxToAtm[j]].vdw;
                  if(test < 0.0F)
                    test = 0.0F;
                  test = test * test;
                }
                if(test < cutoff2) {
                  float weight = cutoff - sqrt1f(test);
                  float *at_col = ColorGet(I->Obj.G, I->AtomInfo[cs->IdxToAtm[j]].color);
                  color[0] += at_col[0] * weight;
                  color[1] += at_col[1] * weight;
                  color[2] += at_col[2] * weight;
                  tot_weight += weight;
                }
                if(test <= nearest) {
                  result = j;
                  nearest = test;
                }
                j = MapNext(map, j);
              }
            }
      } else {
        register int j;
        register float test, *v = cs->Coord;
        for(j = 0; j < cs->NIndex; j++) {
          test = diffsq3f(v, point);
          if(sub_vdw) {
            test = sqrt1f(test);
            test -= I->AtomInfo[cs->IdxToAtm[j]].vdw;
            if(test < 0.0F)
              test = 0.0F;
            test = test * test;
          }
          if(test < cutoff2) {
            float weight = cutoff - sqrt1f(test);
            float *color = ColorGet(I->Obj.G, I->AtomInfo[cs->IdxToAtm[j]].color);
            color[0] += color[0] * weight;
            color[1] += color[1] * weight;
            color[2] += color[2] * weight;
            tot_weight += weight;
          }
          if(test <= nearest) {
            result = j;
            nearest = test;
          }
          v += 3;
        }
      }
      if(result >= 0)
        result = cs->IdxToAtm[result];
    }
  }
  if(dist) {
    if(result >= 0) {
      *dist = sqrt1f(nearest);
      if(tot_weight > 0.0F) {
        color[0] /= tot_weight;
        color[1] /= tot_weight;
        color[2] /= tot_weight;
      }
    } else {
      *dist = -1.0F;
    }
  }
  return result;
}

int ObjectMoleculeGetNearestAtomIndex(ObjectMolecule * I, float *point, float cutoff,
                                      int state, float *dist)
{
  int result = -1;
  register float nearest = -1.0F;
  if(state < 0)
    state = ObjectGetCurrentState(&I->Obj, true);
  if((state >= 0) && (state < I->NCSet)) {
    CoordSet *cs = I->CSet[state];
    if(cs) {
      MapType *map;
      CoordSetUpdateCoord2IdxMap(cs, cutoff);
      nearest = cutoff * cutoff;
      if((map = cs->Coord2Idx)) {
        int a, b, c, d, e, f, j;
        register float test;
        register float *v;
        MapLocus(map, point, &a, &b, &c);
        for(d = a - 1; d <= a + 1; d++)
          for(e = b - 1; e <= b + 1; e++)
            for(f = c - 1; f <= c + 1; f++) {
              j = *(MapFirst(map, d, e, f));
              while(j >= 0) {
                v = cs->Coord + (3 * j);
                test = diffsq3f(v, point);
                if(test <= nearest) {
                  result = j;
                  nearest = test;
                }
                j = MapNext(map, j);
              }
            }
      } else {
        register int j;
        register float test, *v = cs->Coord;
        for(j = 0; j < cs->NIndex; j++) {
          test = diffsq3f(v, point);
          if(test <= nearest) {
            result = j;
            nearest = test;
          }
          v += 3;
        }
      }
      if(result >= 0)
        result = cs->IdxToAtm[result];
    }
  }
  if(dist) {
    if(result >= 0) {
      *dist = sqrt1f(nearest);
    } else {
      *dist = -1.0F;
    }
  }
  return result;
}

int ObjectMoleculeGetPrioritizedOther(int *other, int a1, int a2, int *double_sided)
{
  int a3 = -1;
  int lvl = -1, ck, ck_lvl;
  int offset;
  int ar_count = 0;

  a3 = -1;
  lvl = -1;
  if(a1 >= 0) {
    offset = other[a1];
    if(offset >= 0) {
      while(1) {
        ck = other[offset];
        if(ck != a2) {
          if(ck >= 0) {
            ck_lvl = other[offset + 1];
            if(ck_lvl > lvl) {
              a3 = ck;
              lvl = ck_lvl;
            }
            ar_count += other[offset + 2];
          } else
            break;
        }
        offset += 3;
      }
    }
  }
  if(a2 >= 0) {
    offset = other[a2];
    if(offset >= 0) {
      while(1) {
        ck = other[offset];
        if(ck != a1) {
          if(ck >= 0) {
            ck_lvl = other[offset + 1];
            if(ck_lvl > lvl) {
              a3 = ck;
              lvl = ck_lvl;
            }
            ar_count += other[offset + 2];
          } else
            break;
        }
        offset += 3;
      }
    }
  }

  if(double_sided) {
    if(ar_count == 4)
      *double_sided = true;
    else
      *double_sided = false;
  }
  return a3;
}

int ObjectMoleculeIsAtomBondedToName(ObjectMolecule * obj, int a0, char *name)
{
  int a2, s;
  int bonded = false;

  if(a0 >= 0) {
    s = obj->Neighbor[a0];
    s++;                        /* skip count */
    while(1) {
      a2 = obj->Neighbor[s];
      if(a2 < 0)
        break;
      if(WordMatch(obj->Obj.G, obj->AtomInfo[a2].name, name, true) < 0)
        bonded = true;
      break;
      s += 2;
    }
  }
  return bonded;
}

int ObjectMoleculeAreAtomsBonded2(ObjectMolecule * obj0, int a0, ObjectMolecule * obj1,
                                  int a1)
{
  /* assumes neighbor list is current */

  if(obj0 != obj1)
    return false;
  else {
    int a2, s;

    if(a0 >= 0) {
      s = obj0->Neighbor[a0];
      s++;                      /* skip count */
      while(1) {
        a2 = obj0->Neighbor[s];
        if(a2 < 0)
          break;
        if(a1 == a2)
          return true;
        s += 2;
      }
    }
  }
  return false;
}

int ObjectMoleculeDoesAtomNeighborSele(ObjectMolecule * I, int index, int sele)
{
  int result = false;
  ObjectMoleculeUpdateNeighbors(I);
  if(index < I->NAtom) {
    int a1;
    int n;
    AtomInfoType *ai;

    n = I->Neighbor[index] + 1;
    while(1) {                  /* look for an attached non-hydrogen as a base */
      a1 = I->Neighbor[n];
      n += 2;
      if(a1 < 0)
        break;
      ai = I->AtomInfo + a1;
      if(SelectorIsMember(I->Obj.G, ai->selEntry, sele)) {
        result = true;
        break;
      }
    }
  }
  return result;
}

static void assign_pdb_known_residue(PyMOLGlobals * G, AtomInfoType * ai1,
                                     AtomInfoType * ai2, int *bond_order)
{
  int order = *(bond_order);
  char *name1 = ai1->name;
  char *name2 = ai2->name;
  char *resn1 = ai1->resn;

  /* nasty high-speed hack to get bond valences and formal charges 
     for standard residues */
  if(((!name1[1]) && (!name2[1])) &&
     (((name1[0] == 'C') && (name2[0] == 'O')) ||
      ((name1[0] == 'O') && (name2[0] == 'C')))) {
    order = 2;
  } else if((!name2[1]) && (name2[0] == 'C') && (!strcmp(name1, "OXT"))) {
    ai1->formalCharge = -1;
    ai1->chemFlag = false;
  } else if((!name1[1]) && (name1[0] == 'C') && (!strcmp(name2, "OXT"))) {
    ai2->formalCharge = -1;
    ai2->chemFlag = false;
  } else {
    switch (resn1[0]) {
    case 'A':
      switch (resn1[1]) {
      case 'R':
        switch (resn1[2]) {
        case 'G':              /* ARG... */
          switch (resn1[3]) {
          case 0:
          case 'P':            /*  ARG, ARGP */
            if(!strcmp(name1, "NH1")) {
              ai1->formalCharge = 1;
              ai1->chemFlag = false;
            } else if(!strcmp(name2, "NH1")) {
              ai2->formalCharge = 1;
              ai2->chemFlag = false;
            }
            break;
          }
          if(((!strcmp(name1, "CZ")) && (!strcmp(name2, "NH1"))) ||
             ((!strcmp(name2, "CZ")) && (!strcmp(name1, "NH1"))))
            order = 2;
          break;
        }
        break;
      case 'S':
        switch (resn1[2]) {
        case 'P':              /* ASP... */
          switch (resn1[3]) {
          case 0:
          case 'M':            /* ASP, ASPM minus assumption */
            if(!strcmp(name1, "OD2")) {
              ai1->formalCharge = -1;
              ai1->chemFlag = false;
            } else if(!strcmp(name2, "OD2")) {
              ai2->formalCharge = -1;
              ai2->chemFlag = false;
            }
            break;
          }
          if(((!strcmp(name1, "CG")) && (!strcmp(name2, "OD1"))) ||
             ((!strcmp(name2, "CG")) && (!strcmp(name1, "OD1"))))
            order = 2;
          break;
        case 'N':              /* ASN  */
          if(((!strcmp(name1, "CG")) && (!strcmp(name2, "OD1"))) ||
             ((!strcmp(name2, "CG")) && (!strcmp(name1, "OD1"))))
            order = 2;
          break;
        }
        break;
      case 0:
        if(((!strcmp(name1, "O2P")) || (!strcmp(name1, "OP2")))) {
          ai1->formalCharge = -1;
          ai1->chemFlag = false;
        } else if(((!strcmp(name2, "O2P")) || (!strcmp(name2, "OP2")))) {
          ai2->formalCharge = -1;
          ai2->chemFlag = false;
        }
        if(((!strcmp(name1, "C8")) && (!strcmp(name2, "N7"))) ||
           ((!strcmp(name2, "C8")) && (!strcmp(name1, "N7"))))
          order = 2;
        else if(((!strcmp(name1, "C4")) && (!strcmp(name2, "C5"))) ||
                ((!strcmp(name2, "C4")) && (!strcmp(name1, "C5"))))
          order = 2;

        else if(((!strcmp(name1, "C6")) && (!strcmp(name2, "N1"))) ||
                ((!strcmp(name2, "C6")) && (!strcmp(name1, "N1"))))
          order = 2;
        else if(((!strcmp(name1, "C2")) && (!strcmp(name2, "N3"))) ||
                ((!strcmp(name2, "C2")) && (!strcmp(name1, "N3"))))
          order = 2;
        else
          if(((!strcmp(name1, "P"))
              && (((!strcmp(name2, "O1P")) || (!strcmp(name2, "OP1")))))
             || ((!strcmp(name2, "P"))
                 && (((!strcmp(name1, "O1P")) || (!strcmp(name1, "OP1"))))))
          order = 2;
        break;
      }
      break;
    case 'C':
      if(resn1[1] == 0) {
        if(((!strcmp(name1, "O2P")) || (!strcmp(name1, "OP2")))) {
          ai1->formalCharge = -1;
          ai1->chemFlag = false;
        } else if(((!strcmp(name2, "O2P")) || (!strcmp(name2, "OP2")))) {
          ai2->formalCharge = -1;
          ai2->chemFlag = false;
        }
        if(((!strcmp(name1, "C2")) && (!strcmp(name2, "O2"))) ||
           ((!strcmp(name2, "C2")) && (!strcmp(name1, "O2"))))
          order = 2;
        else if(((!strcmp(name1, "C4")) && (!strcmp(name2, "N3"))) ||
                ((!strcmp(name2, "C4")) && (!strcmp(name1, "N3"))))
          order = 2;

        else if(((!strcmp(name1, "C5")) && (!strcmp(name2, "C6"))) ||
                ((!strcmp(name2, "C5")) && (!strcmp(name1, "C6"))))
          order = 2;
        else
          if(((!strcmp(name1, "P"))
              && (((!strcmp(name2, "O1P")) || (!strcmp(name2, "OP1")))))
             || ((!strcmp(name2, "P"))
                 && (((!strcmp(name1, "O1P")) || (!strcmp(name1, "OP1"))))))
          order = 2;
      }
      break;
    case 'D':                  /* deoxy nucleic acids */
      switch (resn1[1]) {
      case 'A':
        if(resn1[2] == 0) {
          if(((!strcmp(name1, "O2P")) || (!strcmp(name1, "OP2")))) {
            ai1->formalCharge = -1;
            ai1->chemFlag = false;
          } else if(((!strcmp(name2, "O2P")) || (!strcmp(name2, "OP2")))) {
            ai2->formalCharge = -1;
            ai2->chemFlag = false;
          }
          if(((!strcmp(name1, "C8")) && (!strcmp(name2, "N7"))) ||
             ((!strcmp(name2, "C8")) && (!strcmp(name1, "N7"))))
            order = 2;
          else if(((!strcmp(name1, "C4")) && (!strcmp(name2, "C5"))) ||
                  ((!strcmp(name2, "C4")) && (!strcmp(name1, "C5"))))
            order = 2;

          else if(((!strcmp(name1, "C6")) && (!strcmp(name2, "N1"))) ||
                  ((!strcmp(name2, "C6")) && (!strcmp(name1, "N1"))))
            order = 2;
          else if(((!strcmp(name1, "C2")) && (!strcmp(name2, "N3"))) ||
                  ((!strcmp(name2, "C2")) && (!strcmp(name1, "N3"))))
            order = 2;
          else
            if(((!strcmp(name1, "P"))
                && (((!strcmp(name2, "O1P")) || (!strcmp(name2, "OP1")))))
               || ((!strcmp(name2, "P"))
                   && (((!strcmp(name1, "O1P")) || (!strcmp(name1, "OP1"))))))
            order = 2;
        }
        break;
      case 'C':
        if(resn1[2] == 0) {
          if(((!strcmp(name1, "O2P")) || (!strcmp(name1, "OP2")))) {
            ai1->formalCharge = -1;
            ai1->chemFlag = false;
          } else if(((!strcmp(name2, "O2P")) || (!strcmp(name2, "OP2")))) {
            ai2->formalCharge = -1;
            ai2->chemFlag = false;
          }
          if(((!strcmp(name1, "C2")) && (!strcmp(name2, "O2"))) ||
             ((!strcmp(name2, "C2")) && (!strcmp(name1, "O2"))))
            order = 2;
          else if(((!strcmp(name1, "C4")) && (!strcmp(name2, "N3"))) ||
                  ((!strcmp(name2, "C4")) && (!strcmp(name1, "N3"))))
            order = 2;

          else if(((!strcmp(name1, "C5")) && (!strcmp(name2, "C6"))) ||
                  ((!strcmp(name2, "C5")) && (!strcmp(name1, "C6"))))
            order = 2;
          else
            if(((!strcmp(name1, "P"))
                && (((!strcmp(name2, "O1P")) || (!strcmp(name2, "OP1")))))
               || ((!strcmp(name2, "P"))
                   && (((!strcmp(name1, "O1P")) || (!strcmp(name1, "OP1"))))))
            order = 2;
        }
        break;
      case 'T':
        if(resn1[2] == 0) {
          if(((!strcmp(name1, "O2P")) || (!strcmp(name1, "OP2"))))
            ai1->formalCharge = -1;
          else if(((!strcmp(name2, "O2P")) || (!strcmp(name2, "OP2"))))
            ai2->formalCharge = -1;

          if(((!strcmp(name1, "C2")) && (!strcmp(name2, "O2"))) ||
             ((!strcmp(name2, "C2")) && (!strcmp(name1, "O2"))))
            order = 2;
          else if(((!strcmp(name1, "C4")) && (!strcmp(name2, "O4"))) ||
                  ((!strcmp(name2, "C4")) && (!strcmp(name1, "O4"))))
            order = 2;

          else if(((!strcmp(name1, "C5")) && (!strcmp(name2, "C6"))) ||
                  ((!strcmp(name2, "C5")) && (!strcmp(name1, "C6"))))
            order = 2;
          else
            if(((!strcmp(name1, "P"))
                && (((!strcmp(name2, "O1P")) || (!strcmp(name2, "OP1")))))
               || ((!strcmp(name2, "P"))
                   && (((!strcmp(name1, "O1P")) || (!strcmp(name1, "OP1"))))))
            order = 2;
        }
        break;
      case 'G':
        if(resn1[2] == 0) {
          if(((!strcmp(name1, "O2P")) || (!strcmp(name1, "OP2")))) {
            ai1->formalCharge = -1;
            ai1->chemFlag = false;
          } else if(((!strcmp(name2, "O2P")) || (!strcmp(name2, "OP2")))) {
            ai2->formalCharge = -1;
            ai2->chemFlag = false;
          }
          if(((!strcmp(name1, "C6")) && (!strcmp(name2, "O6"))) ||
             ((!strcmp(name2, "C6")) && (!strcmp(name1, "O6"))))
            order = 2;
          else if(((!strcmp(name1, "C2")) && (!strcmp(name2, "N3"))) ||
                  ((!strcmp(name2, "C2")) && (!strcmp(name1, "N3"))))
            order = 2;
          else if(((!strcmp(name1, "C8")) && (!strcmp(name2, "N7"))) ||
                  ((!strcmp(name2, "C8")) && (!strcmp(name1, "N7"))))
            order = 2;
          else if(((!strcmp(name1, "C4")) && (!strcmp(name2, "C5"))) ||
                  ((!strcmp(name2, "C4")) && (!strcmp(name1, "C5"))))
            order = 2;
          else
            if(((!strcmp(name1, "P"))
                && (((!strcmp(name2, "O1P")) || (!strcmp(name2, "OP1")))))
               || ((!strcmp(name2, "P"))
                   && (((!strcmp(name1, "O1P")) || (!strcmp(name1, "OP1"))))))
            order = 2;
        }
        break;
      case 'U':
        if(resn1[2] == 0) {
          if(((!strcmp(name1, "O2P")) || (!strcmp(name1, "OP2")))) {
            ai1->formalCharge = -1;
            ai1->chemFlag = false;
          } else if(((!strcmp(name2, "O2P")) || (!strcmp(name2, "OP2")))) {
            ai2->formalCharge = -1;
            ai2->chemFlag = false;
          }

          if(((!strcmp(name1, "C2")) && (!strcmp(name2, "O2"))) ||
             ((!strcmp(name2, "C2")) && (!strcmp(name1, "O2"))))
            order = 2;
          else if(((!strcmp(name1, "C4")) && (!strcmp(name2, "O4"))) ||
                  ((!strcmp(name2, "C4")) && (!strcmp(name1, "O4"))))
            order = 2;

          else if(((!strcmp(name1, "C5")) && (!strcmp(name2, "C6"))) ||
                  ((!strcmp(name2, "C5")) && (!strcmp(name1, "C6"))))
            order = 2;
          else
            if(((!strcmp(name1, "P"))
                && (((!strcmp(name2, "O1P")) || (!strcmp(name2, "OP1")))))
               || ((!strcmp(name2, "P"))
                   && (((!strcmp(name1, "O1P")) || (!strcmp(name1, "OP1"))))))
            order = 2;
        }
        break;
      }
      break;
    case 'G':
      switch (resn1[1]) {
      case 'L':
        switch (resn1[2]) {
        case 'U':              /* GLU missing GLUN, GLUH, GLH handling */
          switch (resn1[3]) {
          case 0:
          case 'M':            /* minus */
            if(!strcmp(name1, "OE2")) {
              ai1->formalCharge = -1;
              ai1->chemFlag = false;
            } else if(!strcmp(name2, "OE2")) {
              ai2->formalCharge = -1;
              ai2->chemFlag = false;
            }
            break;
          }
          if(((!strcmp(name1, "CD")) && (!strcmp(name2, "OE1"))) ||
             ((!strcmp(name2, "CD")) && (!strcmp(name1, "OE1"))))
            order = 2;
          break;
        case 'N':              /* GLN or GLU */
          if(((!strcmp(name1, "CD")) && (!strcmp(name2, "OE1"))) ||
             ((!strcmp(name2, "CD")) && (!strcmp(name1, "OE1"))))
            order = 2;
          break;
        }
        break;
      case 0:
        if(((!strcmp(name1, "O2P")) || (!strcmp(name1, "OP2")))) {
          ai1->formalCharge = -1;
          ai1->chemFlag = false;
        } else if(((!strcmp(name2, "O2P")) || (!strcmp(name2, "OP2")))) {
          ai2->formalCharge = -1;
          ai2->chemFlag = false;
        }

        if(((!strcmp(name1, "C6")) && (!strcmp(name2, "O6"))) ||
           ((!strcmp(name2, "C6")) && (!strcmp(name1, "O6"))))
          order = 2;
        else if(((!strcmp(name1, "C2")) && (!strcmp(name2, "N3"))) ||
                ((!strcmp(name2, "C2")) && (!strcmp(name1, "N3"))))
          order = 2;
        else if(((!strcmp(name1, "C8")) && (!strcmp(name2, "N7"))) ||
                ((!strcmp(name2, "C8")) && (!strcmp(name1, "N7"))))
          order = 2;
        else if(((!strcmp(name1, "C4")) && (!strcmp(name2, "C5"))) ||
                ((!strcmp(name2, "C4")) && (!strcmp(name1, "C5"))))
          order = 2;
        else
          if(((!strcmp(name1, "P"))
              && (((!strcmp(name2, "O1P")) || (!strcmp(name2, "OP1")))))
             || ((!strcmp(name2, "P"))
                 && (((!strcmp(name1, "O1P")) || (!strcmp(name1, "OP1"))))))
          order = 2;
        break;
      }
      break;
    case 'H':
      switch (resn1[1]) {
      case 'I':
        switch (resn1[2]) {
        case 'P':
          if(!strcmp(name1, "ND1")) {
            ai1->formalCharge = 1;
            ai1->chemFlag = false;
          } else if(!strcmp(name2, "ND1")) {
            ai2->formalCharge = 1;
            ai2->chemFlag = false;
          }
          if(((!strcmp(name1, "CG")) && (!strcmp(name2, "CD2"))) ||
             ((!strcmp(name2, "CG")) && (!strcmp(name1, "CD2"))))
            order = 2;
          else if(((!strcmp(name1, "CE1")) && (!strcmp(name2, "ND1"))) ||
                  ((!strcmp(name2, "CE1")) && (!strcmp(name1, "ND1"))))
            order = 2;
          break;
        case 'S':
          switch (resn1[3]) {
          case 'A':            /* HISA Gromacs */
          case 'D':            /* HISD Quanta */
            if(((!strcmp(name1, "CG")) && (!strcmp(name2, "CD2"))) ||
               ((!strcmp(name2, "CG")) && (!strcmp(name1, "CD2"))))
              order = 2;
            else if(((!strcmp(name1, "CE1")) && (!strcmp(name2, "NE2"))) ||
                    ((!strcmp(name2, "CE1")) && (!strcmp(name1, "NE2"))))
              order = 2;
            break;
          case 0:              /* plain HIS */
          case 'B':            /* HISB Gromacs */
          case 'E':            /* HISE Quanta */
            if(((!strcmp(name1, "CG")) && (!strcmp(name2, "CD2"))) ||
               ((!strcmp(name2, "CG")) && (!strcmp(name1, "CD2"))))
              order = 2;
            else if(((!strcmp(name1, "CE1")) && (!strcmp(name2, "ND1"))) ||
                    ((!strcmp(name2, "CE1")) && (!strcmp(name1, "ND1"))))
              order = 2;
            break;
          case 'H':            /* HISH Gromacs */
          case 'P':            /* HISP Quanta */
            if(!strcmp(name1, "ND1")) {
              ai1->formalCharge = 1;
              ai1->chemFlag = false;
            } else if(!strcmp(name2, "ND1")) {
              ai2->formalCharge = 1;
              ai2->chemFlag = false;
            }
            if(((!strcmp(name1, "CG")) && (!strcmp(name2, "CD2"))) ||
               ((!strcmp(name2, "CG")) && (!strcmp(name1, "CD2"))))
              order = 2;
            else if(((!strcmp(name1, "CE1")) && (!strcmp(name2, "ND1"))) ||
                    ((!strcmp(name2, "CE1")) && (!strcmp(name1, "ND1"))))
              order = 2;
            break;
          }
          break;
        case 'E':              /* HIE */
          if(((!strcmp(name1, "CG")) && (!strcmp(name2, "CD2"))) ||
             ((!strcmp(name2, "CG")) && (!strcmp(name1, "CD2"))))
            order = 2;
          else if(((!strcmp(name1, "CE1")) && (!strcmp(name2, "ND1"))) ||
                  ((!strcmp(name2, "CE1")) && (!strcmp(name1, "ND1"))))
            order = 2;
          break;
        case 'D':              /* HID */
          if(((!strcmp(name1, "CG")) && (!strcmp(name2, "CD2"))) ||
             ((!strcmp(name2, "CG")) && (!strcmp(name1, "CD2"))))
            order = 2;
          else if(((!strcmp(name1, "CE1")) && (!strcmp(name2, "NE2"))) ||
                  ((!strcmp(name2, "CE1")) && (!strcmp(name1, "NE2"))))
            order = 2;
          break;
        }
        break;
      }
      break;
    case 'I':
      if(resn1[1] == 0) {
        if(((!strcmp(name1, "O2P")) || (!strcmp(name1, "OP2")))) {
          ai1->formalCharge = -1;
          ai1->chemFlag = false;
        } else if(((!strcmp(name2, "O2P")) || (!strcmp(name2, "OP2")))) {
          ai2->formalCharge = -1;
          ai2->chemFlag = false;
        }
        if(((!strcmp(name1, "C8")) && (!strcmp(name2, "N7"))) ||
           ((!strcmp(name2, "C8")) && (!strcmp(name1, "N7"))))
          order = 2;
        else if(((!strcmp(name1, "C4")) && (!strcmp(name2, "C5"))) ||
                ((!strcmp(name2, "C4")) && (!strcmp(name1, "C5"))))
          order = 2;

        else if(((!strcmp(name1, "C6")) && (!strcmp(name2, "N1"))) ||
                ((!strcmp(name2, "C6")) && (!strcmp(name1, "N1"))))
          order = 2;
        else if(((!strcmp(name1, "C2")) && (!strcmp(name2, "N3"))) ||
                ((!strcmp(name2, "C2")) && (!strcmp(name1, "N3"))))
          order = 2;
        else
          if(((!strcmp(name1, "P"))
              && (((!strcmp(name2, "O1P")) || (!strcmp(name2, "OP1")))))
             || ((!strcmp(name2, "P"))
                 && (((!strcmp(name1, "O1P")) || (!strcmp(name1, "OP1"))))))
          order = 2;
      }
      break;
    case 'P':
      switch (resn1[1]) {
      case 'H':                /* PHE */
        if(resn1[2] == 'E') {
          if(((!strcmp(name1, "CG")) && (!strcmp(name2, "CD1"))) ||
             ((!strcmp(name2, "CG")) && (!strcmp(name1, "CD1"))))
            order = 2;
          else if(((!strcmp(name1, "CZ")) && (!strcmp(name2, "CE1"))) ||
                  ((!strcmp(name2, "CZ")) && (!strcmp(name1, "CE1"))))
            order = 2;

          else if(((!strcmp(name1, "CE2")) && (!strcmp(name2, "CD2"))) ||
                  ((!strcmp(name2, "CE2")) && (!strcmp(name1, "CD2"))))
            order = 2;
          break;
        }
      }
      break;
    case 'L':
      switch (resn1[1]) {
      case 'Y':
        switch (resn1[2]) {
        case 'S':              /* LYS. */
          switch (resn1[3]) {
          case 0:
          case 'P':            /* LYS, LYSP */
            if(!strcmp(name1, "NZ")) {
              ai1->formalCharge = 1;
              ai1->chemFlag = false;
            } else if(!strcmp(name2, "NZ")) {
              ai2->formalCharge = 1;
              ai2->chemFlag = false;
            }
            break;
          }
          break;
        }
        break;
      }
      break;
    case 'T':
      switch (resn1[1]) {
      case 'Y':                /* TYR */
        if(resn1[2] == 'R') {
          if(((!strcmp(name1, "CG")) && (!strcmp(name2, "CD1"))) ||
             ((!strcmp(name2, "CG")) && (!strcmp(name1, "CD1"))))
            order = 2;
          else if(((!strcmp(name1, "CZ")) && (!strcmp(name2, "CE1"))) ||
                  ((!strcmp(name2, "CZ")) && (!strcmp(name1, "CE1"))))
            order = 2;

          else if(((!strcmp(name1, "CE2")) && (!strcmp(name2, "CD2"))) ||
                  ((!strcmp(name2, "CE2")) && (!strcmp(name1, "CD2"))))
            order = 2;
          break;
        }
        break;
      case 'R':
        if(resn1[2] == 'P') {
          if(((!strcmp(name1, "CG")) && (!strcmp(name2, "CD1"))) ||
             ((!strcmp(name2, "CG")) && (!strcmp(name1, "CD1"))))
            order = 2;
          else if(((!strcmp(name1, "CZ3")) && (!strcmp(name2, "CE3"))) ||
                  ((!strcmp(name2, "CZ3")) && (!strcmp(name1, "CE3"))))
            order = 2;
          else if(((!strcmp(name1, "CZ2")) && (!strcmp(name2, "CH2"))) ||
                  ((!strcmp(name2, "CZ2")) && (!strcmp(name1, "CH2"))))
            order = 2;
          else if(((!strcmp(name1, "CE2")) && (!strcmp(name2, "CD2"))) ||
                  ((!strcmp(name2, "CE2")) && (!strcmp(name1, "CD2"))))
            order = 2;
          break;
        }
        break;
      case 0:
        if(((!strcmp(name1, "O2P")) || (!strcmp(name1, "OP2")))) {
          ai1->formalCharge = -1;
          ai1->chemFlag = false;
        } else if(((!strcmp(name2, "O2P")) || (!strcmp(name2, "OP2")))) {
          ai2->formalCharge = -1;
          ai2->chemFlag = false;
        }

        if(((!strcmp(name1, "C2")) && (!strcmp(name2, "O2"))) ||
           ((!strcmp(name2, "C2")) && (!strcmp(name1, "O2"))))
          order = 2;
        else if(((!strcmp(name1, "C4")) && (!strcmp(name2, "O4"))) ||
                ((!strcmp(name2, "C4")) && (!strcmp(name1, "O4"))))
          order = 2;

        else if(((!strcmp(name1, "C5")) && (!strcmp(name2, "C6"))) ||
                ((!strcmp(name2, "C5")) && (!strcmp(name1, "C6"))))
          order = 2;
        else
          if(((!strcmp(name1, "P"))
              && (((!strcmp(name2, "O1P")) || (!strcmp(name2, "OP1")))))
             || ((!strcmp(name2, "P"))
                 && (((!strcmp(name1, "O1P")) || (!strcmp(name1, "OP1"))))))
          order = 2;
        break;
      }
      break;
    case 'U':
      if(resn1[1] == 0) {
        if(((!strcmp(name1, "O2P")) || (!strcmp(name1, "OP2")))) {
          ai1->formalCharge = -1;
          ai1->chemFlag = false;
        } else if(((!strcmp(name2, "O2P")) || (!strcmp(name2, "OP2")))) {
          ai2->formalCharge = -1;
          ai2->chemFlag = false;
        }
        if(((!strcmp(name1, "C2")) && (!strcmp(name2, "O2"))) ||
           ((!strcmp(name2, "C2")) && (!strcmp(name1, "O2"))))
          order = 2;
        else if(((!strcmp(name1, "C4")) && (!strcmp(name2, "O4"))) ||
                ((!strcmp(name2, "C4")) && (!strcmp(name1, "O4"))))
          order = 2;

        else if(((!strcmp(name1, "C5")) && (!strcmp(name2, "C6"))) ||
                ((!strcmp(name2, "C5")) && (!strcmp(name1, "C6"))))
          order = 2;
        else
          if(((!strcmp(name1, "P"))
              && (((!strcmp(name2, "O1P")) || (!strcmp(name2, "OP1")))))
             || ((!strcmp(name2, "P"))
                 && (((!strcmp(name1, "O1P")) || (!strcmp(name1, "OP1"))))))
          order = 2;
      }
      break;
    }
  }
  *(bond_order) = order;
}

void ObjectMoleculeFixChemistry(ObjectMolecule * I, int sele1, int sele2, int invalidate)
{
  int b;
  int flag = false;
  int s1, s2;
  AtomInfoType *ai1, *ai2;
  int order;
  BondType *bond;
  bond = I->Bond;
  for(b = 0; b < I->NBond; b++) {
    flag = false;
    ai1 = I->AtomInfo + bond->index[0];
    ai2 = I->AtomInfo + bond->index[1];
    s1 = ai1->selEntry;
    s2 = ai2->selEntry;

    if((SelectorIsMember(I->Obj.G, s1, sele1) &&
        SelectorIsMember(I->Obj.G, s2, sele2)) ||
       (SelectorIsMember(I->Obj.G, s2, sele1) && SelectorIsMember(I->Obj.G, s1, sele2))) {
      order = -1;
      if(!ai1->resn[3]) {       /* Standard disconnected PDB residue */
        if(AtomInfoSameResidue(I->Obj.G, ai1, ai2)) {
          assign_pdb_known_residue(I->Obj.G, ai1, ai2, &order);
        }
      }
      if(order > 0) {
        bond->order = order;
        ai1->chemFlag = false;
        ai2->chemFlag = false;
        flag = true;
      } else if(invalidate) {
        ai1->chemFlag = false;
        ai2->chemFlag = false;
        flag = true;
      }
    }
    bond++;
  }
  if(flag) {
    ObjectMoleculeInvalidate(I, cRepAll, cRepInvAll, -1);
    SceneChanged(I->Obj.G);
  }
}

void ObjMolPairwiseInit(ObjMolPairwise * pairwise)
{
  UtilZeroMem((char *) pairwise, sizeof(ObjMolPairwise));
  pairwise->trg_vla = VLAlloc(int, 10);
  pairwise->mbl_vla = VLAlloc(int, 10);
}

void ObjMolPairwisePurge(ObjMolPairwise * pairwise)
{
  VLAFreeP(pairwise->trg_vla);
  VLAFreeP(pairwise->mbl_vla);
}

int ObjectMoleculeConvertIDsToIndices(ObjectMolecule * I, int *id, int n_id)
{
  /* return true if all IDs are unique, false if otherwise */

  int min_id, max_id, range, *lookup = NULL;
  int unique = true;

  /* this routine only works if IDs cover a reasonable range --
     should rewrite using a hash table */

  if(I->NAtom) {

    /* determine range */

    {
      int a, cur_id;
      cur_id = I->AtomInfo[0].id;
      min_id = cur_id;
      max_id = cur_id;
      for(a = 1; a < I->NAtom; a++) {
        cur_id = I->AtomInfo[a].id;
        if(min_id > cur_id)
          min_id = cur_id;
        if(max_id < cur_id)
          max_id = cur_id;
      }
    }

    /* create cross-reference table */

    {
      int a, offset;

      range = max_id - min_id + 1;
      lookup = Calloc(int, range);
      for(a = 0; a < I->NAtom; a++) {
        offset = I->AtomInfo[a].id - min_id;
        if(!lookup[offset])
          lookup[offset] = a + 1;
        else
          unique = false;
      }
    }

    /* iterate through IDs and replace with indices or -1 */

    {
      int i, offset, lkup;

      for(i = 0; i < n_id; i++) {
        offset = id[i] - min_id;
        if((offset >= 0) && (offset < range)) {
          lkup = lookup[offset];
          if(lkup > 0) {
            id[i] = lkup - 1;
          } else {
            id[i] = -1;         /* negative means no match */
          }
        } else
          id[i] = -1;
      }
    }
  }

  FreeP(lookup);
  return unique;

}

static char *check_next_pdb_object(char *p, int skip_to_next)
{
  char *start = p;
  char *line_start;
  char cc[MAXLINELEN];
  while(*p) {
    line_start = p;
    p = ncopy(cc, p, 6);
    if((cc[0] == 'H') &&
       (cc[1] == 'E') &&
       (cc[2] == 'A') && (cc[3] == 'D') && (cc[4] == 'E') && (cc[5] == 'R')) {
      if(skip_to_next)
        return line_start;
      else
        return start;
    } else if(((cc[0] == 'A') &&        /* ATOM */
               (cc[1] == 'T') && (cc[2] == 'O') && (cc[3] == 'M') && (cc[4] == ' ') && (cc[5] == ' ')) || ((cc[0] == 'H') &&    /* HETATM */
                                                                                                           (cc[1] == 'E') && (cc[2] == 'T') && (cc[3] == 'A') && (cc[4] == 'T') && (cc[5] == 'M'))) {
      p = nskip(p, 5);
      ParseNTrim(cc, p, 14);
      /* this is a special workaround for the bogus HETATM entry in PDB 4ZN2:
         HETATM20829              0       0.000   0.000   0.000  1.00  0.00      4ZN20773
         which screws up our PDB test case */
      if(cc[0]) {
        return start;
      }
    } else if(((cc[0] == 'E') &&        /* if we pass over the END of the current PDB file, then reset start */
               (cc[1] == 'N') && (cc[2] == 'D'))) {
      if(strcmp("END", cc) == 0) {      /* END */
        if(skip_to_next) {
          start = line_start;
        }
      }
    }
    p = nextline(p);
  }
  return NULL;
}

static int get_multi_object_status(char *p)
{                               /* expensive -- only call
                                   this if there is no
                                   other way to determine
                                   this information */
  int seen_header = 0;
  char cc[MAXLINELEN];
  while(*p) {
    p = ncopy(cc, p, 6);
    if((cc[0] == 'H') &&
       (cc[1] == 'E') &&
       (cc[2] == 'A') && (cc[3] == 'D') && (cc[4] == 'E') && (cc[5] == 'R')) {
      if(seen_header) {
        return 1;
      } else {
        seen_header = true;
      }
    }
    p = nextline(p);
  }
  return -1;
}

int ObjectMoleculeAutoDisableAtomNameWildcard(ObjectMolecule * I)
{
  PyMOLGlobals *G = I->Obj.G;
  char wildcard = 0;
  int found_wildcard = false;

  {
    char *tmp = SettingGet_s(G, NULL, I->Obj.Setting, cSetting_atom_name_wildcard);
    if(tmp && tmp[0]) {
      wildcard = *tmp;
    } else {
      tmp = SettingGet_s(G, NULL, I->Obj.Setting, cSetting_wildcard);
      if(tmp) {
        wildcard = *tmp;
      }
    }
    if(wildcard == 32)
      wildcard = 0;

  }

  if(wildcard) {
    register int a;
    register char *p, ch;
    register AtomInfoType *ai = I->AtomInfo;

    for(a = 0; a < I->NAtom; a++) {
      p = ai->name;
      while((ch = *(p++))) {
        if(ch == wildcard) {
          found_wildcard = true;
          break;
        }
      }
      ai++;
    }
    if(found_wildcard) {
      ExecutiveSetObjSettingFromString(G, cSetting_atom_name_wildcard, " ",
                                       &I->Obj, -1, true, true);
    }
  }
  return found_wildcard;
}


/*========================================================================*/
#define PDB_MAX_TAGS 64

CoordSet *ObjectMoleculePDBStr2CoordSet(PyMOLGlobals * G,
                                        char *buffer,
                                        AtomInfoType ** atInfoPtr,
                                        char **restart_model,
                                        char *segi_override,
                                        M4XAnnoType * m4x,
                                        char *pdb_name,
                                        char **next_pdb,
                                        PDBInfoRec * info, int quiet, int *model_number)
{

  register char *p;
  int nAtom;
  int a, b, c;
  float *coord = NULL;
  CoordSet *cset = NULL;
  AtomInfoType *atInfo = NULL, *ai;
  int AFlag;
  char SSCode;
  int atomCount;
  int bondFlag = false;
  BondType *bond = NULL, *ii1, *ii2;
  int *idx;
  int nBond = 0;
  int b1, b2, nReal, maxAt;
  CSymmetry *symmetry = NULL;
  int auto_show_lines = (int) SettingGet(G, cSetting_auto_show_lines);
  int auto_show_nonbonded = (int) SettingGet(G, cSetting_auto_show_nonbonded);
  int auto_show_spheres = (int) SettingGet(G, cSetting_auto_show_spheres);
  int reformat_names = (int) SettingGet(G, cSetting_pdb_reformat_names_mode);
  int truncate_resn = SettingGetGlobal_b(G, cSetting_pdb_truncate_residue_name);
  char *tags_in = SettingGetGlobal_s(G, cSetting_pdb_echo_tags), *tag_start[PDB_MAX_TAGS];
  int n_tags = 0;
  int foundNextModelFlag = false;
  int ssFlag = false;
  int ss_resv1 = 0, ss_resv2 = 0;
  ResIdent ss_resi1 = "", ss_resi2 = "";
  unsigned char ss_chain1 = 0, ss_chain2 = 0;
  SSEntry *ss_list = NULL;
  int n_ss = 1;
  int *(ss[256]);               /* one array for each chain identifier */

  char cc[MAXLINELEN], tags[MAXLINELEN];
  char cc_saved, ctmp;
  int index;
  int ignore_pdb_segi = 0;
  int ss_valid, ss_found = false;
  SSEntry *sst;
  int ssi = 0;
  int only_read_one_model = false;
  int ignore_conect = SettingGetGlobal_b(G, cSetting_pdb_ignore_conect);
  int have_bond_order = false;
  int seen_model, in_model = false;
  int seen_conect = false;
  int have_conect = false;
  int is_end_of_object = false;
  int literal_names = SettingGetGlobal_b(G, cSetting_pdb_literal_names);
  int bogus_name_alignment = true;
  AtomName literal_name = "";

  if(tags_in && (!quiet) && (!*restart_model)) {
    char *p = tags;
    strcpy(tags, tags_in);

    while(*p) {
      while(*p == ' ')          /* skip spaces */
        p++;
      if(*p) {
        tag_start[n_tags] = p;
        n_tags++;
        while(*p) {
          if(*p != ',') {
            if(*p == ' ')
              *p = 0;
            p++;
          } else
            break;
        }
        if(*p) {                /* terminate tag */
          *p = 0;
          p++;
        }
      }
    }
  }

  if(literal_names)
    reformat_names = 0;

  ignore_pdb_segi = (int) SettingGet(G, cSetting_ignore_pdb_segi);

  p = buffer;
  nAtom = 0;
  if(atInfoPtr)
    atInfo = *atInfoPtr;

  if(!atInfo)
    ErrFatal(G, "PDBStr2CoordSet", "need atom information record!");    /* failsafe for old version.. */

  if(buffer == *restart_model)
    only_read_one_model = true;
  else if(info && (info->multiplex > 0)) {
    only_read_one_model = true;
    if(!info->multi_object_status) {    /* figure out if this is a multi-object (multi-HEADER) pdb file */
      info->multi_object_status = get_multi_object_status(p);
    }
  }
  /* PASS 1 */
  {
    int seen_end_of_atoms = false;

    *restart_model = NULL;
    while(*p) {
      if(n_tags && !quiet) {
        int skip = false;

        if((p[0] == 'H') &&
           (p[1] == 'E') &&
           (p[2] == 'A') && (p[3] == 'D') && (p[4] == 'E') && (p[5] == 'R')) {
          if(nAtom > 0) {
            /* don't print HEADER until next time */
            skip = true;
          }
        }
        if(!skip) {
          /* fast unrolled string match */
          register int tc = 0;
          register char *q;
          register int same;
          while(tc < n_tags) {
            same = true;
            q = tag_start[tc];
            if(p[0] != q[0])
              same = false;
            else if(p[0] && q[0]) {
              if((p[1] != q[1]) && !((p[1] == ' ') && !q[1]))
                same = false;
              else if(p[1] && q[1]) {
                if((p[2] != q[2]) && !((p[2] == ' ') && !q[2]))
                  same = false;
                else if(p[3] && q[3]) {
                  if((p[3] != q[3]) && !((p[3] == ' ') && !q[3]))
                    same = false;
                  else if(p[4] && q[4]) {
                    if((p[4] != q[4]) && !((p[4] == ' ') && !q[4]))
                      same = false;
                    else if(p[5] && q[5]) {
                      if((p[5] != q[5]) && !((p[5] == ' ') && !q[5]))
                        same = false;
                    }
                  }
                }
              }
            }
            if(same) {
              ParseNTrimRight(cc, p, MAXLINELEN - 1);
              /*              OrthoAddOutput(G," PDB: "); */
              OrthoAddOutput(G, cc);
              OrthoNewLine(G, NULL, true);
            }
            tc++;
          }
        }
      }
      if(((p[0] == 'A') && (p[1] == 'T') && (p[2] == 'O') && (p[3] == 'M')) ||  /* ATOM */
         ((p[0] == 'H') && (p[1] == 'E') && (p[2] == 'T') &&    /* HETATM */
          (p[3] == 'A') && (p[4] == 'T') && (p[5] == 'M') && (!*restart_model))) {
        if(!seen_end_of_atoms)
          nAtom++;
        if(bogus_name_alignment) {
          ncopy(cc, nskip(p, 12), 4);   /* copy the atom field */
          if((cc[0] == 32) && (cc[1] != 32)) {  /* check to see if indentation was followed correctly */
            bogus_name_alignment = false;
          }
        }
      } else if((p[0] == 'H') && (p[1] == 'E') && (p[2] == 'L') && (p[3] == 'I') && (p[4] == 'X'))      /* HELIX */
        ssFlag = true;
      else if((p[0] == 'S') && (p[1] == 'H') && (p[2] == 'E') && (p[3] == 'E') && (p[4] == 'T'))        /* SHEET */
        ssFlag = true;
      else if((p[0] == 'H') &&  /* HEADER */
              (p[1] == 'E') &&
              (p[2] == 'A') && (p[3] == 'D') && (p[4] == 'E') && (p[5] == 'R')) {
        if(nAtom > 0) {         /* if we've already found atom records, then this must be a new pdb */
          (*next_pdb) = p;
          break;
        }
      } else if((p[0] == 'R') && (p[1] == 'E') && (p[2] == 'M') &&      /* REMARK */
                (p[3] == 'A') && (p[4] == 'R') && (p[5] == 'K')) {
        ntrim(cc, p, 30);
        if(strcmp("REMARK    GENERATED BY TRJCONV", cc) == 0)
          if(info)
            info->ignore_header_names = true;
      } else if((p[0] == 'E') && (p[1] == 'N') && (p[2] == 'D') &&      /* ENDMDL */
                (p[3] == 'M') && (p[4] == 'D') && (p[5] == 'L') && (!*restart_model)) {
        *restart_model = nextline(p);
        seen_end_of_atoms = true;
        if(only_read_one_model)
          break;
      } else if((p[0] == 'E') && (p[1] == 'N') && (p[2] == 'D')) {      /* stop parsing after END */
        ntrim(cc, p, 6);
        if(strcmp("END", cc) == 0) {    /* END */
          seen_end_of_atoms = true;
          if(next_pdb) {
            p = nextline(p);
            ncopy(cc, p, 6);
            if(strcmp("HEADER", cc) == 0) {
              (*next_pdb) = p;  /* found another PDB file after this one... */
            } else if(strcmp("ENDMDL", cc) == 0) {
              seen_end_of_atoms = false;
            }
          }
          break;
        }
      } else if((p[0] == 'C') && (p[1] == 'O') && (p[2] == 'N') && (p[3] == 'E') && (p[4] == 'C') && (p[5] == 'T')) {   /* CONECT */
        have_conect = true;
        if(!ignore_conect) 
          bondFlag = true;
      } else if((p[0] == 'U') && (p[1] == 'S') && (p[2] == 'E') &&
                (p[3] == 'R') && (!*restart_model)) {

        /* Metaphorics key now 'USER M4X ', changed from 'USER    ' */
        if((p[4] == ' ') && (p[5] == 'M') && (p[6] == '4') &&
           (p[7] == 'X') && (p[8] == ' ') && m4x) {
          p = nskip(p, 10);
          p = ntrim(cc, p, 6);
          m4x->annotated_flag = true;
          switch (cc[0]) {
          case 'H':
            if(WordMatchExact(G, "HINT", cc, true)) {
              p = nskip(p, 1);
              p = ntrim(cc, p, 6);      /* get context name */
              if(WordMatchExact(G, "ALIGN", cc, true)) {        /* ALIGN is special */
                if(!m4x->align) {
                  m4x->align = Calloc(M4XAlignType, 1);
                  M4XAlignInit(m4x->align);
                  p = nskip(p, 8);
                  p = ntrim(cc, p, 6);  /* get visibility of this structure */
                }
              } else if(WordMatchExact(G, "HIDE", cc, true)) {
                m4x->invisible = 1;
              } else {
                if(!m4x->context) {
                  m4x->context = VLACalloc(M4XContextType, 10);
                }
                if(m4x->context) {
                  int cn;
                  int found = false;

                  /* does context already exist ? */
                  for(cn = 0; cn < m4x->n_context; cn++) {
                    if(WordMatchExact(G, m4x->context[cn].name, cc, true)) {
                      found = true;
                      break;
                    }
                  }

                  /* if not, then create it */
                  if(!found) {
                    cn = m4x->n_context++;
                    VLACheck(m4x->context, M4XContextType, cn);
                    UtilNCopy(m4x->context[cn].name, cc, sizeof(WordType));
                  }

                  while(*cc) {
                    p = nskip(p, 1);
                    p = ntrim(cc, p, 6);
                    switch (cc[0]) {
                    case 'B':
                      if(WordMatchExact(G, "BORDER", cc, true)) {
                        /* ignore PDB CONECT if BORDER present */
                        ignore_conect = true;
                        have_bond_order = true;
                        bondFlag = true;
                      }
                      break;
                    case 'S':
                      if(WordMatchExact(G, "SITE", cc, true)) {
                        if(!m4x->context[cn].site) {
                          m4x->context[cn].site = VLAlloc(int, 50);
                        }
                      }
                      break;
                    case 'L':
                      if(WordMatchExact(G, "LIGAND", cc, true)) {
                        if(!m4x->context[cn].ligand) {
                          m4x->context[cn].ligand = VLAlloc(int, 50);
                        }
                      }
                      break;
                    case 'W':
                      if(WordMatchExact(G, "WATER", cc, true)) {
                        if(!m4x->context[cn].water) {
                          m4x->context[cn].water = VLAlloc(int, 50);
                        }
                      }
                      break;
                    case 'H':
                      if(WordMatchExact(G, "HBOND", cc, true)) {
                        if(!m4x->context[cn].hbond) {
                          m4x->context[cn].hbond = VLAlloc(M4XBondType, 50);
                        }
                      }
                      break;
                    case 'N':
                      if(WordMatchExact(G, "NBOND", cc, true)) {
                        if(!m4x->context[cn].nbond) {
                          m4x->context[cn].nbond = VLAlloc(M4XBondType, 50);
                        }
                      }
                      break;
                    }
                  }
                }
              }
            }
          }
        }
      }
      p = nextline(p);
    }
  }

  *restart_model = NULL;
  coord = VLAlloc(float, 3 * nAtom);

  if(atInfo)
    VLACheck(atInfo, AtomInfoType, nAtom);

  if(bondFlag) {
    nBond = 0;
    bond = VLACalloc(BondType, 6 * nAtom);
  }
  p = buffer;
  PRINTFB(G, FB_ObjectMolecule, FB_Blather)
    " ObjectMoleculeReadPDB: Found %i atoms...\n", nAtom ENDFB(G);

  if(ssFlag) {
    for(a = 0; a <= 255; a++) {
      ss[a] = 0;
    }
    ss_list = VLAlloc(SSEntry, 50);
  }

  a = 0;                        /* WATCHOUT */
  atomCount = 0;

  /* PASS 2 */
  seen_model = false;

  while(*p) {
    /*      printf("%c%c%c%c%c%c\n",p[0],p[1],p[2],p[3],p[4],p[5]); */
    AFlag = false;
    SSCode = 0;
    if((p[0] == 'A') &&         /* ATOM */
       (p[1] == 'T') && (p[2] == 'O') && (p[3] == 'M'))
      AFlag = 1;
    else if((p[0] == 'H') &&    /* HETATM */
            (p[1] == 'E') &&
            (p[2] == 'T') && (p[3] == 'A') && (p[4] == 'T') && (p[5] == 'M'))
      AFlag = 2;
    else if((p[0] == 'H') &&    /* HEADER */
            (p[1] == 'E') &&
            (p[2] == 'A') && (p[3] == 'D') && (p[4] == 'E') && (p[5] == 'R')) {

      if(pdb_name) {
        if(atomCount > 0) {
          /* have we already read atoms?  then this is probably a new PDB file */
          (*next_pdb) = p;      /* found another PDB file after this one... */
          break;
        } else if((!info) || (!info->ignore_header_names)) {
          /* if this isn't an MD trajectory... */
          char *pp;
          pp = nskip(p, 62);    /* is there a four-letter PDB code? */
          pp = ntrim(pdb_name, pp, 4);
          if(pdb_name[0] == 0) {        /* if not, is there a plain name (for MERCK!) */
            p = nskip(p, 6);
            p = ntrim(cc, p, 44);
            UtilNCopy(pdb_name, cc, WordLength);
          } else {
            p = pp;
          }
        }
      }
    } else if((p[0] == 'M') &&  /* MODEL */
              (p[1] == 'O') && (p[2] == 'D') && (p[3] == 'E') && (p[4] == 'L')) {
      if(model_number) {
        int tmp;
        p = nskip(p, 10);
        p = ncopy(cc, p, 5);
        if(sscanf(cc, "%d", &tmp) == 1)
          *model_number = tmp;
      }
      seen_model = true;
      in_model = true;
    } else if((p[0] == 'H') &&  /* HELIX */
              (p[1] == 'E') && (p[2] == 'L') && (p[3] == 'I') && (p[4] == 'X')) {
      ss_valid = true;

      p = nskip(p, 19);
      ss_chain1 = (*p);
      p = nskip(p, 2);
      p = ncopy(cc, p, 4);
      if(!sscanf(cc, "%s", ss_resi1))
        ss_valid = false;
      if(!sscanf(cc, "%d", &ss_resv1))
        ss_valid = false;

      p = nskip(p, 6);
      ss_chain2 = (*p);
      p = nskip(p, 2);
      p = ncopy(cc, p, 4);

      if(!sscanf(cc, "%s", ss_resi2))
        ss_valid = false;
      if(!sscanf(cc, "%d", &ss_resv2))
        ss_valid = false;

      if(ss_valid) {
        PRINTFB(G, FB_ObjectMolecule, FB_Blather)
          " ObjectMolecule: read HELIX %c %s %c %s\n",
          ss_chain1, ss_resi1, ss_chain2, ss_resi2 ENDFB(G);
        SSCode = 'H';
      }

      if(ss_chain1 == ' ')
        ss_chain1 = 0;
      if(ss_chain2 == ' ')
        ss_chain2 = 0;

    } else if((p[0] == 'S') &&  /* SHEET */
              (p[1] == 'H') && (p[2] == 'E') && (p[3] == 'E') && (p[4] == 'T')) {
      ss_valid = true;
      p = nskip(p, 21);
      ss_chain1 = (*p);
      p = nskip(p, 1);
      p = ncopy(cc, p, 4);
      if(!sscanf(cc, "%s", ss_resi1))
        ss_valid = false;
      if(!sscanf(cc, "%d", &ss_resv1))
        ss_valid = false;
      p = nskip(p, 6);
      ss_chain2 = (*p);
      p = nskip(p, 1);
      p = ncopy(cc, p, 4);
      if(!sscanf(cc, "%s", ss_resi2))
        ss_valid = false;
      if(!sscanf(cc, "%d", &ss_resv2))
        ss_valid = false;

      if(ss_valid) {
        PRINTFB(G, FB_ObjectMolecule, FB_Blather)
          " ObjectMolecule: read SHEET %c %s %c %s\n",
          ss_chain1, ss_resi1, ss_chain2, ss_resi2 ENDFB(G);
        SSCode = 'S';
      }

      if(ss_chain1 == ' ')
        ss_chain1 = 0;
      if(ss_chain2 == ' ')
        ss_chain2 = 0;

    } else if((p[0] == 'E') &&  /* ENDMDL */
              (p[1] == 'N') &&
              (p[2] == 'D') && (p[3] == 'M') && (p[4] == 'D') && (p[5] == 'L')) {
      if(*restart_model)
        in_model = false;
      else {
        *restart_model = nextline(p);
        in_model = false;
        if(only_read_one_model) {
          char *pp;
          pp = nextline(p);
          if((pp[0] == 'E') &&  /* END we're going to be starting a new object... */
             (pp[1] == 'N') && (pp[2] == 'D')) {
            (*next_pdb) = check_next_pdb_object(nextline(pp), true);
            if(info && (info->multiplex == 0)) {
              /* multiplex == 0:  FORCED multimodel behavior with concatenated PDB files */
              (*restart_model) = (*next_pdb);
              (*next_pdb) = NULL;
              foundNextModelFlag = true;
              info->multi_object_status = -1;
            } else {
              is_end_of_object = true;
            }
          } else if((pp[0] == 'M') &&   /* not a new object...just a new state (model) */
                    (pp[1] == 'O') && (pp[2] == 'D') && (pp[3] == 'E') && (pp[4] == 'L')) {
            if(info && (info->multiplex > 0)) { /* end object if we're multiplexing */
              (*next_pdb) = check_next_pdb_object(pp, true);
              (*restart_model) = NULL;
            } else
              is_end_of_object = false;
          } else {
            if(pp[0] > 32)      /* more content follows... */
              (*next_pdb) = check_next_pdb_object(pp, true);
            else
              (*next_pdb) = NULL;       /* at end of file */
          }
          break;
        }
      }
    } else if((p[0] == 'E') &&  /* END */
              (p[1] == 'N') && (p[2] == 'D')) {
      ntrim(cc, p, 6);
      if((strcmp("END", cc) == 0) && (!in_model)) {     /* stop parsing here... */
        char *pp;
        pp = nextline(p);       /* ...unless we're in MODEL or next line is itself ENDMDL */
        p = ncopy(cc, p, 6);
        if(!((cc[0] == 'E') && (cc[1] == 'N') && (cc[2] == 'D') && (cc[3] == 'M') && (cc[4] == 'D') && (cc[5] == 'L'))) {       /* NOTE: this test seems unnecessary given strcmp above... */
          if(!*next_pdb) {
            (*next_pdb) = check_next_pdb_object(pp, false);
          }
          if((*next_pdb) && info && (!info->multiplex) && !(*restart_model)) {
            /* multiplex == 0:  FORCED multimodel behavior with concatenated PDB files */
            (*restart_model) = (*next_pdb);
            (*next_pdb) = NULL;
            foundNextModelFlag = true;
            info->multi_object_status = -1;
            is_end_of_object = false;
            break;
          }
          if(*next_pdb) {       /* we've found another object... */
            if(*restart_model)
              is_end_of_object = false; /* however, if we're parsing multi-models, then we're not yet at the end */
            else
              is_end_of_object = true;
            break;
          } else if(!seen_model)
            break;
        }
      }
    } else if((p[0] == 'C') &&  /* CRYST1 */
              (p[1] == 'R') &&
              (p[2] == 'Y') &&
              (p[3] == 'S') && (p[4] == 'T') && (p[5] == '1') && (!*restart_model)) {
      if(!symmetry)
        symmetry = SymmetryNew(G);
      if(symmetry) {
        int symFlag = true;
        PRINTFB(G, FB_ObjectMolecule, FB_Blather)
          " PDBStrToCoordSet: Attempting to read symmetry information\n" ENDFB(G);
        p = nskip(p, 6);
        symFlag = true;
        p = ncopy(cc, p, 9);
        if(sscanf(cc, "%f", &symmetry->Crystal->Dim[0]) != 1)
          symFlag = false;
        p = ncopy(cc, p, 9);
        if(sscanf(cc, "%f", &symmetry->Crystal->Dim[1]) != 1)
          symFlag = false;
        p = ncopy(cc, p, 9);
        if(sscanf(cc, "%f", &symmetry->Crystal->Dim[2]) != 1)
          symFlag = false;
        p = ncopy(cc, p, 7);
        if(sscanf(cc, "%f", &symmetry->Crystal->Angle[0]) != 1)
          symFlag = false;
        p = ncopy(cc, p, 7);
        if(sscanf(cc, "%f", &symmetry->Crystal->Angle[1]) != 1)
          symFlag = false;
        p = ncopy(cc, p, 7);
        if(sscanf(cc, "%f", &symmetry->Crystal->Angle[2]) != 1)
          symFlag = false;
        p = nskip(p, 1);
        p = ncopy(symmetry->SpaceGroup, p, 11);
        UtilCleanStr(symmetry->SpaceGroup);
        p = ncopy(cc, p, 3);
        if(sscanf(cc, "%d", &symmetry->PDBZValue) != 1)
          symmetry->PDBZValue = 1;
        if(!symFlag) {
          ErrMessage(G, "PDBStrToCoordSet", "Error reading CRYST1 record\n");
          SymmetryFree(symmetry);
          symmetry = NULL;
        }
      }
    } else if((p[0] == 'S') &&  /* SCALEn */
              (p[1] == 'C') &&
              (p[2] == 'A') &&
              (p[3] == 'L') && (p[4] == 'E') && (!*restart_model) && info) {
      switch (p[5]) {
      case '1':
      case '2':
      case '3':
        {
          int flag = (p[5] - '1');
          int offset = flag * 4;
          int scale_flag = true;
          p = nskip(p, 10);
          p = ncopy(cc, p, 10);
          if(sscanf(cc, "%f", &info->scale.matrix[offset]) != 1)
            scale_flag = false;
          p = ncopy(cc, p, 10);
          if(sscanf(cc, "%f", &info->scale.matrix[offset + 1]) != 1)
            scale_flag = false;
          p = ncopy(cc, p, 10);
          if(sscanf(cc, "%f", &info->scale.matrix[offset + 2]) != 1)
            scale_flag = false;
          p = nskip(p, 5);
          p = ncopy(cc, p, 10);
          if(sscanf(cc, "%f", &info->scale.matrix[offset + 3]) != 1)
            scale_flag = false;
          if(scale_flag)
            info->scale.flag[flag] = true;
          PRINTFB(G, FB_ObjectMolecule, FB_Blather)
            " PDBStrToCoordSet: SCALE%d %8.4f %8.4f %8.4f %8.4f\n", flag + 1,
            info->scale.matrix[offset],
            info->scale.matrix[offset + 1],
            info->scale.matrix[offset + 2], info->scale.matrix[offset + 3]
            ENDFB(G);
        }
        break;
      }
    } else if((p[0] == 'C') &&  /* CONECT */
              (p[1] == 'O') &&
              (p[2] == 'N') &&
              (p[3] == 'E') &&
              (p[4] == 'C') &&
              (p[5] == 'T') &&
              bondFlag && (!ignore_conect) && ((!*restart_model) || (!in_model))) {
      seen_conect = true;
      p = nskip(p, 6);
      p = ncopy(cc, p, 5);
      if(sscanf(cc, "%d", &b1) == 1)
        while(1) {
          p = ncopy(cc, p, 5);
          if(sscanf(cc, "%d", &b2) != 1)
            break;
          else {
            if((b1 >= 0) && (b2 >= 0) && (b1 != b2)) {  /* IDs must be positive and distinct */
              VLACheck(bond, BondType, nBond);
              if(b1 <= b2) {
                bond[nBond].index[0] = b1;      /* temporarily store the atom indexes */
                bond[nBond].index[1] = b2;
                bond[nBond].order = 1;
                bond[nBond].stereo = 0;
              } else {
                bond[nBond].index[0] = b2;
                bond[nBond].index[1] = b1;
                bond[nBond].order = 1;
                bond[nBond].stereo = 0;
              }
              nBond++;
            }
          }
        }
    } else if((p[0] == 'U') &&  /* USER */
              (p[1] == 'S') && (p[2] == 'E') && (p[3] == 'R') && (!*restart_model)) {
      /* Metaphorics key 'USER M4X ' was 'USER     ' */
      if((p[4] == ' ') &&
         (p[5] == 'M') && (p[6] == '4') && (p[7] == 'X') && (p[8] == ' ') && m4x) {

        int parsed_flag = false;

        p = nskip(p, 10);
        p = ntrim(cc, p, 6);

        /* is this a context name or a USER record? */
        switch (cc[0]) {
        case 'X':
          if(WordMatchExact(G, "XNAME", cc, true)) {    /* object name */
            p = nskip(p, 1);
            p = ntrim(m4x->xname, p, 10);
            if(m4x->xname[0]) {
              m4x->xname_flag = true;
              parsed_flag = true;
            }
          }
          break;
        case 'A':              /* alignment information */
          if(WordMatchExact(G, "ALIGN", cc, true)) {
            if(m4x->align && m4x->align->id_at_point) {
              M4XAlignType *align = m4x->align;
              char target[11];
              int atom_id, point_id;
              float fitness;
              p = nskip(p, 1);
              p = ncopy(cc, p, 6);
              if(sscanf(cc, "%d", &atom_id) == 1) {
                p = nskip(p, 1);
                p = ntrim(target, p, 10);
                if(target[0]) {
                  if(!align->target[0])
                    UtilNCopy(align->target, target, WordLength);
                  if(WordMatchExact(G, align->target, target, true)) {  /* must match the one target allowed */
                    p = nskip(p, 1);
                    p = ncopy(cc, p, 6);
                    if(sscanf(cc, "%d", &point_id) == 1) {
                      p = nskip(p, 1);
                      p = ncopy(cc, p, 6);
                      if(sscanf(cc, "%f", &fitness) == 1) {
                        VLACheck(align->id_at_point, int, point_id);
                        VLACheck(align->fitness, float, point_id);
                        if(point_id >= align->n_point)
                          align->n_point = point_id + 1;
                        align->id_at_point[point_id] = atom_id;
                        align->fitness[point_id] = fitness;
                        /*                        printf("read alignment atom %d to target %s point %d fitness %8.3f\n",
                           atom_id,target,point_id,fitness); */
                      }
                    }
                  }
                }
              }
            }
            parsed_flag = true;
          }
          break;
        }

        if((!parsed_flag) && m4x->context) {    /* named context of some sort... */
          int cn;
          int found = false;

          /* does context already exist ? */
          for(cn = 0; cn < m4x->n_context; cn++) {
            if(WordMatchExact(G, m4x->context[cn].name, cc, true)) {
              found = true;
              break;
            }
          }

          if(found) {

            M4XContextType *cont = m4x->context + cn;

            p = nskip(p, 1);
            p = ntrim(cc, p, 6);
            switch (cc[0]) {
            case 'B':
              if(WordMatchExact(G, "BORDER", cc, true) && bondFlag) {
                int order;

                p = nskip(p, 1);
                p = ncopy(cc, p, 6);
                if(sscanf(cc, "%d", &b1) == 1) {
                  p = nskip(p, 1);
                  p = ncopy(cc, p, 6);
                  if(sscanf(cc, "%d", &b2) == 1) {
                    p = nskip(p, 1);
                    p = ncopy(cc, p, 6);
                    if(sscanf(cc, "%d", &order) == 1) {
                      if((b1 >= 0) && (b2 >= 0)) {      /* IDs must be positive */
                        VLACheck(bond, BondType, nBond);
                        if(b1 <= b2) {
                          bond[nBond].index[0] = b1;    /* temporarily store the atom indexes */
                          bond[nBond].index[1] = b2;
                          bond[nBond].order = order;
                          bond[nBond].stereo = 0;

                        } else {
                          bond[nBond].index[0] = b2;
                          bond[nBond].index[1] = b1;
                          bond[nBond].order = order;
                          bond[nBond].stereo = 0;
                        }
                        nBond++;
                      }
                    }
                  }
                }
              }
              break;
            case 'S':
              if(WordMatchExact(G, "SITE", cc, true)) {
                if(cont->site) {
                  int id;
                  while(*cc) {
                    p = nskip(p, 1);
                    p = ncopy(cc, p, 6);
                    if(sscanf(cc, "%d", &id) == 1) {
                      VLACheck(cont->site, int, cont->n_site);
                      cont->site[cont->n_site++] = id;
                    }
                  }
                }
              }
              break;
            case 'L':
              if(WordMatchExact(G, "LIGAND", cc, true)) {
                if(cont->ligand) {
                  int id;
                  while(*cc) {
                    p = nskip(p, 1);
                    p = ncopy(cc, p, 6);
                    if(sscanf(cc, "%d", &id) == 1) {
                      VLACheck(cont->ligand, int, cont->n_ligand);
                      cont->ligand[cont->n_ligand++] = id;
                    }
                  }
                }
              }
              break;
            case 'W':
              if(WordMatchExact(G, "WATER", cc, true)) {
                if(cont->water) {
                  int id;
                  while(*cc) {
                    p = nskip(p, 1);
                    p = ncopy(cc, p, 6);
                    if(sscanf(cc, "%d", &id) == 1) {
                      VLACheck(cont->water, int, cont->n_water);
                      cont->water[cont->n_water++] = id;
                    }
                  }
                }
              }
              break;
            case 'H':
              if(WordMatchExact(G, "HBOND", cc, true)) {
                if(cont->hbond) {
                  int id1, id2;
                  float strength;
                  p = nskip(p, 1);
                  p = ncopy(cc, p, 6);
                  if(sscanf(cc, "%d", &id1) == 1) {
                    p = nskip(p, 1);
                    p = ncopy(cc, p, 6);
                    if(sscanf(cc, "%d", &id2) == 1) {
                      p = nskip(p, 1);
                      p = ncopy(cc, p, 6);
                      if(sscanf(cc, "%f", &strength) == 1) {
                        VLACheck(cont->hbond, M4XBondType, cont->n_hbond);
                        cont->hbond[cont->n_hbond].atom1 = id1;
                        cont->hbond[cont->n_hbond].atom2 = id2;
                        cont->hbond[cont->n_hbond].strength = strength;
                        cont->n_hbond++;
                      }
                    }
                  }
                }
              }
              break;
            case 'N':
              if(WordMatchExact(G, "NBOND", cc, true)) {
                if(cont->nbond) {
                  int id1, id2;
                  float strength;
                  p = nskip(p, 1);
                  p = ncopy(cc, p, 6);
                  if(sscanf(cc, "%d", &id1) == 1) {
                    p = nskip(p, 1);
                    p = ncopy(cc, p, 6);
                    if(sscanf(cc, "%d", &id2) == 1) {
                      p = nskip(p, 1);
                      p = ncopy(cc, p, 6);
                      if(sscanf(cc, "%f", &strength) == 1) {
                        VLACheck(cont->nbond, M4XBondType, cont->n_nbond);
                        cont->nbond[cont->n_nbond].atom1 = id1;
                        cont->nbond[cont->n_nbond].atom2 = id2;
                        cont->nbond[cont->n_nbond].strength = strength;
                        cont->n_nbond++;
                      }
                    }
                  }
                }
              }
              break;
            }
          }
        }
      }
    } else if((p[0] == 'A') &&  /* ANISOU */
              (p[1] == 'N') &&
              (p[2] == 'I') &&
              (p[3] == 'S') &&
              (p[4] == 'O') && (p[5] == 'U') && (!*restart_model) && (atomCount)) {
      ai = atInfo + atomCount - 1;

      /* TODO: check atom identifier match */

      {
        int dummy;
        p = nskip(p, 6);
        p = ncopy(cc, p, 5);
        if(!sscanf(cc, "%d", &dummy))
          dummy = 0;
        if(dummy == ai->id) {   /* ATOM ID must match */
          p = nskip(p, 17);
          {
            int dummy;
            p = ncopy(cc, p, 7);
            if(sscanf(cc, "%d", &dummy))
              ai->U11 = dummy / 10000.0F;
            p = ncopy(cc, p, 7);
            if(sscanf(cc, "%d", &dummy))
              ai->U22 = dummy / 10000.0F;
            p = ncopy(cc, p, 7);
            if(sscanf(cc, "%d", &dummy))
              ai->U33 = dummy / 10000.0F;
            p = ncopy(cc, p, 7);
            if(sscanf(cc, "%d", &dummy))
              ai->U12 = dummy / 10000.0F;
            p = ncopy(cc, p, 7);
            if(sscanf(cc, "%d", &dummy))
              ai->U13 = dummy / 10000.0F;
            p = ncopy(cc, p, 7);
            if(sscanf(cc, "%d", &dummy))
              ai->U23 = dummy / 10000.0F;
          }
        }
      }
    }

    /* END KEYWORDS */

    /* Secondary structure records */

    if(SSCode) {

      /* pretty confusing how this works... the following efficient (i.e. array-based)
         secondary structure lookup even when there are multiple insertion codes
         and where there may be multiple SS records for the residue using different 
         insertion codes */

      if(!ss[ss_chain1]) {      /* allocate new array for chain if necc. */
        ss[ss_chain1] = Calloc(int, cResvMask + 1);
      }

      sst = NULL;
      for(b = ss_resv1; b <= ss_resv2; b++) {   /* iterate over all residues indicated */
        index = b & cResvMask;
        if(ss[ss_chain1][index])
          sst = NULL;           /* make a unique copy in the event of multiple entries for one resv */
        if(!sst) {
          VLACheck(ss_list, SSEntry, n_ss);
          ssi = n_ss++;
          sst = ss_list + ssi;
          sst->resv1 = ss_resv1;
          sst->resv2 = ss_resv2;
          sst->chain1 = ss_chain1;
          sst->chain2 = ss_chain2;
          sst->type = SSCode;
          strcpy(sst->resi1, ss_resi1);
          strcpy(sst->resi2, ss_resi2);
          ss_found = true;
        }
        sst->next = ss[ss_chain1][index];
        ss[ss_chain1][index] = ssi;
        if(sst->next)
          sst = NULL;           /* force another unique copy */
      }

      if(ss_chain2 != ss_chain1) {      /* handle case where chains are not the same (undefined in PDB spec?) */
        if(!ss[ss_chain2]) {
          ss[ss_chain2] = Calloc(int, cResvMask + 1);
        }
        sst = NULL;
        for(b = ss_resv1; b <= ss_resv2; b++) { /* iterate over all residues indicated */
          index = b & cResvMask;
          if(ss[ss_chain2][index])
            sst = NULL;         /* make a unique copy in the event of multiple entries for one resv */
          if(!sst) {
            VLACheck(ss_list, SSEntry, n_ss);
            ssi = n_ss++;
            sst = ss_list + ssi;
            sst->resv1 = ss_resv1;
            sst->resv2 = ss_resv2;
            sst->chain1 = ss_chain1;
            sst->chain2 = ss_chain2;
            sst->type = SSCode;
            strcpy(sst->resi1, ss_resi1);
            strcpy(sst->resi2, ss_resi2);
          }
          sst->next = ss[ss_chain2][index];
          ss[ss_chain2][index] = ssi;
          if(sst->next)
            sst = NULL;         /* force another unique copy */
        }
      }
    }
    /* Atom records */

    if(AFlag && (!*restart_model)) {
      ai = atInfo + atomCount;
      p = nskip(p, 6);
      p = ncopy(cc, p, 5);
      if(!sscanf(cc, "%d", &ai->id))
        ai->id = 0;
      ai->rank = atomCount;

      p = nskip(p, 1);          /* to 12 */
      p = ncopy(literal_name, p, 4);
      if(literal_names) {
        strcpy(ai->name, literal_name);
      } else {
        ParseNTrim(ai->name, literal_name, 4);
      }

      p = ncopy(cc, p, 1);
      if(*cc == 32)
        ai->alt[0] = 0;
      else {
        ai->alt[0] = *cc;
        ai->alt[1] = 0;
      }

      p = ncopy(cc, p, 4);      /* now allowing for 4-letter residues */
      if(!sscanf(cc, "%s", ai->resn))
        ai->resn[0] = 0;
      else if(truncate_resn)    /* unless specifically disabled */
        ai->resn[3] = 0;

      if(ai->name[0]) {
        int name_len = strlen(ai->name);
        char name[4];
        switch (reformat_names) {
        case 1:                /* pdb compliant: HH12 becomes 2HH1, etc. */
          if(name_len > 3) {
            if((ai->name[0] >= 'A') && ((ai->name[0] <= 'Z')) &&
               (ai->name[3] >= '0') && (ai->name[3] <= '9')) {
              if(!(((ai->name[1] >= 'a') && (ai->name[1] <= 'z')) || ((ai->name[0] == 'C') && (ai->name[1] == 'L')) ||  /* try to be smart about */
                   ((ai->name[0] == 'B') && (ai->name[1] == 'R')) ||    /* distinguishing common atoms */
                   ((ai->name[0] == 'C') && (ai->name[1] == 'A')) ||    /* in all-caps from typical */
                   ((ai->name[0] == 'F') && (ai->name[1] == 'E')) ||    /* nonatomic abbreviations */
                   ((ai->name[0] == 'C') && (ai->name[1] == 'U')) ||
                   ((ai->name[0] == 'N') && (ai->name[1] == 'A')) ||
                   ((ai->name[0] == 'N') && (ai->name[1] == 'I')) ||
                   ((ai->name[0] == 'M') && (ai->name[1] == 'G')) ||
                   ((ai->name[0] == 'M') && (ai->name[1] == 'N')) ||
                   ((ai->name[0] == 'H') && (ai->name[1] == 'G')) ||
                   ((ai->name[0] == 'S') && (ai->name[1] == 'E')) ||
                   ((ai->name[0] == 'S') && (ai->name[1] == 'I')) ||
                   ((ai->name[0] == 'Z') && (ai->name[1] == 'N'))
                 )) {
                ctmp = ai->name[3];
                ai->name[3] = ai->name[2];
                ai->name[2] = ai->name[1];
                ai->name[1] = ai->name[0];
                ai->name[0] = ctmp;
              }
            }
          } else if(name_len == 3) {
            if((ai->name[0] == 'H') &&
               (ai->name[1] >= 'A') && ((ai->name[1] <= 'Z')) &&
               (ai->name[2] >= '0') && (ai->name[2] <= '9')) {
              AtomInfoGetPDB3LetHydroName(G, ai->resn, ai->name, name);
              if(name[0] == ' ')
                strcpy(ai->name, name + 1);
              else
                strcpy(ai->name, name);
            }
          }
          break;
        case 2:                /* amber compliant: 2HH1 becomes HH12 */
        case 3:                /* pdb compliant, but use IUPAC within PyMOL */
          if(ai->name[0]) {
            if((ai->name[0] >= '0') && (ai->name[0] <= '9') &&
               (!((ai->name[1] >= '0') && (ai->name[1] <= '9'))) && (ai->name[1] != 0)) {
              switch (strlen(ai->name)) {
              default:
                break;
              case 2:
                ctmp = ai->name[0];
                ai->name[0] = ai->name[1];
                ai->name[1] = ctmp;
                break;
              case 3:
                ctmp = ai->name[0];
                ai->name[0] = ai->name[1];
                ai->name[1] = ai->name[2];
                ai->name[2] = ctmp;
                break;
              case 4:
                ctmp = ai->name[0];
                ai->name[0] = ai->name[1];
                ai->name[1] = ai->name[2];
                ai->name[2] = ai->name[3];
                ai->name[3] = ctmp;
                break;
              }
              break;
        default:               /* AS IS */
              break;
            }
          }
          break;
        case 4:                /* simply read trim and write back out with 3-letter names starting from the
                                   second column, and four-letter names starting in the first */
          ncopy(cc, ai->name, 44);
          ParseNTrim(ai->name, cc, 4);
          break;
        }
      }

      p = ncopy(cc, p, 1);
      if(*cc == ' ')
        ai->chain[0] = 0;
      else {
        ai->chain[0] = *cc;
        ai->chain[1] = 0;
      }

      p = ncopy(cc, p, 5);      /* we treat insertion records as part of the residue identifier */
      if(!sscanf(cc, "%s", ai->resi))
        ai->resi[0] = 0;
      ai->resv = AtomResvFromResi(ai->resi);

      if(ssFlag) {              /* get secondary structure information (if avail) */

        ss_chain1 = ai->chain[0];
        index = ai->resv & cResvMask;
        if(ss[ss_chain1]) {
          ssi = ss[ss_chain1][index];
          while(ssi) {
            sst = ss_list + ssi;        /* contains shared entry, or unique linked list for each residue */
            /*                printf("%d<=%d<=%d, %s<=%s<=%s ", 
               sst->resv1,ai->resv,sst->resv2,
               sst->resi1,ai->resi,sst->resi2); */
            if(ai->resv >= sst->resv1)
              if(ai->resv <= sst->resv2)
                if((ai->resv != sst->resv1)
                   || (WordCompare(G, ai->resi, sst->resi1, true) >= 0))
                  if((ai->resv != sst->resv2)
                     || (WordCompare(G, ai->resi, sst->resi2, true) <= 0)) {
                    ai->ssType[0] = sst->type;
                    /*                          printf(" Y\n"); */
                    break;
                  }
            /*                printf(" N\n"); */
            ssi = sst->next;
          }
        }

      } else {
        ai->cartoon = cCartoon_tube;
      }
      if((!info) || (!info->is_pqr_file)) {     /* standard PDB file */

        p = nskip(p, 3);
        p = ncopy(cc, p, 8);
        sscanf(cc, "%f", coord + a);
        p = ncopy(cc, p, 8);
        sscanf(cc, "%f", coord + (a + 1));
        p = ncopy(cc, p, 8);
        sscanf(cc, "%f", coord + (a + 2));

        p = ncopy(cc, p, 6);
        if(!sscanf(cc, "%f", &ai->q))
          ai->q = 1.0;

        p = ncopy(cc, p, 6);
        if(!sscanf(cc, "%f", &ai->b))
          ai->b = 0.0;

        p = nskip(p, 6);
        p = ncopy(cc, p, 4);
        if(!ignore_pdb_segi) {
          if(!segi_override[0]) {
            if(!sscanf(cc, "%s", ai->segi))
              ai->segi[0] = 0;
            else {
              cc_saved = cc[3];
              ncopy(cc, p, 4);
              if((cc_saved == '1') &&   /* atom ID overflow? (nonstandard use...)... */
                 (cc[0] == '0') &&
                 (cc[1] == '0') && (cc[2] == '0') && (cc[3] == '0') && atomCount) {
                strcpy(segi_override, (ai - 1)->segi);
                strcpy(ai->segi, (ai - 1)->segi);
              }
            }
          } else {
            strcpy(ai->segi, segi_override);
          }
        } else {
          ai->segi[0] = 0;
        }

        p = ncopy(cc, p, 2);
        if(!sscanf(cc, "%s", ai->elem))
          ai->elem[0] = 0;
        else if(!((((ai->elem[0] >= 'a') && (ai->elem[0] <= 'z')) ||    /* don't get confused by PDB misuse */
                   ((ai->elem[0] >= 'A') && (ai->elem[0] <= 'Z'))) &&
                  (((ai->elem[1] == 0) ||
                    ((ai->elem[1] >= 'a') && (ai->elem[1] <= 'z')) ||
                    ((ai->elem[1] >= 'A') && (ai->elem[1] <= 'Z'))))))
          ai->elem[0] = 0;

        if(!ai->elem[0]) {
          if(((literal_name[0] == ' ') || ((literal_name[0] >= '0') && (literal_name[0] <= '9'))) && (literal_name[1] >= 'A') && (literal_name[1] <= 'Z')) {    /* infer element from name column */
            ai->elem[0] = literal_name[1];
            ai->elem[1] = 0;
          } else if(((literal_name[0] >= 'A') && (literal_name[0] <= 'Z')) && (((literal_name[1] >= 'A') && (literal_name[1] <= 'Z')) || ((literal_name[1] >= 'a') && (literal_name[1] <= 'z')))) {     /* infer element from name column */
            ai->elem[0] = literal_name[0];
            ai->elem[2] = 0;
            if((literal_name[1] >= 'A') && (literal_name[1] <= 'Z')) {  /* second letter is capitalized */
              if(bogus_name_alignment) {
                /* if other atom names aren't properly aligned */
                ai->elem[1] = 0;        /* kill 2nd letter */
              } else if(literal_name[0] == 'H') {
                /* or if this is an ultra-bogus PDB with inconsistent 
                   indendentation, and this is likely a hydrogen */
                ai->elem[1] = 0;        /* kill 2nd letter */
              } else {
                ai->elem[1] = tolower(literal_name[1]);
              }
            } else
              ai->elem[1] = literal_name[1];
          }
        }

        p = ncopy(cc, p, 2);
        if((cc[1] == '-') || (cc[1] == '+')) {
          /* only read formal charge when sign is present */
          char ctmp = cc[0];
          cc[0] = cc[1];
          cc[1] = ctmp;
          if(!sscanf(cc, "%d", &ai->formalCharge))
            ai->formalCharge = 0;
        }

        for(c = 0; c < cRepCnt; c++) {
          ai->visRep[c] = false;
        }

        /* end normal PDB */
      } else if(info && info->is_pqr_file) {
        /* PQR file format...not well defined, but basically PDB
           with charge and radius instead of B and Q.  Right now,
           we insist on PDB column format through the chain ID,
           and then switch over to whitespace delimited parsing 
           for the coordinates, charge, and radius */

        p = ParseWordNumberCopy(cc, p, MAXLINELEN - 1);
        sscanf(cc, "%f", coord + a);
        p = ParseWordNumberCopy(cc, p, MAXLINELEN - 1);
        sscanf(cc, "%f", coord + (a + 1));
        p = ParseWordNumberCopy(cc, p, MAXLINELEN - 1);
        sscanf(cc, "%f", coord + (a + 2));

        p = ParseWordNumberCopy(cc, p, MAXLINELEN - 1);
        if(!sscanf(cc, "%f", &ai->partialCharge))
          ai->partialCharge = 0.0F;

        p = ParseWordNumberCopy(cc, p, MAXLINELEN - 1);
        if(sscanf(cc, "%f", &ai->elec_radius) != 1)
          ai->elec_radius = 0.0F;
      }

      ai->visRep[cRepLine] = auto_show_lines;   /* show lines by default */
      ai->visRep[cRepNonbonded] = auto_show_nonbonded;  /* show lines by default */
      ai->visRep[cRepSphere] = auto_show_spheres;

      if(AFlag == 1)
        ai->hetatm = 0;
      else {
        ai->hetatm = 1;
        ai->flags = cAtomFlag_ignore;
      }

      AtomInfoAssignParameters(G, ai);
      AtomInfoAssignColors(G, ai);

      PRINTFD(G, FB_ObjectMolecule)
        "%s %s %s %s %8.3f %8.3f %8.3f %6.2f %6.2f %s\n",
        ai->name, ai->resn, ai->resi, ai->chain,
        *(coord + a), *(coord + a + 1), *(coord + a + 2), ai->b, ai->q, ai->segi ENDFD;

      if(atomCount < (nAtom - 1)) {     /* safety */
        a += 3;
        atomCount++;
      }
    }
    p = nextline(p);
  }

  /* END PASS 2 */

  if(bondFlag) {
    UtilSortInPlace(G, bond, nBond, sizeof(BondType), (UtilOrderFn *) BondInOrder);
    if(nBond) {
      if(!have_bond_order) {    /* handle PDB bond-order kludge */
        ii1 = bond;
        ii2 = bond + 1;
        nReal = 1;
        ii1->order = 1;
        a = nBond - 1;
        while(a) {
          if((ii1->index[0] == ii2->index[0]) && (ii1->index[1] == ii2->index[1])) {
            ii1->order++;       /* count dup */
          } else {
            ii1++;              /* non-dup, make copy */
            ii1->index[0] = ii2->index[0];
            ii1->index[1] = ii2->index[1];
            ii1->order = ii2->order;
            ii1->stereo = ii2->stereo;
            nReal++;
          }
          ii2++;
          a--;
        }
        nBond = nReal;
      }
      /* now, find atoms we're looking for */

      /* determine ranges */
      maxAt = nAtom;
      ii1 = bond;
      for(a = 0; a < nBond; a++) {
        if(ii1->index[0] > maxAt)
          maxAt = ii1->index[0];
        if(ii1->index[1] > maxAt)
          maxAt = ii1->index[1];
        ii1++;
      }
      for(a = 0; a < nAtom; a++)
        if(maxAt < atInfo[a].id)
          maxAt = atInfo[a].id;
      /* build index */
      maxAt++;
      idx = Alloc(int, maxAt + 1);
      for(a = 0; a < maxAt; a++) {
        idx[a] = -1;
      }
      for(a = 0; a < nAtom; a++)
        idx[atInfo[a].id] = a;
      /* convert indices to bonds */
      ii1 = bond;
      ii2 = bond;
      nReal = 0;
      {
        int unbond_cations = SettingGetGlobal_i(G, cSetting_pdb_unbond_cations);
        int flag;

        for(a = 0; a < nBond; a++) {

          if((ii1->index[0] >= 0) && ((ii1->index[1]) >= 0)) {
            if((idx[ii1->index[0]] >= 0) && (idx[ii1->index[1]] >= 0)) {        /* in case PDB file has bad bonds */
              ii2->index[0] = idx[ii1->index[0]];
              ii2->index[1] = idx[ii1->index[1]];
              ii2->order = ii1->order;
              if((ii2->index[0] >= 0) && (ii2->index[1] >= 0)) {

                if(!have_bond_order) {  /* handle PDB bond order kludge */
                  if(ii1->order <= 2)
                    ii2->order = 1;
                  else if(ii1->order <= 4)
                    ii2->order = 2;
                  else if(ii1->order <= 6)
                    ii2->order = 3;
                  else
                    ii2->order = 4;
                }
                flag = true;
                if(unbond_cations) {
                  if(AtomInfoIsFreeCation(G, atInfo + ii2->index[0]))
                    flag = false;
                  else if(AtomInfoIsFreeCation(G, atInfo + ii2->index[1]))
                    flag = false;
                }
                if(flag) {
                  atInfo[ii2->index[0]].bonded = true;
                  atInfo[ii2->index[1]].bonded = true;
                  nReal++;
                  ii2++;
                }
              }
            }
          }
          ii1++;
        }
      }
      nBond = nReal;
      FreeP(idx);
    }
  }
  if(ss_found && !quiet) {
    PRINTFB(G, FB_ObjectMolecule, FB_Details)
      " ObjectMolecule: Read secondary structure assignments.\n" ENDFB(G);
  }
  if(symmetry && !quiet && (!only_read_one_model)) {
    PRINTFB(G, FB_ObjectMolecule, FB_Details)
      " ObjectMolecule: Read crystal symmetry information.\n" ENDFB(G);
  }
  PRINTFB(G, FB_ObjectMolecule, FB_Blather)
    " PDBStr2CoordSet: Read %d bonds from CONECT records (%p).\n", nBond,
    (void *) bond ENDFB(G);
  cset = CoordSetNew(G);
  cset->NIndex = nAtom;
  cset->Coord = coord;
  cset->TmpBond = bond;
  cset->NTmpBond = nBond;
  if(symmetry)
    cset->Symmetry = symmetry;
  if(atInfoPtr)
    *atInfoPtr = atInfo;

  if((*restart_model) && (!foundNextModelFlag)) {
    /* if plan on need to reading another model into this object, 
       make sure there is another model to read... */
    p = *restart_model;
    while(*p) {
      if((p[0] == 'H') && (p[1] == 'E') && (p[2] == 'A') &&
         (p[3] == 'D') && (p[4] == 'E') && (p[5] == 'R')) {
        /* seeing HEADER means we're off the end of the existing file */
        break;
      }
      if((p[0] == 'M') && (p[1] == 'O') && (p[2] == 'D') &&
         (p[3] == 'E') && (p[4] == 'L') && (p[5] == ' ')) {
        foundNextModelFlag = true;
        break;
      }
      if((p[0] == 'E') && (p[1] == 'N') && (p[2] == 'D') &&
         (p[3] == 'M') && (p[4] == 'D') && (p[5] == 'L')) {
        foundNextModelFlag = true;
        break;
      }
      p = nextline(p);
    }
    if(!foundNextModelFlag) {
      *restart_model = NULL;
    }
  }

  if(ssFlag) {
    for(a = 0; a <= 255; a++) {
      FreeP(ss[a]);
    }
    VLAFreeP(ss_list);
  }

  if(!seen_model)
    *model_number = 1;

  if((*restart_model) && (*next_pdb)) { /* if we're mixing multistate objects and
                                           trajectories, then enforce sanity by 
                                           reading the models first... */
    if(is_end_of_object)
      (*restart_model) = NULL;
    else if((*next_pdb) < (*restart_model))
      (*next_pdb) = NULL;
  }
  return (cset);
}


/*========================================================================*/

void ObjectMoleculeM4XAnnotate(ObjectMolecule * I, M4XAnnoType * m4x, char *script_file,
                               int match_colors, int nbr_sele)
{
  int a;
  WordType name;
  M4XContextType *cont;

  if(m4x) {
    for(a = 0; a < m4x->n_context; a++) {
      cont = m4x->context + a;

      if(cont->site) {
        UtilNCopy(name, I->Obj.Name, sizeof(WordType));
        UtilNConcat(name, "_", sizeof(WordType));
        UtilNConcat(name, cont->name, sizeof(WordType));
        UtilNConcat(name, "_site", sizeof(WordType));
        SelectorSelectByID(I->Obj.G, name, I, cont->site, cont->n_site);
      }
      if(cont->ligand) {
        UtilNCopy(name, I->Obj.Name, sizeof(WordType));
        UtilNConcat(name, "_", sizeof(WordType));
        UtilNConcat(name, cont->name, sizeof(WordType));
        UtilNConcat(name, "_ligand", sizeof(WordType));
        SelectorSelectByID(I->Obj.G, name, I, cont->ligand, cont->n_ligand);
      }
      if(cont->water) {
        UtilNCopy(name, I->Obj.Name, sizeof(WordType));
        UtilNConcat(name, "_", sizeof(WordType));
        UtilNConcat(name, cont->name, sizeof(WordType));
        UtilNConcat(name, "_water", sizeof(WordType));
        SelectorSelectByID(I->Obj.G, name, I, cont->water, cont->n_water);
      }
      if(cont->hbond) {
        ObjectDist *distObj;
        UtilNCopy(name, I->Obj.Name, sizeof(WordType));
        UtilNConcat(name, "_", sizeof(WordType));
        UtilNConcat(name, cont->name, sizeof(WordType));
        UtilNConcat(name, "_hbond", sizeof(WordType));
        ExecutiveDelete(I->Obj.G, name);
        distObj = ObjectDistNewFromM4XBond(I->Obj.G, NULL,
                                           I, cont->hbond, cont->n_hbond, nbr_sele);
        if(match_colors)
          distObj->Obj.Color = I->Obj.Color;
        else
          distObj->Obj.Color = ColorGetIndex(I->Obj.G, "yellow");
        ObjectSetName((CObject *) distObj, name);
        if(distObj)
          ExecutiveManageObject(I->Obj.G, (CObject *) distObj, false, true);
      }

      if(cont->nbond && 0) {
        /*        ObjectDist *distObj; */
        UtilNCopy(name, I->Obj.Name, sizeof(WordType));
        UtilNConcat(name, "_", sizeof(WordType));
        UtilNConcat(name, cont->name, sizeof(WordType));
        UtilNConcat(name, "_nbond", sizeof(WordType));
        ExecutiveDelete(I->Obj.G, name);
        /*        distObj = ObjectDistNewFromM4XBond(I->Obj.G,NULL,
           I,
           cont->nbond,
           cont->n_nbond);
           if(distObj)
           ExecutiveManageObject(I->Obj.G,(CObject*)distObj,false,true); */

        {
          CGO *cgo = NULL;
          ObjectCGO *ocgo;

          cgo = CGONew(I->Obj.G);
          /*
             CGOBegin(cgo,GL_LINES);
             for(a=0;a<op1.nvv1;a++) {
             CGOVertexv(cgo,op2.vv1+(a*3));
             MatrixApplyTTTfn3f(1,v1,op2.ttt,op1.vv1+(a*3));
             CGOVertexv(cgo,v1);
           */
          CGOEnd(cgo);
          CGOStop(cgo);
          ocgo = ObjectCGOFromCGO(I->Obj.G, NULL, cgo, 0);
          if(match_colors)
            ocgo->Obj.Color = I->Obj.Color;
          else
            ocgo->Obj.Color = ColorGetIndex(I->Obj.G, "yellow");
          ObjectSetName((CObject *) ocgo, name);
          ExecutiveDelete(I->Obj.G, name);

          ExecutiveManageObject(I->Obj.G, (CObject *) ocgo, false, true);

          SceneInvalidate(I->Obj.G);
        }

      }

    }
    if(script_file)
      PParse(I->Obj.G, script_file);
  }
}

void ObjectMoleculeInitHBondCriteria(PyMOLGlobals * G, HBondCriteria * hbc)
{
  hbc->maxAngle = SettingGet_f(G, NULL, NULL, cSetting_h_bond_max_angle);
  hbc->maxDistAtMaxAngle = SettingGet_f(G, NULL, NULL, cSetting_h_bond_cutoff_edge);
  hbc->maxDistAtZero = SettingGet_f(G, NULL, NULL, cSetting_h_bond_cutoff_center);
  hbc->power_a = SettingGet_f(G, NULL, NULL, cSetting_h_bond_power_a);
  hbc->power_b = SettingGet_f(G, NULL, NULL, cSetting_h_bond_power_b);
  hbc->cone_dangle =
    (float) cos(PI * 0.5 * SettingGet_f(G, NULL, NULL, cSetting_h_bond_cone) / 180.0F);
  if(hbc->maxDistAtMaxAngle != 0.0F) {
    hbc->factor_a = (float) (0.5 / pow(hbc->maxAngle, hbc->power_a));
    hbc->factor_b = (float) (0.5 / pow(hbc->maxAngle, hbc->power_b));
  }
}

static int ObjectMoleculeTestHBond(float *donToAcc, float *donToH, float *hToAcc,
                                   float *accPlane, HBondCriteria * hbc)
{
  float nDonToAcc[3], nDonToH[3], nAccPlane[3], nHToAcc[3];
  double angle;
  double cutoff;
  double curve;
  double dist;
  float dangle;

/* A ~~ H - D */

  normalize23f(donToAcc, nDonToAcc);
  normalize23f(hToAcc, nHToAcc);
  if(accPlane) {                /* remember, plane need not exist if it's water... */
    normalize23f(accPlane, nAccPlane);
#if 0
    if(dot_product3f(nDonToAcc, nAccPlane) > (-hbc->cone_dangle))       /* don't allow D behind Acceptor plan */
      return 0;
#endif
    if(dot_product3f(nHToAcc, nAccPlane) > (-hbc->cone_dangle)) /* don't allow H behind Acceptor plane */
      return 0;
  }

  normalize23f(donToH, nDonToH);
  normalize23f(donToAcc, nDonToAcc);

  dangle = dot_product3f(nDonToH, nDonToAcc);
  if((dangle < 1.0F) && (dangle > 0.0F))
    angle = 180.0 * acos((double) dangle) / PI;
  else if(dangle > 0.0F)
    angle = 0.0;
  else
    angle = 90.0;

  if(angle > hbc->maxAngle)
    return 0;

  /* interpolate cutoff based on ADH angle */

  if(hbc->maxDistAtMaxAngle != 0.0F) {
    curve = (pow(angle, hbc->power_a) * hbc->factor_a +
             pow(angle, hbc->power_b) * hbc->factor_b);

    cutoff = (hbc->maxDistAtMaxAngle * curve) + (hbc->maxDistAtZero * (1.0 - curve));
  } else {
    cutoff = hbc->maxDistAtZero;
  }

  /*
     printf("angle %8.3f curve %8.3f %8.3f %8.3f %8.3f\n",angle,
     curve,cutoff,hbc->maxDistAtMaxAngle,hbc->maxDistAtZero);
   */

  dist = length3f(donToAcc);

  if(dist > cutoff)
    return 0;
  else
    return 1;

}


/*========================================================================*/

static int ObjectMoleculeFindBestDonorH(ObjectMolecule * I,
                                        int atom,
                                        int state, float *dir, float *best, int *is_real)
{
  int result = 0;
  CoordSet *cs;
  int n, nn;
  int idx;
  int a1;
  float cand[3], cand_dir[3];
  float best_dot = 0.0F, cand_dot;
  float *orig;

  ObjectMoleculeUpdateNeighbors(I);

  if((state >= 0) && (state < I->NCSet) && (cs = I->CSet[state]) && (atom < I->NAtom)) {

    if(I->DiscreteFlag) {
      if(cs == I->DiscreteCSet[atom]) {
        idx = I->DiscreteAtmToIdx[atom];
      } else {
        idx = -1;
      }
    } else {
      idx = cs->AtmToIdx[atom];
    }

    if(idx >= 0) {

      orig = cs->Coord + 3 * idx;

      /*  do we need to add any new hydrogens? */

      n = I->Neighbor[atom];
      nn = I->Neighbor[n++];

      /*      printf("nn %d valence %d %s\n",nn,
         I->AtomInfo[atom].valence,I->AtomInfo[atom].name);
         {
         int i;
         for(i=0;i<nn;i++) {
         printf("%d \n",I->Neighbor[n+2*i]);
         }
         }
       */

      if((nn < I->AtomInfo[atom].valence) || I->AtomInfo[atom].hb_donor) {      /* is there an implicit hydrogen? */
        if(ObjectMoleculeFindOpenValenceVector(I, state, atom, best, dir, -1)) {
          result = true;
          best_dot = dot_product3f(best, dir);
          add3f(orig, best, best);
          if(is_real)
            *is_real = false;
        }
      }
      /* iterate through real hydrogens looking for best match
         with desired direction */

      while(1) {                /* look for an attached non-hydrogen as a base */
        a1 = I->Neighbor[n];
        n += 2;
        if(a1 < 0)
          break;
        if(I->AtomInfo[a1].protons == 1) {      /* hydrogen */
          if(ObjectMoleculeGetAtomVertex(I, state, a1, cand)) { /* present */

            subtract3f(cand, orig, cand_dir);
            normalize3f(cand_dir);
            cand_dot = dot_product3f(cand_dir, dir);
            if(result) {        /* improved */
              if((best_dot < cand_dot) || ((is_real) && (!*is_real))) {
                best_dot = cand_dot;
                copy3f(cand, best);
                if(is_real)
                  *is_real = true;
              }
            } else {            /* first */
              result = true;
              copy3f(cand, best);
              best_dot = cand_dot;
              if(is_real)
                *is_real = true;
            }
          }
        }
      }
    }
  }
  return result;
}


/*========================================================================*/

int ObjectMoleculeGetCheckHBond(int *h_is_real,
                                float *h_crd_ret,
                                ObjectMolecule * don_obj,
                                int don_atom,
                                int don_state,
                                ObjectMolecule * acc_obj,
                                int acc_atom, int acc_state, HBondCriteria * hbc)
{
  int result = 0;
  CoordSet *csD, *csA;
  int idxD, idxA;
  float *vAcc, *vDon;
  float donToAcc[3];
  float donToH[3];
  float bestH[3];
  float hToAcc[3];
  float accPlane[3], *vAccPlane = NULL;

  /* first, check for existence of coordinate sets */

  if((don_state >= 0) &&
     (don_state < don_obj->NCSet) &&
     (csD = don_obj->CSet[don_state]) &&
     (acc_state >= 0) &&
     (acc_state < acc_obj->NCSet) &&
     (csA = acc_obj->CSet[acc_state]) &&
     (don_atom < don_obj->NAtom) && (acc_atom < acc_obj->NAtom)) {

    /* now check for coordinates of these actual atoms */

    if(don_obj->DiscreteFlag) {
      if(csD == don_obj->DiscreteCSet[don_atom]) {
        idxD = don_obj->DiscreteAtmToIdx[don_atom];
      } else {
        idxD = -1;
      }
    } else {
      idxD = csD->AtmToIdx[don_atom];
    }

    if(acc_obj->DiscreteFlag) {
      if(csA == acc_obj->DiscreteCSet[acc_atom]) {
        idxA = acc_obj->DiscreteAtmToIdx[acc_atom];
      } else {
        idxA = -1;
      }
    } else {
      idxA = csA->AtmToIdx[acc_atom];
    }

    if((idxA >= 0) && (idxD >= 0)) {

      /* now get local geometries, including 
         real or virtual hydrogen atom positions */

      vDon = csD->Coord + 3 * idxD;
      vAcc = csA->Coord + 3 * idxA;

      subtract3f(vAcc, vDon, donToAcc);

      if(ObjectMoleculeFindBestDonorH(don_obj,
                                      don_atom, don_state, donToAcc, bestH, h_is_real)) {

        subtract3f(bestH, vDon, donToH);
        subtract3f(vAcc, bestH, hToAcc);

        if(ObjectMoleculeGetAvgHBondVector(acc_obj, acc_atom,
                                           acc_state, accPlane, hToAcc) > 0.1) {
          vAccPlane = &accPlane[0];
        }
#if 0
        {
          float tmp[3];
          CGOLinewidth(DebugCGO, 4.0F);
          CGOBegin(DebugCGO, GL_LINES);

          CGOColor(DebugCGO, 0.0F, 1.0F, 0.0F); /* green */
          CGOVertexv(DebugCGO, vDon);
          normalize23f(donToAcc, tmp);
          add3f(vDon, tmp, tmp);
          CGOVertexv(DebugCGO, tmp);

          CGOColor(DebugCGO, 1.0F, 0.0F, 0.0F); /* red */
          CGOVertexv(DebugCGO, bestH);
          normalize23f(hToAcc, tmp);
          add3f(bestH, tmp, tmp);
          CGOVertexv(DebugCGO, tmp);

          CGOColor(DebugCGO, 0.0F, 1.0F, 1.0F); /* cyan */
          CGOVertexv(DebugCGO, vDon);
          normalize23f(donToH, tmp);
          add3f(vDon, tmp, tmp);
          CGOVertexv(DebugCGO, tmp);

          if(vAccPlane) {
            CGOColor(DebugCGO, 1.0F, 1.0F, 0.0F);       /* yellow */
            CGOVertexv(DebugCGO, vAcc);
            normalize23f(vAccPlane, tmp);
            add3f(vAcc, tmp, tmp);
            CGOVertexv(DebugCGO, tmp);
          }

          CGOEnd(DebugCGO);

        }
#endif
        result = ObjectMoleculeTestHBond(donToAcc, donToH, hToAcc, vAccPlane, hbc);
        if(result && h_crd_ret && h_is_real && *h_is_real)
          copy3f(bestH, h_crd_ret);
      }
    }
  }

  return (result);
}


/*========================================================================*/

float ObjectMoleculeGetMaxVDW(ObjectMolecule * I)
{
  float max_vdw = 0.0F;
  int a;
  AtomInfoType *ai;
  if(I->NAtom) {
    ai = I->AtomInfo;
    for(a = 0; a < I->NAtom; a++) {
      if(max_vdw < ai->vdw)
        max_vdw = ai->vdw;
      ai++;
    }
  }
  return (max_vdw);
}


/*========================================================================*/
#ifndef _PYMOL_NOPY
static PyObject *ObjectMoleculeCSetAsPyList(ObjectMolecule * I)
{
  PyObject *result = NULL;
  int a;
  result = PyList_New(I->NCSet);
  for(a = 0; a < I->NCSet; a++) {
    if(I->CSet[a]) {
      PyList_SetItem(result, a, CoordSetAsPyList(I->CSet[a]));
    } else {
      PyList_SetItem(result, a, PConvAutoNone(Py_None));
    }
  }
  return (PConvAutoNone(result));
}
#endif


/*static PyObject *ObjectMoleculeDiscreteCSetAsPyList(ObjectMolecule *I)
  {
  PyObject *result = NULL;
  return(PConvAutoNone(result));
  }*/

#ifndef _PYMOL_NOPY
static int ObjectMoleculeCSetFromPyList(ObjectMolecule * I, PyObject * list)
{
  int ok = true;
  int a;
  if(ok)
    ok = PyList_Check(list);
  if(ok) {
    VLACheck(I->CSet, CoordSet *, I->NCSet);
    for(a = 0; a < I->NCSet; a++) {
      if(ok)
        ok = CoordSetFromPyList(I->Obj.G, PyList_GetItem(list, a), &I->CSet[a]);
      PRINTFB(I->Obj.G, FB_ObjectMolecule, FB_Debugging)
        " ObjectMoleculeCSetFromPyList: ok %d after CoordSet %d\n", ok, a ENDFB(I->Obj.G);

      if(ok)
        if(I->CSet[a])          /* WLD 030205 */
          I->CSet[a]->Obj = I;
    }
  }
  return (ok);
}
#endif

#ifndef _PYMOL_NOPY
static PyObject *ObjectMoleculeBondAsPyList(ObjectMolecule * I)
{
  PyObject *result = NULL;
  PyObject *bond_list;
  BondType *bond;
  int a;

  result = PyList_New(I->NBond);
  bond = I->Bond;
  for(a = 0; a < I->NBond; a++) {
    bond_list = PyList_New(7);
    PyList_SetItem(bond_list, 0, PyInt_FromLong(bond->index[0]));
    PyList_SetItem(bond_list, 1, PyInt_FromLong(bond->index[1]));
    PyList_SetItem(bond_list, 2, PyInt_FromLong(bond->order));
    PyList_SetItem(bond_list, 3, PyInt_FromLong(bond->id));
    PyList_SetItem(bond_list, 4, PyInt_FromLong(bond->stereo));
    PyList_SetItem(bond_list, 5, PyInt_FromLong(bond->unique_id));
    PyList_SetItem(bond_list, 6, PyInt_FromLong(bond->has_setting));
    PyList_SetItem(result, a, bond_list);
    bond++;
  }

  return (PConvAutoNone(result));
}
#endif

#ifndef _PYMOL_NOPY
static int ObjectMoleculeBondFromPyList(ObjectMolecule * I, PyObject * list)
{
  int ok = true;
  int a;
  int stereo, ll = 0;
  PyObject *bond_list = NULL;
  BondType *bond;
  if(ok)
    ok = PyList_Check(list);
  if(ok)
    ok = ((I->Bond = VLAlloc(BondType, I->NBond)) != NULL);
  bond = I->Bond;
  for(a = 0; a < I->NBond; a++) {
    if(ok)
      bond_list = PyList_GetItem(list, a);
    if(ok)
      ok = PyList_Check(bond_list);
    if(ok)
      ll = PyList_Size(bond_list);
    if(ok)
      ok = PConvPyIntToInt(PyList_GetItem(bond_list, 0), &bond->index[0]);
    if(ok)
      ok = PConvPyIntToInt(PyList_GetItem(bond_list, 1), &bond->index[1]);
    if(ok)
      ok = PConvPyIntToInt(PyList_GetItem(bond_list, 2), &bond->order);
    if(ok)
      ok = PConvPyIntToInt(PyList_GetItem(bond_list, 3), &bond->id);
    if(ok)
      ok = PConvPyIntToInt(PyList_GetItem(bond_list, 4), &stereo);
    if(ok)
      bond->stereo = (short int) stereo;
    if(ok && (ll > 5)) {
      int has_setting;
      if(ok)
        ok = PConvPyIntToInt(PyList_GetItem(bond_list, 5), &bond->unique_id);
      if(ok)
        ok = PConvPyIntToInt(PyList_GetItem(bond_list, 6), &has_setting);
      if(ok)
        bond->has_setting = (short int) has_setting;
      if(ok && bond->unique_id) {       /* reserve existing IDs */
        bond->unique_id = SettingUniqueConvertOldSessionID(I->Obj.G, bond->unique_id);
        AtomInfoReserveUniqueID(I->Obj.G, bond->unique_id);
      }
    }
    bond++;
  }
  PRINTFB(I->Obj.G, FB_ObjectMolecule, FB_Debugging)
    " ObjectMoleculeBondFromPyList: ok %d after restore\n", ok ENDFB(I->Obj.G);

  return (ok);
}
#endif

#ifndef _PYMOL_NOPY
static PyObject *ObjectMoleculeAtomAsPyList(ObjectMolecule * I)
{
  PyObject *result = NULL;
  AtomInfoType *ai;
  int a;

  result = PyList_New(I->NAtom);
  ai = I->AtomInfo;
  for(a = 0; a < I->NAtom; a++) {
    PyList_SetItem(result, a, AtomInfoAsPyList(I->Obj.G, ai));
    ai++;
  }
  return (PConvAutoNone(result));
}
#endif

#ifndef _PYMOL_NOPY
static int ObjectMoleculeAtomFromPyList(ObjectMolecule * I, PyObject * list)
{
  int ok = true;
  int a;
  AtomInfoType *ai;
  if(ok)
    ok = PyList_Check(list);
  VLACheck(I->AtomInfo, AtomInfoType, I->NAtom + 1);
  ai = I->AtomInfo;
  for(a = 0; a < I->NAtom; a++) {
    if(ok)
      ok = AtomInfoFromPyList(I->Obj.G, ai, PyList_GetItem(list, a));
    ai++;
  }
  PRINTFB(I->Obj.G, FB_ObjectMolecule, FB_Debugging)
    " ObjectMoleculeAtomFromPyList: ok %d \n", ok ENDFB(I->Obj.G);
  return (ok);
}
#endif

int ObjectMoleculeNewFromPyList(PyMOLGlobals * G, PyObject * list,
                                ObjectMolecule ** result)
{
#ifdef _PYMOL_NOPY
  return 0;
#else
  int ok = true;
  ObjectMolecule *I = NULL;
  int discrete_flag;
  int ll;
  (*result) = NULL;

  if(ok)
    ok = PyList_Check(list);
  if(ok)
    ll = PyList_Size(list);
  /* TO SUPPORT BACKWARDS COMPATIBILITY...
     Always check ll when adding new PyList_GetItem's */
  if(ok)
    ok = PConvPyIntToInt(PyList_GetItem(list, 8), &discrete_flag);

  I = ObjectMoleculeNew(G, discrete_flag);
  if(ok)
    ok = (I != NULL);

  if(ok)
    ok = ObjectFromPyList(G, PyList_GetItem(list, 0), &I->Obj);
  if(ok)
    ok = PConvPyIntToInt(PyList_GetItem(list, 1), &I->NCSet);
  if(ok)
    ok = PConvPyIntToInt(PyList_GetItem(list, 2), &I->NBond);
  if(ok)
    ok = PConvPyIntToInt(PyList_GetItem(list, 3), &I->NAtom);
  if(ok)
    ok = ObjectMoleculeCSetFromPyList(I, PyList_GetItem(list, 4));
  if(ok)
    ok = CoordSetFromPyList(G, PyList_GetItem(list, 5), &I->CSTmpl);
  if(ok)
    ok = ObjectMoleculeBondFromPyList(I, PyList_GetItem(list, 6));
  if(ok)
    ok = ObjectMoleculeAtomFromPyList(I, PyList_GetItem(list, 7));
  if(ok)
    ok = PConvPyIntToInt(PyList_GetItem(list, 8), &I->DiscreteFlag);
  if(ok)
    ok = PConvPyIntToInt(PyList_GetItem(list, 9), &I->NDiscrete);
  if(ok)
    I->Symmetry = SymmetryNewFromPyList(G, PyList_GetItem(list, 10));
  if(ok)
    ok = PConvPyIntToInt(PyList_GetItem(list, 11), &I->CurCSet);
  if(ok)
    ok = PConvPyIntToInt(PyList_GetItem(list, 12), &I->BondCounter);
  if(ok)
    ok = PConvPyIntToInt(PyList_GetItem(list, 13), &I->AtomCounter);
  if(ok && I->DiscreteFlag) {
    int *dcs = NULL;
    int a, i;
    CoordSet *cs;
    VLACheck(I->DiscreteAtmToIdx, int, I->NDiscrete);
    VLACheck(I->DiscreteCSet, CoordSet *, I->NDiscrete);
    if(ok)
      ok =
        PConvPyListToIntArrayInPlace(PyList_GetItem(list, 14), I->DiscreteAtmToIdx,
                                     I->NDiscrete);
    if(ok)
      ok = PConvPyListToIntArray(PyList_GetItem(list, 15), &dcs);
    if(ok) {

      for(a = 0; a < I->NDiscrete; a++) {
        i = dcs[a];
        I->DiscreteCSet[a] = NULL;
        if(i >= 0) {
          if(i < I->NCSet) {
            cs = I->CSet[i];
            if(cs) {
              I->DiscreteCSet[a] = cs;
            }
          }
        }
      }
    }
    FreeP(dcs);
  }

  ObjectMoleculeInvalidate(I, cRepAll, cRepInvAll, -1);
  if(ok)
    (*result) = I;
  else {
    /* cleanup? */
  }

  return (ok);
#endif
}


/*========================================================================*/
PyObject *ObjectMoleculeAsPyList(ObjectMolecule * I)
{
#ifdef _PYMOL_NOPY
  return NULL;
#else
  PyObject *result = NULL;

  /* first, dump the atoms */

  result = PyList_New(16);
  PyList_SetItem(result, 0, ObjectAsPyList(&I->Obj));
  PyList_SetItem(result, 1, PyInt_FromLong(I->NCSet));
  PyList_SetItem(result, 2, PyInt_FromLong(I->NBond));
  PyList_SetItem(result, 3, PyInt_FromLong(I->NAtom));
  PyList_SetItem(result, 4, ObjectMoleculeCSetAsPyList(I));
  PyList_SetItem(result, 5, CoordSetAsPyList(I->CSTmpl));
  PyList_SetItem(result, 6, ObjectMoleculeBondAsPyList(I));
  PyList_SetItem(result, 7, ObjectMoleculeAtomAsPyList(I));
  PyList_SetItem(result, 8, PyInt_FromLong(I->DiscreteFlag));
  PyList_SetItem(result, 9, PyInt_FromLong(I->NDiscrete));
  PyList_SetItem(result, 10, SymmetryAsPyList(I->Symmetry));
  PyList_SetItem(result, 11, PyInt_FromLong(I->CurCSet));
  PyList_SetItem(result, 12, PyInt_FromLong(I->BondCounter));
  PyList_SetItem(result, 13, PyInt_FromLong(I->AtomCounter));

  if(I->DiscreteFlag) {
    int *dcs;
    int a;
    CoordSet *cs;

    /* get coordinate set indices */

    for(a = 0; a < I->NCSet; a++) {
      cs = I->CSet[a];
      if(cs) {
        cs->tmp_index = a;
      }
    }

    dcs = Alloc(int, I->NDiscrete);

    for(a = 0; a < I->NDiscrete; a++) {
      cs = I->DiscreteCSet[a];
      if(cs)
        dcs[a] = cs->tmp_index;
      else
        dcs[a] = -1;
    }

    PyList_SetItem(result, 14, PConvIntArrayToPyList(I->DiscreteAtmToIdx, I->NDiscrete));
    PyList_SetItem(result, 15, PConvIntArrayToPyList(dcs, I->NDiscrete));
    FreeP(dcs);
  } else {
    PyList_SetItem(result, 14, PConvAutoNone(NULL));
    PyList_SetItem(result, 15, PConvAutoNone(NULL));
  }

#if 0
  CObject Obj;
  struct CoordSet **CSet;
  int NCSet;
  struct CoordSet *CSTmpl;      /* template for trajectories, etc. */
  BondType *Bond;
  AtomInfoType *AtomInfo;
  int NAtom;
  int NBond;
  int DiscreteFlag, NDiscrete;
  int *DiscreteAtmToIdx;
  struct CoordSet **DiscreteCSet;
  int CurCSet;
  int SeleBase;                 /* for internal usage by  selector & only valid during selection process */
  CSymmetry *Symmetry;
  int *Neighbor;
  float *UndoCoord[cUndoMask + 1];
  int UndoState[cUndoMask + 1];
  int UndoNIndex[cUndoMask + 1];
  int UndoIter;
  CGO *UnitCellCGO;
  int BondCounter;
  int AtomCounter;
  struct CSculpt *Sculpt;
#endif

  return (PConvAutoNone(result));
#endif
}


/*========================================================================*/
int ObjectMoleculeConnect(ObjectMolecule * I, BondType ** bond, AtomInfoType * ai,
                          struct CoordSet *cs, int bondSearchMode,
                          int connectModeOverride)
{
#define cMULT 1
  PyMOLGlobals *G = I->Obj.G;
  int a, b, c, d, e, f, i, j;
  int a1, a2;
  float *v1, *v2, dst;
  int maxBond;
  MapType *map;
  int nBond;
  BondType *ii1, *ii2;
  int flag;
  int order;
  AtomInfoType *ai1, *ai2;
  float cutoff_s;
  float cutoff_h;
  float cutoff_v;
  float cutoff;
  float max_cutoff;
  int water_flag;
  int repeat = true;
  int discrete_chains = SettingGetGlobal_i(G, cSetting_pdb_discrete_chains);
  int connect_bonded = SettingGetGlobal_b(G, cSetting_connect_bonded);
  int connect_mode = SettingGetGlobal_i(G, cSetting_connect_mode);
  int unbond_cations = SettingGetGlobal_i(G, cSetting_pdb_unbond_cations);
  cutoff_v = SettingGet(G, cSetting_connect_cutoff);
  cutoff_s = cutoff_v + 0.2F;
  cutoff_h = cutoff_v - 0.2F;
  max_cutoff = cutoff_s;

  if(connectModeOverride >= 0)
    connect_mode = connectModeOverride;

  if(connect_mode == 2) {       /* force use of distance-based connectivity,
                                   ignoring that provided with file */
    bondSearchMode = true;
    cs->NTmpBond = 0;
    FreeP(cs->TmpBond);
  }

  /*  FeedbackMask[FB_ObjectMolecule]=0xFF; */
  nBond = 0;
  maxBond = cs->NIndex * 8;
  (*bond) = VLACalloc(BondType, maxBond);
  while(repeat) {
    repeat = false;

    if(cs->NIndex && bondSearchMode) {  /* &&(!I->DiscreteFlag) WLD 010527 */

      PRINTFB(G, FB_ObjectMolecule, FB_Blather)
        " ObjectMoleculeConnect: Searching for bonds amongst %d coordinates.\n",
        cs->NIndex ENDFB(G);
      if(Feedback(G, FB_ObjectMolecule, FB_Debugging)) {
        for(a = 0; a < cs->NIndex; a++)
          printf(" ObjectMoleculeConnect: coord %d %8.3f %8.3f %8.3f\n",
                 a, cs->Coord[a * 3], cs->Coord[a * 3 + 1], cs->Coord[a * 3 + 2]);
      }

      switch (connect_mode) {
      case 0:                  /* distance-based and explicit (not HETATM to HETATM) */
      case 3:                  /* distance-based and explicit (even HETATM to HETATM) */
      case 2:                  /* distance-based only */  {
          /* distance-based bond location  */
          int violations = 0;
          int *cnt = Alloc(int, cs->NIndex);
          int valcnt;

          for(i = 0; i < cs->NIndex; i++) {
            valcnt = AtomInfoGetExpectedValence(G, ai + cs->IdxToAtm[i]);
            if(valcnt < 0)
              valcnt = 6;
            cnt[i] = valcnt;
          }

          map = MapNew(G, max_cutoff + MAX_VDW, cs->Coord, cs->NIndex, NULL);
          if(map) {
            int dim12 = map->D1D2;
            int dim2 = map->Dim[2];

            for(i = 0; i < cs->NIndex; i++) {
              if(nBond > maxBond)
                break;
              v1 = cs->Coord + (3 * i);

              a1 = cs->IdxToAtm[i];
              ai1 = ai + a1;

              MapLocus(map, v1, &a, &b, &c);
              for(d = a - 1; d <= a + 1; d++) {
                int *j_ptr1 = map->Head + d * dim12 + (b - 1) * dim2;
                for(e = b - 1; e <= b + 1; e++) {
                  int *j_ptr2 = j_ptr1 + c - 1;
                  j_ptr1 += dim2;
                  for(f = c - 1; f <= c + 1; f++) {
                    j = *(j_ptr2++);    /*  *MapFirst(map,d,e,f)); */
                    while(j >= 0) {

                      if(i < j) {
                        v2 = cs->Coord + (3 * j);
                        dst = (float) diff3f(v1, v2);

                        a2 = cs->IdxToAtm[j];
                        ai2 = ai + a2;

                        dst -= ((ai1->vdw + ai2->vdw) / 2);

                        /* quick hack for water detection.  
                           they don't usually don't have CONECT records 
                           and may not be HETATMs though they are supposed to be... */

                        water_flag = false;
                        if(AtomInfoKnownWaterResName(G, ai1->resn))
                          water_flag = true;
                        else if(AtomInfoKnownWaterResName(G, ai2->resn))
                          water_flag = true;

                        /* workaround for hydrogens and sulfurs... */

                        if(ai1->hydrogen || ai2->hydrogen)
                          cutoff = cutoff_h;
                        else if(((ai1->elem[0] == 'S') && (!ai1->elem[1])) ||
                                ((ai2->elem[0] == 'S') && (!ai2->elem[1])))
                          cutoff = cutoff_s;
                        else
                          cutoff = cutoff_v;
                        if((dst <= cutoff) &&
                           (!(ai1->hydrogen && ai2->hydrogen)) &&
                           (water_flag || (!cs->TmpBond) ||
                            ((!(ai1->hetatm && ai2->hetatm) ||
                              (connect_mode == 3)))) &&
                           ((discrete_chains < 1) ||
                            ai1->chain[0] == ai2->chain[0]) &&
                           (connect_bonded || (!(ai1->bonded && ai2->bonded)))) {
                          flag = true;
                          if(water_flag)
                            if(!AtomInfoSameResidue(G, ai1, ai2))
                              flag = false;

                          if(flag) {
                            if(ai1->alt[0] != ai2->alt[0]) {    /* handle alternate conformers */
                              if(ai1->alt[0] && ai2->alt[0])
                                flag = false;   /* don't connect atoms with different, non-NULL
                                                   alternate conformations */
                            } else if(ai1->alt[0] && ai2->alt[0])
                              if(!AtomInfoSameResidue(G, ai1, ai2))
                                if(ai1->alt[0] != ai2->alt[0])
                                  flag = false; /* don't connect different, non-NULL 
                                                   alt conformations in 
                                                   different residues */
                          }

                          if(flag) {
                            if(ai1->alt[0] || ai2->alt[0])
                              if(water_flag)    /* hack to clean up water bonds */
                                if(!AtomInfoSameResidue(G, ai1, ai2))
                                  flag = false;
                          }

                          if(flag && unbond_cations) {
                            if(AtomInfoIsFreeCation(G, ai1))
                              flag = false;
                            else if(AtomInfoIsFreeCation(G, ai2))
                              flag = false;
                          }

                          if(flag) {
                            VLACheck((*bond), BondType, nBond);
                            (*bond)[nBond].index[0] = a1;
                            (*bond)[nBond].index[1] = a2;
                            (*bond)[nBond].stereo = 0;
                            order = 1;
                            if(discrete_chains < 0) {   /* if we're allowing bonds between chains,
                                                           then make sure things don't get out of hand */
                              if(cnt[i] == -1)
                                violations++;
                              if(cnt[j] == -1)
                                violations++;
                              cnt[i]--;
                              cnt[j]--;
                              if(violations > (cs->NIndex >> 3)) {
                                /* if more than 12% of the structure has excessive #'s of bonds... */
                                PRINTFB(G, FB_ObjectMolecule, FB_Blather)
                                  " ObjectMoleculeConnect: Assuming chains are discrete...\n"
                                  ENDFB(G);
                                discrete_chains = 1;
                                repeat = true;
                                goto do_it_again;
                              }
                            }
                            if(ai1->hetatm && (!ai1->resn[3])) {        /* common HETATMs we should know about... */
                              switch (ai1->resn[0]) {
                              case 'M':
                                switch (ai1->resn[1]) {
                                case 'S':
                                  switch (ai1->resn[2]) {
                                  case 'E':
                                    if(((!ai1->name[1]) && (!ai2->name[1])) &&
                                       (((ai1->name[0] == 'C') && (ai2->name[0] == 'O'))
                                        || ((ai1->name[0] == 'O')
                                            && (ai2->name[0] == 'C')))) {
                                      if(AtomInfoSameResidue(G, ai1, ai2)) {
                                        order = 2;
                                      }
                                    }
                                    break;
                                  }
                                  break;
                                }
                                break;
                              }
                            } else if((!ai1->hetatm)) {
                              if(AtomInfoSameResidue(I->Obj.G, ai1, ai2)) {
                                /* Standard disconnected PDB residue */
                                assign_pdb_known_residue(G, ai1, ai2, &order);
                              }
                            }
                            (*bond)[nBond].order = -order;      /* store tentative valence as negative */
                            nBond++;
                          }
                        }
                      }
                      j = MapNext(map, j);
                    }
                  }
                }
              }
            }
          do_it_again:
            MapFree(map);
            FreeP(cnt);
          }
        }
      case 1:                  /* only use explicit connectivity from file (don't do anything) */
        break;
      case 4:                  /* dictionary-based connectivity */
        /* TODO */
        break;
      }
      PRINTFB(G, FB_ObjectMolecule, FB_Blather)
        " ObjectMoleculeConnect: Found %d bonds.\n", nBond ENDFB(G);
      if(Feedback(G, FB_ObjectMolecule, FB_Debugging)) {
        for(a = 0; a < nBond; a++)
          printf(" ObjectMoleculeConnect: bond %d ind0 %d ind1 %d\n",
                 a, (*bond)[a].index[0], (*bond)[a].index[1]);
      }
    }
    if(repeat) {
      nBond = 0;
    }
  }
  if(cs->NTmpBond && cs->TmpBond) {
    int check_conect_all = false;
    int pdb_conect_all = false;
    PRINTFB(G, FB_ObjectMolecule, FB_Blather)
      " ObjectMoleculeConnect: incorporating explicit bonds. %d %d\n",
      nBond, cs->NTmpBond ENDFB(G);
    if((nBond == 0) && (cs->NTmpBond > 0) &&
       bondSearchMode && (connect_mode == 0) && cs->NIndex) {
      /* if we were no bonds were found, and we have explicit connectivity,
       * try to determine if we need to set pdb_conect_mode */
      for(i = 0; i < cs->NIndex; i++) {
        a1 = cs->IdxToAtm[i];
        ai1 = ai + a1;
        if(ai1->bonded && (!ai1->hetatm)) {
          /* apparent PDB ATOM record with explicit bonding... */
          check_conect_all = true;
          break;
        }
      }
    }

    VLACheck((*bond), BondType, (nBond + cs->NTmpBond));
    ii1 = (*bond) + nBond;
    ii2 = cs->TmpBond;
    {
      register int n_atom = I->NAtom;
      for(a = 0; a < cs->NTmpBond; a++) {
        a1 = cs->IdxToAtm[ii2->index[0]];       /* convert bonds from index space */
        a2 = cs->IdxToAtm[ii2->index[1]];       /* to atom space */
        if((a1 >= 0) && (a2 >= 0) && (a1 < n_atom) && (a2 < n_atom)) {
          if(check_conect_all) {
            if((!ai[a1].hetatm) && (!ai[a2].hetatm)) {
              /* found bond between non-HETATMs -- so tell PyMOL to CONECT all ATOMs
               * when writing out a PDB file */
              pdb_conect_all = true;
            }
          }
          ai[a1].bonded = true;
          ai[a2].bonded = true;
          (*ii1) = (*ii2);      /* note this copies owned ids and thus settings etc. */
          ii1->index[0] = a1;
          ii1->index[1] = a2;
          ii1++;
          ii2++;

        }
      }
    }

    nBond = nBond + cs->NTmpBond;
    VLAFreeP(cs->TmpBond);
    cs->NTmpBond = 0;

    if(pdb_conect_all) {
      int dummy;
      if(!SettingGetIfDefined_b(G, I->Obj.Setting, cSetting_pdb_conect_all, &dummy)) {
        CSetting **handle = NULL;
        if(I->Obj.fGetSettingHandle) {
          handle = I->Obj.fGetSettingHandle(&I->Obj, -1);
          if(handle) {
            SettingCheckHandle(G, handle);
            SettingSet_b(*handle, cSetting_pdb_conect_all, true);
          }
        }
      }
    }
  }

  if(cs->NTmpLinkBond && cs->TmpLinkBond) {
    PRINTFB(G, FB_ObjectMolecule, FB_Blather)
      "ObjectMoleculeConnect: incorporating linkage bonds. %d %d\n",
      nBond, cs->NTmpLinkBond ENDFB(G);
    VLACheck((*bond), BondType, (nBond + cs->NTmpLinkBond));
    ii1 = (*bond) + nBond;
    ii2 = cs->TmpLinkBond;
    for(a = 0; a < cs->NTmpLinkBond; a++) {
      a1 = ii2->index[0];       /* first atom is in object */
      a2 = cs->IdxToAtm[ii2->index[1]]; /* second is in the cset */
      ai[a1].bonded = true;
      ai[a2].bonded = true;
      (*ii1) = (*ii2);          /* note this copies owned ids and thus settings etc. */
      ii1->index[0] = a1;
      ii1->index[1] = a2;
      ii1++;
      ii2++;
    }
    nBond = nBond + cs->NTmpLinkBond;
    VLAFreeP(cs->TmpLinkBond);
    cs->NTmpLinkBond = 0;
  }

  PRINTFD(G, FB_ObjectMolecule)
    " ObjectMoleculeConnect: elminating duplicates with %d bonds...\n", nBond ENDFD;

  if(!I->DiscreteFlag) {
    UtilSortInPlace(G, (*bond), nBond, sizeof(BondType), (UtilOrderFn *) BondInOrder);
    if(nBond) {                 /* eliminate duplicates */
      ii1 = (*bond) + 1;
      ii2 = (*bond) + 1;
      a = nBond - 1;
      nBond = 1;
      if(a > 0)
        while(a--) {
          if((ii2->index[0] != (ii1 - 1)->index[0]) ||
             (ii2->index[1] != (ii1 - 1)->index[1])) {
            *(ii1++) = *(ii2++);        /* copy bond */
            nBond++;
          } else {
            if((ii2->order > 0) && ((ii1 - 1)->order < 0))
              (ii1 - 1)->order = ii2->order;    /* use most certain valence */
            ii2++;              /* skip bond */
          }
        }
      VLASize((*bond), BondType, nBond);
    }
  }
  /* restore bond oder positivity */

  ii1 = *bond;
  for(a = 0; a < nBond; a++) {
    if(ii1->order < 0)
      ii1->order = -ii1->order;
    ii1++;
  }

  PRINTFD(G, FB_ObjectMolecule)
    " ObjectMoleculeConnect: leaving with %d bonds...\n", nBond ENDFD;
  return (nBond);
}


/*========================================================================*/
void ObjectMoleculeSort(ObjectMolecule * I)
{                               /* sorts atoms and bonds */
  int *index;
  int *outdex = NULL;
  register int a, b;
  CoordSet *cs, **dcs;
  AtomInfoType *atInfo;
  int *dAtmToIdx;
  if(!I->DiscreteFlag) {        /* currently, discrete objects are never sorted */
    int n_bytes = sizeof(int) * I->NAtom;
    int already_in_order = true;
    int i_NAtom = I->NAtom;
    index = AtomInfoGetSortedIndex(I->Obj.G, &I->Obj, I->AtomInfo, i_NAtom, &outdex);
    for(a = 0; a < i_NAtom; a++) {
      if(index[a] != a) {
        already_in_order = false;
        break;
      }
    }
    if(!already_in_order) {     /* if we aren't already in perfect order */

      for(a = 0; a < I->NBond; a++) {   /* bonds */
        I->Bond[a].index[0] = outdex[I->Bond[a].index[0]];
        I->Bond[a].index[1] = outdex[I->Bond[a].index[1]];
      }

      for(a = -1; a < I->NCSet; a++) {  /* coordinate set mapping */
        if(a < 0) {
          cs = I->CSTmpl;
        } else {
          cs = I->CSet[a];
        }

        if(cs) {
          register int cs_NIndex = cs->NIndex;
          register int *cs_IdxToAtm = cs->IdxToAtm;
          register int *cs_AtmToIdx = cs->AtmToIdx;
          for(b = 0; b < cs_NIndex; b++)
            cs_IdxToAtm[b] = outdex[cs_IdxToAtm[b]];
          if(cs_AtmToIdx) {
            memset(cs_AtmToIdx, -1, n_bytes);
            /*          for(b=0;b<i_NAtom;b++)
               cs_AtmToIdx[b]=-1; */
            for(b = 0; b < cs_NIndex; b++) {
              cs_AtmToIdx[cs_IdxToAtm[b]] = b;
            }
          }
        }
      }

      atInfo = (AtomInfoType *) VLAMalloc(i_NAtom, sizeof(AtomInfoType), 5, true);
      /* autozero here is important */
      for(a = 0; a < i_NAtom; a++)
        atInfo[a] = I->AtomInfo[index[a]];
      VLAFreeP(I->AtomInfo);
      I->AtomInfo = atInfo;

      if(I->DiscreteFlag) {
        dcs = VLAlloc(CoordSet *, i_NAtom);
        dAtmToIdx = VLAlloc(int, i_NAtom);
        for(a = 0; a < i_NAtom; a++) {
          b = index[a];
          dcs[a] = I->DiscreteCSet[b];
          dAtmToIdx[a] = I->DiscreteAtmToIdx[b];
        }
        VLAFreeP(I->DiscreteCSet);
        VLAFreeP(I->DiscreteAtmToIdx);
        I->DiscreteCSet = dcs;
        I->DiscreteAtmToIdx = dAtmToIdx;
      }
    }
    AtomInfoFreeSortedIndexes(I->Obj.G, index, outdex);

    UtilSortInPlace(I->Obj.G, I->Bond, I->NBond, sizeof(BondType),
                    (UtilOrderFn *) BondInOrder);
    /* sort...important! */
    ObjectMoleculeInvalidate(I, cRepAll, cRepInvAtoms, -1);     /* important */

  }
}

Generated by  Doxygen 1.6.0   Back to index