edit ignore

git-svn-id: svn+ssh://gitlab/srv/svn-repos/pdb-redo/trunk@220 a1961a4f-ab94-4bcc-80e8-33b5a54de466

include/cif++/Compound.h

@@ -119,6 +119,8 @@ class Compound
 	bool isWater() const;
 
 	std::vector<ChiralCentre> chiralCentres() const			{ return mChiralCentres; }
+	
+	bool isIsomerOf(const Compound& c) const;
 
   private:
 //	Entity&						mEntity;

include/cif++/Structure.h

@@ -92,6 +92,7 @@ class Atom
 	AtomType type() const;
 
 	Point location() const;
+	void location(Point p);
 
 	const Compound& comp() const;
 	const Entity& ent() const;
@@ -338,6 +339,7 @@ class Structure
 	// Actions
 	void removeAtom(Atom& a);
 	void swapAtoms(Atom& a1, Atom& a2);	// swap the labels for these atoms
+	void moveAtom(Atom& a, Point p);	// move atom to a new location
 	
   private:
 	friend Polymer;

src/Cif++.cpp

@@ -1883,9 +1883,9 @@ void Row::assign(const string& name, const string& value, bool emplacing)
 			if (childCat == nullptr)
 				continue;
 
-#if DEBUG
-cerr << "fixing linked item " << child->mCategory->mName << '.' << child->mTag << endl;
-#endif
+//#if DEBUG
+//cerr << "fixing linked item " << child->mCategory->mName << '.' << child->mTag << endl;
+//#endif
 				
 			auto rows = childCat->find(Key(child->mTag) == oldValue);
 			for (auto& cr: rows)
@@ -2008,9 +2008,9 @@ void Row::swap(const string& name, ItemRow* a, ItemRow* b)
 			if (childCat == nullptr)
 				continue;
 
-#if DEBUG
-cerr << "fixing linked item " << child->mCategory->mName << '.' << child->mTag << endl;
-#endif
+//#if DEBUG
+//cerr << "fixing linked item " << child->mCategory->mName << '.' << child->mTag << endl;
+//#endif
 			if (ai != nullptr)
 			{
 				auto rows = childCat->find(Key(child->mTag) == string(ai->mText));

src/Compound.cpp

@@ -7,6 +7,7 @@
 #include <boost/algorithm/string.hpp>
 #include <boost/filesystem/operations.hpp>
 #include <boost/filesystem/fstream.hpp>
+#include <boost/thread.hpp>
 
 #include "cif++/Cif++.h"
 #include "cif++/Compound.h"
@@ -23,6 +24,8 @@ class CompoundFactory
   public:
 	
 	static CompoundFactory& instance();
+
+	const Compound* get(string id);
 	const Compound* create(string id);
 
   private:
@@ -33,8 +36,165 @@ class CompoundFactory
 
 	fs::path mClibdMon;
 	vector<Compound*> mCompounds;
+	boost::shared_mutex mMutex;
 };
 
