added support for cif loading

tools/structfact/CMakeLists.txt

@@ -20,7 +20,8 @@ set(CMAKE_CXX_STANDARD 17)
 add_definitions(-std=c++17 -fconcepts)
 add_definitions(${Boost_CXX_FLAGS})
 
-include_directories("${PROJECT_SOURCE_DIR}" "${Boost_INCLUDE_DIRS}/.." "../..")
+include_directories("${PROJECT_SOURCE_DIR}" "${Boost_INCLUDE_DIRS}/.." "../.."
+	"../../ext/gemmi/include" "../../ext/gemmi/third_party")
 
 if(BUILD_LIB)
 	set(CMAKE_POSITION_INDEPENDENT_CODE TRUE)

tools/structfact/loadcif.h

+/**
+ * get atom positions from cif
+ * @author Tobias Weber <tweber@ill.fr>
+ * @date Jan-2019
+ * @license GPLv3, see 'LICENSE' file
+ */
+
+#ifndef __LOAD_CIF_H__
+#define __LOAD_CIF_H__
+
+#include <vector>
+#include <string>
+#include <fstream>
+#include <sstream>
+
+#include <boost/algorithm/string.hpp>
+
+#include <gemmi/cif.hpp>
+#include <gemmi/symmetry.hpp>
+
+#include "libs/_cxx20/math_algos.h"
+using namespace m_ops;
+
+
+template<class t_real=double>
+struct Lattice
+{
+	t_real a, b, c;
+	t_real alpha, beta, gamma;
+};
+
+
+/**
+ * gets the symmetry operations from the CIF
+ */
+template<class t_vec, class t_mat, class t_real = typename t_vec::value_type>
+std::vector<t_mat> get_cif_ops(gemmi::cif::Block& block)
+{
+	std::vector<t_mat> ops;
+
+	auto colOps = block.find_values("_symmetry_equiv_pos_as_xyz");
+	for(std::size_t row=0; row<colOps.length(); ++row)
+	{
+		auto op = gemmi::parse_triplet(colOps[row])/*.wrap()*/;
+		auto M = op.float_seitz();
+
+		t_mat mat = m::create<t_mat>({
+			std::get<0>(std::get<0>(M)), std::get<1>(std::get<0>(M)), std::get<2>(std::get<0>(M)), std::get<3>(std::get<0>(M)),
+			std::get<0>(std::get<1>(M)), std::get<1>(std::get<1>(M)), std::get<2>(std::get<1>(M)), std::get<3>(std::get<1>(M)),
+			std::get<0>(std::get<2>(M)), std::get<1>(std::get<2>(M)), std::get<2>(std::get<2>(M)), std::get<3>(std::get<2>(M)),
+			std::get<0>(std::get<3>(M)), std::get<1>(std::get<3>(M)), std::get<2>(std::get<3>(M)), std::get<3>(std::get<3>(M)) });
+
+		ops.emplace_back(std::move(mat));
+	}
+
+	return ops;
+}
+
+
+/**
+ * gets the symmetry operations from the CIF's space group
+ */
+template<class t_vec, class t_mat, class t_real = typename t_vec::value_type>
+std::vector<t_mat> get_cif_sg_ops(gemmi::cif::Block& block)
+{
+	std::vector<t_mat> ops;
+
+	if(auto val = block.find_values("_symmetry_space_group_name_H-M"); val.length())
+	{
+		std::string sgname = boost::trim_copy(val[0]);
+
+		// remove quotation marks
+		if(sgname[0] == '\'' || sgname[0] == '\"')
+			sgname.erase(sgname.begin());
+		if(sgname[sgname.size()-1] == '\'' || sgname[sgname.size()-1] == '\"')
+			sgname.erase(sgname.begin()+sgname.size()-1);
+
+		if(auto sg = gemmi::find_spacegroup_by_name(sgname))
+		{
+			auto symops = sg->operations().all_ops_sorted();
+			for(const auto &op : symops)
+			{
+				auto M = op.float_seitz();
+
+				t_mat mat = m::create<t_mat>({
+					std::get<0>(std::get<0>(M)), std::get<1>(std::get<0>(M)), std::get<2>(std::get<0>(M)), std::get<3>(std::get<0>(M)),
+					std::get<0>(std::get<1>(M)), std::get<1>(std::get<1>(M)), std::get<2>(std::get<1>(M)), std::get<3>(std::get<1>(M)),
+					std::get<0>(std::get<2>(M)), std::get<1>(std::get<2>(M)), std::get<2>(std::get<2>(M)), std::get<3>(std::get<2>(M)),
+					std::get<0>(std::get<3>(M)), std::get<1>(std::get<3>(M)), std::get<2>(std::get<3>(M)), std::get<3>(std::get<3>(M)) });
+
+				ops.emplace_back(std::move(mat));
+			}
+		}
+	}
+
+	return ops;
+}
+
+
+/**
+ * loads the lattice parameters and the atom positions from a CIF
+ */
+template<class t_vec, class t_mat, class t_real = typename t_vec::value_type>
+std::tuple<const char*, std::vector<t_vec>, std::vector<std::vector<t_vec>>, std::vector<std::string>, Lattice<t_real>> 
+load_cif(const std::string& filename)
+{
+	auto ifstr = std::ifstream(filename);
+	if(!ifstr)
+		return std::make_tuple("Cannot open CIF.", std::vector<t_vec>{}, std::vector<std::vector<t_vec>>{}, std::vector<std::string>{}, Lattice{});
+
+	// load CIF
+	auto cif = gemmi::cif::read_istream(ifstr, 4096, filename.c_str());
+
+	if(!cif.blocks.size())
+		return std::make_tuple("No blocks in CIF.", std::vector<t_vec>{}, std::vector<std::vector<t_vec>>{}, std::vector<std::string>{}, Lattice{});
+
+	// get the block
+	/*const*/ auto& block = cif.sole_block();
+
+
+	// lattice
+	t_real a{}, b{}, c{}, alpha{}, beta{}, gamma{};
+	if(auto val = block.find_values("_cell_length_a"); val.length()) std::istringstream{val[0]} >> a;
+	if(auto val = block.find_values("_cell_length_b"); val.length()) std::istringstream{val[0]} >> b;
+	if(auto val = block.find_values("_cell_length_c"); val.length()) std::istringstream{val[0]} >> c;
+	if(auto val = block.find_values("_cell_angle_alpha"); val.length()) std::istringstream{val[0]} >> alpha;
+	if(auto val = block.find_values("_cell_angle_beta"); val.length()) std::istringstream{val[0]} >> beta;
+	if(auto val = block.find_values("_cell_angle_gamma"); val.length()) std::istringstream{val[0]} >> gamma;
+
+	Lattice<t_real> latt{.a=a, .b=b, .c=c, .alpha=alpha, .beta=beta, .gamma=gamma};
+
+
+	// fractional atom positions
+	std::vector<t_vec> atoms;
+	auto tabAtoms = block.find("_atom_site", {"_type_symbol", "_fract_x", "_fract_y", "_fract_z"});
+
+	std::vector<std::string> atomnames;
+	for(std::size_t row=0; row<tabAtoms.length(); ++row)
+	{
+		atomnames.push_back(tabAtoms[row][0]);
+
+		t_real x{}, y{}, z{};
+		std::istringstream{tabAtoms[row][1]} >> x;
+		std::istringstream{tabAtoms[row][2]} >> y;
+		std::istringstream{tabAtoms[row][3]} >> z;
+		atoms.emplace_back(t_vec{{x, y, z}});
+	}
+
+
+	// generate all atoms using symmetry ops
+	std::vector<std::vector<t_vec>> generatedatoms;
+	auto ops = get_cif_ops<t_vec, t_mat, t_real>(block);
+	if(!ops.size()) // if ops are not directly given, use standard ops from space group
+		ops = get_cif_sg_ops<t_vec, t_mat, t_real>(block);
+
+	for(t_vec atom : atoms)
+	{
+		// make homogeneuous 4-vector
+		if(atom.size() == 3) atom.push_back(1);
+
+		std::vector<t_vec> newatoms;
+
+		for(const auto& op : ops)
+		{
+			auto newatom = op*atom;
+			newatom.resize(3);
+
+			for(int i=0; i<3; ++i)
+			{
+				newatom[i] = std::fmod(newatom[i], 1);
+				while(newatom[i] < 0) newatom[i] += 1;
+				while(newatom[i] >= 1) newatom[i] -= 1;
+
+				newatom[i] -= 0.5;
+			}
+
+
+			// position already occupied?
+			if(std::find_if(newatoms.begin(), newatoms.end(), [&newatom](const t_vec& vec)->bool
+			{
+				return m::equals<t_vec>(vec, newatom);
+			}) == newatoms.end())
+			{
+				newatoms.emplace_back(std::move(newatom));
+			}
+		}
+
+		generatedatoms.emplace_back(std::move(newatoms));
+	}
+
+
+	return std::make_tuple(nullptr, atoms, generatedatoms, atomnames, latt);
+}
+
+#endif

