| 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 |
|
*/ |
| 121 |
|
* @return the value of Rosen's B integral |
| 122 |
|
* @note N. Rosen, Phys. Rev., 38 (1931), 255 |
| 123 |
|
*/ |
| 124 |
< |
inline RealType RosenB(int n, RealType alpha) |
| 125 |
< |
{ |
| 126 |
< |
RealType TheSum, Term; |
| 126 |
< |
RealType RosenB_, PSinhRosenA, PCoshRosenA, PHyperRosenA; |
| 127 |
< |
int nu; |
| 128 |
< |
bool IsPositive; |
| 124 |
> |
inline RealType RosenB(int n, RealType alpha) { |
| 125 |
> |
RealType RosenB_; |
| 126 |
> |
|
| 127 |
|
if (alpha != 0.) |
| 128 |
|
{ |
| 129 |
< |
Term = 1.; |
| 130 |
< |
TheSum = 1.; |
| 133 |
< |
IsPositive = true; |
| 129 |
> |
RealType Term = 1.; |
| 130 |
> |
bool IsPositive = true; |
| 131 |
|
|
| 132 |
|
// These two expressions are (up to constant factors) equivalent |
| 133 |
|
// to computing the hyperbolic sine and cosine of a respectively |
| 134 |
|
// The series consists of adding up these terms in an alternating fashion |
| 135 |
< |
PSinhRosenA = exp(alpha) - exp(-alpha); |
| 136 |
< |
PCoshRosenA = -exp(alpha) - exp(-alpha); |
| 137 |
< |
TheSum = PSinhRosenA; |
| 135 |
> |
RealType PSinhRosenA = exp(alpha) - exp(-alpha); |
| 136 |
> |
RealType PCoshRosenA = -exp(alpha) - exp(-alpha); |
| 137 |
> |
RealType TheSum = PSinhRosenA; |
| 138 |
> |
RealType PHyperRosenA; |
| 139 |
|
for (unsigned nu=1; nu<=n; nu++) |
| 140 |
|
{ |
| 141 |
|
if (IsPositive) |
| 203 |
|
*/ |
| 204 |
|
inline RealType sSTOCoulInt(RealType a, RealType b, int m, int n, RealType R) |
| 205 |
|
{ |
| 208 |
– |
int nu, p; |
| 206 |
|
RealType x, K2; |
| 207 |
|
RealType Factor1, Factor2, Term, OneElectronTerm; |
| 208 |
|
RealType eps, epsi; |
| 216 |
|
|
| 217 |
|
// First compute the two-electron component |
| 218 |
|
RealType sSTOCoulInt_ = 0.; |
| 219 |
< |
if (x == 0.) // Pathological case |
| 219 |
> |
if (std::fabs(x) < OpenMD::NumericConstant::epsilon) // Pathological case |
| 220 |
|
{ |
| 221 |
< |
if ((a==b) && (m==n)) |
| 222 |
< |
{ |
| 223 |
< |
for (int nu=0; nu<=2*n-1; nu++) |
| 224 |
< |
{ |
| 225 |
< |
K2 = 0.; |
| 226 |
< |
for (unsigned p=0; p<=2*n+m; p++) K2 += 1. / fact[p]; |
| 227 |
< |
sSTOCoulInt_ += K2 * fact[2*n+m] / fact[m]; |
| 228 |
< |
} |
| 229 |
< |
sSTOCoulInt_ = 2 * a / (n * fact[2*n]) * sSTOCoulInt_; |
| 230 |
< |
} |
| 231 |
< |
else |
| 232 |
< |
{ |
| 233 |
< |
// Not implemented |
| 234 |
< |
cerr << "ERROR, sSTOCoulInt cannot compute from arguments\n"; |
| 235 |
< |
cerr << "a = " << a << "\tb = " << b << "\tm = " << m <<"\tn = " << n << "\t R = " << R << "\n"; |
| 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 |
< |
} |
| 221 |
> |
|
| 222 |
> |
// This solution for the one-center coulomb integrals comes from |
| 223 |
> |
// Yoshiyuki Hase, Computers & Chemistry 9(4), pp. 285-287 (1985). |
| 224 |
> |
|
| 225 |
> |
RealType Term1 = fact[2*m - 1] / pow(2*a, 2*m); |
| 226 |
> |
RealType Term2 = 0.; |
| 227 |
> |
for (int nu = 1; nu <= 2*n; nu++) { |
| 228 |
> |
Term2 += nu * pow(2*b, 2*n - nu) * fact[2*(m+n)-nu-1] / |
| 229 |
> |
(fact[2*n-nu]*2*n * pow(2*(a+b), 2*(m+n)-nu)); |
| 230 |
> |
} |
| 231 |
> |
sSTOCoulInt_ = pow(2*a, 2*m+1) * (Term1 - Term2) / fact[2*m]; |
| 232 |
> |
|
| 233 |
> |
// Original QTPIE code for the one-center case is below. Doesn't |
| 234 |
> |
// appear to generate the correct one-center results. |
| 235 |
> |
// |
| 236 |
> |
// if ((a==b) && (m==n)) |
| 237 |
> |
// { |
| 238 |
> |
// for (int nu=0; nu<=2*n-1; nu++) |
| 239 |
> |
// { |
| 240 |
> |
// K2 = 0.; |
| 241 |
> |
// for (unsigned p=0; p<=2*n+m; p++) K2 += 1. / fact[p]; |
| 242 |
> |
// sSTOCoulInt_ += K2 * fact[2*n+m] / fact[m]; |
| 243 |
> |
// } |
| 244 |
> |
// sSTOCoulInt_ = 2 * a / (n * fact[2*n]) * sSTOCoulInt_; |
| 245 |
> |
// } |
| 246 |
> |
// else |
| 247 |
> |
// { |
| 248 |
> |
// // Not implemented |
| 249 |
> |
// printf("ERROR, sSTOCoulInt cannot compute from arguments\n"); |
| 250 |
> |
// printf("a = %lf b = %lf m = %d n = %d R = %lf\n",a, b, m, n, R); |
| 251 |
> |
// exit(0); |
| 252 |
> |
// } |
| 253 |
> |
|
| 254 |
|
} |
| 255 |
|
else |
| 256 |
|
{ |
| 394 |
|
*/ |
| 395 |
|
inline RealType sSTOCoulIntGrad(RealType a, RealType b, int m, int n, RealType R) |
| 396 |
|
{ |
| 397 |
< |
RealType x, y, z, K2, TheSum; |
| 397 |
> |
RealType x; |
| 398 |
|
// x is the argument of the auxiliary RosenA and RosenB functions |
| 399 |
|
x = 2. * a * R; |
| 400 |
|
|
| 407 |
|
} |
| 408 |
|
else |
| 409 |
|
{ |
| 410 |
+ |
RealType K2, TheSum; |
| 411 |
|
if (a == b) |
| 412 |
|
{ |
| 413 |
|
TheSum = 0.; |
| 433 |
|
{ |
| 434 |
|
// Slater exponents are different |
| 435 |
|
// First calculate some useful arguments |
| 436 |
< |
y = R*(a+b); |
| 437 |
< |
z = R*(a-b); |
| 436 |
> |
RealType y = R*(a+b); |
| 437 |
> |
RealType z = R*(a-b); |
| 438 |
|
TheSum = 0.; |
| 439 |
|
for (int nu=0; nu<=2*n-1; nu++) |
| 440 |
|
{ |
| 479 |
|
* @note Derived in QTPIE research notes, May 15 2007 |
| 480 |
|
*/ |
| 481 |
|
inline RealType sSTOOvIntGrad(RealType a, RealType b, int m, int n, RealType R) |
| 482 |
< |
{ |
| 471 |
< |
RealType w, x, y, z, TheSum; |
| 472 |
< |
|
| 482 |
> |
{ |
| 483 |
|
// Calculate first term |
| 484 |
|
RealType sSTOOvIntGrad_ = (m+n+1.)/R * sSTOOvInt(a, b, m, n, R); |
| 485 |
|
|
| 486 |
|
// Calculate remaining terms; answers depend on exponents |
| 487 |
< |
TheSum = 0.; |
| 488 |
< |
x = a * R; |
| 487 |
> |
RealType TheSum = 0.; |
| 488 |
> |
RealType x = a * R; |
| 489 |
|
if (a == b) |
| 490 |
|
{ |
| 491 |
|
for (int q=0; q<=(m+n)/2; q++) |
| 494 |
|
} |
| 495 |
|
else |
| 496 |
|
{ |
| 497 |
< |
w = b*R; |
| 498 |
< |
y = 0.5*R*(a+b); |
| 499 |
< |
z = 0.5*R*(a-b); |
| 497 |
> |
RealType w = b*R; |
| 498 |
> |
RealType y = 0.5*R*(a+b); |
| 499 |
> |
RealType z = 0.5*R*(a-b); |
| 500 |
|
for (int q=0; q<m+n; q++) |
| 501 |
|
TheSum = TheSum + RosenD(m,n,q) * |
| 502 |
|
((a-b)*RosenB(q+1,z)*RosenA(m+n-q ,y) |
| 508 |
|
|
| 509 |
|
/** |
| 510 |
|
* @brief Calculates a Slater-type orbital exponent based on the hardness parameters |
| 511 |
< |
* @param Hardness: chemical hardness in atomic units |
| 511 |
> |
* @param hardness: chemical hardness in atomic units |
| 512 |
|
* @param n: principal quantum number |
| 513 |
|
* @note Modified for use with OpenMD by Gezelter and Michalka. |
| 514 |
|
*/ |
