backup platonyzer

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

include/cif++/DistanceMap.h

@@ -4,8 +4,6 @@
 
 
 #include <unordered_map>
 #include <unordered_map>
 
 
-//#include <boost/hash/hash.hpp>
-
 #include <clipper/clipper.h>
 #include <clipper/clipper.h>
 
 
 #include "cif++/Structure.h"
 #include "cif++/Structure.h"
@@ -19,7 +17,7 @@ class DistanceMap
 	DistanceMap(const Structure& p, const clipper::Spacegroup& spacegroup, const clipper::Cell& cell,
 	DistanceMap(const Structure& p, const clipper::Spacegroup& spacegroup, const clipper::Cell& cell,
 		float maxDistance);
 		float maxDistance);
 
 
-	// simplified version for subsets of atoms (used in refining e.g.)	
+	// simplified version for subsets of atoms (used in refining e.g.)
 //	DistanceMap(const Structure& p, const std::vector<Atom>& atoms);
 //	DistanceMap(const Structure& p, const std::vector<Atom>& atoms);
 	
 	
 	DistanceMap(const DistanceMap&) = delete;
 	DistanceMap(const DistanceMap&) = delete;
@@ -30,6 +28,12 @@ class DistanceMap
 	std::vector<Atom> near(const Atom& a, float maxDistance = 3.5f) const;
 	std::vector<Atom> near(const Atom& a, float maxDistance = 3.5f) const;
 	std::vector<Atom> near(const Point& p, float maxDistance = 3.5f) const;
 	std::vector<Atom> near(const Point& p, float maxDistance = 3.5f) const;
 
 
+	static clipper::Coord_orth
+		CalculateOffsetForCell(const Structure& p, const clipper::Spacegroup& spacegroup, const clipper::Cell& cell);
+
+	static std::vector<clipper::RTop_orth>
+		AlternativeSites(const clipper::Spacegroup& spacegroup, const clipper::Cell& cell);
+
   private:
   private:
 
 
 	typedef std::map<std::tuple<size_t,size_t>,std::tuple<float,int32_t>> DistMap;
 	typedef std::map<std::tuple<size_t,size_t>,std::tuple<float,int32_t>> DistMap;
@@ -49,4 +53,112 @@ class DistanceMap
 	std::vector<clipper::RTop_orth>			mRtOrth;
 	std::vector<clipper::RTop_orth>			mRtOrth;
 };
 };
 
 
+// --------------------------------------------------------------------
+
+class SymmetryAtomIteratorFactory
+{
+  public:
+	// SymmetryAtomIteratorFactory(const Structure& p);
+	SymmetryAtomIteratorFactory(const Structure& p, const clipper::Spacegroup& spacegroup, const clipper::Cell& cell)
+		: mD(DistanceMap::CalculateOffsetForCell(p, spacegroup, cell))
+		, mRtOrth(DistanceMap::AlternativeSites(spacegroup, cell)) {}
+
+	SymmetryAtomIteratorFactory(const SymmetryAtomIteratorFactory&) = delete;
+	SymmetryAtomIteratorFactory& operator=(const SymmetryAtomIteratorFactory&) = delete;
+
+	class SymmetryAtomIterator : public std::iterator<std::forward_iterator_tag, const Atom>
+	{
+	  public:
+		typedef std::iterator<std::forward_iterator_tag, const Atom>	baseType;
+		typedef typename baseType::pointer								pointer;
+		typedef typename baseType::reference							reference;
+
+		SymmetryAtomIterator(const SymmetryAtomIteratorFactory& factory, const Atom& atom)
+			: m_f(&factory), m_i(0), m_a(atom), m_c(atom) {}
+
+		SymmetryAtomIterator(const SymmetryAtomIteratorFactory& factory, const Atom& atom, int)
+			: SymmetryAtomIterator(factory, atom)
+		{
+			m_i = m_f->mRtOrth.size();
+		}
+
+		SymmetryAtomIterator(const SymmetryAtomIterator& iter)
+			: m_f(iter.m_f), m_i(iter.m_i), m_a(iter.m_a), m_c(iter.m_c) {}
+
+		SymmetryAtomIterator& operator=(const SymmetryAtomIterator& iter)
+		{
+			if (this != &iter)
+			{
+				m_f = iter.m_f;
+				m_i = iter.m_i;
+				m_a = iter.m_a;
+				m_c = iter.m_c;
+			}
+			return *this;
+		}
+
+		reference operator*()		{ return m_c; }
+		pointer operator->()		{ return &m_c; }
+
+		SymmetryAtomIterator operator++()
+		{
+			if (++m_i < m_f->mRtOrth.size())
+				m_c = m_a.symmetryCopy(m_f->mD, m_f->mRtOrth[m_i]);
+			return *this;
+		}
+
+		SymmetryAtomIterator operator++(int)
+		{
+			SymmetryAtomIterator result(*this);
+			this->operator++();
+			return result;
+		}
+
+		bool operator==(const SymmetryAtomIterator& iter) const
+		{
+			return m_f == iter.m_f and m_i == iter.m_i;
+		}
+
+		bool operator!=(const SymmetryAtomIterator& iter) const
+		{
+			return m_f != iter.m_f or m_i != iter.m_i;
+		}
+
+	  private:
+		const SymmetryAtomIteratorFactory* m_f;
+		size_t m_i;
+		Atom m_a, m_c;
+	};
+
+	class SymmetryAtomIteratorRange
+	{
+	  public:
+		SymmetryAtomIteratorRange(const SymmetryAtomIteratorFactory& f, const Atom& a)
+			: m_f(f), m_a(a) {}
+
+		SymmetryAtomIterator begin()
+		{
+			return SymmetryAtomIterator(m_f, m_a);
+		}
+
+		SymmetryAtomIterator end()
+		{
+			return SymmetryAtomIterator(m_f, m_a, 1);
+		}
+
+	  private:
+		const SymmetryAtomIteratorFactory& m_f;
+		Atom m_a;
+	};
+
+	SymmetryAtomIteratorRange operator()(const Atom& a) const
+	{
+		return SymmetryAtomIteratorRange(*this, a);
+	}
+
+  private:
+	Point									mD;			// needed to move atoms to center
+	std::vector<clipper::RTop_orth>			mRtOrth;
+};
+
 }
 }

