cont. of the previous commit (symmetry in SmallStructure)

add '.' to the possible characters in the order string
remove SmallStructure::find_spacegroup()
initialize spacegroup_number to 0

docs/mol.rst

@@ -192,32 +192,42 @@ the space group are read and stored in member variables:
     >>> st.spacegroup_number
     164
 
-and the function `set_spacegroup("SH2")` is automatically
+and the function `set_spacegroup("S.H2")` is automatically
 run to set `spacegroup`:
 
 .. doctest::
 
     >>> st.spacegroup
     <gemmi.SpaceGroup("P -3 m 1")>
 
-`set_spacegroup()` takes one argument, a string in which each character
-specifies what to use for space group determination:
+`set_spacegroup()` takes one argument, a string in which characters
+specify what to use, and in what order, for space group determination:
 
 * `S` = symmetry operations stored in `symops`,
 * `H` = Hall symbol from `spacegroup_hall` (we compare symmetry operations
   encoded in the Hall symbol, not the strings),
 * `1` = H-M symbol; for space groups such as "P n n n" that have two origin
   choices listed in the International Tables, use *Origin Choice 1*,
 * `2` = H-M symbol, with *Origin Choice 2* where applicable,
-* `N` = the space group number.
+* `N` = the space group number,
+* `.` (after S or H) = if the symmetry operations pass sanity checks,
+  stop and use them regardless of whether they correspond to one of
+  the settings tabulated in Gemmi.
 
-The first item that matches one of the 560+ space group settings tabulated
-in Gemmi sets `spacegroup`. To use a different order of items than SH2,
-call set_spacegroup() again:
+If a symbol or operations match one of the 560+ space group settings tabulated
+in Gemmi, `spacegroup` is set to this setting. Otherwise, if `.` is encountered
+and the previous character (`S` or `H`) was evaluated to a valid set of symops,
+it is assumed that these operations were correct: `spacegroup` is left null
+and `cell.images` are set from the list of operations.
+About 350 (out of 500,000+) entries in the COD use such settings.
+Most of them have an unconventional choice of the origin
+(e.g. "P 1 21 1 (a,b,c-1/4)").
+
+To use a different order of items than "S.H2", call set_spacegroup() again:
 
 .. doctest::
 
-    >>> st.set_spacegroup('S1')
+    >>> st.set_spacegroup('H.1')
 
 Errors such as an incorrect format of the symop triplets or of the Hall
 symbol are silently ignored, and the consistency between different items
@@ -248,8 +258,8 @@ would be used to make gemmi::GroupOps from symops::
     GroupOps split_centering_vectors(const std::vector<Op>& ops)
 
 
-SmallStructure <-> Structure
-----------------------------
+without CIF file
+----------------
 
 If your structure is stored in a macromolecular format (PDB, mmCIF)
 you can read it first as macromolecular :ref:`hierarchy <mcra>`

docs/symmetry.rst

