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

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());
}

TEST(HydrogenToolsTest, valencyAdjustment_O)
{
  Molecule mol;
  Atom O = mol.addAtom(8);
  int expectedAdjustment = 2;
  for (int i = 0; i < 8; ++i, --expectedAdjustment) {
    EXPECT_EQ(expectedAdjustment, HydrogenTools::valencyAdjustment(O));
    mol.addBond(mol.addAtom(1), O, 1);
  }
}

TEST(RingPerceiverTest, ethanol)
{
  Molecule molecule;
  molecule.addAtom(6);
  molecule.addAtom(6);
  molecule.addAtom(8);
  molecule.addBond(molecule.atom(0), molecule.atom(1), 1);
  molecule.addBond(molecule.atom(1), molecule.atom(2), 1);

  RingPerceiver perceiver(&molecule);
  std::vector<std::vector<size_t>> rings = perceiver.rings();
  EXPECT_EQ(rings.size(), static_cast<size_t>(0));
}

TEST(UnitCellTest, wrapAtomsToUnitCell)
{
  Molecule mol = createCrystal(static_cast<Real>(3.0),
                               static_cast<Real>(4.0),
                               static_cast<Real>(5.0),
                               static_cast<Real>(90.0),
                               static_cast<Real>(120.0),
                               static_cast<Real>(77.0));
  for (int i = 0; i < 10; ++i)
    mol.addAtom(1).setPosition3d(Vector3::Zero());

  Array<Vector3> fcoords;
  for (int i = 0; i < 10; ++i) {
    fcoords.push_back(Vector3(static_cast<Real>(i + i / 10.),
                              static_cast<Real>(i + 2 * i / 10.),
                              static_cast<Real>(i + 3 * i / 10.)));
  }

  EXPECT_TRUE(CrystalTools::setFractionalCoordinates(mol, fcoords));
  EXPECT_TRUE(CrystalTools::wrapAtomsToUnitCell(mol));
  fcoords.clear();
  EXPECT_TRUE(CrystalTools::fractionalCoordinates(mol, fcoords));
  for (std::vector<Vector3>::const_iterator it = fcoords.begin(),
       itEnd = fcoords.end(); it != itEnd; ++it) {
    EXPECT_GE(it->x(), static_cast<Real>(0.0));
    EXPECT_LE(it->x(), static_cast<Real>(1.0));
    EXPECT_GE(it->y(), static_cast<Real>(0.0));
    EXPECT_LE(it->y(), static_cast<Real>(1.0));
    EXPECT_GE(it->z(), static_cast<Real>(0.0));
    EXPECT_LE(it->z(), static_cast<Real>(1.0));
  }
}

void FileParser::readORCA()
{
    std::string tempString;
    QRegExp rx("", Qt::CaseInsensitive, QRegExp::RegExp2);

    while (1) {
        Molecule *molecule = new Molecule();

        // GOAL TO MATCH -
        //        ---------------------------------
        //        CARTESIAN COORDINATES (ANGSTROEM)
        //        ---------------------------------
        //          N     1.641610   -0.243155    1.175097
        //          C     2.719261   -0.702755    0.549942
        //          C     2.620982   -0.629038   -0.867321
        //          H     3.602940   -1.075919    1.062705
        //          N     1.452945   -0.127890   -1.286957
        //          H     3.405318   -0.967835   -1.537305
        //          N    -1.204954    0.202179    1.259883
        getline(infile, tempString);
        while (tempString.find("CARTESIAN COORDINATES (ANGSTROEM)") == string::npos &&
               infile.eof() == false) {
            getline(infile, tempString);
        }
        if (infile.eof()) {
            break;
        }
#ifdef QT_DEBUG
        std::cout << "readORCA: 'CARTESIAN COORDINATES (ANGSTROEM)' found.\n";
#endif

        // Geometry is reported in Angstroms
        myUnits = Angstrom;

        // Read in atom information
        rx.setPattern("(?:\\s*)(\\w+)(?:\\s+)(-?\\d+\\.\\d+)(?:\\s+)(-?\\d+\\.\\d+)(?:\\s+)(-?\\d+"
                      "\\.\\d+)(?:\\s*)");
        while (1) {
            getline(infile, tempString);
            if (rx.exactMatch(tempString.c_str()) == true) {
                AtomEntry *atom = new AtomEntry;
                atom->Label = rx.cap(1);
                atom->x = rx.cap(2).toDouble();
                atom->y = rx.cap(3).toDouble();
                atom->z = rx.cap(4).toDouble();
                molecule->addAtom(atom);
#ifdef QT_DEBUG
                std::cout << std::setw(5) << atom->Label.toStdString() << " " << std::setw(16)
                          << std::setprecision(10) << atom->x << " " << std::setw(16)
                          << std::setprecision(10) << atom->y << " " << std::setw(16)
                          << std::setprecision(10) << atom->z << std::endl;
#endif
            } else if (tempString.size() == 0) {
                myMoleculeList.push_back(molecule);
                break;
            }
        }
    }
}

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());
}

