backup

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

include/cif++/Symmetry.h

@@ -13,6 +13,7 @@ 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);
 int GetSpacegroupNumber(const std::string& spacegroup);	// alternative for clipper's parsing code
+std::string SpacegroupToxHM(std::string spacegroup);
 
 // --------------------------------------------------------------------
 // To iterate over all near symmetry copies of an atom
@@ -21,10 +22,19 @@ class SymmetryAtomIteratorFactory
 {
   public:
 	// SymmetryAtomIteratorFactory(const Structure& p);
+
 	SymmetryAtomIteratorFactory(const Structure& p, const clipper::Spacegroup& spacegroup, const clipper::Cell& cell)
-		: mD(CalculateOffsetForCell(p, spacegroup, cell))
-		, mRtOrth(AlternativeSites(spacegroup, cell))
+		: mSpacegroupNr(spacegroup.spacegroup_number())
 		, mSpacegroup(spacegroup)
+		, mD(CalculateOffsetForCell(p, spacegroup, cell))
+		, mRtOrth(AlternativeSites(spacegroup, cell))
+		, mCell(cell) {}
+
+	SymmetryAtomIteratorFactory(const Structure& p, int spacegroupNr, const clipper::Cell& cell)
+		: mSpacegroupNr(spacegroupNr)
+		, mSpacegroup(GetCCP4SpacegroupDescr(spacegroupNr))
+		, mD(CalculateOffsetForCell(p, mSpacegroup, cell))
+		, mRtOrth(AlternativeSites(mSpacegroup, cell))
 		, mCell(cell) {}
 
 	SymmetryAtomIteratorFactory(const SymmetryAtomIteratorFactory&) = delete;
@@ -123,9 +133,10 @@ class SymmetryAtomIteratorFactory
 	std::string symop_mmcif(const Atom& a) const;
 
   private:
+	int										mSpacegroupNr;
+	clipper::Spacegroup						mSpacegroup;
 	Point									mD;			// needed to move atoms to center
 	std::vector<clipper::RTop_orth>			mRtOrth;
-	clipper::Spacegroup						mSpacegroup;
 	clipper::Cell							mCell;
 };
 

src/Symmetry.cpp

@@ -237,6 +237,57 @@ int GetSpacegroupNumber(const std::string& spacegroup)
 
 // -----------------------------------------------------------------------
 
