LIBINT  2.1.0-stable
shell.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 Library General Public License, version 2,
7  * as published by the Free Software Foundation.
8  *
9  * This program is distributed in the hope that it will be useful,
10  * but WITHOUT ANY WARRANTY; without even the implied warranty of
11  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12  * GNU General Public License for more details.
13  *
14  * You should have received a copy of the GNU Library General Public License
15  * along with this program. If not, see http://www.gnu.org/licenses/.
16  *
17  */
18 
19 #ifndef _libint2_src_lib_libint_shell_h_
20 #define _libint2_src_lib_libint_shell_h_
21 
22 #include <libint2/cxxstd.h>
23 #if LIBINT2_CPLUSPLUS_STD < 2011
24 # error "libint2/shell.h requires C++11 support"
25 #endif
26 
27 
28 #include <iostream>
29 #include <array>
30 #include <vector>
31 #include <cassert>
32 #include <cmath>
33 
34 #include <libint2.h>
35 
36 namespace libint2 {
37 
38  namespace math {
40  static constexpr std::array<int64_t,21> fac = {{1LL, 1LL, 2LL, 6LL, 24LL, 120LL, 720LL, 5040LL, 40320LL, 362880LL, 3628800LL, 39916800LL,
41  479001600LL, 6227020800LL, 87178291200LL, 1307674368000LL, 20922789888000LL,
42  355687428096000LL, 6402373705728000LL, 121645100408832000LL,
43  2432902008176640000LL}};
45  static constexpr std::array<int64_t,31> df_Kminus1 = {{1LL, 1LL, 1LL, 2LL, 3LL, 8LL, 15LL, 48LL, 105LL, 384LL, 945LL, 3840LL, 10395LL, 46080LL, 135135LL,
46  645120LL, 2027025LL, 10321920LL, 34459425LL, 185794560LL, 654729075LL,
47  3715891200LL, 13749310575LL, 81749606400LL, 316234143225LL, 1961990553600LL,
48  7905853580625LL, 51011754393600LL, 213458046676875LL, 1428329123020800LL,
49  6190283353629375LL}};
51  template <typename Int> int64_t bc(Int i, Int j) {
52  assert(i < fac.size());
53  assert(j < fac.size());
54  assert(i >= j);
55  return fac[i] / (fac[j] * fac[i-j]);
56  }
57  }
58 
60 
78  struct Shell {
79  typedef double real_t;
80 
82  struct Contraction {
83  int l;
84  bool pure;
85  std::vector<real_t> coeff;
86  bool operator==(const Contraction& other) const {
87  return &other == this || (l == other.l && pure == other.pure && coeff == other.coeff);
88  }
89  bool operator!=(const Contraction& other) const {
90  return not this->operator==(other);
91  }
92  size_t cartesian_size() const {
93  return (l + 1) * (l + 2) / 2;
94  }
95  size_t size() const {
96  return pure ? (2 * l + 1) : cartesian_size();
97  }
98  };
99 
100  std::vector<real_t> alpha;
101  std::vector<Contraction> contr;
103  std::vector<real_t> max_ln_coeff;
104 
105  Shell() = default;
106  Shell(const Shell&) = default;
107  // intel does not support "move ctor = default"
108  Shell(Shell&& other) :
109  alpha(std::move(other.alpha)),
110  contr(std::move(other.contr)),
111  O(std::move(other.O)),
112  max_ln_coeff(std::move(other.max_ln_coeff)) {
113  }
114  Shell& operator=(const Shell&) = default;
115  // intel does not support "move asgnmt = default"
116  Shell& operator=(Shell&& other) {
117  alpha = std::move(other.alpha);
118  contr = std::move(other.contr);
119  O = std::move(other.O);
120  max_ln_coeff = std::move(other.max_ln_coeff);
121  return *this;
122  }
123  Shell(std::vector<real_t> _alpha,
124  std::vector<Contraction> _contr,
126  alpha(std::move(_alpha)),
127  contr(std::move(_contr)),
128  O(std::move(_O)) {
129  // embed normalization factors into contraction coefficients
130  renorm();
131  }
132 
133  Shell& move(std::array<real_t, 3> new_origin) {
134  O = std::move(new_origin);
135  return *this;
136  }
137 
138  size_t cartesian_size() const {
139  size_t s = 0;
140  for(const auto& c: contr) { s += c.cartesian_size(); }
141  return s;
142  }
143  size_t size() const {
144  size_t s = 0;
145  for(const auto& c: contr) { s += c.size(); }
146  return s;
147  }
148 
149  size_t ncontr() const { return contr.size(); }
150  size_t nprim() const { return alpha.size(); }
151 
152  bool operator==(const Shell& other) const {
153  return &other == this || (O == other.O && alpha == other.alpha && contr == other.contr);
154  }
155  bool operator!=(const Shell& other) const {
156  return not this->operator==(other);
157  }
158 
159  static char am_symbol(size_t l) {
160  static char lsymb[] = "spdfghikmnoqrtuvwxyz";
161  assert(l<=19);
162  return lsymb[l];
163  }
164  static unsigned short am_symbol_to_l(char am_symbol) {
165  const char AM_SYMBOL = ::toupper(am_symbol);
166  switch (AM_SYMBOL) {
167  case 'S': return 0;
168  case 'P': return 1;
169  case 'D': return 2;
170  case 'F': return 3;
171  case 'G': return 4;
172  case 'H': return 5;
173  case 'I': return 6;
174  case 'K': return 7;
175  case 'M': return 8;
176  case 'N': return 9;
177  case 'O': return 10;
178  case 'Q': return 11;
179  case 'R': return 12;
180  case 'T': return 13;
181  case 'U': return 14;
182  case 'V': return 15;
183  case 'W': return 16;
184  case 'X': return 17;
185  case 'Y': return 18;
186  case 'Z': return 19;
187  default: throw "invalid angular momentum label";
188  }
189  }
190 
192  typedef enum {false_value=0,true_value=1,default_value=2} value_t;
193  defaultable_boolean() : value_(default_value) {}
194  defaultable_boolean(bool v) : value_(static_cast<value_t>(v?1:0)) {}
195  bool is_default() const { return value_ == default_value; }
196  operator bool() const { assert(value_ != default_value); return value_ == true_value; }
197  private:
198  value_t value_;
199  };
200 
206  static bool result{true};
207  if (not flag.is_default()) {
208  result = flag;
209  }
210  return result;
211  }
212 
214  static const Shell& unit() {
215  static const Shell unitshell{make_unit()};
216  return unitshell;
217  }
218 
219  private:
220 
221  // this makes a unit shell
222  struct make_unit{};
223  Shell(make_unit) :
224  alpha{0.0}, // exponent = 0
225  contr{Contraction{0, false, {1.0}}}, // contraction coefficient = 1
226  O{{0.0, 0.0, 0.0}}, // placed at origin
227  max_ln_coeff{0.0} {
228  }
229 
233  void renorm() {
234  using libint2::math::df_Kminus1;
235  const auto sqrt_Pi_cubed = real_t{5.56832799683170784528481798212};
236  const auto np = nprim();
237  for(auto& c: contr) {
238  assert(c.l <= 15); // due to df_Kminus1[] a 64-bit integer type; kinda ridiculous restriction anyway
239  for(auto p=0; p!=np; ++p) {
240  assert(alpha[p] >= 0.0);
241  if (alpha[p] != 0.) {
242  const auto two_alpha = 2 * alpha[p];
243  const auto two_alpha_to_am32 = pow(two_alpha,c.l+1) * sqrt(two_alpha);
244  const auto normalization_factor = sqrt(pow(2,c.l) * two_alpha_to_am32/(sqrt_Pi_cubed * df_Kminus1[2*c.l] ));
245 
246  c.coeff[p] *= normalization_factor;
247  }
248  }
249 
250  // need to force normalization to unity?
251  if (do_enforce_unit_normalization()) {
252  // compute the self-overlap of the , scale coefficients by its inverse square root
253  double norm{0};
254  for(auto p=0; p!=np; ++p) {
255  for(auto q=0; q<=p; ++q) {
256  auto gamma = alpha[p] + alpha[q];
257  norm += (p==q ? 1.0 : 2.0) * df_Kminus1[2*c.l] * sqrt_Pi_cubed * c.coeff[p] * c.coeff[q] /
258  (pow(2,c.l) * pow(gamma,c.l+1) * sqrt(gamma));
259  }
260  }
261  auto normalization_factor = 1.0 / sqrt(norm);
262  for(auto p=0; p!=np; ++p) {
263  c.coeff[p] *= normalization_factor;
264  }
265  }
266 
267  }
268 
269  // update max log coefficients
270  max_ln_coeff.resize(np);
271  for(auto p=0; p!=np; ++p) {
272  real_t max_ln_c = - std::numeric_limits<real_t>::max();
273  for(auto& c: contr) {
274  max_ln_c = std::max(max_ln_c, log(std::abs(c.coeff[p])));
275  }
276  max_ln_coeff[p] = max_ln_c;
277  }
278  }
279 
280  };
281 
282  inline std::ostream& operator<<(std::ostream& os, const Shell& sh) {
283  os << "Shell:( O={" << sh.O[0] << "," << sh.O[1] << "," << sh.O[2] << "}" << std::endl;
284  os << " ";
285  for(const auto& c: sh.contr) {
286  os << " {l=" << c.l << ",sph=" << c.pure << "}";
287  }
288  os << std::endl;
289 
290  for(auto i=0; i<sh.alpha.size(); ++i) {
291  os << " " << sh.alpha[i];
292  for(const auto& c: sh.contr) {
293  os << " " << c.coeff.at(i);
294  }
295  os << std::endl;
296  }
297 
298  return os;
299  }
300 
302  struct ShellPair {
303  typedef Shell::real_t real_t;
304 
305  struct PrimPairData {
306  real_t P[3];
307  real_t K;
308  real_t one_over_gamma;
309  real_t scr;
310  int p1;
311  int p2;
312  };
313 
314  std::vector<PrimPairData> primpairs;
315  real_t AB[3];
316 
317  ShellPair() : primpairs() { for(int i=0; i!=3; ++i) AB[i] = 0.; }
318 
319  ShellPair(size_t max_nprim) : primpairs() {
320  primpairs.reserve(max_nprim*max_nprim);
321  for(int i=0; i!=3; ++i) AB[i] = 0.;
322  }
323 
324  // initializes "expensive" primitive pair data; primitive pairs are
325  void init(const Shell& s1, const Shell& s2, const real_t& ln_prec) {
326 
327  primpairs.clear();
328 
329  const auto& A = s1.O;
330  const auto& B = s2.O;
331  real_t AB2 = 0.;
332  for(int i=0; i!=3; ++i) {
333  AB[i] = A[i] - B[i];
334  AB2 += AB[i]*AB[i];
335  }
336 
337  size_t c = 0;
338  for(size_t p1=0; p1!=s1.alpha.size(); ++p1) {
339  for(size_t p2=0; p2!=s2.alpha.size(); ++p2) {
340 
341  const auto& a1 = s1.alpha[p1];
342  const auto& a2 = s2.alpha[p2];
343  const auto gamma = a1 + a2;
344  const auto oogamma = 1.0 / gamma;
345 
346  const auto rho = a1 * a2 * oogamma;
347  const auto minus_rho_times_AB2 = -rho*AB2;
348  const auto screen_fac = minus_rho_times_AB2 + s1.max_ln_coeff[p1] + s2.max_ln_coeff[p2];
349  if (screen_fac < ln_prec)
350  continue;
351 
352  primpairs.resize(c+1);
353  PrimPairData& p = primpairs[c];
354  p.scr = screen_fac;
355  p.p1 = p1;
356  p.p2 = p2;
357  p.K = exp(minus_rho_times_AB2) * oogamma;
358  p.P[0] = (a1 * A[0] + a2 * B[0]) * oogamma;
359  p.P[1] = (a1 * A[1] + a2 * B[1]) * oogamma;
360  p.P[2] = (a1 * A[2] + a2 * B[2]) * oogamma;
361  p.one_over_gamma = oogamma;
362 
363  ++c;
364  }
365  }
366  }
367 
368  };
369 
370 } // namespace libint2
371 
372 #endif /* _libint2_src_lib_libint_shell_h_ */
std::vector< real_t > alpha
exponents
Definition: shell.h:100
ShellPair pre-computes shell-pair data, primitive pairs are screened to finite precision.
Definition: shell.h:302
generally-contracted Solid-Harmonic/Cartesion Gaussian Shell
Definition: shell.h:78
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:205
Defaults definitions for various parameters assumed by Libint.
Definition: algebra.cc:23
real_t P[3]
Definition: shell.h:306
Definition: stdarray.h:18
std::array< real_t, 3 > O
origin
Definition: shell.h:102
std::vector< real_t > max_ln_coeff
maximum ln of (absolute) contraction coefficient for each primitive
Definition: shell.h:103
Array idential to C++0X arrays.
Definition: stdarray_bits.h:14
contracted Gaussian = angular momentum + sph/cart flag + contraction coefficients ...
Definition: shell.h:82
static const Shell & unit()
Definition: shell.h:214
std::vector< Contraction > contr
contractions
Definition: shell.h:101
Definition: shell.h:305