C++ Molecule::bond方法代码示例

TEST_F(MoleculeTest, removeBond)
{
  Molecule molecule;
  Atom a = molecule.addAtom(1);
  Atom b = molecule.addAtom(1);
  Bond bondAB = molecule.addBond(a, b);
  Atom c = molecule.addAtom(1);
  molecule.addBond(b, c, 2);

  EXPECT_EQ(3, molecule.atomCount());
  EXPECT_EQ(2, molecule.bondCount());
  EXPECT_TRUE(molecule.bond(a, b).isValid());
  EXPECT_TRUE(molecule.bond(b, c).isValid());

  molecule.removeBond(bondAB);

  EXPECT_EQ(3, molecule.atomCount());
  EXPECT_EQ(1, molecule.bondCount());
  EXPECT_FALSE(molecule.bond(a, b).isValid());
  EXPECT_TRUE(molecule.bond(b, c).isValid());

  molecule.clearBonds();

  EXPECT_EQ(0, molecule.bondCount());
}

void BallAndStick::process(const Molecule &molecule,
                           Rendering::GroupNode &node)
{
  // Add a sphere node to contain all of the spheres.
  GeometryNode *geometry = new GeometryNode;
  node.addChild(geometry);
  SphereGeometry *spheres = new SphereGeometry;
  spheres->identifier().molecule = &molecule;
  spheres->identifier().type = Rendering::AtomType;
  geometry->addDrawable(spheres);

  for (Index i = 0; i < molecule.atomCount(); ++i) {
    Core::Atom atom = molecule.atom(i);
    unsigned char atomicNumber = atom.atomicNumber();
    const unsigned char *c = Elements::color(atomicNumber);
    Vector3ub color(c[0], c[1], c[2]);
    spheres->addSphere(atom.position3d().cast<float>(), color,
                       static_cast<float>(Elements::radiusVDW(atomicNumber))
                       * 0.3f);
  }

  float bondRadius = 0.1f;
  CylinderGeometry *cylinders = new CylinderGeometry;
  cylinders->identifier().molecule = &molecule;
  cylinders->identifier().type = Rendering::BondType;
  geometry->addDrawable(cylinders);
  for (Index i = 0; i < molecule.bondCount(); ++i) {
    Core::Bond bond = molecule.bond(i);
    Vector3f pos1 = bond.atom1().position3d().cast<float>();
    Vector3f pos2 = bond.atom2().position3d().cast<float>();
    Vector3ub color1(Elements::color(bond.atom1().atomicNumber()));
    Vector3ub color2(Elements::color(bond.atom2().atomicNumber()));
    Vector3f bondVector = pos2 - pos1;
    float bondLength = bondVector.norm();
    bondVector /= bondLength;
    switch (bond.order()) {
    case 3: {
      Vector3f delta = bondVector.unitOrthogonal() * (2.0f * bondRadius);
      cylinders->addCylinder(pos1 + delta, bondVector, bondLength, bondRadius,
                             color1, color2, i);
      cylinders->addCylinder(pos1 - delta, bondVector, bondLength, bondRadius,
                             color1, color2, i);
    }
    default:
    case 1:
      cylinders->addCylinder(pos1, bondVector, bondLength, bondRadius,
                             color1, color2, i);
      break;
    case 2: {
      Vector3f delta = bondVector.unitOrthogonal() * bondRadius;
      cylinders->addCylinder(pos1 + delta, bondVector, bondLength, bondRadius,
                             color1, color2, i);
      cylinders->addCylinder(pos1 - delta, bondVector, bondLength, bondRadius,
                             color1, color2, i);
    }
    }
  }
}

TEST(CjsonTest, bonds)
{
  CjsonFormat cjson;
  Molecule molecule;
  bool success = cjson.readFile(std::string(AVOGADRO_DATA) +
                                "/data/ethane.cjson", molecule);
  EXPECT_TRUE(success);
  EXPECT_EQ(cjson.error(), "");
  EXPECT_EQ(molecule.data("name").toString(), "Ethane");
  EXPECT_EQ(molecule.atomCount(), static_cast<size_t>(8));
  EXPECT_EQ(molecule.bondCount(), static_cast<size_t>(7));

  Bond bond = molecule.bond(0);
  EXPECT_EQ(bond.atom1().index(), static_cast<size_t>(0));
  EXPECT_EQ(bond.atom2().index(), static_cast<size_t>(1));
  EXPECT_EQ(bond.order(), static_cast<unsigned char>(1));
  bond = molecule.bond(6);
  EXPECT_EQ(bond.atom1().index(), static_cast<size_t>(4));
  EXPECT_EQ(bond.atom2().index(), static_cast<size_t>(7));
  EXPECT_EQ(bond.order(), static_cast<unsigned char>(1));
}

