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

include/cif++/model.hpp

@@ -667,6 +667,8 @@ class branch : public std::vector<sugar>
 	}
 
 	sugar &construct_sugar(const std::string &compound_id);
+	sugar &construct_sugar(const std::string &compound_id, const std::string &atom_id,
+		int linked_sugar_nr, const std::string &linked_atom_id);
 
   private:
 	friend sugar;

include/cif++/point.hpp

@@ -54,7 +54,7 @@ class quaternion_type
   public:
 	using value_type = T;
 
-	constexpr explicit quaternion_type(value_type const &value_a = value_type(), value_type const &value_b = value_type(), value_type const &value_c = value_type(), value_type const &value_d = value_type())
+	constexpr explicit quaternion_type(value_type const &value_a = {}, value_type const &value_b = {}, value_type const &value_c = {}, value_type const &value_d = {})
 		: a(value_a)
 		, b(value_b)
 		, c(value_c)
@@ -305,6 +305,21 @@ class quaternion_type
 	constexpr value_type get_c() const { return c; }
 	constexpr value_type get_d() const { return d; }
 
+	constexpr bool operator==(const quaternion_type &rhs) const
+	{
+		return a == rhs.a and b == rhs.b and c == rhs.c and d == rhs.d;
+	}
+
+	constexpr bool operator!=(const quaternion_type &rhs) const
+	{
+		return a != rhs.a or b != rhs.b or c != rhs.c or d != rhs.d;
+	}
+
+	constexpr operator bool() const
+	{
+		return operator!=({});
+	}
+
   private:
 	value_type a, b, c, d;
 };
@@ -483,6 +498,13 @@ struct point_type
 		m_z = p.get_d();
 	}
 
+	constexpr void rotate(const quaternion &q, point_type pivot)
+	{
+		operator-=(pivot);
+		rotate(q);
+		operator+=(pivot);
+	}
+
 #if HAVE_LIBCLIPPER
 	operator clipper::Coord_orth() const
 	{
@@ -665,6 +687,12 @@ point nudge(point p, float offset);
 quaternion construct_from_angle_axis(float angle, point axis);
 std::tuple<double, point> quaternion_to_angle_axis(quaternion q);
 
+/// @brief Given four points and an angle, return the quaternion required to rotate
+/// point p4 along the p2-p3 axis and around point p3 to obtain the required within
+/// an accuracy of esd
+quaternion construct_for_dihedral_angle(point p1, point p2, point p3, point p4,
+	float angle, float esd);
+
 point centroid(const std::vector<point> &Points);
 point center_points(std::vector<point> &Points);
 

src/model.cpp

@@ -1244,6 +1244,36 @@ sugar &branch::construct_sugar(const std::string &compound_id)
 	return result;
 }
 