bool GaussianOutputFile::readInputOrientationGeometries(vector<Molecule>& molecules)
{
    Atom atom;
    Molecule mol;
    string line;
    bool ok;

    ok = false;

    int i;

    i = 0;
    while(isStringInFile("Input orientation"))
    {
//		if(isStringInFile("Input orientation"))
//		{
        /*
         * Jumping four lines that contains the orientation header
         *
         */
        for(size_t i = 0; i < 4; i++)
            line = readLine();

        /*
         * Reading the atoms
         *
         */
        line = readLine();
        while(line.find("-------") == string::npos)
        {
            atom.clear();
            atom.setName(line.substr(13,16));
            atom.setX(atof(line.substr(33,46).c_str()));
            atom.setY(atof(line.substr(47,59).c_str()));
            atom.setZ(atof(line.substr(60,70).c_str()));
            mol.addAtom(atom);
            line = readLine();
        }

        molecules.push_back(mol);

        mol.clear();

        ok = true;
//		}
    }

    if(ok)
        return true;
    else
        return false;
}

void HydrogenTools::adjustHydrogens(Molecule &molecule, Adjustment adjustment)
{
  // This vector stores indices of hydrogens that need to be removed. Additions
  // are made first, followed by removals to keep indexing sane.
  std::vector<size_t> badHIndices;

  // Temporary container for calls to generateNewHydrogenPositions.
  std::vector<Vector3> newHPos;

  // Convert the adjustment option to a couple of booleans
  bool doAdd(adjustment == Add || adjustment == AddAndRemove);
  bool doRemove(adjustment == Remove || adjustment == AddAndRemove);

  // Limit to only the original atoms:
  const size_t numAtoms = molecule.atomCount();

  // Iterate through all atoms in the molecule, adding hydrogens as needed
  // and building up a list of hydrogens that should be removed.
  for (size_t atomIndex = 0; atomIndex < numAtoms; ++atomIndex) {
    const Atom atom(molecule.atom(atomIndex));
    int hDiff = valencyAdjustment(atom);
    // Add hydrogens:
    if (doAdd && hDiff > 0) {
      newHPos.clear();
      generateNewHydrogenPositions(atom, hDiff, newHPos);
      for (std::vector<Vector3>::const_iterator it = newHPos.begin(),
           itEnd = newHPos.end(); it != itEnd; ++it) {
        Atom newH(molecule.addAtom(1));
        newH.setPosition3d(*it);
        molecule.addBond(atom, newH, 1);
      }
    }
    // Add bad hydrogens to our list of hydrogens to remove:
    else if (doRemove && hDiff < 0) {
      extraHydrogenIndices(atom, -hDiff, badHIndices);
    }
  }

  // Remove dead hydrogens now. Remove them in reverse-index order to keep
  // indexing sane.
  if (doRemove && !badHIndices.empty()) {
    std::sort(badHIndices.begin(), badHIndices.end());
    std::vector<size_t>::iterator newEnd(std::unique(badHIndices.begin(),
                                                     badHIndices.end()));
    badHIndices.resize(std::distance(badHIndices.begin(), newEnd));
    for (std::vector<size_t>::const_reverse_iterator it = badHIndices.rbegin(),
         itEnd = badHIndices.rend(); it != itEnd; ++it) {
      molecule.removeAtom(*it);
    }
  }
}

