C++ Molecule类代码示例

TEST(CjsonTest, atoms)
{
  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));
  Atom atom = molecule.atom(0);
  EXPECT_EQ(atom.atomicNumber(), static_cast<unsigned char>(1));
  atom = molecule.atom(1);
  EXPECT_EQ(atom.atomicNumber(), static_cast<unsigned char>(6));
  EXPECT_EQ(atom.position3d().x(),  0.751621);
  EXPECT_EQ(atom.position3d().y(), -0.022441);
  EXPECT_EQ(atom.position3d().z(), -0.020839);

  atom = molecule.atom(7);
  EXPECT_EQ(atom.atomicNumber(), static_cast<unsigned char>(1));
  EXPECT_EQ(atom.position3d().x(), -1.184988);
  EXPECT_EQ(atom.position3d().y(),  0.004424);
  EXPECT_EQ(atom.position3d().z(), -0.987522);
}

bool AbbreviationExpander::expandAtomAbbreviation (Molecule &mol, int v)
{
   // Determine bonds configuration: number of bonds on the left and on the right
   Vec3f &pos = mol.getAtomXyz(v);
   const Vertex &vertex = mol.getVertex(v);
   int count_left = 0, count_right = 0, count_middle = 0;
   Array<int> left_atoms, right_atoms;
   for (int nei = vertex.neiBegin(); nei != vertex.neiEnd(); nei = vertex.neiNext(nei))
   {
      int nei_atom = vertex.neiVertex(nei);
      Vec3f &nei_pos = mol.getAtomXyz(nei_atom);
      int order = mol.getBondOrder(vertex.neiEdge(nei));
      if (nei_pos.x < pos.x)
      {
         count_left += order;
         left_atoms.push(nei_atom);
      }
      else
      {
         count_right += order;
         right_atoms.push(nei_atom);
      }

      Vec3f diff = nei_pos;
      diff.sub(pos);
      if (fabs(diff.y) > fabs(diff.x))
         count_middle += order;
   }

   int input_order, output_order;

   bool on_left = false, on_right = false, is_single = false;
   if (vertex.degree() == 1)
   {
      if (count_middle)
         on_left = on_right = true;
      if (count_left)
         on_right = true;
      if (count_right)
         on_left = true;
      input_order = std::max(count_left, count_right);
      output_order = 0;
      is_single = true;
   }
   else
   {
      on_right = true;
      input_order = count_left;
      output_order = count_right;
   }

   // Try to expand according to the directions
   Array<Options> options;

   if (on_right && !on_left)
   {
      options.push(Options(RIGHT, RIGHT));
      options.push(Options(LEFT, LEFT));
   }
   else
   {
      if (on_right)
         options.push(Options(RIGHT, RIGHT));
      if (on_left)
      {
         options.push(Options(LEFT, LEFT));
         options.push(Options(LEFT, RIGHT));
         options.push(Options(RIGHT, RIGHT));
         options.push(Options(LEFT, LEFT, 2));
      }
   }

   // Duplicate all the options with ignore case flag
   int opt_cnt = options.size();
   for (int i = 0; i < opt_cnt; i++)
   {
      Options opt = options[i];
      opt.ignore_case = true;
      options.push(opt);
   }

   bool found = false;
   Molecule expanded;
   for (int i = 0; i < options.size(); i++)
   {
      options[i].set(*this);

      if (expand(mol.getPseudoAtom(v), input_order, output_order, expanded))
      {
         found = true;
         break;
      }
   }

   if (!found)
      return false;

   // Merge expanded abbreviation with the source molecule and connect bonds
   Array<int> mapping;
   mol.mergeWithMolecule(expanded, &mapping);
//.........这里部分代码省略.........

static PyObject *set(PyObject *self, PyObject *args) {
  int i, molid, frame;
  PyObject *selected, *val;
  char *attr = 0;

  //
  // get molid, frame, list, attribute, and value
  //
  if (!PyArg_ParseTuple(args, (char *)"iiO!sO!", &molid, &frame,
                        &PyTuple_Type, &selected, 
                        &attr, &PyTuple_Type, &val ))
    return NULL;  // bad args

  // 
  // check that we have been given either one value or one for each selected
  // atom
  //
  int num_selected = PyTuple_Size(selected);
  int tuplesize = PyTuple_Size(val);
  if (tuplesize != 1 && tuplesize != num_selected) {
    PyErr_SetString(PyExc_ValueError, "wrong number of items");
    return NULL; 
  }
 
  //
  // check molecule
  //
  VMDApp *app = get_vmdapp();
  Molecule *mol = app->moleculeList->mol_from_id(molid);
  if (!mol) {
    PyErr_SetString(PyExc_ValueError, "molecule no longer exists");
    return NULL;
  }
  const int num_atoms = mol->nAtoms;

  //
  // Check for a valid attribute
  //
  SymbolTable *table = app->atomSelParser;
  int attrib_index = table->find_attribute(attr);
  if (attrib_index == -1) {
    PyErr_SetString(PyExc_ValueError, "unknown atom attribute");
    return NULL;
  }
  SymbolTableElement *elem = table->fctns.data(attrib_index);
  if (elem->is_a != SymbolTableElement::KEYWORD &&
      elem->is_a != SymbolTableElement::SINGLEWORD) {
    PyErr_SetString(PyExc_ValueError, "attribute is not a keyword or singleword");
    return NULL;
  }
  if (!table->is_changeable(attrib_index)) {
    PyErr_SetString(PyExc_ValueError, "attribute is not modifiable");
    return NULL; 
  }

  // 
  // convert the list of selected atoms into an array of integer flags
  //
  // XXX should check that selected contains valid indices
  int *flgs = new int[num_atoms];
  memset(flgs,0,num_atoms*sizeof(int));
  for (i=0; i<num_selected; i++)
    flgs[PyInt_AsLong(PyTuple_GET_ITEM(selected,i))] = 1;
 
  //  
  // set the data
  //

  // singlewords can never be set, so macro is NULL.
  atomsel_ctxt context(table, mol, frame, NULL);
  if (elem->returns_a == SymbolTableElement::IS_INT) {
    int *list = new int[num_atoms];
    if (tuplesize > 1) {
      int j=0;
      for (int i=0; i<num_atoms; i++) {
        if (flgs[i])
          list[i] = PyInt_AsLong(PyTuple_GET_ITEM(val, j++));
      }
    } else {
      for (int i=0; i<num_atoms; i++) {
        if (flgs[i])
          list[i] = PyInt_AsLong(PyTuple_GET_ITEM(val, 0));
      }
    }
    elem->set_keyword_int(&context, num_atoms, list, flgs);
    delete [] list;

  } else if (elem->returns_a == SymbolTableElement::IS_FLOAT) {
    double *list = new double[num_atoms];
    if (tuplesize > 1) { 
      int j=0;
      for (int i=0; i<num_atoms; i++) { 
        if (flgs[i])
          list[i] = PyFloat_AsDouble(PyTuple_GET_ITEM(val, j++));
      }
    } else {
      for (int i=0; i<num_atoms; i++) {
        if (flgs[i])
          list[i] = PyFloat_AsDouble(PyTuple_GET_ITEM(val, 0));
      }
//.........这里部分代码省略.........

  void NPT::moveA() {
    SimInfo::MoleculeIterator i;
    Molecule::IntegrableObjectIterator  j;
    Molecule* mol;
    StuntDouble* sd;
    Vector3d Tb, ji;
    RealType mass;
    Vector3d vel;
    Vector3d pos;
    Vector3d frc;
    Vector3d sc;
    int index;

    thermostat = snap->getThermostat();
    loadEta();
    
    instaTemp =thermo.getTemperature();
    press = thermo.getPressureTensor();
    instaPress = PhysicalConstants::pressureConvert* (press(0, 0) + press(1, 1) + press(2, 2)) / 3.0;
    instaVol =thermo.getVolume();

    Vector3d  COM = thermo.getCom();

    //evolve velocity half step

    calcVelScale();

    for (mol = info_->beginMolecule(i); mol != NULL; 
         mol = info_->nextMolecule(i)) {

      for (sd = mol->beginIntegrableObject(j); sd != NULL;
	   sd = mol->nextIntegrableObject(j)) {
                
	vel = sd->getVel();
	frc = sd->getFrc();

	mass = sd->getMass();

	getVelScaleA(sc, vel);

	// velocity half step  (use chi from previous step here):

	vel += dt2*PhysicalConstants::energyConvert/mass* frc - dt2*sc;
	sd->setVel(vel);

	if (sd->isDirectional()) {

	  // get and convert the torque to body frame

	  Tb = sd->lab2Body(sd->getTrq());

	  // get the angular momentum, and propagate a half step

	  ji = sd->getJ();

	  ji += dt2*PhysicalConstants::energyConvert * Tb 
            - dt2*thermostat.first* ji;
                
	  rotAlgo_->rotate(sd, ji, dt);

	  sd->setJ(ji);
	}
            
      }
    }
    // evolve chi and eta  half step

    thermostat.first += dt2 * (instaTemp / targetTemp - 1.0) / tt2;
    
    evolveEtaA();

    //calculate the integral of chidt
    thermostat.second += dt2 * thermostat.first;
    
    flucQ_->moveA();


    index = 0;
    for (mol = info_->beginMolecule(i); mol != NULL; 
         mol = info_->nextMolecule(i)) {

      for (sd = mol->beginIntegrableObject(j); sd != NULL;
	   sd = mol->nextIntegrableObject(j)) {

	oldPos[index++] = sd->getPos();            

      }
    }
    
    //the first estimation of r(t+dt) is equal to  r(t)

    for(int k = 0; k < maxIterNum_; k++) {
      index = 0;
      for (mol = info_->beginMolecule(i); mol != NULL; 
           mol = info_->nextMolecule(i)) {

	for (sd = mol->beginIntegrableObject(j); sd != NULL;
	     sd = mol->nextIntegrableObject(j)) {

	  vel = sd->getVel();
//.........这里部分代码省略.........

int main(int argc, char* argv[])
{
    CommandlineParser parpars("RMSDCalculator", "calculate RMSD between ligand poses", VERSION, String(__DATE__), "Analysis");
    parpars.registerMandatoryInputFile("i", "input molecule file");
    parpars.registerMandatoryInputFile("org", "molecule file containing the original ('true') poses");
    parpars.registerOptionalOutputFile("o", "output molecule file");
    parpars.registerFlag("quiet", "by quiet, i.e. do not print progress information");
    String man = "This tool calculates the RMSD between different conformations of the same molecule.\n\nThis tool can be used to evaluate the differences between ligand poses taken from co-crystal structures, e.g. generated by a docking run.\nNote:Molecules may be sorted differently in the two input files; a topology hashkey will be used to match molecules to each other.\n\nOutput of this tool is a molecule file which will for each molecule contain a property-tag 'RMSD' holding the calculated RMSD value.";
    parpars.setToolManual(man);
    parpars.setSupportedFormats("i",MolFileFactory::getSupportedFormats());
    parpars.setSupportedFormats("org",MolFileFactory::getSupportedFormats());
    parpars.setSupportedFormats("o","mol2,sdf,drf");
    parpars.parse(argc, argv);

    // Retrieve coordinates of original poses
    GenericMolFile* original = MolFileFactory::open(parpars.get("org"));
    HashMap<String, list<Vector3> > original_poses;
    for (Molecule* mol = original->read(); mol; delete mol, mol = original->read())
    {
        String topology_hash;
        FlexibleMolecule::generateTopologyHash(mol, topology_hash, true);
        if (original_poses.find(topology_hash) != original_poses.end())
        {
            Log<<"[Warning:] more than one 'original' conformation for a molecule detected. Will use only the first conformation and ignore all other."<<endl;
        }
        else
        {
            list<Vector3> l;
            HashMap<String, list<Vector3> >::iterator map_it = original_poses.insert(make_pair(topology_hash, l)).first;

            for (AtomConstIterator it = mol->beginAtom(); +it; it++)
            {
                if (it->getElement().getSymbol() != "H")
                {
                    map_it->second.push_back(it->getPosition());
                }
            }
        }
    }
    delete original;

    // Retrieve coordinates of input poses and calculate RMSDs
    GenericMolFile* input = MolFileFactory::open(parpars.get("i"));
    GenericMolFile* output = 0;
    String filename = parpars.get("o");
    if (filename != CommandlineParser::NOT_FOUND)
    {
        output = MolFileFactory::open(filename, ios::out, input);
    }

    double average_RMSD = 0;
    int no_mols = 0;
    int no_valid_rmsds = 0;
    bool quiet = (parpars.get("quiet")!=CommandlineParser::NOT_FOUND);

    for (Molecule* mol = input->read(); mol; delete mol, mol = input->read())
    {
        no_mols++;
        String topology_hash;
        FlexibleMolecule::generateTopologyHash(mol, topology_hash, true);

        HashMap<String, list<Vector3> >::iterator map_it = original_poses.find(topology_hash);
        if (map_it == original_poses.end())
        {
            Log<<"[Warning:] no original pose for molecule '"<<mol->getName()<<"' found, its RMSD can thus not be computed."<<endl;
            mol->setProperty("RMSD", "N/A");
        }
        else
        {
            double RMSD = 0;
            list<Vector3>::iterator list_it = map_it->second.begin();
            int no_heavy_atoms = 0;
            AtomConstIterator it = mol->beginAtom();
            for (; +it ; it++)
            {
                if (it->getElement().getSymbol() != "H" && list_it != map_it->second.end())
                {
                    RMSD += pow(it->getPosition().getDistance(*list_it), 2);
                    no_heavy_atoms++;
                    list_it++;
                }
            }
            if (it != mol->endAtom() || list_it != map_it->second.end())
            {
                Log.error()<<"[Error:] Number of heavy atoms of input pose do not match number of heavy atoms of original pose!!"<<endl;
                return 1;
            }
            RMSD = sqrt(RMSD/no_heavy_atoms);
            mol->setProperty("RMSD", RMSD);
            average_RMSD += RMSD;
            no_valid_rmsds++;

            if (!quiet) Log << "RMSD for molecule "<<no_mols<<", '"<<mol->getName()<<"' = "<<RMSD<<endl;
        }

        if (output) *output << *mol;
    }

    average_RMSD /= no_valid_rmsds;

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

bool ChemicalPotential::getDeletion(TMoleculeContainer* cell, double* minco, double* maxco)
{
	 if(this->remainingDeletions.empty()) return false;
	 if(cell->getNumberOfParticles() == 0) return false;
	 
	 unsigned idx = *this->remainingDeletions.begin();
	 this->remainingDeletions.erase(this->remainingDeletions.begin());
	 double tminco[3];
	 double tmaxco[3];
	 if(restrictedControlVolume) for(int d=0; d < 3; d++)
	 {
	    tminco[d] = (minco[d] > control_bottom[d])? minco[d]
	                                              : control_bottom[d];
	    tmaxco[d] = (maxco[d] < control_top[d])? maxco[d]
	                                           : control_top[d];
	 }
	 else for(int d=0; d < 3; d++)
	 {
	    tminco[d] = minco[d];
	    tmaxco[d] = maxco[d];
	 }
	 
	 Molecule* m = cell->begin();
	 int j=0;
	 for(unsigned i=0; (i < idx); i++)
	 {
	    while(( (m->r(0) > tmaxco[0]) || (m->r(1) > tmaxco[1]) ||
	      (m->r(2) > tmaxco[2]) || (m->r(0) < tminco[0]) ||
	      (m->r(1) < tminco[1]) || (m->r(2) < tminco[2]) ||
							(m->componentid() != this->componentid) )
					 && (m != cell->end()))
			{
				 m = cell->next();
				 if(m == cell->end())
	 {
	    if(j == 0) return false;
	    m = cell->begin();
	    j = 0;
	 }
			}
			m = cell->next();
			j++;
			if(m == cell->end())
			{
	 if(j == 0) return false;
	 m = cell->begin();
	 j = 0;
			}
	 }
	 while( (m->r(0) > tmaxco[0]) || (m->r(1) > tmaxco[1]) ||
	        (m->r(2) > tmaxco[2]) || (m->r(0) < tminco[0]) ||
	  (m->r(1) < tminco[1]) || (m->r(2) < tminco[2]) ||
					(m->componentid() != this->componentid) )
	 {
	    m = cell->next();
	    if(m == cell->end())
	    {
	 if(j == 0) return false;
	 m = cell->begin();
			}
	 }
#ifndef NDEBUG
	 global_log->debug() << "ID " << m->id() << " selected for deletion (index " << idx << ")." << std::endl;
#endif
	 assert(m->id() < nextid);
	 return true;
}

  void DensityPlot::process() {
    Molecule* mol;
    RigidBody* rb;
    SimInfo::MoleculeIterator mi;
    Molecule::RigidBodyIterator rbIter;

    DumpReader reader(info_, dumpFilename_);    
    int nFrames = reader.getNFrames();
    for (int i = 0; i < nFrames; i += step_) {
      reader.readFrame(i);
      currentSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot();

      for (mol = info_->beginMolecule(mi); mol != NULL; 
	   mol = info_->nextMolecule(mi)) {
        //change the positions of atoms which belong to the rigidbodies
        for (rb = mol->beginRigidBody(rbIter); rb != NULL; 
	     rb = mol->nextRigidBody(rbIter)) {
          rb->updateAtoms();
        }
        
      }
    
      if (evaluator_.isDynamic()) {
	seleMan_.setSelectionSet(evaluator_.evaluate());
      }

      if (cmEvaluator_.isDynamic()) {
	cmSeleMan_.setSelectionSet(cmEvaluator_.evaluate());
      }

      Vector3d origin = calcNewOrigin();

      Mat3x3d hmat = currentSnapshot_->getHmat();
      RealType slabVolume = deltaR_ * hmat(0, 0) * hmat(1, 1);
      int k; 
      for (StuntDouble* sd = seleMan_.beginSelected(k); sd != NULL; 
	   sd = seleMan_.nextSelected(k)) {


        if (!sd->isAtom()) {
          sprintf( painCave.errMsg, 
		   "Can not calculate electron density if it is not atom\n");
          painCave.severity = OPENMD_ERROR;
          painCave.isFatal = 1;
          simError(); 
        }
            
        Atom* atom = static_cast<Atom*>(sd);
        GenericData* data = atom->getAtomType()->getPropertyByName("nelectron");
        if (data == NULL) {
          sprintf( painCave.errMsg, "Can not find Parameters for nelectron\n");
          painCave.severity = OPENMD_ERROR;
          painCave.isFatal = 1;
          simError(); 
        }
            
        DoubleGenericData* doubleData = dynamic_cast<DoubleGenericData*>(data);
        if (doubleData == NULL) {
          sprintf( painCave.errMsg,
                   "Can not cast GenericData to DoubleGenericData\n");
          painCave.severity = OPENMD_ERROR;
          painCave.isFatal = 1;
          simError();   
        }
            
        RealType nelectron = doubleData->getData();
        LennardJonesAdapter lja = LennardJonesAdapter(atom->getAtomType());
        RealType sigma = lja.getSigma() * 0.5;
        RealType sigma2 = sigma * sigma;
            
        Vector3d pos = sd->getPos() - origin;
        for (int j =0; j < nRBins_; ++j) {
          Vector3d tmp(pos);
          RealType zdist =j * deltaR_ - halfLen_;
          tmp[2] += zdist;
          if (usePeriodicBoundaryConditions_) 
            currentSnapshot_->wrapVector(tmp);
              
          RealType wrappedZdist = tmp.z() + halfLen_;
          if (wrappedZdist < 0.0 || wrappedZdist > len_) {
            continue;
          }
              
          int which = int(wrappedZdist / deltaR_);
          density_[which] += nelectron * exp(-zdist*zdist/(sigma2*2.0)) /(slabVolume* sqrt(2*NumericConstant::PI*sigma*sigma));
              
        }            
      }        
    }
  
    int nProcessed = nFrames /step_;
    std::transform(density_.begin(), density_.end(), density_.begin(), 
		   std::bind2nd(std::divides<RealType>(), nProcessed));  
    writeDensity();
        

  
  }

  void Shake::doConstraint(ConstraintPairFuncPtr func) {
    if (!doShake_) return;

    Molecule* mol;
    SimInfo::MoleculeIterator mi;
    ConstraintElem* consElem;
    Molecule::ConstraintElemIterator cei;
    ConstraintPair* consPair;
    Molecule::ConstraintPairIterator cpi;
    
    for (mol = info_->beginMolecule(mi); mol != NULL; 
         mol = info_->nextMolecule(mi)) {
      for (consElem = mol->beginConstraintElem(cei); consElem != NULL; 
           consElem = mol->nextConstraintElem(cei)) {
	consElem->setMoved(true);
	consElem->setMoving(false);
      }
    }
    
    //main loop of constraint algorithm
    int done = 0;
    int iteration = 0;
    while(!done && iteration < maxConsIteration_){
      done = 1;

      //loop over every constraint pair

      for (mol = info_->beginMolecule(mi); mol != NULL; 
           mol = info_->nextMolecule(mi)) {
	for (consPair = mol->beginConstraintPair(cpi); consPair != NULL; 
             consPair = mol->nextConstraintPair(cpi)) {
          

	  //dispatch constraint algorithm
	  if(consPair->isMoved()) {
	    int exeStatus = (this->*func)(consPair);

	    switch(exeStatus){
	    case consFail:
	      sprintf(painCave.errMsg,
		      "Constraint failure in Shake::constrainA, "
                      "Constraint Fail\n");
	      painCave.isFatal = 1;
	      simError();                             
                            
	      break;
	    case consSuccess:
	      // constrain the pair by moving two elements
	      done = 0;
	      consPair->getConsElem1()->setMoving(true);
	      consPair->getConsElem2()->setMoving(true);
	      break;
	    case consAlready:
	      // current pair is already constrained, do not need to
	      // move the elements
	      break;
	    default:          
	      sprintf(painCave.errMsg, "ConstraintAlgorithm::doConstraint() "
                      "Error: unrecognized status");
	      painCave.isFatal = 1;
	      simError();                           
	      break;
	    }      
	  }
	}
      }//end for(iter->first())

#ifdef IS_MPI
      MPI_Allreduce(MPI_IN_PLACE, &done, 1, MPI_INT, MPI_LAND, MPI_COMM_WORLD);
#endif

      errorCheckPoint();

      for (mol = info_->beginMolecule(mi); mol != NULL; 
           mol = info_->nextMolecule(mi)) {
	for (consElem = mol->beginConstraintElem(cei); consElem != NULL; 
             consElem = mol->nextConstraintElem(cei)) {
	  consElem->setMoved(consElem->getMoving());
	  consElem->setMoving(false);
	}
      }

      iteration++;
    }//end while

    if (!done){
      sprintf(painCave.errMsg,
              "Constraint failure in Shake::constrainA, "
              "too many iterations: %d\n",
              iteration);
      painCave.isFatal = 1;
      simError();    
    }
    
    errorCheckPoint();
  }

   /*
   * Evaluate contributions to accumulators.
   */
   void MdPairEnergyCoefficients::sample(long iStep) 
   {
      if (isAtInterval(iStep)) {

         #ifndef INTER_NOPAIR
//            if (!system().isPairListCurrent()) {
//               system().buildPairList();
//             }

            // First clear neighbor list
            clear();

            // Loop over pairs and construct neighbor list per molecule
            PairIterator iter;
            Atom *atom0Ptr;
            Atom *atom1Ptr;

            for (pairListPtr_->begin(iter); iter.notEnd(); ++iter) {
               iter.getPair(atom0Ptr, atom1Ptr);
               if (selector_.match(*atom0Ptr, *atom1Ptr)) {
                  Pair< Atom * > atomPair;
                  Molecule *moleculePtr;
                  Species *speciesPtr;
                  int speciesId, moleculeId;

                  atomPair[0] = atom0Ptr;
                  atomPair[1] = atom1Ptr;

                  // Store pair in neighbor list of boths atoms' molecules
                  // (do double counting, since pair list counts only unique
                  // pairs)
                  moleculePtr = & atom0Ptr->molecule();
                  speciesPtr =  & moleculePtr->species();
                  speciesId = speciesPtr->id();
                  moleculeId = moleculePtr->id();

                  moleculeNeighbors_[speciesId][moleculeId].
                     append(atomPair);

                  moleculePtr = & atom1Ptr->molecule();
                  speciesPtr =  & moleculePtr->species();
                  speciesId = speciesPtr->id();
                  moleculeId = moleculePtr->id();

                  // store atomPair in reverse order
                  atomPair[0] = atom1Ptr;
                  atomPair[1] = atom0Ptr;

                  moleculeNeighbors_[speciesId][moleculeId].
                     append(atomPair);
               } 
            } // end pair loop

            // Now loop over all molecules
            double moleculePESq = 0.0;
            double pairEnergy=0.0;
            double twoMoleculePESq=0.0;

            double rsq;
            int iSpecies1, iMolecule1;

            for (iSpecies1 = 0; iSpecies1 < nSpecies_; ++iSpecies1) {
               Species *species1Ptr;

               species1Ptr = &system().simulation().species(iSpecies1); 
               for (iMolecule1 = 0; iMolecule1 < species1Ptr->capacity();
                  ++iMolecule1) {
                  int iSpecies2, iMolecule2;

                  double moleculePairEnergy=0.0;
                  
                  // clear array for pair energy with neighboring molecules
                  for (iSpecies2 = 0; iSpecies2 < nSpecies_; ++iSpecies2) {
                     Species *species2Ptr;

                     species2Ptr = &system().simulation().species(iSpecies2);
                     for (iMolecule2 = 0; iMolecule2 < species2Ptr->capacity();
                        ++iMolecule2) {
                        twoMoleculePairEnergy_[iSpecies2][iMolecule2]=0.0;
                     }
                  }

                  // Loop over this molecule's neighbors
                  DSArray< Pair< Atom *> > *neighborsPtr;
                  ConstArrayIterator< Pair< Atom *> > it;

                  neighborsPtr=&moleculeNeighbors_[iSpecies1][iMolecule1];
                  for (neighborsPtr->begin(it); it.notEnd(); ++it) {
                     Atom *atom0Ptr, *atom1Ptr;
                     atom0Ptr = (*it)[0];
                     atom1Ptr = (*it)[1];

                     if (selector_.match(*atom0Ptr,*atom1Ptr)) {
                        double energy;
                        Species *species2Ptr;
                        Molecule *molecule2Ptr;
                        int species2Id,molecule2Id;
//.........这里部分代码省略.........

	bool GAMESSLogFile::write(const Molecule& molecule)
	{
		System S;
		S.insert(*(Molecule*)molecule.create(true));
		return write(S);
	}

void mergeDRFiles(vector<String>& names, string& output_file, Size& best_k, string& e_property, double& score_cutoff, double& score_cuton)
{

	DockResultFile* output = new DockResultFile(output_file, ios::out);
	bool sort_by_scores = 1;
	if (e_property == "") sort_by_scores = 0;

	vector<Result*> new_results;

	/// First of all, copy Result data
	map<Result::Method, Result*> result_map;
	for (Size file = 0; file < names.size(); file++)
	{
		DockResultFile* input = new DockResultFile(names[file]);

		const vector<Result*>* results = input->getResults();
		for (Size i = 0; i < results->size(); i++)
		{
			map<Result::Method, Result*>::iterator it = result_map.find((*results)[i]->getMethod());
			if (it == result_map.end())
			{
				Result* result_copy = new Result(*(*results)[i]);
				if (!sort_by_scores) output->addResult(result_copy);
				else new_results.push_back(result_copy);
				result_map.insert(make_pair(result_copy->getMethod(), result_copy));
			}
			else
			{
				*it->second += *(*results)[i];
			}
		}

		input->close();
		delete input;
	}


	if (e_property != "")
	{
		e_property = "score_"+new_results.back()->getMethodString();
	}

	/// If no sorting is desired, iterate over all input-files and write each input-molecules to output-file
	if (!sort_by_scores)
	{
		output->disableAutomaticResultCreation();

		for (Size file = 0; file < names.size(); file++)
		{
			GenericMolFile* input = MolFileFactory::open(names[file]);

			int mol_no = 0;
			for (Molecule* mol = input->read(); mol; mol = input->read(), mol_no++)
			{
				*output << *mol;
				delete mol;
				Log.level(20) << "\r" << names[file] << " : " << mol_no+1;
				Log.flush();
			}
			Log.level(20)<<endl;
			Log.flush();

			input->close();
			delete input;
		}
	}

	/// If sorting is desired, iterate over all input-files and save each input-molecules to a map.
	/// Then write all FlexibleMolecules in this map to the output file and adapt the Result objects.
	else
	{
		multimap < double, FlexibleMolecule* > compounds; // map containing score and conformation-ID
		set < String > IDs; // IDs of the base-conformations

		for (Size file = 0; file < names.size(); file++)
		{
			DockResultFile* input = new DockResultFile(names[file]);
			int mol_no = 0;
			for (Molecule* mol = input->read(); mol; mol = input->read(), mol_no++)
			{
				if (!mol->hasProperty(e_property))
				{
					Log.level(10) << "Compound " << mol->getName() << " in file " << names[file] << " has no score property. Skipping this compound." << endl;
					for (Size i = 0; i < new_results.size(); i++)
					{
						new_results[i]->erase(input->getCurrentLigand());
					}
					delete mol;
					continue;
				}

				double score = ((String)mol->getProperty(e_property).toString()).toFloat();

				if (score > score_cutoff || score < score_cuton)
				{
					for (Size i = 0; i < new_results.size(); i++)
					{
						new_results[i]->erase(input->getCurrentLigand());
					}
					delete mol;
//.........这里部分代码省略.........

void sortMolecules(vector<String>& names, string& output_file, Size& best_k, string& e_property, double& score_cutoff, double& score_cuton)
{
	multimap<double, Molecule*> compounds;

	for (Size file = 0; file < names.size(); file++)
	{
		GenericMolFile* input = MolFileFactory::open(names[file]);

		int mol_no = 0;
		for (Molecule* mol = input->read(); mol; mol = input->read(), mol_no++)
		{
			if (!mol->hasProperty(e_property))
			{
				Log.level(10) << "Compound " << mol->getName() << " in file " << names[file] << " has no score property. Skipping this compound." << endl;
				delete mol;
				continue;
			}

			double score = ((String)mol->getProperty(e_property).toString()).toFloat();

			if (score > score_cutoff || score < score_cuton)
			{
				delete mol;
				continue;
			}

			if ((compounds.size() < best_k || score < compounds.rbegin()->first))
			{
				compounds.insert(make_pair(score, mol));

				if (compounds.size() > best_k)
				{
					delete compounds.rbegin()->second;
					multimap<double, Molecule*>::iterator it = compounds.end();
					it--;
					compounds.erase(it);
				}
			}
			else
			{
				delete mol;
			}
			Log.level(20) << "\r" << names[file] << " : " << mol_no+1 << flush;
		}
		Log.level(20) << endl;
		Log.flush();

		input->close();
		delete input;
	}

	if (compounds.size() < best_k)
	{
		Log.level(20) << "found " << compounds.size() << " compounds matching the given criteria." << endl;
	}

	GenericMolFile* output = MolFileFactory::open(output_file, ios::out, "mol2.gz");

	for (multimap < double, Molecule* > ::iterator it = compounds.begin();
		it!=compounds.end(); it++)
	{
		*output << *it->second;
		delete it->second;
	}

	output->close();
	delete output;
}

static PyObject *contacts(PyObject *self, PyObject *args) {
  
  int mol1, frame1, mol2, frame2;
  PyObject *selected1, *selected2;
  float cutoff;
  if (!PyArg_ParseTuple(args, (char *)"iiO!iiO!f:atomselection.contacts",
        &mol1, &frame1, &PyTuple_Type, &selected1,
        &mol2, &frame2, &PyTuple_Type, &selected2,
        &cutoff))
    return NULL;
  VMDApp *app = get_vmdapp();
  AtomSel *sel1 = sel_from_py(mol1, frame1, selected1, app);
  AtomSel *sel2 = sel_from_py(mol2, frame2, selected2, app);
  if (!sel1 || !sel2) {
    delete sel1;
    delete sel2;
    return NULL;
  }
  const float *ts1 = sel1->coordinates(app->moleculeList);
  const float *ts2 = sel2->coordinates(app->moleculeList);
  if (!ts1 || !ts2) {
    PyErr_SetString(PyExc_ValueError, "No coordinates in selection");
    delete sel1;
    delete sel2;
    return NULL;
  }
  Molecule *mol = app->moleculeList->mol_from_id(mol1);

  GridSearchPair *pairlist = vmd_gridsearch3(
      ts1, sel1->num_atoms, sel1->on,
      ts2, sel2->num_atoms, sel2->on,
      cutoff, -1, (sel1->num_atoms + sel2->num_atoms) * 27);

  delete sel1;
  delete sel2;
  GridSearchPair *p, *tmp;
  PyObject *list1 = PyList_New(0);
  PyObject *list2 = PyList_New(0);
  PyObject *tmp1;
  PyObject *tmp2;
  for (p=pairlist; p != NULL; p=tmp) {
    // throw out pairs that are already bonded
    MolAtom *a1 = mol->atom(p->ind1);
    if (mol1 != mol2 || !a1->bonded(p->ind2)) {
      // Needed to avoid a memory leak. Append increments the refcount 
      // of whatever gets added to it, but so does PyInt_FromLong.
      // Without a decref, the integers created never have their refcount
      //  go to zero, and you leak memory.
      tmp1 = PyInt_FromLong(p->ind1);
      tmp2 = PyInt_FromLong(p->ind2);
      PyList_Append(list1, tmp1);
      PyList_Append(list2, tmp2);
      Py_DECREF(tmp1);
      Py_DECREF(tmp2);
    }
    tmp = p->next;
    free(p);
  }
  PyObject *result = PyList_New(2);
  PyList_SET_ITEM(result, 0, list1);
  PyList_SET_ITEM(result, 1, list2);
  return result;
}

using std::ios;

Molecule* b;

CHECK(Molecule() throw())
	b = new Molecule;
	TEST_NOT_EQUAL(b, 0)
RESULT											

CHECK(~Molecule() throw())
	delete b;
RESULT

CHECK(Molecule(const Molecule& molecule, bool deep = true) throw())
	Atom a1;
	Molecule m("a"), m2;
	m.append(a1);
	m2 = Molecule(m);
	TEST_EQUAL(m2.getName(), "a")
	TEST_EQUAL(m2.countAtoms(), 1)
RESULT

CHECK(Molecule(const String& name) throw())
	Molecule m("a");
	TEST_EQUAL(m.getName(), "a")
	Molecule m2("");
	TEST_EQUAL(m2.getName(), "")
RESULT

CHECK([EXTRA] clear())
	System s;

void CanonicalSmilesSaver::saveMolecule (Molecule &mol_) const
{
   if (mol_.vertexCount() < 1)
      return;

   QS_DEF(Array<int>, ignored);
   QS_DEF(Array<int>, order);
   QS_DEF(Array<int>, ranks);
   QS_DEF(Molecule, mol);
   int i;

   if (mol_.repeating_units.size() > 0)
      throw Error("can not canonicalize a polymer");

   // Detect hydrogens configuration if aromatic but not ambiguous
   bool found_invalid_h = false;
   for (i = mol_.vertexBegin(); i != mol_.vertexEnd(); i = mol_.vertexNext(i))
   {
      if (mol_.isRSite(i) || mol_.isPseudoAtom(i))
         continue;

      if (mol_.getImplicitH_NoThrow(i, -1) == -1)
         found_invalid_h = true;
   }
   if (found_invalid_h)
   {
      AromaticityOptions options;
      options.method = AromaticityOptions::GENERIC;
      options.unique_dearomatization = true;
      MoleculeDearomatizer::restoreHydrogens(mol_, options);
   }

   mol.clone(mol_, 0, 0);

   // TODO: canonicalize allenes properly
   mol.allene_stereo.clear();

   ignored.clear_resize(mol.vertexEnd());
   ignored.zerofill();

   for (i = mol.vertexBegin(); i < mol.vertexEnd(); i = mol.vertexNext(i))
      if (mol.convertableToImplicitHydrogen(i))
         ignored[i] = 1;

   for (i = mol.edgeBegin(); i != mol.edgeEnd(); i = mol.edgeNext(i))
      if (mol.getBondTopology(i) == TOPOLOGY_RING && mol.cis_trans.getParity(i) != 0)
      {
         // we save cis/trans ring bonds into SMILES, but only those who
         // do not participate in bigger ring systems
         const Edge &edge = mol.getEdge(i);

         if (mol.getAtomRingBondsCount(edge.beg) != 2 ||
             mol.getAtomRingBondsCount(edge.end) != 2)
         {
            mol.cis_trans.setParity(i, 0);
            continue;
         }

         // also, discard the cis-trans bonds that have been converted to aromatic
         const Vertex &beg = mol.getVertex(edge.beg);
         const Vertex &end = mol.getVertex(edge.end);
         bool have_singlebond_beg = false;
         bool have_singlebond_end = false;
         int j;
         
         for (j = beg.neiBegin(); j != beg.neiEnd(); j = beg.neiNext(j))
            if (mol.getBondOrder(beg.neiEdge(j)) == BOND_SINGLE)
               have_singlebond_beg = true;

         for (j = end.neiBegin(); j != end.neiEnd(); j = end.neiNext(j))
            if (mol.getBondOrder(end.neiEdge(j)) == BOND_SINGLE)
               have_singlebond_end = true;

         if (!have_singlebond_beg || !have_singlebond_end)
         {
            mol.cis_trans.setParity(i, 0);
            continue;
         }
      }
         
   MoleculeAutomorphismSearch of;

   of.detect_invalid_cistrans_bonds = find_invalid_stereo;
   of.detect_invalid_stereocenters = find_invalid_stereo;
   of.find_canonical_ordering = true;
   of.ignored_vertices = ignored.ptr();
   of.process(mol);
   of.getCanonicalNumbering(order);

   for (i = mol.edgeBegin(); i != mol.edgeEnd(); i = mol.edgeNext(i))
      if (mol.cis_trans.getParity(i) != 0 && of.invalidCisTransBond(i))
         mol.cis_trans.setParity(i, 0);

   for (i = mol.vertexBegin(); i != mol.vertexEnd(); i = mol.vertexNext(i))
      if (mol.stereocenters.getType(i) > MoleculeStereocenters::ATOM_ANY && of.invalidStereocenter(i))
         mol.stereocenters.remove(i);

   ranks.clear_resize(mol.vertexEnd());

   for (i = 0; i < order.size(); i++)
//.........这里部分代码省略.........