+// --------------------------------------------------------------------
+// Compound helper classes
+
+struct CompoundAtomLess
+{
+	bool operator()(const CompoundAtom& a, const CompoundAtom& b) const
+	{
+		int d = a.id.compare(b.id);
+		if (d == 0)
+			d = a.typeSymbol - b.typeSymbol;
+		return d < 0;
+	}
+};
+
+struct CompoundBondLess
+{
+	bool operator()(const CompoundBond& a, const CompoundBond& b) const
+	{
+		int d = a.atomID[0].compare(b.atomID[0]);
+		if (d == 0)
+			d = a.atomID[1].compare(b.atomID[1]);
+		if (d == 0)
+			d = a.type - b.type;
+		return d < 0;
+	}
+};
+
+// --------------------------------------------------------------------
+// Brute force comparison of two structures, when they are isomers the
+// mapping between the atoms of both is returned
+
+struct Node
+{
+	string									id;
+	AtomType								symbol;
+	vector<tuple<size_t,CompoundBondType>>	links;
+};
+
+// Check to see if the nodes a[iA] and b[iB] are the start of a similar sub structure
+bool SubStructuresAreIsomeric(
+	const vector<Node>& a, const vector<Node>& b, size_t iA, size_t iB,
+	vector<bool> visitedA, vector<bool> visitedB, vector<tuple<string,string>>& outMapping)
+{
+	bool result = false;
+
+	auto& na = a[iA];
+	auto& nb = b[iB];
+	size_t N = na.links.size();
+	
+	if (na.symbol == nb.symbol and nb.links.size() == N)
+	{
+		result = true;
+		
+		visitedA[iA] = true;
+		visitedB[iB] = true;
+		
+		// we now have two sets of links to compare. 
+		// To keep code clean, first create the list of possible permutations
+		
+		vector<size_t> ilb(N);
+		iota(ilb.begin(), ilb.end(), 0);
+		
+		for (;;)
+		{
+			result = true;
+			vector<tuple<string,string>> m = outMapping;
+			
+			for (size_t i = 0; result and i < N; ++i)
+			{
+				size_t lA, lB;
+				CompoundBondType typeA, typeB;
+
+				tie(lA, typeA) = na.links[i];
+				tie(lB, typeB) = nb.links[ilb[i]];
+				
+				if (typeA != typeB or visitedA[lA] != visitedB[lB])
+					result = false;
+				else if (not visitedA[lA])
+					result = SubStructuresAreIsomeric(a, b, lA, lB, visitedA, visitedB, m);
+			}
+			
+			if (result)
+			{
+				swap(m, outMapping);
+				break;
+			}
+
+			if (not next_permutation(ilb.begin(), ilb.end()))
+				break;
+		}
+		
+		if (result)
+			outMapping.emplace_back(na.id, nb.id);
+	}
+	
+	return result;
+}
+
+bool StructuresAreIsomeric(vector<CompoundAtom> atomsA, const vector<CompoundBond>& bondsA,
+	vector<CompoundAtom> atomsB, const vector<CompoundBond>& bondsB,
+	vector<tuple<string,string>>& outMapping)
+{
+	assert(atomsA.size() == atomsB.size());
+	assert(bondsA.size() == bondsB.size());
+
+	vector<Node> a, b;
+	map<string,size_t> ma, mb;
+	
+	for (auto& atomA: atomsA)
+	{
+		ma[atomA.id] = a.size();
+		a.push_back({atomA.id, atomA.typeSymbol});
+	}
+	
+	for (auto& bondA: bondsA)
+	{
+		size_t atom1 = ma.at(bondA.atomID[0]);
+		size_t atom2 = ma.at(bondA.atomID[1]);
+		
+		a[atom1].links.emplace_back(atom2, bondA.type);
+		a[atom2].links.emplace_back(atom1, bondA.type);
+	}
+
+	for (auto& atomB: atomsB)
+	{
+		mb[atomB.id] = b.size();
+		b.push_back({atomB.id, atomB.typeSymbol});
+	}
+
+	for (auto& bondB: bondsB)
+	{
+		size_t atom1 = mb.at(bondB.atomID[0]);
+		size_t atom2 = mb.at(bondB.atomID[1]);
+		
+		b[atom1].links.emplace_back(atom2, bondB.type);
+		b[atom2].links.emplace_back(atom1, bondB.type);
+	}
+
+	size_t N = atomsA.size();
+	
+	bool result = false;
+
+	// try each atom in B to see if it can be traced to be similar to A starting at zero
+	for (size_t ib = 0; ib < N; ++ib)
+	{
+		vector<bool> va(N, false), vb(N, false);
+
+		if (SubStructuresAreIsomeric(a, b, 0, ib, va, vb, outMapping))
+		{
+			result = true;
+			break;
+		}
+	}
+	
+	return result;
+}
 
 // --------------------------------------------------------------------
 // Compound
@@ -160,7 +320,10 @@ CompoundAtom Compound::getAtomById(const string& atomId) const
 
 const Compound* Compound::create(const string& id)
 {
-	return CompoundFactory::instance().create(id);
+	auto result = CompoundFactory::instance().get(id);
+	if (result == nullptr)
+		result = CompoundFactory::instance().create(id);
+	return result;
 }
 
 // --------------------------------------------------------------------
@@ -188,8 +351,31 @@ CompoundFactory& CompoundFactory::instance()
 }
 
 // id is the three letter code
+
+const Compound* CompoundFactory::get(std::string id)
+{
+	boost::shared_lock<boost::shared_mutex> lock(mMutex);
+
+	ba::to_upper(id);
+
+	Compound* result = nullptr;
+	
+	for (auto cmp: mCompounds)
+	{
+		if (cmp->id() == id)
+		{
+			result = cmp;
+			break;
+		}
+	}
+	
+	return result;
+}
+
 const Compound* CompoundFactory::create(std::string id)
 {
+	boost::upgrade_lock<boost::shared_mutex> lock(mMutex);
+
 	ba::to_upper(id);
 
 	Compound* result = nullptr;
@@ -206,6 +392,10 @@ const Compound* CompoundFactory::create(std::string id)
 	if (result == nullptr)
 	{
 		fs::path resFile = mClibdMon / ba::to_lower_copy(id.substr(0, 1)) / (id + ".cif");
+
+		if (not fs::exists(resFile) and	(id == "COM" or id == "CON" or "PRN")) 		// seriously...
+			mClibdMon / ba::to_lower_copy(id.substr(0, 1)) / (id + '_' + id + ".cif");
+
 		fs::ifstream file(resFile);
 		if (file.is_open())
 		{
@@ -246,6 +436,7 @@ const Compound* CompoundFactory::create(std::string id)
 					id, AtomTypeTraits(symbol).type(), energy, charge
 				});
 			}
+			sort(atoms.begin(), atoms.end(), CompoundAtomLess());
 
 			auto& compBonds = cf["comp_" + id]["chem_comp_bond"];
 			
