Introduced cif::crystal

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,6 +230,9 @@ namespace literals
 	}
 } // namespace literals
 
+// --------------------------------------------------------------------
+// The transformation class
+
 class transformation
 {
   public:
@@ -269,8 +248,6 @@ class transformation
 	transformation &operator=(const transformation &) = default;
 	transformation &operator=(transformation &&) = default;
 
-	// point operator()(const cell &c, const point &pt) const;
-
 	point operator()(const point &pt) const
 	{
 		return m_rotation * pt + m_translation;
@@ -292,6 +269,33 @@ class transformation
 };
 
 // --------------------------------------------------------------------
+// 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 +335,45 @@ 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);
+	}
+	
+	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 +385,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

@@ -391,29 +391,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 +453,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 +464,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

@@ -324,18 +324,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 +351,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 +373,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 +398,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 +408,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));
 		}