include/cif++/Structure.h

@@ -83,6 +83,8 @@ class Atom
 	void location(Point p);
 	void location(Point p);
 	
 	
 	Atom symmetryCopy(const Point& d, const clipper::RTop_orth& rt);
 	Atom symmetryCopy(const Point& d, const clipper::RTop_orth& rt);
+	bool isSymmetryCopy() const;
+	std::string symmetry() const;
 	
 	
 	const Compound& comp() const;
 	const Compound& comp() const;
 	bool isWater() const;
 	bool isWater() const;

src/DistanceMap.cpp

@@ -28,14 +28,97 @@ inline ostream& operator<<(ostream& os, const Atom& a)
 
 
 // --------------------------------------------------------------------
 // --------------------------------------------------------------------
 
 
-vector<clipper::RTop_orth> AlternativeSites(const clipper::Spacegroup& spacegroup,
+clipper::Coord_orth DistanceMap::CalculateOffsetForCell(const Structure& p, const clipper::Spacegroup& spacegroup, const clipper::Cell& cell)
+{
+	auto& atoms = p.atoms();
+	size_t dim = atoms.size();
+	
+	vector<clipper::Coord_orth> locations;
+	locations.reserve(dim);
+	
+	// bounding box
+	Point pMin(numeric_limits<float>::max(), numeric_limits<float>::max(), numeric_limits<float>::max()),
+		  pMax(numeric_limits<float>::min(), numeric_limits<float>::min(), numeric_limits<float>::min());
+	
+	for (auto& atom: atoms)
+	{
+		auto p = atom.location();
+		locations.push_back(p);
+
+		if (pMin.mX > p.mX)
+			pMin.mX = p.mX;
+		if (pMin.mY > p.mY)
+			pMin.mY = p.mY;
+		if (pMin.mZ > p.mZ)
+			pMin.mZ = p.mZ;
+
+		if (pMax.mX < p.mX)
+			pMax.mX = p.mX;
+		if (pMax.mY < p.mY)
+			pMax.mY = p.mY;
+		if (pMax.mZ < p.mZ)
+			pMax.mZ = p.mZ;
+	};
+	
+	// correct locations so that the median of x, y and z are inside the cell
+	vector<float> c(dim);
+	auto median = [&]()
+	{
+		return dim % 1 == 0
+			? c[dim / 2]
+			: (c[dim / 2 - 1] + c[dim / 2]) / 2;
+	};
+	
+	transform(locations.begin(), locations.end(), c.begin(), [](auto& l) { return l[0]; });
+	sort(c.begin(), c.end());
+	float mx = median();
+	
+	transform(locations.begin(), locations.end(), c.begin(), [](auto& l) { return l[1]; });
+	sort(c.begin(), c.end());
+	float my = median();
+	
+	transform(locations.begin(), locations.end(), c.begin(), [](auto& l) { return l[2]; });
+	sort(c.begin(), c.end());
+	float mz = median();
+
+	if (cif::VERBOSE > 1)
+		cerr << "median position of atoms: " << Point(mx, my, mz) << endl;
+	
+	auto calculateD = [&](float m, float c)
+	{
+		float d = 0;
+		while (m + d < -(c / 2))
+			d += c;
+		while (m + d > (c / 2))
+			d -= c;
+		return d;
+	};
+
+	Point D;
+
+	D.mX = calculateD(mx, cell.a());
+	D.mY = calculateD(my, cell.b());
+	D.mZ = calculateD(mz, cell.c());
+	
+	if (D.mX != 0 or D.mY != 0 or D.mZ != 0)
+	{
+		if (cif::VERBOSE)
+			cerr << "moving coorinates by " << D.mX << ", " << D.mY << " and " << D.mZ << endl;
+	}
+
+	return D;	
+}
+
+// --------------------------------------------------------------------
+
+vector<clipper::RTop_orth> DistanceMap::AlternativeSites(const clipper::Spacegroup& spacegroup,
 	const clipper::Cell& cell)
 	const clipper::Cell& cell)
 {
 {
 	vector<clipper::RTop_orth> result;
 	vector<clipper::RTop_orth> result;
 	
 	
-	// to make index 0 equal to identity, no 
+	// to make the operation at index 0 equal to identity
 	result.push_back(clipper::RTop_orth::identity());
 	result.push_back(clipper::RTop_orth::identity());
-	
+
 	for (int i = 0; i < spacegroup.num_symops(); ++i)
 	for (int i = 0; i < spacegroup.num_symops(); ++i)
 	{
 	{
 		const auto& symop = spacegroup.symop(i);
 		const auto& symop = spacegroup.symop(i);
@@ -44,10 +127,14 @@ vector<clipper::RTop_orth> AlternativeSites(const clipper::Spacegroup& spacegrou
 			for (int v: { -1, 0, 1})
 			for (int v: { -1, 0, 1})
 				for (int w: { -1, 0, 1})
 				for (int w: { -1, 0, 1})
 				{
 				{
-					result.push_back(
-						clipper::RTop_frac(
+					if (i == 0 and u == 0 and v == 0 and w == 0)
+						continue;
+
+					auto rtop = clipper::RTop_frac(
 							symop.rot(), symop.trn() + clipper::Vec3<>(u, v, w)
 							symop.rot(), symop.trn() + clipper::Vec3<>(u, v, w)
-						).rtop_orth(cell));
+						).rtop_orth(cell);
+
+					result.push_back(move(rtop));
 				}
 				}
 	}
 	}
 	
 	

src/Structure.cpp

@@ -154,6 +154,7 @@ struct AtomImpl
 		, mRefcount(1), mRow(i.mRow), mCompound(i.mCompound)
 		, mRefcount(1), mRow(i.mRow), mCompound(i.mCompound)
 		, mRadius(i.mRadius), mCachedProperties(i.mCachedProperties)
 		, mRadius(i.mRadius), mCachedProperties(i.mCachedProperties)
 		, mSymmetryCopy(i.mSymmetryCopy), mClone(true)
 		, mSymmetryCopy(i.mSymmetryCopy), mClone(true)
+		, mSymmetry(i.mSymmetry)
 	{
 	{
 	}
 	}
 	
 	
@@ -185,6 +186,8 @@ struct AtomImpl
 		mLocation += d;
 		mLocation += d;
 		mLocation = ((clipper::Coord_orth)mLocation).transform(rt);
 		mLocation = ((clipper::Coord_orth)mLocation).transform(rt);
 		mLocation -= d;
 		mLocation -= d;
+
+		mSymmetry = clipper::Symop(rt).format();
 	}
 	}
 
 
 	void prefetch()
 	void prefetch()
@@ -379,6 +382,7 @@ struct AtomImpl
 	
 	
 	bool				mSymmetryCopy = false;
 	bool				mSymmetryCopy = false;
 	bool				mClone = false;
 	bool				mClone = false;
+	string				mSymmetry;
 };
 };
 
 
 //Atom::Atom(const File& f, const string& id)
 //Atom::Atom(const File& f, const string& id)
