continued with bz tool

tools/bz/bz.cpp

@@ -415,6 +415,10 @@ void BZDlg::CalcBZ()
 	std::ostringstream ostr;
 	ostr.precision(g_prec);
 
+	ostr << "# centring symmetry operations" << std::endl;
+	for(const t_mat& op : ops_centr)
+		ostr << op << std::endl;
+
 	for(t_real h=-maxBZ; h<=maxBZ; ++h)
 	{
 		for(t_real k=-maxBZ; k<=maxBZ; ++k)
@@ -422,9 +426,12 @@ void BZDlg::CalcBZ()
 			for(t_real l=-maxBZ; l<=maxBZ; ++l)
 			{
 				auto Q = tl2::create<t_vec>({ h, k, l });
+				if(!is_reflection_allowed<t_mat, t_vec, t_real>(
+					Q, ops_centr, g_eps).first)
+					continue;
+
 				auto Q_invA = m_crystB * Q;
 				auto Qabs_invA = tl2::norm(Q_invA);
-
 			}
 		}
 	}

tools/bz/bz.h

@@ -71,12 +71,6 @@ enum : int
 };
 
 
-struct SymOp
-{
-	t_mat op;
-};
-
-
 class BZDlg : public QDialog
 {
 public:
@@ -133,7 +127,7 @@ protected:
 	void ImportCIF();
 	void GetSymOpsFromSG();
 
-	std::vector<SymOp> GetSymOps(bool only_centring = false) const;
+	std::vector<t_mat> GetSymOps(bool only_centring = false) const;
 	void CalcB(bool bFullRecalc=true);
 	void CalcBZ();
 

tools/bz/bz_ops.cpp

@@ -73,7 +73,7 @@ t_mat BZDlg::StrToOp(const std::string& str)
 {
 	t_mat op = tl2::unit<t_mat>(4);
 
-	std::istringstream istr;
+	std::istringstream istr(str);
 	for(std::size_t row=0; row<op.size1(); ++row)
 		for(std::size_t col=0; col<op.size2(); ++col)
 			istr >> op(row, col);
@@ -379,30 +379,28 @@ void BZDlg::GetSymOpsFromSG()
 /**
  * reads symmetry operations from table
  */
-std::vector<SymOp> BZDlg::GetSymOps(bool only_centring) const
+std::vector<t_mat> BZDlg::GetSymOps(bool only_centring) const
 {
-	std::vector<SymOp> vec;
+	std::vector<t_mat> ops;
 
 	for(int row=0; row<m_symops->rowCount(); ++row)
 	{
-		auto *op = m_symops->item(row, COL_OP);
-
-		if(!op)
+		auto *op_item = m_symops->item(row, COL_OP);
+		if(!op_item)
 		{
 			std::cerr << "Invalid entry in row " << row << "." << std::endl;
 			continue;
 		}
 
-		SymOp symop;
-		symop.op = StrToOp(op->text().toStdString());
+		t_mat op = StrToOp(op_item->text().toStdString());
 
 		bool add_op = true;
 		if(only_centring)
-			add_op = tl2::hom_is_centering<t_mat>(symop.op, g_eps);
+			add_op = tl2::hom_is_centering<t_mat>(op, g_eps);
 
 		if(add_op)
-			vec.emplace_back(std::move(symop));
+			ops.emplace_back(std::move(op));
 	}
 
-	return vec;
+	return ops;
 }

tools/structfact/loadcif.h

@@ -37,11 +37,10 @@
 #include <gemmi/cif.hpp>
 #include <gemmi/symmetry.hpp>
 
-#include "tlibs2/libs/math20.h"
+#include "tlibs2/libs/maths.h"
 using namespace tl2_ops;
 
 
-
 template<class t_real=double>
 struct Lattice
 {
@@ -50,7 +49,6 @@ struct Lattice
 };
 
 
-
 /**
  * removes quotation marks around a string
  */
@@ -67,7 +65,6 @@ void remove_quotes(t_str& str)
 }
 
 
-
 /**
  * gets the atom positions from the CIF
  */
@@ -118,7 +115,6 @@ get_cif_atoms(gemmi::cif::Block& block)
 }
 
 
-
 /**
  * gets the symmetry operations from the CIF
  */
@@ -167,7 +163,6 @@ std::vector<t_mat> get_cif_ops(gemmi::cif::Block& block)
 }
 
 
-
 /**
  * gets the symmetry operations from the CIF's space group
  * (use tl2::equals_all to check if space group operations are the same)
@@ -204,7 +199,6 @@ std::vector<t_mat> get_cif_sg_ops(gemmi::cif::Block& block)
 }
 
 
-
 /**
  * loads the lattice parameters and the atom positions from a CIF
  */
@@ -280,7 +274,6 @@ load_cif(const std::string& filename, t_real eps=1e-6)
 }
 
 
-
 /**
  * gets space group description strings and symmetry operations
  */
@@ -327,7 +320,6 @@ get_sgs(bool bAddNr=true, bool bAddHall=true)
 }
 
 
-
 /**
  * finds all space groups which transform the initial positions into the final ones
  */
@@ -377,4 +369,43 @@ find_matching_sgs(
 }
 
 
+/**
+ * checks for allowed Bragg reflections
+ *
+ * algorithm based on Clipper's HKL_class
+ *   constructor in clipper/core/coords.cpp by K. Cowtan, 2013
+ * @see http://www.ysbl.york.ac.uk/~cowtan/clipper/
+ *
+ * symmetry operation S on position r: R*r + t
+ * F = sum<S>( exp(2*pi*i * (R*r + t)*G) )
+ *   = sum<S>( exp(2*pi*i * ((R*r)*G + t*G)) )
+ *   = sum<S>( exp(2*pi*i * (r*(G*R) + t*G)) )
+ *   = sum<S>( exp(2*pi*i * (r*(G*R)))  *  exp(2*pi*i * (G*t)) )
+ */
+template<class t_mat, class t_vec, class t_real = typename t_vec::value_type,
+	template<class...> class t_cont = std::vector>
+std::pair<bool, std::size_t>
+is_reflection_allowed(const t_vec& Q, const t_cont<t_mat>& symops, t_real eps)
+requires tl2::is_mat<t_mat> && tl2::is_vec<t_vec>
+{
+	for(std::size_t opidx=0; opidx<symops.size(); ++opidx)
+	{
+		const t_mat& mat = symops[opidx];
+		t_mat rot = tl2::submat<t_mat>(mat, 0,0, 3,3); // rotation part of the symop
+		rot = tl2::trans(rot);                         // recip -> transpose
+
+		if(tl2::equals<t_vec>(Q, rot*Q, eps))          // does Q transform into itself
+		{
+			t_vec trans = tl2::create<t_vec>({ mat(0,3), mat(1,3), mat(2,3) });
+
+			// does Q translate to multiples of the lattice vector?
+			if(!tl2::is_integer<t_real>(tl2::inner<t_vec>(trans, Q), eps))
+				return std::make_pair(false, opidx);
+		}
+	}
+
+	return std::make_pair(true, symops.size());
+}
+
+
 #endif