void MoleculeTest::removeAtom()
{
  // Should now to two atoms in the Molecule with ids 0 and 1.
  Atom *a = m_molecule->atom(1);
  m_molecule->removeAtom(a);
  QVERIFY(m_molecule->numAtoms() == 1);
  m_molecule->removeAtom(0ul);
  QVERIFY(m_molecule->numAtoms() == 0);

  // Check behavior of removing an atom that is owned by a different
  // molecule. Should not crash.
  Molecule mol;
  a = mol.addAtom();
  m_molecule->removeAtom(a);
}

bool AbbreviationExpander::tryCarbonChain (TokenChain &tokens, size_t &offset, Molecule &m, AttPoint &attach_to)
{
   if (attach_to.order != 1)
      return false;

   Token &cur = tokens[offset];
   if (cur.type != Token::Element)
      return false;

   if (cur.multiplier == 1 || cur.index != ELEM_C || offset + 1 == tokens.size())
      return false;

   Token &next = tokens[offset + 1];
   // Check CnH(2n+1), CnH2n pattern
   if (next.multiplier > 1 && next.index == ELEM_H)
   {
      bool tail;
      if (next.multiplier == cur.multiplier * 2)
      {
         // Intermediate carbon chain
         tail = false;
      }
      else if (next.multiplier == cur.multiplier * 2 + 1)
      {
         // Terminator carbon chain
         tail = true;
      }
      else
         return false;

      for (int i = 0; i < cur.multiplier; i++)
      {
         int idx = m.addAtom(ELEM_C);
         attachBond(m, attach_to, idx);

         attach_to = AttPoint(idx, 1);
      }
      if (tail)
         attach_to = AttPoint(-1, 0);
      offset += 2;
      return true;
   }
   return false;
}

TEST(HydrogenToolsTest, adjustHydrogens_adjustments)
{
  for (int i = 0; i < 3; ++i) {
    HydrogenTools::Adjustment adjustment;
    std::string expectedFormula;
    switch (i) {
    case 0:
      adjustment = HydrogenTools::Add;
      expectedFormula = "C2H14";
      break;
    case 1:
      adjustment = HydrogenTools::Remove;
      expectedFormula = "C2H5";
      break;
    case 2:
      adjustment = HydrogenTools::AddAndRemove;
      expectedFormula = "C2H8";
      break;
    }

    Molecule mol;
    Atom C1 = mol.addAtom(6); // Overbond this atom
    mol.addBond(C1, mol.addAtom(1));
    mol.addBond(C1, mol.addAtom(1));
    mol.addBond(C1, mol.addAtom(1));
    mol.addBond(C1, mol.addAtom(1));
    mol.addBond(C1, mol.addAtom(1));
    mol.addBond(C1, mol.addAtom(1));
    mol.addBond(C1, mol.addAtom(1));
    mol.addBond(C1, mol.addAtom(1));
    mol.addBond(C1, mol.addAtom(1));
    mol.addBond(C1, mol.addAtom(1));
    Atom C2 = mol.addAtom(6); // Underbond this atom
    mol.addBond(C2, mol.addAtom(1));

    EXPECT_EQ(std::string("C2H11"), mol.formula());
    HydrogenTools::adjustHydrogens(mol, adjustment);
    EXPECT_EQ(expectedFormula, mol.formula());
  }
}

