C++ AtomIterator类代码示例

   /*
   * Zero forces on all local atoms and optionally on ghosts.
   */
   void AtomStorage::zeroForces(bool zeroGhosts)
   {
      int factor = 2;

      // Zero forces for all local atoms
      if (nAtom() > atomCapacity_/factor) {
         atoms_.zeroForces();
         // Optimization to allow sequential access
      } else {
         AtomIterator atomIter;
         begin(atomIter);
         for( ; atomIter.notEnd(); ++atomIter){
            atomIter->force().zero();
         }
      }
   
      // If using reverse communication, zero ghost atoms
      if (zeroGhosts && nGhost()) {
         if (nGhost() > ghostCapacity_/factor) {
            ghosts_.zeroForces();
            // Optimization to allow sequential access
         } else {
            GhostIterator ghostIter;
            begin(ghostIter);
            for( ; ghostIter.notEnd(); ++ghostIter){
               ghostIter->force().zero();
            }
         }
      }

   }

void BondGeometry::calculateCurrentBondGeometry(const Structure& structure) {
    int bondIndex = 0;
    for (AtomIterator iterator = structure.getAtomIterator();
            !iterator.isDone(); iterator.next()) {
        Atom atom = iterator.getCurrentAtom();
        setBondsForAnAtom(atom, bondIndex);
    }
}

	void IMGDock::reassignForces()
	{
		int i = 0;
		for (AtomIterator it = ligand_->beginAtom(); !it.isEnd(); it++)
		{
			it->setForce(old_forces_[i]);
			i++;
		}
	}

	void IMGDock::storeForces()
	{
		int i = 0;
		for (AtomIterator it = ligand_->beginAtom(); !it.isEnd(); it++)
		{
			old_forces_[i] = it->getForce();
			i++;
		}
	}

	void IMGDock::resetRotations()
	{
		current_conformation_ = vector < int > (current_conformation_.size(), 0);
		int atom = 0;
		for (AtomIterator it = ligand_->beginAtom(); +it; it++, atom++)
		{
			it->setPosition(original_atom_positions_[atom]);
		}
	}

	void IMGDock::saveAtomPositions()
	{
		original_atom_positions_.resize(scoring_function_->getNoLigandAtoms());
		Size atom = 0;
		for (AtomIterator it = ligand_->beginAtom(); +it; it++, atom++)
		{
			original_atom_positions_[atom] = it->getPosition();
		}
	}

	void IMGDock::translateLigand(Vector3& v)
	{
		TMatrix4x4<float> M;
		M.setTranslation(v);

		// transform all atoms of the ligand
		for (AtomIterator it = ligand_->beginAtom(); it != ligand_->endAtom(); it++)
		{
			it->setPosition(M*it->getPosition());
		}
	}

	void GridAnalysis::moveProbeGroup_(const Vector3& destination)
	{
		TMatrix4x4<float> M;
		Vector3 translation_vector = destination-center_;
		M.setTranslation(translation_vector);
		for (AtomIterator it = probe_group_.beginAtom(); +it; it++)
		{
			it->setPosition(M*it->getPosition());
		}
		center_ = destination;
	}

	void GridAnalysis::rotateProbeGroup_(int axis, int degree)
	{
		TMatrix4x4<float> M;
		TAngle<float> angle(degree, false);
		Vector3 axis_vector(axis == 0, axis == 1, axis == 2);
		M.setRotation(angle, axis_vector);

		for (AtomIterator it = probe_group_.beginAtom(); it != probe_group_.endAtom(); it++)
		{
			it->setPosition(M*(it->getPosition()-center_)+center_);
		}
	}

	void IMGDock::rotateLigand(int a, int degree)
	{
		TMatrix4x4<float> M;
		TAngle<float> angle(degree, false);
		Vector3 axis(a == 0, a == 1, a == 2);
		M.setRotation(angle, axis);
		const Vector3& origin = scoring_function_->getLigandCenter();

		for (AtomIterator it = ligand_->beginAtom(); it != ligand_->endAtom(); it++)
		{
			it->setPosition(M*(it->getPosition()-origin)+origin);
		}
	}