+sugar &branch::construct_sugar(const std::string &compound_id, const std::string &atom_id,
+	int linked_sugar_nr, const std::string &linked_atom_id)
+{
+	auto &result = construct_sugar(compound_id);
+
+	auto &linked = get_sugar_by_num(linked_sugar_nr);
+	result.set_link(linked.get_atom_by_atom_id(linked_atom_id));
+
+	auto &db = m_structure->get_datablock();
+
+	auto &pdbx_entity_branch_link = db["pdbx_entity_branch_link"];
+	auto linkID = pdbx_entity_branch_link.get_unique_id("");
+
+	db["pdbx_entity_branch_link"].emplace({
+		{ "link_id", linkID },
+		{ "entity_id", get_entity_id() },
+		{ "entity_branch_list_num_1", result.num() }, 
+		{ "comp_id_1", compound_id }, 
+		{ "atom_id_1", atom_id },
+		{ "leaving_atom_id_1", "O1" }, 	/// TODO: Need to fix this!
+		{ "entity_branch_list_num_2", linked.num() }, 
+		{ "comp_id_2", linked.get_compound_id() }, 
+		{ "atom_id_2", linked_atom_id }, 
+		{ "leaving_atom_id_2", "." }, 
+		{ "value_order", "sing" }
+	});
+
+	return result;
+}
+
 std::string branch::name(const sugar &s) const
 {
 	using namespace literals;
@@ -2275,10 +2305,16 @@ std::string structure::create_non_poly(const std::string &entity_id, std::vector
 
 branch &structure::create_branch()
 {
+	auto &entity = m_db["entity"];
 	auto &struct_asym = m_db["struct_asym"];
 
-	std::string asym_id = struct_asym.get_unique_id();
-	auto entity_id = m_db["entity"].get_unique_id("");
+	auto entity_id = entity.get_unique_id("");
+	auto asym_id = struct_asym.get_unique_id();
+
+	entity.emplace({
+		{"id", entity_id},
+		{"type", "branched"}
+	});
 
 	struct_asym.emplace({
 		{"id", asym_id},

src/point.cpp

@@ -355,6 +355,44 @@ point center_points(std::vector<point> &Points)
 	return t;
 }
 
+quaternion construct_for_dihedral_angle(point p1, point p2, point p3, point p4,
+	float angle, float esd)
+{
+	p1 -= p3;
+	p2 -= p3;
+	p4 -= p3;
+	p3 -= p3;
+
+	quaternion q;
+	auto axis = p2;
+
+	float dh = dihedral_angle(p1, p2, p3, p4);
+	for (int iteration = 0; iteration < 100; ++iteration)
+	{
+		float delta = std::fmod(angle - dh, 360);
+
+		if (delta < -180)
+			delta += 360;
+		if (delta > 180)
+			delta -= 360;
+
+		if (std::abs(delta) < esd)
+			break;
+
+		// if (iteration > 0)
+		// 	std::cout << cif::coloured(("iteration " + std::to_string(iteration)).c_str(), cif::scBLUE, cif::scBLACK) << " delta: " << delta << std::endl;
+
+		auto q2 = construct_from_angle_axis(delta, axis);
+		q = iteration == 0 ? q2 : q * q2;
+			
+		p4.rotate(q2);
+
+		dh = dihedral_angle(p1, p2, p3, p4);
+	}
+
+	return q;
+}
+
 point centroid(const std::vector<point> &pts)
 {
 	point result;

src/utilities.cpp

@@ -171,13 +171,15 @@ struct ProgressImpl
 
 void ProgressImpl::Run()
 {
+	using namespace std::literals;
+
 	bool printedAny = false;
 
 	try
 	{
 		for (;;)
 		{
-			std::this_thread::sleep_for(std::chrono::milliseconds(100));
+			std::this_thread::sleep_for(2s);
 
 			std::unique_lock lock(mMutex);
 

test/unit-3d-test.cpp

@@ -87,19 +87,20 @@ BOOST_AUTO_TEST_CASE(t1)
 	// q = Normalize(q);
 
 	// Quaternion q{ 0.1, 0.2, 0.3, 0.4 };
-	cif::quaternion q{0.5, 0.5, 0.5, 0.5};
+	cif::quaternion q{ 0.5, 0.5, 0.5, 0.5 };
 	q = normalize(q);
 
 	const auto &&[angle0, axis0] = cif::quaternion_to_angle_axis(q);
 
 	std::vector<cif::point> p1{
-		{16.979, 13.301, 44.555},
-		{18.150, 13.525, 43.680},
-		{18.656, 14.966, 43.784},
-		{17.890, 15.889, 44.078},
-		{17.678, 13.270, 42.255},
-		{16.248, 13.734, 42.347},
-		{15.762, 13.216, 43.724}};
+		{ 16.979, 13.301, 44.555 },
+		{ 18.150, 13.525, 43.680 },
+		{ 18.656, 14.966, 43.784 },
+		{ 17.890, 15.889, 44.078 },
+		{ 17.678, 13.270, 42.255 },
+		{ 16.248, 13.734, 42.347 },
+		{ 15.762, 13.216, 43.724 }
+	};
 
 	auto p2 = p1;
 
@@ -136,7 +137,7 @@ BOOST_AUTO_TEST_CASE(t2)
 
 	cif::point xp = cif::cross_product(p[1] - p[0], p[2] - p[0]);
 
-	auto q = cif::construct_from_angle_axis(45, xp); //mmcif::Normalize(Quaternion{45 * mmcif::kPI / 180, xp.mX, xp.mY, xp.mZ});
+	auto q = cif::construct_from_angle_axis(45, xp); // mmcif::Normalize(Quaternion{45 * mmcif::kPI / 180, xp.mX, xp.mY, xp.mZ});
 
 	auto &&[angle, axis] = cif::quaternion_to_angle_axis(q);
 
@@ -153,7 +154,7 @@ BOOST_AUTO_TEST_CASE(t3)
 
 	cif::point xp = cif::cross_product(p[1] - p[0], p[2] - p[0]);
 
-	auto q = cif::construct_from_angle_axis(45, xp); //mmcif::Normalize(Quaternion{45 * mmcif::kPI / 180, xp.mX, xp.mY, xp.mZ});
+	auto q = cif::construct_from_angle_axis(45, xp); // mmcif::Normalize(Quaternion{45 * mmcif::kPI / 180, xp.mX, xp.mY, xp.mZ});
 
 	auto v = p[1];
 	v -= p[0];
@@ -166,3 +167,48 @@ BOOST_AUTO_TEST_CASE(t3)
 
 	BOOST_TEST(a == 45, tt::tolerance(0.01));
 }
+
+BOOST_AUTO_TEST_CASE(dh_q_1)
+{
+	struct
+	{
+		float angle;
+		cif::point pts[4];
+	} tests[] = {
+		{ -97.5,
+			{ { 68.8649979, -7.34800005, 54.3769989 },
+				{ 68.1350021, -8.18700027, 53.6489983 },
+				{ 68.7760239, -9.07335377, 52.7140236 },
+				{ 68.9000015, -10.3944235, 53.2217026 } } },
+		{ 80.3,
+			{ { 0.304512024, 0.531184196, 2.25860214 },
+				{ 0.956512451, 0.0321846008, 1.07460022 },
+				{ 0, 0, 0 },
+				{ 0.21336633, -1.09552193, -0.878999829 } } },
+		{ -97.5,
+			{ { 0.088973999, 1.72535372, 1.66297531 },
+				{ -0.641021729, 0.886353493, 0.93497467 },
+				{ 0, 0, 0 },
+				{ 1.29433727, -0.395142615, 0.432300746 } } },
+		{ -97.5,
+			{
+				{ 0.088973999, 1.72535372, 1.66297531 },
+				{ -0.641021729, 0.886353493, 0.93497467 },
+				{ 0, 0, 0 },
+				{ 1.33983064, 0.384027064, -0.275154471 },
+
+			} }
+	};
+
+	for (auto &&[angle, pts] : tests)
+	{
+		auto q = cif::construct_for_dihedral_angle(pts[0], pts[1], pts[2], pts[3], angle, 1);
+
+		pts[3] -= pts[2];
+		pts[3].rotate(q);
+		pts[3] += pts[2];
+
+		auto dh = cif::dihedral_angle(pts[0], pts[1], pts[2], pts[3]);
+		BOOST_TEST(dh == angle, tt::tolerance(0.1f));
+	}
+}
\ No newline at end of file