XtalOptimizer: Option to refine unit cell

This commit is contained in:
2026-01-27 14:57:29 +01:00
parent 626f7abc30
commit 52c3c3f91a
2 changed files with 30 additions and 18 deletions
@@ -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<float>(latt_vec0[0]), static_cast<float>(latt_vec0[1]), static_cast<float>(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) {
@@ -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;