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