21 #ifndef _libint2_src_lib_libint_shell_h_ 22 #define _libint2_src_lib_libint_shell_h_ 24 #include <libint2/util/cxxstd.h> 25 #if LIBINT2_CPLUSPLUS_STD < 2011 26 # error "libint2/shell.h requires C++11 support" 37 #include <libint2/util/small_vector.h> 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,
54 template <
typename Int> int64_t bc(Int i, Int j) {
55 assert(i < Int(fac.size()));
56 assert(j < Int(fac.size()));
58 return fac[i] / (fac[j] * fac[i-j]);
82 typedef double real_t;
88 svector<real_t> coeff;
91 return &other ==
this || (l == other.l && pure == other.pure && coeff == other.coeff);
94 return not this->operator==(other);
96 size_t cartesian_size()
const {
97 return (l + 1) * (l + 2) / 2;
100 return pure ? (2 * l + 1) : cartesian_size();
113 alpha(std::move(other.alpha)),
114 contr(std::move(other.contr)),
115 O(std::move(other.O)),
121 alpha = std::move(other.alpha);
122 contr = std::move(other.contr);
123 O = std::move(other.O);
127 Shell(svector<real_t> _alpha,
128 svector<Contraction> _contr,
130 alpha(std::move(_alpha)),
131 contr(std::move(_contr)),
138 O = std::move(new_origin);
142 size_t cartesian_size()
const {
144 for(
const auto& c:
contr) { s += c.cartesian_size(); }
147 size_t size()
const {
149 for(
const auto& c:
contr) { s += c.size(); }
153 size_t ncontr()
const {
return contr.size(); }
154 size_t nprim()
const {
return alpha.size(); }
156 bool operator==(
const Shell& other)
const {
157 return &other ==
this || (
O == other.O &&
alpha == other.alpha &&
contr == other.contr);
159 bool operator!=(
const Shell& other)
const {
160 return not this->operator==(other);
163 static char am_symbol(
size_t l) {
164 static char lsymb[] =
"spdfghikmnoqrtuvwxyz";
168 static unsigned short am_symbol_to_l(
char am_symbol) {
169 const char AM_SYMBOL = ::toupper(am_symbol);
191 default:
throw "invalid angular momentum label";
196 typedef enum {false_value=0,true_value=1,default_value=2} value_t;
199 bool is_default()
const {
return value_ == default_value; }
200 operator bool()
const { assert(value_ != default_value);
return value_ == true_value; }
210 static bool result{
true};
211 if (not flag.is_default()) {
219 static const Shell unitshell{make_unit()};
226 const auto alpha = this->alpha.at(p);
227 assert(
alpha >= 0.0);
228 const auto l =
contr.at(c).l;
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;
245 contr{{0,
false, {1}}},
254 using libint2::math::df_Kminus1;
256 const auto sqrt_Pi_cubed = real_t{5.56832799683170784528481798212};
257 const auto np = nprim();
258 for(
auto& c:
contr) {
260 for(
auto p=0ul; p!=np; ++p) {
261 assert(
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] ));
267 c.coeff[p] *= normalization_factor;
275 for(
auto p=0ul; p!=np; ++p) {
276 for(decltype(p) q=0ul; q<=p; ++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));
282 auto normalization_factor = 1 / sqrt(norm);
283 for(
auto p=0ul; p!=np; ++p) {
284 c.coeff[p] *= normalization_factor;
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])));
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;
306 for(
const auto& c: sh.contr) {
307 os <<
" {l=" << c.l <<
",sph=" << c.pure <<
"}";
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);
324 typedef Shell::real_t real_t;
329 real_t one_over_gamma;
335 std::vector<PrimPairData> primpairs;
338 ShellPair() : primpairs() {
for(
int i=0; i!=3; ++i) AB[i] = 0.; }
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.;
344 template <
typename Real> ShellPair(
const Shell& s1,
const Shell& s2, Real ln_prec) {
345 init(s1, s2, ln_prec);
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);
358 template <
typename Real>
void init(
const Shell& s1,
const Shell& s2, Real ln_prec) {
362 const auto& A = s1.
O;
363 const auto& B = s2.
O;
365 for(
int i=0; i!=3; ++i) {
371 for(
size_t p1=0; p1!=s1.
alpha.size(); ++p1) {
372 for(
size_t p2=0; p2!=s2.
alpha.size(); ++p2) {
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;
379 const auto rho = a1 * a2 * oogamma;
380 const auto minus_rho_times_AB2 = -rho*AB2;
382 if (screen_fac < ln_prec)
385 primpairs.resize(c+1);
390 p.K = exp(minus_rho_times_AB2) * oogamma;
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;
400 p.one_over_gamma = oogamma;
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