TEST(RWMoleculeTest, MoleculeToRWMolecule)
{
  Molecule mol;
  typedef Molecule::AtomType Atom;
  typedef Molecule::BondType Bond;
  Atom a0 = mol.addAtom(1);
  Atom a1 = mol.addAtom(6);
  Atom a2 = mol.addAtom(9);
  Bond b0 = mol.addBond(a0, a2);
  a1.setPosition3d(Vector3(0, 6, 9));
  b0.setOrder(3);

  RWMolecule rwmol(mol, 0);
  EXPECT_EQ(rwmol.atomCount(), mol.atomCount());
  EXPECT_EQ(rwmol.bondCount(), mol.bondCount());
  EXPECT_EQ(rwmol.atom(2).atomicNumber(), mol.atom(2).atomicNumber());
  EXPECT_EQ(rwmol.bond(0).order(), mol.bond(0).order());
}

  void BondDelegate::initialize()
  {
    Molecule *molecule = m_widget->molecule();

    // remove any existing rows
    if (m_label->childCount())
      model()->removeRows(m_label, 0, m_label->childCount());

    // add the bonds...  
    model()->insertRows(m_label, 0, molecule->numBonds());
    for (int i = 0; i < m_label->childCount(); ++i) {
      ProjectTreeItem *item = m_label->child(i);
      item->setData(0, tr("bond %1").arg(i));
      // set the primitive
      PrimitiveList primitives;
      primitives.append(molecule->bond(i));
      item->setPrimitives(primitives);
    }

  }

   /* 
   * Generate, attempt and accept or reject a Monte Carlo move.
   *
   * A complete move involve the following stepts:
   *
   *    1) choose a molecule-i randomly;
   *    2) find molecule-j of the same type as i, and monomer pairs satisfying
   *       the distance criterion exist;
   *    3) Delete the two rebridging monomers; first molecule-i, then -j.
   *    4) Rebuild molecule-j, then -i. (in the JCP paper, the order is chosen
   *       at random.)
   *
   */
   bool CfbDoubleRebridgeMove::move() 
   {
      bool      found, accept;
      int       iMol, jMol, sign, beginId, nEnd, i;
      Molecule *iPtr;     // pointer to i-th molecule
      Molecule *jPtr;     // pointer to j-th molecule
      int      *bonds;    // bond types of a consequtive blocks of atoms
      Atom     *thisPtr;  // pointer to the current end atom
      Atom     *tempPtr;  // pointer used for swap atom positions
      double    rosenbluth, rosen_r,  rosen_f;
      double    energy, energy_r, energy_f;
      double    po2n, pn2o;
      Vector    swapPos;

      incrementNAttempt();

      // Choose a random direction
      if (random().uniform(0.0, 1.0) > 0.5)
         sign = +1;
      else
         sign = -1;

      // Search for a pair of candidate sites and calculate probability 
      found = forwardScan(sign, iMol, jMol, beginId, po2n);
      if (!found) return found;

      iPtr = &(system().molecule(speciesId_, iMol));
      jPtr = &(system().molecule(speciesId_, jMol));

      // Save positions for atoms to be regrown
      thisPtr = &(iPtr->atom(beginId));
      tempPtr = &(jPtr->atom(beginId));
      for (i = 0; i < nRegrow_; ++i) {
         iOldPos_[i] = thisPtr->position();
         jOldPos_[i] = tempPtr->position();

         thisPtr += sign;
         tempPtr += sign;
      }

      // Prepare for the rebridge move
      accept = false;

      bonds = new int[nRegrow_ + 1];
      if (sign == +1) {
         for (i = 0; i < nRegrow_ + 1; ++i)
            bonds[i] = iPtr->bond(beginId - 1 + nRegrow_ - i).typeId();
      } else {
         for (i = 0; i < nRegrow_ + 1; ++i)
            bonds[i] = iPtr->bond(beginId - nRegrow_ + i).typeId();
      }

      // Deleting atoms: first i-th molecule, then j-th
      rosen_r = 1.0;
      energy_r = 0.0;

      thisPtr = &(iPtr->atom(beginId + (nRegrow_-1)*sign));
      deleteSequence(nRegrow_, sign, thisPtr, bonds, rosenbluth, energy);
      rosen_r *= rosenbluth;
      energy_r += energy;

      thisPtr = &(jPtr->atom(beginId + (nRegrow_-1)*sign));
      deleteSequence(nRegrow_, sign, thisPtr, bonds, rosenbluth, energy);
      rosen_r *= rosenbluth;
      energy_r += energy;

      // Swap the atom positions of the dangling end
      if (sign == +1)
         nEnd = iPtr->nAtom() - beginId - nRegrow_;
      else
         nEnd = beginId + 1 - nRegrow_;

      thisPtr = &(iPtr->atom(beginId + nRegrow_*sign));
      tempPtr = &(jPtr->atom(beginId + nRegrow_*sign));
      for (i = 0; i < nEnd; ++i) {
         swapPos = thisPtr->position();
         //system().moveAtom(*thisPtr, tempPtr->position());
         thisPtr->position() = tempPtr->position();
         #ifndef INTER_NOPAIR
         system().pairPotential().updateAtomCell(*thisPtr);
         #endif
         //system().moveAtom(*tempPtr, swapPos);
         tempPtr->position() = swapPos;
         #ifndef INTER_NOPAIR
         system().pairPotential().updateAtomCell(*tempPtr);
         #endif

//.........这里部分代码省略.........

   /* 
   * Generate, attempt and accept or reject a Monte Carlo move.
   */
   bool CfbHomoReptationMove::move() 
   {
      Vector    oldPos, newPos;
      double    rosen_r,  rosen_f;
      double    energy_r, energy_f;
      Atom     *tailPtr; // pointer to the tail atom (to be removed)
      Atom     *atomPtr; // resetable atom pointer
      Molecule *molPtr;  // pointer to randomly chosen molecule
      int       length, sign, headId, tailId, bondType, i;
      bool      accept;
     
      incrementNAttempt();

      // Choose a molecule at random
      molPtr = &(system().randomMolecule(speciesId_));
      length = molPtr->nAtom();

      // Choose which chain end to regrow
      if (random().uniform(0.0, 1.0) > 0.5) {
         sign = +1;
         headId = length - 1;
         tailId = 0;
      } else {
         sign = -1;
         headId = 0;
         tailId = length - 1;
      }

      // Store current positions of tail
      oldPos = molPtr->atom(tailId).position();
   
      // Delete tail monomers
      tailPtr = &(molPtr->atom(tailId));
      atomPtr = tailPtr + sign;
      if (sign == 1) {
         bondType = molPtr->bond(0).typeId();
      } else {
         bondType = molPtr->bond(length-2).typeId();
      }
      deleteEndAtom(tailPtr, atomPtr, bondType, rosen_r, energy_r);

      #ifndef INTER_NOPAIR
      // Delete from McSystem cell list
      system().pairPotential().deleteAtom(*tailPtr);
      #endif
   
      // Regrow head, using tailPtr to store properties of the new head.
      atomPtr = &(molPtr->atom(headId));
      tailPtr->setTypeId(atomPtr->typeId());
      tailPtr->mask().clear();
      #ifndef MCMD_NOMASKBONDED
      tailPtr->mask().append(*atomPtr);
      #endif
      if (sign == 1) {
         bondType = molPtr->bond(length-2).typeId();
      } else {
         bondType = molPtr->bond(0).typeId();
      }
      addEndAtom(tailPtr, atomPtr, bondType, rosen_f, energy_f);

      // Restore original type and connectivity of tail Atom
      atomPtr = tailPtr + sign;
      tailPtr->setTypeId(atomPtr->typeId());
      tailPtr->mask().clear();
      #ifndef MCMD_NOMASKBONDED
      tailPtr->mask().append(*atomPtr);
      #endif

      // Decide whether to accept or reject
      accept = random().metropolis(rosen_f/rosen_r);
      if (accept) {

         // Increment nAccept, the number accepted moves.
         incrementNAccept();

         // Store new head position
         newPos = tailPtr->position();
   
         tailPtr->position() = atomPtr->position();

         #ifndef INTER_NOPAIR
         // Add back to system cell list
         system().pairPotential().addAtom(*tailPtr);
         #endif
   
         // Shift atom positions towards the head
         for (i=1; i < length - 1; ++i) {
            //system().moveAtom(*atomPtr, (atomPtr+sign)->position());
            atomPtr->position() = (atomPtr+sign)->position();
            #ifndef INTER_NOPAIR
            system().pairPotential().updateAtomCell(*atomPtr);
            #endif
            atomPtr += sign;
         }
   
         // Move head atom to new chosen position
         //system().moveAtom(*atomPtr, newPos);
//.........这里部分代码省略.........

   /* 
   * Generate, attempt and accept or reject a Monte Carlo move.
   *
   * Convention for index
   *       head of molecule                tail of molecule
   *          |                                   |
   *          0   1   2   3 ...                length-1
   *          o---o---o---o---o---o---o---o---o---o
   *              |_______________|___________|
   *                valid beginId    nRegrow
   *
   */
   bool CfbRebridgeMove::move() 
   {
      Molecule *molPtr;   // pointer to randomly chosen molecule
      Atom     *thisPtr;  // pointer to the current end atom
      int       length, sign, beginId, endId, i;
      int       *bonds;
      double    rosen_r,  rosen_f;
      double    energy_r, energy_f;
      bool      accept;

      incrementNAttempt();

      // Choose a molecule at random
      molPtr = &(system().randomMolecule(speciesId_));
      length = molPtr->nAtom();

      // Require that chain length - 2 >= nRegrow_
      if (nRegrow_ > length - 2) {
         UTIL_THROW("nRegrow_  >= chain length");
      }

      // Allocate bonds array
      bonds = new int[nRegrow_ + 1];

      // Choose the beginId of the growing interior bridge
      beginId = random().uniformInt(1, length - nRegrow_);

      // Choose direction: sign = +1 if beginId < endId
      if (random().uniform(0.0, 1.0) > 0.5) {
         sign = +1;
         endId = beginId + (nRegrow_ - 1);
      } else {
         sign = -1;
         endId   = beginId;
         beginId = endId + (nRegrow_ - 1);
      }

      // Store current atomic positions from segment to be regrown
      thisPtr = &(molPtr->atom(beginId));
      for (i = 0; i < nRegrow_; ++i) {
         oldPos_[i] = thisPtr->position();
         thisPtr += sign;
      }

      // Get bond types.
      if (sign == +1) {
         for (i = 0; i < nRegrow_ + 1; ++i)
            bonds[i] = molPtr->bond(endId - i).typeId();
      } else {
         for (i = 0; i < nRegrow_ + 1; ++i)
            bonds[i] = molPtr->bond(endId - 1 + i).typeId();
      }

      // Delete atoms: endId -> beginId
      thisPtr = &(molPtr->atom(endId));
      deleteSequence(nRegrow_, sign, thisPtr, bonds, rosen_r, energy_r);

      // Regrow atoms: endId -> beginId
      thisPtr = &(molPtr->atom(beginId));
      addSequence(nRegrow_, sign, thisPtr, bonds, rosen_f, energy_f);

      // Release bonds array
      delete [] bonds;

      // Decide whether to accept or reject
      accept = random().metropolis(rosen_f/rosen_r);
      if (accept) {

         // Increment counter for accepted moves of this class.
         incrementNAccept();

         // If the move is accepted, keep current positions.

      } else {

         // If the move is rejected, restore old positions
         thisPtr = &(molPtr->atom(beginId));
         for (i = 0; i < nRegrow_; ++i) {
            //system().moveAtom(*thisPtr, oldPos_[i]);
            thisPtr->position() = oldPos_[i];
            #ifndef INTER_NOPAIR
            system().pairPotential().updateAtomCell(*thisPtr);
            #endif
            thisPtr += sign;
         }
   
      }

//.........这里部分代码省略.........

Molecule::bond_type bond_2(const Molecule &mol, const Molecule::atom_type &source, const Molecule::atom_type &target)
{
  return mol.bond(source, target);
}

Molecule::bond_type bond_1(const Molecule &mol, Helium::Index index)
{
  return mol.bond(index);
}