Merge branch 'develop' of github.com:PDB-REDO/libcifpp into develop

CMakeLists.txt

@@ -168,6 +168,8 @@ if(MSVC)
 endif()
 find_package(ZLIB REQUIRED)
 
+find_package(Eigen3 )
+
 include(FindFilesystem)
 list(APPEND CIFPP_REQUIRED_LIBRARIES ${STDCPPFS_LIBRARY})
 
@@ -277,7 +279,7 @@ target_include_directories(cifpp
 	"$<INSTALL_INTERFACE:${CMAKE_INSTALL_INCLUDEDIR}>"
 )
 
-target_link_libraries(cifpp PUBLIC Threads::Threads ZLIB::ZLIB ${CIFPP_REQUIRED_LIBRARIES})
+target_link_libraries(cifpp PUBLIC Threads::Threads ZLIB::ZLIB ${CIFPP_REQUIRED_LIBRARIES} Eigen3::Eigen)
 
 if(CMAKE_CXX_COMPILER_ID STREQUAL "AppleClang")
 	target_link_options(cifpp PRIVATE -undefined dynamic_lookup)
@@ -425,7 +427,7 @@ if(ENABLE_TESTING)
 
 		add_executable(${CIFPP_TEST} ${CIFPP_TEST_SOURCE})
 
-		target_link_libraries(${CIFPP_TEST} PRIVATE Threads::Threads cifpp::cifpp Boost::boost)
+		target_link_libraries(${CIFPP_TEST} PRIVATE Threads::Threads cifpp::cifpp Boost::boost Eigen3::Eigen)
 
 		if(MSVC)
 			# Specify unwind semantics so that MSVC knowns how to handle exceptions

include/cif++/matrix.hpp

@@ -518,7 +518,7 @@ auto eigen(const matrix_expression<M> &mat, bool sort)
 				value_type h = ev[j] - ev[i];
 				value_type t;
 
-				if (a_ij == 0 and h == 0)
+				if (h == 0 and a_ij == 0)
 					t = 1;
 				else if (std::fabs(a_ij) > 1.0e-12 * std::fabs(h))
 				{

include/cif++/point.hpp

@@ -320,7 +320,7 @@ class quaternion_type
 
 	constexpr operator bool() const
 	{
-		return operator!=({});
+		return a != 0 or b != 0 or c != 0 or d != 0;
 	}
 
   private:

include/cif++/symmetry.hpp

@@ -182,34 +182,10 @@ class spacegroup;
 class rtop;
 class sym_op;
 
-// --------------------------------------------------------------------
-// class cell
-
-class cell
-{
-  public:
-	cell(float a, float b, float c, float alpha = 90.f, float beta = 90.f, float gamma = 90.f);
-	cell(const datablock &db);
-
-	float get_a() const { return m_a; }
-	float get_b() const { return m_b; }
-	float get_c() const { return m_c; }
-
-	float get_alpha() const { return m_alpha; }
-	float get_beta() const { return m_beta; }
-	float get_gamma() const { return m_gamma; }
-
-	matrix3x3<float> get_orthogonal_matrix() const { return m_orthogonal; }
-	matrix3x3<float> get_fractional_matrix() const { return m_fractional; }
-
-  private:
-	void init();
-
-	float m_a, m_b, m_c, m_alpha, m_beta, m_gamma;
-	matrix3x3<float> m_orthogonal, m_fractional;
-};
 
 /// @brief A class that encapsulates the symmetry operations as used in PDB files, i.e. a rotational number and a translation vector
+/// The syntax in string format follows the syntax as used in mmCIF files, i.e. rotational number followed by underscore and the
+/// three translations where 5 is no movement.
 struct sym_op
 {
   public:
@@ -254,26 +230,28 @@ namespace literals
 	}
 } // namespace literals
 