@@ -506,6 +510,11 @@ string Atom::labelAsymId() const
 	return mImpl->mAsymID;
 	return mImpl->mAsymID;
 }
 }
 
 
+string Atom::labelAltId() const
+{
+	return mImpl->mAltID;
+}
+
 int Atom::labelSeqId() const
 int Atom::labelSeqId() const
 {
 {
 	return mImpl->mSeqID;
 	return mImpl->mSeqID;
@@ -521,6 +530,16 @@ string Atom::authAtomId() const
 	return property<string>("auth_atom_id");
 	return property<string>("auth_atom_id");
 }
 }
 
 
+string Atom::authAltId() const
+{
+	return property<string>("auth_alt_id");
+}
+
+string Atom::pdbxAuthInsCode() const
+{
+	return property<string>("pdbx_PDB_ins_code");
+}
+
 string Atom::authCompId() const
 string Atom::authCompId() const
 {
 {
 	return property<string>("auth_comp_id");
 	return property<string>("auth_comp_id");
@@ -560,6 +579,16 @@ Atom Atom::symmetryCopy(const Point& d, const clipper::RTop_orth& rt)
 	return Atom(new AtomImpl(*mImpl, d, rt));
 	return Atom(new AtomImpl(*mImpl, d, rt));
 }
 }
 
 
+bool Atom::isSymmetryCopy() const
+{
+	return mImpl->mSymmetryCopy;
+}
+
+string Atom::symmetry() const
+{
+	return mImpl->mSymmetry;
+}
+
 const Compound& Atom::comp() const
 const Compound& Atom::comp() const
 {
 {
 	return mImpl->comp();
 	return mImpl->comp();