C++ PRECONDITION函数代码示例

int
TListViewCtrl::GetTopIndex() const
{
  PRECONDITION(GetHandle());
  return static_cast<int>(SendMessage(LVM_GETTOPINDEX));
}

      // ------------------------------------------------------------------------
      //
      //
      //
      // ------------------------------------------------------------------------
      void addAngles(const ROMol &mol,const AtomicParamVect &params,
                     ForceFields::ForceField *field){
        PRECONDITION(mol.getNumAtoms()==params.size(),"bad parameters");
        PRECONDITION(field,"bad forcefield");
        ROMol::ADJ_ITER nbr1Idx;
        ROMol::ADJ_ITER end1Nbrs;
        ROMol::ADJ_ITER nbr2Idx;
        ROMol::ADJ_ITER end2Nbrs;
        RingInfo *rings=mol.getRingInfo();

        unsigned int nAtoms=mol.getNumAtoms();
        for(unsigned int j=0;j<nAtoms;j++){
          if(!params[j]) continue;
          const Atom *atomJ=mol.getAtomWithIdx(j);
          if(atomJ->getDegree()==1) continue;
          boost::tie(nbr1Idx,end1Nbrs)=mol.getAtomNeighbors(atomJ);
          for (;nbr1Idx!=end1Nbrs;nbr1Idx++) {
            const Atom *atomI=mol[*nbr1Idx].get();
            unsigned int i=atomI->getIdx();
            if(!params[i]) continue;
            boost::tie(nbr2Idx,end2Nbrs)=mol.getAtomNeighbors(atomJ);
            for (;nbr2Idx!=end2Nbrs;nbr2Idx++) {
              if (nbr2Idx<(nbr1Idx+1)) {
                continue;
              }
              const Atom *atomK=mol[*nbr2Idx].get();
              unsigned int k=atomK->getIdx();
              if(!params[k]) continue;
              // skip special cases:
              if( !(atomJ->getHybridization()==Atom::SP3D && atomJ->getDegree()==5) ){
                const Bond *b1 =mol.getBondBetweenAtoms(i,j);
                const Bond *b2 =mol.getBondBetweenAtoms(k,j);
                // FIX: recognize amide bonds here.
                AngleBendContrib *contrib;
                int order=0;
                switch(atomJ->getHybridization()){
                case Atom::SP:
                  order=1;
                  break;
                // the following is a hack to get decent geometries
                // with 3- and 4-membered rings incorporating sp2 atoms
                case Atom::SP2:
                  order=3;
                  // if the central atom is in a ring of size 3
                  if (rings->isAtomInRingOfSize(j, 3)) {
                    // if the central atom and one of the bonded atoms, but not the
                    //  other one are inside a ring, then this angle is between a
                    // ring substituent and a ring edge
                    if ((rings->isAtomInRingOfSize(i, 3) && !rings->isAtomInRingOfSize(k, 3))
                      || (!rings->isAtomInRingOfSize(i, 3) && rings->isAtomInRingOfSize(k, 3))) {
                      order = 30;
                    }
                    // if all atoms are inside the ring, then this is one of ring angles
                    else if (rings->isAtomInRingOfSize(i, 3) && rings->isAtomInRingOfSize(k, 3)) {
                      order = 35;
                    }
                  }
                  // if the central atom is in a ring of size 4
                  else if (rings->isAtomInRingOfSize(j, 4)) {
                    // if the central atom and one of the bonded atoms, but not the
                    //  other one are inside a ring, then this angle is between a
                    // ring substituent and a ring edge
                    if ((rings->isAtomInRingOfSize(i, 4) && !rings->isAtomInRingOfSize(k, 4))
                      || (!rings->isAtomInRingOfSize(i, 4) && rings->isAtomInRingOfSize(k, 4))) {
                      order = 40;
                    }
                    // if all atoms are inside the ring, then this is one of ring angles
                    else if (rings->isAtomInRingOfSize(i, 4) && rings->isAtomInRingOfSize(k, 4)) {
                      order = 45;
                    }
                  }
                  break;
                case Atom::SP3D2:
                  order=4;
                  break;
                default:
                  order=0;
                  break;
                } 
                  
                contrib = new AngleBendContrib(field,i,j,k,
                                               b1->getBondTypeAsDouble(),
                                               b2->getBondTypeAsDouble(),
                                               params[i],params[j],params[k],order);
                field->contribs().push_back(ForceFields::ContribPtr(contrib));
              }
            }
          }
        }
      }

      // ------------------------------------------------------------------------
      //
      //
      //
      // ------------------------------------------------------------------------
      void addTorsions(const ROMol &mol,const AtomicParamVect &params,
                       ForceFields::ForceField *field,
                       std::string torsionBondSmarts){
        PRECONDITION(mol.getNumAtoms()==params.size(),"bad parameters");
        PRECONDITION(field,"bad forcefield");

        // find all of the torsion bonds:
        std::vector<MatchVectType> matchVect;
        ROMol *query=SmartsToMol(torsionBondSmarts);
        TEST_ASSERT(query);
        unsigned int nHits=SubstructMatch(mol,*query,matchVect);
        delete query;

        for(unsigned int i=0; i<nHits; i++){
          MatchVectType match=matchVect[i];
          TEST_ASSERT(match.size()==2);
          int idx1=match[0].second;
          int idx2=match[1].second;
          if(!params[idx1]||!params[idx2]) continue;
          const Bond *bond=mol.getBondBetweenAtoms(idx1,idx2);
          std::vector<TorsionAngleContrib *> contribsHere;
          TEST_ASSERT(bond);
          const Atom *atom1=mol.getAtomWithIdx(idx1);
          const Atom *atom2=mol.getAtomWithIdx(idx2);

          if( (atom1->getHybridization()==Atom::SP2||atom1->getHybridization()==Atom::SP3) &&
              (atom2->getHybridization()==Atom::SP2||atom2->getHybridization()==Atom::SP3) ){
            ROMol::OEDGE_ITER beg1,end1;
            boost::tie(beg1,end1) = mol.getAtomBonds(atom1);
            while(beg1!=end1){
              const Bond *tBond1=mol[*beg1].get();
              if(tBond1!=bond){
                int bIdx = tBond1->getOtherAtomIdx(idx1);
                ROMol::OEDGE_ITER beg2,end2;
                boost::tie(beg2,end2) = mol.getAtomBonds(atom2);
                while(beg2 != end2){
                  const Bond *tBond2=mol[*beg2].get();
                  if(tBond2!=bond && tBond2!=tBond1){
                    int eIdx=tBond2->getOtherAtomIdx(idx2);
                    // make sure this isn't a three-membered ring:
                    if(eIdx != bIdx){
                      // we now have a torsion involving atoms (bonds):
                      //  bIdx - (tBond1) - idx1 - (bond) - idx2 - (tBond2) - eIdx
                      TorsionAngleContrib *contrib;

                      // if either of the end atoms is SP2 hybridized, set a flag
                      // here.  
                      bool hasSP2=false;
                      if(mol.getAtomWithIdx(bIdx)->getHybridization()==Atom::SP2 ||
                         mol.getAtomWithIdx(bIdx)->getHybridization()==Atom::SP2) {
                        hasSP2 = true;
                      }
                      //std::cout << "Torsion: " << bIdx << "-" << idx1 << "-" << idx2 << "-" << eIdx << std::endl;
                      //if(okToIncludeTorsion(mol,bond,bIdx,idx1,idx2,eIdx)){
                        //std::cout << "  INCLUDED" << std::endl;
                        contrib = new TorsionAngleContrib(field,bIdx,idx1,idx2,eIdx,
                                                          bond->getBondTypeAsDouble(),
                                                          atom1->getAtomicNum(),
                                                          atom2->getAtomicNum(),
                                                          atom1->getHybridization(),
                                                          atom2->getHybridization(),
                                                          params[idx1],params[idx2],
                                                          hasSP2);
                        field->contribs().push_back(ForceFields::ContribPtr(contrib));
                        contribsHere.push_back(contrib);
                      //}
                    }
                  }
                  beg2++;
                }
              }
              beg1++;
            }
          }
          // now divide the force constant for each contribution to the torsion energy
          // about this bond by the number of contribs about this bond:
          for(std::vector<TorsionAngleContrib *>::iterator chI=contribsHere.begin();
              chI!=contribsHere.end();++chI){
            (*chI)->scaleForceConstant(contribsHere.size());
          }
        }

      }

