C++ ForceField::calcEnergy方法代码示例

void testIssue239(){
  BOOST_LOG(rdErrorLog) << "-------------------------------------" << std::endl;
  BOOST_LOG(rdErrorLog) << "    Testing Issue239." << std::endl;

  RWMol *mol;
  int needMore;
  ForceFields::ForceField *field;
  double e1,e2;

  std::string pathName=getenv("RDBASE");
  pathName += "/Code/GraphMol/ForceFieldHelpers/UFF/test_data";
  mol = MolFileToMol(pathName+"/Issue239.mol",false);
  TEST_ASSERT(mol);
  MolOps::sanitizeMol(*mol);

  field=UFF::constructForceField(*mol);
  TEST_ASSERT(field);
  field->initialize();
  needMore = field->minimize(200,1e-6,1e-3);
  e1 = field->calcEnergy();
  needMore = field->minimize(200,1e-6,1e-3);
  e2 = field->calcEnergy();
  TEST_ASSERT(feq(e2,e1,0.1));
  
  delete mol;
  delete field;
  

  BOOST_LOG(rdErrorLog) << "  done" << std::endl;
}

void testIssue242()
{
// FIX: Changes to the forcefield (connected to Issue 408) have
// made it so that this particular problem no longer manifests
// in this molecule/embedding. A new test case is needed.
  BOOST_LOG(rdErrorLog) << "-------------------------------------" << std::endl;
  BOOST_LOG(rdErrorLog) << "    Testing Issue242." << std::endl;

  RWMol *mol, *mol2;
  int needMore;
  ForceFields::ForceField *field = 0, *field2 = 0;
  std::string mb1,mb2;
  double e1,e2;

  std::string pathName = getenv("RDBASE");
  pathName += "/Code/GraphMol/ForceFieldHelpers/MMFF/test_data";

  mol = MolFileToMol(pathName + "/Issue242-2.mol");
  TEST_ASSERT(mol);

  mol2 = MolFileToMol(pathName + "/Issue242-2.mol");
  TEST_ASSERT(mol2);

  TEST_ASSERT(DGeomHelpers::EmbedMolecule(*mol2, 30, 2370) >= 0);
  mb1 = MolToMolBlock(*mol);
  mb2 = MolToMolBlock(*mol2);

  //std::cout << mb1 << "------\n";
  //std::cout << mb2 << "------\n";

  field = MMFF::constructForceField(*mol);
  TEST_ASSERT(field);
  field->initialize();
  field2 = MMFF::constructForceField(*mol2);
  TEST_ASSERT(field2);
  field2->initialize();
  e1 = field->calcEnergy();
  e2 = field2->calcEnergy();
  BOOST_LOG(rdInfoLog) << "E1: " << e1 << std::endl;
  BOOST_LOG(rdInfoLog) << "E2: " << e2 << std::endl;
  //TEST_ASSERT(feq(e2,e1,0.1));

  needMore = field->minimize(200,1.0e-4);
  needMore = field2->minimize(200,1.0e-4);
  e1 = field->calcEnergy();
  e2 = field2->calcEnergy();
  BOOST_LOG(rdInfoLog) << "E1: " << e1 << std::endl;
  BOOST_LOG(rdInfoLog) << "E2: " << e2 << std::endl;
  TEST_ASSERT(feq(e2, e1, 0.1));

  delete mol;
  delete mol2;
  delete field;
  delete field2;
  

  BOOST_LOG(rdErrorLog) << "  done" << std::endl;
}

void testUFFMultiThread(){
  BOOST_LOG(rdErrorLog) << "-------------------------------------" << std::endl;
  BOOST_LOG(rdErrorLog) << "    Test UFF multithreading" << std::endl;

  ForceFields::ForceField *field;

  std::string pathName = getenv("RDBASE");
  pathName += "/Code/GraphMol/ForceFieldHelpers/UFF/test_data";
  SDMolSupplier suppl(pathName + "/bulk.sdf");
  std::vector<ROMol *> mols;
  while(!suppl.atEnd()&&mols.size()<100){
    ROMol *mol=0;
    try{
      mol=suppl.next();
    } catch(...){
      continue;
    }
    if(!mol) continue;
    mols.push_back(mol);
  }

  std::cerr<<"generating reference data"<<std::endl;
  std::vector<double> energies(mols.size(),0.0);
  for(unsigned int i=0;i<mols.size();++i){
    ROMol mol(*mols[i]);
    ForceFields::ForceField *field = 0;
    try {
      field = UFF::constructForceField(mol);
    } catch (...) {
      field = 0;
    }
    TEST_ASSERT(field);
    field->initialize();
    int failed = field->minimize(500);
    TEST_ASSERT(!failed);
    energies[i]=field->calcEnergy();
    delete field;
  }
  
  boost::thread_group tg;

  std::cerr<<"processing"<<std::endl;
  unsigned int count=4;
  for(unsigned int i=0;i<count;++i){
    std::cerr<<" launch :"<<i<<std::endl;std::cerr.flush();
    tg.add_thread(new boost::thread(runblock_uff,mols,energies,count,i));
  }
  tg.join_all();
  
  BOOST_FOREACH(ROMol *mol,mols){
    delete mol;
  }
  BOOST_LOG(rdErrorLog) << "  done" << std::endl;
}

