added function calculate_b_aniso_range()

replaces undocumented function get_minimum_b()

docs/mol.rst

@@ -2706,7 +2706,8 @@ so just in case we store it as a string.
   >>> model.name
   '1'
 
-----
+Subchains
+---------
 
 As was discussed before, the PDBx/mmCIF format has also
 a set of parallel identifiers. In particular, it has
@@ -2747,7 +2748,8 @@ The `ResidueSpan` is described in the next section.
 ..
   TODO: find_residue_group, sole_residue, get_all_residue_names
 
-----
+Helper functions
+----------------
 
 .. _model_count_atom:
 
@@ -2765,8 +2767,6 @@ In Python, `Model` has also methods for often needed calculations:
   <gemmi.Position(-5.7572, 16.4099, 2.88299)>
   >>> model.has_hydrogen()
   False
-  >>> model.calculate_b_iso_range()
-  (7.670000076293945, 46.880001068115234)
 
 The first two function can take a :ref:`Selection <selections>`
 as an argument. For example, we can count sulfur atoms with:
@@ -2778,6 +2778,18 @@ as an argument. For example, we can count sulfur atoms with:
   >>> model.count_occupancies(gemmi.Selection('[S]'))
   2.0
 
+Two functions calculate the range of ADP (B-factor) values in the model.
+One function considers only isotropic B values, while the other uses
+minimum and maximum eigenvalues of anisotropic ADPs. For atoms lacking
+ANISOU records, it falls back to the isotropic B-factor.
+
+.. doctest::
+
+  >>> model.calculate_b_iso_range()  # doctest: +ELLIPSIS
+  (7.67000..., 46.88000...)
+  >>> model.calculate_b_aniso_range()  # doctest: +ELLIPSIS
+  (3.523999..., 122.568275...)
+
 In C++, the same functionality is provided by templated functions
 from `gemmi/calculate.hpp` and `gemmi/select.hpp`.
 These functions (in C++) can be applied not only

include/gemmi/calculate.hpp

@@ -6,6 +6,7 @@
 #define GEMMI_CALCULATE_HPP_
 
 #include <array>
+#include <algorithm>  // for std::min, std::minmax
 #include "model.hpp"
 
 namespace gemmi {
@@ -73,6 +74,28 @@ template<> inline std::pair<float,float> calculate_b_iso_range(const Atom& atom)
   return {atom.b_iso, atom.b_iso};
 }
 
+/// uses min/max eigenvalues of Baniso, or Biso if B-factor is isotropic
+inline std::pair<double, double> calculate_b_aniso_range(const Model& model) {
+  std::pair<double, double> range{INFINITY, -INFINITY};
+  for (const Chain& chain : model.chains)
+    for (const Residue& residue : chain.residues)
+      for (const Atom& atom : residue.atoms) {
+        if (atom.occ == 0)
+          continue;
+        if (atom.aniso.nonzero()) {
+          std::array<double,3> eig = atom.aniso.calculate_eigenvalues();
+          auto u = std::minmax({eig[0], eig[1], eig[2]});
+          range.first = std::min(range.first, u.first * u_to_b());
+          range.second = std::max(range.second, u.second * u_to_b());
+        } else {
+          range.first = std::min(range.first, (double) atom.b_iso);
+          range.second = std::max(range.second, (double) atom.b_iso);
+        }
+      }
+  return range;
+}
+
+
 template<class T> void expand_box(const T& obj, Box<Position>& box) {
   for (const auto& child : obj.children())
     expand_box(child, box);

include/gemmi/dencalc.hpp

@@ -10,6 +10,7 @@
 #include "formfact.hpp" // for ExpSum
 #include "grid.hpp"     // for Grid
 #include "model.hpp"    // for Structure, ...
+#include "calculate.hpp" // for calculate_b_aniso_range
 
 namespace gemmi {
 
@@ -66,23 +67,6 @@ Real it92_radius_approx(Real b) {
   return (8.5f + 0.075f * b) / (2.4f + 0.0045f * b);
 }
 
-inline double get_minimum_b(const Model& model) {
-  double b_min = 1000.;
-  for (const Chain& chain : model.chains)
-    for (const Residue& residue : chain.residues)
-      for (const Atom& atom : residue.atoms) {
-        if (atom.occ == 0) continue;
-        double b = atom.b_iso;
-        if (atom.aniso.nonzero()) {
-          std::array<double,3> eig = atom.aniso.calculate_eigenvalues();
-          b = std::min(std::min(eig[0], eig[1]), eig[2]) * u_to_b();
-        }
-        if (b < b_min)
-          b_min = b;
-      }
-  return b_min;
-}
-
 // Usual usage:
 // - set d_min and optionally also other parameters,
 // - set addends to f' values for your wavelength (see fprime.hpp)
@@ -112,7 +96,7 @@ struct DensityCalculator {
     double spacing = requested_grid_spacing();
     if (spacing <= 0)
       spacing = std::min(std::min(grid.spacing[0], grid.spacing[1]), grid.spacing[2]);
-    double b_min = get_minimum_b(model);
+    double b_min = calculate_b_aniso_range(model).first;
     blur = std::max(u_to_b() / 1.1 * sq(spacing) - b_min, 0.);
   }
 

prog/sfcalc.cpp

@@ -661,9 +661,10 @@ void process_with_table(bool use_st, gemmi::Structure& st, const gemmi::SmallStr
         // and on the atomic cutoff radius, and probably it would be too
         // hard to estimate. Here we use the same formula as in Refmac.
         dencalc.set_refmac_compatible_blur(st.models[0]);
-        if (p.options[Verbose])
-          fprintf(stderr, "B_min=%g, B_add=%g\n",
-                  gemmi::get_minimum_b(st.models[0]), dencalc.blur);
+        if (p.options[Verbose]) {
+          double b_min = gemmi::calculate_b_aniso_range(st.models[0]).first;
+          fprintf(stderr, "B_min=%g, B_add=%g\n", b_min, dencalc.blur);
+        }
       }
       gemmi::AtomicRadiiSet radii_choice = gemmi::AtomicRadiiSet::VanDerWaals;
       if (p.options[RadiiSet]) {

python/mol.cpp

@@ -269,6 +269,7 @@ void add_mol(py::module& m) {
         return calculate_center_of_mass(self).get();
     })
     .def("calculate_b_iso_range", &calculate_b_iso_range<Model>)
+    .def("calculate_b_aniso_range", &calculate_b_aniso_range)
     .def("transform_pos_and_adp", transform_pos_and_adp<Model>, py::arg("tr"))
     .def("split_chains_by_segments", &split_chains_by_segments)
     .def("clone", [](const Model& self) { return new Model(self); })