| 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 |
| 85 |
|
* \frac{\alpha^\nu}{\nu!} |
| 86 |
|
* \f] |
| 87 |
|
* @param n - principal quantum number |
| 88 |
< |
* @param alpha - Slater exponent |
| 88 |
> |
* @param a - Slater exponent |
| 89 |
|
* @return the value of Rosen's A integral |
| 90 |
|
* @note N. Rosen, Phys. Rev., 38 (1931), 255 |
| 91 |
|
*/ |
| 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; |
| 125 |
|
{ |
| 126 |
|
RealType TheSum, Term; |
| 127 |
|
RealType RosenB_, PSinhRosenA, PCoshRosenA, PHyperRosenA; |
| 125 |
– |
int nu; |
| 128 |
|
bool IsPositive; |
| 129 |
|
if (alpha != 0.) |
| 130 |
|
{ |
| 174 |
|
* @note N. Rosen, Phys. Rev., 38 (1931), 255 |
| 175 |
|
*/ |
| 176 |
|
inline RealType RosenD(int m, int n, int p) |
| 175 |
– |
// [?ERROR?] can't return int |
| 177 |
|
{ |
| 177 |
– |
printf("RosenD\n"); |
| 178 |
– |
exit(0); |
| 178 |
|
if (m+n+p > maxFact) |
| 179 |
|
{ |
| 180 |
|
printf("Error, arguments exceed maximum factorial computed %d > %d\n", m+n+p, maxFact); |
| 205 |
|
*/ |
| 206 |
|
inline RealType sSTOCoulInt(RealType a, RealType b, int m, int n, RealType R) |
| 207 |
|
{ |
| 209 |
– |
int nu, p; |
| 208 |
|
RealType x, K2; |
| 209 |
|
RealType Factor1, Factor2, Term, OneElectronTerm; |
| 210 |
|
RealType eps, epsi; |
| 218 |
|
|
| 219 |
|
// First compute the two-electron component |
| 220 |
|
RealType sSTOCoulInt_ = 0.; |
| 221 |
< |
if (x == 0.) // Pathological case |
| 221 |
> |
if (std::fabs(x) < OpenMD::NumericConstant::epsilon) // Pathological case |
| 222 |
|
{ |
| 223 |
< |
if ((a==b) && (m==n)) |
| 224 |
< |
{ |
| 225 |
< |
for (int nu=0; nu<=2*n-1; nu++) |
| 226 |
< |
{ |
| 227 |
< |
K2 = 0.; |
| 228 |
< |
for (unsigned p=0; p<=2*n+m; p++) K2 += 1. / fact[p]; |
| 229 |
< |
sSTOCoulInt_ += K2 * fact[2*n+m] / fact[m]; |
| 230 |
< |
} |
| 231 |
< |
sSTOCoulInt_ = 2 * a / (n * fact[2*n]) * sSTOCoulInt_; |
| 232 |
< |
} |
| 233 |
< |
else |
| 234 |
< |
{ |
| 235 |
< |
// Not implemented |
| 236 |
< |
printf("ERROR, sSTOCoulInt cannot compute from arguments\n"); |
| 237 |
< |
printf("a = %lf b = %lf m = %d n = %d R = %lf\n",a, b, m, n, R); |
| 238 |
< |
exit(0); |
| 239 |
< |
} |
| 223 |
> |
|
| 224 |
> |
// This solution for the one-center coulomb integrals comes from |
| 225 |
> |
// Yoshiyuki Hase, Computers & Chemistry 9(4), pp. 285-287 (1985). |
| 226 |
> |
|
| 227 |
> |
RealType Term1 = fact[2*m - 1] / pow(2*a, 2*m); |
| 228 |
> |
RealType Term2 = 0.; |
| 229 |
> |
for (int nu = 1; nu <= 2*n; nu++) { |
| 230 |
> |
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)); |
| 231 |
> |
} |
| 232 |
> |
sSTOCoulInt_ = pow(2*a, 2*m+1) * (Term1 - Term2) / fact[2*m]; |
| 233 |
> |
|
| 234 |
> |
// Original QTPIE code for the one-center case is below. Doesn't |
| 235 |
> |
// appear to generate the correct one-center results. |
| 236 |
> |
// |
| 237 |
> |
// if ((a==b) && (m==n)) |
| 238 |
> |
// { |
| 239 |
> |
// for (int nu=0; nu<=2*n-1; nu++) |
| 240 |
> |
// { |
| 241 |
> |
// K2 = 0.; |
| 242 |
> |
// for (unsigned p=0; p<=2*n+m; p++) K2 += 1. / fact[p]; |
| 243 |
> |
// sSTOCoulInt_ += K2 * fact[2*n+m] / fact[m]; |
| 244 |
> |
// } |
| 245 |
> |
// sSTOCoulInt_ = 2 * a / (n * fact[2*n]) * sSTOCoulInt_; |
| 246 |
> |
// } |
| 247 |
> |
// else |
| 248 |
> |
// { |
| 249 |
> |
// // Not implemented |
| 250 |
> |
// printf("ERROR, sSTOCoulInt cannot compute from arguments\n"); |
| 251 |
> |
// printf("a = %lf b = %lf m = %d n = %d R = %lf\n",a, b, m, n, R); |
| 252 |
> |
// exit(0); |
| 253 |
> |
// } |
| 254 |
> |
|
| 255 |
|
} |
| 256 |
|
else |
| 257 |
|
{ |
| 510 |
|
|
| 511 |
|
/** |
| 512 |
|
* @brief Calculates a Slater-type orbital exponent based on the hardness parameters |
| 513 |
< |
* @param Hardness: chemical hardness in atomic units |
| 513 |
> |
* @param hardness: chemical hardness in atomic units |
| 514 |
|
* @param n: principal quantum number |
| 515 |
|
* @note Modified for use with OpenMD by Gezelter and Michalka. |
| 516 |
|
*/ |
