Loading...
Searching...
No Matches
seams::steinhardt Namespace Reference

Functions

void minImage (double &dx, double &dy, double &dz, double bx, double by, double bz)
double normLM (int orderL, int absM)
double legendreAmp (int orderL, int absM, const double *s, const double *c)
void ylmAll (int orderL, double theta, double phi, double *out)
void qlmAddBond (int orderL, double dx, double dy, double dz, double *qlmInterleaved, int row, int nComp, int &nUsed)
void qlmOneAtom (int iatom, int orderL, const double *xyz, const int *offsets, const int *cols, double bx, double by, double bz, double *qlmInterleaved)
void qlmOneAtomDr (int iatom, int orderL, const double *dr, const int *offsets, const int *cols, double *qlmInterleaved)
void qlOneAtom (int iatom, int orderL, const double *qlmInterleaved, const int *offsets, const int *cols, double *ql, double *qlBar)

Function Documentation

◆ legendreAmp()

double seams::steinhardt::legendreAmp ( int orderL,
int absM,
const double * s,
const double * c )
inline

Definition at line 122 of file steinhardt_device.hpp.

123 {
124 if (orderL == 3) {
125 switch (absM) {
126 case 0:
127 return 5.0 * c[3] - 3.0 * c[1];
128 case 1:
129 return s[1] * (5.0 * c[2] - 1.0);
130 case 2:
131 return s[2] * c[1];
132 case 3:
133 return s[3];
134 default:
135 return 0.0;
136 }
137 }
138 if (orderL == 4) {
139 switch (absM) {
140 case 0:
141 return 35.0 * c[4] - 30.0 * c[2] + 3.0;
142 case 1:
143 return s[1] * (7.0 * c[3] - 3.0 * c[1]);
144 case 2:
145 return s[2] * (7.0 * c[2] - 1.0);
146 case 3:
147 return s[3] * c[1];
148 case 4:
149 return s[4];
150 default:
151 return 0.0;
152 }
153 }
154 if (orderL == 6) {
155 switch (absM) {
156 case 0:
157 return 231.0 * c[6] - 315.0 * c[4] + 105.0 * c[2] - 5.0;
158 case 1:
159 return s[1] * (33.0 * c[5] - 30.0 * c[3] + 5.0 * c[1]);
160 case 2:
161 return s[2] * (33.0 * c[4] - 18.0 * c[2] + 1.0);
162 case 3:
163 return s[3] * (11.0 * c[3] - 3.0 * c[1]);
164 case 4:
165 return s[4] * (11.0 * c[2] - 1.0);
166 case 5:
167 return s[5] * c[1];
168 case 6:
169 return s[6];
170 default:
171 return 0.0;
172 }
173 }
174 if (orderL == 8) {
175 switch (absM) {
176 case 0:
177 return (35.0 + c[2] * (-1260.0 + c[2] * (6930.0 + c[2] * (-12012.0 + c[2] * 6435.0)))) / 128.0;
178 case 1:
179 return s[1] * c[1] * (-315.0 + c[2] * (3465.0 + c[2] * (-9009.0 + c[2] * 6435.0))) / 16.0;
180 case 2:
181 return s[2] * (-315.0 + c[2] * (10395.0 + c[2] * (-45045.0 + c[2] * 45045.0))) / 16.0;
182 case 3:
183 return s[3] * c[1] * (10395.0 + c[2] * (-90090.0 + c[2] * 135135.0)) / 8.0;
184 case 4:
185 return s[4] * (10395.0 + c[2] * (-270270.0 + c[2] * 675675.0)) / 8.0;
186 case 5:
187 return s[5] * c[1] * (-135135.0 + c[2] * 675675.0) / 2.0;
188 case 6:
189 return s[6] * (-135135.0 + c[2] * 2027025.0) / 2.0;
190 case 7:
191 return s[7] * c[1] * 2027025.0;
192 case 8:
193 return s[8] * 2027025.0;
194 default:
195 return 0.0;
196 }
197 }
198 return 0.0;
199}