tools/structfact/structfact.cpp

@@ -19,15 +19,19 @@
 
 #include <iostream>
 #include <fstream>
+#include <random>
+#include <chrono>
 
 #include <boost/version.hpp>
 #include <boost/config.hpp>
 #include <boost/property_tree/ptree.hpp>
 #include <boost/property_tree/xml_parser.hpp>
+#include <boost/algorithm/string.hpp>
 #include <boost/algorithm/string/replace.hpp>
 namespace algo = boost::algorithm;
 namespace pt = boost::property_tree;
 
+#include "loadcif.h"
 #include "libs/algos.h"
 #include "libs/helper.h"
 #include "libs/_cxx20/math_algos.h"
@@ -110,6 +114,7 @@ StructFactDlg::StructFactDlg(QWidget* pParent) : QDialog{pParent},
 		QToolButton *pTabBtnDown = new QToolButton(m_nucleipanel);
 		QToolButton *pTabBtnLoad = new QToolButton(m_nucleipanel);
 		QToolButton *pTabBtnSave = new QToolButton(m_nucleipanel);
+		QToolButton *pTabBtnImportCIF = new QToolButton(m_nucleipanel);
 		QToolButton *pTabBtn3DView = new QToolButton(m_nucleipanel);
 
 		m_nuclei->setSizePolicy(QSizePolicy{QSizePolicy::Expanding, QSizePolicy::Expanding});
