C++ OBMol::SetConformer方法代码示例

void runConformerComparisons(OBMol &moleculeA, OBMol &moleculeB, bool verbose=false) { // A is the target and B is the reference
    vector< vector<double> > coordAs, coordBs, comAs, comBs, eVectAs, eVectBs;
    vector<double> VDWsA, VDWsB, massesA, massesB, currentPCACoordA, currentSDCoordA, bestCoordsA;
    vector< vector<double> > atomMatchScoringTable;

    double molecularWeightA = moleculeA.GetMolWt(), molecularWeightB = moleculeB.GetMolWt();
    double bestVolumeOverlap = -(numeric_limits<double>::max)();
    int bestJ = -1, bestI = -1, stepCount = 0;

    generateVDWRadiusListFromMolecule(VDWsA, moleculeA);
    generateVDWRadiusListFromMolecule(VDWsB, moleculeB);
    generateAtomicMassesListFromMolecule(massesA, moleculeA);
    generateAtomicMassesListFromMolecule(massesB, moleculeB);
    generateCoordsMatrixFromMoleculeConformers(coordAs, moleculeA);
    generateCoordsMatrixFromMoleculeConformers(coordBs, moleculeB);
    getMoleculeConformerCenterCoords(comAs, moleculeA);
    getMoleculeConformerCenterCoords(comBs, moleculeB);
    generateAtomMatchScoringTableFromTwoMolecules(atomMatchScoringTable, moleculeA, moleculeB);
    eVectAs.resize( moleculeA.NumConformers() ), eVectBs.resize( moleculeB.NumConformers() );
    for (unsigned int k=0; k < moleculeA.NumConformers(); k++) getPCAEigenMatrix(eVectAs[k], coordAs[k], massesA, molecularWeightA);
    for (unsigned int k=0; k < moleculeB.NumConformers(); k++) getPCAEigenMatrix(eVectBs[k], coordBs[k], massesB, molecularWeightB);
    cout << "Finished setting up data; running search...\n";

    for (unsigned int j=0; j < moleculeB.NumConformers(); j++) {
        for (unsigned int i=0; i < moleculeA.NumConformers(); i++) {
            PCAEngine(currentPCACoordA, coordAs[i], coordBs[j], eVectAs[i], eVectBs[j], comAs[i], comBs[j], VDWsA, VDWsB, atomMatchScoringTable);
            double currentVolumeOverlap = steepestDescentEngine(currentSDCoordA, currentPCACoordA, coordBs[j], VDWsA, VDWsB, comAs[i], atomMatchScoringTable, 1.0, 10.0 * M_PI / 180.0);
            if (currentVolumeOverlap > bestVolumeOverlap) {
                bestVolumeOverlap = currentVolumeOverlap;
                bestCoordsA = currentSDCoordA;
                bestI = i;
                bestJ = j;
            }
            RANKS_AND_COORDS.push_back(RanksAndCoords(currentVolumeOverlap, i, j, currentSDCoordA));
            cout << "ROUND " << ++stepCount << " A#" << i << " and B#" << j << " = " << currentVolumeOverlap << endl;
        }
    }
    cout << "\nThe best overlap is between conformer A#" << bestI << " and B#" << bestJ << ", which, after PCA followed by Steepest Descent, produces a volume overlap of " << bestVolumeOverlap << endl;

    //saveCoordsMatrixToMolecule(moleculeA, bestCoordsA);

    addConformerToMolecule(moleculeA, bestCoordsA);
    moleculeA.SetConformer(moleculeA.NumConformers() - 1);
    moleculeB.SetConformer(bestJ);
    MOLECULE_CONFORMER = bestJ;
}

  bool OBMoleculeFormat::DoOutputOptions(OBBase* pOb, OBConversion* pConv)
  {
    if(pConv->IsOption("addoutindex", OBConversion::GENOPTIONS)) {
      stringstream ss;
      ss << pOb->GetTitle() << " " << pConv->GetOutputIndex();
      pOb->SetTitle(ss.str().c_str());
    }

    OBMol* pmol = dynamic_cast<OBMol*> (pOb);
    if(pmol) {
      if(pConv->IsOption("writeconformers", OBConversion::GENOPTIONS)) {
        //The last conformer is written in the calling function
        unsigned int c = 0;
        for (; c < pmol->NumConformers()-1; ++c) {
          pmol->SetConformer(c);
          if(!pConv->GetOutFormat()->WriteMolecule(pmol, pConv))
            break;
        }
        pmol->SetConformer(c);
      }
    }
    return true;
  }

  void OpConfab::Run(OBConversion* pConv, OBMol* pmol)
  {
    OBMol mol = *pmol;
    
    N++;
    cout << "**Molecule " << N << endl << "..title = " << mol.GetTitle() << endl;
    cout << "..number of rotatable bonds = " << mol.NumRotors() << endl;
    mol.AddHydrogens();
    bool success = pff->Setup(mol);
    if (!success) {
      cout << "!!Cannot set up forcefield for this molecule\n"
           << "!!Skipping\n" << endl;
      return;
    }
    pff->DiverseConfGen(rmsd_cutoff, conf_cutoff, energy_cutoff, verbose);

    pff->GetConformers(mol);
    int nconfs = include_original ? mol.NumConformers() : mol.NumConformers() - 1;
    unsigned int c = include_original ? 0 : 1;

    // If mol.NumRotors is 0 and originals have not been included, then nconfs
    // may be 0. Here, if nconfs is 0, we include the original input conformer
    if (nconfs == 0) {
      nconfs = mol.NumConformers();
      c = 0;
    }

    cout << "..generated " << nconfs << " conformers" << endl;

    for (; c < mol.NumConformers(); ++c) {
      mol.SetConformer(c);
      if(!pConv->GetOutFormat()->WriteMolecule(&mol, pConv))
        break;
    }
    cout << endl;

  }