◆ minImage()

void seams::steinhardt::minImage ( double & dx,
double & dy,
double & dz,
double bx,
double by,
double bz )
inline

Definition at line 21 of file steinhardt_device.hpp.

22 {
23 if (dx < -0.5 * bx) {
24 dx += bx;
25 }
26 if (dx >= 0.5 * bx) {
27 dx -= bx;
28 }
29 if (dy < -0.5 * by) {
30 dy += by;
31 }
32 if (dy >= 0.5 * by) {
33 dy -= by;
34 }
35 if (dz < -0.5 * bz) {
36 dz += bz;
37 }
38 if (dz >= 0.5 * bz) {
39 dz -= bz;
40 }
41}

◆ normLM()

double seams::steinhardt::normLM ( int orderL,
int absM )
inline

Definition at line 43 of file steinhardt_device.hpp.

43 {
44 constexpr double pi = 3.14159265358979323846;
45 if (orderL == 3) {
46 switch (absM) {
47 case 0:
48 return 0.25 * std::sqrt(7.0 / pi);
49 case 1:
50 return 0.125 * std::sqrt(21.0 / pi);
51 case 2:
52 return 0.25 * std::sqrt(105.0 / (2.0 * pi));
53 case 3:
54 return 0.125 * std::sqrt(35.0 / pi);
55 default:
56 return 0.0;
57 }
58 }
59 if (orderL == 4) {
60 switch (absM) {
61 case 0:
62 return 0.1875 * std::sqrt(1.0 / pi);
63 case 1:
64 return 0.375 * std::sqrt(5.0 / pi);
65 case 2:
66 return 0.375 * std::sqrt(5.0 / (2.0 * pi));
67 case 3:
68 return 0.375 * std::sqrt(35.0 / pi);
69 case 4:
70 return 0.1875 * std::sqrt(35.0 / (2.0 * pi));
71 default:
72 return 0.0;
73 }
74 }
75 if (orderL == 6) {
76 switch (absM) {
77 case 0:
78 return 0.03125 * std::sqrt(13.0 / pi);
79 case 1:
80 return 0.0625 * std::sqrt(273.0 / (2.0 * pi));
81 case 2:
82 return 0.015625 * std::sqrt(1365.0 / pi);
83 case 3:
84 return 0.03125 * std::sqrt(1365.0 / pi);
85 case 4:
86 return 0.09375 * std::sqrt(91.0 / (2.0 * pi));
87 case 5:
88 return 0.09375 * std::sqrt(1001.0 / pi);
89 case 6:
90 return 0.015625 * std::sqrt(3003.0 / pi);
91 default:
92 return 0.0;
93 }
94 }
95 if (orderL == 8) {
96 switch (absM) {
97 case 0:
98 return std::sqrt(17.0 / (4.0 * pi));
99 case 1:
100 return std::sqrt(17.0 / (288.0 * pi));
101 case 2:
102 return std::sqrt(17.0 / (20160.0 * pi));
103 case 3:
104 return std::sqrt(17.0 / (1330560.0 * pi));
105 case 4:
106 return std::sqrt(17.0 / (79833600.0 * pi));
107 case 5:
108 return std::sqrt(17.0 / (4151347200.0 * pi));
109 case 6:
110 return std::sqrt(17.0 / (174356582400.0 * pi));
111 case 7:
112 return std::sqrt(17.0 / (5230697472000.0 * pi));
113 case 8:
114 return std::sqrt(17.0 / (83691159552000.0 * pi));
115 default:
116 return 0.0;
117 }
118 }
119 return 0.0;
120}

◆ qlmAddBond()

void seams::steinhardt::qlmAddBond ( int orderL,
double dx,
double dy,
double dz,
double * qlmInterleaved,
int row,
int nComp,
int & nUsed )
inline

Definition at line 245 of file steinhardt_device.hpp.