bool compareMolecules(Molecule& mol1, Molecule& mol2)
{
	if(mol1.countAtoms()!=mol2.countAtoms())
		return false;
	if(mol1.countBonds()!=mol2.countBonds())
		return false;

	AtomIterator ai;
	vector<Vector3> pos1;
	vector<float> q1;
	BALL_FOREACH_ATOM(mol1, ai)
	{
		pos1.push_back(ai->getPosition());
		q1.push_back(ai->getCharge());
	}

   /*
   * Record current positions of all local atoms and lock storage.
   */
   void AtomStorage::makeSnapshot()
   {
      // Precondition
      if (!isCartesian()) {
         UTIL_THROW("Error: Coordinates not Cartesian in makeSnapshot");
      }
 
      AtomIterator iter;
      int i = 0;
      for (begin(iter); iter.notEnd(); ++iter) {
         snapshot_[i] = iter->position();
         ++i;
      }
      locked_ = true;
   }

	void DockingAlgorithm::mapLigandOntoReferenceLigand()
	{
		if (!scoring_function_)
		{
			Log.error()<<"Error DockingAlgorithm::mapLigandOntoReferenceLigand() : ScoringFunction not set!"<<endl;
			return;
		}

		AtomContainer* ligand = scoring_function_->getLigand();

		if (!ligand)
		{
			Log.error()<<"Error DockingAlgorithm::mapLigandOntoReferenceLigand() : Ligand not set!"<<endl;
			return;
		}
		if (!reference_ligand_)
		{
			Log.error()<<"Error DockingAlgorithm::mapLigandOntoReferenceLigand() : Reference ligand not set!"<<endl;
			return;
		}

		Timer timer;
		timer.start();

		double lower_bound = 2;
		double upper_bound = 5;
		double tolerance = 1;
		Size no_matched_atoms = 0;
		double rmsd = 0;

		Matrix4x4 T = mapCompounds(*ligand, *reference_ligand_, no_matched_atoms, rmsd, upper_bound, lower_bound, tolerance);

		for (AtomIterator it = ligand->beginAtom(); +it; it++)
		{
			it->setPosition(T*it->getPosition());
		}

		timer.stop();
		Log.level(10)<<"superposing ligand: "<<timer.getClockTime()<<" seconds"<<endl;
	}

   /*
   *  First half of two-step velocity-Verlet integrator. 
   */
   void NphIntegrator::integrateStep1()
   {
      Vector dv;
      AtomIterator atomIter;

      Simulation& sys = simulation();
      sys.computeVirialStress();
      sys.computeKineticStress();
      sys.computeKineticEnergy();

      if (sys.domain().isMaster()) {
         P_target_ = simulation().boundaryEnsemble().pressure();
         T_kinetic_ = sys.kineticEnergy()*2.0/ndof_;
         Tensor stress = sys.virialStress();
         stress += sys.kineticStress();

         P_curr_diag_ = Vector(stress(0, 0), stress(1,1), stress(2,2));
         double P_curr = (1.0/3.0)*stress.trace();

         double mtk_term = (1.0/2.0)*dt_*T_kinetic_/W_;

         // Advance barostat
         double V = sys.boundary().volume();
         if (mode_ == Cubic) {
            nu_[0] += (1.0/2.0)*dt_*V/W_*(P_curr - P_target_) + mtk_term;
            nu_[1] = nu_[2] = nu_[0];
         } else if (mode_ == Tetragonal) {
            nu_[0] += (1.0/2.0)*dt_*V/W_*(P_curr_diag_[0] - P_target_) + mtk_term;
            nu_[1] += (1.0/2.0)*dt_*V/W_*((1.0/2.0)*(P_curr_diag_[1]+P_curr_diag_[2]) - P_target_) + mtk_term;
            nu_[2] = nu_[1];
         } else if (mode_  == Orthorhombic) {
            nu_[0] += (1.0/2.0)*dt_*V/W_*(P_curr_diag_[0] - P_target_) + mtk_term;
            nu_[1] += (1.0/2.0)*dt_*V/W_*(P_curr_diag_[1] - P_target_) + mtk_term;
            nu_[2] += (1.0/2.0)*dt_*V/W_*(P_curr_diag_[2] - P_target_) + mtk_term;
         }

      }

      #ifdef UTIL_MPI
      bcast(domain().communicator(), nu_,0);
      #endif

      // Precompute loop invariant quantities
      mtk_term_2_ = (nu_[0]+nu_[1]+nu_[2])/ndof_;
      Vector v_fac = Vector((1.0/4.0)*(nu_[0]+mtk_term_2_),
                            (1.0/4.0)*(nu_[1]+mtk_term_2_),
                            (1.0/4.0)*(nu_[2]+mtk_term_2_));
      exp_v_fac_ = Vector(exp(-v_fac[0]*dt_),
                          exp(-v_fac[1]*dt_),
                          exp(-v_fac[2]*dt_));
      Vector exp_v_fac_2 = Vector(exp(-(2.0*v_fac[0])*dt_),
                                  exp(-(2.0*v_fac[1])*dt_),
                                  exp(-(2.0*v_fac[2])*dt_));
      Vector r_fac = Vector((1.0/2.0)*nu_[0],
                            (1.0/2.0)*nu_[1],
                            (1.0/2.0)*nu_[2]);
      Vector exp_r_fac = Vector(exp(r_fac[0]*dt_),
                                exp(r_fac[1]*dt_),
                                exp(r_fac[2]*dt_));

      // Coefficients of sinh(x)/x = a_0 + a_2*x^2 + a_4*x^4 + a_6*x^6 + a_8*x^8 + a_10*x^10
      const double a[] = {double(1.0), double(1.0/6.0), double(1.0/120.0), 
                          double(1.0/5040.0), double(1.0/362880.0), double(1.0/39916800.0)};

      Vector arg_v = Vector(v_fac[0]*dt_, v_fac[1]*dt_, v_fac[2]*dt_);
      Vector arg_r = Vector(r_fac[0]*dt_, r_fac[1]*dt_, r_fac[2]*dt_);

      sinhx_fac_v_ = Vector(0.0,0.0,0.0);
      Vector sinhx_fac_r = Vector(0.0,0.0,0.0);

      Vector term_v = Vector(1.0,1.0,1.0);
      Vector term_r = Vector(1.0,1.0,1.0);

      for (unsigned int i = 0; i < 6; i++) {
         sinhx_fac_v_ += Vector(a[0]*term_v[0],
                                a[1]*term_v[1],
                                a[2]*term_v[2]);
         sinhx_fac_r += Vector(a[0]*term_r[0],
                               a[1]*term_r[1],
                               a[2]*term_r[2]);
         term_v = Vector(term_v[0] * arg_v[0] * arg_v[0],
                         term_v[1] * arg_v[1] * arg_v[1],
                         term_v[2] * arg_v[2] * arg_v[2]);
         term_r = Vector(term_r[0] * arg_r[0] * arg_r[0],
                         term_r[1] * arg_r[1] * arg_r[1],
                         term_r[2] * arg_r[2] * arg_r[2]);
      }

      // 1st half of NPH
      Vector vtmp;
      double prefactor; // = 0.5*dt/mass
      atomStorage().begin(atomIter);
      for ( ; atomIter.notEnd(); ++atomIter) {
         prefactor = prefactors_[atomIter->typeId()];

         dv.multiply(atomIter->force(), prefactor);
         dv[0] = dv[0] * exp_v_fac_[0]*sinhx_fac_v_[0];
//.........这里部分代码省略.........

		void EditMode::addStructure_()
		{
			QAction* action = dynamic_cast<QAction*>(sender());
			if (action == 0) return;

			String name = ascii(action->text());

			scene_->deselect();

			if (!fragment_db_)
			{
				fragment_db_ = new FragmentDB("fragments/Editing-Fragments.db");
			}

			list<AtomContainer*> containers = scene_->getContainers();
			if (containers.size() == 0) return;

			Residue* residue = fragment_db_->getResidueCopy(name);
			if (residue == 0)
			{
				residue = fragment_db_->getResidueCopy(name + "-Skeleton");
				if (residue == 0) return;
			}

			Vector3 p;
			Size nr = 0;
			AtomIterator ait = residue->beginAtom();
			for (;+ait; ++ait)
			{
				p += ait->getPosition();
				nr++;
			}

			if (nr == 0)
			{
				BALLVIEW_DEBUG
					delete residue;
				return;
			}

			p /= (float) nr;

			QPoint pos = scene_->mapFromGlobal(mouse_pos_new_);

			Matrix4x4 m;
			Vector3 x = scene_->mapViewportTo3D(pos.x(), pos.y());
			TransformationProcessor tf;

			Vector3 vv = scene_->getStage()->getCamera().getViewVector();
			float l = vv.getLength();
			if (!Maths::isZero(l)) vv /= l;
			Vector3 axis = Vector3(1,0,0) % vv;
			if (axis.getSquareLength() != 0)
			{
				Angle a = vv.getAngle(Vector3(1,0,0));
				m.setRotation(a, axis);
				tf.setTransformation(m);
				residue->apply(tf);
			}

			m.setTranslation(x - p);
			tf.setTransformation(m);
			residue->apply(tf);

			AtomContainer* s = *containers.begin();
			if (RTTI::isKindOf<System>(*s))
			{
				System* system = (System*) s;
				Molecule* mol = system->getMolecule(0);
				if (mol == 0)
				{
					mol = new Molecule();
					system->insert(*mol);
				}
				s = mol;
			}

			s->insert(*residue);
			scene_->getMainControl()->deselectCompositeRecursive(s, true);
			scene_->getMainControl()->selectCompositeRecursive(residue, true);
			scene_->getMainControl()->update(*s);
			scene_->notify(new NewSelectionMessage);
//			setMode(MOVE__MODE);
		}