symopnr

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

include/cif++/Symmetry.h

@@ -12,8 +12,6 @@ namespace mmcif
 
 clipper::Coord_orth CalculateOffsetForCell(const Structure& p, const clipper::Spacegroup& spacegroup, const clipper::Cell& cell);
 std::vector<clipper::RTop_orth> AlternativeSites(const clipper::Spacegroup& spacegroup, const clipper::Cell& cell);
-int32_t GetRotationalIndexNumber(const clipper::Spacegroup& spacegroup, const clipper::Cell& cell,
-	const clipper::RTop_orth& rt);
 
 // --------------------------------------------------------------------
 // To iterate over all near symmetry copies of an atom

src/Symmetry.cpp

@@ -130,6 +130,70 @@ vector<clipper::RTop_orth> AlternativeSites(const clipper::Spacegroup& spacegrou
 	return result;
 }
 
+// --------------------------------------------------------------------
+// Unfortunately, clipper has a different numbering scheme than PDB
+// for rotation numbers. So we created a table to map those.
+// Perhaps a bit over the top, but hey....
+
+#include "SymOpTable_data.cpp"
+
+int32_t GetRotationalIndexNumber(int spacegroup, const clipper::RTop_frac& rt)
+{
+	auto& rot = rt.rot();
+	auto& trn = rt.trn();
+
+	auto rte = [&rot](int i, int j) { return static_cast<int8_t>(lrint(rot(i, j))); };
+
+	SymopData k
+	{
+		rte(0, 0), rte(0, 1), rte(0, 2),
+		rte(1, 0), rte(1, 1), rte(1, 2),
+		rte(2, 0), rte(2, 1), rte(2, 2)
+	};
+
+	for (int i = 0; i < 3; ++i)
+	{
+		int n = lrint(trn[i] * 24);
+		int d = 24;
+
+		if (n == 0 or n == 24)
+			continue;		// is 0, 0 in our table
+
+		for (int i = 5; i > 1; --i)
+			if (n % i == 0 and d % i == 0)
+			{
+				n /= i;
+				d /= i;
+			}
+
+		switch (i)
+		{
+			case 0: k.trn_0_0 = n; k.trn_0_1 = d; break;
+			case 1: k.trn_1_0 = n; k.trn_1_1 = d; break;
+			case 2: k.trn_2_0 = n; k.trn_2_1 = d; break;
+		}
+	}
+
+	const size_t N = sizeof(kSymopNrTable) / sizeof(SymopDataBlock);
+	size_t L = 0, R = N - 1;
+	while (L <= R)
+	{
+		size_t i = (L + R) / 2;
+		if (kSymopNrTable[i].spacegroupNr < spacegroup)
+			L = i + 1;
+		else
+			R = i - 1;
+	}
+
+	for (size_t i = L; i < N and kSymopNrTable[i].spacegroupNr == spacegroup; ++i)
+	{
+		if (kSymopNrTable[i].rt.iv == k.iv)
+			return kSymopNrTable[i].rotationalNr;
+	}
+
+	throw runtime_error("Symmetry operation was not found in table, cannot find rotational number");
+}
+
 // -----------------------------------------------------------------------
 
 string SymmetryAtomIteratorFactory::symop_mmcif(const Atom& a) const
@@ -163,6 +227,8 @@ string SymmetryAtomIteratorFactory::symop_mmcif(const Atom& a) const
 
 							auto rtop_f = rtop.rtop_frac(mCell);
 
