From e1c661e48db1e1d869742c1060dda47af08de8ab Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Wed, 13 May 2026 13:32:25 +0200 Subject: [PATCH] Revert "Gemmi: Include through FetchContent full gemmi library (not limited cpp/hpp files)" This reverts commit b06dfc8357e9706d96f55eb0991e48aaeabdcae2. --- CMakeLists.txt | 11 +- common/CMakeLists.txt | 2 +- docs/SOFTWARE.md | 2 +- image_analysis/CMakeLists.txt | 2 +- symmetry/CMakeLists.txt | 2 + symmetry/LICENSE.txt | 373 ++++++++++ symmetry/gemmi/cellred.hpp | 406 +++++++++++ symmetry/gemmi/fail.hpp | 93 +++ symmetry/gemmi/math.hpp | 458 +++++++++++++ symmetry/gemmi/symmetry.hpp | 1044 ++++++++++++++++++++++++++++ symmetry/gemmi/unitcell.hpp | 618 +++++++++++++++++ symmetry/symmetry.cpp | 1215 +++++++++++++++++++++++++++++++++ 12 files changed, 4214 insertions(+), 12 deletions(-) create mode 100644 symmetry/CMakeLists.txt create mode 100644 symmetry/LICENSE.txt create mode 100644 symmetry/gemmi/cellred.hpp create mode 100644 symmetry/gemmi/fail.hpp create mode 100644 symmetry/gemmi/math.hpp create mode 100644 symmetry/gemmi/symmetry.hpp create mode 100644 symmetry/gemmi/unitcell.hpp create mode 100644 symmetry/symmetry.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index 9beada6d0..c49e06444 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -95,13 +95,6 @@ SET(BUILD_SHARED_LIBS OFF CACHE BOOL "" FORCE) SET(HTTPLIB_USE_NON_BLOCKING_GETADDRINFO OFF CACHE BOOL "" FORCE) SET(HTTPLIB_REQUIRE_ZLIB ON CACHE BOOL "" FORCE) -FetchContent_Declare( - gemmi - GIT_REPOSITORY https://github.com/fleon-psi/gemmi - GIT_TAG d6dcc1f57eedf7ba34a7d2d2ed283075113040bf - EXCLUDE_FROM_ALL -) - FetchContent_Declare( spdlog GIT_REPOSITORY https://github.com/gabime/spdlog.git @@ -149,7 +142,7 @@ FetchContent_Declare( EXCLUDE_FROM_ALL ) -FetchContent_MakeAvailable(zstd sls_detector_package catch2 hdf5 spdlog httplib gemmi) +FetchContent_MakeAvailable(zstd sls_detector_package catch2 hdf5 spdlog httplib) ADD_SUBDIRECTORY(jungfrau) ADD_SUBDIRECTORY(compression) @@ -160,7 +153,7 @@ ADD_SUBDIRECTORY(reader) ADD_SUBDIRECTORY(detector_control) ADD_SUBDIRECTORY(image_puller) ADD_SUBDIRECTORY(preview) -#ADD_SUBDIRECTORY(symmetry) +ADD_SUBDIRECTORY(symmetry) ADD_SUBDIRECTORY(xds-plugin) IF (JFJOCH_WRITER_ONLY) diff --git a/common/CMakeLists.txt b/common/CMakeLists.txt index d9b1d8b23..a15a652bd 100644 --- a/common/CMakeLists.txt +++ b/common/CMakeLists.txt @@ -130,7 +130,7 @@ ADD_LIBRARY(JFJochCommon STATIC ScalingSettings.h ) -TARGET_LINK_LIBRARIES(JFJochCommon JFJochLogger Compression JFCalibration gemmi_cpp Threads::Threads -lrt ) +TARGET_LINK_LIBRARIES(JFJochCommon JFJochLogger Compression JFCalibration gemmi Threads::Threads -lrt ) TARGET_LINK_LIBRARIES(JFJochZMQ "$") diff --git a/docs/SOFTWARE.md b/docs/SOFTWARE.md index bddc39521..d64671e8e 100644 --- a/docs/SOFTWARE.md +++ b/docs/SOFTWARE.md @@ -33,7 +33,6 @@ Automatically downloaded by CMake and statically linked: * Catch2 testing library - see [github.com/catchorg/Catch2](https://github.com/catchorg/Catch2) * Ceres Solver library for least square optimization - see [http://ceres-solver.org/] * Spdlog logging library - see [github.com/gabime/spdlog](https://github.com/gabime/spdlog) -* GEMMI library by Global Phasing - see [github.com/project-gemmi/gemmi](https://github.com/project-gemmi/gemmi) Please follow the link provided above to check for LICENSE file. Building code with dependencies above requires access from the build system to github.com. Directly included in the repository: @@ -45,5 +44,6 @@ Directly included in the repository: * LZ4 compression by Y.Collet - see [github.com/lz4/lz4](https://github.com/lz4/lz4) * ZeroMQ library (through slsDetectorPackage) - see [github.com/zeromq/libzmq](https://github.com/zeromq/libzmq) * Base64 decoder/encoder - see [gist.github.com/tomykaira](https://gist.github.com/tomykaira/f0fd86b6c73063283afe550bc5d77594) +* GEMMI library by Global Phasing - see [github.com/project-gemmi/gemmi](https://github.com/project-gemmi/gemmi) For license check LICENSE file in respective directory diff --git a/image_analysis/CMakeLists.txt b/image_analysis/CMakeLists.txt index 6bb525187..394dbaed7 100644 --- a/image_analysis/CMakeLists.txt +++ b/image_analysis/CMakeLists.txt @@ -44,4 +44,4 @@ ADD_SUBDIRECTORY(scale_merge) ADD_SUBDIRECTORY(image_preprocessing) ADD_SUBDIRECTORY(azint) -TARGET_LINK_LIBRARIES(JFJochImageAnalysis JFJochAzIntEngine JFJochImagePreprocessing JFJochBraggPrediction JFJochBraggIntegration JFJochLatticeSearch JFJochIndexing JFJochSpotFinding JFJochCommon JFJochGeomRefinement JFJochScaleMerge) +TARGET_LINK_LIBRARIES(JFJochImageAnalysis JFJochAzIntEngine JFJochImagePreprocessing JFJochBraggPrediction JFJochBraggIntegration JFJochLatticeSearch JFJochIndexing JFJochSpotFinding JFJochCommon JFJochGeomRefinement JFJochScaleMerge gemmi) diff --git a/symmetry/CMakeLists.txt b/symmetry/CMakeLists.txt new file mode 100644 index 000000000..fed3f7925 --- /dev/null +++ b/symmetry/CMakeLists.txt @@ -0,0 +1,2 @@ +ADD_LIBRARY(gemmi STATIC symmetry.cpp gemmi/symmetry.hpp gemmi/fail.hpp) +TARGET_INCLUDE_DIRECTORIES(gemmi PUBLIC .) \ No newline at end of file diff --git a/symmetry/LICENSE.txt b/symmetry/LICENSE.txt new file mode 100644 index 000000000..14e2f777f --- /dev/null +++ b/symmetry/LICENSE.txt @@ -0,0 +1,373 @@ +Mozilla Public License Version 2.0 +================================== + +1. Definitions +-------------- + +1.1. "Contributor" + means each individual or legal entity that creates, contributes to + the creation of, or owns Covered Software. + +1.2. "Contributor Version" + means the combination of the Contributions of others (if any) used + by a Contributor and that particular Contributor's Contribution. + +1.3. "Contribution" + means Covered Software of a particular Contributor. + +1.4. "Covered Software" + means Source Code Form to which the initial Contributor has attached + the notice in Exhibit A, the Executable Form of such Source Code + Form, and Modifications of such Source Code Form, in each case + including portions thereof. + +1.5. "Incompatible With Secondary Licenses" + means + + (a) that the initial Contributor has attached the notice described + in Exhibit B to the Covered Software; or + + (b) that the Covered Software was made available under the terms of + version 1.1 or earlier of the License, but not also under the + terms of a Secondary License. + +1.6. "Executable Form" + means any form of the work other than Source Code Form. + +1.7. "Larger Work" + means a work that combines Covered Software with other material, in + a separate file or files, that is not Covered Software. + +1.8. "License" + means this document. + +1.9. "Licensable" + means having the right to grant, to the maximum extent possible, + whether at the time of the initial grant or subsequently, any and + all of the rights conveyed by this License. + +1.10. "Modifications" + means any of the following: + + (a) any file in Source Code Form that results from an addition to, + deletion from, or modification of the contents of Covered + Software; or + + (b) any new file in Source Code Form that contains any Covered + Software. + +1.11. "Patent Claims" of a Contributor + means any patent claim(s), including without limitation, method, + process, and apparatus claims, in any patent Licensable by such + Contributor that would be infringed, but for the grant of the + License, by the making, using, selling, offering for sale, having + made, import, or transfer of either its Contributions or its + Contributor Version. + +1.12. "Secondary License" + means either the GNU General Public License, Version 2.0, the GNU + Lesser General Public License, Version 2.1, the GNU Affero General + Public License, Version 3.0, or any later versions of those + licenses. + +1.13. "Source Code Form" + means the form of the work preferred for making modifications. + +1.14. "You" (or "Your") + means an individual or a legal entity exercising rights under this + License. For legal entities, "You" includes any entity that + controls, is controlled by, or is under common control with You. For + purposes of this definition, "control" means (a) the power, direct + or indirect, to cause the direction or management of such entity, + whether by contract or otherwise, or (b) ownership of more than + fifty percent (50%) of the outstanding shares or beneficial + ownership of such entity. + +2. License Grants and Conditions +-------------------------------- + +2.1. Grants + +Each Contributor hereby grants You a world-wide, royalty-free, +non-exclusive license: + +(a) under intellectual property rights (other than patent or trademark) + Licensable by such Contributor to use, reproduce, make available, + modify, display, perform, distribute, and otherwise exploit its + Contributions, either on an unmodified basis, with Modifications, or + as part of a Larger Work; and + +(b) under Patent Claims of such Contributor to make, use, sell, offer + for sale, have made, import, and otherwise transfer either its + Contributions or its Contributor Version. + +2.2. Effective Date + +The licenses granted in Section 2.1 with respect to any Contribution +become effective for each Contribution on the date the Contributor first +distributes such Contribution. + +2.3. Limitations on Grant Scope + +The licenses granted in this Section 2 are the only rights granted under +this License. No additional rights or licenses will be implied from the +distribution or licensing of Covered Software under this License. +Notwithstanding Section 2.1(b) above, no patent license is granted by a +Contributor: + +(a) for any code that a Contributor has removed from Covered Software; + or + +(b) for infringements caused by: (i) Your and any other third party's + modifications of Covered Software, or (ii) the combination of its + Contributions with other software (except as part of its Contributor + Version); or + +(c) under Patent Claims infringed by Covered Software in the absence of + its Contributions. + +This License does not grant any rights in the trademarks, service marks, +or logos of any Contributor (except as may be necessary to comply with +the notice requirements in Section 3.4). + +2.4. Subsequent Licenses + +No Contributor makes additional grants as a result of Your choice to +distribute the Covered Software under a subsequent version of this +License (see Section 10.2) or under the terms of a Secondary License (if +permitted under the terms of Section 3.3). + +2.5. Representation + +Each Contributor represents that the Contributor believes its +Contributions are its original creation(s) or it has sufficient rights +to grant the rights to its Contributions conveyed by this License. + +2.6. Fair Use + +This License is not intended to limit any rights You have under +applicable copyright doctrines of fair use, fair dealing, or other +equivalents. + +2.7. Conditions + +Sections 3.1, 3.2, 3.3, and 3.4 are conditions of the licenses granted +in Section 2.1. + +3. Responsibilities +------------------- + +3.1. Distribution of Source Form + +All distribution of Covered Software in Source Code Form, including any +Modifications that You create or to which You contribute, must be under +the terms of this License. You must inform recipients that the Source +Code Form of the Covered Software is governed by the terms of this +License, and how they can obtain a copy of this License. You may not +attempt to alter or restrict the recipients' rights in the Source Code +Form. + +3.2. Distribution of Executable Form + +If You distribute Covered Software in Executable Form then: + +(a) such Covered Software must also be made available in Source Code + Form, as described in Section 3.1, and You must inform recipients of + the Executable Form how they can obtain a copy of such Source Code + Form by reasonable means in a timely manner, at a charge no more + than the cost of distribution to the recipient; and + +(b) You may distribute such Executable Form under the terms of this + License, or sublicense it under different terms, provided that the + license for the Executable Form does not attempt to limit or alter + the recipients' rights in the Source Code Form under this License. + +3.3. Distribution of a Larger Work + +You may create and distribute a Larger Work under terms of Your choice, +provided that You also comply with the requirements of this License for +the Covered Software. If the Larger Work is a combination of Covered +Software with a work governed by one or more Secondary Licenses, and the +Covered Software is not Incompatible With Secondary Licenses, this +License permits You to additionally distribute such Covered Software +under the terms of such Secondary License(s), so that the recipient of +the Larger Work may, at their option, further distribute the Covered +Software under the terms of either this License or such Secondary +License(s). + +3.4. Notices + +You may not remove or alter the substance of any license notices +(including copyright notices, patent notices, disclaimers of warranty, +or limitations of liability) contained within the Source Code Form of +the Covered Software, except that You may alter any license notices to +the extent required to remedy known factual inaccuracies. + +3.5. Application of Additional Terms + +You may choose to offer, and to charge a fee for, warranty, support, +indemnity or liability obligations to one or more recipients of Covered +Software. However, You may do so only on Your own behalf, and not on +behalf of any Contributor. You must make it absolutely clear that any +such warranty, support, indemnity, or liability obligation is offered by +You alone, and You hereby agree to indemnify every Contributor for any +liability incurred by such Contributor as a result of warranty, support, +indemnity or liability terms You offer. You may include additional +disclaimers of warranty and limitations of liability specific to any +jurisdiction. + +4. Inability to Comply Due to Statute or Regulation +--------------------------------------------------- + +If it is impossible for You to comply with any of the terms of this +License with respect to some or all of the Covered Software due to +statute, judicial order, or regulation then You must: (a) comply with +the terms of this License to the maximum extent possible; and (b) +describe the limitations and the code they affect. Such description must +be placed in a text file included with all distributions of the Covered +Software under this License. Except to the extent prohibited by statute +or regulation, such description must be sufficiently detailed for a +recipient of ordinary skill to be able to understand it. + +5. Termination +-------------- + +5.1. The rights granted under this License will terminate automatically +if You fail to comply with any of its terms. However, if You become +compliant, then the rights granted under this License from a particular +Contributor are reinstated (a) provisionally, unless and until such +Contributor explicitly and finally terminates Your grants, and (b) on an +ongoing basis, if such Contributor fails to notify You of the +non-compliance by some reasonable means prior to 60 days after You have +come back into compliance. Moreover, Your grants from a particular +Contributor are reinstated on an ongoing basis if such Contributor +notifies You of the non-compliance by some reasonable means, this is the +first time You have received notice of non-compliance with this License +from such Contributor, and You become compliant prior to 30 days after +Your receipt of the notice. + +5.2. If You initiate litigation against any entity by asserting a patent +infringement claim (excluding declaratory judgment actions, +counter-claims, and cross-claims) alleging that a Contributor Version +directly or indirectly infringes any patent, then the rights granted to +You by any and all Contributors for the Covered Software under Section +2.1 of this License shall terminate. + +5.3. In the event of termination under Sections 5.1 or 5.2 above, all +end user license agreements (excluding distributors and resellers) which +have been validly granted by You or Your distributors under this License +prior to termination shall survive termination. + +************************************************************************ +* * +* 6. Disclaimer of Warranty * +* ------------------------- * +* * +* Covered Software is provided under this License on an "as is" * +* basis, without warranty of any kind, either expressed, implied, or * +* statutory, including, without limitation, warranties that the * +* Covered Software is free of defects, merchantable, fit for a * +* particular purpose or non-infringing. The entire risk as to the * +* quality and performance of the Covered Software is with You. * +* Should any Covered Software prove defective in any respect, You * +* (not any Contributor) assume the cost of any necessary servicing, * +* repair, or correction. This disclaimer of warranty constitutes an * +* essential part of this License. No use of any Covered Software is * +* authorized under this License except under this disclaimer. * +* * +************************************************************************ + +************************************************************************ +* * +* 7. Limitation of Liability * +* -------------------------- * +* * +* Under no circumstances and under no legal theory, whether tort * +* (including negligence), contract, or otherwise, shall any * +* Contributor, or anyone who distributes Covered Software as * +* permitted above, be liable to You for any direct, indirect, * +* special, incidental, or consequential damages of any character * +* including, without limitation, damages for lost profits, loss of * +* goodwill, work stoppage, computer failure or malfunction, or any * +* and all other commercial damages or losses, even if such party * +* shall have been informed of the possibility of such damages. This * +* limitation of liability shall not apply to liability for death or * +* personal injury resulting from such party's negligence to the * +* extent applicable law prohibits such limitation. Some * +* jurisdictions do not allow the exclusion or limitation of * +* incidental or consequential damages, so this exclusion and * +* limitation may not apply to You. * +* * +************************************************************************ + +8. Litigation +------------- + +Any litigation relating to this License may be brought only in the +courts of a jurisdiction where the defendant maintains its principal +place of business and such litigation shall be governed by laws of that +jurisdiction, without reference to its conflict-of-law provisions. +Nothing in this Section shall prevent a party's ability to bring +cross-claims or counter-claims. + +9. Miscellaneous +---------------- + +This License represents the complete agreement concerning the subject +matter hereof. If any provision of this License is held to be +unenforceable, such provision shall be reformed only to the extent +necessary to make it enforceable. Any law or regulation which provides +that the language of a contract shall be construed against the drafter +shall not be used to construe this License against a Contributor. + +10. Versions of the License +--------------------------- + +10.1. New Versions + +Mozilla Foundation is the license steward. Except as provided in Section +10.3, no one other than the license steward has the right to modify or +publish new versions of this License. Each version will be given a +distinguishing version number. + +10.2. Effect of New Versions + +You may distribute the Covered Software under the terms of the version +of the License under which You originally received the Covered Software, +or under the terms of any subsequent version published by the license +steward. + +10.3. Modified Versions + +If you create software not governed by this License, and you want to +create a new license for such software, you may create and use a +modified version of this License if you rename the license and remove +any references to the name of the license steward (except to note that +such modified license differs from this License). + +10.4. Distributing Source Code Form that is Incompatible With Secondary +Licenses + +If You choose to distribute Source Code Form that is Incompatible With +Secondary Licenses under the terms of this version of the License, the +notice described in Exhibit B of this License must be attached. + +Exhibit A - Source Code Form License Notice +------------------------------------------- + + This Source Code Form is subject to the terms of the Mozilla Public + License, v. 2.0. If a copy of the MPL was not distributed with this + file, You can obtain one at http://mozilla.org/MPL/2.0/. + +If it is not possible or desirable to put the notice in a particular +file, then You may include the notice in a location (such as a LICENSE +file in a relevant directory) where a recipient would be likely to look +for such a notice. + +You may add additional accurate notices of copyright ownership. + +Exhibit B - "Incompatible With Secondary Licenses" Notice +--------------------------------------------------------- + + This Source Code Form is "Incompatible With Secondary Licenses", as + defined by the Mozilla Public License, v. 2.0. diff --git a/symmetry/gemmi/cellred.hpp b/symmetry/gemmi/cellred.hpp new file mode 100644 index 000000000..7174f7766 --- /dev/null +++ b/symmetry/gemmi/cellred.hpp @@ -0,0 +1,406 @@ +// Copyright 2021 Global Phasing Ltd. +// +// Unit cell reductions: Buerger, Niggli, Selling-Delaunay. + +#ifndef GEMMI_CELLRED_HPP_ +#define GEMMI_CELLRED_HPP_ + +#include +#include +#include // for unique_ptr +#include "math.hpp" // for deg +#include "symmetry.hpp" // for Op +#include "unitcell.hpp" // for UnitCell + +namespace gemmi { + +struct SellingVector; + +// GruberVector contains G6 vector (G for Gruber) and cell reduction algorithms. +// Originally, in B. Gruber, Acta Cryst. A29, 433 (1973), the vector was called +// "characteristic" of a lattice/cell. +// Functions that take epsilon as a parameter use it for comparisons, +// as proposed in Grosse-Kunstleve et al, Acta Cryst. (2004) A60, 1. +struct GruberVector { + // a.a b.b c.c 2b.c 2a.c 2a.b + double A, B, C, xi, eta, zeta; // the 1973 paper uses names A B C ξ η ζ + std::unique_ptr change_of_basis; // we use only Op::Rot + + // m - orthogonalization matrix of a primitive cell + explicit GruberVector(const Mat33& m) + : A(m.column_dot(0,0)), + B(m.column_dot(1,1)), + C(m.column_dot(2,2)), + xi(2 * m.column_dot(1,2)), + eta(2 * m.column_dot(0,2)), + zeta(2 * m.column_dot(0,1)) {} + + explicit GruberVector(const std::array& g6) + : A(g6[0]), B(g6[1]), C(g6[2]), xi(g6[3]), eta(g6[4]), zeta(g6[5]) {} + + GruberVector(const UnitCell& u, char centring, bool track_change_of_basis=false) + : GruberVector(u.primitive_orth_matrix(centring)) { + if (track_change_of_basis) + set_change_of_basis(Op{centred_to_primitive(centring), {0,0,0}, 'x'}); + } + + GruberVector(const UnitCell& u, const SpaceGroup* sg, bool track_change_of_basis=false) + : GruberVector(u, sg ? sg->centring_type() : 'P', track_change_of_basis) {} + + void set_change_of_basis(const Op& op) { change_of_basis.reset(new Op(op)); } + + std::array parameters() const { return {A, B, C, xi, eta, zeta}; } + std::array cell_parameters() const { + // inverse of UnitCell::g6() + double a = std::sqrt(A); + double b = std::sqrt(B); + double c = std::sqrt(C); + return {a, b, c, + deg(std::acos(xi/(2*b*c))), + deg(std::acos(eta/(2*a*c))), + deg(std::acos(zeta/(2*a*b)))}; + } + UnitCell get_cell() const { return UnitCell(cell_parameters()); } + + SellingVector selling() const; + + bool is_normalized() const { + // eq(3) from Gruber 1973 + return A <= B && B <= C && + (A != B || std::abs(xi) <= std::abs(eta)) && + (B != C || std::abs(eta) <= std::abs(zeta)) && + (xi > 0) == (eta > 0) && (xi > 0) == (zeta > 0); + } + + bool is_buerger(double epsilon=1e-9) const { + return is_normalized() && + // eq (4) from Gruber 1973 + std::abs(xi) <= B + epsilon && + std::abs(eta) <= A + epsilon && + std::abs(zeta) <= A + epsilon; + } + + // Algorithm N from Gruber (1973). + // Returns branch taken in N3. + void normalize(double eps=1e-9) { + auto step_N1 = [&]() { + if (A - B > eps || (A - B >= -eps && std::abs(xi) > std::abs(eta) + eps)) { // N1 + std::swap(A, B); + std::swap(xi, eta); + if (change_of_basis) + swap_columns_and_negate(0, 1); + } + }; + step_N1(); + if (B - C > eps || (B - C >= -eps && std::abs(eta) > std::abs(zeta) + eps)) { // N2 + std::swap(B, C); + std::swap(eta, zeta); + if (change_of_basis) + swap_columns_and_negate(1, 2); + // To make it faster, instead of "go to the point N1" we repeat N1 once + // (which is equivalent - three swaps are sufficient to reorder ABC). + step_N1(); + } + // N3 + // xi * eta * zeta > 0 <=> positive count is 1 or 3 and no zeros + int pos_count = (xi > eps) + (eta > eps) + (zeta > eps); + int nonneg_count = (xi >= -eps) + (eta >= -eps) + (zeta >= -eps); + double sgn = (pos_count == nonneg_count && pos_count % 2 == 1) ? 1 : -1; + if (change_of_basis) { + if (sgn * xi < -eps) negate_column(0); + if (sgn * eta < -eps) negate_column(1); + if (sgn * zeta < -eps) negate_column(2); + if (pos_count != nonneg_count && pos_count % 2 == 1) + negate_column(std::fabs(zeta) <= eps ? 2 : + std::fabs(eta) <= eps ? 1 : 0); + } + xi = std::copysign(xi, sgn); + eta = std::copysign(eta, sgn); + zeta = std::copysign(zeta, sgn); + } + + // Algorithm B from Gruber (1973). + // Returns true if no change was needed. + bool buerger_step() { + if (std::abs(xi) > B) { // B2 + double j = std::floor(0.5*xi/B + 0.5); + C += j * (j*B - xi); + xi -= 2 * j * B; + eta -= j * zeta; + } else if (std::abs(eta) > A) { // B3 + double j = std::floor(0.5*eta/A + 0.5); + C += j * (j*A - eta); + xi -= j * zeta; + eta -= 2 * j * A; + } else if (std::abs(zeta) > A) { // B4 + double j = std::floor(0.5*zeta/A + 0.5); + B += j * (j*A - zeta); + xi -= j * eta; + zeta -= 2 * j * A; + } else if (xi + eta + zeta + A + B < 0) { // B5 + double j = std::floor(0.5 * (xi + eta) / (A + B + zeta) + 0.5); + C += j * (j * (A + B + zeta) - (xi + eta)); + xi -= j * (2*B + zeta); + eta -= j * (2*A + zeta); + } else { + return true; + } + return false; + } + + // Returns number of iterations. + int buerger_reduce() { + int n = 0; + double prev_sum = -1; + int stall_count = 0; + for (;;) { + normalize(); + // In rare cases numerical errors push the algorithm into infinite loop, + // as described in Grosse-Kunstleve et al, Acta Cryst. (2004) A60, 1. + // Ad-hoc solution: stop if a+b+c is stalled for 5 iterations. + if (++n > 8) { // don't waste time during the first few iterations + double sum = std::sqrt(A) + std::sqrt(B) + std::sqrt(C); + if (std::abs(sum - prev_sum) < sum * 1e-6) { + if (++stall_count == 5) + break; + } else { + stall_count = 0; + } + prev_sum = sum; + } + if (buerger_step()) + break; + } + return n; + } + + // To be called after normalize() or is_normalized(). + // Returns true if it already was Niggli cell. + // Algorithm from Krivy & Gruber, Acta Cryst. (1976) A32, 297. + bool niggli_step(double epsilon=1e-9) { + if (std::abs(xi) > B + epsilon || // step 5. from Krivy-Gruber (1976) + (xi >= B - epsilon && 2 * eta < zeta - epsilon) || + (xi <= -(B - epsilon) && zeta < -epsilon)) { + double sign_xi = xi >= 0 ? 1 : -1; + C += B - xi * sign_xi; + eta -= zeta * sign_xi; + xi -= 2 * B * sign_xi; + if (change_of_basis) + add_column(1, 2, -int(sign_xi)); + } else if (std::abs(eta) > A + epsilon || // step 6. + (eta >= A - epsilon && 2 * xi < zeta - epsilon) || + (eta <= -(A - epsilon) && zeta < -epsilon)) { + double sign_eta = eta >= 0 ? 1 : -1; + C += A - eta * sign_eta; + xi -= zeta * sign_eta; + eta -= 2 * A * sign_eta; + if (change_of_basis) + add_column(0, 2, -int(sign_eta)); + } else if (std::abs(zeta) > A + epsilon || // step 7. + (zeta >= A - epsilon && 2 * xi < eta - epsilon) || + (zeta <= -(A - epsilon) && eta < -epsilon)) { + double sign_zeta = zeta >= 0 ? 1 : -1; + B += A - zeta * sign_zeta; + xi -= eta * sign_zeta; + zeta -= 2 * A * sign_zeta; + if (change_of_basis) + add_column(0, 1, -int(sign_zeta)); + } else if (xi + eta + zeta + A + B < -epsilon || // step 8. + (xi + eta + zeta + A + B <= epsilon && 2 * (A + eta) + zeta > epsilon)) { + C += A + B + xi + eta + zeta; + xi += 2 * B + zeta; + eta += 2 * A + zeta; + if (change_of_basis) { + add_column(0, 2, 1); + add_column(1, 2, 1); + } + } else { + return true; + } + return false; + } + + // Returns number of iterations. + int niggli_reduce(double epsilon=1e-9, int iteration_limit=100) { + int n = 0; + for (;;) { + normalize(epsilon); + if (++n == iteration_limit || niggli_step(epsilon)) + break; + } + return n; + } + + bool is_niggli(double epsilon=1e-9) const { + return is_normalized() && GruberVector(parameters()).niggli_step(epsilon); + } + +private: + void swap_columns_and_negate(int i, int j) { + for (auto& r : change_of_basis->rot) + std::swap(r[i], r[j]); + for (auto& r : change_of_basis->rot) + for (auto& v : r) + v = -v; + } + void negate_column(int i) { + for (auto& r : change_of_basis->rot) + r[i] = -r[i]; + } + void add_column(int pos, int dest, int sign) { + for (auto& r : change_of_basis->rot) + r[dest] += sign * r[pos]; + } +}; + + +// Selling-Delaunay reduction. Based on: +// - chapter "Delaunay reduction and standardization" in +// International Tables for Crystallography vol. A (2016), sec. 3.1.2.3. +// https://onlinelibrary.wiley.com/iucr/itc/Ac/ch3o1v0001/ +// - Patterson & Love (1957), Acta Cryst. 10, 111, +// "Remarks on the Delaunay reduction", doi:10.1107/s0365110x57000328 +// - Andrews et al (2019), Acta Cryst. A75, 115, +// "Selling reduction versus Niggli reduction for crystallographic lattices". +struct SellingVector { + // b.c a.c a.b a.d b.d c.d + std::array s; + + explicit SellingVector(const std::array& s_) : s(s_) {} + + explicit SellingVector(const Mat33& orth) { + Vec3 b[4]; + for (int i = 0; i < 3; ++i) + b[i] = orth.column_copy(i); + b[3]= -b[0] - b[1] - b[2]; + s[0] = b[1].dot(b[2]); + s[1] = b[0].dot(b[2]); + s[2] = b[0].dot(b[1]); + s[3] = b[0].dot(b[3]); + s[4] = b[1].dot(b[3]); + s[5] = b[2].dot(b[3]); + } + + SellingVector(const UnitCell& u, char centring) + : SellingVector(u.primitive_orth_matrix(centring)) {} + SellingVector(const UnitCell& u, const SpaceGroup* sg) + : SellingVector(u, sg ? sg->centring_type() : 'P') {} + + // The reduction minimizes the sum b_i^2 which is equal to -2 sum s_i. + double sum_b_squared() const { + return -2 * (s[0] + s[1] + s[2] + s[3] + s[4] + s[5]); + } + + bool is_reduced(double eps=1e-9) const { + return std::all_of(s.begin(), s.end(), [eps](double x) { return x <= eps; }); + } + + bool reduce_step(double eps=1e-9) { + //printf(" s = %g %g %g %g %g %g sum=%g\n", + // s[0], s[1], s[2], s[3], s[4], s[5], sum_b_squared()); + const int table[6][5] = { + // When negating s[n] we need to apply operations from table[n]: + // 2 x add, subtract, 2 x swap&add + {2, 4, 3, 1, 5}, // 0 + {2, 3, 4, 0, 5}, // 1 + {1, 3, 5, 0, 4}, // 2 + {1, 2, 0, 4, 5}, // 3 + {0, 2, 1, 3, 5}, // 4 + {0, 1, 2, 3, 4}, // 5 + }; + + double max_s = eps; + int max_s_pos = -1; + for (int i = 0; i < 6; ++i) + if (s[i] > max_s) { + max_s = s[i]; + max_s_pos = i; + } + if (max_s_pos < 0) + return false; + const int (&indices)[5] = table[max_s_pos]; + s[max_s_pos] = -max_s; + s[indices[0]] += max_s; + s[indices[1]] += max_s; + s[indices[2]] -= max_s; + std::swap(s[indices[3]], s[indices[4]]); + s[indices[3]] += max_s; + s[indices[4]] += max_s; + //printf(" s[%d]=%g sum: %g\n", max_s_pos, max_s, sum_b_squared()); + return true; + } + + // Returns number of iterations. + int reduce(double eps=1e-9, int iteration_limit=100) { + int n = 0; + while (++n != iteration_limit) + if (!reduce_step(eps)) + break; + return n; + } + + std::array g6_parameters() const { + return {-s[1]-s[2]-s[3], -s[0]-s[2]-s[4], -s[0]-s[1]-s[5], 2*s[0], 2*s[1], 2*s[2]}; + } + + GruberVector gruber() const { return GruberVector(g6_parameters()); } + + // Swap values to make a <= b <= c <= d + void sort(double eps=1e-9) { + double abcd_sq_neg[4] = { + // -a^2, -b^2, -c^2, -d^2 (negated - to be sorted in descending order) + s[1]+s[2]+s[3], s[0]+s[2]+s[4], s[0]+s[1]+s[5], s[3]+s[4]+s[5] + }; + // First, make sure that d >= a,b,c (therefore -d^2 <= -a^2,...). + int min_idx = 3; + for (int i = 0; i < 3; ++i) + if (abcd_sq_neg[i] < abcd_sq_neg[min_idx] - eps) + min_idx = i; + switch (min_idx) { + case 0: // a <-> d + std::swap(s[1], s[5]); + std::swap(s[2], s[4]); + break; + case 1: // b <-> d + std::swap(s[0], s[5]); + std::swap(s[2], s[3]); + break; + case 2: // c <-> d + std::swap(s[0], s[4]); + std::swap(s[1], s[3]); + break; + } + // we could stop here and not care about the order of a,b,c. + std::swap(abcd_sq_neg[min_idx], abcd_sq_neg[3]); + if (abcd_sq_neg[0] < abcd_sq_neg[1] - eps) { // a <-> b + std::swap(s[0], s[1]); + std::swap(s[3], s[4]); + std::swap(abcd_sq_neg[0], abcd_sq_neg[1]); + } + if (abcd_sq_neg[1] < abcd_sq_neg[2] - eps) { // b <-> c + std::swap(s[1], s[2]); + std::swap(s[4], s[5]); + std::swap(abcd_sq_neg[1], abcd_sq_neg[2]); + } + if (abcd_sq_neg[0] < abcd_sq_neg[1] - eps) { // a <-> b + std::swap(s[0], s[1]); + std::swap(s[3], s[4]); + //std::swap(abcd_sq_neg[0], abcd_sq_neg[1]); + } + } + + std::array cell_parameters() const { + return gruber().cell_parameters(); + } + UnitCell get_cell() const { return UnitCell(cell_parameters()); } +}; + +inline SellingVector GruberVector::selling() const { + double s0 = 0.5 * xi; + double s1 = 0.5 * eta; + double s2 = 0.5 * zeta; + return SellingVector({s0, s1, s2, -A - s1 - s2, -B - s0 - s2, -C - s0 - s1}); +} + +} // namespace gemmi +#endif diff --git a/symmetry/gemmi/fail.hpp b/symmetry/gemmi/fail.hpp new file mode 100644 index 000000000..10596385b --- /dev/null +++ b/symmetry/gemmi/fail.hpp @@ -0,0 +1,93 @@ +// Copyright 2017 Global Phasing Ltd. +// +// fail(), unreachable() and __declspec/__attribute__ macros + +#ifndef GEMMI_FAIL_HPP_ +#define GEMMI_FAIL_HPP_ + +#include // for errno +#include // for runtime_error +#include // for system_error +#include +#include // for forward + +#ifdef __INTEL_COMPILER +// warning #2196: routine is both "inline" and "noinline" +# pragma warning disable 2196 +#endif +#if defined(__GNUG__) && !defined(__clang__) +# pragma GCC diagnostic push +# pragma GCC diagnostic ignored "-Wattributes" +#endif + +#if defined(__GNUC__) || defined(__clang__) +# define GEMMI_COLD __attribute__((cold)) +#elif defined(_MSC_VER) +# define GEMMI_COLD __declspec(noinline) +#else +# define GEMMI_COLD __attribute__((noinline)) +#endif + +#if __cplusplus >= 202002L || _MSVC_LANG >= 202002L +# define GEMMI_LIKELY(x) (x) [[likely]] +# define GEMMI_UNLIKELY(x) (x) [[unlikely]] +#elif defined(__GNUC__) || defined(__clang__) +# define GEMMI_LIKELY(x) (__builtin_expect(!!(x), 1)) +# define GEMMI_UNLIKELY(x) (__builtin_expect(!!(x), 0)) +#else +# define GEMMI_LIKELY(x) (x) +# define GEMMI_UNLIKELY(x) (x) +#endif + +#if defined(_WIN32) +# if defined(GEMMI_SHARED) +# if defined(GEMMI_BUILD) +# define GEMMI_DLL __declspec(dllexport) +# else +# define GEMMI_DLL __declspec(dllimport) +# endif // GEMMI_BUILD +# else +# define GEMMI_DLL +# endif // GEMMI_SHARED +#else +# define GEMMI_DLL __attribute__((visibility("default"))) +#endif + +namespace gemmi { + +[[noreturn]] +inline void fail(const std::string& msg) { throw std::runtime_error(msg); } + +template [[noreturn]] +void fail(std::string&& str, T&& arg1, Args&&... args) { + str += arg1; + fail(std::move(str), std::forward(args)...); +} + +[[noreturn]] +inline GEMMI_COLD void fail(const char* msg) { throw std::runtime_error(msg); } + +[[noreturn]] +inline GEMMI_COLD void sys_fail(const std::string& msg) { + throw std::system_error(errno, std::system_category(), msg); +} +[[noreturn]] +inline GEMMI_COLD void sys_fail(const char* msg) { + throw std::system_error(errno, std::system_category(), msg); +} + +// unreachable() is used to silence GCC -Wreturn-type and hint the compiler +[[noreturn]] inline void unreachable() { +#if defined(__GNUC__) || defined(__clang__) + __builtin_unreachable(); +#elif defined(_MSC_VER) + __assume(0); +#endif +} + +#if defined(__GNUG__) && !defined(__clang__) +# pragma GCC diagnostic pop +#endif + +} // namespace gemmi +#endif diff --git a/symmetry/gemmi/math.hpp b/symmetry/gemmi/math.hpp new file mode 100644 index 000000000..87c10516d --- /dev/null +++ b/symmetry/gemmi/math.hpp @@ -0,0 +1,458 @@ +// Copyright 2018 Global Phasing Ltd. +// +// Math utilities. 3D linear algebra. + +#ifndef GEMMI_MATH_HPP_ +#define GEMMI_MATH_HPP_ + +#include // for fabs, cos, sqrt, round +#include // for min +#include +#include // for out_of_range +#include // for enable_if, is_integral + +namespace gemmi { + +constexpr double pi() { return 3.1415926535897932384626433832795029; } + +// The value used in converting between energy[eV] and wavelength[Angstrom]. +// $ units -d15 'h * c / eV / angstrom' +constexpr double hc() { return 12398.4197386209; } + +// The Bohr radius (a0) in Angstroms. +constexpr double bohrradius() { return 0.529177210903; } + +// for Mott-Bethe factor +constexpr double mott_bethe_const() { return 1. / (2 * pi() * pi() * bohrradius()); } + +// Used in conversion of ADPs (atomic displacement parameters). +constexpr double u_to_b() { return 8 * pi() * pi(); } + +constexpr double deg(double angle) { return 180.0 / pi() * angle; } +constexpr double rad(double angle) { return pi() / 180.0 * angle; } + +constexpr float sq(float x) { return x * x; } +constexpr double sq(double x) { return x * x; } + +inline double log_cosh(double x) { + // cosh(x) would overflow for x > 710.5, so we calculate: + // ln(cosh(x)) = ln(e^x + e^-x) - ln(2) = ln(e^x * (1 + e^-2x)) - ln(2) + x = std::abs(x); + return x - std::log(2) + std::log1p(std::exp(-2 * x)); +} + +inline int iround(double d) { return static_cast(std::round(d)); } + +inline double angle_abs_diff(double a, double b, double full=360.0) { + double d = std::fabs(a - b); + if (d > full) + d -= std::floor(d / full) * full; + return std::min(d, full - d); +} + +// similar to C++17 std::clamp() +template constexpr T clamp(T v, T lo, T hi) { + return std::min(std::max(v, lo), hi); +} + +template +struct Vec3_ { + Real x, y, z; + + Vec3_() : x(0), y(0), z(0) {} + Vec3_(Real x_, Real y_, Real z_) : x(x_), y(y_), z(z_) {} + explicit Vec3_(std::array h) : x(h[0]), y(h[1]), z(h[2]) {} + + Real& at(int i) { + switch (i) { + case 0: return x; + case 1: return y; + case 2: return z; + default: throw std::out_of_range("Vec3 index must be 0, 1 or 2."); + } + } + Real at(int i) const { return const_cast(this)->at(i); } + + Vec3_ operator-() const { return {-x, -y, -z}; } + Vec3_ operator-(const Vec3_& o) const { return {x-o.x, y-o.y, z-o.z}; } + Vec3_ operator+(const Vec3_& o) const { return {x+o.x, y+o.y, z+o.z}; } + Vec3_ operator*(Real d) const { return {x*d, y*d, z*d}; } + Vec3_ operator/(Real d) const { return *this * (1.0/d); } + Vec3_& operator-=(const Vec3_& o) { *this = *this - o; return *this; } + Vec3_& operator+=(const Vec3_& o) { *this = *this + o; return *this; } + Vec3_& operator*=(Real d) { *this = *this * d; return *this; } + Vec3_& operator/=(Real d) { return operator*=(1.0/d); } + + Vec3_ negated() const { return {-x, -y, -z}; } + Real dot(const Vec3_& o) const { return x*o.x + y*o.y + z*o.z; } + Vec3_ cross(const Vec3_& o) const { + return {y*o.z - z*o.y, z*o.x - x*o.z, x*o.y - y*o.x}; + } + Real length_sq() const { return x * x + y * y + z * z; } + Real length() const { return std::sqrt(length_sq()); } + Vec3_ changed_magnitude(Real m) const { return operator*(m / length()); } + Vec3_ normalized() const { return changed_magnitude(1.0); } + Real dist_sq(const Vec3_& o) const { return (*this - o).length_sq(); } + Real dist(const Vec3_& o) const { return std::sqrt(dist_sq(o)); } + Real cos_angle(const Vec3_& o) const { + return dot(o) / std::sqrt(length_sq() * o.length_sq()); + } + Real angle(const Vec3_& o) const { + return std::acos(clamp(cos_angle(o), -1., 1.)); + } + bool approx(const Vec3_& o, Real epsilon) const { + return std::fabs(x - o.x) <= epsilon && + std::fabs(y - o.y) <= epsilon && + std::fabs(z - o.z) <= epsilon; + } + bool has_nan() const { + return std::isnan(x) || std::isnan(y) || std::isnan(z); + } +}; + +using Vec3 = Vec3_; +using Vec3f = Vec3_; + +inline Vec3 operator*(double d, const Vec3& v) { return v * d; } + +/// Rodrigues' rotation formula: rotate vector v about given axis of rotation +/// (which must be a unit vector) by given angle (in radians). +inline Vec3 rotate_about_axis(const Vec3& v, const Vec3& axis, double theta) { + double sin_theta = std::sin(theta); + double cos_theta = std::cos(theta); + return v * cos_theta + axis.cross(v) * sin_theta + + axis * (axis.dot(v) * (1 - cos_theta)); +} + +struct Mat33 { + double a[3][3] = { {1.,0.,0.}, {0.,1.,0.}, {0.,0.,1.} }; + + // make it accessible with ".a" + typedef double row_t[3]; + const row_t& operator[](int i) const { return a[i]; } + row_t& operator[](int i) { return a[i]; } + + Mat33() = default; + explicit Mat33(double d) : a{{d, d, d}, {d, d, d}, {d, d, d}} {} + Mat33(double a1, double a2, double a3, double b1, double b2, double b3, + double c1, double c2, double c3) + : a{{a1, a2, a3}, {b1, b2, b3}, {c1, c2, c3}} {} + + static Mat33 from_columns(const Vec3& c1, const Vec3& c2, const Vec3& c3) { + return Mat33(c1.x, c2.x, c3.x, c1.y, c2.y, c3.y, c1.z, c2.z, c3.z); + } + + Vec3 row_copy(int i) const { + if (i < 0 || i > 2) + throw std::out_of_range("Mat33 row index must be 0, 1 or 2."); + return Vec3(a[i][0], a[i][1], a[i][2]); + } + + Vec3 column_copy(int i) const { + if (i < 0 || i > 2) + throw std::out_of_range("Mat33 column index must be 0, 1 or 2."); + return Vec3(a[0][i], a[1][i], a[2][i]); + } + + Mat33 operator+(const Mat33& b) const { + return Mat33(a[0][0] + b[0][0], a[0][1] + b[0][1], a[0][2] + b[0][2], + a[1][0] + b[1][0], a[1][1] + b[1][1], a[1][2] + b[1][2], + a[2][0] + b[2][0], a[2][1] + b[2][1], a[2][2] + b[2][2]); + } + Mat33 operator-(const Mat33& b) const { + return Mat33(a[0][0] - b[0][0], a[0][1] - b[0][1], a[0][2] - b[0][2], + a[1][0] - b[1][0], a[1][1] - b[1][1], a[1][2] - b[1][2], + a[2][0] - b[2][0], a[2][1] - b[2][1], a[2][2] - b[2][2]); + } + + Vec3 multiply(const Vec3& p) const { + return {a[0][0] * p.x + a[0][1] * p.y + a[0][2] * p.z, + a[1][0] * p.x + a[1][1] * p.y + a[1][2] * p.z, + a[2][0] * p.x + a[2][1] * p.y + a[2][2] * p.z}; + } + Vec3 left_multiply(const Vec3& p) const { + return {a[0][0] * p.x + a[1][0] * p.y + a[2][0] * p.z, + a[0][1] * p.x + a[1][1] * p.y + a[2][1] * p.z, + a[0][2] * p.x + a[1][2] * p.y + a[2][2] * p.z}; + } + // p has elements from the main diagonal of a 3x3 diagonal matrix + Mat33 multiply_by_diagonal(const Vec3& p) const { + return Mat33(a[0][0] * p.x, a[0][1] * p.y, a[0][2] * p.z, + a[1][0] * p.x, a[1][1] * p.y, a[1][2] * p.z, + a[2][0] * p.x, a[2][1] * p.y, a[2][2] * p.z); + } + Mat33 multiply(const Mat33& b) const { + Mat33 r; + for (int i = 0; i != 3; ++i) + for (int j = 0; j != 3; ++j) + r[i][j] = a[i][0] * b[0][j] + a[i][1] * b[1][j] + a[i][2] * b[2][j]; + return r; + } + Mat33 transpose() const { + return Mat33(a[0][0], a[1][0], a[2][0], + a[0][1], a[1][1], a[2][1], + a[0][2], a[1][2], a[2][2]); + } + double trace() const { return a[0][0] + a[1][1] + a[2][2]; } + + bool approx(const Mat33& other, double epsilon) const { + for (int i = 0; i < 3; ++i) + for (int j = 0; j < 3; ++j) + if (std::fabs(a[i][j] - other.a[i][j]) > epsilon) + return false; + return true; + } + bool has_nan() const { + for (int i = 0; i < 3; ++i) + for (int j = 0; j < 3; ++j) + if (std::isnan(a[i][j])) + return true; + return false; + } + + double determinant() const { + return a[0][0] * (a[1][1]*a[2][2] - a[2][1]*a[1][2]) + + a[0][1] * (a[1][2]*a[2][0] - a[2][2]*a[1][0]) + + a[0][2] * (a[1][0]*a[2][1] - a[2][0]*a[1][1]); + } + Mat33 inverse() const { + Mat33 inv; + double inv_det = 1.0 / determinant(); + inv[0][0] = inv_det * (a[1][1] * a[2][2] - a[2][1] * a[1][2]); + inv[0][1] = inv_det * (a[0][2] * a[2][1] - a[0][1] * a[2][2]); + inv[0][2] = inv_det * (a[0][1] * a[1][2] - a[0][2] * a[1][1]); + inv[1][0] = inv_det * (a[1][2] * a[2][0] - a[1][0] * a[2][2]); + inv[1][1] = inv_det * (a[0][0] * a[2][2] - a[0][2] * a[2][0]); + inv[1][2] = inv_det * (a[1][0] * a[0][2] - a[0][0] * a[1][2]); + inv[2][0] = inv_det * (a[1][0] * a[2][1] - a[2][0] * a[1][1]); + inv[2][1] = inv_det * (a[2][0] * a[0][1] - a[0][0] * a[2][1]); + inv[2][2] = inv_det * (a[0][0] * a[1][1] - a[1][0] * a[0][1]); + return inv; + } + bool is_identity() const { + return a[0][0] == 1 && a[0][1] == 0 && a[0][2] == 0 && + a[1][0] == 0 && a[1][1] == 1 && a[1][2] == 0 && + a[2][0] == 0 && a[2][1] == 0 && a[2][2] == 1; + } + + double column_dot(int i, int j) const { + return a[0][i] * a[0][j] + a[1][i] * a[1][j] + a[2][i] * a[2][j]; + } + + bool is_upper_triangular() const { + return a[1][0] == 0 && a[2][0] == 0 && a[2][1] == 0; + } +}; + +struct UpperTriangularMat33 { + double a11 = 0, a12 = 0, a13 = 0; + double a22 = 0, a23 = 0; + double a33 = 0; + UpperTriangularMat33() = default; + UpperTriangularMat33& operator=(const Mat33& m) { + if (m.is_upper_triangular()) { + a11 = m[0][0]; + a12 = m[0][1]; + a13 = m[0][2]; + a22 = m[1][1]; + a23 = m[1][2]; + a33 = m[2][2]; + } else { + a11 = a12 = a13 = a22 = a23 = a33 = NAN; + } + return *this; + } + Vec3 multiply(const Vec3& p) const { + return {a11 * p.x + a12 * p.y + a13 * p.z, + a22 * p.y + a23 * p.z, + a33 * p.z}; + } +}; + +// Symmetric matrix 3x3. Used primarily for an ADP tensor. +template struct SMat33 { + T u11, u22, u33, u12, u13, u23; + + // The PDB ANISOU record has the above order, but in a different context + // (such as metric tensor) the order of Voigt notation may be preferred. + std::array elements_pdb() const { return {{u11, u22, u33, u12, u13, u23}}; } + std::array elements_voigt() const { return {{u11, u22, u33, u23, u13, u12}}; } + + Mat33 as_mat33() const { + return Mat33(u11, u12, u13, u12, u22, u23, u13, u23, u33); + } + + // the arguments i and j must be in [0,2], i.e. 0, 1 or 2. + T& unchecked_ref(int i, int j) { + T* ptrs[9] = {&u11, &u12, &u13, &u12, &u22, &u23, &u13, &u23, &u33}; + return *ptrs[3 * i + j]; + } + + T trace() const { return u11 + u22 + u33; } + bool nonzero() const { return trace() != 0; } + + bool all_zero() const { + return u11 == 0 && u22 == 0 && u33 == 0 && u12 == 0 && u13 == 0 && u23 == 0; + } + + void scale(T s) const { + u11 *= s; u22 *= s; u33 *= s; u12 *= s; u13 *= s; u23 *= s; + } + + template + SMat33 scaled(Real s) const { + return SMat33{u11*s, u22*s, u33*s, u12*s, u13*s, u23*s}; + } + + // returns U + kI + SMat33 added_kI(T k) const { + return {u11+k, u22+k, u33+k, u12, u13, u23}; + } + + // returns squared norm r^T U r where U is this matrix and vector r is arg + template + auto r_u_r(const Vec3_& r) const -> decltype(r.x+u11) { + return r.x * r.x * u11 + r.y * r.y * u22 + r.z * r.z * u33 + + 2 * (r.x * r.y * u12 + r.x * r.z * u13 + r.y * r.z * u23); + } + double r_u_r(const std::array& h) const { + // it's faster to first convert ints to doubles (Vec3) + return r_u_r(Vec3(h)); + } + + Vec3 multiply(const Vec3& p) const { + return {u11 * p.x + u12 * p.y + u13 * p.z, + u12 * p.x + u22 * p.y + u23 * p.z, + u13 * p.x + u23 * p.y + u33 * p.z}; + } + + SMat33 operator-(const SMat33& o) const { + return {u11-o.u11, u22-o.u22, u33-o.u33, u12-o.u12, u13-o.u13, u23-o.u23}; + } + SMat33 operator+(const SMat33& o) const { + return {u11+o.u11, u22+o.u22, u33+o.u33, u12+o.u12, u13+o.u13, u23+o.u23}; + } + + // return M U M^T + template + SMat33 transformed_by(const Mat33& m) const { + // slightly faster than m.multiply(as_mat33()).multiply(m.transpose()); + auto elem = [&](int i, int j) { + return static_cast( + m[i][0] * (m[j][0] * u11 + m[j][1] * u12 + m[j][2] * u13) + + m[i][1] * (m[j][0] * u12 + m[j][1] * u22 + m[j][2] * u23) + + m[i][2] * (m[j][0] * u13 + m[j][1] * u23 + m[j][2] * u33)); + }; + return SMat33{elem(0, 0), elem(1, 1), elem(2, 2), + elem(0, 1), elem(0, 2), elem(1, 2)}; + } + + T determinant() const { + return u11 * (u22*u33 - u23*u23) + + u12 * (u23*u13 - u33*u12) + + u13 * (u12*u23 - u13*u22); + } + + SMat33 inverse_(T det) const { + SMat33 inv; + T inv_det = 1.0f / det; + inv.u11 = inv_det * (u22 * u33 - u23 * u23); + inv.u22 = inv_det * (u11 * u33 - u13 * u13); + inv.u33 = inv_det * (u11 * u22 - u12 * u12); + inv.u12 = inv_det * (u13 * u23 - u12 * u33); + inv.u13 = inv_det * (u12 * u23 - u13 * u22); + inv.u23 = inv_det * (u12 * u13 - u11 * u23); + return inv; + } + SMat33 inverse() const { + return inverse_(determinant()); + } + + /// Based on https://en.wikipedia.org/wiki/Eigenvalue_algorithm + /// To calculate both eigenvalues and eigenvectors use eig3.hpp + std::array calculate_eigenvalues() const { + double p1 = u12*u12 + u13*u13 + u23*u23; + if (p1 == 0) + return {{u11, u22, u33}}; + double q = (1./3.) * trace(); + SMat33 b{u11 - q, u22 - q, u33 - q, u12, u13, u23}; + double p2 = sq(b.u11) + sq(b.u22) + sq(b.u33) + 2 * p1; + double p = std::sqrt((1./6.) * p2); + double r = b.determinant() / ((1./3.) * p2 * p); + double phi = 0; + if (r <= -1) + phi = (1./3.) * pi(); + else if (r < 1) + phi = (1./3.) * std::acos(r); + double eig1 = q + 2 * p * std::cos(phi); + double eig3 = q + 2 * p * std::cos(phi + 2./3.*pi()); + return {{eig1, 3 * q - eig1 - eig3, eig3}}; + } +}; + +struct Transform { + Mat33 mat; + Vec3 vec; + + Transform inverse() const { + Mat33 minv = mat.inverse(); + return {minv, minv.multiply(vec).negated()}; + } + + Vec3 apply(const Vec3& x) const { return mat.multiply(x) + vec; } + + Transform combine(const Transform& b) const { + return {mat.multiply(b.mat), vec + mat.multiply(b.vec)}; + } + + bool is_identity() const { + return mat.is_identity() && vec.x == 0. && vec.y == 0. && vec.z == 0.; + } + void set_identity() { mat = Mat33(); vec = Vec3(); } + + bool has_nan() const { + return mat.has_nan() || vec.has_nan(); + } + + bool approx(const Transform& o, double epsilon) const { + return mat.approx(o.mat, epsilon) && vec.approx(o.vec, epsilon); + } +}; + +template +struct Box { + Pos minimum = Pos(INFINITY, INFINITY, INFINITY); + Pos maximum = Pos(-INFINITY, -INFINITY, -INFINITY); + void extend(const Pos& p) { + if (p.x < minimum.x) minimum.x = p.x; + if (p.y < minimum.y) minimum.y = p.y; + if (p.z < minimum.z) minimum.z = p.z; + if (p.x > maximum.x) maximum.x = p.x; + if (p.y > maximum.y) maximum.y = p.y; + if (p.z > maximum.z) maximum.z = p.z; + } + Pos get_size() const { return maximum - minimum; } + void add_margins(const Pos& p) { minimum -= p; maximum += p; } + void add_margin(double m) { add_margins(Pos(m, m, m)); } +}; + +// internally used functions +namespace impl { +// MSVC is missing isnan(IntegralType), so we define is_nan as a replacement +template +typename std::enable_if::value, bool>::type +is_nan(T) { return false; } +template +typename std::enable_if::value, bool>::type +is_nan(T a) { return std::isnan(a); } + +template +typename std::enable_if::value, bool>::type +is_same(T a, T b) { return a == b; } +template +typename std::enable_if::value, bool>::type +is_same(T a, T b) { return std::isnan(b) ? std::isnan(a) : a == b; } +} // namespace impl + +} // namespace gemmi +#endif diff --git a/symmetry/gemmi/symmetry.hpp b/symmetry/gemmi/symmetry.hpp new file mode 100644 index 000000000..203324fdd --- /dev/null +++ b/symmetry/gemmi/symmetry.hpp @@ -0,0 +1,1044 @@ +// Copyright 2017-2019 Global Phasing Ltd. +// +// Crystallographic Symmetry. Space Groups. Coordinate Triplets. +// +// If this is all that you need from Gemmi you can just copy this file, +// src/symmetry.cpp fail.hpp and LICENSE.txt to your project. + +#ifndef GEMMI_SYMMETRY_HPP_ +#define GEMMI_SYMMETRY_HPP_ + +#include // for strtol, abs +#include +#include // for sort, remove +#include // for hash +#include // for invalid_argument +#include +#include // for tie +#include + +#include "fail.hpp" // for fail, unreachable + +namespace gemmi { + +// OP + +// Op is a symmetry operation, or a change-of-basic transformation, +// or a different operation of similar kind. +// Both "rotation" matrix and translation vector are fractional, with DEN +// used as the denominator. +struct GEMMI_DLL Op { + static constexpr int DEN = 24; // 24 to handle 1/8 in change-of-basis + typedef std::array, 3> Rot; + typedef std::array Tran; + + Rot rot; + Tran tran; + char notation = ' '; + + bool is_hkl() const { return notation == 'h'; } + + Op as_hkl() const { + return is_hkl() ? *this : Op{rot, {0,0,0}, 'h'}; + } + Op as_xyz() const { + return is_hkl() ? Op{rot, {0,0,0}, 'x'} : *this; + } + + std::string triplet(char style=' ') const; + + Op inverse() const; + + Op::Tran wrapped_tran() const { + Op::Tran t = tran; + for (int i = 0; i != 3; ++i) { + if (t[i] >= DEN) // elements need to be in [0,DEN) + t[i] %= DEN; + else if (t[i] < 0) + t[i] = ((t[i] + 1) % DEN) + DEN - 1; + } + return t; + } + + // If the translation points outside of the unit cell, wrap it. + Op& wrap() { + tran = wrapped_tran(); + return *this; + } + + Op& translate(const Tran& a) { + for (int i = 0; i != 3; ++i) + tran[i] += a[i]; + return *this; + } + + Op translated(const Tran& a) const { return Op(*this).translate(a); } + + Op add_centering(const Tran& a) const { return translated(a).wrap(); } + + Rot negated_rot() const { + return {{{-rot[0][0], -rot[0][1], -rot[0][2]}, + {-rot[1][0], -rot[1][1], -rot[1][2]}, + {-rot[2][0], -rot[2][1], -rot[2][2]}}}; + } + + static Rot transpose(const Rot& rot) { + return {{{rot[0][0], rot[1][0], rot[2][0]}, + {rot[0][1], rot[1][1], rot[2][1]}, + {rot[0][2], rot[1][2], rot[2][2]}}}; + } + Rot transposed_rot() const { return transpose(rot); } + + // DEN^3 for rotation, -DEN^3 for rotoinversion + int det_rot() const { + return rot[0][0] * (rot[1][1] * rot[2][2] - rot[1][2] * rot[2][1]) + - rot[0][1] * (rot[1][0] * rot[2][2] - rot[1][2] * rot[2][0]) + + rot[0][2] * (rot[1][0] * rot[2][1] - rot[1][1] * rot[2][0]); + } + + // Rotation-part type based on Table 1 in RWGK, Acta Cryst. A55, 383 (1999) + int rot_type() const { + int det = det_rot(); + int tr_den = rot[0][0] + rot[1][1] + rot[2][2]; + int tr = tr_den / DEN; + const int table[] = {0, 0, 2, 3, 4, 6, 1}; + if (std::abs(det) == DEN * DEN * DEN && tr * DEN == tr_den && std::abs(tr) <= 3) + return det > 0 ? table[3 + tr] : -table[3 - tr]; + return 0; + } + + Op combine(const Op& b) const { + if (is_hkl() != b.is_hkl()) + fail("can't combine real- and reciprocal-space Op"); + Op r; + for (int i = 0; i != 3; ++i) { + r.tran[i] = tran[i] * Op::DEN; + for (int j = 0; j != 3; ++j) { + r.rot[i][j] = (rot[i][0] * b.rot[0][j] + + rot[i][1] * b.rot[1][j] + + rot[i][2] * b.rot[2][j]) / Op::DEN; + r.tran[i] += rot[i][j] * b.tran[j]; + } + r.tran[i] /= Op::DEN; + } + r.notation = notation; + return r; + } + + std::array apply_to_xyz(const std::array& xyz) const { + if (is_hkl()) + fail("can't apply reciprocal-space Op to xyz"); + std::array out; + for (int i = 0; i != 3; ++i) + out[i] = (rot[i][0] * xyz[0] + rot[i][1] * xyz[1] + rot[i][2] * xyz[2] + + tran[i]) / Op::DEN; + return out; + } + + // Miller is defined in the same way in namespace gemmi in unitcell.hpp + using Miller = std::array; + + Miller apply_to_hkl_without_division(const Miller& hkl) const { + Miller r; + for (int i = 0; i != 3; ++i) + r[i] = (rot[0][i] * hkl[0] + rot[1][i] * hkl[1] + rot[2][i] * hkl[2]); + return r; + } + static Miller divide_hkl_by_DEN(const Miller& hkl) { + return {{ hkl[0] / DEN, hkl[1] / DEN, hkl[2] / DEN }}; + } + Miller apply_to_hkl(const Miller& hkl) const { + return divide_hkl_by_DEN(apply_to_hkl_without_division(hkl)); + } + + double phase_shift(const Miller& hkl) const { + constexpr double mult = -2 * 3.1415926535897932384626433832795 / Op::DEN; + return mult * (hkl[0] * tran[0] + hkl[1] * tran[1] + hkl[2] * tran[2]); + } + + std::array, 4> int_seitz() const { + std::array, 4> t; + for (int i = 0; i < 3; ++i) + t[i] = { rot[i][0], rot[i][1], rot[i][2], tran[i] }; + t[3] = { 0, 0, 0, 1 }; + return t; + } + + std::array, 4> float_seitz() const { + std::array, 4> t; + double m = 1.0 / Op::DEN; + for (int i = 0; i < 3; ++i) + t[i] = { m * rot[i][0], m * rot[i][1], m * rot[i][2], m * tran[i] }; + t[3] = { 0., 0., 0., 1. }; + return t; + } + + static constexpr Op identity() { + return {{{{DEN,0,0}, {0,DEN,0}, {0,0,DEN}}}, {0,0,0}, ' '}; + } + static constexpr Op::Rot inversion_rot() { + return {{{-DEN,0,0}, {0,-DEN,0}, {0,0,-DEN}}}; + } + bool operator<(const Op& rhs) const { + return std::tie(rot, tran) < std::tie(rhs.rot, rhs.tran); + } +}; + +inline bool operator==(const Op& a, const Op& b) { + return a.rot == b.rot && a.tran == b.tran; +} +inline bool operator!=(const Op& a, const Op& b) { return !(a == b); } + +inline Op operator*(const Op& a, const Op& b) { return a.combine(b).wrap(); } +inline Op& operator*=(Op& a, const Op& b) { a = a * b; return a; } + +inline Op Op::inverse() const { + int detr = det_rot(); + if (detr == 0) + fail("cannot invert matrix: " + Op{rot, {0,0,0}, notation}.triplet()); + int d2 = Op::DEN * Op::DEN; + Op inv; + inv.rot[0][0] = d2 * (rot[1][1] * rot[2][2] - rot[2][1] * rot[1][2]) / detr; + inv.rot[0][1] = d2 * (rot[0][2] * rot[2][1] - rot[0][1] * rot[2][2]) / detr; + inv.rot[0][2] = d2 * (rot[0][1] * rot[1][2] - rot[0][2] * rot[1][1]) / detr; + inv.rot[1][0] = d2 * (rot[1][2] * rot[2][0] - rot[1][0] * rot[2][2]) / detr; + inv.rot[1][1] = d2 * (rot[0][0] * rot[2][2] - rot[0][2] * rot[2][0]) / detr; + inv.rot[1][2] = d2 * (rot[1][0] * rot[0][2] - rot[0][0] * rot[1][2]) / detr; + inv.rot[2][0] = d2 * (rot[1][0] * rot[2][1] - rot[2][0] * rot[1][1]) / detr; + inv.rot[2][1] = d2 * (rot[2][0] * rot[0][1] - rot[0][0] * rot[2][1]) / detr; + inv.rot[2][2] = d2 * (rot[0][0] * rot[1][1] - rot[1][0] * rot[0][1]) / detr; + for (int i = 0; i != 3; ++i) + inv.tran[i] = (-tran[0] * inv.rot[i][0] + -tran[1] * inv.rot[i][1] + -tran[2] * inv.rot[i][2]) / Op::DEN; + inv.notation = notation; + return inv; +} + +// inverse of Op::float_seitz() +GEMMI_DLL Op seitz_to_op(const std::array, 4>& t); + +// helper function for use in AsuBrick::str() +GEMMI_DLL void append_op_fraction(std::string& s, int w); + +// TRIPLET -> OP +GEMMI_DLL std::array parse_triplet_part(const std::string& s, char& notation, + double* decimal_fract=nullptr); +GEMMI_DLL Op parse_triplet(const std::string& s, char notation=' '); + +// GROUPS OF OPERATIONS + +// corresponds to Table A1.4.2.2 in ITfC vol.B (edition 2010) +inline std::vector centring_vectors(char centring_type) { + constexpr int h = Op::DEN / 2; + constexpr int t = Op::DEN / 3; + constexpr int d = 2 * t; + // note: find_centering() depends on the order of operations in vector + switch (centring_type & ~0x20) { + case 'P': return {{0, 0, 0}}; + case 'A': return {{0, 0, 0}, {0, h, h}}; + case 'B': return {{0, 0, 0}, {h, 0, h}}; + case 'C': return {{0, 0, 0}, {h, h, 0}}; + case 'I': return {{0, 0, 0}, {h, h, h}}; + case 'R': return {{0, 0, 0}, {d, t, t}, {t, d, d}}; + // hall_symbols.html has no H, ITfC 2010 has no S and T + case 'H': return {{0, 0, 0}, {d, t, 0}, {t, d, 0}}; + case 'S': return {{0, 0, 0}, {t, t, d}, {d, t, d}}; + case 'T': return {{0, 0, 0}, {t, d, t}, {d, t, d}}; + case 'F': return {{0, 0, 0}, {0, h, h}, {h, 0, h}, {h, h, 0}}; + default: fail("not a centring type: ", centring_type); + } +} + + +struct GroupOps { + std::vector sym_ops; + std::vector cen_ops; + + int order() const { return static_cast(sym_ops.size()*cen_ops.size()); } + + void add_missing_elements(); + void add_missing_elements_part2(const std::vector& gen, + size_t max_size, bool ignore_bad_gen); + + bool add_inversion() { + size_t init_size = sym_ops.size(); + sym_ops.reserve(2 * init_size); + for (const Op& op : sym_ops) { + Op::Rot neg = op.negated_rot(); + if (find_by_rotation(neg)) { + sym_ops.resize(init_size); + return false; + } + sym_ops.push_back({neg, op.tran, op.notation}); + } + return true; + } + + char find_centering() const { + if (cen_ops.size() == 1 && cen_ops[0] == Op::Tran{0, 0, 0}) + return 'P'; + std::vector trans = cen_ops; + std::sort(trans.begin(), trans.end()); + for (char c : {'A', 'B', 'C', 'I', 'F', 'R', 'H', 'S', 'T'}) { + std::vector c_vectors = centring_vectors(c); + if (c == 'R' || c == 'H') // these two are returned not sorted + std::swap(c_vectors[1], c_vectors[2]); + if (trans == c_vectors) + return c; + } + return 0; + } + + Op* find_by_rotation(const Op::Rot& r) { + for (Op& op : sym_ops) + if (op.rot == r) + return &op; + return nullptr; + } + + const Op* find_by_rotation(const Op::Rot& r) const { + return const_cast(this)->find_by_rotation(r); + } + + bool is_centrosymmetric() const { + return find_by_rotation(Op::inversion_rot()) != nullptr; + } + + bool is_reflection_centric(const Op::Miller& hkl) const { + Op::Miller mhkl = {{-Op::DEN * hkl[0], -Op::DEN * hkl[1], -Op::DEN * hkl[2]}}; + for (const Op& op : sym_ops) + if (op.apply_to_hkl_without_division(hkl) == mhkl) + return true; + return false; + } + + int epsilon_factor_without_centering(const Op::Miller& hkl) const { + Op::Miller denh = {{Op::DEN * hkl[0], Op::DEN * hkl[1], Op::DEN * hkl[2]}}; + int epsilon = 0; + for (const Op& op : sym_ops) + if (op.apply_to_hkl_without_division(hkl) == denh) + ++epsilon; + return epsilon; + } + int epsilon_factor(const Op::Miller& hkl) const { + return epsilon_factor_without_centering(hkl) * (int) cen_ops.size(); + } + + static bool has_phase_shift(const Op::Tran& c, const Op::Miller& hkl) { + return (hkl[0] * c[0] + hkl[1] * c[1] + hkl[2] * c[2]) % Op::DEN != 0; + } + + bool is_systematically_absent(const Op::Miller& hkl) const { + for (auto i = cen_ops.begin() + 1; i != cen_ops.end(); ++i) + if (has_phase_shift(*i, hkl)) + return true; + Op::Miller denh = {{Op::DEN * hkl[0], Op::DEN * hkl[1], Op::DEN * hkl[2]}}; + for (auto op = sym_ops.begin() + 1; op != sym_ops.end(); ++op) + if (op->apply_to_hkl_without_division(hkl) == denh) { + for (const Op::Tran& c : cen_ops) + if (has_phase_shift({{op->tran[0] + c[0], + op->tran[1] + c[1], + op->tran[2] + c[2]}}, hkl)) + return true; + } + return false; + } + + void change_basis_impl(const Op& cob, const Op& inv) { + if (sym_ops.empty() || cen_ops.empty()) + return; + + // Apply change-of-basis to sym_ops. + // Ignore the first item in sym_ops -- it's identity. + for (auto op = sym_ops.begin() + 1; op != sym_ops.end(); ++op) + *op = cob.combine(*op).combine(inv).wrap(); + + // The number of centering vectors may be different. + // As an ad-hoc method (not proved to be robust) add lattice points + // from a super-cell. + int idet = inv.det_rot() / (Op::DEN * Op::DEN * Op::DEN); + if (idet > 1) { + std::vector new_cen_ops; + new_cen_ops.reserve(cen_ops.size() * idet * idet * idet); + for (int i = 0; i < idet; ++i) + for (int j = 0; j < idet; ++j) + for (int k = 0; k < idet; ++k) + for (Op::Tran& cen : cen_ops) + new_cen_ops.push_back({i * Op::DEN + cen[0], + j * Op::DEN + cen[1], + k * Op::DEN + cen[2]}); + cen_ops.swap(new_cen_ops); + } + + // Apply change-of-basis to centering vectors + Op cvec = Op::identity(); + for (auto tr = cen_ops.begin() + 1; tr != cen_ops.end(); ++tr) { + cvec.tran = *tr; + *tr = cob.combine(cvec).combine(inv).wrap().tran; + } + + // Remove redundant centering vectors. + for (int i = static_cast(cen_ops.size()) - 1; i > 0; --i) + for (int j = i - 1; j >= 0; --j) + if (cen_ops[i] == cen_ops[j]) { + cen_ops.erase(cen_ops.begin() + i); + break; + } + } + + void change_basis_forward(const Op& cob) { change_basis_impl(cob, cob.inverse()); } + void change_basis_backward(const Op& inv) { change_basis_impl(inv.inverse(), inv); } + + std::vector all_ops_sorted() const { + std::vector ops; + ops.reserve(sym_ops.size() * cen_ops.size()); + for (const Op& so : sym_ops) + for (const Op::Tran& co : cen_ops) + ops.push_back(so.add_centering(co)); + std::sort(ops.begin(), ops.end()); + return ops; + } + + Op get_op(int n) const { + int n_cen = n / (int) sym_ops.size(); + int n_sym = n % (int) sym_ops.size(); + return sym_ops.at(n_sym).add_centering(cen_ops.at(n_cen)); + } + + bool is_same_as(const GroupOps& other) const { + if (cen_ops.size() != other.cen_ops.size() || + sym_ops.size() != other.sym_ops.size()) + return false; + return all_ops_sorted() == other.all_ops_sorted(); + } + + bool has_same_centring(const GroupOps& other) const { + if (cen_ops.size() != other.cen_ops.size()) + return false; + if (std::is_sorted(cen_ops.begin(), cen_ops.end()) && + std::is_sorted(other.cen_ops.begin(), other.cen_ops.end())) + return cen_ops == other.cen_ops; + std::vector v1 = cen_ops; + std::vector v2 = other.cen_ops; + std::sort(v1.begin(), v1.end()); + std::sort(v2.begin(), v2.end()); + return v1 == v2; + } + + bool has_same_rotations(const GroupOps& other) const { + if (sym_ops.size() != other.sym_ops.size()) + return false; + auto sorted_rotations = [](const GroupOps& g) { + std::vector r(g.sym_ops.size()); + for (size_t i = 0; i != r.size(); ++i) + r[i] = g.sym_ops[i].rot; + std::sort(r.begin(), r.end()); + return r; + }; + return sorted_rotations(*this) == sorted_rotations(other); + } + + // minimal multiplicity for real-space grid in each direction + // examples: 1,2,1 for P21, 1,1,6 for P61 + std::array find_grid_factors() const { + const int T = Op::DEN; + int r[3] = {T, T, T}; + for (Op op : *this) + for (int i = 0; i != 3; ++i) + if (op.tran[i] != 0 && op.tran[i] < r[i]) + r[i] = op.tran[i]; + return {T / r[0], T / r[1], T / r[2]}; + } + + bool are_directions_symmetry_related(int u, int v) const { + for (const Op& op : sym_ops) + if (op.rot[u][v] != 0) + return true; + return false; + } + + // remove translation part of sym_ops + GroupOps derive_symmorphic() const { + GroupOps r(*this); + for (Op& op : r.sym_ops) + op.tran[0] = op.tran[1] = op.tran[2] = 0; + return r; + } + + struct Iter { + const GroupOps& gops; + int n_sym, n_cen; + void operator++() { + if (++n_sym == (int) gops.sym_ops.size()) { + ++n_cen; + n_sym = 0; + } + } + Op operator*() const { + return gops.sym_ops.at(n_sym).translated(gops.cen_ops.at(n_cen)).wrap(); + } + bool operator==(const Iter& other) const { + return n_sym == other.n_sym && n_cen == other.n_cen; + } + bool operator!=(const Iter& other) const { return !(*this == other); } + }; + + Iter begin() const { return {*this, 0, 0}; } + Iter end() const { return {*this, 0, (int) cen_ops.size()}; } +}; + +inline void GroupOps::add_missing_elements() { + // We always keep identity as sym_ops[0]. + if (sym_ops.empty() || sym_ops[0] != Op::identity()) + fail("oops"); + if (sym_ops.size() == 1) + return; + constexpr size_t max_size = 1024; + // Below we assume that all centring vectors are already known (in cen_ops) + // so when checking for a new element we compare only the 3x3 matrix. + // Dimino's algorithm. https://physics.stackexchange.com/a/351400/95713 + std::vector gen(sym_ops.begin() + 1, sym_ops.end()); + sym_ops.resize(2); + const Op::Rot idrot = Op::identity().rot; + for (Op g = sym_ops[1] * sym_ops[1]; g.rot != idrot; g *= sym_ops[1]) { + sym_ops.push_back(g); + if (sym_ops.size() > max_size) + fail("Too many elements in the group - bad generators"); + } + // the rest is in separate function b/c it's reused in twin.hpp + add_missing_elements_part2(gen, max_size, false); +} + +inline void GroupOps::add_missing_elements_part2(const std::vector& gen, + size_t max_size, bool ignore_bad_gen) { + for (size_t i = 1; i < gen.size(); ++i) { + std::vector coset_repr(1, Op::identity()); + size_t init_size = sym_ops.size(); + for (;;) { + size_t len = coset_repr.size(); + for (size_t j = 0; j != len; ++j) { + for (size_t n = 0; n != i + 1; ++n) { + Op sg = gen[n] * coset_repr[j]; + if (find_by_rotation(sg.rot) == nullptr) { + sym_ops.push_back(sg); + for (size_t k = 1; k != init_size; ++k) + sym_ops.push_back(sg * sym_ops[k]); + coset_repr.push_back(sg); + } + } + } + if (len == coset_repr.size()) + break; + if (sym_ops.size() > max_size) { + if (!ignore_bad_gen) + fail("Too many elements in the group - bad generators"); + // ignore this generator and continue with the next one + sym_ops.resize(init_size); + break; + } + } + } +} + +// Create GroupOps from Ops by separating centering vectors +inline GroupOps split_centering_vectors(const std::vector& ops) { + const Op identity = Op::identity(); + GroupOps go; + go.sym_ops.push_back(identity); + for (const Op& op : ops) + if (Op* old_op = go.find_by_rotation(op.rot)) { + Op::Tran tran = op.wrapped_tran(); + if (op.rot == identity.rot) // pure shift + go.cen_ops.push_back(tran); + if (tran == identity.tran) // or rather |tran| < |old_op->tran| ? + old_op->tran = op.tran; + } else { + go.sym_ops.push_back(op); + } + return go; +} + +GEMMI_DLL GroupOps generators_from_hall(const char* hall); + +inline GroupOps symops_from_hall(const char* hall) { + GroupOps ops = generators_from_hall(hall); + ops.add_missing_elements(); + return ops; +} + +// CRYSTAL SYSTEMS, POINT GROUPS AND LAUE CLASSES + +enum class CrystalSystem : unsigned char { + Triclinic=0, Monoclinic, Orthorhombic, Tetragonal, Trigonal, Hexagonal, Cubic +}; + +inline const char* crystal_system_str(CrystalSystem system) { + static const char* names[7] = { + "triclinic", "monoclinic", "orthorhombic", "tetragonal", + "trigonal", "hexagonal", "cubic" + }; + return names[static_cast(system)]; +} + +enum class PointGroup : unsigned char { + C1=0, Ci, C2, Cs, C2h, D2, C2v, D2h, C4, S4, C4h, D4, C4v, D2d, D4h, C3, + C3i, D3, C3v, D3d, C6, C3h, C6h, D6, C6v, D3h, D6h, T, Th, O, Td, Oh +}; + +inline const char* point_group_hm(PointGroup pg) { + static const char hm_pointgroup_names[32][6] = { + "1", "-1", "2", "m", "2/m", "222", "mm2", "mmm", + "4", "-4", "4/m", "422", "4mm", "-42m", "4/mmm", "3", + "-3", "32", "3m", "-3m", "6", "-6", "6/m", "622", + "6mm", "-62m", "6/mmm", "23", "m-3", "432", "-43m", "m-3m", + }; + return hm_pointgroup_names[static_cast(pg)]; +} + +// http://reference.iucr.org/dictionary/Laue_class +enum class Laue : unsigned char { + L1=0, L2m, Lmmm, L4m, L4mmm, L3, L3m, L6m, L6mmm, Lm3, Lm3m +}; + +inline Laue pointgroup_to_laue(PointGroup pg) { + static const Laue laue[32] = { + Laue::L1, Laue::L1, + Laue::L2m, Laue::L2m, Laue::L2m, + Laue::Lmmm, Laue::Lmmm, Laue::Lmmm, + Laue::L4m, Laue::L4m, Laue::L4m, + Laue::L4mmm, Laue::L4mmm, Laue::L4mmm, Laue::L4mmm, + Laue::L3, Laue::L3, + Laue::L3m, Laue::L3m, Laue::L3m, + Laue::L6m, Laue::L6m, Laue::L6m, + Laue::L6mmm, Laue::L6mmm, Laue::L6mmm, Laue::L6mmm, + Laue::Lm3, Laue::Lm3, + Laue::Lm3m, Laue::Lm3m, Laue::Lm3m, + }; + return laue[static_cast(pg)]; +} + +// return centrosymmetric pointgroup from the Laue class +inline PointGroup laue_to_pointgroup(Laue laue) { + static const PointGroup pg[11] = { + PointGroup::Ci, PointGroup::C2h, PointGroup::D2h, PointGroup::C4h, + PointGroup::D4h, PointGroup::C3i, PointGroup::D3d, PointGroup::C6h, + PointGroup::D6h, PointGroup::Th, PointGroup::Oh + }; + return pg[static_cast(laue)]; +} + +inline const char* laue_class_str(Laue laue) { + return point_group_hm(laue_to_pointgroup(laue)); +} + +inline CrystalSystem crystal_system(Laue laue) { + static const CrystalSystem crystal_systems[11] = { + CrystalSystem::Triclinic, + CrystalSystem::Monoclinic, + CrystalSystem::Orthorhombic, + CrystalSystem::Tetragonal, CrystalSystem::Tetragonal, + CrystalSystem::Trigonal, CrystalSystem::Trigonal, + CrystalSystem::Hexagonal, CrystalSystem::Hexagonal, + CrystalSystem::Cubic, CrystalSystem::Cubic + }; + return crystal_systems[static_cast(laue)]; +} + +inline CrystalSystem crystal_system(PointGroup pg) { + return crystal_system(pointgroup_to_laue(pg)); +} + +inline unsigned char point_group_index_and_category(int space_group_number) { + // 0x20=Sohncke, 0x40=enantiomorphic, 0x80=symmorphic + enum : unsigned char { S=0x20, E=(0x20|0x40), Y=0x80, Z=(0x20|0x80) }; + static const unsigned char indices[230] = { + 0|Z, 1|Y, 2|Z, 2|S, 2|Z, 3|Y, 3, 3|Y, 3, 4|Y, // 1-10 + 4, 4|Y, 4, 4, 4, 5|Z, 5|S, 5|S, 5|S, 5|S, // 11-20 + 5|Z, 5|Z, 5|Z, 5|S, 6|Y, 6, 6, 6, 6, 6, // 21-30 + 6, 6, 6, 6, 6|Y, 6, 6, 6|Y, 6, 6, // 31-40 + 6, 6|Y, 6, 6|Y, 6, 6, 7|Y, 7, 7, 7, // 41-50 + 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, // 51-60 + 7, 7, 7, 7, 7|Y, 7, 7, 7, 7|Y, 7, // 61-70 + 7|Y, 7, 7, 7, 8|Z, 8|E, 8|S, 8|E, 8|Z, 8|S, // 71-80 + 9|Y, 9|Y, 10|Y, 10, 10, 10, 10|Y, 10, 11|Z, 11|S, // 81-90 + 11|E, 11|E, 11|S, 11|S, 11|E, 11|E, 11|Z, 11|S, 12|Y, 12, // 91-100 + 12, 12, 12, 12, 12, 12, 12|Y, 12, 12, 12, // 101-110 + 13|Y, 13, 13, 13, 13|Y, 13, 13, 13, 13|Y, 13, // 111-120 + 13|Y, 13, 14|Y, 14, 14, 14, 14, 14, 14, 14, // 121-130 + 14, 14, 14, 14, 14, 14, 14, 14, 14|Y, 14, // 131-140 + 14, 14, 15|Z, 15|E, 15|E, 15|Z, 16|Y, 16|Y, 17|Z, 17|Z, // 141-150 + 17|E, 17|E, 17|E, 17|E, 17|Z, 18|Y, 18|Y, 18, 18, 18|Y, // 151-160 + 18, 19|Y, 19, 19|Y, 19, 19|Y, 19, 20|Z, 20|E, 20|E, // 161-170 + 20|E, 20|E, 20|S, 21|Y, 22|Y, 22, 23|Z, 23|E, 23|E, 23|E, // 171-180 + 23|E, 23|S, 24|Y, 24, 24, 24, 25|Y, 25, 25|Y, 25, // 181-190 + 26|Y, 26, 26, 26, 27|Z, 27|Z, 27|Z, 27|S, 27|S, 28|Y, // 191-200 + 28, 28|Y, 28, 28|Y, 28, 28, 29|Z, 29|S, 29|Z, 29|S, // 201-210 + 29|Z, 29|E, 29|E, 29|S, 30|Y, 30|Y, 30|Y, 30, 30, 30, // 211-220 + 31|Y, 31, 31, 31, 31|Y, 31, 31, 31, 31|Y, 31 // 221-230 + }; + return indices[space_group_number-1]; +} + +inline PointGroup point_group(int space_group_number) { + auto n = point_group_index_and_category(space_group_number); + return static_cast(n & 0x1f); +} + +// true for 65 Sohncke (non-enantiogenic) space groups +inline bool is_sohncke(int space_group_number) { + return (point_group_index_and_category(space_group_number) & 0x20) != 0; +} + +// true for 22 space groups (11 enantiomorphic pairs) +inline bool is_enantiomorphic(int space_group_number) { + return (point_group_index_and_category(space_group_number) & 0x40) != 0; +} + +// true for 73 space groups +inline bool is_symmorphic(int space_group_number) { + return (point_group_index_and_category(space_group_number) & 0x80) != 0; +} + +/// Inversion center of the Euclidean normalizer that is not at the origin of +/// reference settings. Returns (0,0,0) if absent. Based on tables in ch. 3.5 +/// of ITA (2016) doi:10.1107/97809553602060000933 (column "Inversion through +/// a centre at"). +inline Op::Tran nonzero_inversion_center(int space_group_number) { + constexpr int D = Op::DEN; + switch (space_group_number) { + case 43: return {D/8, D/8, 0}; + case 80: return {D/4, 0, 0}; + case 98: return {D/4, 0, D/8}; + case 109: return {D/4, 0, 0}; + case 110: return {D/4, 0, 0}; + case 122: return {D/4, 0, D/8}; + case 210: return {D/8, D/8, D/8}; + default: return {0, 0, 0}; + } +} + +GEMMI_DLL const char* get_basisop(int basisop_idx); + + +// Returns a change-of-basis operator for centred -> primitive transformation. +// The same operator as inverse of z2p_op in sgtbx. +inline Op::Rot centred_to_primitive(char centring_type) { + constexpr int D = Op::DEN; + constexpr int H = Op::DEN / 2; + constexpr int T = Op::DEN / 3; + switch (centring_type) { + case 'P': return {{{D,0,0}, {0,D,0}, {0,0,D}}}; + case 'A': return {{{-D,0,0}, {0,-H,H}, {0,H,H}}}; + case 'B': return {{{-H,0,H}, {0,-D,0}, {H,0,H}}}; + case 'C': return {{{H,H,0}, {H,-H,0}, {0,0,-D}}}; + case 'I': return {{{-H,H,H}, {H,-H,H}, {H,H,-H}}}; + case 'R': return {{{2*T,-T,-T}, {T,T,-2*T}, {T,T,T}}}; + case 'H': return {{{2*T,-T,0}, {T,T,0}, {0,0,D}}}; // not used normally + case 'F': return {{{0,H,H}, {H,0,H}, {H,H,0}}}; + default: fail("not a centring type: ", centring_type); + } +} + + +// LIST OF CRYSTALLOGRAPHIC SPACE GROUPS + +struct SpaceGroup { // typically 44 bytes + int number; + int ccp4; + char hm[11]; // Hermann-Mauguin (international) notation + char ext; + char qualifier[5]; + char hall[15]; + int basisop_idx; + + std::string xhm() const { + std::string ret = hm; + if (ext) { + ret += ':'; + ret += ext; + } + return ret; + } + + char centring_type() const { return ext == 'R' ? 'P' : hm[0]; } + + // (old) CCP4 spacegroup names start with H for hexagonal setting + char ccp4_lattice_type() const { return ext == 'H' ? 'H' : hm[0]; } + + // P 1 2 1 -> P2, but P 1 1 2 -> P112. R 3:H -> H3. + std::string short_name() const { + std::string s(hm); + size_t len = s.size(); + if (len > 6 && s[2] == '1' && s[len - 2] == ' ' && s[len - 1] == '1') + s = s[0] + s.substr(4, len - 4 - 2); + if (ext == 'H') + s[0] = 'H'; + s.erase(std::remove(s.begin(), s.end(), ' '), s.end()); + return s; + } + + // As explained in Phenix newsletter CCN_2011_01.pdf#page=12 + // the PDB uses own, non-standard symbols for rhombohedral space groups. + std::string pdb_name() const { + std::string s; + s += ccp4_lattice_type(); + s += hm+1; + return s; + } + + bool is_sohncke() const { return gemmi::is_sohncke(number); } + bool is_enantiomorphic() const { return gemmi::is_enantiomorphic(number); } + bool is_symmorphic() const { return gemmi::is_symmorphic(number); } + PointGroup point_group() const { return gemmi::point_group(number); } + const char* point_group_hm() const { + return gemmi::point_group_hm(point_group()); + } + Laue laue_class() const { return pointgroup_to_laue(point_group()); } + const char* laue_str() const { return laue_class_str(laue_class()); } + CrystalSystem crystal_system() const { + return gemmi::crystal_system(point_group()); + } + const char* crystal_system_str() const { + return gemmi::crystal_system_str(crystal_system()); + } + bool is_centrosymmetric() const { + return laue_to_pointgroup(laue_class()) == point_group(); + } + + /// returns 'a', 'b' or 'c' for monoclinic SG, '\0' otherwise + char monoclinic_unique_axis() const { + if (crystal_system() == CrystalSystem::Monoclinic) + return qualifier[qualifier[0] == '-' ? 1 : 0]; + return '\0'; + } + + const char* basisop_str() const { return get_basisop(basisop_idx); } + Op basisop() const { return parse_triplet(basisop_str()); } + bool is_reference_setting() const { return basisop_idx == 0; } + + Op centred_to_primitive() const { + return {gemmi::centred_to_primitive(centring_type()), {0,0,0}, 'x'}; + } + + /// Returns change-of-hand operator. Compatible with similar sgtbx function. + Op change_of_hand_op() const { + if (is_centrosymmetric()) + return Op::identity(); + Op::Tran t = nonzero_inversion_center(number); + Op op{Op::inversion_rot(), {2*t[0], 2*t[1], 2*t[2]}, 'x'}; + if (!is_reference_setting()) { + Op b = basisop(); + op = b.combine(op).combine(b.inverse()); + } + return op; + } + + GroupOps operations() const { return symops_from_hall(hall); } +}; + +struct SpaceGroupAltName { + char hm[11]; + char ext; + int pos; +}; + +struct GEMMI_DLL spacegroup_tables { + static const SpaceGroup main[564]; + static const SpaceGroupAltName alt_names[28]; + static const unsigned char ccp4_hkl_asu[230]; +}; + +inline const SpaceGroup* find_spacegroup_by_number(int ccp4) noexcept { + if (ccp4 == 0) + return &spacegroup_tables::main[0]; + for (const SpaceGroup& sg : spacegroup_tables::main) + if (sg.ccp4 == ccp4) + return &sg; + return nullptr; +} + +inline const SpaceGroup& get_spacegroup_by_number(int ccp4) { + const SpaceGroup* sg = find_spacegroup_by_number(ccp4); + if (sg == nullptr) + throw std::invalid_argument("Invalid space-group number: " + + std::to_string(ccp4)); + return *sg; +} + +inline const SpaceGroup& get_spacegroup_reference_setting(int number) { + for (const SpaceGroup& sg : spacegroup_tables::main) + if (sg.number == number && sg.is_reference_setting()) + return sg; + throw std::invalid_argument("Invalid space-group number: " + + std::to_string(number)); +} + +/// If angles alpha and gamma are provided, they are used to +/// distinguish hexagonal and rhombohedral settings (e.g. for "R 3"). +/// \param prefer can specify preferred H/R settings and 1/2 origin choice. +/// For example, prefer="2H" means the origin choice 2 and hexagonal +/// settings. The default is "1H". +GEMMI_DLL const SpaceGroup* find_spacegroup_by_name(std::string name, + double alpha=0., double gamma=0., + const char* prefer=nullptr); + +inline const SpaceGroup& get_spacegroup_by_name(const std::string& name) { + const SpaceGroup* sg = find_spacegroup_by_name(name); + if (sg == nullptr) + throw std::invalid_argument("Unknown space-group name: " + name); + return *sg; +} + +inline const SpaceGroup& get_spacegroup_p1() { + return spacegroup_tables::main[0]; +} + +inline const SpaceGroup* find_spacegroup_by_ops(const GroupOps& gops) { + char c = gops.find_centering(); + for (const SpaceGroup& sg : spacegroup_tables::main) + if ((c == sg.hall[0] || c == sg.hall[1]) && + gops.is_same_as(sg.operations())) + return &sg; + return nullptr; +} + +// Reciprocal space asu (asymmetric unit). +// The same 12 choices of ASU as in CCP4 symlib and cctbx. +struct ReciprocalAsu { + int idx; + Op::Rot rot{}; // value-initialized only to avoid -Wmaybe-uninitialized + bool is_ref; + + ReciprocalAsu(const SpaceGroup* sg, bool tnt=false) { + if (sg == nullptr) + fail("Missing space group"); + idx = spacegroup_tables::ccp4_hkl_asu[sg->number - 1]; + if (tnt) { + idx += 10; + is_ref = true; // TNT ASU is given wrt current (not standard) settings + } else { + is_ref = sg->is_reference_setting(); + if (!is_ref) + rot = sg->basisop().rot; + } + } + + bool is_in(const Op::Miller& hkl) const { + if (is_ref) + return is_in_reference_setting(hkl[0], hkl[1], hkl[2]); + Op::Miller r; + for (int i = 0; i != 3; ++i) + r[i] = rot[0][i] * hkl[0] + rot[1][i] * hkl[1] + rot[2][i] * hkl[2]; + return is_in_reference_setting(r[0], r[1], r[2]); + } + + bool is_in_reference_setting(int h, int k, int l) const { + switch (idx) { + // 0-9: CCP4 hkl asu, 10-19: TNT hkl asu + case 0: return l>0 || (l==0 && (h>0 || (h==0 && k>=0))); + case 1: return k>=0 && (l>0 || (l==0 && h>=0)); + case 12: // orthorhombic-D + case 2: return h>=0 && k>=0 && l>=0; + case 3: return l>=0 && ((h>=0 && k>0) || (h==0 && k==0)); + case 14: // tetragonal-D, hexagonal-D + case 4: return h>=k && k>=0 && l>=0; + case 5: return (h>=0 && k>0) || (h==0 && k==0 && l>=0); + case 16: // trigonal-D P312 + case 6: return h>=k && k>=0 && (k>0 || l>=0); + case 17: // trigonal-D P321 + case 7: return h>=k && k>=0 && (h>k || l>=0); + case 8: return h>=0 && ((l>=h && k>h) || (l==h && k==h)); + case 9: return k>=l && l>=h && h>=0; + case 10: return k>0 || (k==0 && (h>0 || (h==0 && l>=0))); // triclinic + case 11: return k>=0 && (h>0 || (h==0 && l>=0)); // monoclinic-B + case 13: return l>=0 && ((k>=0 && h>0) || (h==0 && k==0)); // tetragonal-C, hexagonal-C + case 15: return (k>=0 && h>0) || (h==0 && k==0 && l>=0); // trigonal-C + case 18: return k>=0 && l>=0 && ((h>k && h>l) || (h==k && h>=l)); // cubic-T + case 19: return h>=k && k>=l && l>=0; // cubic-O + } + unreachable(); + } + + const char* condition_str() const { + switch (idx) { + case 0: return "l>0 or (l=0 and (h>0 or (h=0 and k>=0)))"; + case 1: return "k>=0 and (l>0 or (l=0 and h>=0))"; + case 12: + case 2: return "h>=0 and k>=0 and l>=0"; + case 3: return "l>=0 and ((h>=0 and k>0) or (h=0 and k=0))"; + case 14: + case 4: return "h>=k and k>=0 and l>=0"; + case 5: return "(h>=0 and k>0) or (h=0 and k=0 and l>=0)"; + case 16: + case 6: return "h>=k and k>=0 and (k>0 or l>=0)"; + case 17: + case 7: return "h>=k and k>=0 and (h>k or l>=0)"; + case 8: return "h>=0 and ((l>=h and k>h) or (l=h and k=h))"; + case 9: return "k>=l and l>=h and h>=0"; + case 10: return "k>0 or (k==0 and (h>0 or (h=0 and l>=0)))"; + case 11: return "k>=0 and (h>0 or (h=0 and l>=0))"; + case 13: return "l>=0 and ((k>=0 and h>0) or (h=0 and k==0))"; + case 15: return "(k>=0 and h>0) or (h=0 and k==0 and l>=0)"; + case 18: return "k>=0 and l>=0 and ((h>k and h>l) or (h=k and h>=l))"; + case 19: return "h>=k and k>=l and l>=0"; + } + unreachable(); + } + + /// Returns hkl in asu and MTZ ISYM - 2*n-1 for reflections in the positive + /// asu (I+ of a Friedel pair), 2*n for reflections in the negative asu (I-). + std::pair to_asu(const Op::Miller& hkl, const std::vector& sym_ops) const { + int isym = 0; + for (const Op& op : sym_ops) { + ++isym; + Op::Miller new_hkl = op.apply_to_hkl_without_division(hkl); + if (is_in(new_hkl)) + return {Op::divide_hkl_by_DEN(new_hkl), isym}; + ++isym; + Op::Miller negated_new_hkl{{-new_hkl[0], -new_hkl[1], -new_hkl[2]}}; + if (is_in(negated_new_hkl)) + return {Op::divide_hkl_by_DEN(negated_new_hkl), isym}; + } + fail("Oops, maybe inconsistent GroupOps?"); + } + + std::pair to_asu(const Op::Miller& hkl, const GroupOps& gops) const { + return to_asu(hkl, gops.sym_ops); + } + + /// Similar to to_asu(), but the second returned value is sign: true for + or centric + std::pair to_asu_sign(const Op::Miller& hkl, const GroupOps& gops) const { + std::pair neg = {{0,0,0}, true}; + for (const Op& op : gops.sym_ops) { + Op::Miller new_hkl = op.apply_to_hkl_without_division(hkl); + if (is_in(new_hkl)) + return {Op::divide_hkl_by_DEN(new_hkl), true}; + Op::Miller negated_new_hkl{{-new_hkl[0], -new_hkl[1], -new_hkl[2]}}; + if (is_in(negated_new_hkl)) + // don't return it yet, because for centric reflection we prefer (+) + neg = {Op::divide_hkl_by_DEN(negated_new_hkl), false}; + } + if (neg.second) + fail("Oops, maybe inconsistent GroupOps?"); + return neg; + } +}; + +} // namespace gemmi + +namespace std { +template<> struct hash { + size_t operator()(const gemmi::Op& op) const { + size_t h = 0; + for (int i = 0; i != 3; ++i) + for (int j = 0; j != 3; ++j) + h = (h << 2) ^ (op.rot[i][j] + 1); + for (int i = 0; i != 3; ++i) + h = (h << 5) ^ op.tran[i]; + return h; + } +}; +} // namespace std + +#endif diff --git a/symmetry/gemmi/unitcell.hpp b/symmetry/gemmi/unitcell.hpp new file mode 100644 index 000000000..25bb8b467 --- /dev/null +++ b/symmetry/gemmi/unitcell.hpp @@ -0,0 +1,618 @@ +// Copyright 2017 Global Phasing Ltd. +// +// Unit cell. + +#ifndef GEMMI_UNITCELL_HPP_ +#define GEMMI_UNITCELL_HPP_ + +#include +#include // for cos, sin, sqrt, floor, NAN +#include +#include "math.hpp" +#include "fail.hpp" // for fail +#include "symmetry.hpp" // for Op, SpaceGroup + +namespace gemmi { + +inline Mat33 rot_as_mat33(const Op::Rot& rot) { + double mult = 1.0 / Op::DEN; + return Mat33(mult * rot[0][0], mult * rot[0][1], mult * rot[0][2], + mult * rot[1][0], mult * rot[1][1], mult * rot[1][2], + mult * rot[2][0], mult * rot[2][1], mult * rot[2][2]); +} +inline Mat33 rot_as_mat33(const Op& op) { return rot_as_mat33(op.rot); } + + +inline Vec3 tran_as_vec3(const Op& op) { + double mult = 1.0 / Op::DEN; + return Vec3(mult * op.tran[0], mult * op.tran[1], mult * op.tran[2]); +} + +/// Coordinates in Angstroms - orthogonal (Cartesian) coordinates. +struct Position : Vec3 { + using Vec3::Vec3; + Position() = default; + explicit Position(const Vec3& v) : Vec3(v) {} + Position operator-() const { return Position(Vec3::operator-()); } + Position operator-(const Position& o) const { return Position(Vec3::operator-(o)); } + Position operator+(const Position& o) const { return Position(Vec3::operator+(o)); } + Position operator*(double d) const { return Position(Vec3::operator*(d)); } + Position operator/(double d) const { return Position(Vec3::operator/(d)); } + Position& operator-=(const Position& o) { *this = *this - o; return *this; } + Position& operator+=(const Position& o) { *this = *this + o; return *this; } + Position& operator*=(double d) { *this = *this * d; return *this; } + Position& operator/=(double d) { return operator*=(1.0/d); } +}; + +inline Position operator*(double d, const Position& v) { return v * d; } + +/// Fractional coordinates. +struct Fractional : Vec3 { + using Vec3::Vec3; + Fractional() = default; + explicit Fractional(const Vec3& v) : Vec3(v) {} + Fractional operator-(const Fractional& o) const { + return Fractional(Vec3::operator-(o)); + } + Fractional operator+(const Fractional& o) const { + return Fractional(Vec3::operator+(o)); + } + Fractional wrap_to_unit() const { + return {x - std::floor(x), y - std::floor(y), z - std::floor(z)}; + } + Fractional wrap_to_zero() const { + return {x - std::round(x), y - std::round(y), z - std::round(z)}; + } + Fractional round() const { + return {std::round(x), std::round(y), std::round(z)}; + } + void move_toward_zero_by_one() { + if (x > 0.5) x -= 1.0; else if (x < -0.5) x += 1.0; + if (y > 0.5) y -= 1.0; else if (y < -0.5) y += 1.0; + if (z > 0.5) z -= 1.0; else if (z < -0.5) z += 1.0; + } +}; + +enum class Asu : unsigned char { Same, Different, Any }; + +/// Result of find_nearest_image +struct NearestImage { + double dist_sq; + int pbc_shift[3] = { 0, 0, 0 }; + int sym_idx = 0; + + double dist() const { return std::sqrt(dist_sq); } + bool same_asu() const { + return pbc_shift[0] == 0 && pbc_shift[1] == 0 && pbc_shift[2] == 0 && sym_idx == 0; + } + + /// Returns a string such as 1555 or 1_555. + std::string symmetry_code(bool underscore) const { + std::string s = std::to_string(sym_idx + 1); + if (underscore) + s += '_'; + if (unsigned(5 + pbc_shift[0]) <= 9 && + unsigned(5 + pbc_shift[1]) <= 9 && + unsigned(5 + pbc_shift[2]) <= 9) { // normal, quick path + for (int shift : pbc_shift) + s += char('5' + shift); + } else { // problematic, non-standard path + for (int i = 0; i < 3; ++i) { + if (i != 0 && underscore) + s += '_'; + s += std::to_string(5 + pbc_shift[i]); + } + } + return s; + } +}; + + +/// Like Transform, but apply() arg is Fractional (not Vec3 - for type safety). +struct FTransform : Transform { + FTransform() = default; + FTransform(const Transform& t) : Transform(t) {} + Fractional apply(const Fractional& p) const { + return Fractional(Transform::apply(p)); + } +}; + +/// Non-crystallographic symmetry operation (such as in the MTRIXn record) +struct NcsOp { + std::string id; + bool given; + Transform tr; + Position apply(const Position& p) const { return Position(tr.apply(p)); } +}; + +/// A synonym for convenient passing of hkl. +using Miller = std::array; + +struct MillerHash { + std::size_t operator()(const Miller& hkl) const noexcept { + return std::size_t((hkl[0] * 1024 + hkl[1]) * 1024 + hkl[2]); // NOLINT misplaced cast + } +}; + +struct UnitCellParameters { + double a = 1.0, b = 1.0, c = 1.0; + double alpha = 90.0, beta = 90.0, gamma = 90.0; + + UnitCellParameters() = default; + explicit UnitCellParameters(const double (&par)[6]) { + a = par[0]; b = par[1]; c = par[2]; alpha = par[3]; beta = par[4]; gamma = par[5]; + } + explicit UnitCellParameters(const std::array& par) { + a = par[0]; b = par[1]; c = par[2]; alpha = par[3]; beta = par[4]; gamma = par[5]; + } + + bool operator==(const UnitCellParameters& o) const { + return a == o.a && b == o.b && c == o.c && + alpha == o.alpha && beta == o.beta && gamma == o.gamma; + } + bool operator!=(const UnitCellParameters& o) const { return !operator==(o); } + + bool approx(const UnitCellParameters& o, double epsilon) const { + auto eq = [&](double x, double y) { return std::fabs(x - y) < epsilon; }; + return eq(a, o.a) && eq(b, o.b) && eq(c, o.c) && + eq(alpha, o.alpha) && eq(beta, o.beta) && eq(gamma, o.gamma); + } +}; + +/// Unit cell. Contains cell parameters as well as pre-calculated +/// orthogonalization and fractionalization matrices, volume, and more. +/// Contains symmetry operations (incl. NCS) if they were set from outside. +struct UnitCell : UnitCellParameters { + UnitCell() = default; + UnitCell(double a_, double b_, double c_, + double alpha_, double beta_, double gamma_) { + set(a_, b_, c_, alpha_, beta_, gamma_); + } + UnitCell(const std::array& v) { set_from_array(v); } + + Transform orth; + Transform frac; + double volume = 1.0; + /// reciprocal parameters a*, b*, c*, alpha*, beta*, gamma* + double ar = 1.0, br = 1.0, cr = 1.0; + double cos_alphar = 0.0, cos_betar = 0.0, cos_gammar = 0.0; + bool explicit_matrices = false; + short cs_count = 0; // crystallographic symmetries except identity + std::vector images; // symmetry operations + + // Non-crystalline (for example NMR) structures are supposed to use fake + // unit cell 1x1x1, but sometimes they don't. A number of non-crystalline + // entries in the PDB has incorrectly set unit cell or fract. matrix, + // that is why we check both. + bool is_crystal() const { return a != 1.0 && frac.mat[0][0] != 1.0; } + + // compare lengths using relative tolerance rel, angles using tolerance deg + bool is_similar(const UnitCell& o, double rel, double deg) const { + auto siml = [&](double x, double y) { return std::fabs(x - y) < rel * std::max(x, y); }; + auto sima = [&](double x, double y) { return std::fabs(x - y) < deg; }; + return siml(a, o.a) && siml(b, o.b) && siml(c, o.c) && + sima(alpha, o.alpha) && sima(beta, o.beta) && sima(gamma, o.gamma); + } + + void calculate_properties() { + // ensure exact values for right angles + double cos_alpha = alpha == 90. ? 0. : std::cos(rad(alpha)); + double cos_beta = beta == 90. ? 0. : std::cos(rad(beta)); + double cos_gamma = gamma == 90. ? 0. : std::cos(rad(gamma)); + double sin_alpha = alpha == 90. ? 1. : std::sin(rad(alpha)); + double sin_beta = beta == 90. ? 1. : std::sin(rad(beta)); + double sin_gamma = gamma == 90. ? 1. : std::sin(rad(gamma)); + if (sin_alpha == 0 || sin_beta == 0 || sin_gamma == 0) + fail("Impossible angle - N*180deg."); + + // volume - formula from Giacovazzo p.62 + volume = a * b * c * std::sqrt(1 - cos_alpha * cos_alpha + - cos_beta * cos_beta - cos_gamma * cos_gamma + + 2 * cos_alpha * cos_beta * cos_gamma); + + // reciprocal parameters a*, b*, ... (Giacovazzo, p. 64) + ar = b * c * sin_alpha / volume; + br = a * c * sin_beta / volume; + cr = a * b * sin_gamma / volume; + double cos_alphar_sin_beta = (cos_beta * cos_gamma - cos_alpha) / sin_gamma; + cos_alphar = cos_alphar_sin_beta / sin_beta; + //cos_alphar = (cos_beta * cos_gamma - cos_alpha) / (sin_beta * sin_gamma); + cos_betar = (cos_alpha * cos_gamma - cos_beta) / (sin_alpha * sin_gamma); + cos_gammar = (cos_alpha * cos_beta - cos_gamma) / (sin_alpha * sin_beta); + + if (explicit_matrices) + return; + + // The orthogonalization matrix we use is described in ITfC B p.262: + // "An alternative mode of orthogonalization, used by the Protein + // Data Bank and most programs, is to align the a1 axis of the unit + // cell with the Cartesian X_1 axis, and to align the a*_3 axis with the + // Cartesian X_3 axis." + double sin_alphar = std::sqrt(1.0 - cos_alphar * cos_alphar); + orth.mat = {a, b * cos_gamma, c * cos_beta, + 0., b * sin_gamma, -c * cos_alphar_sin_beta, + 0., 0. , c * sin_beta * sin_alphar}; + orth.vec = {0., 0., 0.}; + + double o12 = -cos_gamma / (sin_gamma * a); + double o13 = -(cos_gamma * cos_alphar_sin_beta + cos_beta * sin_gamma) + / (sin_alphar * sin_beta * sin_gamma * a); + double o23 = cos_alphar / (sin_alphar * sin_gamma * b); + frac.mat = {1 / a, o12, o13, + 0., 1 / orth.mat[1][1], o23, + 0., 0., 1 / orth.mat[2][2]}; + frac.vec = {0., 0., 0.}; + } + + double cos_alpha() const { return alpha == 90. ? 0. : std::cos(rad(alpha)); } + + /// B matrix following convention from Busing & Levy (1967), not from cctbx. + /// Cf. https://dials.github.io/documentation/conventions.html + Mat33 calculate_matrix_B() const { + double sin_gammar = std::sqrt(1 - cos_gammar * cos_gammar); + double sin_betar = std::sqrt(1 - cos_betar * cos_betar); + return Mat33(ar, br * cos_gammar, cr * cos_betar, + 0., br * sin_gammar, -cr * sin_betar * cos_alpha(), + 0., 0., 1.0 / c); + } + + /// The equivalent isotropic displacement factor. + /// Based on Fischer & Tillmanns (1988). Acta Cryst. C44, 775-776. + /// The argument is a non-orthogonalized tensor U, + /// i.e. the one from SmallStructure::Site, but not from Atom. + double calculate_u_eq(const SMat33& ani) const { + double aar = a * ar; + double bbr = b * br; + double ccr = c * cr; + // it could be optimized using orth.mat[0][1] and orth.mat[0][2] + double cos_beta = beta == 90. ? 0. : std::cos(rad(beta)); + double cos_gamma = gamma == 90. ? 0. : std::cos(rad(gamma)); + return 1/3. * (sq(aar) * ani.u11 + sq(bbr) * ani.u22 + sq(ccr) * ani.u33 + + 2 * (aar * bbr * cos_gamma * ani.u12 + + aar * ccr * cos_beta * ani.u13 + + bbr * ccr * cos_alpha() * ani.u23)); + } + + void set_matrices_from_fract(const Transform& f) { + // mmCIF _atom_sites.fract_transf_* and PDB SCALEn records usually contain + // fewer significant digits than the unit cell parameters, and sometimes are + // just wrong. Use them only if we seem to have non-standard crystal frame. + if (f.mat.approx(frac.mat, 1e-4) && f.vec.approx(frac.vec, 1e-6)) + return; + // The SCALE record is sometimes incorrect. Here we only catch cases + // when CRYST1 is set as for non-crystal and SCALE is very suspicious. + if (frac.mat[0][0] == 1.0 && (f.mat[0][0] == 0.0 || f.mat[0][0] > 1.0)) + return; + frac = f; + orth = f.inverse(); + explicit_matrices = true; + } + + void set(double a_, double b_, double c_, + double alpha_, double beta_, double gamma_) { + if (gamma_ == 0.0) // ignore empty/partial CRYST1 (example: 3iyp) + return; + a = a_; + b = b_; + c = c_; + alpha = alpha_; + beta = beta_; + gamma = gamma_; + calculate_properties(); + } + + void set_from_parameters(const UnitCellParameters& p) { + set(p.a, p.b, p.c, p.alpha, p.beta, p.gamma); + } + + void set_from_array(const std::array& v) { set(v[0], v[1], v[2], v[3], v[4], v[5]); } + + void set_from_vectors(const Vec3& va, const Vec3& vb, const Vec3& vc) { + set(va.length(), vb.length(), vc.length(), + deg(vb.angle(vc)), deg(vc.angle(va)), deg(va.angle(vb))); + } + + UnitCell changed_basis_backward(const Op& op, bool set_images) { + Mat33 mat = orth.mat.multiply(rot_as_mat33(op)); + UnitCell new_cell; + new_cell.set_from_vectors(mat.column_copy(0), + mat.column_copy(1), + mat.column_copy(2)); + if (set_images && !images.empty()) { + new_cell.images.reserve(images.size()); + Transform tr{rot_as_mat33(op), tran_as_vec3(op)}; + Transform tr_inv = tr.inverse(); + for (const FTransform& im : images) + new_cell.images.push_back(tr.combine(im).combine(tr_inv)); + } + return new_cell; + } + + UnitCell changed_basis_forward(const Op& op, bool set_images) { + return changed_basis_backward(op.inverse(), set_images); + } + + bool is_compatible_with_groupops(const GroupOps& gops, double eps=1e-3) const { + std::array metric = metric_tensor().elements_voigt(); + for (const Op& op : gops.sym_ops) { + Mat33 m = orth.mat.multiply(rot_as_mat33(op)); + std::array other = {{ + m.column_dot(0,0), m.column_dot(1,1), m.column_dot(2,2), + m.column_dot(1,2), m.column_dot(0,2), m.column_dot(0,1) + }}; + for (int i = 0; i < 6; ++i) + if (std::fabs(metric[i] - other[i]) > eps) + return false; + } + return true; + } + + bool is_compatible_with_spacegroup(const SpaceGroup* sg, double eps=1e-3) const { + return sg ? is_compatible_with_groupops(sg->operations(), eps) : false; + } + + void set_cell_images_from_groupops(const GroupOps& group_ops) { + images.clear(); + cs_count = (short) group_ops.order() - 1; + images.reserve(cs_count); + 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; + } + } + + void add_ncs_images_to_cs_images(const std::vector& ncs) { + assert(cs_count == (short) images.size()); + for (const NcsOp& ncs_op : ncs) + if (!ncs_op.given) { + // We need it to operates on fractional, not orthogonal coordinates. + FTransform f = frac.combine(ncs_op.tr.combine(orth)); + images.push_back(f); + for (int i = 0; i < cs_count; ++i) + images.push_back(images[i].combine(f)); + } + } + + std::vector get_ncs_transforms() const { + std::vector ncs; + for (size_t n = cs_count; n < images.size(); n += cs_count + 1) + ncs.push_back(images[n]); + return ncs; + } + + Position orthogonalize(const Fractional& f) const { + return Position(orth.apply(f)); + } + Fractional fractionalize(const Position& o) const { + return Fractional(frac.apply(o)); + } + + /// orthogonalize_difference(a-b) == orthogonalize(a) - orthogonalize(b) + // The shift (fract.vec) can be non-zero in non-standard settings, + // just do not apply it here. + Position orthogonalize_difference(const Fractional& delta) const { + return Position(orth.mat.multiply(delta)); + } + /// the inverse of orthogonalize_difference + Fractional fractionalize_difference(const Position& delta) const { + return Fractional(frac.mat.multiply(delta)); + } + + /// Returns box containing fractional box (a cuboid in fractional + /// coordinates can be a parallelepiped in Cartesian coordinates). + Box orthogonalize_box(const Box& f) const { + Box r; + r.minimum = orthogonalize(f.minimum); + r.maximum = orthogonalize(f.maximum); + if (alpha != 90. || beta == 90. || gamma == 90.) { + r.extend(orthogonalize({f.minimum.x, f.minimum.y, f.maximum.z})); + r.extend(orthogonalize({f.minimum.x, f.maximum.y, f.maximum.z})); + r.extend(orthogonalize({f.minimum.x, f.maximum.y, f.minimum.z})); + r.extend(orthogonalize({f.maximum.x, f.maximum.y, f.minimum.z})); + r.extend(orthogonalize({f.maximum.x, f.minimum.y, f.minimum.z})); + r.extend(orthogonalize({f.maximum.x, f.minimum.y, f.maximum.z})); + } + return r; + } + + Transform orthogonalize_transform(const FTransform& ftr) const { + return orth.combine(ftr.combine(frac)); + } + Transform op_as_transform(const Op& op) const { + return orthogonalize_transform(Transform{rot_as_mat33(op), tran_as_vec3(op)}); + } + + double distance_sq(const Fractional& pos1, const Fractional& pos2) const { + Fractional diff = (pos1 - pos2).wrap_to_zero(); + return orthogonalize_difference(diff).length_sq(); + } + double distance_sq(const Position& pos1, const Position& pos2) const { + return distance_sq(fractionalize(pos1), fractionalize(pos2)); + } + + double volume_per_image() const { + return is_crystal() ? volume / (1 + images.size()) : NAN; + } + + // Helper function. PBC = periodic boundary conditions. + bool search_pbc_images(Fractional&& diff, NearestImage& image) const { + int neg_shift[3] = {0, 0, 0}; + if (is_crystal()) { + for (int j = 0; j < 3; ++j) + neg_shift[j] = iround(diff.at(j)); + diff.x -= neg_shift[0]; + diff.y -= neg_shift[1]; + diff.z -= neg_shift[2]; + } + Position orth_diff = orthogonalize_difference(diff); + double dsq = orth_diff.length_sq(); + if (dsq < image.dist_sq) { + image.dist_sq = dsq; + for (int j = 0; j < 3; ++j) + image.pbc_shift[j] = -neg_shift[j]; + return true; + } + return false; + } + + NearestImage find_nearest_image(const Position& ref, const Position& pos, Asu asu) const { + NearestImage image; + if (asu == Asu::Different) + image.dist_sq = INFINITY; + else + image.dist_sq = ref.dist_sq(pos); + if (asu == Asu::Same) + return image; + Fractional fpos = fractionalize(pos); + Fractional fref = fractionalize(ref); + search_pbc_images(fpos - fref, image); + if (asu == Asu::Different && + image.pbc_shift[0] == 0 && image.pbc_shift[1] == 0 && image.pbc_shift[2] == 0) + image.dist_sq = INFINITY; + for (int n = 0; n != static_cast(images.size()); ++n) + if (search_pbc_images(images[n].apply(fpos) - fref, image)) + image.sym_idx = n + 1; + return image; + } + + void apply_transform(Fractional& fpos, int image_idx, bool inverse) const { + if (image_idx > 0) { + const FTransform& t = images.at(image_idx - 1); + if (!inverse) + fpos = t.apply(fpos); + else + fpos = FTransform(t.inverse()).apply(fpos); + } + } + + NearestImage find_nearest_pbc_image(const Fractional& fref, Fractional fpos, + int image_idx=0) const { + NearestImage sym_image; + sym_image.dist_sq = INFINITY; + sym_image.sym_idx = image_idx; + apply_transform(fpos, image_idx, false); + search_pbc_images(fpos - fref, sym_image); + return sym_image; + } + NearestImage find_nearest_pbc_image(const Position& ref, const Position& pos, + int image_idx=0) const { + return find_nearest_pbc_image(fractionalize(ref), fractionalize(pos), image_idx); + } + + std::vector find_nearest_pbc_images(const Fractional& fref, double dist, + const Fractional& fpos, int image_idx) const { + std::vector results; + NearestImage im = find_nearest_pbc_image(fref, fpos, image_idx); + int sh[3] = {im.pbc_shift[0], im.pbc_shift[1], im.pbc_shift[2]}; + for (im.pbc_shift[0] = sh[0]-1; im.pbc_shift[0] <= sh[0]+1; ++im.pbc_shift[0]) + for (im.pbc_shift[1] = sh[1]-1; im.pbc_shift[1] <= sh[1]+1; ++im.pbc_shift[1]) + for (im.pbc_shift[2] = sh[2]-1; im.pbc_shift[2] <= sh[2]+1; ++im.pbc_shift[2]) { + Fractional shift(im.pbc_shift[0], im.pbc_shift[1], im.pbc_shift[2]); + im.dist_sq = orthogonalize_difference(fpos - fref + shift).length_sq(); + if (im.dist_sq <= sq(dist)) + results.push_back(im); + } + return results; + } + + Position orthogonalize_in_pbc(const Position& ref, + const Fractional& fpos) const { + Fractional fref = fractionalize(ref); + return orthogonalize_difference((fpos - fref).wrap_to_zero()) + ref; + } + + Position find_nearest_pbc_position(const Position& ref, const Position& pos, + int image_idx, bool inverse=false) const { + Fractional fpos = fractionalize(pos); + apply_transform(fpos, image_idx, inverse); + return orthogonalize_in_pbc(ref, fpos); + } + + // apply NearestImage symmetry to fpos + Fractional fract_image(const NearestImage& im, Fractional fpos) { + apply_transform(fpos, im.sym_idx, false); + return fpos + Fractional(im.pbc_shift[0], im.pbc_shift[1], im.pbc_shift[2]); + } + + /// Counts nearby symmetry mates (0 = none, 3 = 4-fold axis, etc). + /// \pre is_crystal() + int is_special_position(const Fractional& fpos, double max_dist) const { + const double max_dist_sq = max_dist * max_dist; + int n = 0; + for (const FTransform& image : images) { + Fractional fdiff = (image.apply(fpos) - fpos).wrap_to_zero(); + if (orthogonalize_difference(fdiff).length_sq() < max_dist_sq) + ++n; + } + return n; + } + int is_special_position(const Position& pos, double max_dist = 0.8) const { + return is_special_position(fractionalize(pos), max_dist); + } + + /// Calculate 1/d^2 for specified hkl reflection. + /// 1/d^2 = (2*sin(theta)/lambda)^2 + // The indices are integers, but they may be stored as floating-point + // numbers (MTZ format) so we use double to avoid conversions. + double calculate_1_d2_double(double h, double k, double l) const { + double arh = ar * h; + double brk = br * k; + double crl = cr * l; + return arh * arh + brk * brk + crl * crl + 2 * (arh * brk * cos_gammar + + arh * crl * cos_betar + + brk * crl * cos_alphar); + } + double calculate_1_d2(const Miller& hkl) const { + return calculate_1_d2_double(hkl[0], hkl[1], hkl[2]); + } + + /// Calculate d-spacing. + /// d = lambda/(2*sin(theta)) + double calculate_d(const Miller& hkl) const { + return 1.0 / std::sqrt(calculate_1_d2(hkl)); + } + + /// Calculate (sin(theta)/lambda)^2 = d*^2/4 + double calculate_stol_sq(const Miller& hkl) const { + return 0.25 * calculate_1_d2(hkl); + } + + /// https://dictionary.iucr.org/Metric_tensor + SMat33 metric_tensor() const { + // the order in SMat33 is ... m12 m13 m23 -> a.a b.b c.c a.b a.c b.c + return {a*a, b*b, c*c, a*orth.mat[0][1], a*orth.mat[0][2], b*c*cos_alpha()}; + } + + SMat33 reciprocal_metric_tensor() const { + return {ar*ar, br*br, cr*cr, ar*br*cos_gammar, ar*cr*cos_betar, br*cr*cos_alphar}; + } + + /// Returns reciprocal unit cell. + UnitCell reciprocal() const { + auto acosd = [](double x) { return deg(std::acos(x)); }; + return UnitCell(ar, br, cr, + acosd(cos_alphar), acosd(cos_betar), acosd(cos_gammar)); + } + + Miller get_hkl_limits(double dmin) const { + return {{int(a / dmin), int(b / dmin), int(c / dmin)}}; + } + + Mat33 primitive_orth_matrix(char centring_type) const { + if (centring_type == 'P') + return orth.mat; + Mat33 c2p = rot_as_mat33(centred_to_primitive(centring_type)); + return orth.mat.multiply(c2p); + } +}; + +} // namespace gemmi +#endif diff --git a/symmetry/symmetry.cpp b/symmetry/symmetry.cpp new file mode 100644 index 000000000..05b5d1c07 --- /dev/null +++ b/symmetry/symmetry.cpp @@ -0,0 +1,1215 @@ +// Copyright Global Phasing Ltd. + +#include +#include // for fabs +#include // for memchr, strchr + +static const char* skip_space(const char* p) { + if (p) + while (*p == ' ' || *p == '\t' || *p == '_') // '_' can be used as space + ++p; + return p; +} + +namespace gemmi { + +// TRIPLET -> OP + +// param only can be set to 'h', 'x', 'a' or ' ' (any), to limit accepted characters. +// decimal_fract is useful only for non-crystallographic ops (such as x+0.12) +std::array parse_triplet_part(const std::string& s, char& notation, double* decimal_fract) { + constexpr char a_ = 'a' & ~3; + constexpr char h_ = 'h' & ~3; + constexpr char x_ = 'x' & ~3; + static const signed char letter2index[] = + // a b c d e f g h i j k l + { a_+0, a_+1, a_+2, 0, 0, 0, 0, h_+0, 0, 0, h_+1, h_+2, + // m n o p q r s t u v w x y z + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, x_+0, x_+1, x_+2 }; + auto interpret_letter = [&](char c) { + size_t idx = size_t((c | 0x20) - 'a'); // "|0x20" = to lower + if (idx >= sizeof(letter2index) || letter2index[idx] == 0) + fail("unexpected character '", c, "' in: ", s); + auto value = letter2index[idx]; + int detected_notation = value & ~3; + if ((notation | 0x20) == ' ') + notation = detected_notation; + else if (((notation | 0x20) & ~3) != detected_notation) + fail("Unexpected notation (letter set) in: ", s); + return value & 3; + }; + + std::array r = { 0, 0, 0, 0 }; + int num = Op::DEN; + const char* c = s.c_str(); + while (*(c = skip_space(c))) { + if (*c == '+' || *c == '-') { + num = (*c == '+' ? Op::DEN : -Op::DEN); + c = skip_space(++c); + } + if (num == 0) + fail("wrong or unsupported triplet format: " + s); + int r_idx; + int den = 1; + double fract = 0; + if ((*c >= '0' && *c <= '9') || *c == '.') { + // syntax examples in this branch: "1", "-1/2", "+2*x", "1/2 * b" + char* endptr; + int n = std::strtol(c, &endptr, 10); + // some COD CIFs have decimal fractions ("-x+0.25", ".5+Y", "1.25000-y") + if (*endptr == '.') { + // avoiding strtod() etc which is locale-dependent + fract = n; + for (double denom = 0.1; *++endptr >= '0' && *endptr <= '9'; denom *= 0.1) + fract += int(*endptr - '0') * denom; + double rounded = std::round(fract * num); + if (!decimal_fract) { + if (std::fabs(rounded - fract * num) > 0.05) + fail("unexpected number in a symmetry triplet part: " + s); + num = int(rounded); + } + } else { + num *= n; + } + if (*endptr == '/') + den = std::strtol(endptr + 1, &endptr, 10); + if (*endptr == '*') { + c = skip_space(endptr + 1); + r_idx = interpret_letter(*c); + ++c; + } else { + c = endptr; + r_idx = 3; + } + } else { + // syntax examples in this branch: "x", "+a", "-k/3" + r_idx = interpret_letter(*c); + c = skip_space(++c); + if (*c == '/') { + char* endptr; + den = std::strtol(c + 1, &endptr, 10); + c = endptr; + } + } + if (den != 1) { + if (den <= 0 || Op::DEN % den != 0 || fract != 0) + fail("Wrong denominator " + std::to_string(den) + " in: " + s); + num /= den; + } + r[r_idx] += num; + if (decimal_fract) + decimal_fract[r_idx] = num > 0 ? fract : -fract; + num = 0; + } + if (num != 0) + fail("trailing sign in: " + s); + return r; +} + +Op parse_triplet(const std::string& s, char notation) { + if (std::count(s.begin(), s.end(), ',') != 2) + fail("expected exactly two commas in triplet"); + size_t comma1 = s.find(','); + size_t comma2 = s.find(',', comma1 + 1); + char save_notation = notation; + notation = (notation | 0x20) & ~3; + if (notation != 'x' && notation != 'h' && notation != '`' && notation != ' ') // '`' == a' & ~3 + fail("parse_triplet(): unexpected notation='", save_notation, "'"); + auto a = parse_triplet_part(s.substr(0, comma1), notation); + auto b = parse_triplet_part(s.substr(comma1 + 1, comma2 - (comma1 + 1)), notation); + auto c = parse_triplet_part(s.substr(comma2 + 1), notation); + Op::Rot rot = {{{a[0], a[1], a[2]}, {b[0], b[1], b[2]}, {c[0], c[1], c[2]}}}; + Op::Tran tran = {a[3], b[3], c[3]}; + if (notation == 'h') { + if (tran != Op::Tran{0, 0, 0}) + fail("parse_triplet(): reciprocal-space Op cannot have translation: ", s); + rot = Op::transpose(rot); + } + return { rot, tran, notation }; +} + + +// OP -> TRIPLET + +namespace { + +// much faster than s += std::to_string(n) for n in 0 ... 99 +void append_small_number(std::string& s, int n) { + if (n < 0 || n >= 100) { + s += std::to_string(n); + } else if (n < 10) { + s += char('0' + n); + } else { // 10 ... 99 + int tens = n / 10; + s += char('0' + tens); + s += char('0' + n - 10 * tens); + } +} + +void append_sign_of(std::string& s, int n) { + if (n < 0) + s += '-'; + else if (!s.empty()) + s += '+'; +} + +// append w/DEN fraction reduced to the lowest terms +std::pair get_op_fraction(int w) { + // Op::DEN == 24 == 2 * 2 * 2 * 3 + int denom = 1; + for (int i = 0; i != 3; ++i) + if (w % 2 == 0) // 2, 2, 2 + w /= 2; + else + denom *= 2; + if (w % 3 == 0) // 3 + w /= 3; + else + denom *= 3; + return {w, denom}; +} + +void append_fraction(std::string& s, std::pair frac) { + append_small_number(s, frac.first); + if (frac.second != 1) { + s += '/'; + append_small_number(s, frac.second); + } +} + +std::string make_triplet_part(const std::array& xyz, int w, char style) { + std::string s; + const char* letters = "xyz hkl abc XYZ HKL ABC"; + switch((style | 0x20) & ~3) { // |0x20 converts to lower case + case 'h': letters += 4; break; + case '`': letters += 8; break; // 'a', because 'a'&~3 == 0x60 == '`' + } + if (!(style & 0x20)) // not lower + letters += 12; + for (int i = 0; i != 3; ++i) + if (xyz[i] != 0) { + append_sign_of(s, xyz[i]); + int a = std::abs(xyz[i]); + if (a != Op::DEN) { + std::pair frac = get_op_fraction(a); + if (frac.first == 1) { // e.g. "x/3" + s += letters[i]; + s += '/'; + append_small_number(s, frac.second); + } else { // e.g. "2/3*x" + append_fraction(s, frac); + s += '*'; + s += letters[i]; + } + } else { + s += letters[i]; + } + } + if (w != 0) { + append_sign_of(s, w); + std::pair frac = get_op_fraction(std::abs(w)); + append_fraction(s, frac); + } + return s; +} + +} // anonymous namespace + +Op seitz_to_op(const std::array, 4>& t) { + static_assert(Op::DEN == 24, ""); + auto check_round = [](double d) { + double r = std::round(d * Op::DEN); + if (std::fabs(r - d * Op::DEN) > 0.05) + fail("all numbers in Seitz matrix must be equal Z/24"); + return static_cast(r); + }; + Op op; + if (std::fabs(t[3][0]) + std::fabs(t[3][1]) + std::fabs(t[3][2]) + + std::fabs(t[3][3] - 1) > 1e-3) + fail("the last row in Seitz matrix must be [0 0 0 1]"); + for (int i = 0; i < 3; ++i) { + for (int j = 0; j < 3; ++j) + op.rot[i][j] = check_round(t[i][j]); + op.tran[i] = check_round(t[i][3]); + } + op.notation = 'x'; + return op; +} + +void append_op_fraction(std::string& s, int w) { + append_fraction(s, get_op_fraction(w)); +} + +std::string Op::triplet(char style) const { + if (style == ' ') + style = (notation & ~0x20) ? notation : 'x'; + char lower_style = (style | 0x20) & ~3; + if (lower_style == 'h' && !is_hkl()) + fail("triplet(): can't write real-space triplet as hkl"); + if (lower_style != 'h' && is_hkl()) + fail("triplet(): can't write reciprocal-space triplet as xyz"); + // 'x'==0x78, 'h'==0x68, 'a'==0x61, so 'a'&~3 == 0x60 == '`' + if (lower_style != 'x' && lower_style != 'h' && lower_style != '`') + fail("unexpected triplet style: '", style, "'"); + // parse_triplet() transposes hkl ops such as l,h,k + auto r = !is_hkl()? rot : transposed_rot(); + return make_triplet_part(r[0], tran[0], style) + + "," + make_triplet_part(r[1], tran[1], style) + + "," + make_triplet_part(r[2], tran[2], style); +} + + +// INTERPRETING HALL SYMBOLS +// based on both ITfC vol.B ch.1.4 (2010) +// and http://cci.lbl.gov/sginfo/hall_symbols.html + +// matrices for Nz from Table 3 and 4 from hall_symbols.html +namespace { +Op::Rot hall_rotation_z(int N) { + constexpr int d = Op::DEN; + switch (N) { + case 1: return {{{d,0,0}, {0,d,0}, {0,0,d}}}; + case 2: return {{{-d,0,0}, {0,-d,0}, {0,0,d}}}; + case 3: return {{{0,-d,0}, {d,-d,0}, {0,0,d}}}; + case 4: return {{{0,-d,0}, {d,0,0}, {0,0,d}}}; + case 6: return {{{d,-d,0}, {d,0,0}, {0,0,d}}}; + case '\'': return {{{0,-d,0},{-d,0,0}, {0,0,-d}}}; + case '"': return {{{0,d,0}, { d,0,0}, {0,0,-d}}}; + case '*': return {{{0,0,d}, { d,0,0}, {0,d,0}}}; + default: fail("incorrect axis definition"); + } +} +Op::Tran hall_translation_from_symbol(char symbol) { + constexpr int h = Op::DEN / 2; + constexpr int q = Op::DEN / 4; + switch (symbol) { + case 'a': return {h, 0, 0}; + case 'b': return {0, h, 0}; + case 'c': return {0, 0, h}; + case 'n': return {h, h, h}; + case 'u': return {q, 0, 0}; + case 'v': return {0, q, 0}; + case 'w': return {0, 0, q}; + case 'd': return {q, q, q}; + default: fail(std::string("unknown symbol: ") + symbol); + } +} + +Op hall_matrix_symbol(const char* start, const char* end, int pos, int& prev) { + Op op = Op::identity(); + bool neg = (*start == '-'); + const char* p = (neg ? start + 1 : start); + if (*p < '1' || *p == '5' || *p > '6') + fail("wrong n-fold order notation: " + std::string(start, end)); + int N = *p++ - '0'; + int fractional_tran = 0; + char principal_axis = '\0'; + char diagonal_axis = '\0'; + for (; p < end; ++p) { + if (*p >= '1' && *p <= '5') { + if (fractional_tran != '\0') + fail("two numeric subscripts"); + fractional_tran = *p - '0'; + } else if (*p == '\'' || *p == '"' || *p == '*') { + if (N != (*p == '*' ? 3 : 2)) + fail("wrong symbol: " + std::string(start, end)); + diagonal_axis = *p; + } else if (*p == 'x' || *p == 'y' || *p == 'z') { + principal_axis = *p; + } else { + op.translate(hall_translation_from_symbol(*p)); + } + } + // fill in implicit values + if (!principal_axis && !diagonal_axis) { + if (pos == 1) { + principal_axis = 'z'; + } else if (pos == 2 && N == 2) { + if (prev == 2 || prev == 4) + principal_axis = 'x'; + else if (prev == 3 || prev == 6) + diagonal_axis = '\''; + } else if (pos == 3 && N == 3) { + diagonal_axis = '*'; + } else if (N != 1) { + fail("missing axis"); + } + } + // get the operation + op.rot = hall_rotation_z(diagonal_axis ? diagonal_axis : N); + if (neg) + op.rot = op.negated_rot(); + auto alter_order = [](const Op::Rot& r, int i, int j, int k) { + return Op::Rot{{ {r[i][i], r[i][j], r[i][k]}, + {r[j][i], r[j][j], r[j][k]}, + {r[k][i], r[k][j], r[k][k]} }}; + }; + if (principal_axis == 'x') + op.rot = alter_order(op.rot, 2, 0, 1); + else if (principal_axis == 'y') + op.rot = alter_order(op.rot, 1, 2, 0); + if (fractional_tran) + op.tran[principal_axis - 'x'] += Op::DEN / N * fractional_tran; + prev = N; + return op; +} + +// Parses either short (0 0 1) or long notation (x,y,z+1/12) +// but without multipliers (such as 1/2x) to keep things simple for now. +Op parse_hall_change_of_basis(const char* start, const char* end) { + if (std::memchr(start, ',', end - start) != nullptr) // long symbol + return parse_triplet(std::string(start, end)); + // short symbol (0 0 1) + Op cob = Op::identity(); + char* endptr; + for (int i = 0; i != 3; ++i) { + cob.tran[i] = std::strtol(start, &endptr, 10) % 12 * (Op::DEN / 12); + start = endptr; + } + if (endptr != end) + fail("unexpected change-of-basis format: " + std::string(start, end)); + return cob; +} +} // anonymous namespace + +GroupOps generators_from_hall(const char* hall) { + auto find_blank = [](const char* p) { + while (*p != '\0' && *p != ' ' && *p != '\t' && *p != '_') // '_' == ' ' + ++p; + return p; + }; + if (hall == nullptr) + fail("null"); + hall = skip_space(hall); + GroupOps ops; + ops.sym_ops.emplace_back(Op::identity()); + bool centrosym = (hall[0] == '-'); + const char* lat = skip_space(centrosym ? hall + 1 : hall); + if (!lat) + fail("not a hall symbol: " + std::string(hall)); + ops.cen_ops = centring_vectors(*lat); + int counter = 0; + int prev = 0; + const char* part = skip_space(lat + 1); + while (*part != '\0' && *part != '(') { + const char* space = find_blank(part); + ++counter; + if (part[0] != '1' || (part[1] != ' ' && part[1] != '\0')) { + Op op = hall_matrix_symbol(part, space, counter, prev); + ops.sym_ops.emplace_back(op); + } + part = skip_space(space); + } + if (centrosym) + ops.sym_ops.push_back({Op::identity().negated_rot(), {0,0,0}, 'x'}); + if (*part == '(') { + const char* rb = std::strchr(part, ')'); + if (!rb) + fail("missing ')': " + std::string(hall)); + if (ops.sym_ops.empty()) + fail("misplaced translation: " + std::string(hall)); + ops.change_basis_forward(parse_hall_change_of_basis(part + 1, rb)); + + if (*skip_space(find_blank(rb + 1)) != '\0') + fail("unexpected characters after ')': " + std::string(hall)); + } + return ops; +} + + +const SpaceGroup spacegroup_tables::main[564] = { + // This table was generated by tools/gen_sg_table.py. + // First 530 entries in the same order as in SgInfo, sgtbx and ITB. + // Note: spacegroup 68 has three duplicates with different H-M names. + { 1, 1, "P 1" , 0, "", "P 1" , 0 }, // 0 + { 2, 2, "P -1" , 0, "", "-P 1" , 0 }, // 1 + { 3, 3, "P 1 2 1" , 0, "b", "P 2y" , 0 }, // 2 + { 3, 1003, "P 1 1 2" , 0, "c", "P 2" , 1 }, // 3 + { 3, 0, "P 2 1 1" , 0, "a", "P 2x" , 2 }, // 4 + { 4, 4, "P 1 21 1" , 0, "b", "P 2yb" , 0 }, // 5 + { 4, 1004, "P 1 1 21" , 0, "c", "P 2c" , 1 }, // 6 + { 4, 0, "P 21 1 1" , 0, "a", "P 2xa" , 2 }, // 7 + { 5, 5, "C 1 2 1" , 0, "b1", "C 2y" , 0 }, // 8 + { 5, 2005, "A 1 2 1" , 0, "b2", "A 2y" , 3 }, // 9 + { 5, 4005, "I 1 2 1" , 0, "b3", "I 2y" , 4 }, // 10 + { 5, 0, "A 1 1 2" , 0, "c1", "A 2" , 1 }, // 11 + { 5, 1005, "B 1 1 2" , 0, "c2", "B 2" , 5 }, // 12 + { 5, 0, "I 1 1 2" , 0, "c3", "I 2" , 6 }, // 13 + { 5, 0, "B 2 1 1" , 0, "a1", "B 2x" , 2 }, // 14 + { 5, 0, "C 2 1 1" , 0, "a2", "C 2x" , 7 }, // 15 + { 5, 0, "I 2 1 1" , 0, "a3", "I 2x" , 8 }, // 16 + { 6, 6, "P 1 m 1" , 0, "b", "P -2y" , 0 }, // 17 + { 6, 1006, "P 1 1 m" , 0, "c", "P -2" , 1 }, // 18 + { 6, 0, "P m 1 1" , 0, "a", "P -2x" , 2 }, // 19 + { 7, 7, "P 1 c 1" , 0, "b1", "P -2yc" , 0 }, // 20 + { 7, 0, "P 1 n 1" , 0, "b2", "P -2yac" , 9 }, // 21 + { 7, 0, "P 1 a 1" , 0, "b3", "P -2ya" , 3 }, // 22 + { 7, 0, "P 1 1 a" , 0, "c1", "P -2a" , 1 }, // 23 + { 7, 0, "P 1 1 n" , 0, "c2", "P -2ab" , 10}, // 24 + { 7, 1007, "P 1 1 b" , 0, "c3", "P -2b" , 5 }, // 25 + { 7, 0, "P b 1 1" , 0, "a1", "P -2xb" , 2 }, // 26 + { 7, 0, "P n 1 1" , 0, "a2", "P -2xbc" , 11}, // 27 + { 7, 0, "P c 1 1" , 0, "a3", "P -2xc" , 7 }, // 28 + { 8, 8, "C 1 m 1" , 0, "b1", "C -2y" , 0 }, // 29 + { 8, 0, "A 1 m 1" , 0, "b2", "A -2y" , 3 }, // 30 + { 8, 0, "I 1 m 1" , 0, "b3", "I -2y" , 4 }, // 31 + { 8, 0, "A 1 1 m" , 0, "c1", "A -2" , 1 }, // 32 + { 8, 1008, "B 1 1 m" , 0, "c2", "B -2" , 5 }, // 33 + { 8, 0, "I 1 1 m" , 0, "c3", "I -2" , 6 }, // 34 + { 8, 0, "B m 1 1" , 0, "a1", "B -2x" , 2 }, // 35 + { 8, 0, "C m 1 1" , 0, "a2", "C -2x" , 7 }, // 36 + { 8, 0, "I m 1 1" , 0, "a3", "I -2x" , 8 }, // 37 + { 9, 9, "C 1 c 1" , 0, "b1", "C -2yc" , 0 }, // 38 + { 9, 0, "A 1 n 1" , 0, "b2", "A -2yab" , 12}, // 39 + { 9, 0, "I 1 a 1" , 0, "b3", "I -2ya" , 13}, // 40 + { 9, 0, "A 1 a 1" , 0, "-b1", "A -2ya" , 3 }, // 41 + { 9, 0, "C 1 n 1" , 0, "-b2", "C -2yac" , 14}, // 42 + { 9, 0, "I 1 c 1" , 0, "-b3", "I -2yc" , 4 }, // 43 + { 9, 0, "A 1 1 a" , 0, "c1", "A -2a" , 1 }, // 44 + { 9, 0, "B 1 1 n" , 0, "c2", "B -2ab" , 15}, // 45 + { 9, 0, "I 1 1 b" , 0, "c3", "I -2b" , 16}, // 46 + { 9, 1009, "B 1 1 b" , 0, "-c1", "B -2b" , 5 }, // 47 + { 9, 0, "A 1 1 n" , 0, "-c2", "A -2ab" , 10}, // 48 + { 9, 0, "I 1 1 a" , 0, "-c3", "I -2a" , 6 }, // 49 + { 9, 0, "B b 1 1" , 0, "a1", "B -2xb" , 2 }, // 50 + { 9, 0, "C n 1 1" , 0, "a2", "C -2xac" , 17}, // 51 + { 9, 0, "I c 1 1" , 0, "a3", "I -2xc" , 18}, // 52 + { 9, 0, "C c 1 1" , 0, "-a1", "C -2xc" , 7 }, // 53 + { 9, 0, "B n 1 1" , 0, "-a2", "B -2xab" , 11}, // 54 + { 9, 0, "I b 1 1" , 0, "-a3", "I -2xb" , 8 }, // 55 + { 10, 10, "P 1 2/m 1" , 0, "b", "-P 2y" , 0 }, // 56 + { 10, 1010, "P 1 1 2/m" , 0, "c", "-P 2" , 1 }, // 57 + { 10, 0, "P 2/m 1 1" , 0, "a", "-P 2x" , 2 }, // 58 + { 11, 11, "P 1 21/m 1", 0, "b", "-P 2yb" , 0 }, // 59 + { 11, 1011, "P 1 1 21/m", 0, "c", "-P 2c" , 1 }, // 60 + { 11, 0, "P 21/m 1 1", 0, "a", "-P 2xa" , 2 }, // 61 + { 12, 12, "C 1 2/m 1" , 0, "b1", "-C 2y" , 0 }, // 62 + { 12, 0, "A 1 2/m 1" , 0, "b2", "-A 2y" , 3 }, // 63 + { 12, 0, "I 1 2/m 1" , 0, "b3", "-I 2y" , 4 }, // 64 + { 12, 0, "A 1 1 2/m" , 0, "c1", "-A 2" , 1 }, // 65 + { 12, 1012, "B 1 1 2/m" , 0, "c2", "-B 2" , 5 }, // 66 + { 12, 0, "I 1 1 2/m" , 0, "c3", "-I 2" , 6 }, // 67 + { 12, 0, "B 2/m 1 1" , 0, "a1", "-B 2x" , 2 }, // 68 + { 12, 0, "C 2/m 1 1" , 0, "a2", "-C 2x" , 7 }, // 69 + { 12, 0, "I 2/m 1 1" , 0, "a3", "-I 2x" , 8 }, // 70 + { 13, 13, "P 1 2/c 1" , 0, "b1", "-P 2yc" , 0 }, // 71 + { 13, 0, "P 1 2/n 1" , 0, "b2", "-P 2yac" , 9 }, // 72 + { 13, 0, "P 1 2/a 1" , 0, "b3", "-P 2ya" , 3 }, // 73 + { 13, 0, "P 1 1 2/a" , 0, "c1", "-P 2a" , 1 }, // 74 + { 13, 0, "P 1 1 2/n" , 0, "c2", "-P 2ab" , 10}, // 75 + { 13, 1013, "P 1 1 2/b" , 0, "c3", "-P 2b" , 5 }, // 76 + { 13, 0, "P 2/b 1 1" , 0, "a1", "-P 2xb" , 2 }, // 77 + { 13, 0, "P 2/n 1 1" , 0, "a2", "-P 2xbc" , 11}, // 78 + { 13, 0, "P 2/c 1 1" , 0, "a3", "-P 2xc" , 7 }, // 79 + { 14, 14, "P 1 21/c 1", 0, "b1", "-P 2ybc" , 0 }, // 80 + { 14, 2014, "P 1 21/n 1", 0, "b2", "-P 2yn" , 9 }, // 81 + { 14, 3014, "P 1 21/a 1", 0, "b3", "-P 2yab" , 3 }, // 82 + { 14, 0, "P 1 1 21/a", 0, "c1", "-P 2ac" , 1 }, // 83 + { 14, 0, "P 1 1 21/n", 0, "c2", "-P 2n" , 10}, // 84 + { 14, 1014, "P 1 1 21/b", 0, "c3", "-P 2bc" , 5 }, // 85 + { 14, 0, "P 21/b 1 1", 0, "a1", "-P 2xab" , 2 }, // 86 + { 14, 0, "P 21/n 1 1", 0, "a2", "-P 2xn" , 11}, // 87 + { 14, 0, "P 21/c 1 1", 0, "a3", "-P 2xac" , 7 }, // 88 + { 15, 15, "C 1 2/c 1" , 0, "b1", "-C 2yc" , 0 }, // 89 + { 15, 0, "A 1 2/n 1" , 0, "b2", "-A 2yab" , 12}, // 90 + { 15, 0, "I 1 2/a 1" , 0, "b3", "-I 2ya" , 13}, // 91 + { 15, 0, "A 1 2/a 1" , 0, "-b1", "-A 2ya" , 3 }, // 92 + { 15, 0, "C 1 2/n 1" , 0, "-b2", "-C 2yac" , 19}, // 93 + { 15, 0, "I 1 2/c 1" , 0, "-b3", "-I 2yc" , 4 }, // 94 + { 15, 0, "A 1 1 2/a" , 0, "c1", "-A 2a" , 1 }, // 95 + { 15, 0, "B 1 1 2/n" , 0, "c2", "-B 2ab" , 15}, // 96 + { 15, 0, "I 1 1 2/b" , 0, "c3", "-I 2b" , 16}, // 97 + { 15, 1015, "B 1 1 2/b" , 0, "-c1", "-B 2b" , 5 }, // 98 + { 15, 0, "A 1 1 2/n" , 0, "-c2", "-A 2ab" , 10}, // 99 + { 15, 0, "I 1 1 2/a" , 0, "-c3", "-I 2a" , 6 }, // 100 + { 15, 0, "B 2/b 1 1" , 0, "a1", "-B 2xb" , 2 }, // 101 + { 15, 0, "C 2/n 1 1" , 0, "a2", "-C 2xac" , 17}, // 102 + { 15, 0, "I 2/c 1 1" , 0, "a3", "-I 2xc" , 18}, // 103 + { 15, 0, "C 2/c 1 1" , 0, "-a1", "-C 2xc" , 7 }, // 104 + { 15, 0, "B 2/n 1 1" , 0, "-a2", "-B 2xab" , 11}, // 105 + { 15, 0, "I 2/b 1 1" , 0, "-a3", "-I 2xb" , 8 }, // 106 + { 16, 16, "P 2 2 2" , 0, "", "P 2 2" , 0 }, // 107 + { 17, 17, "P 2 2 21" , 0, "", "P 2c 2" , 0 }, // 108 + { 17, 1017, "P 21 2 2" , 0, "cab", "P 2a 2a" , 1 }, // 109 + { 17, 2017, "P 2 21 2" , 0, "bca", "P 2 2b" , 2 }, // 110 + { 18, 18, "P 21 21 2" , 0, "", "P 2 2ab" , 0 }, // 111 + { 18, 3018, "P 2 21 21" , 0, "cab", "P 2bc 2" , 1 }, // 112 + { 18, 2018, "P 21 2 21" , 0, "bca", "P 2ac 2ac" , 2 }, // 113 + { 19, 19, "P 21 21 21", 0, "", "P 2ac 2ab" , 0 }, // 114 + { 20, 20, "C 2 2 21" , 0, "", "C 2c 2" , 0 }, // 115 + { 20, 0, "A 21 2 2" , 0, "cab", "A 2a 2a" , 1 }, // 116 + { 20, 0, "B 2 21 2" , 0, "bca", "B 2 2b" , 2 }, // 117 + { 21, 21, "C 2 2 2" , 0, "", "C 2 2" , 0 }, // 118 + { 21, 0, "A 2 2 2" , 0, "cab", "A 2 2" , 1 }, // 119 + { 21, 0, "B 2 2 2" , 0, "bca", "B 2 2" , 2 }, // 120 + { 22, 22, "F 2 2 2" , 0, "", "F 2 2" , 0 }, // 121 + { 23, 23, "I 2 2 2" , 0, "", "I 2 2" , 0 }, // 122 + { 24, 24, "I 21 21 21", 0, "", "I 2b 2c" , 0 }, // 123 + { 25, 25, "P m m 2" , 0, "", "P 2 -2" , 0 }, // 124 + { 25, 0, "P 2 m m" , 0, "cab", "P -2 2" , 1 }, // 125 + { 25, 0, "P m 2 m" , 0, "bca", "P -2 -2" , 2 }, // 126 + { 26, 26, "P m c 21" , 0, "", "P 2c -2" , 0 }, // 127 + { 26, 0, "P c m 21" , 0, "ba-c", "P 2c -2c" , 7 }, // 128 + { 26, 0, "P 21 m a" , 0, "cab", "P -2a 2a" , 1 }, // 129 + { 26, 0, "P 21 a m" , 0, "-cba", "P -2 2a" , 3 }, // 130 + { 26, 0, "P b 21 m" , 0, "bca", "P -2 -2b" , 2 }, // 131 + { 26, 0, "P m 21 b" , 0, "a-cb", "P -2b -2" , 5 }, // 132 + { 27, 27, "P c c 2" , 0, "", "P 2 -2c" , 0 }, // 133 + { 27, 0, "P 2 a a" , 0, "cab", "P -2a 2" , 1 }, // 134 + { 27, 0, "P b 2 b" , 0, "bca", "P -2b -2b" , 2 }, // 135 + { 28, 28, "P m a 2" , 0, "", "P 2 -2a" , 0 }, // 136 + { 28, 0, "P b m 2" , 0, "ba-c", "P 2 -2b" , 7 }, // 137 + { 28, 0, "P 2 m b" , 0, "cab", "P -2b 2" , 1 }, // 138 + { 28, 0, "P 2 c m" , 0, "-cba", "P -2c 2" , 3 }, // 139 + { 28, 0, "P c 2 m" , 0, "bca", "P -2c -2c" , 2 }, // 140 + { 28, 0, "P m 2 a" , 0, "a-cb", "P -2a -2a" , 5 }, // 141 + { 29, 29, "P c a 21" , 0, "", "P 2c -2ac" , 0 }, // 142 + { 29, 0, "P b c 21" , 0, "ba-c", "P 2c -2b" , 7 }, // 143 + { 29, 0, "P 21 a b" , 0, "cab", "P -2b 2a" , 1 }, // 144 + { 29, 0, "P 21 c a" , 0, "-cba", "P -2ac 2a" , 3 }, // 145 + { 29, 0, "P c 21 b" , 0, "bca", "P -2bc -2c" , 2 }, // 146 + { 29, 0, "P b 21 a" , 0, "a-cb", "P -2a -2ab" , 5 }, // 147 + { 30, 30, "P n c 2" , 0, "", "P 2 -2bc" , 0 }, // 148 + { 30, 0, "P c n 2" , 0, "ba-c", "P 2 -2ac" , 7 }, // 149 + { 30, 0, "P 2 n a" , 0, "cab", "P -2ac 2" , 1 }, // 150 + { 30, 0, "P 2 a n" , 0, "-cba", "P -2ab 2" , 3 }, // 151 + { 30, 0, "P b 2 n" , 0, "bca", "P -2ab -2ab" , 2 }, // 152 + { 30, 0, "P n 2 b" , 0, "a-cb", "P -2bc -2bc" , 5 }, // 153 + { 31, 31, "P m n 21" , 0, "", "P 2ac -2" , 0 }, // 154 + { 31, 0, "P n m 21" , 0, "ba-c", "P 2bc -2bc" , 7 }, // 155 + { 31, 0, "P 21 m n" , 0, "cab", "P -2ab 2ab" , 1 }, // 156 + { 31, 0, "P 21 n m" , 0, "-cba", "P -2 2ac" , 3 }, // 157 + { 31, 0, "P n 21 m" , 0, "bca", "P -2 -2bc" , 2 }, // 158 + { 31, 0, "P m 21 n" , 0, "a-cb", "P -2ab -2" , 5 }, // 159 + { 32, 32, "P b a 2" , 0, "", "P 2 -2ab" , 0 }, // 160 + { 32, 0, "P 2 c b" , 0, "cab", "P -2bc 2" , 1 }, // 161 + { 32, 0, "P c 2 a" , 0, "bca", "P -2ac -2ac" , 2 }, // 162 + { 33, 33, "P n a 21" , 0, "", "P 2c -2n" , 0 }, // 163 + { 33, 0, "P b n 21" , 0, "ba-c", "P 2c -2ab" , 7 }, // 164 + { 33, 0, "P 21 n b" , 0, "cab", "P -2bc 2a" , 1 }, // 165 + { 33, 0, "P 21 c n" , 0, "-cba", "P -2n 2a" , 3 }, // 166 + { 33, 0, "P c 21 n" , 0, "bca", "P -2n -2ac" , 2 }, // 167 + { 33, 0, "P n 21 a" , 0, "a-cb", "P -2ac -2n" , 5 }, // 168 + { 34, 34, "P n n 2" , 0, "", "P 2 -2n" , 0 }, // 169 + { 34, 0, "P 2 n n" , 0, "cab", "P -2n 2" , 1 }, // 170 + { 34, 0, "P n 2 n" , 0, "bca", "P -2n -2n" , 2 }, // 171 + { 35, 35, "C m m 2" , 0, "", "C 2 -2" , 0 }, // 172 + { 35, 0, "A 2 m m" , 0, "cab", "A -2 2" , 1 }, // 173 + { 35, 0, "B m 2 m" , 0, "bca", "B -2 -2" , 2 }, // 174 + { 36, 36, "C m c 21" , 0, "", "C 2c -2" , 0 }, // 175 + { 36, 0, "C c m 21" , 0, "ba-c", "C 2c -2c" , 7 }, // 176 + { 36, 0, "A 21 m a" , 0, "cab", "A -2a 2a" , 1 }, // 177 + { 36, 0, "A 21 a m" , 0, "-cba", "A -2 2a" , 3 }, // 178 + { 36, 0, "B b 21 m" , 0, "bca", "B -2 -2b" , 2 }, // 179 + { 36, 0, "B m 21 b" , 0, "a-cb", "B -2b -2" , 5 }, // 180 + { 37, 37, "C c c 2" , 0, "", "C 2 -2c" , 0 }, // 181 + { 37, 0, "A 2 a a" , 0, "cab", "A -2a 2" , 1 }, // 182 + { 37, 0, "B b 2 b" , 0, "bca", "B -2b -2b" , 2 }, // 183 + { 38, 38, "A m m 2" , 0, "", "A 2 -2" , 0 }, // 184 + { 38, 0, "B m m 2" , 0, "ba-c", "B 2 -2" , 7 }, // 185 + { 38, 0, "B 2 m m" , 0, "cab", "B -2 2" , 1 }, // 186 + { 38, 0, "C 2 m m" , 0, "-cba", "C -2 2" , 3 }, // 187 + { 38, 0, "C m 2 m" , 0, "bca", "C -2 -2" , 2 }, // 188 + { 38, 0, "A m 2 m" , 0, "a-cb", "A -2 -2" , 5 }, // 189 + { 39, 39, "A b m 2" , 0, "", "A 2 -2b" , 0 }, // 190 + { 39, 0, "B m a 2" , 0, "ba-c", "B 2 -2a" , 7 }, // 191 + { 39, 0, "B 2 c m" , 0, "cab", "B -2a 2" , 1 }, // 192 + { 39, 0, "C 2 m b" , 0, "-cba", "C -2a 2" , 3 }, // 193 + { 39, 0, "C m 2 a" , 0, "bca", "C -2a -2a" , 2 }, // 194 + { 39, 0, "A c 2 m" , 0, "a-cb", "A -2b -2b" , 5 }, // 195 + { 40, 40, "A m a 2" , 0, "", "A 2 -2a" , 0 }, // 196 + { 40, 0, "B b m 2" , 0, "ba-c", "B 2 -2b" , 7 }, // 197 + { 40, 0, "B 2 m b" , 0, "cab", "B -2b 2" , 1 }, // 198 + { 40, 0, "C 2 c m" , 0, "-cba", "C -2c 2" , 3 }, // 199 + { 40, 0, "C c 2 m" , 0, "bca", "C -2c -2c" , 2 }, // 200 + { 40, 0, "A m 2 a" , 0, "a-cb", "A -2a -2a" , 5 }, // 201 + { 41, 41, "A b a 2" , 0, "", "A 2 -2ab" , 0 }, // 202 + { 41, 0, "B b a 2" , 0, "ba-c", "B 2 -2ab" , 7 }, // 203 + { 41, 0, "B 2 c b" , 0, "cab", "B -2ab 2" , 1 }, // 204 + { 41, 0, "C 2 c b" , 0, "-cba", "C -2ac 2" , 3 }, // 205 + { 41, 0, "C c 2 a" , 0, "bca", "C -2ac -2ac" , 2 }, // 206 + { 41, 0, "A c 2 a" , 0, "a-cb", "A -2ab -2ab" , 5 }, // 207 + { 42, 42, "F m m 2" , 0, "", "F 2 -2" , 0 }, // 208 + { 42, 0, "F 2 m m" , 0, "cab", "F -2 2" , 1 }, // 209 + { 42, 0, "F m 2 m" , 0, "bca", "F -2 -2" , 2 }, // 210 + { 43, 43, "F d d 2" , 0, "", "F 2 -2d" , 0 }, // 211 + { 43, 0, "F 2 d d" , 0, "cab", "F -2d 2" , 1 }, // 212 + { 43, 0, "F d 2 d" , 0, "bca", "F -2d -2d" , 2 }, // 213 + { 44, 44, "I m m 2" , 0, "", "I 2 -2" , 0 }, // 214 + { 44, 0, "I 2 m m" , 0, "cab", "I -2 2" , 1 }, // 215 + { 44, 0, "I m 2 m" , 0, "bca", "I -2 -2" , 2 }, // 216 + { 45, 45, "I b a 2" , 0, "", "I 2 -2c" , 0 }, // 217 + { 45, 0, "I 2 c b" , 0, "cab", "I -2a 2" , 1 }, // 218 + { 45, 0, "I c 2 a" , 0, "bca", "I -2b -2b" , 2 }, // 219 + { 46, 46, "I m a 2" , 0, "", "I 2 -2a" , 0 }, // 220 + { 46, 0, "I b m 2" , 0, "ba-c", "I 2 -2b" , 7 }, // 221 + { 46, 0, "I 2 m b" , 0, "cab", "I -2b 2" , 1 }, // 222 + { 46, 0, "I 2 c m" , 0, "-cba", "I -2c 2" , 3 }, // 223 + { 46, 0, "I c 2 m" , 0, "bca", "I -2c -2c" , 2 }, // 224 + { 46, 0, "I m 2 a" , 0, "a-cb", "I -2a -2a" , 5 }, // 225 + { 47, 47, "P m m m" , 0, "", "-P 2 2" , 0 }, // 226 + { 48, 48, "P n n n" , '1', "", "P 2 2 -1n" , 20}, // 227 + { 48, 0, "P n n n" , '2', "", "-P 2ab 2bc" , 0 }, // 228 + { 49, 49, "P c c m" , 0, "", "-P 2 2c" , 0 }, // 229 + { 49, 0, "P m a a" , 0, "cab", "-P 2a 2" , 1 }, // 230 + { 49, 0, "P b m b" , 0, "bca", "-P 2b 2b" , 2 }, // 231 + { 50, 50, "P b a n" , '1', "", "P 2 2 -1ab" , 21}, // 232 + { 50, 0, "P b a n" , '2', "", "-P 2ab 2b" , 0 }, // 233 + { 50, 0, "P n c b" , '1', "cab", "P 2 2 -1bc" , 22}, // 234 + { 50, 0, "P n c b" , '2', "cab", "-P 2b 2bc" , 1 }, // 235 + { 50, 0, "P c n a" , '1', "bca", "P 2 2 -1ac" , 23}, // 236 + { 50, 0, "P c n a" , '2', "bca", "-P 2a 2c" , 2 }, // 237 + { 51, 51, "P m m a" , 0, "", "-P 2a 2a" , 0 }, // 238 + { 51, 0, "P m m b" , 0, "ba-c", "-P 2b 2" , 7 }, // 239 + { 51, 0, "P b m m" , 0, "cab", "-P 2 2b" , 1 }, // 240 + { 51, 0, "P c m m" , 0, "-cba", "-P 2c 2c" , 3 }, // 241 + { 51, 0, "P m c m" , 0, "bca", "-P 2c 2" , 2 }, // 242 + { 51, 0, "P m a m" , 0, "a-cb", "-P 2 2a" , 5 }, // 243 + { 52, 52, "P n n a" , 0, "", "-P 2a 2bc" , 0 }, // 244 + { 52, 0, "P n n b" , 0, "ba-c", "-P 2b 2n" , 7 }, // 245 + { 52, 0, "P b n n" , 0, "cab", "-P 2n 2b" , 1 }, // 246 + { 52, 0, "P c n n" , 0, "-cba", "-P 2ab 2c" , 3 }, // 247 + { 52, 0, "P n c n" , 0, "bca", "-P 2ab 2n" , 2 }, // 248 + { 52, 0, "P n a n" , 0, "a-cb", "-P 2n 2bc" , 5 }, // 249 + { 53, 53, "P m n a" , 0, "", "-P 2ac 2" , 0 }, // 250 + { 53, 0, "P n m b" , 0, "ba-c", "-P 2bc 2bc" , 7 }, // 251 + { 53, 0, "P b m n" , 0, "cab", "-P 2ab 2ab" , 1 }, // 252 + { 53, 0, "P c n m" , 0, "-cba", "-P 2 2ac" , 3 }, // 253 + { 53, 0, "P n c m" , 0, "bca", "-P 2 2bc" , 2 }, // 254 + { 53, 0, "P m a n" , 0, "a-cb", "-P 2ab 2" , 5 }, // 255 + { 54, 54, "P c c a" , 0, "", "-P 2a 2ac" , 0 }, // 256 + { 54, 0, "P c c b" , 0, "ba-c", "-P 2b 2c" , 7 }, // 257 + { 54, 0, "P b a a" , 0, "cab", "-P 2a 2b" , 1 }, // 258 + { 54, 0, "P c a a" , 0, "-cba", "-P 2ac 2c" , 3 }, // 259 + { 54, 0, "P b c b" , 0, "bca", "-P 2bc 2b" , 2 }, // 260 + { 54, 0, "P b a b" , 0, "a-cb", "-P 2b 2ab" , 5 }, // 261 + { 55, 55, "P b a m" , 0, "", "-P 2 2ab" , 0 }, // 262 + { 55, 0, "P m c b" , 0, "cab", "-P 2bc 2" , 1 }, // 263 + { 55, 0, "P c m a" , 0, "bca", "-P 2ac 2ac" , 2 }, // 264 + { 56, 56, "P c c n" , 0, "", "-P 2ab 2ac" , 0 }, // 265 + { 56, 0, "P n a a" , 0, "cab", "-P 2ac 2bc" , 1 }, // 266 + { 56, 0, "P b n b" , 0, "bca", "-P 2bc 2ab" , 2 }, // 267 + { 57, 57, "P b c m" , 0, "", "-P 2c 2b" , 0 }, // 268 + { 57, 0, "P c a m" , 0, "ba-c", "-P 2c 2ac" , 7 }, // 269 + { 57, 0, "P m c a" , 0, "cab", "-P 2ac 2a" , 1 }, // 270 + { 57, 0, "P m a b" , 0, "-cba", "-P 2b 2a" , 3 }, // 271 + { 57, 0, "P b m a" , 0, "bca", "-P 2a 2ab" , 2 }, // 272 + { 57, 0, "P c m b" , 0, "a-cb", "-P 2bc 2c" , 5 }, // 273 + { 58, 58, "P n n m" , 0, "", "-P 2 2n" , 0 }, // 274 + { 58, 0, "P m n n" , 0, "cab", "-P 2n 2" , 1 }, // 275 + { 58, 0, "P n m n" , 0, "bca", "-P 2n 2n" , 2 }, // 276 + { 59, 59, "P m m n" , '1', "", "P 2 2ab -1ab" , 21}, // 277 + { 59, 1059, "P m m n" , '2', "", "-P 2ab 2a" , 0 }, // 278 + { 59, 0, "P n m m" , '1', "cab", "P 2bc 2 -1bc" , 22}, // 279 + { 59, 0, "P n m m" , '2', "cab", "-P 2c 2bc" , 1 }, // 280 + { 59, 0, "P m n m" , '1', "bca", "P 2ac 2ac -1ac", 23}, // 281 + { 59, 0, "P m n m" , '2', "bca", "-P 2c 2a" , 2 }, // 282 + { 60, 60, "P b c n" , 0, "", "-P 2n 2ab" , 0 }, // 283 + { 60, 0, "P c a n" , 0, "ba-c", "-P 2n 2c" , 7 }, // 284 + { 60, 0, "P n c a" , 0, "cab", "-P 2a 2n" , 1 }, // 285 + { 60, 0, "P n a b" , 0, "-cba", "-P 2bc 2n" , 3 }, // 286 + { 60, 0, "P b n a" , 0, "bca", "-P 2ac 2b" , 2 }, // 287 + { 60, 0, "P c n b" , 0, "a-cb", "-P 2b 2ac" , 5 }, // 288 + { 61, 61, "P b c a" , 0, "", "-P 2ac 2ab" , 0 }, // 289 + { 61, 0, "P c a b" , 0, "ba-c", "-P 2bc 2ac" , 3 }, // 290 + { 62, 62, "P n m a" , 0, "", "-P 2ac 2n" , 0 }, // 291 + { 62, 0, "P m n b" , 0, "ba-c", "-P 2bc 2a" , 7 }, // 292 + { 62, 0, "P b n m" , 0, "cab", "-P 2c 2ab" , 1 }, // 293 + { 62, 0, "P c m n" , 0, "-cba", "-P 2n 2ac" , 3 }, // 294 + { 62, 0, "P m c n" , 0, "bca", "-P 2n 2a" , 2 }, // 295 + { 62, 0, "P n a m" , 0, "a-cb", "-P 2c 2n" , 5 }, // 296 + { 63, 63, "C m c m" , 0, "", "-C 2c 2" , 0 }, // 297 + { 63, 0, "C c m m" , 0, "ba-c", "-C 2c 2c" , 7 }, // 298 + { 63, 0, "A m m a" , 0, "cab", "-A 2a 2a" , 1 }, // 299 + { 63, 0, "A m a m" , 0, "-cba", "-A 2 2a" , 3 }, // 300 + { 63, 0, "B b m m" , 0, "bca", "-B 2 2b" , 2 }, // 301 + { 63, 0, "B m m b" , 0, "a-cb", "-B 2b 2" , 5 }, // 302 + { 64, 64, "C m c a" , 0, "", "-C 2ac 2" , 0 }, // 303 + { 64, 0, "C c m b" , 0, "ba-c", "-C 2ac 2ac" , 7 }, // 304 + { 64, 0, "A b m a" , 0, "cab", "-A 2ab 2ab" , 1 }, // 305 + { 64, 0, "A c a m" , 0, "-cba", "-A 2 2ab" , 3 }, // 306 + { 64, 0, "B b c m" , 0, "bca", "-B 2 2ab" , 2 }, // 307 + { 64, 0, "B m a b" , 0, "a-cb", "-B 2ab 2" , 5 }, // 308 + { 65, 65, "C m m m" , 0, "", "-C 2 2" , 0 }, // 309 + { 65, 0, "A m m m" , 0, "cab", "-A 2 2" , 1 }, // 310 + { 65, 0, "B m m m" , 0, "bca", "-B 2 2" , 2 }, // 311 + { 66, 66, "C c c m" , 0, "", "-C 2 2c" , 0 }, // 312 + { 66, 0, "A m a a" , 0, "cab", "-A 2a 2" , 1 }, // 313 + { 66, 0, "B b m b" , 0, "bca", "-B 2b 2b" , 2 }, // 314 + { 67, 67, "C m m a" , 0, "", "-C 2a 2" , 0 }, // 315 + { 67, 0, "C m m b" , 0, "ba-c", "-C 2a 2a" , 14}, // 316 + { 67, 0, "A b m m" , 0, "cab", "-A 2b 2b" , 1 }, // 317 + { 67, 0, "A c m m" , 0, "-cba", "-A 2 2b" , 3 }, // 318 + { 67, 0, "B m c m" , 0, "bca", "-B 2 2a" , 2 }, // 319 + { 67, 0, "B m a m" , 0, "a-cb", "-B 2a 2" , 5 }, // 320 + { 68, 68, "C c c a" , '1', "", "C 2 2 -1ac" , 24}, // 321 + { 68, 0, "C c c a" , '2', "", "-C 2a 2ac" , 0 }, // 322 + { 68, 0, "C c c b" , '1', "ba-c", "C 2 2 -1ac" , 24}, // 323 (==321) + { 68, 0, "C c c b" , '2', "ba-c", "-C 2a 2c" , 21}, // 324 + { 68, 0, "A b a a" , '1', "cab", "A 2 2 -1ab" , 25}, // 325 + { 68, 0, "A b a a" , '2', "cab", "-A 2a 2b" , 1 }, // 326 + { 68, 0, "A c a a" , '1', "-cba", "A 2 2 -1ab" , 25}, // 327 (==325) + { 68, 0, "A c a a" , '2', "-cba", "-A 2ab 2b" , 3 }, // 328 + { 68, 0, "B b c b" , '1', "bca", "B 2 2 -1ab" , 26}, // 329 + { 68, 0, "B b c b" , '2', "bca", "-B 2ab 2b" , 2 }, // 330 + { 68, 0, "B b a b" , '1', "a-cb", "B 2 2 -1ab" , 26}, // 331 (==329) + { 68, 0, "B b a b" , '2', "a-cb", "-B 2b 2ab" , 5 }, // 332 + { 69, 69, "F m m m" , 0, "", "-F 2 2" , 0 }, // 333 + { 70, 70, "F d d d" , '1', "", "F 2 2 -1d" , 27}, // 334 + { 70, 0, "F d d d" , '2', "", "-F 2uv 2vw" , 0 }, // 335 + { 71, 71, "I m m m" , 0, "", "-I 2 2" , 0 }, // 336 + { 72, 72, "I b a m" , 0, "", "-I 2 2c" , 0 }, // 337 + { 72, 0, "I m c b" , 0, "cab", "-I 2a 2" , 1 }, // 338 + { 72, 0, "I c m a" , 0, "bca", "-I 2b 2b" , 2 }, // 339 + { 73, 73, "I b c a" , 0, "", "-I 2b 2c" , 0 }, // 340 + { 73, 0, "I c a b" , 0, "ba-c", "-I 2a 2b" , 28}, // 341 + { 74, 74, "I m m a" , 0, "", "-I 2b 2" , 0 }, // 342 + { 74, 0, "I m m b" , 0, "ba-c", "-I 2a 2a" , 28}, // 343 + { 74, 0, "I b m m" , 0, "cab", "-I 2c 2c" , 1 }, // 344 + { 74, 0, "I c m m" , 0, "-cba", "-I 2 2b" , 3 }, // 345 + { 74, 0, "I m c m" , 0, "bca", "-I 2 2a" , 2 }, // 346 + { 74, 0, "I m a m" , 0, "a-cb", "-I 2c 2" , 5 }, // 347 + { 75, 75, "P 4" , 0, "", "P 4" , 0 }, // 348 + { 76, 76, "P 41" , 0, "", "P 4w" , 0 }, // 349 + { 77, 77, "P 42" , 0, "", "P 4c" , 0 }, // 350 + { 78, 78, "P 43" , 0, "", "P 4cw" , 0 }, // 351 + { 79, 79, "I 4" , 0, "", "I 4" , 0 }, // 352 + { 80, 80, "I 41" , 0, "", "I 4bw" , 0 }, // 353 + { 81, 81, "P -4" , 0, "", "P -4" , 0 }, // 354 + { 82, 82, "I -4" , 0, "", "I -4" , 0 }, // 355 + { 83, 83, "P 4/m" , 0, "", "-P 4" , 0 }, // 356 + { 84, 84, "P 42/m" , 0, "", "-P 4c" , 0 }, // 357 + { 85, 85, "P 4/n" , '1', "", "P 4ab -1ab" , 29}, // 358 + { 85, 0, "P 4/n" , '2', "", "-P 4a" , 0 }, // 359 + { 86, 86, "P 42/n" , '1', "", "P 4n -1n" , 30}, // 360 + { 86, 0, "P 42/n" , '2', "", "-P 4bc" , 0 }, // 361 + { 87, 87, "I 4/m" , 0, "", "-I 4" , 0 }, // 362 + { 88, 88, "I 41/a" , '1', "", "I 4bw -1bw" , 31}, // 363 + { 88, 0, "I 41/a" , '2', "", "-I 4ad" , 0 }, // 364 + { 89, 89, "P 4 2 2" , 0, "", "P 4 2" , 0 }, // 365 + { 90, 90, "P 4 21 2" , 0, "", "P 4ab 2ab" , 0 }, // 366 + { 91, 91, "P 41 2 2" , 0, "", "P 4w 2c" , 0 }, // 367 + { 92, 92, "P 41 21 2" , 0, "", "P 4abw 2nw" , 0 }, // 368 + { 93, 93, "P 42 2 2" , 0, "", "P 4c 2" , 0 }, // 369 + { 94, 94, "P 42 21 2" , 0, "", "P 4n 2n" , 0 }, // 370 + { 95, 95, "P 43 2 2" , 0, "", "P 4cw 2c" , 0 }, // 371 + { 96, 96, "P 43 21 2" , 0, "", "P 4nw 2abw" , 0 }, // 372 + { 97, 97, "I 4 2 2" , 0, "", "I 4 2" , 0 }, // 373 + { 98, 98, "I 41 2 2" , 0, "", "I 4bw 2bw" , 0 }, // 374 + { 99, 99, "P 4 m m" , 0, "", "P 4 -2" , 0 }, // 375 + {100, 100, "P 4 b m" , 0, "", "P 4 -2ab" , 0 }, // 376 + {101, 101, "P 42 c m" , 0, "", "P 4c -2c" , 0 }, // 377 + {102, 102, "P 42 n m" , 0, "", "P 4n -2n" , 0 }, // 378 + {103, 103, "P 4 c c" , 0, "", "P 4 -2c" , 0 }, // 379 + {104, 104, "P 4 n c" , 0, "", "P 4 -2n" , 0 }, // 380 + {105, 105, "P 42 m c" , 0, "", "P 4c -2" , 0 }, // 381 + {106, 106, "P 42 b c" , 0, "", "P 4c -2ab" , 0 }, // 382 + {107, 107, "I 4 m m" , 0, "", "I 4 -2" , 0 }, // 383 + {108, 108, "I 4 c m" , 0, "", "I 4 -2c" , 0 }, // 384 + {109, 109, "I 41 m d" , 0, "", "I 4bw -2" , 0 }, // 385 + {110, 110, "I 41 c d" , 0, "", "I 4bw -2c" , 0 }, // 386 + {111, 111, "P -4 2 m" , 0, "", "P -4 2" , 0 }, // 387 + {112, 112, "P -4 2 c" , 0, "", "P -4 2c" , 0 }, // 388 + {113, 113, "P -4 21 m" , 0, "", "P -4 2ab" , 0 }, // 389 + {114, 114, "P -4 21 c" , 0, "", "P -4 2n" , 0 }, // 390 + {115, 115, "P -4 m 2" , 0, "", "P -4 -2" , 0 }, // 391 + {116, 116, "P -4 c 2" , 0, "", "P -4 -2c" , 0 }, // 392 + {117, 117, "P -4 b 2" , 0, "", "P -4 -2ab" , 0 }, // 393 + {118, 118, "P -4 n 2" , 0, "", "P -4 -2n" , 0 }, // 394 + {119, 119, "I -4 m 2" , 0, "", "I -4 -2" , 0 }, // 395 + {120, 120, "I -4 c 2" , 0, "", "I -4 -2c" , 0 }, // 396 + {121, 121, "I -4 2 m" , 0, "", "I -4 2" , 0 }, // 397 + {122, 122, "I -4 2 d" , 0, "", "I -4 2bw" , 0 }, // 398 + {123, 123, "P 4/m m m" , 0, "", "-P 4 2" , 0 }, // 399 + {124, 124, "P 4/m c c" , 0, "", "-P 4 2c" , 0 }, // 400 + {125, 125, "P 4/n b m" , '1', "", "P 4 2 -1ab" , 21}, // 401 + {125, 0, "P 4/n b m" , '2', "", "-P 4a 2b" , 0 }, // 402 + {126, 126, "P 4/n n c" , '1', "", "P 4 2 -1n" , 20}, // 403 + {126, 0, "P 4/n n c" , '2', "", "-P 4a 2bc" , 0 }, // 404 + {127, 127, "P 4/m b m" , 0, "", "-P 4 2ab" , 0 }, // 405 + {128, 128, "P 4/m n c" , 0, "", "-P 4 2n" , 0 }, // 406 + {129, 129, "P 4/n m m" , '1', "", "P 4ab 2ab -1ab", 29}, // 407 + {129, 0, "P 4/n m m" , '2', "", "-P 4a 2a" , 0 }, // 408 + {130, 130, "P 4/n c c" , '1', "", "P 4ab 2n -1ab" , 29}, // 409 + {130, 0, "P 4/n c c" , '2', "", "-P 4a 2ac" , 0 }, // 410 + {131, 131, "P 42/m m c", 0, "", "-P 4c 2" , 0 }, // 411 + {132, 132, "P 42/m c m", 0, "", "-P 4c 2c" , 0 }, // 412 + {133, 133, "P 42/n b c", '1', "", "P 4n 2c -1n" , 32}, // 413 + {133, 0, "P 42/n b c", '2', "", "-P 4ac 2b" , 0 }, // 414 + {134, 134, "P 42/n n m", '1', "", "P 4n 2 -1n" , 33}, // 415 + {134, 0, "P 42/n n m", '2', "", "-P 4ac 2bc" , 0 }, // 416 + {135, 135, "P 42/m b c", 0, "", "-P 4c 2ab" , 0 }, // 417 + {136, 136, "P 42/m n m", 0, "", "-P 4n 2n" , 0 }, // 418 + {137, 137, "P 42/n m c", '1', "", "P 4n 2n -1n" , 32}, // 419 + {137, 0, "P 42/n m c", '2', "", "-P 4ac 2a" , 0 }, // 420 + {138, 138, "P 42/n c m", '1', "", "P 4n 2ab -1n" , 33}, // 421 + {138, 0, "P 42/n c m", '2', "", "-P 4ac 2ac" , 0 }, // 422 + {139, 139, "I 4/m m m" , 0, "", "-I 4 2" , 0 }, // 423 + {140, 140, "I 4/m c m" , 0, "", "-I 4 2c" , 0 }, // 424 + {141, 141, "I 41/a m d", '1', "", "I 4bw 2bw -1bw", 34}, // 425 + {141, 0, "I 41/a m d", '2', "", "-I 4bd 2" , 0 }, // 426 + {142, 142, "I 41/a c d", '1', "", "I 4bw 2aw -1bw", 35}, // 427 + {142, 0, "I 41/a c d", '2', "", "-I 4bd 2c" , 0 }, // 428 + {143, 143, "P 3" , 0, "", "P 3" , 0 }, // 429 + {144, 144, "P 31" , 0, "", "P 31" , 0 }, // 430 + {145, 145, "P 32" , 0, "", "P 32" , 0 }, // 431 + {146, 146, "R 3" , 'H', "", "R 3" , 0 }, // 432 + {146, 1146, "R 3" , 'R', "", "P 3*" , 36}, // 433 + {147, 147, "P -3" , 0, "", "-P 3" , 0 }, // 434 + {148, 148, "R -3" , 'H', "", "-R 3" , 0 }, // 435 + {148, 1148, "R -3" , 'R', "", "-P 3*" , 36}, // 436 + {149, 149, "P 3 1 2" , 0, "", "P 3 2" , 0 }, // 437 + {150, 150, "P 3 2 1" , 0, "", "P 3 2\"" , 0 }, // 438 + {151, 151, "P 31 1 2" , 0, "", "P 31 2 (0 0 4)", 0 }, // 439 + {152, 152, "P 31 2 1" , 0, "", "P 31 2\"" , 0 }, // 440 + {153, 153, "P 32 1 2" , 0, "", "P 32 2 (0 0 2)", 0 }, // 441 + {154, 154, "P 32 2 1" , 0, "", "P 32 2\"" , 0 }, // 442 + {155, 155, "R 3 2" , 'H', "", "R 3 2\"" , 0 }, // 443 + {155, 1155, "R 3 2" , 'R', "", "P 3* 2" , 36}, // 444 + {156, 156, "P 3 m 1" , 0, "", "P 3 -2\"" , 0 }, // 445 + {157, 157, "P 3 1 m" , 0, "", "P 3 -2" , 0 }, // 446 + {158, 158, "P 3 c 1" , 0, "", "P 3 -2\"c" , 0 }, // 447 + {159, 159, "P 3 1 c" , 0, "", "P 3 -2c" , 0 }, // 448 + {160, 160, "R 3 m" , 'H', "", "R 3 -2\"" , 0 }, // 449 + {160, 1160, "R 3 m" , 'R', "", "P 3* -2" , 36}, // 450 + {161, 161, "R 3 c" , 'H', "", "R 3 -2\"c" , 0 }, // 451 + {161, 1161, "R 3 c" , 'R', "", "P 3* -2n" , 36}, // 452 + {162, 162, "P -3 1 m" , 0, "", "-P 3 2" , 0 }, // 453 + {163, 163, "P -3 1 c" , 0, "", "-P 3 2c" , 0 }, // 454 + {164, 164, "P -3 m 1" , 0, "", "-P 3 2\"" , 0 }, // 455 + {165, 165, "P -3 c 1" , 0, "", "-P 3 2\"c" , 0 }, // 456 + {166, 166, "R -3 m" , 'H', "", "-R 3 2\"" , 0 }, // 457 + {166, 1166, "R -3 m" , 'R', "", "-P 3* 2" , 36}, // 458 + {167, 167, "R -3 c" , 'H', "", "-R 3 2\"c" , 0 }, // 459 + {167, 1167, "R -3 c" , 'R', "", "-P 3* 2n" , 36}, // 460 + {168, 168, "P 6" , 0, "", "P 6" , 0 }, // 461 + {169, 169, "P 61" , 0, "", "P 61" , 0 }, // 462 + {170, 170, "P 65" , 0, "", "P 65" , 0 }, // 463 + {171, 171, "P 62" , 0, "", "P 62" , 0 }, // 464 + {172, 172, "P 64" , 0, "", "P 64" , 0 }, // 465 + {173, 173, "P 63" , 0, "", "P 6c" , 0 }, // 466 + {174, 174, "P -6" , 0, "", "P -6" , 0 }, // 467 + {175, 175, "P 6/m" , 0, "", "-P 6" , 0 }, // 468 + {176, 176, "P 63/m" , 0, "", "-P 6c" , 0 }, // 469 + {177, 177, "P 6 2 2" , 0, "", "P 6 2" , 0 }, // 470 + {178, 178, "P 61 2 2" , 0, "", "P 61 2 (0 0 5)", 0 }, // 471 + {179, 179, "P 65 2 2" , 0, "", "P 65 2 (0 0 1)", 0 }, // 472 + {180, 180, "P 62 2 2" , 0, "", "P 62 2 (0 0 4)", 0 }, // 473 + {181, 181, "P 64 2 2" , 0, "", "P 64 2 (0 0 2)", 0 }, // 474 + {182, 182, "P 63 2 2" , 0, "", "P 6c 2c" , 0 }, // 475 + {183, 183, "P 6 m m" , 0, "", "P 6 -2" , 0 }, // 476 + {184, 184, "P 6 c c" , 0, "", "P 6 -2c" , 0 }, // 477 + {185, 185, "P 63 c m" , 0, "", "P 6c -2" , 0 }, // 478 + {186, 186, "P 63 m c" , 0, "", "P 6c -2c" , 0 }, // 479 + {187, 187, "P -6 m 2" , 0, "", "P -6 2" , 0 }, // 480 + {188, 188, "P -6 c 2" , 0, "", "P -6c 2" , 0 }, // 481 + {189, 189, "P -6 2 m" , 0, "", "P -6 -2" , 0 }, // 482 + {190, 190, "P -6 2 c" , 0, "", "P -6c -2c" , 0 }, // 483 + {191, 191, "P 6/m m m" , 0, "", "-P 6 2" , 0 }, // 484 + {192, 192, "P 6/m c c" , 0, "", "-P 6 2c" , 0 }, // 485 + {193, 193, "P 63/m c m", 0, "", "-P 6c 2" , 0 }, // 486 + {194, 194, "P 63/m m c", 0, "", "-P 6c 2c" , 0 }, // 487 + {195, 195, "P 2 3" , 0, "", "P 2 2 3" , 0 }, // 488 + {196, 196, "F 2 3" , 0, "", "F 2 2 3" , 0 }, // 489 + {197, 197, "I 2 3" , 0, "", "I 2 2 3" , 0 }, // 490 + {198, 198, "P 21 3" , 0, "", "P 2ac 2ab 3" , 0 }, // 491 + {199, 199, "I 21 3" , 0, "", "I 2b 2c 3" , 0 }, // 492 + {200, 200, "P m -3" , 0, "", "-P 2 2 3" , 0 }, // 493 + {201, 201, "P n -3" , '1', "", "P 2 2 3 -1n" , 20}, // 494 + {201, 0, "P n -3" , '2', "", "-P 2ab 2bc 3" , 0 }, // 495 + {202, 202, "F m -3" , 0, "", "-F 2 2 3" , 0 }, // 496 + {203, 203, "F d -3" , '1', "", "F 2 2 3 -1d" , 27}, // 497 + {203, 0, "F d -3" , '2', "", "-F 2uv 2vw 3" , 0 }, // 498 + {204, 204, "I m -3" , 0, "", "-I 2 2 3" , 0 }, // 499 + {205, 205, "P a -3" , 0, "", "-P 2ac 2ab 3" , 0 }, // 500 + {206, 206, "I a -3" , 0, "", "-I 2b 2c 3" , 0 }, // 501 + {207, 207, "P 4 3 2" , 0, "", "P 4 2 3" , 0 }, // 502 + {208, 208, "P 42 3 2" , 0, "", "P 4n 2 3" , 0 }, // 503 + {209, 209, "F 4 3 2" , 0, "", "F 4 2 3" , 0 }, // 504 + {210, 210, "F 41 3 2" , 0, "", "F 4d 2 3" , 0 }, // 505 + {211, 211, "I 4 3 2" , 0, "", "I 4 2 3" , 0 }, // 506 + {212, 212, "P 43 3 2" , 0, "", "P 4acd 2ab 3" , 0 }, // 507 + {213, 213, "P 41 3 2" , 0, "", "P 4bd 2ab 3" , 0 }, // 508 + {214, 214, "I 41 3 2" , 0, "", "I 4bd 2c 3" , 0 }, // 509 + {215, 215, "P -4 3 m" , 0, "", "P -4 2 3" , 0 }, // 510 + {216, 216, "F -4 3 m" , 0, "", "F -4 2 3" , 0 }, // 511 + {217, 217, "I -4 3 m" , 0, "", "I -4 2 3" , 0 }, // 512 + {218, 218, "P -4 3 n" , 0, "", "P -4n 2 3" , 0 }, // 513 + {219, 219, "F -4 3 c" , 0, "", "F -4a 2 3" , 0 }, // 514 + {220, 220, "I -4 3 d" , 0, "", "I -4bd 2c 3" , 0 }, // 515 + {221, 221, "P m -3 m" , 0, "", "-P 4 2 3" , 0 }, // 516 + {222, 222, "P n -3 n" , '1', "", "P 4 2 3 -1n" , 20}, // 517 + {222, 0, "P n -3 n" , '2', "", "-P 4a 2bc 3" , 0 }, // 518 + {223, 223, "P m -3 n" , 0, "", "-P 4n 2 3" , 0 }, // 519 + {224, 224, "P n -3 m" , '1', "", "P 4n 2 3 -1n" , 30}, // 520 + {224, 0, "P n -3 m" , '2', "", "-P 4bc 2bc 3" , 0 }, // 521 + {225, 225, "F m -3 m" , 0, "", "-F 4 2 3" , 0 }, // 522 + {226, 226, "F m -3 c" , 0, "", "-F 4a 2 3" , 0 }, // 523 + {227, 227, "F d -3 m" , '1', "", "F 4d 2 3 -1d" , 27}, // 524 + {227, 0, "F d -3 m" , '2', "", "-F 4vw 2vw 3" , 0 }, // 525 + {228, 228, "F d -3 c" , '1', "", "F 4d 2 3 -1ad" , 37}, // 526 + {228, 0, "F d -3 c" , '2', "", "-F 4ud 2vw 3" , 0 }, // 527 + {229, 229, "I m -3 m" , 0, "", "-I 4 2 3" , 0 }, // 528 + {230, 230, "I a -3 d" , 0, "", "-I 4bd 2c 3" , 0 }, // 529 + // And extra entries from syminfo.lib + { 5, 5005, "I 1 21 1" , 0, "b4", "I 2yb" , 38}, // 530 + { 5, 3005, "C 1 21 1" , 0, "b5", "C 2yb" , 14}, // 531 + { 18, 1018, "P 21212(a)", 0, "", "P 2ab 2a" , 14}, // 532 + { 20, 1020, "C 2 2 21a)", 0, "", "C 2ac 2" , 39}, // 533 + { 21, 1021, "C 2 2 2a" , 0, "", "C 2ab 2b" , 14}, // 534 + { 22, 1022, "F 2 2 2a" , 0, "", "F 2 2c" , 40}, // 535 + { 23, 1023, "I 2 2 2a" , 0, "", "I 2ab 2bc" , 33}, // 536 + { 94, 1094, "P 42 21 2a", 0, "", "P 4bc 2a" , 20}, // 537 + {197, 1197, "I 2 3a" , 0, "", "I 2ab 2bc 3" , 30}, // 538 + // And extra entries from Crystallographic Space Group Diagrams and Tables + // http://img.chem.ucl.ac.uk/sgp/ + // We want to have all entries from Open Babel and PDB. + // If available, Hall symbols are taken from + // https://cci.lbl.gov/cctbx/multiple_cell.html + // triclinic - enlarged unit cells + { 1, 0, "A 1" , 0, "", "A 1" , 41}, // 539 + { 1, 0, "B 1" , 0, "", "B 1" , 42}, // 540 + { 1, 0, "C 1" , 0, "", "C 1" , 43}, // 541 + { 1, 0, "F 1" , 0, "", "F 1" , 44}, // 542 + { 1, 0, "I 1" , 0, "", "I 1" , 45}, // 543 + { 2, 0, "A -1" , 0, "", "-A 1" , 41}, // 544 + { 2, 0, "B -1" , 0, "", "-B 1" , 42}, // 545 + { 2, 0, "C -1" , 0, "", "-C 1" , 43}, // 546 + { 2, 0, "F -1" , 0, "", "-F 1" , 44}, // 547 + { 2, 0, "I -1" , 0, "", "-I 1" , 45}, // 548 + // monoclinic (qualifiers such as "b1" are assigned arbitrary unique numbers) + { 3, 0, "B 1 2 1" , 0, "b1", "B 2y" , 46}, // 549 + { 3, 0, "C 1 1 2" , 0, "c1", "C 2" , 47}, // 550 + { 4, 0, "B 1 21 1" , 0, "b1", "B 2yb" , 46}, // 551 + { 4, 0, "C 1 1 21" , 0, "c2", "C 2c" , 47}, // 552 + { 5, 0, "F 1 2 1" , 0, "b6", "F 2y" , 48}, // 553 + { 8, 0, "F 1 m 1" , 0, "b4", "F -2y" , 48}, // 554 + { 9, 0, "F 1 d 1" , 0, "b4", "F -2yuw" , 49}, // 555 + { 12, 0, "F 1 2/m 1" , 0, "b4", "-F 2y" , 48}, // 556 + // orthorhombic + { 64, 0, "A b a m" , 0, "", "-A 2 2ab" , 3 }, // 557 (==306) + // tetragonal - enlarged C- and F-centred unit cells + { 89, 0, "C 4 2 2" , 0, "", "C 4 2" , 50}, // 558 + { 90, 0, "C 4 2 21" , 0, "", "C 4a 2" , 50}, // 559 + { 97, 0, "F 4 2 2" , 0, "", "F 4 2" , 50}, // 560 + {115, 0, "C -4 2 m" , 0, "", "C -4 2" , 50}, // 561 + {117, 0, "C -4 2 b" , 0, "", "C -4 2ya" , 50}, // 562 + {139, 0, "F 4/m m m" , 0, "", "-F 4 2" , 50}, // 563 +}; + +const SpaceGroupAltName spacegroup_tables::alt_names[28] = { + // In 1990's ITfC vol.A changed some of the standard names, introducing + // symbol 'e'. sgtbx interprets these new symbols with option ad_hoc_1992. + // spglib uses only the new symbols. + {"A e m 2", 0, 190}, // A b m 2 + {"B m e 2", 0, 191}, // B m a 2 + {"B 2 e m", 0, 192}, // B 2 c m + {"C 2 m e", 0, 193}, // C 2 m b + {"C m 2 e", 0, 194}, // C m 2 a + {"A e 2 m", 0, 195}, // A c 2 m + {"A e a 2", 0, 202}, // A b a 2 + {"B b e 2", 0, 203}, // B b a 2 + {"B 2 e b", 0, 204}, // B 2 c b + {"C 2 c e", 0, 205}, // C 2 c b + {"C c 2 e", 0, 206}, // C c 2 a + {"A e 2 a", 0, 207}, // A c 2 a + {"C m c e", 0, 303}, // C m c a + {"C c m e", 0, 304}, // C c m b + {"A e m a", 0, 305}, // A b m a + {"A e a m", 0, 306}, // A c a m + {"B b e m", 0, 307}, // B b c m + {"B m e b", 0, 308}, // B m a b + {"C m m e", 0, 315}, // C m m a + {"A e m m", 0, 317}, // A b m m + {"B m e m", 0, 319}, // B m c m + {"C c c e", '1', 321}, // C c c a + {"C c c e", '2', 322}, // C c c a + {"A e a a", '1', 325}, // A b a a + {"A e a a", '2', 326}, // A b a a + {"B b e b", '1', 329}, // B b c b + {"B b e b", '2', 330}, // B b c b + // help with parsing of unusual setting names that are present in the PDB + {"P 21 21 2a", 0, 532}, // P 21212(a) +}; + +// This table was generated by tools/gen_reciprocal_asu.py. +const unsigned char spacegroup_tables::ccp4_hkl_asu[230] = { + 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, + 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, + 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, + 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, + 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, + 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 5, 6, 7, 6, 7, 6, 7, 7, 7, + 6, 7, 6, 7, 7, 6, 6, 7, 7, 7, 7, 3, 3, 3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, + 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 9, 9, + 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9 +}; + +// Generated by tools/gen_sg_table.py. +const char* get_basisop(int basisop_idx) { + static const char* basisops[51] = { + "x,y,z", // 0 + "z,x,y", // 1 + "y,z,x", // 2 + "z,y,-x", // 3 + "x,y,-x+z", // 4 + "-x,z,y", // 5 + "-x+z,x,y", // 6 + "y,-x,z", // 7 + "y,-x+z,x", // 8 + "x-z,y,z", // 9 + "z,x-z,y", // 10 + "y,z,x-z", // 11 + "z,y,-x+z", // 12 + "x+z,y,-x", // 13 + "x+1/4,y+1/4,z", // 14 + "-x+z,z,y", // 15 + "-x,x+z,y", // 16 + "y,-x+z,z", // 17 + "y,-x,x+z", // 18 + "x+1/4,y-1/4,z", // 19 + "x-1/4,y-1/4,z-1/4", // 20 + "x-1/4,y-1/4,z", // 21 + "z,x-1/4,y-1/4", // 22 + "y-1/4,z,x-1/4", // 23 + "x-1/2,y-1/4,z+1/4", // 24 + "z+1/4,x-1/2,y-1/4", // 25 + "y-1/4,z+1/4,x-1/2", // 26 + "x+1/8,y+1/8,z+1/8", // 27 + "x+1/4,y-1/4,z+1/4", // 28 + "x-1/4,y+1/4,z", // 29 + "x+1/4,y+1/4,z+1/4", // 30 + "x,y+1/4,z+1/8", // 31 + "x-1/4,y+1/4,z+1/4", // 32 + "x-1/4,y+1/4,z-1/4", // 33 + "x-1/2,y+1/4,z+1/8", // 34 + "x-1/2,y+1/4,z-3/8", // 35 + "-y+z,x+z,-x+y+z", // 36 + "x-1/8,y-1/8,z-1/8", // 37 + "x+1/4,y+1/4,-x+z-1/4", // 38 + "x+1/4,y,z", // 39 + "x,y,z+1/4", // 40 + "-x,-y/2+z/2,y/2+z/2", // 41 + "-x/2+z/2,-y,x/2+z/2", // 42 + "x/2+y/2,x/2-y/2,-z", // 43 + "y/2+z/2,x/2+z/2,x/2+y/2", // 44 + "-x/2+y/2+z/2,x/2-y/2+z/2,x/2+y/2-z/2", // 45 + "x/2,y,-x/2+z", // 46 + "-x/2+z,x/2,y", // 47 + "x-z/2,y,z/2", // 48 + "x+z/2,y,z/2", // 49 + "x/2+y/2,-x/2+y/2,z", // 50 + }; + return basisops[basisop_idx]; +} + +const SpaceGroup* find_spacegroup_by_name(std::string name, double alpha, double gamma, + const char* prefer) { + bool prefer_2 = false; + bool prefer_R = false; + if (prefer) + for (const char* p = prefer; *p != '\0'; ++p) { + if (*p == '2') + prefer_2 = true; + else if (*p == 'R') + prefer_R = true; + else if (*p != '1' && *p != 'H') + throw std::invalid_argument("find_spacegroup_by_name(): invalid arg 'prefer'"); + } + const char* p = skip_space(name.c_str()); + if (*p >= '0' && *p <= '9') { // handle numbers + char *endptr; + long n = std::strtol(p, &endptr, 10); + return *endptr == '\0' ? find_spacegroup_by_number(n) : nullptr; + } + char first = *p & ~0x20; // to uppercase + if (first == '\0') + return nullptr; + if (first == 'H') + first = 'R'; + p = skip_space(p+1); + size_t start = p - name.c_str(); + // change letters to lower case, except the letter after : + for (size_t i = start; i < name.size(); ++i) { + if (name[i] >= 'A' && name[i] <= 'Z') + name[i] |= 0x20; // to lowercase + else if (name[i] == ':') + while (++i < name.size()) + if (name[i] >= 'a' && name[i] <= 'z') + name[i] &= ~0x20; // to uppercase + } + // allow names ending with R or H, such as R3R instead of R3:R + if (name.back() == 'h' || name.back() == 'r') { + name.back() &= ~0x20; // to uppercase + name.insert(name.end() - 1, ':'); + } + // The string that const char* p points to was just modified. + // This confuses some compilers (GCC 4.8), so let's re-assign p. + p = name.c_str() + start; + + for (const SpaceGroup& sg : spacegroup_tables::main) + if (sg.hm[0] == first) { + if (sg.hm[2] == *p) { + const char* a = skip_space(p + 1); + const char* b = skip_space(sg.hm + 3); + // In IT 1935 and 1952, symbols of centrosymmetric, cubic space groups + // 200-206 and 221-230 had symbol 3 (not -3), e.g. Pm3 instead of Pm-3, + // as listed in Table 3.3.3.1 in ITfC (2016) vol. A, p.788. + while ((*a == *b && *b != '\0') || + (*a == '3' && *b == '-' && b == sg.hm + 4 && *++b == '3')) { + a = skip_space(a+1); + b = skip_space(b+1); + } + if (*b == '\0') { + if (*a == '\0') { + // Change hexagonal settings to rhombohedral if the unit cell + // angles are more consistent with the latter. + // We have possible ambiguity in the hexagonal crystal family. + // For instance, "R 3" may mean "R 3:H" (hexagonal setting) or + // "R 3:R" (rhombohedral setting). The :H symbols come first + // in the table and are used by default. The ratio gamma:alpha + // is 120:90 in the hexagonal system and 1:1 in rhombohedral. + // We assume that the 'R' entry follows directly the 'H' entry. + if (sg.ext == 'H' && (alpha == 0. ? prefer_R : gamma < 1.125 * alpha)) + return &sg + 1; + // Similarly, the origin choice #2 follows directly #1. + if (sg.ext == '1' && prefer_2) + return &sg + 1; + return &sg; + } + if (*a == ':' && *skip_space(a+1) == sg.ext) + return &sg; + } + } else if (sg.hm[2] == '1' && sg.hm[3] == ' ') { + // check monoclinic short names, matching P2 to "P 1 2 1"; + // as an exception "B 2" == "B 1 1 2" (like in the PDB) + const char* b = sg.hm + 4; + if (*b != '1' || (first == 'B' && *++b == ' ' && *++b != '1')) { + char end = (b == sg.hm + 4 ? ' ' : '\0'); + const char* a = skip_space(p); + while (*a == *b && *b != end) { + ++a; + ++b; + } + if (*skip_space(a) == '\0' && *b == end) + return &sg; + } + } + } + for (const SpaceGroupAltName& sg : spacegroup_tables::alt_names) + if (sg.hm[0] == first && sg.hm[2] == *p) { + const char* a = skip_space(p + 1); + const char* b = skip_space(sg.hm + 3); + while (*a == *b && *b != '\0') { + a = skip_space(a+1); + b = skip_space(b+1); + } + if (*b == '\0' && + (*a == '\0' || (*a == ':' && *skip_space(a+1) == sg.ext))) + return &spacegroup_tables::main[sg.pos]; + } + return nullptr; +} + +} // namespace gemmi +