1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
|
/**
* @file Tools for manipulating SymOp bonds
*
* (c) Schrodinger, Inc.
*/
#include "SymOpTools.h"
#include "CoordSet.h"
#include "SymOp.h"
#include "pymol/algorithm.h"
#include <algorithm>
namespace pymol
{
/**
* Find all symmetry operations of a PBC (periodic boundary condition) bond.
*
* Guesses the correct symmetry operations based on the distance to
* all possible symmetry mates. This will include all images that are within 5%
* of the shortest distance (if tolerance is 1.05).
*
* @param cs Coordinate set which provides atom coordinates and symmetry
* information
* @param atm1 Atom index "from"
* @param atm2 Atom index "to"
* @param tolerance Distance tolerance (>= 1.0)
* @return List of symmetry operations and distances for atm2, given that atm1
* has symmetry operation 1_555.
*/
std::vector<std::pair<SymOp, float>> find_bond_symops(
CoordSet const& cs, unsigned atm1, unsigned atm2, float tolerance)
{
assert(tolerance >= 1);
auto const* sym = cs.getSymmetry();
if (!sym) {
return {};
}
int const idx[] = {
cs.atmToIdx(atm1),
cs.atmToIdx(atm2),
};
if (idx[0] == -1 || idx[1] == -1) {
return {};
}
auto const n_symmat = sym->getNSymMat();
float const* const v_fr = cs.coordPtr(idx[0]);
float v_to_buf[3];
float dist_min = FLT_MAX;
float dist_min_tol = FLT_MAX;
std::vector<std::pair<SymOp, float>> symops;
for (auto symop = SymOp("1_444"); symop.x < 2; ++symop.x) {
for (symop.y = -1; symop.y < 2; ++symop.y) {
for (symop.z = -1; symop.z < 2; ++symop.z) {
for (symop.index = 0; symop.index < n_symmat; ++symop.index) {
// atom which binds to its own symmetry mate must not be 1_555
if (idx[0] == idx[1] && !symop) {
continue;
}
auto const v_to = cs.coordPtrSym(idx[1], symop, v_to_buf);
assert(v_to);
auto const dist = diffsq3f(v_fr, v_to);
if (dist_min_tol < dist) {
continue;
}
if (dist_min > dist) {
dist_min = dist;
dist_min_tol = dist_min * tolerance * tolerance;
using cr_t = decltype(symops)::const_reference;
pymol::erase_if(symops,
[dist_min_tol](cr_t sd) { return dist_min_tol < sd.second; });
}
symops.emplace_back(symop, dist);
}
}
}
}
return symops;
}
} // namespace pymol
|
