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

	void UCK::getGraph(vector<String>& v, PairVector& e, const Molecule& mol)
	{
		weight_ = 0.0;
		Size count = 0;
		vector<pair<String, Size> >* mol_name;
		mol_name = new vector<pair<String, Size> >;
		bool found_atom = false;

		for(AtomConstIterator atit1 = mol.beginAtom(); atit1 != mol.endAtom(); ++atit1)
		{
			if(ignore_hydrogens_ && atit1->getElement()==PTE[1]) continue;

			// find chemical formula
			for(Size i = 0; i != mol_name->size(); ++i)
			{
				if((*mol_name)[i].first == atit1->getElement().getSymbol())	// increase number of already existing molecules
				{
					(*mol_name)[i].second++;
					found_atom = true;
					break;
				}
			}
			
			if(!found_atom)	// add current atom to formula, if it doesn't exist
			{
				mol_name->push_back(make_pair(atit1->getElement().getSymbol(),1));
				found_atom = false;
			}
			found_atom = false;
			
			weight_ += atit1->getElement().getAtomicWeight();
			v.push_back(atit1->getElement().getSymbol());	// add atom-name to label-list
			Size dest = 0;
			// find bonds from current atom to all other atoms and store them in e
			for(AtomConstIterator atit2 = mol.beginAtom(); atit2 != mol.endAtom(); ++atit2)
			{
				if(ignore_hydrogens_ && atit2->getElement()==PTE[1]) continue;

				if(atit1->getBond(*atit2) != 0)
				{
					e.push_back(make_pair(count, dest));
				}
				++dest;
			}
			++count;
		}

		sort(mol_name->begin(), mol_name->end());		// sort vector mol_name in order to get the lexicographically ordered
		for(Size i = 0; i != mol_name->size(); ++i)	// chemical formula
		{
			formula_ += ((*mol_name)[i].first)+(String)(*mol_name)[i].second;
		}
			
		delete mol_name;
		return;
	}

		Molecule* RGroupAssembler::generateNextMolecule()
	{
		if (!not_finished_) return NULL;
			Molecule* mol = new Molecule;
		System* scaffold_copy = new System;
		if (scaffolds_[scaffold_counter_] != NULL)
		{
			*scaffold_copy = *scaffolds_[scaffold_counter_];
		}
		else
		{
			return 0;
		}
		mol->insert(*scaffold_copy);
			for (AtomIterator atom_it = mol->beginAtom(); +atom_it; )
		{
			String atom_name = atom_it->getName();
			if (isPlacemark_(atom_name))
			{
				atom_it++;
				appendMoiety_(mol, moieties_[atom_name][moiety_counter_[atom_name]], atom_name);
			}
			else
			{
				atom_it++;
			}
		}
			not_finished_ = incrementCounter_();
		return mol;
	}

  DistanceFinder::DistanceFinder(SimInfo* info) : info_(info) {
    nObjects_.push_back(info_->getNGlobalAtoms()+info_->getNGlobalRigidBodies());
    nObjects_.push_back(info_->getNGlobalBonds());
    nObjects_.push_back(info_->getNGlobalBends());
    nObjects_.push_back(info_->getNGlobalTorsions());
    nObjects_.push_back(info_->getNGlobalInversions());

    stuntdoubles_.resize(nObjects_[STUNTDOUBLE]);
    bonds_.resize(nObjects_[BOND]);
    bends_.resize(nObjects_[BEND]);
    torsions_.resize(nObjects_[TORSION]);
    inversions_.resize(nObjects_[INVERSION]);
    
    SimInfo::MoleculeIterator mi;
    Molecule::AtomIterator ai;
    Molecule::RigidBodyIterator rbIter;
    Molecule::BondIterator bondIter;
    Molecule::BendIterator bendIter;
    Molecule::TorsionIterator torsionIter;
    Molecule::InversionIterator inversionIter;

    Molecule* mol;
    Atom* atom;
    RigidBody* rb;
    Bond* bond;
    Bend* bend;
    Torsion* torsion;
    Inversion* inversion;    
    
    for (mol = info_->beginMolecule(mi); mol != NULL; 
         mol = info_->nextMolecule(mi)) {
        
      for(atom = mol->beginAtom(ai); atom != NULL; 
          atom = mol->nextAtom(ai)) {
	stuntdoubles_[atom->getGlobalIndex()] = atom;
      }
      for (rb = mol->beginRigidBody(rbIter); rb != NULL; 
           rb = mol->nextRigidBody(rbIter)) {
	stuntdoubles_[rb->getGlobalIndex()] = rb;
      }
      for (bond = mol->beginBond(bondIter); bond != NULL; 
           bond = mol->nextBond(bondIter)) {
        bonds_[bond->getGlobalIndex()] = bond;
      }   
      for (bend = mol->beginBend(bendIter); bend != NULL; 
           bend = mol->nextBend(bendIter)) {
        bends_[bend->getGlobalIndex()] = bend;
      }   
      for (torsion = mol->beginTorsion(torsionIter); torsion != NULL; 
           torsion = mol->nextTorsion(torsionIter)) {
        torsions_[torsion->getGlobalIndex()] = torsion;
      }   
      for (inversion = mol->beginInversion(inversionIter); inversion != NULL; 
           inversion = mol->nextInversion(inversionIter)) {
        inversions_[inversion->getGlobalIndex()] = inversion;
      }   

    }
  }

  void IndexFinder::init() {

    SimInfo::MoleculeIterator mi;
    Molecule::AtomIterator ai;
    Molecule::RigidBodyIterator rbIter;
    Molecule::BondIterator bondIter;
    Molecule::BendIterator bendIter;
    Molecule::TorsionIterator torsionIter;
    Molecule::InversionIterator inversionIter;

    Molecule* mol;
    Atom* atom;
    RigidBody* rb;
    Bond* bond;
    Bend* bend;
    Torsion* torsion;
    Inversion* inversion;    
    
    for (mol = info_->beginMolecule(mi); mol != NULL; 
         mol = info_->nextMolecule(mi)) {
           
      SelectionSet ss(nObjects_);
      ss.bitsets_[MOLECULE].setBitOn(mol->getGlobalIndex());

      for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
	ss.bitsets_[STUNTDOUBLE].setBitOn(atom->getGlobalIndex());
      }
      for (rb = mol->beginRigidBody(rbIter); rb != NULL; 
           rb = mol->nextRigidBody(rbIter)) {
        ss.bitsets_[STUNTDOUBLE].setBitOn(rb->getGlobalIndex());
      }
      for (bond = mol->beginBond(bondIter); bond != NULL; 
           bond = mol->nextBond(bondIter)) {
        ss.bitsets_[BOND].setBitOn(bond->getGlobalIndex());
      }   
      for (bend = mol->beginBend(bendIter); bend != NULL; 
           bend = mol->nextBend(bendIter)) {
        ss.bitsets_[BEND].setBitOn(bend->getGlobalIndex());
      }   
      for (torsion = mol->beginTorsion(torsionIter); torsion != NULL; 
           torsion = mol->nextTorsion(torsionIter)) {
        ss.bitsets_[TORSION].setBitOn(torsion->getGlobalIndex());
      }   
      for (inversion = mol->beginInversion(inversionIter); inversion != NULL; 
           inversion = mol->nextInversion(inversionIter)) {
        ss.bitsets_[INVERSION].setBitOn(inversion->getGlobalIndex());
      }   

      selectionSets_[mol->getGlobalIndex()] = ss;
    }
  }

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;

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

  void BondAngleDistribution::process() {
    Molecule* mol;
    Atom* atom;
    RigidBody* rb;
    int myIndex;
    SimInfo::MoleculeIterator mi;
    Molecule::RigidBodyIterator rbIter;
    Molecule::AtomIterator ai;
    StuntDouble* sd;
    Vector3d vec;
    std::vector<Vector3d> bondvec;
    RealType r;    
    int nBonds;    
    int i;

    bool usePeriodicBoundaryConditions_ = info_->getSimParams()->getUsePeriodicBoundaryConditions();
    
    DumpReader reader(info_, dumpFilename_);    
    int nFrames = reader.getNFrames();
    frameCounter_ = 0;
    
    nTotBonds_ = 0;
    
    for (int istep = 0; istep < nFrames; istep += step_) {
      reader.readFrame(istep);
      frameCounter_++;
      currentSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot();
      
      if (evaluator_.isDynamic()) {
        seleMan_.setSelectionSet(evaluator_.evaluate());
      }

      // update the positions of atoms which belong to the rigidbodies
      
      for (mol = info_->beginMolecule(mi); mol != NULL; 
           mol = info_->nextMolecule(mi)) {
        for (rb = mol->beginRigidBody(rbIter); rb != NULL; 
             rb = mol->nextRigidBody(rbIter)) {
          rb->updateAtoms();
        }        
      }           
            
      // outer loop is over the selected StuntDoubles:

      for (sd = seleMan_.beginSelected(i); sd != NULL; 
           sd = seleMan_.nextSelected(i)) {

        myIndex = sd->getGlobalIndex();
        nBonds = 0;
        bondvec.clear();
        
        // inner loop is over all other atoms in the system:
        
        for (mol = info_->beginMolecule(mi); mol != NULL; 
             mol = info_->nextMolecule(mi)) {
          for (atom = mol->beginAtom(ai); atom != NULL; 
               atom = mol->nextAtom(ai)) {

            if (atom->getGlobalIndex() != myIndex) {

              vec = sd->getPos() - atom->getPos();       

              if (usePeriodicBoundaryConditions_) 
                currentSnapshot_->wrapVector(vec);
              
              // Calculate "bonds" and make a pair list 
              
              r = vec.length();
              
              // Check to see if neighbor is in bond cutoff 
              
              if (r < rCut_) { 
                // Add neighbor to bond list's
                bondvec.push_back(vec);
                nBonds++;
                nTotBonds_++;
              }  
            }
          }
          
          
          for (int i = 0; i < nBonds-1; i++ ){
            Vector3d vec1 = bondvec[i];
            vec1.normalize();
            for(int j = i+1; j < nBonds; j++){
              Vector3d vec2 = bondvec[j];
              
              vec2.normalize();
	      
              RealType theta = acos(dot(vec1,vec2))*180.0/NumericConstant::PI;
              
              
              if (theta > 180.0){
                theta = 360.0 - theta;
              }
              int whichBin = int(theta/deltaTheta_);
              
              histogram_[whichBin] += 2;
            }
          }           
//.........这里部分代码省略.........

	bool KCFFile::write(const Molecule& molecule)
		throw(File::CannotWrite)
	{
		if (!isOpen() || getOpenMode() != std::ios::out)
		{
			throw File::CannotWrite(__FILE__, __LINE__, name_);
		}

		// An alias for simplicity's sake...
		std::ostream& os(getFileStream());
		
		// Write ENTRY block
		// number of blanks????  properties are not read, written??? Which ones are there?
		os << ENTRY_TAG << "      " << molecule.getName() << std::endl;
		
		static char buffer[BALL_MAX_LINE_LENGTH];

		// Write NODE block
		// How to create the KEGG atom types? How many blanks?
		// This is not specified in the KCF format description, so we use what we can
    // deduce from example files.
		// First line gets the NODE tag
		os << NODE_TAG << "      " << molecule.countAtoms() << "\n"; 
		Size count = 1;
		AtomConstIterator ai(molecule.beginAtom());
		std::map<const Atom*, Position> atom_to_index;
		for (; +ai; ++ai, ++count)
		{
			// Write the atom line.
			// Blanks????
			String type = ai->getTypeName();
			String comment;
			
			// Make sure the type is in the set of KEGG types????
			// Blanks?
			sprintf(buffer, "             %d %s %s %6.4f %6.4f %s\n", 
							count, type.c_str(), ai->getElement().getSymbol().c_str(), 
							ai->getPosition().x, ai->getPosition().y, comment.c_str());
			os << buffer;
			
			// Remember the index of the current atom to map atom
			// pointers back to indices for the EDGE section.
			atom_to_index[&*ai] = count;
		}
		
		// Write EDGE block. Walk over all bonds to do so.
		// Blanks????
		os << "EDGE    " << molecule.countBonds() << "\n";
		count = 1;
		for (ai = molecule.beginAtom(); +ai; ++ai)
		{
			for (Atom::BondConstIterator bi(ai->beginBond()); +bi; ++bi)
			{
				Position index1 = atom_to_index[bi->getFirstAtom()];
				Position index2 = atom_to_index[bi->getSecondAtom()];
				String comment;
		
				// Write every bond just once				
				if (bi->getFirstAtom() == &*ai)
				{
					sprintf(buffer, "          %4d %4d %4d %1d%s\n", 
									count, index1, index2, bi->getOrder(), comment.c_str());
					os << buffer;
					++count;
				}
			}
		}
		
		// Write the DELIMITER block
		os << DELIMITER_TAG << std::endl;
		
		return true;
	}

//.........这里部分代码省略.........
	else
	{
		configuration_string += "use_polar_radius_rules " + polar_radius_rule_file_name + "\n";

		tmp = path.find(polar_radius_rule_file_name);
		if (tmp != "")
		{
			polar_radius_rule_file_name = tmp;
		}
		if (verbosity > 0)
		{
			Log.info() << "Using " << polar_radius_rule_file_name << " as polar radius rule file" << endl;
		}

		INIFile radius_ini(polar_radius_rule_file_name);
		radius_ini.read();
		RadiusRuleProcessor radius_rules(radius_ini);
		system.apply(radius_rules);
	}

	if (verbosity > 8)
	{
		Log.info() << "system statistics:" << endl;
		Log.info() << "molecules: " << system.countMolecules() << endl;
		Log.info() << "proteins:  " << system.countProteins() << endl;
		Log.info() << "fragments: " << system.countFragments() << endl;
		Log.info() << "atoms:     " << system.countAtoms() << endl;
	}


	// check for uassigned type names
	AtomIterator it;
	int count = 0;
	for (it = system.beginAtom(); +it; ++it)
	{
		count++;
		if (it->getElement().getSymbol() == "?")
		{
			Log.info() << "Got symbol \"?\": " << it->getFullName() << endl;
		}
		if (it->getCharge() > 8.0)
		{
			Log.error() << "Mysterious charge: " << it->getCharge() << "\t"
				<< it->getFullName() << "\t" << count << "\t" << it->getElement().getSymbol() << endl;
		}
	}

	// scale charges
	if (verbosity > 0)
	{
		Log.info() << "Scaling receptor charges by " << receptor_charge_scaling_factor << endl;
	}
	if (parpars.has("s"))
	{
		for (it = protein->beginAtom(); +it; ++it)
		{
			it->setCharge(it->getCharge() * receptor_charge_scaling_factor);
		}
	}

	if (verbosity > 0)
	{
		Log.info() << "Scaling ligand charges by " << ligand_charge_scaling_factor << endl;
	}
	if (parpars.has("t"))
	{