@@ -119,6 +124,7 @@ StructFactDlg::StructFactDlg(QWidget* pParent) : QDialog{pParent},
 		pTabBtnDown->setSizePolicy(QSizePolicy{QSizePolicy::Expanding, QSizePolicy::Fixed});
 		pTabBtnLoad->setSizePolicy(QSizePolicy{QSizePolicy::Expanding, QSizePolicy::Fixed});
 		pTabBtnSave->setSizePolicy(QSizePolicy{QSizePolicy::Expanding, QSizePolicy::Fixed});
+		pTabBtnImportCIF->setSizePolicy(QSizePolicy{QSizePolicy::Expanding, QSizePolicy::Fixed});
 		pTabBtn3DView->setSizePolicy(QSizePolicy{QSizePolicy::Expanding, QSizePolicy::Fixed});
 
 		pTabBtnAdd->setText("Add Nucleus");
@@ -127,6 +133,7 @@ StructFactDlg::StructFactDlg(QWidget* pParent) : QDialog{pParent},
 		pTabBtnDown->setText("Move Nuclei Down");
 		pTabBtnLoad->setText("Load...");
 		pTabBtnSave->setText("Save...");
+		pTabBtnImportCIF->setText("Import CIF...");
 		pTabBtn3DView->setText("3D View...");
 
 
@@ -150,6 +157,7 @@ StructFactDlg::StructFactDlg(QWidget* pParent) : QDialog{pParent},
 		pTabGrid->addWidget(pTabBtnDown, y,3,1,1);
 		pTabGrid->addWidget(pTabBtnLoad, ++y,0,1,1);
 		pTabGrid->addWidget(pTabBtnSave, y,1,1,1);