void testSFIssue2378119() {
  BOOST_LOG(rdErrorLog) << "-------------------------------------" << std::endl;
  BOOST_LOG(rdErrorLog) << "    Testing SFIssue2378119." << std::endl;

  std::string pathName = getenv("RDBASE");
  pathName += "/Code/GraphMol/ForceFieldHelpers/MMFF/test_data";
  {
    RWMol *mol = MolFileToMol(pathName + "/Issue2378119.mol");
    TEST_ASSERT(mol);
    ForceFields::ForceField *field = MMFF::constructForceField(*mol);
    TEST_ASSERT(field);
    field->initialize();
    double e1 = field->calcEnergy();
    TEST_ASSERT(e1 > 0.0 && e1 < 1.0e12);

    int needMore = field->minimize(200, 1.0e-6, 1.0e-3);
    TEST_ASSERT(!needMore);
    double e2 = field->calcEnergy();
    TEST_ASSERT(e2 < e1);

    delete mol;
    delete field;
  }
  {
    RWMol *mol = MolFileToMol(pathName + "/Issue2378119.2.mol");
    TEST_ASSERT(mol);
    ForceFields::ForceField *field = MMFF::constructForceField(*mol);
    TEST_ASSERT(field);
    field->initialize();
    double e1 = field->calcEnergy();
    TEST_ASSERT(e1 == 0.0);

    int needMore = field->minimize(200, 1.0e-6, 1.0e-3);
    TEST_ASSERT(!needMore);
    double e2 = field->calcEnergy();
    TEST_ASSERT(e2 == e1);

    delete mol;
    delete field;
  }
  {
    RWMol *mol = MolFileToMol(pathName + "/Issue2378119.2.mol");
    TEST_ASSERT(mol);
    ForceFields::ForceField *field =
        MMFF::constructForceField(*mol, 100.0, -1, false);
    TEST_ASSERT(field);
    field->initialize();
    double e1 = field->calcEnergy();
    TEST_ASSERT(e1 > 0.0 && e1 < 1.0e12);

    int needMore = field->minimize(200, 1.0e-6, 1.0e-3);
    TEST_ASSERT(!needMore);
    double e2 = field->calcEnergy();
    TEST_ASSERT(e2 < e1);

    delete mol;
    delete field;
  }
  BOOST_LOG(rdErrorLog) << "  done" << std::endl;
}