bool GaussianOutputFile::readFirstInputOrientationGeometry(Molecule& mol)
{
    Atom atom;
    string line;

    if(!rewind())
        return false;

    if(isStringInFile("Input orientation"))
    {
        /*
         * Jumping four lines that contains the orientation header
         *
         */
        for(size_t i = 0; i < 4; i++)
            line = readLine();

        /*
         * Reading the atoms
         *
         */
        line = readLine();
        while(line.find("-------") == string::npos)
        {
            atom.clear();
            atom.setName(line.substr(13,16));
            atom.setX(atof(line.substr(33,46).c_str()));
            atom.setY(atof(line.substr(47,59).c_str()));
            atom.setZ(atof(line.substr(60,70).c_str()));
            mol.addAtom(atom);
            line = readLine();


        }

    }
    else
        return false;

    return true;
}

void MoleculeLayoutGraphSmart::saveDebug ()
{
   int i;
   Molecule mol;
   QS_DEF(Array<int>, mapping);

   mapping.clear_resize(vertexEnd());

   for (i = vertexBegin(); i < vertexEnd(); i = vertexNext(i))
   {
      if (getVertexType(i) == ELEMENT_NOT_DRAWN)
         continue;

      mapping[i] = mol.addAtom(ELEM_C);
      mol.setAtomXyz(mapping[i], getPos(i).x, getPos(i).y, 0);
   }

   for (i = edgeBegin(); i < edgeEnd(); i = edgeNext(i))
   {
      if (getEdgeType(i) == ELEMENT_NOT_DRAWN)
         continue;

      const Edge &edge = getEdge(i);

      mol.addBond(mapping[edge.beg], mapping[edge.end], BOND_SINGLE);
   }

   static int id = 0;
   char out_name[100];

   sprintf_s(out_name, "D:\\mf\\draw\\trace_my\\%03d.mol", id);

   FileOutput fo(out_name);
   MolfileSaver ms(fo);

   ms.saveMolecule(mol);

   id++;
}

//.........这里部分代码省略.........
      first_atom = false;

      if (bond != 0)
         bond->end = _atoms.size() - 1;

      memset(&atom, 0, sizeof(_AtomDesc));
      atom.hydrogens = -1;
      atom.valence = -1;
      atom.pseudo_atom_idx = -1;
      atom.rsite = false;

      if (code > 0 && (code < ELEM_MAX || code == CMF_PSEUDOATOM || code == CMF_RSITE || code == CMF_RSITE_EXT))
      {
         if (!_readAtom(code, atom, _atoms.size() - 1))
            break;
         continue;
      }

      if (!_getNextCode(code))
         break;

   } while (true); 

   // if have internal decoder, finish it