+		pTabGrid->addWidget(pTabBtnImportCIF, y,2,1,1);
 		pTabGrid->addWidget(pTabBtn3DView, y,3,1,1);
 
 		auto sep1 = new QFrame(m_nucleipanel); sep1->setFrameStyle(QFrame::HLine);
@@ -190,6 +198,7 @@ StructFactDlg::StructFactDlg(QWidget* pParent) : QDialog{pParent},
 		connect(pTabBtnDown, &QToolButton::clicked, this, &StructFactDlg::MoveTabItemDown);
 		connect(pTabBtnLoad, &QToolButton::clicked, this, &StructFactDlg::Load);
 		connect(pTabBtnSave, &QToolButton::clicked, this, &StructFactDlg::Save);
+		connect(pTabBtnImportCIF, &QToolButton::clicked, this, &StructFactDlg::ImportCIF);
 		connect(pTabBtn3DView, &QToolButton::clicked, this, [this]()
 		{
 			// plot widget
@@ -801,6 +810,73 @@ void StructFactDlg::Save()
 	std::ofstream ofstr{filename.toStdString()};
 	pt::write_xml(ofstr, node, pt::xml_writer_make_settings('\t', 1, std::string{"utf-8"}));
 }
+
+
+void StructFactDlg::ImportCIF()
+{
+	QString dirLast = m_sett->value("dir", "").toString();
+	QString filename = QFileDialog::getOpenFileName(this, "Import CIF", dirLast, "CIF Files (*.cif *.CIF)");
+	if(filename=="" || !QFile::exists(filename))
+		return;
+	m_sett->setValue("dir", QFileInfo(filename).path());
+
+	auto [errstr, atoms, generatedatoms, atomnames, lattice] = load_cif<t_vec, t_mat>(filename.toStdString());
+	if(errstr)
+	{
+		QMessageBox::critical(this, "Structure Factors", errstr);
+		return;
+	}
+
+
+	// clear old nuclei
+	DelTabItem(true);
+
+	// lattice
+	{
+		std::ostringstream ostr; ostr << lattice.a;
+		m_editA->setText(ostr.str().c_str());
+	}
+	{
+		std::ostringstream ostr; ostr << lattice.b;
+		m_editB->setText(ostr.str().c_str());
+	}
+	{
+		std::ostringstream ostr; ostr << lattice.c;
+		m_editC->setText(ostr.str().c_str());
+	}
+	{
+		std::ostringstream ostr; ostr << lattice.alpha;
+		m_editAlpha->setText(ostr.str().c_str());
+	}
+	{
+		std::ostringstream ostr; ostr << lattice.beta;
+		m_editBeta->setText(ostr.str().c_str());
+	}
+	{
+		std::ostringstream ostr; ostr << lattice.gamma;
+		m_editGamma->setText(ostr.str().c_str());
+	}
+
+
+	// atoms
+	std::mt19937 gen{tl2::epoch<unsigned int>()};
+	for(std::size_t atomnum=0; atomnum<atoms.size(); ++atomnum)
+	{
+		// random colour
+		std::ostringstream ostrcol;
+		std::uniform_int_distribution<int> dist{0, 255};
+		ostrcol << "#" << std::hex << std::setw(2) << std::setfill('0') << dist(gen) 
+			<< std::setw(2) << std::setfill('0') << dist(gen) 
+			<< std::setw(2) << std::setfill('0') << dist(gen);
+
+		for(std::size_t symnr=0; symnr<generatedatoms[atomnum].size(); ++symnr)
+		{
+			AddTabItem(-1, atomnames[atomnum], 0, 0, 
+				generatedatoms[atomnum][symnr][0],  generatedatoms[atomnum][symnr][1], generatedatoms[atomnum][symnr][2], 
+				1, ostrcol.str());
+		}
+	}
+}
 // ----------------------------------------------------------------------------
 
 

tools/structfact/structfact.h

@@ -118,6 +118,7 @@ protected:
 
 	void Load();
 	void Save();
+	void ImportCIF();
 
 	std::vector<NuclPos> GetNuclei() const;
 	void Calc();