LIBINT  2.6.0
cgshellinfo.h
1 /*
2  * Copyright (C) 2004-2019 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 
28 namespace libint2 {
29 
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
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 
446  // 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
Defaults definitions for various parameters assumed by Libint.
Definition: algebra.cc:24
Definition: cgshellinfo.h:439
CartesianShellNormalization
Normalization convention for Cartesian Gaussian shells.
Definition: cgshellinfo.h:31
Definition: cgshellinfo.h:437
Definition: cgshellinfo.h:63
static void compute(int(&cartindex)[lmax+1][lmax+1][lmax+1])
Definition: cgshellinfo.h:82
provides ordering maps for up to angular momentum lmax and ordering specified by CGShellOrderingSpec
Definition: cgshellinfo.h:464