space group test

tools/cif2xml/CMakeLists.txt

@@ -13,7 +13,7 @@ set(CMAKE_VERBOSE_MAKEFILE TRUE)
 find_package(Boost REQUIRED)
 
 
-set(CMAKE_CXX_STANDARD 20)
+set(CMAKE_CXX_STANDARD 17)
 add_definitions(-std=c++2a -fconcepts)
 add_definitions(${Boost_CXX_FLAGS})
 
@@ -27,6 +27,11 @@ include_directories(
 
 add_executable(cif2xml
 	cif2xml.cpp ../structfact/loadcif.h)
-target_link_libraries(cif2xml 
+target_link_libraries(cif2xml
 #	-static-libstdc++ -static-libgcc
 )
+
+
+add_executable(findsg
+	findsg.cpp ../structfact/loadcif.h)
+target_link_libraries(findsg)

tools/cif2xml/findsg.cpp

+/**
+ * finds a matching space group
+ * @author Tobias Weber <tweber@ill.fr>
+ * @date jan-2020
+ * @license GPLv3, see 'LICENSE' file
+ */
+
+#include "../structfact/loadcif.h"
+#include "libs/_cxx20/math_algos.h"
+
+#include <gemmi/version.hpp>
+#include <iostream>
+#include <memory>
+
+
+using t_real = double;
+using t_vec = std::vector<t_real>;
+using t_mat = m::mat<t_real, std::vector>;
+
+constexpr t_real g_eps = 1e-6;
+constexpr int g_prec = 6;
+
+
+
+/**
+ * find matching spacegroup
+ */
+std::vector<std::tuple<int, std::string, std::vector<t_mat>>> find_sgs(
+	const std::vector<t_vec>& vecInit, const std::vector<t_vec>& vecFinal)
+{
+	return find_matching_sgs<t_vec, t_mat, t_real>(vecInit, vecFinal);
+}
+
+
+
+/**
+ * entry point
+ */
+int main(int argc, char** argv)
+{
+	// TODO: read input data
+
+	// test data
+	std::vector<t_vec> vecInit
+	{{
+		m::create<t_vec>({ 0.1,  0.1,  0.1}),
+		m::create<t_vec>({ 0.4, -0.1, -0.4}),
+		m::create<t_vec>({-0.1, -0.4,  0.4}),
+		m::create<t_vec>({-0.4,  0.4, -0.1}),
+	}};
+	std::vector<t_vec> vecFinal
+	{{
+		m::create<t_vec>({ 0.1,  0.1,  0.1}),
+		m::create<t_vec>({ 0.4, -0.1, -0.4}),
+		m::create<t_vec>({-0.1, -0.4,  0.4}),
+		m::create<t_vec>({-0.4,  0.4, -0.1}),
+	}};
+
+
+	std::cout << "Full set of positions to match:\n";
+	std::size_t ctr = 1;
+	for(const auto& pos : vecFinal)
+		std::cout << "\t(" << ctr++ << ") " << pos << "\n";
+	std::cout << std::endl;
+
+	while(1)
+	{
+		std::cout << "\n--------------------------------------------------------------------------------\n";
+		std::cout << "Base set of positions:\n";
+		ctr = 1;
+		for(const auto& pos : vecInit)
+			std::cout << "\t(" << ctr++ << ") " << pos << "\n";
+		std::cout << std::endl;
+
+		auto matchingSGs = find_sgs(vecInit, vecFinal);
+
+		std::cout << "Matching space groups:\n";
+		ctr = 1;
+		for(const auto& sg : matchingSGs)
+			std::cout << "\t(" << ctr++ << ") " << std::get<1>(sg) << "\n";
+		std::cout << "--------------------------------------------------------------------------------\n";
+		std::cout << std::endl;
+
+		vecInit.pop_back();
+		if(vecInit.size() == 0)
+			break;
+	}
+
+	return 0;
+}

tools/moldyn/CMakeLists.txt

@@ -21,7 +21,7 @@ find_package(Qt5 REQUIRED COMPONENTS Core Gui Widgets OpenGL)
 set(CMAKE_AUTOUIC TRUE)
 set(CMAKE_AUTOMOC TRUE)
 
-set(CMAKE_CXX_STANDARD 20)
+set(CMAKE_CXX_STANDARD 17)
 add_definitions(-std=c++2a -fconcepts)
 add_definitions(${Boost_CXX_FLAGS})
 add_definitions(-DBOOST_SYSTEM_NO_DEPRECATED -DBOOST_ERROR_CODE_HEADER_ONLY)

tools/structfact/loadcif.h

@@ -264,7 +264,11 @@ load_cif(const std::string& filename, t_real eps=1e-6)
  * gets space group description strings and symmetry operations
  */
 template<class t_mat, class t_real = typename t_mat::value_type>
-std::vector<std::tuple<int, std::string, std::vector<t_mat>>>
+std::vector<std::tuple<
+	int,				// sg number 
+	std::string, 		// description
+	std::vector<t_mat>	// symops
+	>>
 get_sgs(bool bAddNr=true, bool bAddHall=true)
 {
 	std::vector<std::tuple<int, std::string, std::vector<t_mat>>> sgs;
@@ -302,4 +306,54 @@ get_sgs(bool bAddNr=true, bool bAddHall=true)
 }
 
 
+
+/**
+ * finds all space groups which transform the initial positions into the final ones
+ */
+template<class t_vec, class t_mat, class t_real = typename t_mat::value_type>
+std::vector<std::tuple<
+	int,				// sg number 
+	std::string, 		// description
+	std::vector<t_mat>	// symops
+	>>
+find_matching_sgs(
+	const std::vector<t_vec>& posInit, const std::vector<t_vec>& _posFinal, 
+	t_real eps=1e-6)
+{
+	std::vector<t_vec> posFinal = m::keep_atoms_in_uc<t_vec, t_real>(_posFinal);
+
+
+	std::vector<std::tuple<int, std::string, std::vector<t_mat>>> matchingSGs;
+	auto sgs = get_sgs<t_mat, t_real>();
+
+	// iterate spacegroups
+	for(const auto& [sgNum, sgName, sgOps] : sgs)
+	{
+		// generate symmetry-equivalent positions
+		std::vector<t_vec> generatedpos;
+
+		for(const t_vec& pos : posInit)
+		{
+			std::vector<t_vec> newpos = m::apply_ops_hom<t_vec, t_mat, t_real>(pos, sgOps, eps);
+			generatedpos.insert(generatedpos.end(), newpos.begin(), newpos.end());
+		}
+
+		//for(const auto& thepos : generatedpos)
+		//	std::cout << thepos << std::endl;
+
+		// filter multiple occupancies in generatedpos
+		generatedpos = m::remove_duplicates<t_vec>(generatedpos, eps);
+
+
+		// no match
+		if(!m::equals_all<t_vec>(generatedpos, posFinal, eps, 3))
+			continue;
+
+		matchingSGs.emplace_back(std::make_tuple(sgNum, sgName, sgOps));
+	}
+
+	return matchingSGs;
+}
+
+
 #endif