void testGitHubIssue62() {
  BOOST_LOG(rdErrorLog) << "-------------------------------------" << std::endl;
  BOOST_LOG(rdErrorLog) << "    Testing GitHubIssue62." << std::endl;

  std::string pathName=getenv("RDBASE");
  pathName += "/Code/GraphMol/ForceFieldHelpers/UFF/test_data";
  {
    double energyValues[] =
      { 38.687, 174.698, 337.986, 115.248,
        2.482, 1.918, 10.165, 98.469, 39.078,
        267.236, 15.747, 202.121, 205.539,
        20.044, 218.986, 79.627 };
    SmilesMolSupplier smiSupplier(pathName + "/Issue62.smi");
    SDWriter *sdfWriter = new SDWriter(pathName + "/Issue62.sdf");
    for (unsigned int i = 0; i < smiSupplier.length(); ++i) {
      ROMol *mol = MolOps::addHs(*(smiSupplier[i]));
      TEST_ASSERT(mol);
      std::string molName = "";
      if (mol->hasProp(common_properties::_Name)) {
        mol->getProp(common_properties::_Name, molName);
      }
      DGeomHelpers::EmbedMolecule(*mol);
      ForceFields::ForceField *field = UFF::constructForceField(*mol);
      TEST_ASSERT(field);
      field->initialize();
      int needMore = field->minimize(200, 1.e-6, 1.e-3);
      TEST_ASSERT(!needMore);
      sdfWriter->write(*mol);
      double e = field->calcEnergy();
      BOOST_LOG(rdErrorLog) << molName << " " << e << std::endl;
      TEST_ASSERT(fabs(e - energyValues[i]) < 1.);
    }
    sdfWriter->close();
    BOOST_LOG(rdErrorLog) << "  done" << std::endl;
  }
}

  PyObject *embedBoundsMatrix(python::object boundsMatArg,int maxIters=10,
                              bool randomizeOnFailure=false,int numZeroFail=2,
                              python::list weights=python::list(),
                              int randomSeed=-1){
    PyObject *boundsMatObj = boundsMatArg.ptr();
    if(!PyArray_Check(boundsMatObj))
      throw_value_error("Argument isn't an array");
    
    PyArrayObject *boundsMat=reinterpret_cast<PyArrayObject *>(boundsMatObj);
    // get the dimensions of the array
    unsigned int nrows = boundsMat->dimensions[0];
    unsigned int ncols = boundsMat->dimensions[1];
    if(nrows!=ncols)
      throw_value_error("The array has to be square");
    if(nrows<=0)
      throw_value_error("The array has to have a nonzero size");
    if (boundsMat->descr->type_num != PyArray_DOUBLE)
      throw_value_error("Only double arrays are currently supported");

    unsigned int dSize = nrows*nrows;
    double *cData = new double[dSize];
    double *inData = reinterpret_cast<double *>(boundsMat->data);
    memcpy(static_cast<void *>(cData), 
           static_cast<const void *>(inData),
           dSize*sizeof(double));

    DistGeom::BoundsMatrix::DATA_SPTR sdata(cData);
    DistGeom::BoundsMatrix bm(nrows,sdata);

    RDGeom::Point3D *positions=new RDGeom::Point3D[nrows];
    std::vector<RDGeom::Point *> posPtrs;
    for (unsigned int i = 0; i < nrows; i++) {
      posPtrs.push_back(&positions[i]);
    }

    RDNumeric::DoubleSymmMatrix distMat(nrows, 0.0);

    // ---- ---- ---- ---- ---- ---- ---- ---- ---- 
    // start the embedding:
    bool gotCoords=false;
    for(int iter=0;iter<maxIters && !gotCoords;iter++){
      // pick a random distance matrix
      DistGeom::pickRandomDistMat(bm,distMat,randomSeed);

      // and embed it:
      gotCoords=DistGeom::computeInitialCoords(distMat,posPtrs,randomizeOnFailure,
                                               numZeroFail,randomSeed);

      // update the seed:
      if(randomSeed>=0) randomSeed+=iter*999;
    }

    if(gotCoords){
      std::map<std::pair<int,int>,double> weightMap;
      unsigned int nElems=PySequence_Size(weights.ptr());
      for(unsigned int entryIdx=0;entryIdx<nElems;entryIdx++){
        PyObject *entry=PySequence_GetItem(weights.ptr(),entryIdx);
        if(!PySequence_Check(entry) || PySequence_Size(entry)!=3){
          throw_value_error("weights argument must be a sequence of 3-sequences");
        }
        int idx1=PyInt_AsLong(PySequence_GetItem(entry,0));
        int idx2=PyInt_AsLong(PySequence_GetItem(entry,1));
        double w=PyFloat_AsDouble(PySequence_GetItem(entry,2));
        weightMap[std::make_pair(idx1,idx2)]=w;
      }
      DistGeom::VECT_CHIRALSET csets;
      ForceFields::ForceField *field = DistGeom::constructForceField(bm,posPtrs,csets,0.0, 0.0,
                                                                     &weightMap);
      CHECK_INVARIANT(field,"could not build dgeom force field");
      field->initialize();
      if(field->calcEnergy()>1e-5){
        int needMore=1;
        while(needMore){
          needMore=field->minimize();
        }
      }
      delete field;
    } else {
      throw_value_error("could not embed matrix");
    }

    // ---- ---- ---- ---- ---- ---- ---- ---- ---- 
    // construct the results matrix:
    npy_intp dims[2];
    dims[0] = nrows;
    dims[1] = 3;
    PyArrayObject *res = (PyArrayObject *)PyArray_SimpleNew(2,dims,NPY_DOUBLE);
    double *resData=reinterpret_cast<double *>(res->data);
    for(unsigned int i=0;i<nrows;i++){
      unsigned int iTab=i*3;
      for (unsigned int j = 0; j < 3; ++j) {
        resData[iTab + j]=positions[i][j]; //.x;
      }
    }
    delete [] positions;

    return PyArray_Return(res);
  }