/*   if (_decoder_obj.get() != 0)
      _decoder_obj->finish(); */

   /* Reading finished, filling molecule */

   int i;

   for (i = 0; i < _atoms.size(); i++)
   {
      mol.addAtom(_atoms[i].label);

      if (_atoms[i].pseudo_atom_idx >= 0)
         mol.setPseudoAtom(i, _pseudo_labels.at(_atoms[i].pseudo_atom_idx));

      if (_atoms[i].rsite_bits > 0)
         mol.setRSiteBits(i, _atoms[i].rsite_bits);

      mol.setAtomCharge(i, _atoms[i].charge);
      mol.setAtomIsotope(i, _atoms[i].isotope);
      if (_atoms[i].hydrogens >= 0)
         mol.setImplicitH(i, _atoms[i].hydrogens);
      mol.setAtomRadical(i, _atoms[i].radical);

      if (_atoms[i].highlighted)
         mol.highlightAtom(i);
   }

   for (i = 0; i < _bonds.size(); i++)
   {
      int type = _bonds[i].type;
      int beg = _bonds[i].beg;
      int end = _bonds[i].end;
      int tmp;

      if (_bonds[i].swap)
         __swap(beg, end, tmp);

      int idx = mol.addBond_Silent(beg, end, type);
      
      if (_bonds[i].in_ring)
         mol.setEdgeTopology(idx, TOPOLOGY_RING);
      else

void FileParser::readQchem31()
{
    std::string tempString;
    QRegExp rx("", Qt::CaseInsensitive, QRegExp::RegExp2);

    while (1) {
        Molecule *molecule = new Molecule();

        // GOAL TO MATCH -
        // Optimization Cycle:   1
        //
        //                     Coordinates (Angstroms)
        //   ATOM              X           Y           Z
        //  1  H           6.264255   -1.214463    0.001562
        //  2  N           0.992577    0.405079   -0.000713
        //  3  C           3.640581    0.699759   -0.000285
        //  4  H           5.692879    1.292403   -0.000233
        //  5  N           1.103169   -1.927005    0.000687
        //  6  C          -3.986741   -0.820146    0.000147
        //  7  C          -4.603474    0.385503    0.001374
        //  8  N           4.158387   -1.540250    0.001279
        //  9  C           5.256835   -0.819403    0.001066
        // 10  H          -5.684216    0.487120    0.002358
        // 11  N           2.932990    1.843103   -0.001138
        // 12  H          -4.383326    2.467645    0.002362
        // 13  H           0.959209    2.465806   -0.001902
        // 14  H           1.658401   -2.773307    0.001857
        // 15  C           1.728674   -0.733977    0.000168
        // 16  N           5.011372    0.541687    0.000124
        // 17  H          -4.452049   -2.729414    0.878462
        // 18  C           3.136252   -0.602051    0.000394
        // 19  H          -0.845397    0.417252   -0.001289
        // 20  H          -4.453040   -2.728702   -0.879157
        // 21  H          -5.802705   -1.976283    0.000726
        // 22  C          -4.719039   -2.131115    0.000042
        // 23  C           1.619390    1.601695   -0.001312
        // 24  H           0.079495   -1.982939    0.000872
        // 25  N          -3.902148    1.574898    0.001480
        // 26  C          -2.511246    1.648434    0.000286
        // 27  O          -1.910441    2.717732    0.000430
        // 28  N          -1.895228    0.410754   -0.001007
        // 29  C          -2.522920   -0.833703   -0.001123
        // 30  O          -1.848461   -1.874351   -0.002325

        getline(infile, tempString);
        while (tempString.find("Optimization Cycle:") == string::npos && infile.eof() == false) {
            getline(infile, tempString);
        }
        if (infile.eof()) {
            break;
        }
#ifdef QT_DEBUG
        std::cout << "readQchem31: 'Optimization Cycle:' found.\n";
#endif

        getline(infile, tempString);
        getline(infile, tempString);
        getline(infile, tempString);

        // qchem31 is reported in Angstroms
        myUnits = Angstrom;

        // Read in atom information
        rx.setPattern("(?:\\s*)(?:\\d+)(?:\\s+)(\\w+)(?:\\s+)(-?\\d+\\.\\d+)(?:\\s+)(-?\\d+\\.\\d+)"
                      "(?:\\s+)(-?\\d+\\.\\d+)");
        while (1) {
            getline(infile, tempString);
            if (rx.exactMatch(tempString.c_str()) == true) {
                AtomEntry *atom = new AtomEntry;
                atom->Label = rx.cap(1);
                atom->x = rx.cap(2).toDouble();
                atom->y = rx.cap(3).toDouble();
                atom->z = rx.cap(4).toDouble();
                molecule->addAtom(atom);
#ifdef QT_DEBUG
                std::cout << std::setw(5) << atom->Label.toStdString() << " " << std::setw(16)
                          << std::setprecision(10) << atom->x << " " << std::setw(16)
                          << std::setprecision(10) << atom->y << " " << std::setw(16)
                          << std::setprecision(10) << atom->z << std::endl;
#endif
            } else {
                myMoleculeList.push_back(molecule);
                break;
            }
        }
    }
}