//------------------------------------------------------------------------------ // This file is part of the OpenStructure project <www.openstructure.org> // // Copyright (C) 2008-2011 by the OpenStructure authors // // This library is free software; you can redistribute it and/or modify it under // the terms of the GNU Lesser General Public License as published by the Free // Software Foundation; either version 3.0 of the License, or (at your option) // any later version. // This library is distributed in the hope that it will be useful, but WITHOUT // ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more // details. // // You should have received a copy of the GNU Lesser General Public License // along with this library; if not, write to the Free Software Foundation, Inc., // 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA //------------------------------------------------------------------------------ #ifndef _MSC_VER #include <sys/time.h> #endif #include <map> #include <boost/shared_ptr.hpp> #include <boost/version.hpp> #include <boost/version.hpp> #if BOOST_VERSION>=103900 #include <boost/smart_ptr/make_shared.hpp> #define MAKE_SHARED_AVAILABLE 1 #else #define MAKE_SHARED_AVAILABLE 0 #endif #include <boost/logic/tribool.hpp> #include <ost/geom/geom.hh> #include <ost/log.hh> #include <ost/profile.hh> #include "atom_impl.hh" #include "entity_impl.hh" #include <ost/mol/entity_visitor.hh> #include <ost/mol/entity_observer.hh> #include "atom_impl.hh" #include "residue_impl.hh" #include "chain_impl.hh" #include "connector_impl.hh" #include "torsion_impl.hh" #include "fragment_impl.hh" #include "query_impl.hh" #include <ost/mol/mol.hh> #include <ost/mol/impl/residue_impl.hh> #include <ost/mol/bond_table.hh> #include <ost/mol/atom_handle.hh> #include <ost/mol/not_connected_error.hh> #include <ost/mol/query_state.hh> #include <ost/profile.hh> #include <ost/integrity_error.hh> #include "torsion_impl.hh" using boost::logic::tribool; using boost::logic::indeterminate; namespace ost { namespace mol { namespace impl { EntityImpl::EntityImpl(): atom_map_(), chain_list_(), connector_map_(), torsion_map_(), transformation_matrix_(), inverse_transformation_matrix_(), identity_transf_(true), atom_organizer_(5.0), fragment_list_(), observer_map_(), xcs_editor_count_(0), ics_editor_count_(0), dirty_flags_(DisableICS), name_(""), next_index_(0L), default_query_flags_(0) { } int EntityImpl::GetAtomCount() const { return static_cast<int>(atom_map_.size()); } int EntityImpl::GetBondCount() const { return static_cast<int>(connector_map_.size()); } mol::BondHandleList EntityImpl::GetBondList() const { mol::BondHandleList bhl; mol::impl::ConnectorImplMap::const_iterator i; for (i=connector_map_.begin();i!=connector_map_.end();++i) { bhl.push_back(mol::BondHandle((i->second))); } return bhl; } int EntityImpl::GetResidueCount() const { int count=0; for (ChainImplList::const_iterator i=chain_list_.begin(), e=chain_list_.end(); i!=e; ++i) { count+=(*i)->GetResidueCount(); } return count; } EntityImplPtr EntityImpl::Copy() { Profile prof("EntityImpl::Copy"); #if MAKE_SHARED_AVAILABLE EntityImplPtr ent_p=boost::make_shared<EntityImpl>(); #else EntityImplPtr ent_p(new EntityImpl); #endif this->DoCopy(ent_p); return ent_p; } ChainImplPtr EntityImpl::InsertChain(const ChainImplPtr& chain) { ChainImplPtr dst_chain=this->InsertChain(chain->GetName()); dst_chain->Assign(*chain.get()); return dst_chain; } void EntityImpl::ReplicateHierarchy(EntityImplPtr dest) { for (ChainImplList::const_iterator i=chain_list_.begin(), e1=chain_list_.end(); i!=e1; ++i) { ChainImplPtr src_chain=*i; ChainImplPtr dst_chain=dest->InsertChain(src_chain); // copy generic properties dst_chain->Assign(*src_chain.get()); for (ResidueImplList::iterator j=src_chain->GetResidueList().begin(), e2=src_chain->GetResidueList().end(); j!=e2; ++j) { ResidueImplPtr src_res=*j; ResidueImplPtr dst_res=dst_chain->AppendResidue(src_res); for (AtomImplList::iterator k=src_res->GetAtomList().begin(), e3=src_res->GetAtomList().end(); k!=e3; ++k) { AtomImplPtr src_atom=*k; AtomImplPtr dst_atom=dst_res->InsertAtom(src_atom); dst_atom->SetIndex(src_atom->GetIndex()); } } } } namespace { AtomImplPtr lookup_atom(const AtomImplPtr& atom, EntityImplPtr dest, int atom_index=-1) { if (atom_index==-1) { ResidueImplPtr res=atom->GetResidue(); AtomImplList::const_iterator i=std::find(res->GetAtomList().begin(), res->GetAtomList().end(), atom); if (i==res->GetAtomList().end()) return AtomImplPtr(); atom_index=i-res->GetAtomList().begin(); } int res_index=atom->GetResidue()->GetIndex(); ChainImplPtr chain=dest->FindChain(atom->GetResidue()->GetChain()->GetName()); return ((chain->GetResidueList())[res_index]->GetAtomList())[atom_index]; } } void EntityImpl::DoCopyBondsAndTorsions(EntityImplPtr dest) { for (ChainImplList::const_iterator i=chain_list_.begin(), i2=dest->chain_list_.begin(), e1=chain_list_.end(); i!=e1; ++i, ++i2) { ChainImplPtr src_chain=*i; ChainImplPtr dst_chain=*i2; for (ResidueImplList::iterator j=src_chain->GetResidueList().begin(), j2=dst_chain->GetResidueList().begin(), e2=src_chain->GetResidueList().end(); j!=e2; ++j, ++j2) { ResidueImplPtr src_res=*j, dst_res=*j2; int res_index=j-src_chain->GetResidueList().begin(); for (AtomImplList::iterator k=src_res->GetAtomList().begin(), k2=dst_res->GetAtomList().begin(), e3=src_res->GetAtomList().end(); k!=e3; ++k, ++k2) { AtomImplPtr src_atom=*k; AtomImplPtr dst_atom=*k2; const ConnectorImplList& sec_conn=src_atom->GetSecondaryConnectors(); for (ConnectorImplList::const_iterator l=sec_conn.begin(), e4=sec_conn.end(); l!=e4; ++l) { ConnectorImplP conn=*l; if (conn->GetFirst()==src_atom) { AtomImplPtr dst_atom2=lookup_atom(conn->GetSecond(), dest); assert(dst_atom2); dest->Connect(dst_atom, dst_atom2, 0, 0, 0, conn->GetBondOrder()); } } } const TorsionImplList& torsions=src_res->GetTorsionList(); for (TorsionImplList::const_iterator k=torsions.begin(), e3=torsions.end(); k!=e3; ++k) { TorsionImplP t=*k; // only add torsion if all involved atoms have a residue index // smaller/equal to our residue index if (t->GetFourth()->GetResidue()->GetIndex()>res_index || t->GetThird()->GetResidue()->GetIndex()>res_index || t->GetSecond()->GetResidue()->GetIndex()>res_index || t->GetFirst()->GetResidue()->GetIndex()>res_index) { continue; } dest->AddTorsion(t->GetName(), lookup_atom(t->GetFirst(), dest), lookup_atom(t->GetSecond(), dest), lookup_atom(t->GetThird(), dest), lookup_atom(t->GetFourth(), dest)); } } } dest->MarkTraceDirty(); } void EntityImpl::DoCopy(EntityImplPtr dest) { // first copy the chain - residue - atom hierarchy before replicating the // bond network and the torsions dest->SetName(this->GetName()); dest->next_index_=0; this->ReplicateHierarchy(dest); this->DoCopyBondsAndTorsions(dest); } int EntityImpl::GetChainCount() const { return static_cast<int>(chain_list_.size()); } const geom::Mat4& EntityImpl::GetTransfMatrix() const { return transformation_matrix_; } const geom::Mat4& EntityImpl::GetInvTransfMatrix() const { return inverse_transformation_matrix_; } bool EntityImpl::IsTransfIdentity() const { return identity_transf_; } EntityImpl::~EntityImpl() { // notify all observers of pending destruct EntityObserverMap::iterator it; for (it=observer_map_.begin();it!=observer_map_.end();++it) { it->second->OnDestroy(); } } geom::AlignedCuboid EntityImpl::GetBounds() const { if (this->GetAtomCount()>0) { geom::Vec3 mmin, mmax; AtomImplMap::const_iterator it=atom_map_.begin(); mmin=mmax=it->second->TransformedPos(); for (++it; it!=atom_map_.end();++it) { mmin=geom::Min(mmin,it->second->TransformedPos()); mmax=geom::Max(mmax,it->second->TransformedPos()); } return geom::AlignedCuboid(mmin, mmax); } else { return geom::AlignedCuboid(geom::Vec3(), geom::Vec3()); } } geom::Vec3 EntityImpl::GetCenterOfAtoms() const { geom::Vec3 center; if (this->GetAtomCount()>0) { for (AtomImplMap::const_iterator it = atom_map_.begin(); it!=atom_map_.end();++it) { center+=it->second->TransformedPos(); } center/=static_cast<Real>(atom_map_.size()); } return center; } geom::Vec3 EntityImpl::GetCenterOfMass() const { geom::Vec3 center; Real mass = this->GetMass(); if (this->GetAtomCount()>0 && mass>0) { for(AtomImplMap::const_iterator it = atom_map_.begin();it!=atom_map_.end();++it) { center+=it->second->TransformedPos()*it->second->GetMass(); } center/=mass; } return center; } Real EntityImpl::GetMass() const { double mass=0.0; for (ChainImplList::const_iterator i=chain_list_.begin(), e=chain_list_.end(); i!=e; ++i) { ChainImplPtr chain=*i; mass+=chain->GetMass(); } return mass; } AtomImplPtr EntityImpl::CreateAtom(const ResidueImplPtr& rp, const String& name, const geom::Vec3& pos, const String& ele) { #if MAKE_SHARED_AVAILABLE AtomImplPtr ap=boost::make_shared<AtomImpl>(shared_from_this(), rp, name, pos, ele,next_index_++); #else AtomImplPtr ap(new AtomImpl(shared_from_this(), rp, name, pos, ele, next_index_++)); #endif if (!identity_transf_) { geom::Vec3 transformed_pos = geom::Vec3(transformation_matrix_*geom::Vec4(pos)); ap->TransformedPos()=transformed_pos; atom_organizer_.Add(ap,transformed_pos); } else { ap->TransformedPos()=pos; atom_organizer_.Add(ap,pos); } atom_map_.insert(AtomImplMap::value_type(ap.get(),ap)); return ap; } void EntityImpl::DeleteAtom(const AtomImplPtr& atom) { atom_map_.erase(atom.get()); atom_organizer_.Remove(atom); } ResidueImplPtr EntityImpl::CreateResidue(const ChainImplPtr& cp, const ResNum& n, const ResidueKey& k) { #if MAKE_SHARED_AVAILABLE return boost::make_shared<ResidueImpl>(shared_from_this(), cp, n, k); #else return ResidueImplPtr(new ResidueImpl(shared_from_this(), cp, n, k)); #endif } ChainImplPtr EntityImpl::InsertChain(const String& cname) { for (ChainImplList::iterator i=chain_list_.begin(), e=chain_list_.end(); i!=e; ++i) { if ((*i)->GetName()==cname) { throw IntegrityError("Can't insert chain. A chain with name '"+cname+ "' already exists"); } } #if MAKE_SHARED_AVAILABLE ChainImplPtr cp=boost::make_shared<ChainImpl>(shared_from_this(), cname); #else ChainImplPtr cp(new ChainImpl(shared_from_this(), cname)); #endif chain_list_.push_back(cp); return cp; } ConnectorImplP EntityImpl::Connect(const AtomImplPtr& first, const AtomImplPtr& second, Real len, Real theta, Real phi, unsigned char bond_order) { if (!first || !second) return ConnectorImplP(); // check for already existing connection if(first->GetPrimaryConnector() && first->GetPrimaryConnector()->GetFirst()==second) { LOG_DEBUG("connect: returning existing bond (1)"); return first->GetPrimaryConnector(); } else if(second->GetPrimaryConnector() && second->GetPrimaryConnector()->GetFirst()==first) { LOG_DEBUG("connect: returning existing bond (2)"); return second->GetPrimaryConnector(); } else { ConnectorImplList clist = first->GetSecondaryConnectors(); for(ConnectorImplList::const_iterator it=clist.begin();it!=clist.end();++it) { if((*it)->GetSecond()==second) { LOG_DEBUG("connect: returning existing bond (3)"); return *it; } } clist = second->GetSecondaryConnectors(); for(ConnectorImplList::const_iterator it=clist.begin();it!=clist.end();++it) { if((*it)->GetSecond()==first) { LOG_DEBUG("connect: returning existing bond (4)"); return *it; } } } // no primary connectors are generated, this is left for the directionality // trace #if MAKE_SHARED_AVAILABLE ConnectorImplP bp=boost::make_shared<ConnectorImpl>(shared_from_this(), first, second, len, theta, phi, bond_order); assert(bp->GetFirst()); assert(bp->GetSecond()); #else ConnectorImplP bp(new ConnectorImpl(shared_from_this(), first, second, len, theta, phi, bond_order)); #endif first->AddSecondaryConnector(bp); second->AddSecondaryConnector(bp); LOG_DEBUG("adding new connector " << bp); connector_map_.insert(ConnectorImplMap::value_type(bp.get(),bp)); return bp; } TorsionImplP EntityImpl::AddTorsion(const String& name, const AtomImplPtr& a1, const AtomImplPtr& a2, const AtomImplPtr& a3, const AtomImplPtr& a4) { if (ConnectorExists(a2, a3) && ConnectorExists(a1, a2) && ConnectorExists(a3, a4)) { #if MAKE_SHARED_AVAILABLE TorsionImplP impl=boost::make_shared<TorsionImpl>(name, a1, a2, a3, a4); #else TorsionImplP impl(new TorsionImpl(name, a1, a2, a3, a4)); #endif torsion_map_.insert(TorsionImplMap::value_type(impl.get(), impl)); std::set<ResidueImpl*> involved_residues; involved_residues.insert(a1->GetResidue().get()); involved_residues.insert(a2->GetResidue().get()); involved_residues.insert(a3->GetResidue().get()); involved_residues.insert(a4->GetResidue().get()); std::set<ResidueImpl*>::const_iterator i=involved_residues.begin(); for (; i!=involved_residues.end(); ++i) { (*i)->AddTorsion(impl); } return impl; } return TorsionImplP(); } bool EntityImpl::SetAngle(const AtomImplPtr& a1, const AtomImplPtr& a2, const AtomImplPtr& a3, Real angle) { ConnectorImplP c12=GetConnector(a1, a2); ConnectorImplP c23=GetConnector(a2, a3); if (c12==c23) return false; geom::Vec3 v12, v23; if (c12 && c12->GetFirst() && c12->GetSecond() && c23 && c23->GetFirst() && c23->GetSecond()) { if (c12==a2->GetPrimaryConnector()) { v12=geom::Vec3(0, 0, -1); v23=c23->GetDir(); } else if (c23==a2->GetPrimaryConnector()) { v12=c12->GetDir(); v23=geom::Vec3( 0, 0, -1); } else { v12=c12->GetDir(); v23=c23->GetDir(); } geom::Vec3 normal; if (std::abs(Dot(v12, v23)) > 0.999) { //The two vectors are colinear. This means we have the freedom to choose //the normal we want to rotate around, as long as it is orthogonal to the // two vectors. normal=OrthogonalVector(v12); } else { normal=geom::Cross(v12, v23); } geom::Mat3 rot=AxisRotation(normal, angle); if (c23==a2->GetPrimaryConnector()) { assert(0 && "not implemented yet."); } else { c23->SetDir(rot*v12); } return true; } return false; } Real EntityImpl::GetAngleXCS(const AtomImplPtr& a1, const AtomImplPtr& a2, const AtomImplPtr& a3) const { ConnectorImplP c12=GetConnector(a1, a2); ConnectorImplP c23=GetConnector(a2, a3); if (c12 && c12->GetFirst() && c12->GetSecond()) { if (c23 && c23->GetFirst() && c23->GetSecond()) { return Angle(a2->TransformedPos()-a1->TransformedPos(), a2->TransformedPos()-a3->TransformedPos()); } else { AtomHandle ah2(a2), ah3(a3); throw NotConnectedError(ah2, ah3); } } else { AtomHandle ah1(a1), ah3(a3); throw NotConnectedError(ah1, ah3); } } Real EntityImpl::GetAngleICS(const AtomImplPtr& a1, const AtomImplPtr& a2, const AtomImplPtr& a3) const { ConnectorImplP c12=GetConnector(a1, a2); ConnectorImplP c23=GetConnector(a2, a3); if (c12 && c12->GetFirst() && c12->GetSecond()) { if (c23 && c23->GetFirst() && c23->GetSecond()) { // all that remains is to calculate the angle of the direction vectors. // in case one of the connector is a2's primary connector we set it's // direction vector to (0, 0, -1). if (c12==a2->GetPrimaryConnector()) { return Angle(geom::Vec3(0,0,-1), c23->GetDir()); } else if (c23==a2->GetPrimaryConnector()) { return Angle(c12->GetDir(), geom::Vec3(0,0,-1)); } else { // both are secondary connectors. return Angle(c12->GetDir(),c23->GetDir()); } } else { throw NotConnectedError(AtomHandle(a2), AtomHandle(a3)); } } else { throw NotConnectedError(AtomHandle(a1), AtomHandle(a2)); } } Real EntityImpl::GetAngle(const AtomImplPtr& a1, const AtomImplPtr& a2, const AtomImplPtr& a3) const { if (dirty_flags_ & DirtyXCS) { return this->GetAngleICS(a1, a2, a3); } else { return this->GetAngleXCS(a1, a2, a3); } } namespace { unsigned long ai_ptr_diff(const AtomImplPtr& aip1, const AtomImplPtr& aip2) { return std::abs(reinterpret_cast<long>(aip1.get())-reinterpret_cast<long>(aip2.get())); } } /* Trace and verify directionality, as encoded by the first and second atom in each connector. The hints from the original order as given by Connect are taken into account. */ void EntityImpl::TraceDirectionality() { fragment_list_.clear(); Profile profile_trace("trace directionality"); // some reseting on the atom level for(AtomImplMap::iterator it=atom_map_.begin();it!=atom_map_.end();++it) { it->second->ClearDirectionality(); it->second->SetVisited(false); it->second->SetTraced(false); } /* the beginning of a fragment is given by two conditions, namely that it does not point to a previous atom, and that none of the connectors have it as second. */ // mark all second atoms for(ConnectorImplMap::iterator it=connector_map_.begin(); it!=connector_map_.end(); ++it) { it->second->GetSecond()->SetVisited(true); } // search atoms that fulfill above criteria, and start recursive trace with them unsigned int traced_atom_count=0; unsigned int traced_atom_count_prev=atom_map_.size()+1; while (true) { for(AtomImplMap::iterator it=atom_map_.begin();it!=atom_map_.end();++it) { if (!it->second->IsVisited() && !it->second->HasPrevious()) { it->second->SetVisited(true); LOG_DEBUG("dir: new fragment"); #if MAKE_SHARED_AVAILABLE FragmentImplP frag=boost::make_shared<FragmentImpl>(it->second); #else FragmentImplP frag( new FragmentImpl(it->second)); #endif fragment_list_.push_back(frag); it->second->TraceDirectionality(frag,ConnectorImplP(), 0, traced_atom_count); } } if(traced_atom_count<atom_map_.size()) { LOG_DEBUG("entering closed loops search, since only " << traced_atom_count << " from " << atom_map_.size() << " atoms were traced"); /* identify closed loops that prohibit the trace from completing */ AtomImplPtr aip1; AtomImplPtr aip2; // step 1: find the connector with the largest atom order gap // (this is a heuristic guess for the best atom to isolate) for (ConnectorImplMap::iterator it=connector_map_.begin(); it!=connector_map_.end(); ++it) { if(!it->second->GetFirst()->HasPrevious() && !it->second->GetSecond()->HasPrevious() ) { // candidates for loop removal if(aip1) { if(ai_ptr_diff(it->second->GetFirst(),it->second->GetSecond())>ai_ptr_diff(aip1,aip2)) { aip1=it->second->GetFirst(); aip2=it->second->GetSecond(); } } else { aip1=it->second->GetFirst(); aip2=it->second->GetSecond(); } } } LOG_DEBUG("found pair [" << aip1 << "] [" << aip2 << "]"); // step 2: swap all connectors pointing to aip2 // thus removing its 'visited' state for(ConnectorImplMap::iterator it=connector_map_.begin(); it!=connector_map_.end(); ++it) { if(it->second->GetSecond()==aip2) { it->second->Switch(); it->second->GetSecond()->SetVisited(true); } } aip2->SetVisited(false); // and now back to the recursive trace } else { // finished trace break; } // infinite loop safeguard if(traced_atom_count == traced_atom_count_prev) { LOG_VERBOSE("Encountered unbreakable locked state during directionality trace"); break; } traced_atom_count_prev=traced_atom_count; } LOG_VERBOSE("Directionality trace completed with " << fragment_list_.size() << " fragment(s)"); } bool EntityImpl::IsXCSDirty() const { return (dirty_flags_ & DirtyXCS)!=0; } void EntityImpl::UpdateFromICS() { LOG_DEBUG("updating external from internal coordinate system"); Profile prof_conv("update external from internal coordinate system"); // traverse through all atoms for(AtomImplMap::iterator it=atom_map_.begin();it!=atom_map_.end();++it) { it->second->SetVisited(false); } for (FragmentImplList::iterator it=fragment_list_.begin(); it!=fragment_list_.end();++it) { LOG_DEBUG("entering root of tree for ICS update"); assert(*it); AtomImplPtr ap=(*it)->GetAtom(); // recursive update, using position of first atom /*geom::Mat4 pmat=(*it)->GetTMat(); pmat=tmat*pmat; ap->SetICSMat(pmat);*/ ap->UpdateFromICS(); LOG_DEBUG("done with tree traversal for ICS update"); } LOG_DEBUG("refreshing all connector positions"); // refresh all connector positions /* for (ConnectorImplMap::iterator it=connector_map_.begin(); it!=connector_map_.end();++it) { it->second->UpdatePosition(); }*/ this->UpdateTransformedPos(); } void EntityImpl::UpdateFromXCS() { LOG_DEBUG("rebuilding internal from external coordinate system"); Profile prof_conv("updating internal from external coordinate system"); if(fragment_list_.empty()) { // this should have been handled by conopology after the // initial connectivity assignments - if not, then do it here Profile profile_trace("finished directionality trace"); LOG_DEBUG("warning: no fragment list, re-running directionality trace"); TraceDirectionality(); } if(Logger::Instance().GetVerbosityLevel()>Logger::DEBUG) { LOG_TRACE("dumping directionality"); for(AtomImplMap::iterator it=atom_map_.begin();it!=atom_map_.end();++it) { LOG_TRACE(" " << it->second << ":"); ConnectorImplList clist = it->second->GetSecondaryConnectors(); for(ConnectorImplList::const_iterator cit=clist.begin();cit!=clist.end();++cit) { LOG_TRACE(" " << (*cit)->GetFirst() << " -> " << (*cit)->GetSecond()); } } } for(AtomImplMap::iterator it=atom_map_.begin();it!=atom_map_.end();++it) { it->second->SetVisited(false); } for (FragmentImplList::iterator it=fragment_list_.begin(); it!=fragment_list_.end();++it) { LOG_DEBUG("entering fragment tree traversal for XCS update"); FragmentImplP frag= *it; AtomImplPtr atom = frag->GetAtom(); { atom->UpdateFromXCS(); } } } void EntityImpl::Apply(EntityVisitor& v) { LOG_TRACE("visitor @" << &v << " visiting entity impl @" << this); v.OnEntry(); for(ChainImplList::iterator it = chain_list_.begin();it!=chain_list_.end();++it) { (*it)->Apply(v); } for(ConnectorImplMap::iterator it = connector_map_.begin();it!=connector_map_.end();++it) { it->second->Apply(v); } for(TorsionImplMap::iterator it = torsion_map_.begin();it!=torsion_map_.end();++it) { it->second->Apply(v); } v.OnExit(); } void EntityImpl::ApplyTransform(const geom::Mat4 transfmat) { geom::Mat4 updated_transformation_matrix=transfmat*transformation_matrix_; this->SetTransform(updated_transformation_matrix); } void EntityImpl::SetTransform(const geom::Mat4 transfmat) { transformation_matrix_=transfmat; inverse_transformation_matrix_=Invert(transformation_matrix_); identity_transf_ = (transformation_matrix_==geom::Mat4()); this->UpdateTransformedPos(); this->MarkOrganizerDirty(); } void EntityImpl::AttachObserver(const EntityObserverPtr& o) { observer_map_.insert(EntityObserverMap::value_type(o.get(),o)); } void EntityImpl::DetachObserver(const EntityObserverPtr& o) { EntityObserverMap::iterator it=observer_map_.find(o.get()); } void EntityImpl::NotifyObserver() { for(EntityObserverMap::iterator it=observer_map_.begin();it!=observer_map_.end();++it) { it->second->OnUpdate(); } } void EntityImpl::Swap(EntityImpl& impl) { atom_map_.swap(impl.atom_map_); chain_list_.swap(impl.chain_list_); connector_map_.swap(impl.connector_map_); torsion_map_.swap(impl.torsion_map_); std::swap(transformation_matrix_,impl.transformation_matrix_); atom_organizer_.Swap(impl.atom_organizer_); fragment_list_.swap(impl.fragment_list_); observer_map_.swap(impl.observer_map_); } TorsionImplP EntityImpl::FindTorsion(const AtomImplPtr& a1, const AtomImplPtr& a2, const AtomImplPtr& a3, const AtomImplPtr& a4) const { TorsionImplMap::const_iterator i; for (i=torsion_map_.begin(); i!=torsion_map_.end(); ++i) { TorsionImplP t = i->second; if (t->Matches(a1, a2, a3, a4)) { return t; } } return TorsionImplP(); } AtomImplList EntityImpl::FindWithin(const geom::Vec3& pos, Real radius) const { SpatialAtomOrganizer::ItemList alist = atom_organizer_.FindWithin(pos,radius); return alist; } namespace { typedef BondTable<AtomView> BondTableType; void update_bond_table(BondTableType& bond_table, const AtomImplPtr& ai, const AtomView& view) { bond_table.Update(BondHandle(ai->GetPrimaryConnector()), view); const ConnectorImplList& clist = ai->GetSecondaryConnectors(); ConnectorImplList::const_iterator it; for (it=clist.begin(); it!=clist.end(); ++it) { bond_table.Update(BondHandle(*it), view); } } } template<bool always_true> EntityView EntityImpl::do_selection(const EntityHandle& eh, const Query& query, QueryFlags flags) const { BondTableType bond_table; LOG_DEBUG("entering do selection"); size_t chain_count = 0, residue_count = 0, atom_count = 0; bool c_added, r_added; LOG_DEBUG("creating view"); EntityView view(eh); LOG_DEBUG("creating query state"); EntityHandle myself(const_cast<impl::EntityImpl*>(this)->shared_from_this()); QueryState qs(query.CreateQueryState(myself)); LOG_DEBUG("entering chain loop"); for (ChainImplList::const_iterator ch_it=chain_list_.begin(); ch_it!=chain_list_.end();++ch_it) { LOG_DEBUG("checking chain " << (*ch_it)->GetName()); c_added = false; tribool c = always_true ? tribool(true) : qs.EvalChain(*ch_it); if (c == true) { LOG_DEBUG("chain is selected"); // Include all residues const ChainImplPtr& ci=*ch_it; ++chain_count; ChainView chain=view.AddChain(ci); ResidueImplList::const_iterator re_it = ci->GetResidueList().begin(); for ( ;re_it!=ci->GetResidueList().end();++re_it) { LOG_DEBUG(" adding residue " << (*re_it)->GetNumber()); ResidueImplPtr& ri = const_cast<ResidueImplPtr&>(*re_it); ++residue_count; ResidueView res =chain.AddResidue(ri); AtomImplList::const_iterator at_it = ri->GetAtomList().begin(); for (;at_it != ri->GetAtomList().end(); ++at_it) { LOG_DEBUG(" adding atom " << (*at_it)->GetName()); ++atom_count; if (flags & QueryFlag::NO_BONDS) { res.AddAtom(*at_it); } else { update_bond_table(bond_table, *at_it, res.AddAtom(*at_it)); } } } } else if (indeterminate(c)) { // Test residues r_added = false; const ChainImplPtr& ci = *ch_it; ChainView chain; ResidueImplList::const_iterator re_it = ci->GetResidueList().begin(); for ( ;re_it!=ci->GetResidueList().end();++re_it) { LOG_DEBUG(" checking residue " << (*re_it)->GetNumber()); tribool r = qs.EvalResidue(*re_it); ResidueImplPtr& ri = const_cast<ResidueImplPtr&>(*re_it); if (r == true) { LOG_DEBUG(" residue is selected"); // Include all atoms ResidueImplPtr& ri = const_cast<ResidueImplPtr&>(*re_it); AtomImplList::const_iterator at_it = ri->GetAtomList().begin(); if (!c_added) { c_added=true; ++chain_count; chain=view.AddChain(ci); } ++residue_count; ResidueView res=chain.AddResidue(ri); for (;at_it != ri->GetAtomList().end(); ++at_it) { LOG_DEBUG(" adding atom " << (*at_it)->GetName()); ++atom_count; if (flags & QueryFlag::NO_BONDS) { res.AddAtom(*at_it); } else { update_bond_table(bond_table, *at_it, res.AddAtom(*at_it)); } } } else if (indeterminate(r)) { // Test atoms r_added = false; ResidueView res; AtomImplList::const_iterator at_it = ri->GetAtomList().begin(); for (;at_it != ri->GetAtomList().end(); ++at_it) { tribool a = qs.EvalAtom(*at_it); if (indeterminate(a) || a==true) { if (!c_added) { c_added = true; ++chain_count; chain = view.AddChain(ci); } if (!r_added) { r_added = true; ++residue_count; res = chain.AddResidue(ri); } if (flags & QueryFlag::MATCH_RESIDUES) { AtomImplList::const_iterator at_it2 = ri->GetAtomList().begin(); for (;at_it2 != ri->GetAtomList().end(); ++at_it2) { LOG_DEBUG(" adding atom " << (*at_it2)->GetQualifiedName()); ++atom_count; if (flags & QueryFlag::NO_BONDS) { res.AddAtom(*at_it2); } else { update_bond_table(bond_table, *at_it2, res.AddAtom(*at_it2)); } } break; } else { LOG_DEBUG(" adding atom " << (*at_it)->GetName()); ++atom_count; if (flags & QueryFlag::NO_BONDS) { res.AddAtom(*at_it); } else { update_bond_table(bond_table, *at_it, res.AddAtom(*at_it)); } } } qs.Reset(Prop::ATOM); } } qs.Reset(Prop::RESIDUE); } } if(!always_true) { qs.Reset(Prop::CHAIN); } } int bond_count=0; if (!(flags & QueryFlag::NO_BONDS)) { LOG_DEBUG("adding bonds"); typedef BondTableType::MapType::const_iterator ConstIt; for (ConstIt it=bond_table.bonds.begin();it!=bond_table.bonds.end(); ++it) { if ((QueryFlag::EXCLUSIVE_BONDS & flags) || it->second.IsComplete()) { BondHandle bh=it->second.bond; view.AddBond(bh); bond_count++; } } } LOG_VERBOSE("selected " << chain_count << " chain(s), " << residue_count << " residue(s), " << atom_count << " atom(s)" << " " << bond_count << " bond(s)"); return view; } EntityView EntityImpl::Select(const EntityHandle& h, const Query& q) const { return do_selection<false>(h, q, default_query_flags_); } EntityView EntityImpl::Select(const EntityHandle& h, const Query& q, QueryFlags flags) const { return do_selection<false>(h, q, flags); } EntityView EntityImpl::CreateFullView(const EntityHandle& h) const { return do_selection<true>(h, Query(), 0); } ChainImplPtr EntityImpl::FindChain(const String& name) const { ChainImplList::const_iterator i; for(i=chain_list_.begin(); i!=chain_list_.end();++i) { if ((*i)->GetName()==name) return *i; } return ChainImplPtr(); } ResidueImplPtr EntityImpl::FindResidue(const String& chain_name, const ResNum& residue) const { ChainImplPtr chain=this->FindChain(chain_name); if (!chain) return ResidueImplPtr(); return chain->FindResidue(residue); } AtomImplPtr EntityImpl::FindAtom(const String& chain_name, const ResNum& residue, const String& atom_name) const { ResidueImplPtr res=this->FindResidue(chain_name, residue); if (!res) return AtomImplPtr(); return res->FindAtom(atom_name); } void EntityImpl::DeleteFromConnMap(const ConnectorImplP& conn) { connector_map_.erase(conn.get()); } TorsionImplMap& EntityImpl::GetTorsionMap() { return torsion_map_; } void EntityImpl::DeleteChain(const ChainImplPtr& chain) { if (chain && chain->GetEntity().get()==this) { chain->DeleteAllResidues(); chain_list_.erase(this->GetChain(chain->GetName())); } } void EntityImpl::UpdateOrganizer() { atom_organizer_.Clear(); for (AtomImplMap::const_iterator i=atom_map_.begin(), e=atom_map_.end(); i!=e; ++i) { atom_organizer_.Add(i->second, i->second->TransformedPos()); } dirty_flags_&=~DirtyOrganizer; } void EntityImpl::IncXCSEditorCount() { xcs_editor_count_++; } void EntityImpl::UpdateXCSIfNeeded() { if (dirty_flags_ & DirtyTrace) { this->TraceDirectionality(); dirty_flags_&=~DirtyTrace; } if (dirty_flags_ & DirtyXCS) { LOG_VERBOSE("XCS marked dirty. Updating"); this->UpdateFromICS(); dirty_flags_&=~DirtyXCS; } } void EntityImpl::MarkXCSDirty() { dirty_flags_|=DirtyXCS|DirtyOrganizer; } void EntityImpl::MarkICSDirty() { dirty_flags_|=DirtyICS; } void EntityImpl::MarkTraceDirty() { dirty_flags_|=DirtyTrace; } void EntityImpl::MarkOrganizerDirty() { dirty_flags_|=DirtyOrganizer; } void EntityImpl::EnableICS() { dirty_flags_=(dirty_flags_ & ~DisableICS) | DirtyTrace; this->UpdateICSIfNeeded(); } void EntityImpl::UpdateICSIfNeeded() { if (dirty_flags_ & DisableICS) { return; } if (dirty_flags_ & DirtyTrace) { LOG_VERBOSE("rerunning directionality trace"); this->TraceDirectionality(); this->UpdateFromXCS(); dirty_flags_&=~DirtyTrace; } else if (dirty_flags_ & DirtyICS) { LOG_VERBOSE("updating internal from external coordinates"); this->UpdateFromXCS(); dirty_flags_&=~DirtyICS; } } void EntityImpl::DecXCSEditorCount() { assert(xcs_editor_count_>0); xcs_editor_count_--; if (xcs_editor_count_==0) { this->UpdateICSIfNeeded(); this->UpdateOrganizerIfNeeded(); } } void EntityImpl::UpdateOrganizerIfNeeded() { if (dirty_flags_ & DirtyOrganizer) { LOG_DEBUG("atom organizer marked as dirty. updating"); this->UpdateOrganizer(); dirty_flags_&=~DirtyOrganizer; } } void EntityImpl::IncICSEditorCount() { ics_editor_count_++; } bool EntityImpl::HasICS() const { return !(dirty_flags_ & DisableICS); } void EntityImpl::DecICSEditorCount() { assert(ics_editor_count_>0); ics_editor_count_--; if (ics_editor_count_==0) { this->UpdateXCSIfNeeded(); this->UpdateOrganizerIfNeeded(); } } impl::ChainImplList::iterator EntityImpl::GetChain(const String& name) { impl::ChainImplList& cc=this->GetChainList(); for (impl::ChainImplList::iterator i=cc.begin(), e=cc.end(); i!=e; ++i) { if ((*i)->GetName()==name) { return i; } } return cc.end(); } pointer_it<ChainImplPtr> EntityImpl::GetChainIter(const String& name) { impl::ChainImplList& cc=this->GetChainList(); for (size_t i=0; i<cc.size(); ++i) { if (cc[i]->GetName()==name) { return &cc.front()+i; } } return pointer_it<ChainImplPtr>(NULL); } void EntityImpl::RenameChain(ChainImplPtr chain, const String& new_name) { //ChainImplList::iterator i; ChainImplPtr ch=this->FindChain(new_name); if ((ch) && (ch != chain)) { throw IntegrityError("unable to rename chain '"+chain->GetName()+ "' to '"+new_name+"', since there is already a chain " "with that name"); } chain->SetName(new_name); } void EntityImpl::UpdateTransformedPos(){ for(AtomImplMap::iterator it = atom_map_.begin();it!=atom_map_.end();++it) { it->second->TransformedPos()=geom::Vec3(transformation_matrix_*geom::Vec4(it->second->OriginalPos())); } } const String& EntityImpl::GetName() const { return name_; } void EntityImpl::SetName(const String& ent_name){ name_=ent_name; } void EntityImpl::ReorderAllResidues() { for(ChainImplList::iterator cit=chain_list_.begin();cit!=chain_list_.end();++cit) { (*cit)->ReorderResidues(); } } void EntityImpl::RenumberAllResidues(int start, bool keep_spacing) { for(ChainImplList::iterator cit=chain_list_.begin();cit!=chain_list_.end();++cit) { (*cit)->RenumberAllResidues(start, keep_spacing); } } }}} // ns