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