Skip to content
GitLab
Verified Commit 290cd478 authored by Tobias WEBER's avatar Tobias WEBER
Browse files

bz tool: added cut plane

parent c599e326
Hide whitespace changes
Inline Side-by-side
tools/bz/CMakeLists.txt
View file @ 290cd478
......@@ -69,7 +69,7 @@ include_directories(
add_executable(takin_bz
bz.cpp bz_file.cpp
bz.cpp bz_file.cpp bz_calc.cpp
bz_plot.cpp bz_main.cpp
bz_ops.cpp bz.h globals.h
plot_cut.cpp plot_cut.h
......
tools/bz/bz.cpp
View file @ 290cd478
......@@ -43,16 +43,8 @@
namespace algo = boost::algorithm;
#include "../structfact/loadcif.h"
#include "tlibs2/libs/phys.h"
#include "tlibs2/libs/algos.h"
#include "tlibs2/libs/qt/helper.h"
#include "pathslib/libs/voronoi.h"
using namespace tl2_ops;
// ----------------------------------------------------------------------------
BZDlg::BZDlg(QWidget* pParent) : QDialog{pParent},
m_sett{new QSettings{"takin", "bz"}}
{
......@@ -79,7 +71,7 @@ BZDlg::BZDlg(QWidget* pParent) : QDialog{pParent},
m_symops->verticalHeader()->setVisible(false);
m_symops->setColumnCount(NUM_COLS);
m_symops->setHorizontalHeaderItem(COL_OP, new QTableWidgetItem{"Symmetry Operation"});
m_symops->setHorizontalHeaderItem(COL_OP, new QTableWidgetItem{"Symmetry Operations"});
m_symops->setColumnWidth(COL_OP, 500);
QToolButton *btnAdd = new QToolButton(symopspanel);
......@@ -198,8 +190,12 @@ BZDlg::BZDlg(QWidget* pParent) : QDialog{pParent},
(*cut)->setMinimum(-99);
(*cut)->setMaximum(99);
(*cut)->setDecimals(2);
(*cut)->setSingleStep(1);
(*cut)->setSingleStep(0.1);
(*cut)->setValue(0);
connect(*cut, static_cast<void (QDoubleSpinBox::*)(double)>
(&QDoubleSpinBox::valueChanged),
this, &BZDlg::CalcBZCut);
}
m_cutZ->setValue(1);
......@@ -211,10 +207,7 @@ BZDlg::BZDlg(QWidget* pParent) : QDialog{pParent},
pGrid->addWidget(new QLabel("Plane Offset:"), 2,0, 1,1);
pGrid->addWidget(m_cutD, 2,1, 1,1);
// signals
// TODO
tabs->addTab(cutspanel, "Cuts");
tabs->addTab(cutspanel, "Cut");
}
......@@ -392,170 +385,3 @@ void BZDlg::closeEvent(QCloseEvent *)
m_sett->setValue("geo_3dview", m_dlgPlot->saveGeometry());
}
}
// ----------------------------------------------------------------------------
/**
* calculate crystal B matrix
*/
void BZDlg::CalcB(bool bFullRecalc)
{
if(m_ignoreCalc)
return;
t_real a,b,c, alpha,beta,gamma;
std::istringstream{m_editA->text().toStdString()} >> a;
std::istringstream{m_editB->text().toStdString()} >> b;
std::istringstream{m_editC->text().toStdString()} >> c;
std::istringstream{m_editAlpha->text().toStdString()} >> alpha;
std::istringstream{m_editBeta->text().toStdString()} >> beta;
std::istringstream{m_editGamma->text().toStdString()} >> gamma;
if(tl2::equals<t_real>(a, 0., g_eps) || a <= 0. ||
tl2::equals<t_real>(b, 0., g_eps) || b <= 0. ||
tl2::equals<t_real>(c, 0., g_eps) || c <= 0. ||
tl2::equals<t_real>(alpha, 0., g_eps) || alpha <= 0. ||
tl2::equals<t_real>(beta, 0., g_eps) || beta <= 0. ||
tl2::equals<t_real>(gamma, 0., g_eps) || gamma <= 0.)
{
QMessageBox::critical(this, "Brillouin Zones",
"Error: Invalid lattice.");
return;
}
t_mat crystB = tl2::B_matrix<t_mat>(a, b, c,
alpha/180.*tl2::pi<t_real>,
beta/180.*tl2::pi<t_real>,
gamma/180.*tl2::pi<t_real>);
bool ok = true;
t_mat crystA = tl2::unit<t_mat>(3);
std::tie(crystA, ok) = tl2::inv(crystB);
if(!ok)
{
QMessageBox::critical(this, "Brillouin Zones",
"Error: Cannot invert B matrix.");
return;
}
m_crystA = crystA * t_real(2)*tl2::pi<t_real>;
m_crystB = crystB;
if(m_plot)
{
t_mat_gl matA{m_crystA};
m_plot->GetRenderer()->SetBTrafo(m_crystB, &matA);
}
if(bFullRecalc)
CalcBZ();
}
/**
* calculate brillouin zone
*/
void BZDlg::CalcBZ()
{
if(m_ignoreCalc)
return;
const auto maxBZ = m_maxBZ->value();
const auto ops_centr = GetSymOps(true);
std::ostringstream ostr;
ostr.precision(g_prec);
#ifdef DEBUG
ostr << "# centring symmetry operations" << std::endl;
for(const t_mat& op : ops_centr)
ostr << op << std::endl;
#endif
std::vector<t_vec> Qs_invA;
Qs_invA.reserve((2*maxBZ+1)*(2*maxBZ+1)*(2*maxBZ+1));
std::size_t idx000 = 0;
for(t_real h=-maxBZ; h<=maxBZ; ++h)
{
for(t_real k=-maxBZ; k<=maxBZ; ++k)
{
for(t_real l=-maxBZ; l<=maxBZ; ++l)
{
t_vec 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;
if(tl2::equals_0(Q, g_eps))
idx000 = Qs_invA.size();
t_vec Q_invA = m_crystB * Q;
t_real Qabs_invA = tl2::norm(Q_invA);
Qs_invA.emplace_back(std::move(Q_invA));
}
}
}
// calculate voronoi diagram
auto [voronoi, triags, neighbours] =
geo::calc_delaunay(3, Qs_invA, false, false, idx000);
ClearPlot();
#ifdef DEBUG
std::ofstream ofstrSites("sites.dat");
std::cout << "cat sites.dat | qvoronoi s p Fv QV" << idx000 << std::endl;
ofstrSites << "3 " << Qs_invA.size() << std::endl;
for(const t_vec& Q : Qs_invA)
{
//PlotAddBraggPeak(Q);
ofstrSites << Q[0] << " " << Q[1] << " " << Q[2] << std::endl;
}
#endif
// add gamma point
PlotAddBraggPeak(Qs_invA[idx000]);
// add voronoi vertices forming the vertices of the BZ
ostr << "\n# Brillouin zone vertices" << std::endl;
for(std::size_t idx=0; idx<voronoi.size(); ++idx)
{
t_vec& voro = voronoi[idx];
tl2::set_eps_0(voro, g_eps);
PlotAddVoronoiVertex(voro);
ostr << "vertex " << idx << ": " << voro << std::endl;
//for(std::size_t nidx : neighbours[idx])
// ostr << "\tneighbour index: " << nidx << std::endl;
}
// calculate the faces of the BZ
auto [bz_verts, bz_triags, bz_neighbours] =
geo::calc_delaunay(3, voronoi, true, false);
std::vector<t_vec> bz_all_triags;
ostr << "\n# Brillouin zone polygons" << std::endl;
for(std::size_t idx_triag=0; idx_triag<bz_triags.size(); ++idx_triag)
{
const auto& triag = bz_triags[idx_triag];
ostr << "polygon " << idx_triag << ": " << std::endl;
for(std::size_t idx_vert=0; idx_vert<triag.size(); ++idx_vert)
{
bz_all_triags.push_back(triag[idx_vert]);
ostr << "\tvertex: " << triag[idx_vert] << std::endl;
}
}
PlotAddTriangles(bz_all_triags);
// brillouin zone description
m_bz->setPlainText(ostr.str().c_str());
}
// ----------------------------------------------------------------------------
tools/bz/bz.h
View file @ 290cd478
......@@ -74,8 +74,12 @@ protected:
QDialog *m_dlgPlot = nullptr;
std::shared_ptr<tl2::GlPlot> m_plot;
std::size_t m_sphere = 0;
std::size_t m_plane = 0;
QLabel *m_labelGlInfos[4] = { nullptr, nullptr, nullptr, nullptr };
QLabel *m_status3D = nullptr;
QCheckBox *m_plot_coordcross = nullptr;
QCheckBox *m_plot_labels = nullptr;
QCheckBox *m_plot_plane = nullptr;
// symops panel
QLineEdit *m_editA = nullptr;
......@@ -128,14 +132,20 @@ protected:
std::vector<t_mat> GetSymOps(bool only_centring = false) const;
void CalcB(bool bFullRecalc=true);
void CalcBZ();
void CalcBZCut();
void ClearPlot();
void PlotAddVoronoiVertex(const t_vec& pos);
void PlotAddBraggPeak(const t_vec& pos);
void PlotAddTriangles(const std::vector<t_vec>& vecs);
void PlotSetPlane(const t_vec& norm, t_real d);
void Set3DStatusMsg(const std::string& msg);
void ShowBZPlot();
void PlotShowCoordCross(bool show);
void PlotShowLabels(bool show);
void PlotShowPlane(bool show);
void PlotMouseDown(bool left, bool mid, bool right);
void PlotMouseUp(bool left, bool mid, bool right);
void PickerIntersection(const t_vec3_gl* pos,
......
tools/bz/bz_calc.cpp 0 → 100644
View file @ 290cd478
/**
* brillouin zone tool
* @author Tobias Weber <tweber@ill.fr>
* @date May-2022
* @license GPLv3, see 'LICENSE' file
*
* ----------------------------------------------------------------------------
* mag-core (part of the Takin software suite)
* Copyright (C) 2018-2022 Tobias WEBER (Institut Laue-Langevin (ILL),
* Grenoble, France).
* "misc" project
* Copyright (C) 2017-2021 Tobias WEBER (privately developed).
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, version 3 of the License.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
* ----------------------------------------------------------------------------
*/
#include "bz.h"
#include <QtWidgets/QMessageBox>
#include <iostream>
#include <sstream>
#include "../structfact/loadcif.h"
#include "tlibs2/libs/phys.h"
#include "tlibs2/libs/algos.h"
#include "tlibs2/libs/qt/helper.h"
#include "pathslib/libs/voronoi.h"
using namespace tl2_ops;
/**
* calculate crystal B matrix
*/
void BZDlg::CalcB(bool bFullRecalc)
{
if(m_ignoreCalc)
return;
t_real a,b,c, alpha,beta,gamma;
std::istringstream{m_editA->text().toStdString()} >> a;
std::istringstream{m_editB->text().toStdString()} >> b;
std::istringstream{m_editC->text().toStdString()} >> c;
std::istringstream{m_editAlpha->text().toStdString()} >> alpha;
std::istringstream{m_editBeta->text().toStdString()} >> beta;
std::istringstream{m_editGamma->text().toStdString()} >> gamma;
if(tl2::equals<t_real>(a, 0., g_eps) || a <= 0. ||
tl2::equals<t_real>(b, 0., g_eps) || b <= 0. ||
tl2::equals<t_real>(c, 0., g_eps) || c <= 0. ||
tl2::equals<t_real>(alpha, 0., g_eps) || alpha <= 0. ||
tl2::equals<t_real>(beta, 0., g_eps) || beta <= 0. ||
tl2::equals<t_real>(gamma, 0., g_eps) || gamma <= 0.)
{
QMessageBox::critical(this, "Brillouin Zones",
"Error: Invalid lattice.");
return;
}
t_mat crystB = tl2::B_matrix<t_mat>(a, b, c,
alpha/180.*tl2::pi<t_real>,
beta/180.*tl2::pi<t_real>,
gamma/180.*tl2::pi<t_real>);
bool ok = true;
t_mat crystA = tl2::unit<t_mat>(3);
std::tie(crystA, ok) = tl2::inv(crystB);
if(!ok)
{
QMessageBox::critical(this, "Brillouin Zones",
"Error: Cannot invert B matrix.");
return;
}
m_crystA = crystA * t_real(2)*tl2::pi<t_real>;
m_crystB = crystB;
if(m_plot)
{
t_mat_gl matA{m_crystA};
m_plot->GetRenderer()->SetBTrafo(m_crystB, &matA);
}
if(bFullRecalc)
CalcBZ();
}
/**
* calculate brillouin zone
*/
void BZDlg::CalcBZ()
{
if(m_ignoreCalc)
return;
const auto maxBZ = m_maxBZ->value();
const auto ops_centr = GetSymOps(true);
std::ostringstream ostr;
ostr.precision(g_prec);
#ifdef DEBUG
ostr << "# centring symmetry operations" << std::endl;
for(const t_mat& op : ops_centr)
ostr << op << std::endl;
#endif
std::vector<t_vec> Qs_invA;
Qs_invA.reserve((2*maxBZ+1)*(2*maxBZ+1)*(2*maxBZ+1));
std::size_t idx000 = 0;
for(t_real h=-maxBZ; h<=maxBZ; ++h)
{
for(t_real k=-maxBZ; k<=maxBZ; ++k)
{
for(t_real l=-maxBZ; l<=maxBZ; ++l)
{
t_vec 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;
if(tl2::equals_0(Q, g_eps))
idx000 = Qs_invA.size();
t_vec Q_invA = m_crystB * Q;
t_real Qabs_invA = tl2::norm(Q_invA);
Qs_invA.emplace_back(std::move(Q_invA));
}
}
}
// calculate voronoi diagram
auto [voronoi, _triags, _neighbours] =
geo::calc_delaunay(3, Qs_invA, false, false, idx000);
voronoi = tl2::remove_duplicates(voronoi, g_eps);
ClearPlot();
#ifdef DEBUG
std::ofstream ofstrSites("sites.dat");
std::cout << "cat sites.dat | qvoronoi s p Fv QV" << idx000 << std::endl;
ofstrSites << "3 " << Qs_invA.size() << std::endl;
for(const t_vec& Q : Qs_invA)
{
//PlotAddBraggPeak(Q);
ofstrSites << Q[0] << " " << Q[1] << " " << Q[2] << std::endl;
}
#endif
// add gamma point
PlotAddBraggPeak(Qs_invA[idx000]);
// add voronoi vertices forming the vertices of the BZ
ostr << "# Brillouin zone vertices" << std::endl;
for(std::size_t idx=0; idx<voronoi.size(); ++idx)
{
t_vec& voro = voronoi[idx];
tl2::set_eps_0(voro, g_eps);
PlotAddVoronoiVertex(voro);
ostr << "vertex " << idx << ": " << voro << std::endl;
}
// calculate the faces of the BZ
auto [bz_verts, bz_triags, bz_neighbours] =
geo::calc_delaunay(3, voronoi, true, false);
std::vector<t_vec> bz_all_triags;
ostr << "\n# Brillouin zone polygons" << std::endl;
for(std::size_t idx_triag=0; idx_triag<bz_triags.size(); ++idx_triag)
{
auto& triag = bz_triags[idx_triag];
ostr << "polygon " << idx_triag << ": " << std::endl;
for(std::size_t idx_vert=0; idx_vert<triag.size(); ++idx_vert)
{
t_vec& vert = triag[idx_vert];
set_eps_0(vert, g_eps);
// find index of vert among voronoi vertices
std::ptrdiff_t voroidx = -1;
if(auto voro_iter = std::find_if(voronoi.begin(), voronoi.end(),
[&vert](const t_vec& vec) -> bool
{
return tl2::equals<t_vec>(vec, vert, g_eps);
}); voro_iter != voronoi.end())
{
voroidx = voro_iter - voronoi.begin();
}
bz_all_triags.push_back(vert);
ostr << "\tvertex " << voroidx << ": "
<< vert << std::endl;
}
}
PlotAddTriangles(bz_all_triags);
// brillouin zone description
m_bz->setPlainText(ostr.str().c_str());
}
/**
* calculate brillouin zone cut
*/
void BZDlg::CalcBZCut()
{
t_real nx = m_cutX->value();
t_real ny = m_cutY->value();
t_real nz = m_cutZ->value();
t_real d = m_cutD->value();
t_vec norm = tl2::create<t_vec>({ nx, ny, nz });
norm /= tl2::norm<t_vec>(norm);
PlotSetPlane(norm, d);
}
tools/bz/bz_plot.cpp
View file @ 290cd478
......@@ -46,6 +46,7 @@ void BZDlg::ShowBZPlot()
m_plot = std::make_shared<tl2::GlPlot>(this);
m_plot->GetRenderer()->SetRestrictCamTheta(false);
m_plot->GetRenderer()->SetCull(false);
m_plot->GetRenderer()->SetBlend(true);
m_plot->GetRenderer()->SetLight(0, tl2::create<t_vec3_gl>({ 5, 5, 5 }));
m_plot->GetRenderer()->SetLight(1, tl2::create<t_vec3_gl>({ -5, -5, -5 }));
m_plot->GetRenderer()->SetCoordMax(1.);
......@@ -57,6 +58,13 @@ void BZDlg::ShowBZPlot()
//comboCoordSys->addItem("Fractional Units (rlu)");
//comboCoordSys->addItem("Lab Units (A)");
m_plot_coordcross = new QCheckBox("Show Coordinates", this);
m_plot_labels = new QCheckBox("Show Labels", this);
m_plot_plane = new QCheckBox("Show Plane", this);
m_plot_coordcross->setChecked(true);
m_plot_labels->setChecked(true);
m_plot_plane->setChecked(true);
m_status3D = new QLabel(this);
m_plot->setSizePolicy(QSizePolicy{QSizePolicy::Expanding, QSizePolicy::Expanding});
......@@ -65,15 +73,21 @@ void BZDlg::ShowBZPlot()
auto grid = new QGridLayout(m_dlgPlot);
grid->setSpacing(2);
grid->setContentsMargins(4,4,4,4);
grid->addWidget(m_plot.get(), 0,0,1,2);
grid->addWidget(m_plot.get(), 0,0,1,3);
//grid->addWidget(labCoordSys, 1,0,1,1);
//grid->addWidget(comboCoordSys, 1,1,1,1);
grid->addWidget(m_status3D, 2,0,1,2);
grid->addWidget(m_plot_coordcross, 1,0,1,1);
grid->addWidget(m_plot_labels, 1,1,1,1);
grid->addWidget(m_plot_plane, 1,2,1,1);
grid->addWidget(m_status3D, 2,0,1,3);
connect(m_plot.get(), &tl2::GlPlot::AfterGLInitialisation, this, &BZDlg::AfterGLInitialisation);
connect(m_plot->GetRenderer(), &tl2::GlPlotRenderer::PickerIntersection, this, &BZDlg::PickerIntersection);
connect(m_plot.get(), &tl2::GlPlot::MouseDown, this, &BZDlg::PlotMouseDown);
connect(m_plot.get(), &tl2::GlPlot::MouseUp, this, &BZDlg::PlotMouseUp);
connect(m_plot_coordcross, &QCheckBox::toggled, this, &BZDlg::PlotShowCoordCross);
connect(m_plot_labels, &QCheckBox::toggled, this, &BZDlg::PlotShowLabels);
connect(m_plot_plane, &QCheckBox::toggled, this, &BZDlg::PlotShowPlane);
/*connect(comboCoordSys, static_cast<void (QComboBox::*)(int)>(&QComboBox::currentIndexChanged), this, [this](int val)
{
if(this->m_plot)
......@@ -93,6 +107,45 @@ void BZDlg::ShowBZPlot()
}
/**
* show or hide the coordinate system
*/
void BZDlg::PlotShowCoordCross(bool show)
{
if(!m_plot) return;
if(auto obj = m_plot->GetRenderer()->GetCoordCross(); obj)
{
m_plot->GetRenderer()->SetObjectVisible(*obj, show);
m_plot->update();
}
}
/**
* show or hide the object labels
*/
void BZDlg::PlotShowLabels(bool show)
{
if(!m_plot) return;
m_plot->GetRenderer()->SetLabelsVisible(show);
m_plot->update();
}
/**
* show or hide the BZ cut plane
*/
void BZDlg::PlotShowPlane(bool show)
{
if(!m_plot) return;
m_plot->GetRenderer()->SetObjectVisible(m_plane, show);
m_plot->update();
}
/**
* add a voronoi vertex to the plot
*/
......@@ -192,6 +245,26 @@ void BZDlg::PlotAddTriangles(const std::vector<t_vec>& _vecs)
}
/**
* set the brillouin zone cut plane
*/
void BZDlg::PlotSetPlane(const t_vec& _norm, t_real d)
{
if(!m_plot) return;
t_vec3_gl norm = tl2::convert<t_vec3_gl>(_norm);
t_vec3_gl norm_old = tl2::create<t_vec3_gl>({ 0, 0, 1 });
t_vec3_gl rot_vec = tl2::create<t_vec3_gl>({ 1, 0, 0 });
t_vec3_gl offs = d * norm;
t_mat_gl rot = tl2::hom_rotation<t_mat_gl, t_vec3_gl>(norm_old, norm, &rot_vec);
t_mat_gl trans = tl2::hom_translation<t_mat_gl>(offs[0], offs[1], offs[2]);
m_plot->GetRenderer()->SetObjectMatrix(m_plane, trans*rot);
m_plot->update();
}