247 {
248 const double r2 = dx * dx + dy * dy + dz * dz;
249 if (r2 == 0.0) {
250 return;
251 }
252 const double r = std::sqrt(r2);
253 const double phi = std::atan2(dx, dy);
254 const double theta = std::acos(dz / r);
255 double ylm[34];
256 ylmAll(orderL, theta, phi, ylm);
257 for (int m = 0; m < nComp; m++) {
258 qlmInterleaved[2 * (row + m)] += ylm[2 * m];
259 qlmInterleaved[2 * (row + m) + 1] += ylm[2 * m + 1];
260 }
261 nUsed++;
262}
seams::steinhardt::ylmAll
void ylmAll(int orderL, double theta, double phi, double *out)
Definition steinhardt_device.hpp:202

◆ qlmOneAtom()

void seams::steinhardt::qlmOneAtom ( int iatom,
int orderL,
const double * xyz,
const int * offsets,
const int * cols,
double bx,
double by,
double bz,
double * qlmInterleaved )
inline

Definition at line 266 of file steinhardt_device.hpp.

268 {
269 (void)bx;
270 (void)by;
271 (void)bz;
272 const int nComp = 2 * orderL + 1;
273 const int iOff = 3 * iatom;
274 const int row = iatom * nComp;
275 for (int m = 0; m < nComp; m++) {
276 qlmInterleaved[2 * (row + m)] = 0.0;
277 qlmInterleaved[2 * (row + m) + 1] = 0.0;
278 }
279 const int j0 = offsets[iatom];
280 const int j1 = offsets[iatom + 1];
281 int nUsed = 0;
282 for (int p = j0; p < j1; p++) {
283 const int jatom = cols[p];
284 const double dx = xyz[iOff] - xyz[3 * jatom];
285 const double dy = xyz[iOff + 1] - xyz[3 * jatom + 1];
286 const double dz = xyz[iOff + 2] - xyz[3 * jatom + 2];
287 qlmAddBond(orderL, dx, dy, dz, qlmInterleaved, row, nComp, nUsed);
288 }
289 if (nUsed == 0) {
290 return;
291 }
292 const double inv = 1.0 / static_cast<double>(nUsed);
293 for (int m = 0; m < nComp; m++) {
294 qlmInterleaved[2 * (row + m)] *= inv;
295 qlmInterleaved[2 * (row + m) + 1] *= inv;
296 }
297}
seams::steinhardt::qlmAddBond
void qlmAddBond(int orderL, double dx, double dy, double dz, double *qlmInterleaved, int row, int nComp, int &nUsed)
Definition steinhardt_device.hpp:245

◆ qlmOneAtomDr()

void seams::steinhardt::qlmOneAtomDr ( int iatom,
int orderL,
const double * dr,
const int * offsets,
const int * cols,
double * qlmInterleaved )
inline

Definition at line 301 of file steinhardt_device.hpp.

303 {
304 const int nComp = 2 * orderL + 1;
305 const int row = iatom * nComp;
306 for (int m = 0; m < nComp; m++) {
307 qlmInterleaved[2 * (row + m)] = 0.0;
308 qlmInterleaved[2 * (row + m) + 1] = 0.0;
309 }
310 const int j0 = offsets[iatom];
311 const int j1 = offsets[iatom + 1];
312 int nUsed = 0;
313 for (int p = j0; p < j1; p++) {
314 qlmAddBond(orderL, dr[3 * p], dr[3 * p + 1], dr[3 * p + 2],
315 qlmInterleaved, row, nComp, nUsed);
316 }
317 if (nUsed == 0) {
318 return;
319 }
320 const double inv = 1.0 / static_cast<double>(nUsed);
321 for (int m = 0; m < nComp; m++) {
322 qlmInterleaved[2 * (row + m)] *= inv;
323 qlmInterleaved[2 * (row + m) + 1] *= inv;
324 }
325}

◆ qlOneAtom()

void seams::steinhardt::qlOneAtom ( int iatom,
int orderL,
const double * qlmInterleaved,
const int * offsets,
const int * cols,
double * ql,
double * qlBar )
inline

Definition at line 327 of file steinhardt_device.hpp.