double calcBondRestLength(const MMFFBond *mmffBondParams)
{
    PRECONDITION(mmffBondParams, "bond parameters not found");

    return mmffBondParams->r0;
}

unsigned int Atom::getDegree() const {
  PRECONDITION(dp_mol,"degree not defined for atoms not associated with molecules");
  return getOwningMol().getAtomDegree(this);
}

//
/// Informs the sheet that information in a sheet has changed. The sheet enables the
/// 'Apply' button.
//
void
TPropertySheet::PageChanged(const TPropertyPage& pg)
{
  PRECONDITION(HPROPSHEETPAGE(pg));
  SendMessage(PSM_CHANGED, TParam1(pg.GetHandle()));
}

void OMFile::removeAllDefaultEncodings(void)
{
  TRACE("OMFile::removeAllDefaultEncodings()");
  PRECONDITION( "Valid default encoding map", _defaultEncodings );
  _defaultEncodings->clear();
}

  MatchVectType findFuncGroupsOnMol(const ROMol &mol, 
				    const FragCatParams *params,
				    INT_VECT &fgBonds) {
    PRECONDITION(params,"bad params");

    fgBonds.clear();
    
    std::pair<int, int> amat;
    MatchVectType aidFgrps;
    std::vector<MatchVectType> fgpMatches;
    std::vector<MatchVectType>::const_iterator mati;
    MatchVectType::const_iterator mi;
    int aid;
    //const ROMol *fgrp;

    INT_VECT_CI bi;
    aidFgrps.clear();
    
    int fid = 0;
    const MOL_SPTR_VECT &fgrps = params->getFuncGroups();
    MOL_SPTR_VECT::const_iterator fgci;
    
    for (fgci = fgrps.begin(); fgci != fgrps.end(); fgci++) {
      const ROMol *fgrp = fgci->get();
      std::string fname;
      (*fgci)->getProp(common_properties::_Name, fname);
      //std::cout << "Groups number: " << fname << "\n";
      //(*fgci)->debugMol(std::cout);
      //mol->debugMol(std::cout);
      // match this functional group onto the molecule
      SubstructMatch(mol, *fgrp, fgpMatches);

      // loop over all the matches we get for this fgroup
      for (mati = fgpMatches.begin(); mati != fgpMatches.end(); mati++) {
	//FIX: we will assume that the first atom in fgrp is always the connection
	// atom
	amat = mati->front();
	aid = amat.second; //FIX: is this correct - the second entry in the pair is the atom ID from mol

	// grab the list of atom Ids from mol that match the functional group 
	INT_VECT bondIds, maids;
	for (mi = mati->begin(); mi != mati->end(); mi++) {
	  maids.push_back(mi->second);
	}

	// create a list of bond IDs from these atom ID 
	// these are the bond in mol that are part of portion that matches the 
	// functional group
	bondIds = Subgraphs::bondListFromAtomList(mol, maids);
	
	// now check if all these bonds have been covered as part of larger 
	// functional group that was dealt with earlier
	// FIX: obviously we assume here that the function groups in params 
	// come in decreasing order of size.
	bool allDone = true;
	for (bi = bondIds.begin(); bi != bondIds.end(); bi++) {
	  if (std::find(fgBonds.begin(), fgBonds.end(), (*bi)) == fgBonds.end()) {
	    allDone = false;
	    fgBonds.push_back(*bi);
	  }
	}
	
	if (!allDone) {
	  // this functional group mapping onto mol is not part of a larger func
	  // group mapping so record it
	  aidFgrps.push_back(std::pair<int, int>(aid, fid));
	}
      }
      fid++;

    }

    
    return aidFgrps;
  }