+// --------------------------------------------------------------------
+// The transformation class
+
 class transformation
 {
   public:
 	transformation(const symop_data &data);
-	transformation(const matrix3x3<float> &r, const cif::point &t)
-		: m_rotation(r)
-		, m_translation(t)
-	{
-	}
+	transformation(const matrix3x3<float> &r, const cif::point &t);
 
 	transformation(const transformation &) = default;
 	transformation(transformation &&) = default;
 	transformation &operator=(const transformation &) = default;
 	transformation &operator=(transformation &&) = default;
 
-	// point operator()(const cell &c, const point &pt) const;
-
-	point operator()(const point &pt) const
+	point operator()(point pt) const
 	{
-		return m_rotation * pt + m_translation;
+		if (m_q)
+			pt.rotate(m_q);
+		else
+			pt = m_rotation * pt;
+
+		return pt + m_translation;
 	}
 
 	friend transformation operator*(const transformation &lhs, const transformation &rhs);
@@ -287,11 +265,46 @@ class transformation
 	friend class spacegroup;
 
   private:
+
+	// Most rotation matrices provided by the International Tables
+	// are really rotation matrices, in those cases we can construct
+	// a quaternion. Unfortunately, that doesn't work for all of them
+
+	void try_create_quaternion();
+
 	matrix3x3<float> m_rotation;
+	quaternion m_q;
 	point m_translation;
 };
 
 // --------------------------------------------------------------------
+// class cell
+
+class cell
+{
+  public:
+	cell(float a, float b, float c, float alpha = 90.f, float beta = 90.f, float gamma = 90.f);
+	cell(const datablock &db);
+
+	float get_a() const { return m_a; }
+	float get_b() const { return m_b; }
+	float get_c() const { return m_c; }
+
+	float get_alpha() const { return m_alpha; }
+	float get_beta() const { return m_beta; }
+	float get_gamma() const { return m_gamma; }
+
+	matrix3x3<float> get_orthogonal_matrix() const { return m_orthogonal; }
+	matrix3x3<float> get_fractional_matrix() const { return m_fractional; }
+
+  private:
+	void init();
+
+	float m_a, m_b, m_c, m_alpha, m_beta, m_gamma;
+	matrix3x3<float> m_orthogonal, m_fractional;
+};
+
+// --------------------------------------------------------------------
 
 int get_space_group_number(const datablock &db);
 int get_space_group_number(std::string_view spacegroup);
@@ -331,15 +344,50 @@ class spacegroup : public std::vector<transformation>
 };
 
 // --------------------------------------------------------------------
-// Symmetry operations on points
+// A crystal combines a cell and a spacegroup.
 
-template <typename T>
-inline point_type<T> operator*(const matrix3x3<T> &m, const point_type<T> &pt)
+class crystal
 {
-	return point_type(m(0, 0) * pt.m_x + m(0, 1) * pt.m_y + m(0, 2) * pt.m_z,
-		m(1, 0) * pt.m_x + m(1, 1) * pt.m_y + m(1, 2) * pt.m_z,
-		m(2, 0) * pt.m_x + m(2, 1) * pt.m_y + m(2, 2) * pt.m_z);
-}
+  public:
+	crystal(const datablock &db)
+		: m_cell(db)
+		, m_spacegroup(db)
+	{
+	}
+
+	crystal(const cell &c, const spacegroup &sg)
+		: m_cell(c)
+		, m_spacegroup(sg)
+	{
+	}
+
+	crystal(const crystal &) = default;
+	crystal(crystal &&) = default;
+	crystal &operator=(const crystal &) = default;
+	crystal &operator=(crystal &&) = default;
+
+	const cell &get_cell() const { return m_cell; }
+	const spacegroup &get_spacegroup() const { return m_spacegroup; }
+
+	point symmetry_copy(const point &pt, sym_op symop) const
+	{
+		return m_spacegroup(pt, m_cell, symop);
+	}
+
+	point inverse_symmetry_copy(const point &pt, sym_op symop) const
+	{
+		return m_spacegroup.inverse(pt, m_cell, symop);
+	}
+	
+	std::tuple<float,point,sym_op> closest_symmetry_copy(point a, point b) const;
+
+  private:
+	cell m_cell;
+	spacegroup m_spacegroup;
+};
+
+// --------------------------------------------------------------------
+// Symmetry operations on points
 
 inline point orthogonal(const point &pt, const cell &c)
 {
@@ -351,35 +399,6 @@ inline point fractional(const point &pt, const cell &c)
 	return c.get_fractional_matrix() * pt;
 }
 
-inline transformation orthogonal(const transformation &t, const cell &c)
-{
-	return transformation(c.get_orthogonal_matrix(), {}) * t * transformation(c.get_fractional_matrix(), {});
-}
-
-inline transformation fractional(const transformation &t, const cell &c)
-{
-	return transformation(c.get_fractional_matrix(), {}) * t * transformation(c.get_orthogonal_matrix(), {});
-}
-
 // --------------------------------------------------------------------
 