@@ -272,8 +463,12 @@ const Compound* CompoundFactory::create(std::string id)
 					b.type = singleBond;
 				}
 				
+				if (b.atomID[0] > b.atomID[1])
+					swap(b.atomID[0], b.atomID[1]);
+				
 				bonds.push_back(b);
 			}
+			sort(bonds.begin(), bonds.end(), CompoundBondLess());
 
 			auto& compChir = cf["comp_" + id]["chem_comp_chir"];
 			
@@ -303,9 +498,11 @@ const Compound* CompoundFactory::create(std::string id)
 				
 				chiralCentres.push_back(cc);
 			}
-			
+
 			result = new Compound(id, name, group, move(atoms), move(bonds),
 				move(chiralCentres));
+
+			boost::upgrade_to_unique_lock<boost::shared_mutex> uniqueLock(lock);
 			mCompounds.push_back(result);
 		}
 	}
@@ -337,4 +534,77 @@ float Compound::atomBondValue(const string& atomId_1, const string& atomId_2) co
 	return i != mBonds.end() ? i->distance : 0;
 }
 
+bool Compound::isIsomerOf(const Compound& c) const
+{
+	bool result = false;
+	
+	for (;;)
+	{
+		// easy tests first
+		if (mId == c.mId)
+		{
+			result = true;
+			break;
+		}
+		
+		if (mAtoms.size() != c.mAtoms.size())
+			break;
+		
+		if (mBonds.size() != c.mBonds.size())
+			break;
+		
+		if (mChiralCentres.size() != c.mChiralCentres.size())
+			break;
+
+		// same number of atoms of each type?
+		map<AtomType,int> aTypeCount, bTypeCount;
+		
+		bool sameAtomNames = true;
+		for (size_t i = 0; i < mAtoms.size(); ++i)
+		{
+			auto& a = mAtoms[i];
+			auto& b = c.mAtoms[i];
+			
+			aTypeCount[a.typeSymbol] += 1;
+			bTypeCount[b.typeSymbol] += 1;
+			
+			if (a.id != b.id or a.typeSymbol != b.typeSymbol)
+				sameAtomNames = false;
+		}
+		
+		if (not sameAtomNames and aTypeCount != bTypeCount)
+			break;
+		
+		bool sameBonds = sameAtomNames;
+		for (size_t i = 0; sameBonds and i < mBonds.size(); ++i)
+		{
+			sameBonds =
+				mBonds[i].atomID[0] == c.mBonds[i].atomID[0] and
+				mBonds[i].atomID[1] == c.mBonds[i].atomID[1] and
+				mBonds[i].type == c.mBonds[i].type;
+		}
+		
+		if (sameBonds)
+		{
+			result = true;
+			break;
+		}
+		
+		// implement rest of tests
+		
+		vector<tuple<string,string>> mapping;
+		result = StructuresAreIsomeric(mAtoms, mBonds, c.mAtoms, c.mBonds, mapping);
+		
+		if (VERBOSE and result)
+		{
+			for (auto& m: mapping)
+				cerr << "  " << get<0>(m) << " => " << get<1>(m) << endl;
+		}
+		
+		break;	
+	}
+	
+	return result;
+}
+
 }

src/Structure.cpp

@@ -239,6 +239,15 @@ struct AtomImpl
 			delete this;
 	}
 	
+	void moveTo(const Point& p)
+	{
+		mRow["Cartn_x"] = p.getX();
+		mRow["Cartn_y"] = p.getY();
+		mRow["Cartn_z"] = p.getZ();
+
+		mLocation = p;
+	}
+	
 	const Compound& comp()
 	{
 		if (mCompound == nullptr)
@@ -426,6 +435,11 @@ Point Atom::location() const
 	return mImpl->mLocation;
 }
 
+void Atom::location(Point p)
+{
+	mImpl->moveTo(p);
+}
+
 const Compound& Atom::comp() const
 {
 	return mImpl->comp();
@@ -812,6 +826,7 @@ struct StructureImpl
 	
 	void removeAtom(Atom& a);
 	void swapAtoms(Atom& a1, Atom& a2);
+	void moveAtom(Atom& a, Point p);
 
 	File*			mFile;
 	uint32			mModelNr;
@@ -862,6 +877,11 @@ void StructureImpl::swapAtoms(Atom& a1, Atom& a2)
 	swap(l3, l4);
 }
 
+void StructureImpl::moveAtom(Atom& a, Point p)
+{
+	a.location(p);
+}
+
 Structure::Structure(File& f, uint32 modelNr)
 	: mImpl(new StructureImpl(*this, f, modelNr))
 {
@@ -1065,4 +1085,9 @@ void Structure::swapAtoms(Atom& a1, Atom& a2)
 	mImpl->swapAtoms(a1, a2);
 }
 
+void Structure::moveAtom(Atom& a, Point p)
+{
+	mImpl->moveAtom(a, p);
+}
+
 }