void constructFragmenterBondTypes(
    std::istream *inStream,
    const std::map<unsigned int, std::string> &atomTypes,
    std::vector<FragmenterBondType> &defs, const std::string &comment,
    bool validate, bool labelByConnector) {
  PRECONDITION(inStream, "no stream");
  defs.clear();
  defs.resize(0);

  unsigned int line = 0;
  while (!inStream->eof()) {
    ++line;
    std::string tempStr = getLine(inStream);
    if (tempStr == "" || tempStr.find(comment) == 0) continue;
    std::vector<std::string> tokens;
    boost::split(tokens, tempStr, boost::is_any_of(" \t"),
                 boost::token_compress_on);
    if (tokens.size() < 3) {
      BOOST_LOG(rdWarningLog) << "line " << line << " is too short"
                              << std::endl;
      continue;
    }
    unsigned int idx1 = boost::lexical_cast<unsigned int>(tokens[0]);
    if (atomTypes.find(idx1) == atomTypes.end()) {
      BOOST_LOG(rdWarningLog) << "atom type #" << idx1 << " not recognized."
                              << std::endl;
      continue;
    }
    unsigned int idx2 = boost::lexical_cast<unsigned int>(tokens[1]);
    if (atomTypes.find(idx2) == atomTypes.end()) {
      BOOST_LOG(rdWarningLog) << "atom type #" << idx2 << " not recognized."
                              << std::endl;
      continue;
    }
    std::string sma1 = atomTypes.find(idx1)->second;
    std::string sma2 = atomTypes.find(idx2)->second;
    std::string smarts = "[$(" + sma1 + ")]" + tokens[2] + "[$(" + sma2 + ")]";
    ROMol *p = SmartsToMol(smarts);
    if (validate) {
      if (!p) {
        BOOST_LOG(rdWarningLog) << "cannot convert SMARTS " << smarts
                                << " to molecule at line " << line << std::endl;
        continue;
      }
    }
    FragmenterBondType fbt;
    fbt.atom1Type = idx1;
    fbt.atom2Type = idx2;
    if (labelByConnector) {
      fbt.atom1Label = idx1;
      fbt.atom2Label = idx2;
    } else {
      fbt.atom1Label = idx2;
      fbt.atom2Label = idx1;
    }
    if (p) {
      // for the purposes of replacing the bond, we'll use just the first
      // character to set the bond type (if we recognize it):
      switch (tokens[2][0]) {
        case '-':
          fbt.bondType = Bond::SINGLE;
          break;
        case '=':
          fbt.bondType = Bond::DOUBLE;
          break;
        case '#':
          fbt.bondType = Bond::TRIPLE;
          break;
        case ':':
          fbt.bondType = Bond::AROMATIC;
          break;
        default:
          fbt.bondType = p->getBondWithIdx(0)->getBondType();
      }
      fbt.query = ROMOL_SPTR(p);
    } else {
      fbt.bondType = Bond::UNSPECIFIED;
      fbt.query = ROMOL_SPTR();
    }
    defs.push_back(fbt);
  }
}