int main(int argc,char *argv[])
{
  if (argc != 2)
    {
      cout << "Usage: obrotamer <file>" << endl;
      cout << "  Outputs the number of rotable bonds and generate a random"
           << " rotamer" << endl;
      return(-1);
    }

  ifstream ifs(argv[1]);
  if (!ifs)
    {
      cerr << "Error! Cannot read input file!" << endl;
      return(-1);
    }
  
  OBConversion conv(&ifs, &cout);
  OBFormat* pFormat;
  
  pFormat = conv.FormatFromExt(argv[1]);
  if ( pFormat == NULL )
    {
      cerr << "Error! Cannot read file format!" << endl;
      return(-1);
    }
  
  // Finally, we can do some work!
  OBMol mol;
  if (! conv.SetInAndOutFormats(pFormat, pFormat))
    {
      cerr << "Error! File format isn't loaded" << endl;
      return (-1);
    }
  
  OBRotorList rl;
  OBRotamerList rotamers;
  int *rotorKey = NULL;
  OBRandom rand;
  rand.TimeSeed();

  while(ifs.peek() != EOF && ifs.good())
    {
      mol.Clear();
      conv.Read(&mol);

      rl.Setup(mol);
      
      cerr << " Number of rotatable bonds: " << rl.Size() << endl;

      // indexed from 1, rotorKey[0] = 0
      std::vector<int> rotorKey(rl.Size() + 1, 0);

      // each entry represents the configuration of a rotor
      // e.g. indexes into OBRotor::GetResolution()
      //       (the different angles to sample from the OBRotorRules database)
      OBRotorIterator ri;
      OBRotor *rotor = rl.BeginRotor(ri);
      for (int i = 1; i < rl.Size() + 1; ++i, rotor = rl.NextRotor(ri))
        rotorKey[i] = rand.NextInt() % rotor->GetResolution().size();

      rotamers.SetBaseCoordinateSets(mol);
      rotamers.Setup(mol, rl);
      rotamers.AddRotamer(rotorKey);

      // This is more useful when you have added many rotamers
      // rather than the one in this example
      rotamers.ExpandConformerList(mol, mol.GetConformers());
      
      // Change the molecule conformation -- from 0 .. rotamers.NumRotamers()
      mol.SetConformer(0);
      conv.Write(&mol);
    } // while reading molecules
  
  return(0);
}