LIBINT 2.9.0
cgshellinfo.h
1/*
2 * Copyright (C) 2004-2024 Edward F. Valeev
3 *
4 * This file is part of Libint library.
5 *
6 * Libint library is free software: you can redistribute it and/or modify
7 * it under the terms of the GNU Lesser General Public License as published by
8 * the Free Software Foundation, either version 3 of the License, or
9 * (at your option) any later version.
10 *
11 * Libint library is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 * GNU Lesser General Public License for more details.
15 *
16 * You should have received a copy of the GNU Lesser General Public License
17 * along with Libint library. If not, see <http://www.gnu.org/licenses/>.
18 *
19 */
20
21#ifndef _libint2_src_bin_libint_cgshellinfo_h_
22#define _libint2_src_bin_libint_cgshellinfo_h_
23
24#include <algorithm>
25#include <cassert>
26#include <utility>
27
28namespace libint2 {
29
31enum class CartesianShellNormalization {
32 standard, // same normalization factor for every function in shell (default)
33 uniform // different normalization factors so that each function has same
34 // norm
35};
36
37namespace detail {
38inline int notxyz(int a, int b) {
39 assert(a != b);
40 int amax = std::max(a, b);
41 int amin = std::min(a, b);
42 if (amin == 0 && amax == 1) return 2;
43 if (amin == 0 && amax == 2) return 1;
44 if (amin == 1 && amax == 2) return 0;
45 assert(false); // unreachable
46 return -1;
47}
48
49inline std::pair<int, int> notxyz(int a) {
50 switch (a) {
51 case 0:
52 return std::make_pair(1, 2);
53 break;
54 case 1:
55 return std::make_pair(0, 2);
56 break;
57 case 2:
58 return std::make_pair(0, 1);
59 break;
60 }
61 assert(false); // unreachable
62 return std::make_pair(-1, -1);
63}
64} // namespace detail
65
66template <CGShellOrdering Ord, unsigned int lmax>
67struct CGShellOrderingGenerator {
68 static void compute(int (&cartindex)[lmax + 1][lmax + 1][lmax + 1]);
69};
70template <unsigned int lmax>
71struct CGShellOrderingGenerator<CGShellOrdering_Standard, lmax> {
72 static void compute(int (&cartindex)[lmax + 1][lmax + 1][lmax + 1]) {
73 // see cgshell_ordering.h
74 for (unsigned int am = 0; am <= lmax; ++am) {
75 int count = 0;
76 for (int i = (am); (i) >= 0; (i)--) {
77 for (int j = (am) - (i); (j) >= 0; (j)--, ++count) {
78 cartindex[am][i][j] = count;
79 }
80 }
81 }
82 }
83};
84template <unsigned int lmax>
85struct CGShellOrderingGenerator<CGShellOrdering_GAMESS, lmax> {
86 static void compute(int (&cartindex)[lmax + 1][lmax + 1][lmax + 1]) {
87 for (unsigned int am = 0; am <= lmax; ++am) {
88 if (am == 0) {
89 cartindex[0][0][0] = 0;
90 continue;
91 }
94 if (am == 4) {
95 cartindex[4][4][0] = 0;
96 cartindex[4][0][4] = 1;
97 cartindex[4][0][0] = 2;
98 cartindex[4][3][1] = 3;
99 cartindex[4][3][0] = 4;
100 cartindex[4][1][3] = 5;
101 cartindex[4][0][3] = 6;
102 cartindex[4][1][0] = 7;
103 cartindex[4][0][1] = 8;
104 cartindex[4][2][2] = 9;
105 cartindex[4][2][0] = 10;
106 cartindex[4][0][2] = 11;
107 cartindex[4][2][1] = 12;
108 cartindex[4][1][2] = 13;
109 cartindex[4][1][1] = 14;
110 continue;
111 }
112
113 unsigned int current_index = 0;
114 unsigned int qn[3] = {0, 0, 0};
115
116 const int ammin = ((int)am + 2) / 3;
117 for (int am1 = am; am1 >= ammin; --am1) {
118 for (int xyz1 = 0; xyz1 < 3; ++xyz1) {
119 qn[xyz1] = am1;
120
121 // distribute the remaining quanta according to the following rules
122
123 // "nothing to distribute" is a special case
124 if (am - am1 == 0) {
125 std::pair<int, int> xyz(detail::notxyz(xyz1));
126 qn[xyz.first] = 0;
127 qn[xyz.second] = 0;
128 cartindex[am][qn[0]][qn[1]] = current_index;
129 ++current_index;
130 } else {
131 int am23 = (int)am - qn[xyz1];
132 const int maxam23 = std::min((int)qn[xyz1], am23);
133 const int minam23 = (am23 + 1) / 2;
134 for (int am2 = maxam23; am2 >= minam23; --am2) {
135 const int xyz2min = (am2 == qn[xyz1]) ? xyz1 + 1 : 0;
136 for (int xyz2 = xyz2min; xyz2 < 3; ++xyz2) {
137 if (xyz1 == xyz2) continue;
138 qn[xyz2] = am2;
139 const int xyz3 = detail::notxyz(xyz1, xyz2);
140 qn[xyz3] = am23 - am2;
141 if ((qn[xyz3] == qn[xyz1] && xyz3 < xyz1) ||
142 (qn[xyz3] == qn[xyz2] && xyz3 < xyz2))
143 continue;
144 {
145 cartindex[am][qn[0]][qn[1]] = current_index;
146 ++current_index;
147 }
148 }
149 }
150 }
151 }
152 }
153 }
154 }
155};
156template <unsigned int lmax>
157struct CGShellOrderingGenerator<CGShellOrdering_ORCA, lmax> {
158 static void compute(int (&cartindex)[lmax + 1][lmax + 1][lmax + 1]) {
159 // identical to GAMESS for s through g functions
160 // identical to standard for h and higher
161 CGShellOrderingGenerator<CGShellOrdering_Standard, lmax>::compute(
162 cartindex);
163
164 for (unsigned int am = 0; am <= std::min(4u, lmax); ++am) {
165 if (am == 0) {
166 cartindex[0][0][0] = 0;
167 continue;
168 }
169 if (am == 1) {
170 cartindex[1][0][0] = 0;
171 cartindex[1][1][0] = 1;
172 cartindex[1][0][1] = 2;
173 continue;
174 }
175 if (am == 2) {
176 cartindex[2][2][0] = 0;
177 cartindex[2][0][2] = 1;
178 cartindex[2][0][0] = 2;
179 cartindex[2][1][1] = 3;
180 cartindex[2][1][0] = 4;
181 cartindex[2][0][1] = 5;
182 continue;
183 }
184 if (am == 3) {
185 cartindex[3][3][0] = 0;
186 cartindex[3][0][3] = 1;
187 cartindex[3][0][0] = 2;
188 cartindex[3][2][1] = 3;
189 cartindex[3][2][0] = 4;
190 cartindex[3][1][2] = 5;
191 cartindex[3][0][2] = 6;
192 cartindex[3][1][0] = 7;
193 cartindex[3][0][1] = 8;
194 cartindex[3][1][1] = 9;
195 continue;
196 }
197 if (am == 4) {
198 cartindex[4][4][0] = 0;
199 cartindex[4][0][4] = 1;
200 cartindex[4][0][0] = 2;
201 cartindex[4][3][1] = 3;
202 cartindex[4][3][0] = 4;
203 cartindex[4][1][3] = 5;
204 cartindex[4][0][3] = 6;
205 cartindex[4][1][0] = 7;
206 cartindex[4][0][1] = 8;
207 cartindex[4][2][2] = 9;
208 cartindex[4][2][0] = 10;
209 cartindex[4][0][2] = 11;
210 cartindex[4][2][1] = 12;
211 cartindex[4][1][2] = 13;
212 cartindex[4][1][1] = 14;
213 continue;
214 }
215 /*
216 if (am == 5) {
217 cartindex[5][5][0] = 0;
218 cartindex[5][4][1] = 1;
219 cartindex[5][4][0] = 2;
220 cartindex[5][3][2] = 3;
221 cartindex[5][3][1] = 4;
222 cartindex[5][3][0] = 5;
223 cartindex[5][2][3] = 6;
224 cartindex[5][2][2] = 7;
225 cartindex[5][2][1] = 8;
226 cartindex[5][2][0] = 9;
227 cartindex[5][1][4] = 10;
228 cartindex[5][1][3] = 11;
229 cartindex[5][1][2] = 12;
230 cartindex[5][1][1] = 13;
231 cartindex[5][1][0] = 14;
232 cartindex[5][0][5] = 15;
233 cartindex[5][0][4] = 16;
234 cartindex[5][0][3] = 17;
235 cartindex[5][0][2] = 18;
236 cartindex[5][0][1] = 19;
237 cartindex[5][0][0] = 20;
238 continue;
239 }
240 if (am == 6) {
241 cartindex[6][6][0] = 0;
242 cartindex[6][5][1] = 1;
243 cartindex[6][5][0] = 2;
244 cartindex[6][4][2] = 3;
245 cartindex[6][4][1] = 4;
246 cartindex[6][4][0] = 5;
247 cartindex[6][3][3] = 6;
248 cartindex[6][3][2] = 7;
249 cartindex[6][3][1] = 8;
250 cartindex[6][3][0] = 9;
251 cartindex[6][2][4] = 10;
252 cartindex[6][2][3] = 11;
253 cartindex[6][2][2] = 12;
254 cartindex[6][2][1] = 13;
255 cartindex[6][2][0] = 14;
256 cartindex[6][1][5] = 15;
257 cartindex[6][1][4] = 16;
258 cartindex[6][1][3] = 17;
259 cartindex[6][1][2] = 18;
260 cartindex[6][1][1] = 19;
261 cartindex[6][1][0] = 20;
262 cartindex[6][0][6] = 21;
263 cartindex[6][0][5] = 22;
264 cartindex[6][0][4] = 23;
265 cartindex[6][0][3] = 24;
266 cartindex[6][0][2] = 25;
267 cartindex[6][0][1] = 26;
268 cartindex[6][0][0] = 27;
269 continue;
270 }
271 if (am == 7) {
272 cartindex[7][7][0] = 0;
273 cartindex[7][6][1] = 1;
274 cartindex[7][6][0] = 2;
275 cartindex[7][5][2] = 3;
276 cartindex[7][5][1] = 4;
277 cartindex[7][5][0] = 5;
278 cartindex[7][4][3] = 6;
279 cartindex[7][4][2] = 7;
280 cartindex[7][4][1] = 8;
281 cartindex[7][4][0] = 9;
282 cartindex[7][3][4] = 10;
283 cartindex[7][3][3] = 11;
284 cartindex[7][3][2] = 12;
285 cartindex[7][3][1] = 13;
286 cartindex[7][3][0] = 14;
287 cartindex[7][2][5] = 15;
288 cartindex[7][2][4] = 16;
289 cartindex[7][2][3] = 17;
290 cartindex[7][2][2] = 18;
291 cartindex[7][2][1] = 19;
292 cartindex[7][2][0] = 20;
293 cartindex[7][1][6] = 21;
294 cartindex[7][1][5] = 22;
295 cartindex[7][1][4] = 23;
296 cartindex[7][1][3] = 24;
297 cartindex[7][1][2] = 25;
298 cartindex[7][1][1] = 26;
299 cartindex[7][1][0] = 27;
300 cartindex[7][0][7] = 28;
301 cartindex[7][0][6] = 29;
302 cartindex[7][0][5] = 30;
303 cartindex[7][0][4] = 31;
304 cartindex[7][0][3] = 32;
305 cartindex[7][0][2] = 33;
306 cartindex[7][0][1] = 34;
307 cartindex[7][0][0] = 35;
308 continue;
309 }
310 if (am == 8) {
311 cartindex[8][8][0] = 0;
312 cartindex[8][7][1] = 1;
313 cartindex[8][7][0] = 2;
314 cartindex[8][6][2] = 3;
315 cartindex[8][6][1] = 4;
316 cartindex[8][6][0] = 5;
317 cartindex[8][5][3] = 6;
318 cartindex[8][5][2] = 7;
319 cartindex[8][5][1] = 8;
320 cartindex[8][5][0] = 9;
321 cartindex[8][4][4] = 10;
322 cartindex[8][4][3] = 11;
323 cartindex[8][4][2] = 12;
324 cartindex[8][4][1] = 13;
325 cartindex[8][4][0] = 14;
326 cartindex[8][3][5] = 15;
327 cartindex[8][3][4] = 16;
328 cartindex[8][3][3] = 17;
329 cartindex[8][3][2] = 18;
330 cartindex[8][3][1] = 19;
331 cartindex[8][3][0] = 20;
332 cartindex[8][2][6] = 21;
333 cartindex[8][2][5] = 22;
334 cartindex[8][2][4] = 23;
335 cartindex[8][2][3] = 24;
336 cartindex[8][2][2] = 25;
337 cartindex[8][2][1] = 26;
338 cartindex[8][2][0] = 27;
339 cartindex[8][1][7] = 28;
340 cartindex[8][1][6] = 29;
341 cartindex[8][1][5] = 30;
342 cartindex[8][1][4] = 31;
343 cartindex[8][1][3] = 32;
344 cartindex[8][1][2] = 33;
345 cartindex[8][1][1] = 34;
346 cartindex[8][1][0] = 35;
347 cartindex[8][0][8] = 36;
348 cartindex[8][0][7] = 37;
349 cartindex[8][0][6] = 38;
350 cartindex[8][0][5] = 39;
351 cartindex[8][0][4] = 40;
352 cartindex[8][0][3] = 41;
353 cartindex[8][0][2] = 42;
354 cartindex[8][0][1] = 43;
355 cartindex[8][0][0] = 44;
356 continue;
357 }
358 */
359 }
360
361#if 0
362 for(int l=0; l<=lmax; ++l) {
363 for(int i=0; i<=l; ++i) {
364 for(int j=0; j<=l-i; ++j) {
365 std::cout << "CGShellOrderingGenerator<CGShellOrdering_ORCA>: cartindex[" << l << "]["
366 << i << "][" << j << "] = " << cartindex[l][i][j] << std::endl;
367 }
368 }
369 }
370#endif
371 }
372};
373
374// see http://www.cmbi.ru.nl/molden/molden_format.html
375template <unsigned int lmax>
376struct CGShellOrderingGenerator<CGShellOrdering_MOLDEN, lmax> {
377 static void compute(int (&cartindex)[lmax + 1][lmax + 1][lmax + 1]) {
378 for (unsigned int am = 0; am <= 4u; ++am) {
379 if (am == 0) {
380 cartindex[0][0][0] = 0;
381 continue;
382 }
383 if (am == 1) {
384 cartindex[1][1][0] = 0;
385 cartindex[1][0][1] = 1;
386 cartindex[1][0][0] = 2;
387 continue;
388 }
389 if (am == 2) {
390 cartindex[2][2][0] = 0;
391 cartindex[2][0][2] = 1;
392 cartindex[2][0][0] = 2;
393 cartindex[2][1][1] = 3;
394 cartindex[2][1][0] = 4;
395 cartindex[2][0][1] = 5;
396 continue;
397 }
398 if (am == 3) {
399 cartindex[3][3][0] = 0;
400 cartindex[3][0][3] = 1;
401 cartindex[3][0][0] = 2;
402 cartindex[3][1][2] = 3;
403 cartindex[3][2][1] = 4;
404 cartindex[3][2][0] = 5;
405 cartindex[3][1][0] = 6;
406 cartindex[3][0][1] = 7;
407 cartindex[3][0][2] = 8;
408 cartindex[3][1][1] = 9;
409 continue;
410 }
411 if (am == 4) {
412 cartindex[4][4][0] = 0;
413 cartindex[4][0][4] = 1;
414 cartindex[4][0][0] = 2;
415 cartindex[4][3][1] = 3;
416 cartindex[4][3][0] = 4;
417 cartindex[4][1][3] = 5;
418 cartindex[4][0][3] = 6;
419 cartindex[4][1][0] = 7;
420 cartindex[4][0][1] = 8;
421 cartindex[4][2][2] = 9;
422 cartindex[4][2][0] = 10;
423 cartindex[4][0][2] = 11;
424 cartindex[4][2][1] = 12;
425 cartindex[4][1][2] = 13;
426 cartindex[4][1][1] = 14;
427 continue;
428 }
429 }
430 }
431};
432
433template <CGShellOrdering Ord, unsigned int lmax = LIBINT_CARTGAUSS_MAX_AM>
434struct CGShellOrderingData {
435 struct Triple {
436 Triple() : i(0), j(0), k(0) {}
437 Triple(int ii, int jj, int kk) : i(ii), j(jj), k(kk) {}
438 int i, j, k;
439 };
440
441 CGShellOrderingData() {
442 // compute cartindex
443 CGShellOrderingGenerator<Ord, lmax>::compute(cartindex);
444 // then use it to compute cartindex_to_ijk
445 for (unsigned int l = 0; l <= lmax; ++l) {
446 for (unsigned int i = 0; i <= l; ++i) {
447 for (unsigned int j = 0; j <= l - i; ++j) {
448 const int c = cartindex[l][i][j];
449 cartindex_to_ijk[l][c] = Triple(i, j, l - i - j);
450 }
451 }
452 }
453 }
454
455 int cartindex[lmax + 1][lmax + 1][lmax + 1];
456 Triple cartindex_to_ijk[lmax + 1][(lmax + 1) * (lmax + 2) / 2];
457};
458
461template <typename OrderingData>
462struct CGShellInfo {
463 // computes cartindex xyz from am i j
464 static int cartindex(unsigned int am, int i, int j) {
465 return data_.cartindex[am][i][j];
466 }
467 // computes i j k from cartindex xyz
468 template <typename I1, typename I2, typename I3>
469 static void cartindex_to_ijk(I1 am, I2 xyz, I3& i, I3& j, I3& k) {
470 const typename OrderingData::Triple& ijk = data_.cartindex_to_ijk[am][xyz];
471 i = ijk.i;
472 j = ijk.j;
473 k = ijk.k;
474 }
475
476 private:
477 static OrderingData data_;
478};
479
480template <typename OrderingData>
481OrderingData CGShellInfo<OrderingData>::data_;
482
483}; // namespace libint2
484
485#endif // header guard
libint2
Defaults definitions for various parameters assumed by Libint.
Definition algebra.cc:24
libint2::CartesianShellNormalization
CartesianShellNormalization
Normalization convention for Cartesian Gaussian shells.
Definition cgshellinfo.h:31
libint2::CGShellInfo
provides ordering maps for up to angular momentum lmax and ordering specified by CGShellOrderingSpec
Definition cgshellinfo.h:462
libint2::CGShellOrderingData::Triple
Definition cgshellinfo.h:435
libint2::CGShellOrderingData
Definition cgshellinfo.h:434
libint2::CGShellOrderingGenerator< CGShellOrdering_GAMESS, lmax >::compute
static void compute(int(&cartindex)[lmax+1][lmax+1][lmax+1])
Definition cgshellinfo.h:86
libint2::CGShellOrderingGenerator
Definition cgshellinfo.h:67
Generated by doxygen 1.10.0