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 ylmAllTrig (int orderL, double sinT, double cosT, double cphi, double sphi, double *out)
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 252 of file steinhardt_device.hpp.

254 {
255 const double r2 = dx * dx + dy * dy + dz * dz;
256 if (r2 == 0.0) {
257 return;
258 }
259 const double r = std::sqrt(r2);
260 const double invr = 1.0 / r;
261 double cosT = dz * invr;
262 if (cosT > 1.0) {
263 cosT = 1.0;
264 } else if (cosT < -1.0) {
265 cosT = -1.0;
266 }
267 const double rho2 = dx * dx + dy * dy;
268 double sinT = 0.0;
269 double cphi = 1.0;
270 double sphi = 0.0;
271 if (rho2 != 0.0) {
272 const double rho = std::sqrt(rho2);
273 sinT = rho * invr;
274 const double invrho = 1.0 / rho;
275 cphi = dy * invrho;
276 sphi = dx * invrho;
277 }
278 double ylm[34];
279 ylmAllTrig(orderL, sinT, cosT, cphi, sphi, ylm);
280 for (int m = 0; m < nComp; m++) {
281 qlmInterleaved[2 * (row + m)] += ylm[2 * m];
282 qlmInterleaved[2 * (row + m) + 1] += ylm[2 * m + 1];
283 }
284 nUsed++;
285}
seams::steinhardt::ylmAllTrig
void ylmAllTrig(int orderL, double sinT, double cosT, double cphi, double sphi, double *out)
Definition steinhardt_device.hpp:205

◆ 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 289 of file steinhardt_device.hpp.

291 {
292 (void)bx;
293 (void)by;
294 (void)bz;
295 const int nComp = 2 * orderL + 1;
296 const int iOff = 3 * iatom;
297 const int row = iatom * nComp;
298 for (int m = 0; m < nComp; m++) {
299 qlmInterleaved[2 * (row + m)] = 0.0;
300 qlmInterleaved[2 * (row + m) + 1] = 0.0;
301 }
302 const int j0 = offsets[iatom];
303 const int j1 = offsets[iatom + 1];
304 int nUsed = 0;
305 for (int p = j0; p < j1; p++) {
306 const int jatom = cols[p];
307 const double dx = xyz[iOff] - xyz[3 * jatom];
308 const double dy = xyz[iOff + 1] - xyz[3 * jatom + 1];
309 const double dz = xyz[iOff + 2] - xyz[3 * jatom + 2];
310 qlmAddBond(orderL, dx, dy, dz, qlmInterleaved, row, nComp, nUsed);
311 }
312 if (nUsed == 0) {
313 return;
314 }
315 const double inv = 1.0 / static_cast<double>(nUsed);
316 for (int m = 0; m < nComp; m++) {
317 qlmInterleaved[2 * (row + m)] *= inv;
318 qlmInterleaved[2 * (row + m) + 1] *= inv;
319 }
320}
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:252

◆ qlmOneAtomDr()

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

Definition at line 324 of file steinhardt_device.hpp.

326 {
327 const int nComp = 2 * orderL + 1;
328 const int row = iatom * nComp;
329 for (int m = 0; m < nComp; m++) {
330 qlmInterleaved[2 * (row + m)] = 0.0;
331 qlmInterleaved[2 * (row + m) + 1] = 0.0;
332 }
333 const int j0 = offsets[iatom];
334 const int j1 = offsets[iatom + 1];
335 int nUsed = 0;
336 for (int p = j0; p < j1; p++) {
337 qlmAddBond(orderL, dr[3 * p], dr[3 * p + 1], dr[3 * p + 2],
338 qlmInterleaved, row, nComp, nUsed);
339 }
340 if (nUsed == 0) {
341 return;
342 }
343 const double inv = 1.0 / static_cast<double>(nUsed);
344 for (int m = 0; m < nComp; m++) {
345 qlmInterleaved[2 * (row + m)] *= inv;
346 qlmInterleaved[2 * (row + m) + 1] *= inv;
347 }
348}

◆ 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 350 of file steinhardt_device.hpp.

352 {
353 constexpr double pi = 3.14159265358979323846;
354 const int nComp = 2 * orderL + 1;
355 const double prefactor = 4.0 * pi / (2.0 * orderL + 1.0);
356 const int row = iatom * nComp;
357 double sumLocal = 0.0;
358 for (int m = 0; m < nComp; m++) {
359 const double re = qlmInterleaved[2 * (row + m)];
360 const double im = qlmInterleaved[2 * (row + m) + 1];
361 sumLocal += re * re + im * im;
362 }
363 ql[iatom] = std::sqrt(prefactor * sumLocal);
364
365 const int j0 = offsets[iatom];
366 const int j1 = offsets[iatom + 1];
367 if (j0 == j1) {
368 qlBar[iatom] = ql[iatom];
369 return;
370 }
371 double barRe[17];
372 double barIm[17];
373 for (int m = 0; m < nComp; m++) {
374 barRe[m] = qlmInterleaved[2 * (row + m)];
375 barIm[m] = qlmInterleaved[2 * (row + m) + 1];
376 }
377 int nContrib = 1;
378 for (int p = j0; p < j1; p++) {
379 const int jatom = cols[p];
380 const int jRow = jatom * nComp;
381 for (int m = 0; m < nComp; m++) {
382 barRe[m] += qlmInterleaved[2 * (jRow + m)];
383 barIm[m] += qlmInterleaved[2 * (jRow + m) + 1];
384 }
385 nContrib++;
386 }
387 const double inv = 1.0 / static_cast<double>(nContrib);
388 double sumBar = 0.0;
389 for (int m = 0; m < nComp; m++) {
390 const double re = barRe[m] * inv;
391 const double im = barIm[m] * inv;
392 sumBar += re * re + im * im;
393 }
394 qlBar[iatom] = std::sqrt(prefactor * sumBar);
395}

◆ ylmAll()

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

Definition at line 244 of file steinhardt_device.hpp.

244 {
245 ylmAllTrig(orderL, std::sin(theta), std::cos(theta), std::cos(phi),
246 std::sin(phi), out);
247}

◆ ylmAllTrig()

void seams::steinhardt::ylmAllTrig ( int orderL,
double sinT,
double cosT,
double cphi,
double sphi,
double * out )
inline

Definition at line 205 of file steinhardt_device.hpp.

206 {
207 const int nComp = 2 * orderL + 1;
208 for (int k = 0; k < nComp; k++) {
209 out[2 * k] = 0.0;
210 out[2 * k + 1] = 0.0;
211 }
212
213 double s[9];
214 double c[9];
215 double pr[9];
216 double pi[9];
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