diff --git a/modules/io/src/mol/sdf_writer.cc b/modules/io/src/mol/sdf_writer.cc index 9f3154e9b99f91a47354a7dcb2f9ccf31cd1b758..40c0660153853bc5edbfeb8e94384194a8daf63b 100644 --- a/modules/io/src/mol/sdf_writer.cc +++ b/modules/io/src/mol/sdf_writer.cc @@ -27,6 +27,7 @@ #include <ost/mol/chain_view.hh> #include <ost/mol/bond_handle.hh> #include <boost/regex.hpp> +#include <boost/bind.hpp> namespace ost { namespace io { @@ -38,7 +39,7 @@ namespace { public: SDFAtomWriter(std::ostream& ostream, std::map<long, int>& atom_indices) : ostr_(ostream), atom_indices_(atom_indices), counter_(0) { - atom_indices_.clear(); + atom_indices_.clear(); } private: public: @@ -60,23 +61,61 @@ namespace { class SDFBondWriter : public mol::EntityViewVisitor { public: - SDFBondWriter(std::ostream& ostream, std::map<long, int>& atom_indices) + SDFBondWriter(std::ostream& ostream, + const std::map<long, int>& atom_indices) : ostr_(ostream), atom_indices_(atom_indices), counter_(0) { } private: + // compare two atoms according to their indices (used for sorting) + bool CompareAtomIdx(const mol::AtomView& first, + const mol::AtomView& second) { + std::map<long, int>::const_iterator aidx_first( + atom_indices_.find(first.GetHashCode())); + std::map<long, int>::const_iterator aidx_second( + atom_indices_.find(second.GetHashCode())); + + if(aidx_first==atom_indices_.end() || aidx_second==atom_indices_.end()) { + throw IOException("Cannot write bond: atom idx not found for sorting"); + } + return (aidx_first->second < aidx_second->second); + } + public: virtual bool VisitAtom(const mol::AtomView& atom) { - counter_++; + ++counter_; // current atom index + + // get all neighboring atoms and sort them according to their atom index mol::AtomViewList atoms = atom.GetBondPartners(); - mol::AtomViewList::iterator atom_iter = atoms.begin(); - for(; atom_iter != atoms.end(); ++atom_iter) { - int atom_index = atom_indices_.find((*atom_iter).GetHashCode())->second; - if(atom_index > counter_) { - int type = 1; - mol::BondHandle bond = atom.GetHandle().FindBondToAtom(atom_iter->GetHandle()); - if(bond.IsValid()) type = bond.GetBondOrder(); + std::sort(atoms.begin(), atoms.end(), bind(&SDFBondWriter::CompareAtomIdx, + this, _1, _2)); + + // iterate all neighboring atoms and print bonds to all atoms with index + // larger than current atom index + for(mol::AtomViewList::iterator atom_iter = atoms.begin(); + atom_iter != atoms.end(); ++atom_iter) { + std::map<long, int>::const_iterator aidx( + atom_indices_.find((*atom_iter).GetHashCode())); + + // check if index was found + if(aidx==atom_indices_.end()) { + throw IOException("Cannot write bond between " + + atom.GetQualifiedName() + " and " + + atom_iter->GetQualifiedName() + + ": atom index not found"); + } + + // only print bonds to atoms with larger index than current index + if(aidx->second > counter_) { + mol::BondHandle bond(atom.GetHandle().FindBondToAtom( + atom_iter->GetHandle())); + if(!bond.IsValid()) { + throw IOException("Bond is invalid between " + + atom.GetQualifiedName() + " and " + + atom_iter->GetQualifiedName()); + } + int type = bond.GetBondOrder(); ostr_ << format("%3i") % counter_ - << format("%3i") % atom_index + << format("%3i") % aidx->second << format("%3i") % type << " 0 0 0" << std::endl; @@ -87,17 +126,17 @@ namespace { private: std::ostream& ostr_; - std::map<long, int>& atom_indices_; + const std::map<long, int>& atom_indices_; int counter_; }; } SDFWriter::SDFWriter(std::ostream& ostream) - : outfile_(), ostr_(ostream), counter_(0) { + : outfile_(), ostr_(ostream), counter_(0), atom_indices_() { } SDFWriter::SDFWriter(const String& filename) - : outfile_(filename.c_str()), ostr_(outfile_), counter_(0) { + : outfile_(filename.c_str()), ostr_(outfile_), counter_(0), atom_indices_() { } SDFWriter::SDFWriter(const boost::filesystem::path& filename): @@ -106,7 +145,7 @@ SDFWriter::SDFWriter(const boost::filesystem::path& filename): #else outfile_(filename.file_string().c_str()), #endif - ostr_(outfile_), counter_(0) { + ostr_(outfile_), counter_(0), atom_indices_() { } void SDFWriter::Write(const mol::EntityView& ent) {