@@ -8,7 +8,7 @@ In C++: `#include <gemmi/symmetry.hpp>` (it's header-only).
 Space group table
 =================
 
-Gemmi tabulates 550+ settings of the 230 crystallographic space groups.
+Gemmi tabulates 560+ settings of the 230 crystallographic space groups.
 Each entry includes:
 
 * `number` -- space group number (1-230),

include/gemmi/small.hpp

@@ -19,6 +19,22 @@
 
 namespace gemmi {
 
+inline bool is_complete(const GroupOps& gops) {
+  for (Op op1 : gops.sym_ops)
+    for (Op op2 : gops.sym_ops)
+      if (gops.find_by_rotation((op1 * op2).rot) == nullptr)
+        return false;
+  return true;
+}
+
+inline std::vector<Op> triplets_to_ops(const std::vector<std::string>& symops) {
+  std::vector<Op> ops;
+  ops.reserve(symops.size());
+  for (const std::string& xyz : symops)
+    ops.push_back(parse_triplet(xyz));
+  return ops;
+}
+
 struct SmallStructure {
   struct Site {
     std::string label;
@@ -58,43 +74,42 @@ struct SmallStructure {
   const SpaceGroup* spacegroup = nullptr;
   std::string spacegroup_hm;
   std::string spacegroup_hall;
-  int spacegroup_number;
+  int spacegroup_number = 0;
   std::vector<std::string> symops;
   std::vector<Site> sites;
   std::vector<AtomType> atom_types;
   double wavelength = 0.; // the first wavelength if multiple
 
   std::vector<Site> get_all_unit_cell_sites() const;
 
-  // deprecated, use directly spacegroup
-  const SpaceGroup* find_spacegroup() const { return spacegroup; }
-
-  const SpaceGroup* find_spacegroup_from_symops() const {
-    if (symops.empty())
-      return nullptr;
-    std::vector<Op> ops;
-    ops.reserve(symops.size());
-    for (const std::string& xyz : symops)
-      ops.push_back(parse_triplet(xyz));
-    GroupOps gops = split_centering_vectors(ops);
-    return find_spacegroup_by_ops(gops);
-  }
-
   void set_spacegroup(const char* order) {
     spacegroup = nullptr;
     if (order)
       for (const char* c = order; *c != '\0' && spacegroup == nullptr; ++c) {
         try {
-          spacegroup = get_spacegroup_from(*c);
+          GroupOps gops;
+          spacegroup = get_spacegroup_from(*c, gops);
+          if (!spacegroup && *(c+1) == '.') {
+            // If symops don't correspond to tabulated settings,
+            // we can't set spacegroup, but we can set UnitCell::images.
+            if (gops.order() == (int) symops.size() && is_complete(gops)) {
+              cell.set_cell_images_from_groupops(gops);
+              return;
+            }
+            ++c;
+          }
         } catch (std::exception&) {}
       }
     setup_cell_images();
   }
 
-  const SpaceGroup* get_spacegroup_from(char c) const {
+  const SpaceGroup* get_spacegroup_from(char c, GroupOps& gops) const {
     switch (lower(c)) {
       case 's':
-        return find_spacegroup_from_symops();
+        if (symops.empty())
+          return nullptr;
+        gops = split_centering_vectors(triplets_to_ops(symops));
+        return find_spacegroup_by_ops(gops);
       case 'h':
         if (spacegroup_hall.empty())
           return nullptr;
@@ -119,7 +134,16 @@ struct SmallStructure {
     std::string err;
     if (!symops.empty())
       try {
-        auto sg = find_spacegroup_from_symops();
+        std::vector<Op> ops = triplets_to_ops(symops);
+        for (Op& op : ops)
+          op.wrap();
+        std::sort(ops.begin(), ops.end());
+        GroupOps gops = split_centering_vectors(ops);
+        if (!is_complete(gops))
+          cat_to(err, "symops list is incomplete or incorrect\n");
+        else if (gops.all_ops_sorted() != ops)
+          cat_to(err, "symops list is incorrect or incomplete or redundant\n");
+        auto sg = find_spacegroup_by_ops(gops);
         if (!sg)
           cat_to(err, "space group from symops not found in the table\n");
         else if (sg != spacegroup)
@@ -182,7 +206,7 @@ struct SmallStructure {
   }
 
   void setup_cell_images() {
-    cell.set_cell_images_from_spacegroup(find_spacegroup());
+    cell.set_cell_images_from_spacegroup(spacegroup);
   }
 };
 

include/gemmi/smcif.hpp

@@ -52,7 +52,7 @@ SmallStructure make_small_structure_from_block(const cif::Block& block_) {
       st.spacegroup_number = cif::as_int(*val, 0);
       break;
     }
-  st.set_spacegroup("SH2");
+  st.set_spacegroup("S.H2");
 
   enum { kLabel, kSymbol, kX, kY, kZ, kUiso, kBiso, kOcc, kDisorderGroup };
   cif::Table atom_table = block.find("_atom_site_",

include/gemmi/unitcell.hpp

@@ -335,18 +335,21 @@ struct UnitCell {
     return sg ? is_compatible_with_groupops(sg->operations(), eps) : false;
   }
 
-  void set_cell_images_from_spacegroup(const SpaceGroup* sg) {
+  void set_cell_images_from_groupops(const GroupOps& group_ops) {
     images.clear();
-    cs_count = 0;
-    if (!sg)
-      return;
-    GroupOps group_ops = sg->operations();
     cs_count = (short) group_ops.order() - 1;
     images.reserve(cs_count);
-    for (Op op : group_ops) {
-      if (op == Op::identity())
-        continue;
-      images.push_back(Transform{rot_as_mat33(op), tran_as_vec3(op)});
+    for (Op op : group_ops)
+      if (op != Op::identity())
+        images.push_back(Transform{rot_as_mat33(op), tran_as_vec3(op)});
+  }
+
+  void set_cell_images_from_spacegroup(const SpaceGroup* sg) {
+    if (sg) {
+      set_cell_images_from_groupops(sg->operations());
+    } else {
+      images.clear();
+      cs_count = 0;
     }
   }
 
@@ -356,9 +359,9 @@ struct UnitCell {
       if (!ncs_op.given) {
         // We need it to operates on fractional, not orthogonal coordinates.
         FTransform f = frac.combine(ncs_op.tr.combine(orth));
-        images.emplace_back(f);
+        images.push_back(f);
         for (int i = 0; i < cs_count; ++i)
-          images.emplace_back(images[i].combine(f));
+          images.push_back(images[i].combine(f));
       }
   }
 

prog/sfcalc.cpp

@@ -380,7 +380,7 @@ void print_structure_factors_sm(const gemmi::SmallStructure& small,
   int max_h = int(max_1_d / small.cell.ar);
   int max_k = int(max_1_d / small.cell.br);
   int max_l = int(max_1_d / small.cell.cr);
-  const gemmi::SpaceGroup* sg = small.find_spacegroup();
+  const gemmi::SpaceGroup* sg = small.spacegroup;
   if (!sg)
     sg = &gemmi::get_spacegroup_p1();
   gemmi::ReciprocalAsu asu(sg);