Loading...
Searching...
No Matches
lattice.hpp
Go to the documentation of this file.
1#pragma once
2
3/**
4 * Lattice vector generation from Bravais lattice index.
5 *
6 * Provides get_lattice_vectors_from_bravais_index, which converts a standard
7 * Bravais lattice index and cell dimension parameters into an explicit 3x3
8 * matrix of lattice vectors in units of the lattice parameter (alat).
9 */
10
11#include <cmath>
12#include <stdexcept>
13#include <string>
14
15#include "layer_0_core/types/constants.hpp"
16#include "layer_0_core/types/matrix_3x3.hpp"
17#include "layer_0_core/types/type_aliases.hpp"
18
19namespace q3 {
20
21/**
22 * Build lattice vectors from cell parameters a, b, c (any consistent length
23 * unit) and angles alpha, beta, gamma (degrees).
24 *
25 * Standard triclinic construction: a1 along x, a2 in the xy plane
26 * (a2 = b(cos γ, sin γ, 0)), a3 completing the cell. Rows of the returned
27 * matrix are a1, a2, a3 in the same unit as the input lengths.
28 *
29 * Throws std::runtime_error for degenerate angle combinations (sin γ ≈ 0 or a
30 * zero-volume cell).
31 */
32inline Matrix3x3 get_lattice_vectors_from_cell_parameters(Real length_a,
33 Real length_b,
34 Real length_c,
35 Real alpha_degrees,
36 Real beta_degrees,
37 Real gamma_degrees) {
38 Real cos_alpha = std::cos(alpha_degrees * constants::degrees_to_radians);
39 Real cos_beta = std::cos(beta_degrees * constants::degrees_to_radians);
40 Real cos_gamma = std::cos(gamma_degrees * constants::degrees_to_radians);
41 Real sin_gamma = std::sin(gamma_degrees * constants::degrees_to_radians);
42 if (std::abs(sin_gamma) < 1e-12)
43 throw std::runtime_error("Invalid lattice angle gamma: " + std::to_string(gamma_degrees));
44
45 Real c_direction_x = cos_beta;
46 Real c_direction_y = (cos_alpha - cos_beta * cos_gamma) / sin_gamma;
47 Real c_direction_z_squared =
48 1.0 - c_direction_x * c_direction_x - c_direction_y * c_direction_y;
49 if (c_direction_z_squared <= 0.0)
50 throw std::runtime_error("Invalid lattice angles: cell has no volume");
51
52 Matrix3x3 lattice_vectors;
53 lattice_vectors[0] = {length_a, 0.0, 0.0};
54 lattice_vectors[1] = {length_b * cos_gamma, length_b * sin_gamma, 0.0};
55 lattice_vectors[2] = {length_c * c_direction_x,
56 length_c * c_direction_y,
57 length_c * std::sqrt(c_direction_z_squared)};
58 return lattice_vectors;
59}
60
61/**
62 * Compute lattice vectors from a Bravais lattice index and cell dimension parameters.
63 *
64 * The index encoding, and the primitive vectors each index stands for, are the `ibrav`
65 * convention published in Quantum ESPRESSO's pw.x input documentation (Doc/INPUT_PW, the
66 * `ibrav` variable). They are reproduced because they are what a QE input file *means*:
67 * positions given in crystal units are coordinates in exactly these vectors, so a
68 * different -- equally valid -- primitive cell for the same Bravais lattice would read
69 * the same file as a different structure. There is no wider standard that fixes them;
70 * other codes number and orient their cells differently. The one-line comment on each
71 * `case` below restates that documentation's entry for its index, in the dimensionless
72 * form this function returns.
73 *
74 * The returned matrix has rows a1, a2, a3 in units of alat (cell_dimension_1).
75 *
76 * Supported bravais_lattice_index values: 0, 1, 2, 3, -3, 4, 5, -5, 6, 7, 8, 9, 10, 11, 12,
77 * -12, 14. The convention also defines 13, -13, -9 and 91; Q3 does not build those cells
78 * and refuses them rather than substituting a neighbouring one.
79 *
80 * Throws std::runtime_error for unsupported bravais_lattice_index values.
81 */
82inline Matrix3x3 get_lattice_vectors_from_bravais_index(int bravais_lattice_index,
83 Real cell_dimension_1,
84 Real cell_dimension_2 = 0,
85 Real cell_dimension_3 = 0,
86 Real cell_dimension_4 = 0,
87 Real cell_dimension_5 = 0,
88 Real cell_dimension_6 = 0) {
89 Matrix3x3 lattice_vectors;
90 Real ratio_c_over_a = cell_dimension_3;
91 Real ratio_b_over_a = cell_dimension_2;
92
93 switch (bravais_lattice_index) {
94 case 0:
95 /**
96 * Free lattice: caller must set lattice vectors separately.
97 * Initialize to identity as placeholder.
98 */
99 lattice_vectors[0] = {1.0, 0.0, 0.0};
100 lattice_vectors[1] = {0.0, 1.0, 0.0};
101 lattice_vectors[2] = {0.0, 0.0, 1.0};
102 break;
103
104 case 1:
105 // Simple cubic: a1=a(1,0,0), a2=a(0,1,0), a3=a(0,0,1)
106 lattice_vectors[0] = {1.0, 0.0, 0.0};
107 lattice_vectors[1] = {0.0, 1.0, 0.0};
108 lattice_vectors[2] = {0.0, 0.0, 1.0};
109 break;
110
111 case 2:
112 // FCC: a1 = a/2 (-1,0,1), a2 = a/2 (0,1,1), a3 = a/2 (-1,1,0)
113 lattice_vectors[0] = {-0.5, 0.0, 0.5};
114 lattice_vectors[1] = {0.0, 0.5, 0.5};
115 lattice_vectors[2] = {-0.5, 0.5, 0.0};
116 break;
117
118 case 3:
119 // BCC: a1 = a/2 (1,1,1), a2 = a/2 (-1,1,1), a3 = a/2 (-1,-1,1)
120 lattice_vectors[0] = {0.5, 0.5, 0.5};
121 lattice_vectors[1] = {-0.5, 0.5, 0.5};
122 lattice_vectors[2] = {-0.5, -0.5, 0.5};
123 break;
124
125 case -3:
126 // BCC alt: a1 = a/2 (-1,1,1), a2 = a/2 (1,-1,1), a3 = a/2 (1,1,-1)
127 lattice_vectors[0] = {-0.5, 0.5, 0.5};
128 lattice_vectors[1] = {0.5, -0.5, 0.5};
129 lattice_vectors[2] = {0.5, 0.5, -0.5};
130 break;
131
132 case 4: {
133 // Hexagonal: a1 = a(1,0,0), a2 = a(-1/2,√3/2,0), a3 = a(0,0,c/a)
134 Real half_sqrt3 = std::sqrt(3.0) / 2.0;
135 lattice_vectors[0] = {1.0, 0.0, 0.0};
136 lattice_vectors[1] = {-0.5, half_sqrt3, 0.0};
137 lattice_vectors[2] = {0.0, 0.0, ratio_c_over_a};
138 break;
139 }
140
141 case 5: {
142 // Trigonal R (3-fold axis c): cell_dimension_4=cos(alpha)
143 Real cos_alpha = cell_dimension_4;
144 Real vector_component_x = std::sqrt((1.0 - cos_alpha) / 2.0);
145 Real vector_component_y = std::sqrt((1.0 - cos_alpha) / 6.0);
146 Real vector_component_z = std::sqrt((1.0 + 2.0 * cos_alpha) / 3.0);
147 lattice_vectors[0] = {vector_component_x, -vector_component_y, vector_component_z};
148 lattice_vectors[1] = {0.0, 2.0 * vector_component_y, vector_component_z};
149 lattice_vectors[2] = {-vector_component_x, -vector_component_y, vector_component_z};
150 break;
151 }
152
153 case -5: {
154 // Trigonal R (3-fold axis <111>): cell_dimension_4=cos(alpha)
155 Real cos_alpha = cell_dimension_4;
156 Real vector_component_y = std::sqrt((1.0 - cos_alpha) / 6.0);
157 Real vector_component_z = std::sqrt((1.0 + 2.0 * cos_alpha) / 3.0);
158 Real vector_component_u =
159 vector_component_z - 2.0 * std::sqrt(2.0) * vector_component_y;
160 Real vector_component_v = vector_component_z + std::sqrt(2.0) * vector_component_y;
161 lattice_vectors[0] = {vector_component_u / std::sqrt(3.0),
162 vector_component_v / std::sqrt(3.0),
163 vector_component_v / std::sqrt(3.0)};
164 lattice_vectors[1] = {vector_component_v / std::sqrt(3.0),
165 vector_component_u / std::sqrt(3.0),
166 vector_component_v / std::sqrt(3.0)};
167 lattice_vectors[2] = {vector_component_v / std::sqrt(3.0),
168 vector_component_v / std::sqrt(3.0),
169 vector_component_u / std::sqrt(3.0)};
170 break;
171 }
172
173 case 6:
174 // Tetragonal P: a1 = a(1,0,0), a2 = a(0,1,0), a3 = a(0,0,c/a)
175 lattice_vectors[0] = {1.0, 0.0, 0.0};
176 lattice_vectors[1] = {0.0, 1.0, 0.0};
177 lattice_vectors[2] = {0.0, 0.0, ratio_c_over_a};
178 break;
179
180 case 7:
181 // Body-centered tetragonal: a1=a/2(1,-1,c/a), a2=a/2(1,1,c/a), a3=a/2(-1,-1,c/a)
182 lattice_vectors[0] = {0.5, -0.5, 0.5 * ratio_c_over_a};
183 lattice_vectors[1] = {0.5, 0.5, 0.5 * ratio_c_over_a};
184 lattice_vectors[2] = {-0.5, -0.5, 0.5 * ratio_c_over_a};
185 break;
186
187 case 8:
188 // Orthorhombic P: a1=a(1,0,0), a2=a(0,b/a,0), a3=a(0,0,c/a)
189 lattice_vectors[0] = {1.0, 0.0, 0.0};
190 lattice_vectors[1] = {0.0, ratio_b_over_a, 0.0};
191 lattice_vectors[2] = {0.0, 0.0, ratio_c_over_a};
192 break;
193
194 case 9:
195 /**
196 * Base-centered orthorhombic (bco):
197 * a1=a/2(1,b/a,0), a2=a/2(-1,b/a,0), a3=a(0,0,c/a)
198 */
199 lattice_vectors[0] = {0.5, 0.5 * ratio_b_over_a, 0.0};
200 lattice_vectors[1] = {-0.5, 0.5 * ratio_b_over_a, 0.0};
201 lattice_vectors[2] = {0.0, 0.0, ratio_c_over_a};
202 break;
203
204 case 10:
205 /**
206 * Face-centered orthorhombic:
207 * a1=a/2(1,0,c/a), a2=a/2(1,b/a,0), a3=a/2(0,b/a,c/a)
208 */
209 lattice_vectors[0] = {0.5, 0.0, 0.5 * ratio_c_over_a};
210 lattice_vectors[1] = {0.5, 0.5 * ratio_b_over_a, 0.0};
211 lattice_vectors[2] = {0.0, 0.5 * ratio_b_over_a, 0.5 * ratio_c_over_a};
212 break;
213
214 case 11:
215 /**
216 * Body-centered orthorhombic:
217 * a1=a/2(1,b/a,c/a), a2=a/2(-1,b/a,c/a), a3=a/2(-1,-b/a,c/a)
218 */
219 lattice_vectors[0] = {0.5, 0.5 * ratio_b_over_a, 0.5 * ratio_c_over_a};
220 lattice_vectors[1] = {-0.5, 0.5 * ratio_b_over_a, 0.5 * ratio_c_over_a};
221 lattice_vectors[2] = {-0.5, -0.5 * ratio_b_over_a, 0.5 * ratio_c_over_a};
222 break;
223
224 case 12: {
225 /**
226 * Monoclinic P (unique axis c): cell_dimension_4=cos(ab)
227 * a1=a(1,0,0), a2=a*b/a(cos(ab),sin(ab),0), a3=a(0,0,c/a)
228 */
229 Real cos_ab = cell_dimension_4;
230 Real sin_ab = std::sqrt(1.0 - cos_ab * cos_ab);
231 lattice_vectors[0] = {1.0, 0.0, 0.0};
232 lattice_vectors[1] = {ratio_b_over_a * cos_ab, ratio_b_over_a * sin_ab, 0.0};
233 lattice_vectors[2] = {0.0, 0.0, ratio_c_over_a};
234 break;
235 }
236
237 case -12: {
238 // Monoclinic P (unique axis b): cell_dimension_5=cos(ac)
239 Real cos_ac = cell_dimension_5;
240 Real sin_ac = std::sqrt(1.0 - cos_ac * cos_ac);
241 lattice_vectors[0] = {1.0, 0.0, 0.0};
242 lattice_vectors[1] = {0.0, ratio_b_over_a, 0.0};
243 lattice_vectors[2] = {ratio_c_over_a * cos_ac, 0.0, ratio_c_over_a * sin_ac};
244 break;
245 }
246
247 case 14: {
248 /**
249 * Triclinic: cell_dimension_4=cos(bc), cell_dimension_5=cos(ac),
250 * cell_dimension_6=cos(ab)
251 */
252 Real cos_bc = cell_dimension_4;
253 Real cos_ac = cell_dimension_5;
254 Real cos_ab = cell_dimension_6;
255 Real sin_ab = std::sqrt(1.0 - cos_ab * cos_ab);
256 lattice_vectors[0] = {1.0, 0.0, 0.0};
257 lattice_vectors[1] = {ratio_b_over_a * cos_ab, ratio_b_over_a * sin_ab, 0.0};
258 Real vector_component_c_x = ratio_c_over_a * cos_ac;
259 Real vector_component_c_y = ratio_c_over_a * (cos_bc - cos_ac * cos_ab) / sin_ab;
260 Real vector_component_c_z =
261 ratio_c_over_a *
262 std::sqrt(1.0 + 2.0 * cos_bc * cos_ac * cos_ab - cos_bc * cos_bc - cos_ac * cos_ac -
263 cos_ab * cos_ab) /
264 sin_ab;
265 lattice_vectors[2] = {vector_component_c_x, vector_component_c_y, vector_component_c_z};
266 break;
267 }
268
269 default:
270 throw std::runtime_error(
271 "Unsupported bravais_lattice_index=" + std::to_string(bravais_lattice_index) +
272 ". Supported: 0 (free), 1 (SC), 2 (FCC), 3/-3 (BCC), 4 (hex), 5/-5 (trig R), "
273 "6 (tet), 7 (bct), 8 (ort), 9 (bco), 10 (fco), 11 (bort), "
274 "12/-12 (mono), 14 (tri)");
275 }
276
277 return lattice_vectors;
278}
279
280} // namespace q3
matrix_3x3.hpp
q3::constants::degrees_to_radians
constexpr Real degrees_to_radians
Definition constants.hpp:29
q3::get_lattice_vectors_from_cell_parameters
Matrix3x3 get_lattice_vectors_from_cell_parameters(Real length_a, Real length_b, Real length_c, Real alpha_degrees, Real beta_degrees, Real gamma_degrees)
Definition lattice.hpp:32
q3::get_lattice_vectors_from_bravais_index
Matrix3x3 get_lattice_vectors_from_bravais_index(int bravais_lattice_index, Real cell_dimension_1, Real cell_dimension_2=0, Real cell_dimension_3=0, Real cell_dimension_4=0, Real cell_dimension_5=0, Real cell_dimension_6=0)
Definition lattice.hpp:82
q3::Real
double Real
Definition type_aliases.hpp:20
q3::Matrix3x3
3×3 matrix stored row-major: row i = mat[i]
Definition matrix_3x3.hpp:15
type_aliases.hpp