| 35 | }; |
| 36 | |
| 37 | ExtendedSystem extend_system( |
| 38 | py::array_t<double> positions, |
| 39 | py::array_t<int> atomic_numbers, |
| 40 | py::array_t<double> cell, |
| 41 | py::array_t<bool> pbc, |
| 42 | double cutoff) |
| 43 | { |
| 44 | if (cutoff < 0) { |
| 45 | throw invalid_argument("Cutoff must be positive."); |
| 46 | } |
| 47 | // Determine the upper limit of how many copies we need in each cell vector |
| 48 | // direction. We take as many copies as needed to reach the radial cutoff. |
| 49 | // Notice that we need to use vectors that are perpendicular to the cell |
| 50 | // vectors to ensure that the correct atoms are included for non-cubic |
| 51 | // cells. |
| 52 | auto positions_u = positions.unchecked<2>(); |
| 53 | auto atomic_numbers_u = atomic_numbers.unchecked<1>(); |
| 54 | auto cell_u = cell.unchecked<2>(); |
| 55 | auto pbc_u = pbc.unchecked<1>(); |
| 56 | vector<double> a = {cell_u(0, 0), cell_u(0, 1), cell_u(0, 2)}; |
| 57 | vector<double> b = {cell_u(1, 0), cell_u(1, 1), cell_u(1, 2)}; |
| 58 | vector<double> c = {cell_u(2, 0), cell_u(2, 1), cell_u(2, 2)}; |
| 59 | |
| 60 | vector<double> p1 = cross(b, c); |
| 61 | vector<double> p2 = cross(c, a); |
| 62 | vector<double> p3 = cross(a, b); |
| 63 | |
| 64 | // Projections of basis vectors onto perpendicular vectors. |
| 65 | double p1_coeff = dot(a, p1) / dot(p1, p1); |
| 66 | double p2_coeff = dot(b, p2) / dot(p2, p2); |
| 67 | double p3_coeff = dot(c, p3) / dot(p3, p3); |
| 68 | for(double &x : p1) { x *= p1_coeff; } |
| 69 | for(double &x : p2) { x *= p2_coeff; } |
| 70 | for(double &x : p3) { x *= p3_coeff; } |
| 71 | vector<vector<double>> vectors = {p1, p2, p3}; |
| 72 | |
| 73 | // Figure out how many copies to take per basis vector. Determined by how |
| 74 | // many perpendicular projections fit into the cutoff distance. |
| 75 | vector<int> n_copies_axis(3); |
| 76 | vector<vector<int>> multipliers; |
| 77 | for (int i=0; i < 3; ++i) { |
| 78 | if (pbc_u(i)) { |
| 79 | double length = norm(vectors[i]); |
| 80 | double factor = cutoff/length; |
| 81 | int multiplier = (int)ceil(factor); |
| 82 | n_copies_axis[i] = multiplier; |
| 83 | |
| 84 | // Store multipliers explicitly into a list in an order that keeps the |
| 85 | // original system in the same place both in space and in index. |
| 86 | vector<int> multiples; |
| 87 | for (int j=0; j < multiplier + 1; ++j) { |
| 88 | multiples.push_back(j); |
| 89 | } |
| 90 | for (int j=-multiplier; j < 0; ++j) { |
| 91 | multiples.push_back(j); |
| 92 | } |
| 93 | multipliers.push_back(multiples); |
| 94 | } else { |
no test coverage detected