+							int rnr = GetRotationalIndexNumber(mSpacegroup.spacegroup_number(), rtop_f);
+
 							uint32_t t[3] = {
 								static_cast<uint32_t>(5 + static_cast<int>(rint(rtop_f.trn()[0]))),
 								static_cast<uint32_t>(5 + static_cast<int>(rint(rtop_f.trn()[1]))),
@@ -172,7 +238,7 @@ string SymmetryAtomIteratorFactory::symop_mmcif(const Atom& a) const
 							if (t[0] > 9 or t[1] > 9 or t[2] > 9)
 								throw runtime_error("Symmetry operation has an out-of-range translation.");
 
-							result += to_string(1 + i) + "_"
+							result += to_string(rnr) + "_"
 								   + to_string(t[0])
 								   + to_string(t[1])
 								   + to_string(t[2]);

tools/symop-map-generator.cpp

@@ -14,8 +14,6 @@ clang++ -I ~/my-clipper/include -L ~/my-clipper/lib -o symop-map-generator symop
 
 #include <cstdlib>
 
-#include <clipper/clipper.h>
-
 using namespace std;
 
 std::regex kNameRx(R"(^(\d+) +(\d+) +(\d+) +(\S+) +(\S+) +(\S+) +'([^']+)'( +'([^']+)')?(?: +!.+)?$)");
@@ -172,26 +170,20 @@ int main()
 		ifstream file(CLIBD + "/symop.lib"s);
 		if (not file.is_open())
 			throw runtime_error("Could not open symop.lib file");
-		
-		string line;
-		clipper::Spacegroup spacegroup;
 
-		enum class State { Initial, InSpacegroup, Error } state = State::Initial;
-		int symopnr = 0;
+		// --------------------------------------------------------------------
+		// unfortunately, the data in symop.lib is not sorted, so we have to
+		// store the data in memory before writing out.
 
-		cout << R"(
-// This file was generated from $CLIBD/symop.lib
+		vector<tuple<int,int,string,array<int,15>>> data;
 
-struct SymopNrTable
-{
-	uint8_t spacegroupNr;
-	uint8_t	rotationalNr:7;
-	int8_t rot[3][3]:2;
-	uint8_t trn[3][2]:4;
-};
+		// --------------------------------------------------------------------
+		
+		string line;
+		string spacegroupName;
+		int spacegroupNr, symopnr;
 
-kSymopNrTable[] = {
-)" << endl;
+		enum class State { Initial, InSpacegroup, Error } state = State::Initial;
 
 		while (getline(file, line))
 		{
@@ -210,12 +202,8 @@ kSymopNrTable[] = {
 						try
 						{
 							SymopParser p;
-							auto symop = p.parse(line);
-
-							cout << "    { " << setw(3) << spacegroup.spacegroup_number() << '\t' << symopnr++;
-							for (auto i: symop)
-								cout << ',' << setw(2) << i;
-							cout << " }," << endl;
+							data.emplace_back(spacegroupNr, symopnr, spacegroupName, p.parse(line));
+							++symopnr;
 						}
 						catch (const exception& e)
 						{
@@ -223,7 +211,6 @@ kSymopNrTable[] = {
 								 << e.what() << endl;
 						}
 
-
 						continue;
 					}
 					// fall through
@@ -237,40 +224,73 @@ kSymopNrTable[] = {
 						throw runtime_error("Name line does not match regular expression");
 					}
 
-					int nr = stoi(m[1].str());
-					if (nr > 230)
-					{
-						state = State::Error;
-						continue;
-					}
-
-					try
-					{
-						spacegroup = clipper::Spacegroup(clipper::Spgr_descr(stoi(m[1].str())));
-						symopnr = 1;
-						
-						cout << "    // " << spacegroup.symbol_hm() << endl;
+					spacegroupNr = stoi(m[1]);
+					spacegroupName = m[7];
+					symopnr = 1;
 
-						if (spacegroup.num_symops() != stoi(m[2]))
-						{
-							cerr << line << endl
-								 << "Num of symops do not match: " << spacegroup.num_symops() << " <> " << m[2] << endl;
-							state = State::Error;
-							continue;
-						}
-
-						state = State::InSpacegroup;
-					}
-					catch (const clipper::Message_fatal& e)
-					{
-						cerr << line << endl
-							 << e.text() << endl;
-						state = State::Error;
-					}
+					state = State::InSpacegroup;
 					break;
 				}
 			}
 		}
+
+		// --------------------------------------------------------------------
+
+		sort(data.begin(), data.end());
+
+		// --------------------------------------------------------------------
+
+		cout << R"(// This file was generated from $CLIBD/symop.lib
+
+union SymopData
+{
+	struct
+	{
+		int rot_0_0:2;
+		int rot_0_1:2;
+		int rot_0_2:2;
+		int rot_1_0:2;
+		int rot_1_1:2;
+		int rot_1_2:2;
+		int rot_2_0:2;
+		int rot_2_1:2;
+		int rot_2_2:2;
+		unsigned int trn_0_0:3;
+		unsigned int trn_0_1:3;
+		unsigned int trn_1_0:3;
+		unsigned int trn_1_1:3;
+		unsigned int trn_2_0:3;
+		unsigned int trn_2_1:3;
+	};
+	uint64_t iv:36;
+};
+
+struct SymopDataBlock
+{
+	uint16_t spacegroupNr;
+	uint8_t rotationalNr;
+	SymopData rt;
+} kSymopNrTable[] = {
+)" << endl;
+
+		spacegroupNr = 0;
+		for (auto& sd: data)
+		{
+			int sp, o;
+			string n;
+			tie(sp, o, n, ignore) = sd;
+
+			if (sp > spacegroupNr)
+				cout << "    // " << n << endl;
+			spacegroupNr = sp;
+
+			cout << "    { " << setw(3) << sp
+					<< ", " << setw(3) << o << ", { ";
+			for (auto i: get<3>(sd))
+				cout << setw(2) << i << ',';
+			cout << " } }," << endl;
+		}
+
 		cout << "};" << endl;
 	}
 	catch (const exception& ex)