mx/Jungfraujoch
1
1
Fork 0
Code Issues Pull Requests 1 Actions 1 Packages Projects Releases 53 Wiki Activity

Gemmi: Include through FetchContent full gemmi library (not limited cpp/hpp files)
Build Packages / Unit tests (push) Failing after 6m56s
Details
Build Packages / build:rpm (rocky8_nocuda) (push) Failing after 7m55s
Details
Build Packages / build:rpm (rocky9_nocuda) (push) Failing after 8m37s
Details
Build Packages / build:rpm (ubuntu2404_nocuda) (push) Failing after 8m50s
Details
Build Packages / build:rpm (ubuntu2204_nocuda) (push) Failing after 9m1s
Details
Build Packages / build:rpm (rocky8_sls9) (push) Failing after 10m40s
Details
Build Packages / build:rpm (rocky9_sls9) (push) Failing after 11m29s
Details
Build Packages / build:rpm (rocky8) (push) Failing after 6m56s
Details
Build Packages / Generate python client (push) Successful in 1m33s
Details
Build Packages / build:rpm (rocky9) (push) Failing after 8m10s
Details
Build Packages / Create release (push) Skipped
Details
Build Packages / Build documentation (push) Successful in 1m44s
Details
Build Packages / build:rpm (ubuntu2204) (push) Failing after 8m33s
Details
Build Packages / DIALS test (push) Failing after 8m15s
Details
Build Packages / build:rpm (ubuntu2404) (push) Failing after 9m30s
Details
Build Packages / XDS test (neggia plugin) (push) Failing after 6m57s
Details
Build Packages / XDS test (JFJoch plugin) (push) Failing after 7m48s
Details
Build Packages / XDS test (durin plugin) (push) Failing after 7m57s
Details

Browse Source
This commit is contained in:
leonarski_f
2026-05-13 13:13:17 +02:00
parent d5b928fa73
commit b06dfc8357
12 changed files with 12 additions and 4214 deletions
Download Patch File Download Diff File Expand all files Collapse all files
CMakeLists.txt
+9 -2
View File
@@ -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)
common/CMakeLists.txt
+1 -1
View File
@@ -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 "$<BUILD_INTERFACE:libzmq-static>")
docs/SOFTWARE.md
+1 -1
View File
@@ -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
image_analysis/CMakeLists.txt
+1 -1
View File
@@ -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)
symmetry/CMakeLists.txt
-2
View File
@@ -1,2 +0,0 @@
ADD_LIBRARY(gemmi STATIC symmetry.cpp gemmi/symmetry.hpp gemmi/fail.hpp)
TARGET_INCLUDE_DIRECTORIES(gemmi PUBLIC .)
symmetry/LICENSE.txt
-373
View File
@@ -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.
symmetry/gemmi/cellred.hpp
-406
View File
@@ -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 <cmath>
#include <array>
#include <memory> // 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<Op> 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<double,6>& 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<double,6> parameters() const { return {A, B, C, xi, eta, zeta}; }
std::array<double,6> 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<double,6> s;
explicit SellingVector(const std::array<double,6>& 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<double,6> 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<double,6> 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
symmetry/gemmi/fail.hpp
-93
View File
@@ -1,93 +0,0 @@
// Copyright 2017 Global Phasing Ltd.
//
// fail(), unreachable() and __declspec/__attribute__ macros
#ifndef GEMMI_FAIL_HPP_
#define GEMMI_FAIL_HPP_
#include <cerrno> // for errno
#include <stdexcept> // for runtime_error
#include <system_error> // for system_error
#include <string>
#include <utility> // 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<typename T, typename... Args> [[noreturn]]
void fail(std::string&& str, T&& arg1, Args&&... args) {
str += arg1;
fail(std::move(str), std::forward<Args>(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
symmetry/gemmi/math.hpp
-458
View File
@@ -1,458 +0,0 @@
// Copyright 2018 Global Phasing Ltd.
//
// Math utilities. 3D linear algebra.
#ifndef GEMMI_MATH_HPP_
#define GEMMI_MATH_HPP_
#include <cmath> // for fabs, cos, sqrt, round
#include <algorithm> // for min
#include <array>
#include <stdexcept> // for out_of_range
#include <type_traits> // 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<int>(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<class T> constexpr T clamp(T v, T lo, T hi) {
return std::min(std::max(v, lo), hi);
}
template <typename Real>
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<int, 3> 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<Vec3_*>(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_<double>;
using Vec3f = Vec3_<float>;
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<typename T> 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<T, 6> elements_pdb() const { return {{u11, u22, u33, u12, u13, u23}}; }
std::array<T, 6> 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<typename Real>
SMat33<Real> scaled(Real s) const {
return SMat33<Real>{u11*s, u22*s, u33*s, u12*s, u13*s, u23*s};
}
// returns U + kI
SMat33<T> 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<typename VT>
auto r_u_r(const Vec3_<VT>& 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<int,3>& 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<typename Real=double>
SMat33<Real> 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<Real>(
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<Real>{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<double, 3> calculate_eigenvalues() const {
double p1 = u12*u12 + u13*u13 + u23*u23;
if (p1 == 0)
return {{u11, u22, u33}};
double q = (1./3.) * trace();
SMat33<double> 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<typename Pos>
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 T>
typename std::enable_if<std::is_integral<T>::value, bool>::type
is_nan(T) { return false; }
template<typename T>
typename std::enable_if<std::is_floating_point<T>::value, bool>::type
is_nan(T a) { return std::isnan(a); }
template<typename T>
typename std::enable_if<std::is_integral<T>::value, bool>::type
is_same(T a, T b) { return a == b; }
template<typename T>
typename std::enable_if<std::is_floating_point<T>::value, bool>::type
is_same(T a, T b) { return std::isnan(b) ? std::isnan(a) : a == b; }
} // namespace impl
} // namespace gemmi
#endif
symmetry/gemmi/symmetry.hpp
-1044
View File
File diff suppressed because it is too large Load Diff
symmetry/gemmi/unitcell.hpp
-618
View File
@@ -1,618 +0,0 @@
// Copyright 2017 Global Phasing Ltd.
//
// Unit cell.
#ifndef GEMMI_UNITCELL_HPP_
#define GEMMI_UNITCELL_HPP_
#include <cassert>
#include <cmath> // for cos, sin, sqrt, floor, NAN
#include <vector>
#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<int, 3>;
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<double,6>& 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<double, 6>& 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<FTransform> 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<double>& 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<double,6>& 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<double,6> metric = metric_tensor().elements_voigt();
for (const Op& op : gops.sym_ops) {
Mat33 m = orth.mat.multiply(rot_as_mat33(op));
std::array<double,6> 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<NcsOp>& 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<FTransform> get_ncs_transforms() const {
std::vector<FTransform> 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<Position> orthogonalize_box(const Box<Fractional>& f) const {
Box<Position> 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<int>(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<NearestImage> find_nearest_pbc_images(const Fractional& fref, double dist,
const Fractional& fpos, int image_idx) const {
std::vector<NearestImage> 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<double> 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<double> 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
symmetry/symmetry.cpp
-1215
View File
File diff suppressed because it is too large Load Diff