-/// @brief Return the symmetry copy of a point based on spacegroup \a sg, cell \a c and symmetry operation \a symop
-/// @param pt The point to transform
-/// @param sg The spacegroup
-/// @param c The cell
-/// @param symop The symmetry operation
-/// @return Newly calculated point
-inline point symmetry_copy(const point &pt, const spacegroup &sg, const cell &c, sym_op symop)
-{
-	return sg(pt, c, symop);
-}
-
-/// @brief Return the point and symmetry operation required to move point \a b as close as possible to point \a a
-/// @param sg The spacegroup
-/// @param c The cell
-/// @param a The point that acts as reference
-/// @param b The point that needs to be moved
-/// @return The calculated distance between the new point and \a a plus the symmetry operation required to operate on \a b
-std::tuple<float,point,sym_op> closest_symmetry_copy(const spacegroup &sg, const cell &c, point a, point b);
-
 } // namespace cif

src/symmetry.cpp

@@ -32,6 +32,8 @@
 
 #include "symop_table_data.hpp"
 
+#include <Eigen/Eigenvalues>
+
 namespace cif
 {
 
@@ -129,11 +131,60 @@ transformation::transformation(const symop_data &data)
 	m_rotation(2, 1) = d[7];
 	m_rotation(2, 2) = d[8];
 
+	try_create_quaternion();
+
 	m_translation.m_x = d[9] == 0 ? 0 : 1.0 * d[9] / d[10];
 	m_translation.m_y = d[11] == 0 ? 0 : 1.0 * d[11] / d[12];
 	m_translation.m_z = d[13] == 0 ? 0 : 1.0 * d[13] / d[14];
 }
 
+transformation::transformation(const matrix3x3<float> &r, const cif::point &t)
+	: m_rotation(r)
+	, m_translation(t)
+{
+	try_create_quaternion();
+}
+
+void transformation::try_create_quaternion()
+{
+	float Qxx = m_rotation(0, 0);
+	float Qxy = m_rotation(0, 1);
+	float Qxz = m_rotation(0, 2);
+	float Qyx = m_rotation(1, 0);
+	float Qyy = m_rotation(1, 1);
+	float Qyz = m_rotation(1, 2);
+	float Qzx = m_rotation(2, 0);
+	float Qzy = m_rotation(2, 1);
+	float Qzz = m_rotation(2, 2);
+
+	Eigen::Matrix4f em;
+
+	em << Qxx - Qyy - Qzz, Qyx + Qxy, Qzx + Qxz, Qzy - Qyz,
+			Qyx + Qxy, Qyy - Qxx - Qzz, Qzy + Qyz, Qxz - Qzx,
+			Qzx + Qxz, Qzy + Qyz, Qzz - Qxx - Qyy, Qyx - Qxy,
+			Qzy - Qyz, Qxz - Qzx, Qyx - Qxy, Qxx + Qyy + Qzz;
+
+	Eigen::EigenSolver<Eigen::Matrix4f> es(em / 3);
+
+	auto ev = es.eigenvalues();
+
+	for (size_t j = 0; j < 4; ++j)
+	{
+		if (std::abs(ev[j].real() - 1) > 0.01)
+			continue;
+		
+		auto col = es.eigenvectors().col(j);
+
+		m_q = normalize(cif::quaternion{
+			static_cast<float>(col(3).real()),
+			static_cast<float>(col(0).real()),
+			static_cast<float>(col(1).real()),
+			static_cast<float>(col(2).real()) });
+
+		break;
+	}
+}
+
 transformation operator*(const transformation &lhs, const transformation &rhs)
 {
 	auto r = lhs.m_rotation * rhs.m_rotation;
@@ -141,8 +192,6 @@ transformation operator*(const transformation &lhs, const transformation &rhs)
 	t = t + lhs.m_translation;
 
 	return transformation(r, t);
-
-	// return transformation(lhs.m_rotation * rhs.m_rotation, lhs.m_rotation * rhs.m_translation + lhs.m_translation);
 }
 
 transformation inverse(const transformation &t)
@@ -391,29 +440,29 @@ int get_space_group_number(const datablock &db)
 
 // --------------------------------------------------------------------
 
