for now, require eigen3

add inverse_symmetry_copy to crystal

CMakeLists.txt

@@ -168,7 +168,7 @@ if(MSVC)
 endif()
 find_package(ZLIB REQUIRED)
 
-find_package(Eigen3)
+find_package(Eigen3 )
 
 include(FindFilesystem)
 list(APPEND CIFPP_REQUIRED_LIBRARIES ${STDCPPFS_LIBRARY})
@@ -279,9 +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 PRIVATE Eigen3::Eigen)
+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)

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

@@ -237,23 +237,22 @@ 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 point &pt) const;
+	point operator()(point pt) const
+	{
+		if (m_q)
+			pt.rotate(m_q);
+		else
+			pt = m_rotation * pt;
 
-	// point operator()(const point &pt) const
-	// {
-	// 	return m_rotation * pt + m_translation;
-	// }
+		return pt + m_translation;
+	}
 
 	friend transformation operator*(const transformation &lhs, const transformation &rhs);
 	friend transformation inverse(const transformation &t);
@@ -266,6 +265,13 @@ 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;
@@ -367,6 +373,11 @@ class crystal
 	{
 		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;
 

src/symmetry.cpp

@@ -121,15 +121,41 @@ transformation::transformation(const symop_data &data)
 {
 	const auto &d = data.data();
 
-	float Qxx = m_rotation(0, 0) = d[0];
-	float Qxy = m_rotation(0, 1) = d[1];
-	float Qxz = m_rotation(0, 2) = d[2];
-	float Qyx = m_rotation(1, 0) = d[3];
-	float Qyy = m_rotation(1, 1) = d[4];
-	float Qyz = m_rotation(1, 2) = d[5];
-	float Qzx = m_rotation(2, 0) = d[6];
-	float Qzy = m_rotation(2, 1) = d[7];
-	float Qzz = m_rotation(2, 2) = d[8];
+	m_rotation(0, 0) = d[0];
+	m_rotation(0, 1) = d[1];
+	m_rotation(0, 2) = d[2];
+	m_rotation(1, 0) = d[3];
+	m_rotation(1, 1) = d[4];
+	m_rotation(1, 2) = d[5];
+	m_rotation(2, 0) = d[6];
+	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;
 
@@ -157,22 +183,6 @@ transformation::transformation(const symop_data &data)
 
 		break;
 	}
-
-	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];
-}
-
-point transformation::operator()(const point &pt) const
-{
-	cif::point p = pt;
-
-	if (m_q)
-		p.rotate(m_q);
-	else
-		p = m_rotation * pt;
-		
-	return p + m_translation;
 }
 
 transformation operator*(const transformation &lhs, const transformation &rhs)
@@ -182,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)