From b06dfc8357e9706d96f55eb0991e48aaeabdcae2 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Wed, 13 May 2026 13:13:17 +0200 Subject: [PATCH] Gemmi: Include through FetchContent full gemmi library (not limited cpp/hpp files) --- 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, 12 insertions(+), 4214 deletions(-) delete mode 100644 symmetry/CMakeLists.txt delete mode 100644 symmetry/LICENSE.txt delete mode 100644 symmetry/gemmi/cellred.hpp delete mode 100644 symmetry/gemmi/fail.hpp delete mode 100644 symmetry/gemmi/math.hpp delete mode 100644 symmetry/gemmi/symmetry.hpp delete mode 100644 symmetry/gemmi/unitcell.hpp delete mode 100644 symmetry/symmetry.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index c49e0644..9beada6d 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -95,6 +95,13 @@ 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 @@ -142,7 +149,7 @@ FetchContent_Declare( EXCLUDE_FROM_ALL ) -FetchContent_MakeAvailable(zstd sls_detector_package catch2 hdf5 spdlog httplib) +FetchContent_MakeAvailable(zstd sls_detector_package catch2 hdf5 spdlog httplib gemmi) ADD_SUBDIRECTORY(jungfrau) ADD_SUBDIRECTORY(compression) @@ -153,7 +160,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 a15a652b..d9b1d8b2 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 Threads::Threads -lrt ) +TARGET_LINK_LIBRARIES(JFJochCommon JFJochLogger Compression JFCalibration gemmi_cpp Threads::Threads -lrt ) TARGET_LINK_LIBRARIES(JFJochZMQ "$") diff --git a/docs/SOFTWARE.md b/docs/SOFTWARE.md index d64671e8..bddc3952 100644 --- a/docs/SOFTWARE.md +++ b/docs/SOFTWARE.md @@ -33,6 +33,7 @@ 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: @@ -44,6 +45,5 @@ 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 394dbaed..6bb52518 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 gemmi) +TARGET_LINK_LIBRARIES(JFJochImageAnalysis JFJochAzIntEngine JFJochImagePreprocessing JFJochBraggPrediction JFJochBraggIntegration JFJochLatticeSearch JFJochIndexing JFJochSpotFinding JFJochCommon JFJochGeomRefinement JFJochScaleMerge) diff --git a/symmetry/CMakeLists.txt b/symmetry/CMakeLists.txt deleted file mode 100644 index fed3f792..00000000 --- a/symmetry/CMakeLists.txt +++ /dev/null @@ -1,2 +0,0 @@ -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 deleted file mode 100644 index 14e2f777..00000000 --- a/symmetry/LICENSE.txt +++ /dev/null @@ -1,373 +0,0 @@ -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 deleted file mode 100644 index 7174f776..00000000 --- a/symmetry/gemmi/cellred.hpp +++ /dev/null @@ -1,406 +0,0 @@ -// 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 deleted file mode 100644 index 10596385..00000000 --- a/symmetry/gemmi/fail.hpp +++ /dev/null @@ -1,93 +0,0 @@ -// 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 deleted file mode 100644 index 87c10516..00000000 --- a/symmetry/gemmi/math.hpp +++ /dev/null @@ -1,458 +0,0 @@ -// 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 deleted file mode 100644 index 203324fd..00000000 --- a/symmetry/gemmi/symmetry.hpp +++ /dev/null @@ -1,1044 +0,0 @@ -// 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 deleted file mode 100644 index 25bb8b46..00000000 --- a/symmetry/gemmi/unitcell.hpp +++ /dev/null @@ -1,618 +0,0 @@ -// 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 deleted file mode 100644 index 05b5d1c0..00000000 --- a/symmetry/symmetry.cpp +++ /dev/null @@ -1,1215 +0,0 @@ -// 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 -