+std::string SpacegroupToxHM(std::string spacegroup)
+{
+	if (spacegroup == "P 1-")
+		spacegroup = "P -1";
+	else if (spacegroup == "P 21 21 2 A")
+		spacegroup = "P 21 21 2 (a)";
+	else if (spacegroup.empty())
+		throw runtime_error("No spacegroup, cannot continue");
+
+	int nr = GetSpacegroupNumber(spacegroup);
+
+	return kSpaceGroups[nr].xHM;
+}
+
+clipper::Spgr_descr GetCCP4SpacegroupDescr(int nr)
+{
+	const size_t N = sizeof(kSymopNrTable) / sizeof(SymopDataBlock);
+	int32_t L = 0, R = static_cast<int32_t>(N - 1);
+	while (L <= R)
+	{
+		int32_t i = (L + R) / 2;
+		if (kSymopNrTable[i].spacegroupNr < nr)
+			L = i + 1;
+		else
+			R = i - 1;
+	}
+
+	if (L < 0 or L >= N)
+		throw runtime_error("Invalid spacegroup number");
+
+	clipper::Spgr_descr::Symop_codes symop_codes;
+
+	for (size_t i = L; i < N and kSymopNrTable[i].spacegroupNr == nr; ++i)
+	{
+		auto& symop = kSymopNrTable[i];
+
+		clipper::ftype m[4][4] = {
+			{ symop.rt.rot_0_0, symop.rt.rot_0_1, symop.rt.rot_0_2, static_cast<float>(symop.rt.trn_0_0) / symop.rt.trn_0_1 },
+			{ symop.rt.rot_1_0, symop.rt.rot_1_1, symop.rt.rot_1_2, static_cast<float>(symop.rt.trn_1_0) / symop.rt.trn_1_1 },
+			{ symop.rt.rot_2_0, symop.rt.rot_2_1, symop.rt.rot_2_2, static_cast<float>(symop.rt.trn_2_0) / symop.rt.trn_2_1 },
+			{ 0, 0, 0, 1 }
+		};
+
+		symop_codes.emplace_back(clipper::Symop(m));
+	}
+
+	return clipper::Spgr_descr(symop_codes);
+}
+
+// -----------------------------------------------------------------------
+
 string SymmetryAtomIteratorFactory::symop_mmcif(const Atom& a) const
 {
 	string result;
@@ -268,7 +319,7 @@ string SymmetryAtomIteratorFactory::symop_mmcif(const Atom& a) const
 
 							auto rtop_f = rtop.rtop_frac(mCell);
 
-							int rnr = GetRotationalIndexNumber(mSpacegroup.spacegroup_number(), rtop_f);
+							int rnr = GetRotationalIndexNumber(mSpacegroupNr, rtop_f);
 
 							uint32_t t[3] = {
 								static_cast<uint32_t>(5 + static_cast<int>(rint(rtop_f.trn()[0]))),

tools/symop-map-generator.cpp

@@ -11,6 +11,7 @@ clang++ -I ~/my-clipper/include -L ~/my-clipper/lib -o symop-map-generator symop
 #include <iomanip>
 #include <fstream>
 #include <regex>
+#include <map>
 
 #include <cstdlib>
 
@@ -156,6 +157,60 @@ class SymopParser
 	int m_trn[3][2] = {};
 };
 
+map<string,string> get_xhm_map(string path)
+{
+	ifstream file(path);
+	if (not file.is_open())
+		throw runtime_error("Could not open syminfo.lib file");
+
+	enum class State { skip, spacegroup } state = State::skip;
+
+	string line;
+	string old, xHM;
+
+	const regex rx(R"(^symbol +(xHM|old) +'(.*)'$)");
+
+	map<string,string> result;
+
+	while (getline(file, line))
+	{
+		switch (state)
+		{
+			case State::skip:
+				if (line == "begin_spacegroup")
+					state = State::spacegroup;
+				break;
+			
+			case State::spacegroup:
+			{
+				smatch m;
+				if (regex_match(line, m, rx))
+				{
+cerr << "match " << m[1] << " => " << m[2] << endl;
+
+					if (m[1] == "old")
+						old = m[2];
+					else if (m[1] == "xHM")
+						xHM = m[2];
+				}
+				else if (line == "end_spacegroup")
+				{
+					if (not (old.empty() or xHM.empty()))
+						result.emplace(old, xHM);
+					
+					state = State::skip;
+					old.clear();
+					xHM.clear();
+				}
+				break;
+			}
+		}
+	}
+
+	return result;
+}
+
+
 int main()
 {
 	try
@@ -236,6 +291,10 @@ int main()
 			}
 		}
 
+		// -----------------------------------------------------------------------
+		
+		auto xhm_map = get_xhm_map(CLIBD + "/syminfo.lib"s);
+
 		// --------------------------------------------------------------------
 
 		sort(data.begin(), data.end());
@@ -248,13 +307,20 @@ int main()
 struct Spacegroup
 {
 	const char* name;
+	const char* xHM;
 	int nr;
 } kSpaceGroups[] =
 {
 )";
 
-		for (auto& [name, nr]: spacegroups)
-			cout << "\t{ \"" << name << "\", " << nr << " }," << endl;
+		for (auto [name, nr]: spacegroups)
+		{
+			name = '"' + name + '"' + string(20 - name.length(), ' ');
+			string xHM = xhm_map[name];
+			xHM = '"' + xHM + '"' + string(20 - xHM.length(), ' ');
+
+			cout << "\t{ " << name << ", " << xHM << ", " << nr << " }," << endl;
+		}
 
 cout << R"(
 };