//
/// Returns the class name of the underlying control associated with
/// the TColumnHeader object.
/// \note The logic used depends on the availability of native
///       Common Control support. In the case where OWL provides
///       the underlying support, we'll specify a TColumnHeader -specific
///       classname although that's not necessary [it eases debugging]
//
TWindow::TGetClassNameReturnType
TColumnHeader::GetClassName()
{
  PRECONDITION(TCommCtrl::IsAvailable());
  return WC_HEADER;
}

void LeaveRawTokenMode (void)
/* Leave raw token mode. */
{
    PRECONDITION (RawMode > 0);
    --RawMode;
}

double calcBondForceConstant(const MMFFBond *mmffBondParams)
{
    PRECONDITION(mmffBondParams, "bond parameters not found");

    return mmffBondParams->kb;
}

std::string GetBondSmiles(const Bond *bond, int atomToLeftIdx, bool doKekule,
                          bool allBondsExplicit) {
  PRECONDITION(bond, "bad bond");
  if (atomToLeftIdx < 0) atomToLeftIdx = bond->getBeginAtomIdx();

  std::string res = "";
  bool aromatic = false;
  if (!doKekule && (bond->getBondType() == Bond::SINGLE ||
                    bond->getBondType() == Bond::DOUBLE ||
                    bond->getBondType() == Bond::AROMATIC)) {
    Atom *a1, *a2;
    a1 = bond->getOwningMol().getAtomWithIdx(atomToLeftIdx);
    a2 = bond->getOwningMol().getAtomWithIdx(
        bond->getOtherAtomIdx(atomToLeftIdx));
    if ((a1->getIsAromatic() && a2->getIsAromatic()) &&
        (a1->getAtomicNum() || a2->getAtomicNum()))
      aromatic = true;
  }

  Bond::BondDir dir = bond->getBondDir();

  if (bond->hasProp(common_properties::_TraversalRingClosureBond)) {
    // std::cerr<<"FLIP: "<<bond->getIdx()<<"
    // "<<bond->getBeginAtomIdx()<<"-"<<bond->getEndAtomIdx()<<std::endl;
    // if(dir==Bond::ENDDOWNRIGHT) dir=Bond::ENDUPRIGHT;
    // else if(dir==Bond::ENDUPRIGHT) dir=Bond::ENDDOWNRIGHT;
    bond->clearProp(common_properties::_TraversalRingClosureBond);
  }

  switch (bond->getBondType()) {
    case Bond::SINGLE:
      if (dir != Bond::NONE && dir != Bond::UNKNOWN) {
        switch (dir) {
          case Bond::ENDDOWNRIGHT:
            if (bond->getOwningMol().hasProp(common_properties::_doIsoSmiles))
              res = "\\";
            break;
          case Bond::ENDUPRIGHT:
            if (bond->getOwningMol().hasProp(common_properties::_doIsoSmiles))
              res = "/";
            break;
          default:
            break;
        }
      } else {
        // if the bond is marked as aromatic and the two atoms
        //  are aromatic, we need no marker (this arises in kekulized
        //  molecules).
        // FIX: we should be able to dump kekulized smiles
        //   currently this is possible by removing all
        //   isAromatic flags, but there should maybe be another way
        if (allBondsExplicit)
          res = "-";
        else if (aromatic && !bond->getIsAromatic())
          res = "-";
      }
      break;
    case Bond::DOUBLE:
      // see note above
      if (!aromatic || !bond->getIsAromatic()) res = "=";
      break;
    case Bond::TRIPLE:
      res = "#";
      break;
    case Bond::AROMATIC:
      if (dir != Bond::NONE && dir != Bond::UNKNOWN) {
        switch (dir) {
          case Bond::ENDDOWNRIGHT:
            if (bond->getOwningMol().hasProp(common_properties::_doIsoSmiles))
              res = "\\";
            break;
          case Bond::ENDUPRIGHT:
            if (bond->getOwningMol().hasProp(common_properties::_doIsoSmiles))
              res = "/";
            break;
          default:
            break;
        }
      } else if (allBondsExplicit || !aromatic) {
        res = ":";
      }
      break;
    case Bond::DATIVE:
      if (atomToLeftIdx >= 0 &&
          bond->getBeginAtomIdx() == static_cast<unsigned int>(atomToLeftIdx))
        res = ">";
      else
        res = "<";
      break;
    default:
      res = "~";
  }
  return res;
}