329 {
330 constexpr double pi = 3.14159265358979323846;
331 const int nComp = 2 * orderL + 1;
332 const double prefactor = 4.0 * pi / (2.0 * orderL + 1.0);
333 const int row = iatom * nComp;
334 double sumLocal = 0.0;
335 for (int m = 0; m < nComp; m++) {
336 const double re = qlmInterleaved[2 * (row + m)];
337 const double im = qlmInterleaved[2 * (row + m) + 1];
338 sumLocal += re * re + im * im;
339 }
340 ql[iatom] = std::sqrt(prefactor * sumLocal);
341
342 const int j0 = offsets[iatom];
343 const int j1 = offsets[iatom + 1];
344 if (j0 == j1) {
345 qlBar[iatom] = ql[iatom];
346 return;
347 }
348 double barRe[17];
349 double barIm[17];
350 for (int m = 0; m < nComp; m++) {
351 barRe[m] = qlmInterleaved[2 * (row + m)];
352 barIm[m] = qlmInterleaved[2 * (row + m) + 1];
353 }
354 int nContrib = 1;
355 for (int p = j0; p < j1; p++) {
356 const int jatom = cols[p];
357 const int jRow = jatom * nComp;
358 for (int m = 0; m < nComp; m++) {
359 barRe[m] += qlmInterleaved[2 * (jRow + m)];
360 barIm[m] += qlmInterleaved[2 * (jRow + m) + 1];
361 }
362 nContrib++;
363 }
364 const double inv = 1.0 / static_cast<double>(nContrib);
365 double sumBar = 0.0;
366 for (int m = 0; m < nComp; m++) {
367 const double re = barRe[m] * inv;
368 const double im = barIm[m] * inv;
369 sumBar += re * re + im * im;
370 }
371 qlBar[iatom] = std::sqrt(prefactor * sumBar);
372}

◆ ylmAll()

void seams::steinhardt::ylmAll ( int orderL,
double theta,
double phi,
double * out )
inline

Definition at line 202 of file steinhardt_device.hpp.

202 {
203 const int nComp = 2 * orderL + 1;
204 for (int k = 0; k < nComp; k++) {
205 out[2 * k] = 0.0;
206 out[2 * k + 1] = 0.0;
207 }
208
209 double s[9];
210 double c[9];
211 double pr[9];
212 double pi[9];
213 const double sinT = std::sin(theta);
214 const double cosT = std::cos(theta);
215 const double cphi = std::cos(phi);
216 const double sphi = std::sin(phi);
217 s[0] = 1.0;
218 c[0] = 1.0;
219 pr[0] = 1.0;
220 pi[0] = 0.0;
221 for (int k = 1; k <= orderL; k++) {
222 s[k] = s[k - 1] * sinT;
223 c[k] = c[k - 1] * cosT;
224 pr[k] = pr[k - 1] * cphi - pi[k - 1] * sphi;
225 pi[k] = pr[k - 1] * sphi + pi[k - 1] * cphi;
226 }
227
228 for (int absM = 0; absM <= orderL; absM++) {
229 const double amp = normLM(orderL, absM) * legendreAmp(orderL, absM, s, c);
230 // Y_{l,-m} = amp * conj(e^{i m phi})
231 const double nre = amp * pr[absM];
232 const double nim = -amp * pi[absM];
233 const int ineg = orderL - absM;
234 out[2 * ineg] = nre;
235 out[2 * ineg + 1] = nim;
236 // Y_{l,+m} = (-1)^m amp * e^{i m phi}
237 const double sign = (absM % 2 == 0) ? 1.0 : -1.0;
238 const int ipos = orderL + absM;
239 out[2 * ipos] = sign * amp * pr[absM];
240 out[2 * ipos + 1] = sign * amp * pi[absM];
241 }
242}
seams::steinhardt::legendreAmp
double legendreAmp(int orderL, int absM, const double *s, const double *c)
Definition steinhardt_device.hpp:122
seams::steinhardt::normLM
double normLM(int orderL, int absM)
Definition steinhardt_device.hpp:43