| 61 |
|
|
| 62 |
|
#include "config.h" |
| 63 |
|
#include <cmath> |
| 64 |
+ |
#include <cstdlib> |
| 65 |
+ |
#include <iostream> |
| 66 |
|
#include "math/Factorials.hpp" |
| 67 |
+ |
#include "utils/NumericConstant.hpp" |
| 68 |
|
|
| 69 |
|
#ifndef NONBONDED_SLATERINTEGRALS_HPP |
| 70 |
|
#define NONBONDED_SLATERINTEGRALS_HPP |
| 96 |
|
{ |
| 97 |
|
RealType Term = 1.; |
| 98 |
|
RosenA_ = Term; |
| 99 |
< |
for (unsigned nu=1; nu<=n; nu++) |
| 99 |
> |
for (int nu=1; nu<=n; nu++) |
| 100 |
|
{ |
| 101 |
|
Term *= a / nu; |
| 102 |
|
RosenA_ += Term; |
| 175 |
|
* @note N. Rosen, Phys. Rev., 38 (1931), 255 |
| 176 |
|
*/ |
| 177 |
|
inline RealType RosenD(int m, int n, int p) |
| 175 |
– |
// [?ERROR?] can't return int |
| 178 |
|
{ |
| 177 |
– |
printf("RosenD\n"); |
| 178 |
– |
exit(0); |
| 179 |
|
if (m+n+p > maxFact) |
| 180 |
|
{ |
| 181 |
|
printf("Error, arguments exceed maximum factorial computed %d > %d\n", m+n+p, maxFact); |
| 220 |
|
|
| 221 |
|
// First compute the two-electron component |
| 222 |
|
RealType sSTOCoulInt_ = 0.; |
| 223 |
< |
if (x == 0.) // Pathological case |
| 223 |
> |
if (std::fabs(x) < OpenMD::NumericConstant::epsilon) // Pathological case |
| 224 |
|
{ |
| 225 |
< |
if ((a==b) && (m==n)) |
| 226 |
< |
{ |
| 227 |
< |
for (int nu=0; nu<=2*n-1; nu++) |
| 228 |
< |
{ |
| 229 |
< |
K2 = 0.; |
| 230 |
< |
for (unsigned p=0; p<=2*n+m; p++) K2 += 1. / fact[p]; |
| 231 |
< |
sSTOCoulInt_ += K2 * fact[2*n+m] / fact[m]; |
| 232 |
< |
} |
| 233 |
< |
sSTOCoulInt_ = 2 * a / (n * fact[2*n]) * sSTOCoulInt_; |
| 234 |
< |
} |
| 235 |
< |
else |
| 236 |
< |
{ |
| 237 |
< |
// Not implemented |
| 238 |
< |
printf("ERROR, sSTOCoulInt cannot compute from arguments\n"); |
| 239 |
< |
printf("a = %lf b = %lf m = %d n = %d R = %lf\n",a, b, m, n, R); |
| 240 |
< |
exit(0); |
| 241 |
< |
} |
| 225 |
> |
|
| 226 |
> |
// This solution for the one-center coulomb integrals comes from |
| 227 |
> |
// Yoshiyuki Hase, Computers & Chemistry 9(4), pp. 285-287 (1985). |
| 228 |
> |
|
| 229 |
> |
RealType Term1 = fact[2*m - 1] / pow(2*a, 2*m); |
| 230 |
> |
RealType Term2 = 0.; |
| 231 |
> |
for (int nu = 1; nu <= 2*n; nu++) { |
| 232 |
> |
Term2 += nu * pow(2*b, 2*n - nu) * fact[2*(m+n)-nu-1] / (fact[2*n-nu]*2*n * pow(2*(a+b), 2*(m+n)-nu)); |
| 233 |
> |
} |
| 234 |
> |
sSTOCoulInt_ = pow(2*a, 2*m+1) * (Term1 - Term2) / fact[2*m]; |
| 235 |
> |
|
| 236 |
> |
// Original QTPIE code for the one-center case is below. Doesn't |
| 237 |
> |
// appear to generate the correct one-center results. |
| 238 |
> |
// |
| 239 |
> |
// if ((a==b) && (m==n)) |
| 240 |
> |
// { |
| 241 |
> |
// for (int nu=0; nu<=2*n-1; nu++) |
| 242 |
> |
// { |
| 243 |
> |
// K2 = 0.; |
| 244 |
> |
// for (unsigned p=0; p<=2*n+m; p++) K2 += 1. / fact[p]; |
| 245 |
> |
// sSTOCoulInt_ += K2 * fact[2*n+m] / fact[m]; |
| 246 |
> |
// } |
| 247 |
> |
// sSTOCoulInt_ = 2 * a / (n * fact[2*n]) * sSTOCoulInt_; |
| 248 |
> |
// } |
| 249 |
> |
// else |
| 250 |
> |
// { |
| 251 |
> |
// // Not implemented |
| 252 |
> |
// printf("ERROR, sSTOCoulInt cannot compute from arguments\n"); |
| 253 |
> |
// printf("a = %lf b = %lf m = %d n = %d R = %lf\n",a, b, m, n, R); |
| 254 |
> |
// exit(0); |
| 255 |
> |
// } |
| 256 |
> |
|
| 257 |
|
} |
| 258 |
|
else |
| 259 |
|
{ |