std::string FragmentSmilesConstruct(
    ROMol &mol, int atomIdx, std::vector<Canon::AtomColors> &colors,
    const UINT_VECT &ranks, bool doKekule, bool canonical,
    bool doIsomericSmiles, bool allBondsExplicit, bool allHsExplicit,
    std::vector<unsigned int> &atomOrdering,
    const boost::dynamic_bitset<> *bondsInPlay = 0,
    const std::vector<std::string> *atomSymbols = 0,
    const std::vector<std::string> *bondSymbols = 0) {
  PRECONDITION(!bondsInPlay || bondsInPlay->size() >= mol.getNumBonds(),
               "bad bondsInPlay");
  PRECONDITION(!atomSymbols || atomSymbols->size() >= mol.getNumAtoms(),
               "bad atomSymbols");
  PRECONDITION(!bondSymbols || bondSymbols->size() >= mol.getNumBonds(),
               "bad bondSymbols");

  Canon::MolStack molStack;
  // try to prevent excessive reallocation
  molStack.reserve(mol.getNumAtoms() + mol.getNumBonds());
  std::stringstream res;

  std::map<int, int> ringClosureMap;
  int ringIdx, closureVal;
  if (!canonical) mol.setProp(common_properties::_StereochemDone, 1);
  std::list<unsigned int> ringClosuresToErase;

  Canon::canonicalizeFragment(mol, atomIdx, colors, ranks, molStack,
                              bondsInPlay, bondSymbols, doIsomericSmiles);
  Bond *bond = 0;
  BOOST_FOREACH (Canon::MolStackElem mSE, molStack) {
    switch (mSE.type) {
      case Canon::MOL_STACK_ATOM:
        if (!ringClosuresToErase.empty()) {
          BOOST_FOREACH (unsigned int rclosure, ringClosuresToErase) {
            ringClosureMap.erase(rclosure);
          }
          ringClosuresToErase.clear();
        }
        // std::cout<<"\t\tAtom: "<<mSE.obj.atom->getIdx()<<std::endl;
        if (!atomSymbols) {
          res << GetAtomSmiles(mSE.obj.atom, doKekule, bond, allHsExplicit);
        } else {
          res << (*atomSymbols)[mSE.obj.atom->getIdx()];
        }
        atomOrdering.push_back(mSE.obj.atom->getIdx());

        break;
      case Canon::MOL_STACK_BOND:
        bond = mSE.obj.bond;
        // std::cout<<"\t\tBond: "<<bond->getIdx()<<std::endl;
        if (!bondSymbols) {
          res << GetBondSmiles(bond, mSE.number, doKekule, allBondsExplicit);
        } else {
          res << (*bondSymbols)[bond->getIdx()];
        }
        break;
      case Canon::MOL_STACK_RING:
        ringIdx = mSE.number;
        // std::cout<<"\t\tRing: "<<ringIdx;
        if (ringClosureMap.count(ringIdx)) {
          // the index is already in the map ->
          //   we're closing a ring, so grab
          //   the index and then delete the value:
          closureVal = ringClosureMap[ringIdx];
          // ringClosureMap.erase(ringIdx);
          ringClosuresToErase.push_back(ringIdx);
        } else {
          // we're opening a new ring, find the index for it:
          closureVal = 1;
          bool done = false;
          // EFF: there's got to be a more efficient way to do this
          while (!done) {
            std::map<int, int>::iterator mapIt;
            for (mapIt = ringClosureMap.begin(); mapIt != ringClosureMap.end();
                 mapIt++) {
              if (mapIt->second == closureVal) break;
            }
            if (mapIt == ringClosureMap.end()) {
              done = true;
            } else {
              closureVal += 1;
            }
          }
          ringClosureMap[ringIdx] = closureVal;
        }
        if (closureVal >= 10) {
          res << "%";
        }
        // std::cerr << " > " << closureVal <<std::endl;
        res << closureVal;
        break;
      case Canon::MOL_STACK_BRANCH_OPEN:
        res << "(";
        break;
      case Canon::MOL_STACK_BRANCH_CLOSE:
        res << ")";
        break;
      default:
        break;
    }

      // ------------------------------------------------------------------------
      //
      //
      //
      // ------------------------------------------------------------------------
      void addTrigonalBipyramidAngles(const Atom *atom,const ROMol &mol, int confId,
                                      const AtomicParamVect &params,
                                      ForceFields::ForceField *field){
        PRECONDITION(atom,"bad atom");
        PRECONDITION(atom->getHybridization()==Atom::SP3D,"bad hybridization");
        PRECONDITION(atom->getDegree()==5,"bad degree");
        PRECONDITION(mol.getNumAtoms()==params.size(),"bad parameters");
        PRECONDITION(field,"bad forcefield");

        const Bond *ax1=0,*ax2=0;
        const Bond *eq1=0,*eq2=0,*eq3=0;

        const Conformer &conf = mol.getConformer(confId);
        //------------------------------------------------------------
        // identify the axial and equatorial bonds:
        double mostNeg=100.0;
        ROMol::OEDGE_ITER beg1,end1;
        boost::tie(beg1,end1) = mol.getAtomBonds(atom);
        unsigned int aid = atom->getIdx();
        while(beg1!=end1){
          const Bond *bond1=mol[*beg1].get();
          unsigned int oaid = bond1->getOtherAtomIdx(aid);
          RDGeom::Point3D v1=conf.getAtomPos(aid).directionVector(conf.getAtomPos(oaid));
                  
          ROMol::OEDGE_ITER beg2,end2;
          boost::tie(beg2,end2) = mol.getAtomBonds(atom);
          while(beg2 != end2){
            const Bond *bond2=mol[*beg2].get();
            if(bond2->getIdx() > bond1->getIdx()){
              unsigned int oaid2 = bond2->getOtherAtomIdx(aid);
              RDGeom::Point3D v2=conf.getAtomPos(aid).directionVector(conf.getAtomPos(oaid2));
              double dot=v1.dotProduct(v2);
              if(dot<mostNeg){
                mostNeg = dot;
                ax1 = bond1;
                ax2 = bond2;
              }
            }
            ++beg2;
          }
          ++beg1;
        }
        CHECK_INVARIANT(ax1,"axial bond not found");
        CHECK_INVARIANT(ax2,"axial bond not found");
	
        boost::tie(beg1,end1) = mol.getAtomBonds(atom);
        while(beg1!=end1){
          const Bond *bond=mol[*beg1].get();
	  ++beg1;
	  if(bond==ax1 || bond==ax2) continue;
	  if(!eq1) eq1=bond;
	  else if(!eq2) eq2=bond;
	  else if(!eq3) eq3=bond;
	}

        CHECK_INVARIANT(eq1,"equatorial bond not found");
        CHECK_INVARIANT(eq2,"equatorial bond not found");
        CHECK_INVARIANT(eq3,"equatorial bond not found");

        
        //------------------------------------------------------------
        // alright, add the angles:
        AngleBendContrib *contrib;
        int atomIdx=atom->getIdx();
        int i,j;

        // Axial-Axial
        i=ax1->getOtherAtomIdx(atomIdx);
        j=ax2->getOtherAtomIdx(atomIdx);
        if(params[i]&&params[j]){
          contrib = new AngleBendContrib(field,i,atomIdx,j,
                                         ax1->getBondTypeAsDouble(),
                                         ax2->getBondTypeAsDouble(),
                                         params[i],params[atomIdx],params[j],2);
          field->contribs().push_back(ForceFields::ContribPtr(contrib));
        }        
        // Equatorial-Equatorial
        i=eq1->getOtherAtomIdx(atomIdx);
        j=eq2->getOtherAtomIdx(atomIdx);
        if(params[i]&&params[j]){
          contrib = new AngleBendContrib(field,i,atomIdx,j,
                                         eq1->getBondTypeAsDouble(),
                                         eq2->getBondTypeAsDouble(),
                                         params[i],params[atomIdx],params[j],3);
          field->contribs().push_back(ForceFields::ContribPtr(contrib));
        }
        i=eq1->getOtherAtomIdx(atomIdx);
        j=eq3->getOtherAtomIdx(atomIdx);
        if(params[i]&&params[j]){
          contrib = new AngleBendContrib(field,i,atomIdx,j,
                                         eq1->getBondTypeAsDouble(),
                                         eq3->getBondTypeAsDouble(),
                                         params[i],params[atomIdx],params[j],3);
          field->contribs().push_back(ForceFields::ContribPtr(contrib));
        }
//.........这里部分代码省略.........