From 52c3c3f91ac4cdae58361a5e1b6e679be70fbfaa Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Tue, 27 Jan 2026 14:57:29 +0100 Subject: [PATCH] XtalOptimizer: Option to refine unit cell --- .../geom_refinement/XtalOptimizer.cpp | 47 ++++++++++++------- .../geom_refinement/XtalOptimizer.h | 1 + 2 files changed, 30 insertions(+), 18 deletions(-) diff --git a/image_analysis/geom_refinement/XtalOptimizer.cpp b/image_analysis/geom_refinement/XtalOptimizer.cpp index f1e775df..ccf2eecc 100644 --- a/image_analysis/geom_refinement/XtalOptimizer.cpp +++ b/image_analysis/geom_refinement/XtalOptimizer.cpp @@ -500,25 +500,30 @@ bool XtalOptimizerInternal(XtalOptimizerData &data, } } - // Parameter bounds - // Lengths - for (int i = 0; i < 3; ++i) { - problem.SetParameterLowerBound(latt_vec1, i, data.min_length_A); - problem.SetParameterUpperBound(latt_vec1, i, data.max_length_A); - } - - if (data.crystal_system == gemmi::CrystalSystem::Monoclinic) { - const double beta_lo = std::max(1e-6, M_PI * (data.min_angle_deg / 180.0)); - const double beta_hi = std::min(M_PI - 1e-6, M_PI * (data.max_angle_deg / 180.0)); - problem.SetParameterLowerBound(latt_vec2, 0, beta_lo); - problem.SetParameterUpperBound(latt_vec2, 0, beta_hi); - } else if (data.crystal_system == gemmi::CrystalSystem::Triclinic) { - // α, β, γ bounds (radians) - const double alo = M_PI * (data.min_angle_deg / 180.0); - const double ahi = M_PI * (data.max_angle_deg / 180.0); + if (!data.refine_unit_cell) { + problem.SetParameterBlockConstant(latt_vec1); + problem.SetParameterBlockConstant(latt_vec2); + } else { + // Parameter bounds + // Lengths for (int i = 0; i < 3; ++i) { - problem.SetParameterLowerBound(latt_vec2, i, alo); - problem.SetParameterUpperBound(latt_vec2, i, ahi); + problem.SetParameterLowerBound(latt_vec1, i, data.min_length_A); + problem.SetParameterUpperBound(latt_vec1, i, data.max_length_A); + } + + if (data.crystal_system == gemmi::CrystalSystem::Monoclinic) { + const double beta_lo = std::max(1e-6, M_PI * (data.min_angle_deg / 180.0)); + const double beta_hi = std::min(M_PI - 1e-6, M_PI * (data.max_angle_deg / 180.0)); + problem.SetParameterLowerBound(latt_vec2, 0, beta_lo); + problem.SetParameterUpperBound(latt_vec2, 0, beta_hi); + } else if (data.crystal_system == gemmi::CrystalSystem::Triclinic) { + // α, β, γ bounds (radians) + const double alo = M_PI * (data.min_angle_deg / 180.0); + const double ahi = M_PI * (data.max_angle_deg / 180.0); + for (int i = 0; i < 3; ++i) { + problem.SetParameterLowerBound(latt_vec2, i, alo); + problem.SetParameterUpperBound(latt_vec2, i, ahi); + } } } @@ -545,6 +550,12 @@ bool XtalOptimizerInternal(XtalOptimizerData &data, if (data.refine_detector_angles) data.geom.PoniRot1_rad(detector_rot[0]).PoniRot2_rad(detector_rot[1]); + if (!data.refine_unit_cell) { + Coord rot_vector{static_cast(latt_vec0[0]), static_cast(latt_vec0[1]), static_cast(latt_vec0[2])}; + Coord rot_vector_norm = rot_vector.Normalize(); + std::cout << "X " << rot_vector.Length() * 180.0 / M_PI << " " << rot_vector_norm << std::endl; + } + if (data.crystal_system == gemmi::CrystalSystem::Orthorhombic) data.latt = AngleAxisAndLengthsToLattice(latt_vec0, latt_vec1, false); else if (data.crystal_system == gemmi::CrystalSystem::Tetragonal) { diff --git a/image_analysis/geom_refinement/XtalOptimizer.h b/image_analysis/geom_refinement/XtalOptimizer.h index 38c53cd3..1854459c 100644 --- a/image_analysis/geom_refinement/XtalOptimizer.h +++ b/image_analysis/geom_refinement/XtalOptimizer.h @@ -26,6 +26,7 @@ struct XtalOptimizerData { bool refine_beam_center = true; bool refine_distance_mm = false; bool refine_detector_angles = false; + bool refine_unit_cell = true; // This refines unit cell size + angles - orientation is always refined bool index_ice_rings = true;