LIBINT  2.6.0
shell.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_lib_libint_shell_h_
22 #define _libint2_src_lib_libint_shell_h_
23 
24 #include <libint2/util/cxxstd.h>
25 #if LIBINT2_CPLUSPLUS_STD < 2011
26 # error "libint2/shell.h requires C++11 support"
27 #endif
28 
29 #include <iostream>
30 #include <array>
31 #include <vector>
32 #include <cassert>
33 #include <cmath>
34 
35 #include <libint2.h>
36 
37 #include <libint2/util/small_vector.h>
38 
39 namespace libint2 {
40 
41  namespace math {
43  static constexpr std::array<int64_t,21> fac = {{1LL, 1LL, 2LL, 6LL, 24LL, 120LL, 720LL, 5040LL, 40320LL, 362880LL, 3628800LL, 39916800LL,
44  479001600LL, 6227020800LL, 87178291200LL, 1307674368000LL, 20922789888000LL,
45  355687428096000LL, 6402373705728000LL, 121645100408832000LL,
46  2432902008176640000LL}};
48  static constexpr std::array<int64_t,31> df_Kminus1 = {{1LL, 1LL, 1LL, 2LL, 3LL, 8LL, 15LL, 48LL, 105LL, 384LL, 945LL, 3840LL, 10395LL, 46080LL, 135135LL,
49  645120LL, 2027025LL, 10321920LL, 34459425LL, 185794560LL, 654729075LL,
50  3715891200LL, 13749310575LL, 81749606400LL, 316234143225LL, 1961990553600LL,
51  7905853580625LL, 51011754393600LL, 213458046676875LL, 1428329123020800LL,
52  6190283353629375LL}};
54  template <typename Int> int64_t bc(Int i, Int j) {
55  assert(i < Int(fac.size()));
56  assert(j < Int(fac.size()));
57  assert(i >= j);
58  return fac[i] / (fac[j] * fac[i-j]);
59  }
60  }
61 
63 
81  struct Shell {
82  typedef double real_t;
83 
85  struct Contraction {
86  int l;
87  bool pure;
88  svector<real_t> coeff;
89 
90  bool operator==(const Contraction& other) const {
91  return &other == this || (l == other.l && pure == other.pure && coeff == other.coeff);
92  }
93  bool operator!=(const Contraction& other) const {
94  return not this->operator==(other);
95  }
96  size_t cartesian_size() const {
97  return (l + 1) * (l + 2) / 2;
98  }
99  size_t size() const {
100  return pure ? (2 * l + 1) : cartesian_size();
101  }
102  };
103 
104  svector<real_t> alpha;
105  svector<Contraction> contr;
107  svector<real_t> max_ln_coeff;
108 
109  Shell() = default;
110  Shell(const Shell&) = default;
111  // intel does not support "move ctor = default"
112  Shell(Shell&& other) noexcept :
113  alpha(std::move(other.alpha)),
114  contr(std::move(other.contr)),
115  O(std::move(other.O)),
116  max_ln_coeff(std::move(other.max_ln_coeff)) {
117  }
118  Shell& operator=(const Shell&) = default;
119  // intel does not support "move asgnmt = default"
120  Shell& operator=(Shell&& other) noexcept {
121  alpha = std::move(other.alpha);
122  contr = std::move(other.contr);
123  O = std::move(other.O);
124  max_ln_coeff = std::move(other.max_ln_coeff);
125  return *this;
126  }
127  Shell(svector<real_t> _alpha,
128  svector<Contraction> _contr,
130  alpha(std::move(_alpha)),
131  contr(std::move(_contr)),
132  O(std::move(_O)) {
133  // embed normalization factors into contraction coefficients
134  renorm();
135  }
136 
137  Shell& move(std::array<real_t, 3> new_origin) {
138  O = std::move(new_origin);
139  return *this;
140  }
141 
142  size_t cartesian_size() const {
143  size_t s = 0;
144  for(const auto& c: contr) { s += c.cartesian_size(); }
145  return s;
146  }
147  size_t size() const {
148  size_t s = 0;
149  for(const auto& c: contr) { s += c.size(); }
150  return s;
151  }
152 
153  size_t ncontr() const { return contr.size(); }
154  size_t nprim() const { return alpha.size(); }
155 
156  bool operator==(const Shell& other) const {
157  return &other == this || (O == other.O && alpha == other.alpha && contr == other.contr);
158  }
159  bool operator!=(const Shell& other) const {
160  return not this->operator==(other);
161  }
162 
163  static char am_symbol(size_t l) {
164  static char lsymb[] = "spdfghikmnoqrtuvwxyz";
165  assert(l<=19);
166  return lsymb[l];
167  }
168  static unsigned short am_symbol_to_l(char am_symbol) {
169  const char AM_SYMBOL = ::toupper(am_symbol);
170  switch (AM_SYMBOL) {
171  case 'S': return 0;
172  case 'P': return 1;
173  case 'D': return 2;
174  case 'F': return 3;
175  case 'G': return 4;
176  case 'H': return 5;
177  case 'I': return 6;
178  case 'K': return 7;
179  case 'M': return 8;
180  case 'N': return 9;
181  case 'O': return 10;
182  case 'Q': return 11;
183  case 'R': return 12;
184  case 'T': return 13;
185  case 'U': return 14;
186  case 'V': return 15;
187  case 'W': return 16;
188  case 'X': return 17;
189  case 'Y': return 18;
190  case 'Z': return 19;
191  default: throw "invalid angular momentum label";
192  }
193  }
194 
196  typedef enum {false_value=0,true_value=1,default_value=2} value_t;
197  defaultable_boolean() : value_(default_value) {}
198  defaultable_boolean(bool v) : value_(static_cast<value_t>(v?1:0)) {}
199  bool is_default() const { return value_ == default_value; }
200  operator bool() const { assert(value_ != default_value); return value_ == true_value; }
201  private:
202  value_t value_;
203  };
204 
210  static bool result{true};
211  if (not flag.is_default()) {
212  result = flag;
213  }
214  return result;
215  }
216 
218  static const Shell& unit() {
219  static const Shell unitshell{make_unit()};
220  return unitshell;
221  }
222 
225  real_t coeff_normalized(size_t c, size_t p) const {
226  const auto alpha = this->alpha.at(p);
227  assert(alpha >= 0.0);
228  const auto l = contr.at(c).l;
229  assert(l <= 15); // due to df_Kminus1[] a 64-bit integer type; kinda ridiculous restriction anyway
230 
231  using libint2::math::df_Kminus1;
232  const auto sqrt_Pi_cubed = real_t{5.56832799683170784528481798212};
233  const auto two_alpha = 2 * alpha;
234  const auto two_alpha_to_am32 = pow(two_alpha,l+1) * sqrt(two_alpha);
235  const auto one_over_N = sqrt((sqrt_Pi_cubed * df_Kminus1[2*l] )/(pow(2,l) * two_alpha_to_am32));
236  return contr.at(c).coeff[p] * one_over_N;
237  }
238 
239  private:
240 
241  // this makes a unit shell
242  struct make_unit{};
243  Shell(make_unit) :
244  alpha{0.0}, // exponent = 0
245  contr{{0, false, {1}}}, // contraction coefficient = 1
246  O{{0, 0, 0}}, // placed at origin
247  max_ln_coeff{0} {
248  }
249 
253  void renorm() {
254  using libint2::math::df_Kminus1;
255  using std::pow;
256  const auto sqrt_Pi_cubed = real_t{5.56832799683170784528481798212};
257  const auto np = nprim();
258  for(auto& c: contr) {
259  assert(c.l <= 15); // due to df_Kminus1[] a 64-bit integer type; kinda ridiculous restriction anyway
260  for(auto p=0ul; p!=np; ++p) {
261  assert(alpha[p] >= 0);
262  if (alpha[p] != 0) {
263  const auto two_alpha = 2 * alpha[p];
264  const auto two_alpha_to_am32 = pow(two_alpha,c.l+1) * sqrt(two_alpha);
265  const auto normalization_factor = sqrt(pow(2,c.l) * two_alpha_to_am32/(sqrt_Pi_cubed * df_Kminus1[2*c.l] ));
266 
267  c.coeff[p] *= normalization_factor;
268  }
269  }
270 
271  // need to force normalization to unity?
273  // compute the self-overlap of the , scale coefficients by its inverse square root
274  double norm{0};
275  for(auto p=0ul; p!=np; ++p) {
276  for(decltype(p) q=0ul; q<=p; ++q) {
277  auto gamma = alpha[p] + alpha[q];
278  norm += (p==q ? 1 : 2) * df_Kminus1[2*c.l] * sqrt_Pi_cubed * c.coeff[p] * c.coeff[q] /
279  (pow(2,c.l) * pow(gamma,c.l+1) * sqrt(gamma));
280  }
281  }
282  auto normalization_factor = 1 / sqrt(norm);
283  for(auto p=0ul; p!=np; ++p) {
284  c.coeff[p] *= normalization_factor;
285  }
286  }
287 
288  }
289 
290  // update max log coefficients
291  max_ln_coeff.resize(np);
292  for(auto p=0ul; p!=np; ++p) {
293  real_t max_ln_c = - std::numeric_limits<real_t>::max();
294  for(auto& c: contr) {
295  max_ln_c = std::max(max_ln_c, std::log(std::abs(c.coeff[p])));
296  }
297  max_ln_coeff[p] = max_ln_c;
298  }
299  }
300 
301  };
302 
303  inline std::ostream& operator<<(std::ostream& os, const Shell& sh) {
304  os << "Shell:( O={" << sh.O[0] << "," << sh.O[1] << "," << sh.O[2] << "}" << std::endl;
305  os << " ";
306  for(const auto& c: sh.contr) {
307  os << " {l=" << c.l << ",sph=" << c.pure << "}";
308  }
309  os << std::endl;
310 
311  for(auto i=0ul; i<sh.alpha.size(); ++i) {
312  os << " " << sh.alpha[i];
313  for(const auto& c: sh.contr) {
314  os << " " << c.coeff.at(i);
315  }
316  os << std::endl;
317  }
318 
319  return os;
320  }
321 
323  struct ShellPair {
324  typedef Shell::real_t real_t;
325 
326  struct PrimPairData {
327  real_t P[3];
328  real_t K;
329  real_t one_over_gamma;
330  real_t scr;
331  int p1;
332  int p2;
333  };
334 
335  std::vector<PrimPairData> primpairs;
336  real_t AB[3];
337 
338  ShellPair() : primpairs() { for(int i=0; i!=3; ++i) AB[i] = 0.; }
339 
340  ShellPair(size_t max_nprim) : primpairs() {
341  primpairs.reserve(max_nprim*max_nprim);
342  for(int i=0; i!=3; ++i) AB[i] = 0.;
343  }
344  template <typename Real> ShellPair(const Shell& s1, const Shell& s2, Real ln_prec) {
345  init(s1, s2, ln_prec);
346  }
347 
348  void resize(std::size_t max_nprim) {
349  const auto max_nprim2 = max_nprim * max_nprim;
350  if (max_nprim * max_nprim > primpairs.size())
351  primpairs.resize(max_nprim2);
352  }
353 
358  template <typename Real> void init(const Shell& s1, const Shell& s2, Real ln_prec) {
359 
360  primpairs.clear();
361 
362  const auto& A = s1.O;
363  const auto& B = s2.O;
364  real_t AB2 = 0.;
365  for(int i=0; i!=3; ++i) {
366  AB[i] = A[i] - B[i];
367  AB2 += AB[i]*AB[i];
368  }
369 
370  size_t c = 0;
371  for(size_t p1=0; p1!=s1.alpha.size(); ++p1) {
372  for(size_t p2=0; p2!=s2.alpha.size(); ++p2) {
373 
374  const auto& a1 = s1.alpha[p1];
375  const auto& a2 = s2.alpha[p2];
376  const auto gamma = a1 + a2;
377  const auto oogamma = 1 / gamma;
378 
379  const auto rho = a1 * a2 * oogamma;
380  const auto minus_rho_times_AB2 = -rho*AB2;
381  const auto screen_fac = minus_rho_times_AB2 + s1.max_ln_coeff[p1] + s2.max_ln_coeff[p2];
382  if (screen_fac < ln_prec)
383  continue;
384 
385  primpairs.resize(c+1);
386  PrimPairData& p = primpairs[c];
387  p.scr = screen_fac;
388  p.p1 = p1;
389  p.p2 = p2;
390  p.K = exp(minus_rho_times_AB2) * oogamma;
391  if (AB2 == 0.) { // this buys a bit more precision
392  p.P[0] = A[0];
393  p.P[1] = A[1];
394  p.P[2] = A[2];
395  } else {
396  p.P[0] = (a1 * A[0] + a2 * B[0]) * oogamma;
397  p.P[1] = (a1 * A[1] + a2 * B[1]) * oogamma;
398  p.P[2] = (a1 * A[2] + a2 * B[2]) * oogamma;
399  }
400  p.one_over_gamma = oogamma;
401 
402  ++c;
403  }
404  }
405  }
406 
407  };
408 
409 } // namespace libint2
410 
411 #endif /* _libint2_src_lib_libint_shell_h_ */
ShellPair pre-computes shell-pair data, primitive pairs are screened to finite precision.
Definition: shell.h:323
generally-contracted Solid-Harmonic/Cartesion Gaussian Shell
Definition: shell.h:81
static bool do_enforce_unit_normalization(defaultable_boolean flag=defaultable_boolean())
sets and/or reports whether the auto-renormalization to unity is set if called without arguments,...
Definition: shell.h:209
Defaults definitions for various parameters assumed by Libint.
Definition: algebra.cc:24
svector< Contraction > contr
contractions
Definition: shell.h:105
real_t P[3]
Definition: shell.h:327
void init(const Shell &s1, const Shell &s2, Real ln_prec)
initializes "expensive" primitive pair data; a pair of primitives with exponents located at whose m...
Definition: shell.h:358
std::array< real_t, 3 > O
origin
Definition: shell.h:106
svector< real_t > max_ln_coeff
maximum ln of (absolute) contraction coefficient for each primitive
Definition: shell.h:107
Array idential to C++0X arrays.
Definition: stdarray_bits.h:14
real_t coeff_normalized(size_t c, size_t p) const
Definition: shell.h:225
contracted Gaussian = angular momentum + sph/cart flag + contraction coefficients
Definition: shell.h:85
svector< real_t > alpha
exponents
Definition: shell.h:104
static const Shell & unit()
Definition: shell.h:218
Definition: shell.h:326