-std::tuple<float,point,sym_op> closest_symmetry_copy(const spacegroup &sg, const cell &c, point a, point b)
+std::tuple<float,point,sym_op> crystal::closest_symmetry_copy(point a, point b) const
 {
-	if (c.get_a() == 0 or c.get_b() == 0 or c.get_c() == 0)
+	if (m_cell.get_a() == 0 or m_cell.get_b() == 0 or m_cell.get_c() == 0)
 		throw std::runtime_error("Invalid cell, contains a dimension that is zero");
 
 	point result_fsb;
 	float result_d = std::numeric_limits<float>::max();
 	sym_op result_s;
 
-	auto fa = fractional(a, c);
-	auto fb = fractional(b, c);
+	auto fa = fractional(a, m_cell);
+	auto fb = fractional(b, m_cell);
 
 	auto o = offsetToOriginFractional(fa);
 
 	fa = fa + o;
 	fb = fb + o;
 
-	a = orthogonal(fa, c);
+	a = orthogonal(fa, m_cell);
 
-	for (size_t i = 0; i < sg.size(); ++i)
+	for (size_t i = 0; i < m_spacegroup.size(); ++i)
 	{
 		sym_op s(i + 1);
-		auto &t = sg[i];
+		auto &t = m_spacegroup[i];
 
 		auto fsb = t(fb);
 
@@ -453,8 +502,9 @@ std::tuple<float,point,sym_op> closest_symmetry_copy(const spacegroup &sg, const
 			s.m_tc += 1;			
 		}
 
-		auto p = orthogonal(fsb, c);
+		auto p = orthogonal(fsb, m_cell);
 		auto dsq = distance_squared(a, p);
+
 		if (result_d > dsq)
 		{
 			result_d = dsq;
@@ -463,7 +513,7 @@ std::tuple<float,point,sym_op> closest_symmetry_copy(const spacegroup &sg, const
 		}
 	}
 
-	auto p = orthogonal(result_fsb - o, c);
+	auto p = orthogonal(result_fsb - o, m_cell);
 
 	return { std::sqrt(result_d), p, result_s };
 }

test/unit-3d-test.cpp

@@ -31,6 +31,8 @@
 
 #include <cif++.hpp>
 
+#include <Eigen/Eigenvalues>
+
 namespace tt = boost::test_tools;
 namespace utf = boost::unit_test;
 
@@ -251,6 +253,72 @@ BOOST_AUTO_TEST_CASE(dh_q_1)
 
 // --------------------------------------------------------------------
 
+
+BOOST_AUTO_TEST_CASE(m2q_0, *utf::tolerance(0.001f))
+{
+	for (size_t i = 0; i < cif::kSymopNrTableSize; ++i)
+	{
+		auto d = cif::kSymopNrTable[i].symop().data();
+
+		cif::matrix3x3<float> rot;
+		float Qxx = rot(0, 0) = d[0];
+		float Qxy = rot(0, 1) = d[1];
+		float Qxz = rot(0, 2) = d[2];
+		float Qyx = rot(1, 0) = d[3];
+		float Qyy = rot(1, 1) = d[4];
+		float Qyz = rot(1, 2) = d[5];
+		float Qzx = rot(2, 0) = d[6];
+		float Qzy = rot(2, 1) = d[7];
+		float Qzz = rot(2, 2) = d[8];
+
+		Eigen::Matrix4f em;
+
+		em << Qxx - Qyy - Qzz, Qyx + Qxy, Qzx + Qxz, Qzy - Qyz,
+		      Qyx + Qxy, Qyy - Qxx - Qzz, Qzy + Qyz, Qxz - Qzx,
+			  Qzx + Qxz, Qzy + Qyz, Qzz - Qxx - Qyy, Qyx - Qxy,
+			  Qzy - Qyz, Qxz - Qzx, Qyx - Qxy, Qxx + Qyy + Qzz;
+
+		Eigen::EigenSolver<Eigen::Matrix4f> es(em / 3);
+
+		auto ev = es.eigenvalues();
+
+		size_t bestJ = 0;
+		float bestEV = -1;
+
+		for (size_t j = 0; j < 4; ++j)
+		{
+			if (bestEV < ev[j].real())
+			{
+				bestEV = ev[j].real();
+				bestJ = j;
+			}
+		}
+
+		if (std::abs(bestEV - 1) > 0.01)
+			continue; // not a rotation matrix
+
+		auto col = es.eigenvectors().col(bestJ);
+
+		auto q = normalize(cif::quaternion{
+			static_cast<float>(col(3).real()),
+			static_cast<float>(col(0).real()),
+			static_cast<float>(col(1).real()),
+			static_cast<float>(col(2).real()) });
+		
+		cif::point p1{ 1, 1, 1 };
+		cif::point p2 = p1;
+		p2.rotate(q);
+
+		cif::point p3 = rot * p1;
+
+		BOOST_TEST(p2.m_x == p3.m_x);
+		BOOST_TEST(p2.m_y == p3.m_y);
+		BOOST_TEST(p2.m_z == p3.m_z);
+	}
+}
+
+// --------------------------------------------------------------------
+
 BOOST_AUTO_TEST_CASE(symm_1)
 {
 	cif::cell c(10, 10, 10);
@@ -324,18 +392,17 @@ BOOST_AUTO_TEST_CASE(symm_4wvp_1, *utf::tolerance(0.1f))
 	auto &db = f.front();
 	cif::mm::structure s(db);
 
-	cif::spacegroup sg(db);
-	cif::cell c(db);
+	cif::crystal c(db);
 
 	cif::point p{ -78.722, 98.528,  11.994 };
 	auto a = s.get_residue("A", 10, "").get_atom_by_atom_id("O");
 
-	auto sp1 = cif::symmetry_copy(a.get_location(), sg, c, "2_565"_symop);
+	auto sp1 = c.symmetry_copy(a.get_location(), "2_565"_symop);
 	BOOST_TEST(sp1.m_x == p.m_x);
 	BOOST_TEST(sp1.m_y == p.m_y);
 	BOOST_TEST(sp1.m_z == p.m_z);
 
-	const auto &[d, sp2, so] = cif::closest_symmetry_copy(sg, c, p, a.get_location());
+	const auto &[d, sp2, so] = c.closest_symmetry_copy(p, a.get_location());
 
 	BOOST_TEST(d < 1);
 
@@ -352,8 +419,7 @@ BOOST_AUTO_TEST_CASE(symm_2bi3_1, *utf::tolerance(0.1f))
 	auto &db = f.front();
 	cif::mm::structure s(db);
 
-	cif::spacegroup sg(db);
-	cif::cell c(db);
+	cif::crystal c(db);
 
 	auto struct_conn = db["struct_conn"];
 	for (const auto &[
@@ -375,14 +441,14 @@ BOOST_AUTO_TEST_CASE(symm_2bi3_1, *utf::tolerance(0.1f))
 		auto a1 = r1.get_atom_by_atom_id(atomid1);
 		auto a2 = r2.get_atom_by_atom_id(atomid2);
 
-		auto sa1 = symmetry_copy(a1.get_location(), sg, c, cif::sym_op(symm1));
-		auto sa2 = symmetry_copy(a2.get_location(), sg, c, cif::sym_op(symm2));
+		auto sa1 = c.symmetry_copy(a1.get_location(), cif::sym_op(symm1));
+		auto sa2 = c.symmetry_copy(a2.get_location(), cif::sym_op(symm2));
 
 		BOOST_TEST(cif::distance(sa1, sa2) == dist);
 
 		auto pa1 = a1.get_location();
 
-		const auto &[d, p, so] = cif::closest_symmetry_copy(sg, c, pa1, a2.get_location());
+		const auto &[d, p, so] = c.closest_symmetry_copy(pa1, a2.get_location());
 
 		BOOST_TEST(p.m_x == sa2.m_x);
 		BOOST_TEST(p.m_y == sa2.m_y);
@@ -400,8 +466,7 @@ BOOST_AUTO_TEST_CASE(symm_3bwh_1, *utf::tolerance(0.1f))
 
 	auto &db = f.front();
 
-	cif::spacegroup sg(db);
-	cif::cell c(db);
+	cif::crystal c(db);
 	cif::mm::structure s(db);
 
 	for (auto a1 : s.atoms())
@@ -411,7 +476,7 @@ BOOST_AUTO_TEST_CASE(symm_3bwh_1, *utf::tolerance(0.1f))
 			if (a1 == a2)
 				continue;
 			
-			const auto&[ d, p, so ] = cif::closest_symmetry_copy(sg, c, a1.get_location(), a2.get_location());
+			const auto&[ d, p, so ] = c.closest_symmetry_copy(a1.get_location(), a2.get_location());
 
 			BOOST_TEST(d == distance(a1.get_location(), p));
 		}