19 #ifndef _libint2_src_lib_libint_shell_h_ 20 #define _libint2_src_lib_libint_shell_h_ 22 #include <libint2/cxxstd.h> 23 #if LIBINT2_CPLUSPLUS_STD < 2011 24 # error "libint2/shell.h requires C++11 support" 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,
51 template <
typename Int> int64_t bc(Int i, Int j) {
52 assert(i < fac.size());
53 assert(j < fac.size());
55 return fac[i] / (fac[j] * fac[i-j]);
79 typedef double real_t;
85 std::vector<real_t> coeff;
87 return &other ==
this || (l == other.l && pure == other.pure && coeff == other.coeff);
90 return not this->operator==(other);
92 size_t cartesian_size()
const {
93 return (l + 1) * (l + 2) / 2;
96 return pure ? (2 * l + 1) : cartesian_size();
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)) {
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);
123 Shell(std::vector<real_t> _alpha,
124 std::vector<Contraction> _contr,
126 alpha(std::move(_alpha)),
127 contr(std::move(_contr)),
134 O = std::move(new_origin);
138 size_t cartesian_size()
const {
140 for(
const auto& c: contr) { s += c.cartesian_size(); }
143 size_t size()
const {
145 for(
const auto& c: contr) { s += c.size(); }
149 size_t ncontr()
const {
return contr.size(); }
150 size_t nprim()
const {
return alpha.size(); }
152 bool operator==(
const Shell& other)
const {
153 return &other ==
this || (O == other.
O && alpha == other.
alpha && contr == other.
contr);
155 bool operator!=(
const Shell& other)
const {
156 return not this->operator==(other);
159 static char am_symbol(
size_t l) {
160 static char lsymb[] =
"spdfghikmnoqrtuvwxyz";
164 static unsigned short am_symbol_to_l(
char am_symbol) {
165 const char AM_SYMBOL = ::toupper(am_symbol);
187 default:
throw "invalid angular momentum label";
192 typedef enum {false_value=0,true_value=1,default_value=2} value_t;
195 bool is_default()
const {
return value_ == default_value; }
196 operator bool()
const { assert(value_ != default_value);
return value_ == true_value; }
206 static bool result{
true};
207 if (not flag.is_default()) {
215 static const Shell unitshell{make_unit()};
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) {
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] ));
246 c.coeff[p] *= normalization_factor;
251 if (do_enforce_unit_normalization()) {
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));
261 auto normalization_factor = 1.0 / sqrt(norm);
262 for(
auto p=0; p!=np; ++p) {
263 c.coeff[p] *= normalization_factor;
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])));
276 max_ln_coeff[p] = max_ln_c;
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;
285 for(
const auto& c: sh.
contr) {
286 os <<
" {l=" << c.l <<
",sph=" << c.pure <<
"}";
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);
303 typedef Shell::real_t real_t;
308 real_t one_over_gamma;
314 std::vector<PrimPairData> primpairs;
317 ShellPair() : primpairs() {
for(
int i=0; i!=3; ++i) AB[i] = 0.; }
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.;
325 void init(
const Shell& s1,
const Shell& s2,
const real_t& ln_prec) {
329 const auto& A = s1.
O;
330 const auto& B = s2.
O;
332 for(
int i=0; i!=3; ++i) {
338 for(
size_t p1=0; p1!=s1.
alpha.size(); ++p1) {
339 for(
size_t p2=0; p2!=s2.
alpha.size(); ++p2) {
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;
346 const auto rho = a1 * a2 * oogamma;
347 const auto minus_rho_times_AB2 = -rho*AB2;
349 if (screen_fac < ln_prec)
352 primpairs.resize(c+1);
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;
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