OpenMD 3.0
Molecular Dynamics in the Open
Loading...
Searching...
No Matches
SlaterIntegrals.hpp
1/***************************************************************************
2 * This program is free software; you can redistribute it and/or modify *
3 * it under the terms of the GNU General Public License as published by *
4 * the Free Software Foundation; either version 3 of the License, or *
5 * (at your option) any later version. *
6 * *
7 * This program is distributed in the hope that it will be useful, *
8 * but WITHOUT ANY WARRANTY; without even the implied warranty of *
9 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
10 * GNU General Public License for more details. *
11 * *
12 * You should have received a copy of the GNU General Public License *
13 * along with this program; if not, see <http://www.gnu.org/licenses/>. *
14 ***************************************************************************/
15
16/***
17 * This file was imported from qtpie, found at http://code.google.com/p/qtpie
18 * and was modified minimally for use in OpenMD.
19 *
20 * No author attribution was found in the code, but it presumably is
21 * the work of J. Chen and Todd J. Martinez.
22 *
23 * QTPIE (charge transfer with polarization current equilibration) is
24 * a new charge model, similar to other charge models like QEq,
25 * fluc-q, EEM or ABEEM. Unlike other existing charge models, however,
26 * it is capable of describing both charge transfer and polarization
27 * phenomena. It is also unique for its ability to describe
28 * intermolecular charge transfer at reasonable computational cost.
29 *
30 * Good references to cite when using this code are:
31 *
32 * J. Chen and T. J. Martinez, "QTPIE: Charge transfer with
33 * polarization current equalization. A fluctuating charge model with
34 * correct asymptotics", Chemical Physics Letters 438 (2007),
35 * 315-320. DOI: 10.1016/j.cplett.2007.02.065. Erratum: ibid, 463
36 * (2008), 288. DOI: 10.1016/j.cplett.2008.08.060
37 *
38 * J. Chen, D. Hundertmark and T. J. Martinez, "A unified theoretical
39 * framework for fluctuating-charge models in atom-space and in
40 * bond-space", Journal of Chemical Physics 129 (2008), 214113. DOI:
41 * 10.1063/1.3021400.
42 *
43 * J. Chen and T. J. Martinez, "Charge conservation in
44 * electronegativity equalization and its implications for the
45 * electrostatic properties of fluctuating-charge models", Journal of
46 * Chemical Physics 131 (2009), 044114. DOI: 10.1063/1.3183167
47 *
48 * J. Chen and T. J. Martinez, "The dissociation catastrophe in
49 * fluctuating-charge models and its implications for the concept of
50 * atomic electronegativity", Progress in Theoretical Chemistry and
51 * Physics, to appear. arXiv:0812.1543
52 *
53 * J. Chen, "Theory and applications of fluctuating-charge models",
54 * PhD (chemical physics) thesis, University of Illinois at
55 * Urbana-Champaign, Department of Chemistry, 2009.
56 *
57 * J. Chen and T. J. Martinez, "Size-extensive polarizabilities with
58 * intermolecular charge transfer in a fluctuating-charge model", in
59 * preparation. arXiv:0812.1544
60 */
61
62#include <config.h>
63
64#include <algorithm>
65#include <cmath>
66#include <cstdlib>
67#include <iostream>
68#include <limits>
69
70#include "math/Factorials.hpp"
71#include "utils/Constants.hpp"
72
73#ifndef NONBONDED_SLATERINTEGRALS_HPP
74#define NONBONDED_SLATERINTEGRALS_HPP
75
76template<typename T>
77inline T sqr(T t) {
78 return t * t;
79}
80template<typename T>
81inline T mod(T x, T m) {
82 return x < 0 ? m - 1 - ((-x) - 1) % m : x % m;
83}
84
85// #include "parameters.h"
86
87/**
88 * @brief Computes Rosen's Guillimer-Zener function A
89 * Computes Rosen's A integral, an auxiliary quantity needed to
90 * compute integrals involving Slater-type orbitals of s symmetry.
91 * \f[
92 * A_n(\alpha) = \int_1^\infty x^n e^{-\alpha x}dx
93 * = \frac{n! e^{-\alpha}}{\alpha^{n+1}}\sum_{\nu=0}^n
94 * \frac{\alpha^\nu}{\nu!}
95 * \f]
96 * @param n - principal quantum number
97 * @param a - Slater exponent
98 * @return the value of Rosen's A integral
99 * @note N. Rosen, Phys. Rev., 38 (1931), 255
100 */
101inline RealType RosenA(int n, RealType a) {
102 RealType RosenA_ = 0.;
103 if (a != 0.) {
104 RealType Term = 1.;
105 RosenA_ = Term;
106 for (int nu = 1; nu <= n; nu++) {
107 Term *= a / nu;
108 RosenA_ += Term;
109 }
110 RosenA_ = (RosenA_ / Term) * (exp(-a) / a);
111 }
112 return RosenA_;
113}
114
115/**
116 * @brief Computes Rosen's Guillimer-Zener function B
117 * Computes Rosen's B integral, an auxiliary quantity needed to
118 * compute integrals involving Slater-type orbitals of s symmetry.
119 * \f[
120 * B_n(\alpha) = \int_{-1}^1 x^n e^{-\alpha x} dx
121 * = \frac{n!}{\alpha^{n+1}}
122 * \sum_{\nu=0}^n \frac{e^\alpha(-\alpha)^\nu
123 * - e^{-\alpha} \alpha^\nu}{\nu!}
124 * \f]
125 * @param n - principal quantum number
126 * @param alpha - Slater exponent
127 * @return the value of Rosen's B integral
128 * @note N. Rosen, Phys. Rev., 38 (1931), 255
129 */
130inline RealType RosenB(int n, RealType alpha) {
131 RealType RosenB_;
132
133 if (alpha != 0.) {
134 RealType Term = 1.;
135 bool IsPositive = true;
136
137 // These two expressions are (up to constant factors) equivalent
138 // to computing the hyperbolic sine and cosine of a respectively
139 // The series consists of adding up these terms in an alternating fashion
140 RealType PSinhRosenA = exp(alpha) - exp(-alpha);
141 RealType PCoshRosenA = -exp(alpha) - exp(-alpha);
142 RealType TheSum = PSinhRosenA;
143 RealType PHyperRosenA;
144 for (unsigned nu = 1; nu <= (unsigned)n; nu++) {
145 if (IsPositive) {
146 PHyperRosenA = PCoshRosenA;
147 IsPositive = false;
148 } else // term to add should be negative
149 {
150 PHyperRosenA = PSinhRosenA;
151 IsPositive = true;
152 }
153 Term *= alpha / nu;
154 TheSum += Term * PHyperRosenA;
155 }
156 RosenB_ = TheSum / (alpha * Term);
157 } else // pathological case of a=0
158 {
159 printf("WARNING, a = 0 in RosenB\n");
160 RosenB_ = (1. - pow(-1., n)) / (n + 1.);
161 }
162 return RosenB_;
163}
164
165/** @brief Computes Rosen's D combinatorial factor
166 * Computes Rosen's D factor, an auxiliary quantity needed to
167 * compute integrals involving Slater-type orbitals of s symmetry.
168 * \f[
169 * RosenD^{mn}_p = \sum_k (-1)^k \frac{m! n!}
170 * {(p-k)!(m-p+k)!(n-k)!k!}
171 * \f]
172 * @return the value of Rosen's D factor
173 * @note N. Rosen, Phys. Rev., 38 (1931), 255
174 */
175inline RealType RosenD(int m, int n, int p) {
176 if (m + n + p > maxFact) {
177 printf("Error, arguments exceed maximum factorial computed %d > %d\n",
178 m + n + p, maxFact);
179 ::exit(0);
180 }
181
182 RealType RosenD_ = 0;
183 for (int k = max(p - m, 0); k <= min(n, p); k++) {
184 if (mod(k, 2) == 0)
185 RosenD_ += (fact[m] / (fact[p - k] * fact[m - p + k])) *
186 (fact[n] / (fact[n - k] * fact[k]));
187 else
188 RosenD_ -= (fact[m] / (fact[p - k] * fact[m - p + k])) *
189 (fact[n] / (fact[n - k] * fact[k]));
190 }
191 return RosenD_;
192}
193
194/** @brief Computes Coulomb integral analytically over s-type STOs
195 * Computes the two-center Coulomb integral over Slater-type orbitals of s
196 * symmetry.
197 * @param a : Slater zeta exponent of first atom in inverse Bohr (au)
198 * @param b : Slater zeta exponent of second atom in inverse Bohr (au)
199 * @param m : principal quantum number of first atom
200 * @param n : principal quantum number of second atom
201 * @param R : internuclear distance in atomic units (bohr)
202 * @return value of the Coulomb potential energy integral
203 * @note N. Rosen, Phys. Rev., 38 (1931), 255
204 * @note In Rosen's paper, this integral is known as K2.
205 */
206inline RealType sSTOCoulInt(RealType a, RealType b, int m, int n, RealType R) {
207 RealType x, K2;
208 RealType Factor1, Factor2, Term, OneElectronTerm;
209 RealType eps, epsi;
210
211 // To speed up calculation, we terminate loop once contributions
212 // to integral fall below the bound, epsilon
213 RealType epsilon = 0.;
214
215 // x is the argument of the auxiliary RosenA and RosenB functions
216 x = 2. * a * R;
217
218 // First compute the two-electron component
219 RealType sSTOCoulInt_ = 0.;
220 if (std::fabs(x) <
221 std::numeric_limits<RealType>::epsilon()) // Pathological case
222 {
223 // This solution for the one-center coulomb integrals comes from
224 // Yoshiyuki Hase, Computers & Chemistry 9(4), pp. 285-287 (1985).
225
226 RealType Term1 = fact[2 * m - 1] / pow(2 * a, 2 * m);
227 RealType Term2 = 0.;
228 for (int nu = 1; nu <= 2 * n; nu++) {
229 Term2 += nu * pow(2 * b, 2 * n - nu) * fact[2 * (m + n) - nu - 1] /
230 (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 } else {
256 OneElectronTerm = 1. / R +
257 pow(x, 2 * m) / (fact[2 * m] * R) *
258 ((x - 2 * m) * RosenA(2 * m - 1, x) - exp(-x)) +
259 sSTOCoulInt_;
260 eps = epsilon / OneElectronTerm;
261 if (a == b) {
262 // Apply Rosen (48)
263 Factor1 = -a * pow(a * R, 2 * m) / (n * fact[2 * m]);
264 for (int nu = 0; nu <= 2 * n - 1; nu++) {
265 Factor2 = (2. * n - nu) / fact[nu] * pow(a * R, nu);
266 epsi = eps / fabs(Factor1 * Factor2);
267 K2 = 0.;
268 for (int p = 0; p <= m + (nu - 1) / 2; p++) {
269 Term = RosenD(2 * m - 1, nu, 2 * p) / (2. * p + 1.) *
270 RosenA(2 * m + nu - 1 - 2 * p, x);
271 K2 += Term;
272 if ((Term > 0) && (Term < epsi)) goto label1;
273 }
274 sSTOCoulInt_ += K2 * Factor2;
275 }
276 label1:
277 sSTOCoulInt_ *= Factor1;
278 } else {
279 Factor1 = -a * pow(a * R, 2 * m) / (2. * n * fact[2 * m]);
280 epsi = eps / fabs(Factor1);
281 if (b == 0.)
282 printf("WARNING: b = 0 in sSTOCoulInt\n");
283 else {
284 // Apply Rosen (54)
285 for (int nu = 0; nu <= 2 * n - 1; nu++) {
286 K2 = 0;
287 for (int p = 0; p <= 2 * m + nu - 1; p++)
288 K2 = K2 + RosenD(2 * m - 1, nu, p) * RosenB(p, R * (a - b)) *
289 RosenA(2 * m + nu - 1 - p, R * (a + b));
290 Term = K2 * (2 * n - nu) / fact[nu] * pow(b * R, nu);
291 sSTOCoulInt_ += Term;
292 if (fabs(Term) < epsi) goto label2;
293 }
294 label2:
295 sSTOCoulInt_ *= Factor1;
296 }
297 }
298 // Now add the one-electron term from Rosen (47) = Rosen (53)
299 sSTOCoulInt_ += OneElectronTerm;
300 }
301 return sSTOCoulInt_;
302}
303
304/**
305 * @brief Computes overlap integral analytically over s-type STOs
306 * Computes the overlap integral over two
307 * Slater-type orbitals of s symmetry.
308 * @param a : Slater zeta exponent of first atom in inverse Bohr (au)
309 * @param b : Slater zeta exponent of second atom in inverse Bohr (au)
310 * @param m : principal quantum number of first atom
311 * @param n : principal quantum number of second atom
312 * @param R : internuclear distance in atomic units (bohr)
313 * @return the value of the sSTOOvInt integral
314 * @note N. Rosen, Phys. Rev., 38 (1931), 255
315 * @note In the Rosen paper, this integral is known as I.
316 */
317inline RealType sSTOOvInt(RealType a, RealType b, int m, int n, RealType R) {
318 RealType Factor, Term, eps;
319
320 // To speed up calculation, we terminate loop once contributions
321 // to integral fall below the bound, epsilon
322 RealType epsilon = 0.;
323 RealType sSTOOvInt_ = 0.;
324
325 if (a == b) {
326 Factor = pow(a * R, m + n + 1) / sqrt(fact[2 * m] * fact[2 * n]);
327 eps = epsilon / fabs(Factor);
328 for (int q = 0; q <= (m + n) / 2; q++) {
329 Term = RosenD(m, n, 2 * q) / (2. * q + 1.) * RosenA(m + n - 2 * q, a * R);
330 sSTOOvInt_ += Term;
331 if (fabs(Term) < eps) exit(0);
332 }
333 sSTOOvInt_ *= Factor;
334 } else {
335 Factor = 0.5 * pow(a * R, m + 0.5) * pow(b * R, n + 0.5) /
336 sqrt(fact[2 * m] * fact[2 * n]);
337 eps = epsilon / fabs(Factor);
338 for (int q = 0; q <= m + n; q++) {
339 Term = RosenD(m, n, q) * RosenB(q, R / 2. * (a - b)) *
340 RosenA(m + n - q, R / 2. * (a + b));
341 sSTOOvInt_ += Term;
342 if (fabs(Term) < eps) exit(0);
343 }
344 sSTOOvInt_ *= Factor;
345 }
346 return sSTOOvInt_;
347}
348
349/**
350 * @brief Computes kinetic energy integral analytically over s-type STOs
351 * Computes the overlap integral over two Slater-type orbitals of s symmetry.
352 * @param a : Slater zeta exponent of first atom in inverse Bohr (au)
353 * @param b : Slater zeta exponent of second atom in inverse Bohr (au)
354 * @param m : principal quantum number of first atom
355 * @param n : principal quantum number of second atom
356 * @param R : internuclear distance in atomic units (bohr)
357 * @return the value of the kinetic energy integral
358 * @note N. Rosen, Phys. Rev., 38 (1931), 255
359 * @note untested
360 */
361inline RealType KinInt(RealType a, RealType b, int m, int n, RealType R) {
362 RealType KinInt_ = -0.5 * b * b * sSTOOvInt(a, b, m, n, R);
363 if (n > 0) {
364 KinInt_ +=
365 b * b * pow(2 * b / (2 * b - 1), 0.5) * sSTOOvInt(a, b, m, n - 1, R);
366 if (n > 1)
367 KinInt_ += pow(n * (n - 1) / ((n - 0.5) * (n - 1.5)), 0.5) *
368 sSTOOvInt(a, b, m, n - 2, R);
369 }
370 return KinInt_;
371}
372
373/**
374 * @brief Computes derivative of Coulomb integral with respect to the
375 * interatomic distance Computes the two-center Coulomb integral over
376 * Slater-type orbitals of s symmetry.
377 * @param a: Slater zeta exponent of first atom in inverse Bohr (au)
378 * @param b: Slater zeta exponent of second atom in inverse Bohr (au)
379 * @param m: principal quantum number of first atom
380 * @param n: principal quantum number of second atom
381 * @param R: internuclear distance in atomic units (bohr)
382 * @return the derivative of the Coulomb potential energy integral
383 * @note Derived in QTPIE research notes, May 15 2007
384 */
385inline RealType sSTOCoulIntGrad(RealType a, RealType b, int m, int n,
386 RealType R) {
387 RealType x;
388 // x is the argument of the auxiliary RosenA and RosenB functions
389 x = 2. * a * R;
390
391 // First compute the two-electron component
392 RealType sSTOCoulIntGrad_ = 0.;
393 if (x == 0) // Pathological case
394 {
395 printf("WARNING: argument given to sSTOCoulIntGrad is 0\n");
396 printf("a = %lf R= %lf\n", a, R);
397 } else {
398 RealType K2, TheSum;
399 if (a == b) {
400 TheSum = 0.;
401 for (int nu = 0; nu <= 2 * (n - 1); nu++) {
402 K2 = 0.;
403 for (int p = 0; p <= (m + nu) / 2; p++)
404 K2 += RosenD(2 * m - 1, nu + 1, 2 * p) / (2 * p + 1.) *
405 RosenA(2 * m + nu - 1 - 2 * p, x);
406 TheSum += (2 * n - nu - 1) / fact[nu] * pow(a * R, nu) * K2;
407 }
408 sSTOCoulIntGrad_ =
409 -pow(a, 2 * m + 2) * pow(R, 2 * m) / (n * fact[2 * m]) * TheSum;
410 TheSum = 0.;
411 for (int nu = 0; nu <= 2 * n - 1; nu++) {
412 K2 = 0.;
413 for (int p = 0; p <= (m + nu - 1) / 2; p++)
414 K2 += RosenD(2 * m - 1, nu, 2 * p) / (2 * p + 1.) *
415 RosenA(2 * m + nu - 2 * p, x);
416 TheSum += (2 * n - nu) / fact[nu] * pow(a * R, nu) * K2;
417 }
418 sSTOCoulIntGrad_ +=
419 2 * pow(a, 2 * m + 2) * pow(R, 2 * m) / (n * fact[2 * m]) * TheSum;
420 } else {
421 // Slater exponents are different
422 // First calculate some useful arguments
423 RealType y = R * (a + b);
424 RealType z = R * (a - b);
425 TheSum = 0.;
426 for (int nu = 0; nu <= 2 * n - 1; nu++) {
427 K2 = 0.;
428 for (int p = 0; p <= 2 * m + nu; p++)
429 K2 += RosenD(2 * m - 1, nu + 1, p) * RosenB(p, z) *
430 RosenA(2 * m + nu - p, y);
431 TheSum += (2 * n - nu - 1) / fact[nu] * pow(b * R, nu) * K2;
432 }
433 sSTOCoulIntGrad_ = -b * pow(a, 2 * m + 1) * pow(R, 2 * m) /
434 (2 * n * fact[2 * m]) * TheSum;
435 TheSum = 0.;
436 for (int nu = 0; nu <= 2 * n; nu++) {
437 K2 = 0.;
438 for (int p = 0; p <= 2 * m - 1 + nu; p++)
439 K2 += RosenD(2 * m - 1, nu, p) *
440 ((a - b) * RosenB(p + 1, z) * RosenA(2 * m + nu - p - 1, y) +
441 (a + b) * RosenB(p, z) * RosenA(2 * m + nu - p, y));
442 TheSum += (2 * n - nu) / fact[nu] * pow(b * R, nu) * K2;
443 }
444 sSTOCoulIntGrad_ +=
445 pow(a, 2 * m + 1) * pow(R, 2 * m) / (2 * n * fact[2 * m]) * TheSum;
446 }
447 // Now add one-electron terms and common term
448 sSTOCoulIntGrad_ = sSTOCoulIntGrad_ - (2. * m + 1.) / sqr(R) +
449 2. * m / R * sSTOCoulInt(a, b, m, n, R) +
450 pow(x, 2 * m) / (fact[2 * m] * sqr(R)) *
451 ((2. * m + 1.) * exp(-x) +
452 2. * m * (1. + 2. * m - x) * RosenA(2 * m - 1, x));
453 }
454 return sSTOCoulIntGrad_;
455}
456
457/**
458 * @brief Computes gradient of overlap integral with respect to the interatomic
459 * diatance Computes the derivative of the overlap integral over two Slater-type
460 * orbitals of s symmetry.
461 * @param a: Slater zeta exponent of first atom in inverse Bohr (au)
462 * @param b: Slater zeta exponent of second atom in inverse Bohr (au)
463 * @param m: principal quantum number of first atom
464 * @param n: principal quantum number of second atom
465 * @param R: internuclear distance in atomic units (bohr)
466 * @return the derivative of the sSTOOvInt integral
467 * @note Derived in QTPIE research notes, May 15 2007
468 */
469inline RealType sSTOOvIntGrad(RealType a, RealType b, int m, int n,
470 RealType R) {
471 // Calculate first term
472 RealType sSTOOvIntGrad_ = (m + n + 1.) / R * sSTOOvInt(a, b, m, n, R);
473
474 // Calculate remaining terms; answers depend on exponents
475 RealType TheSum = 0.;
476 RealType x = a * R;
477 if (a == b) {
478 for (int q = 0; q <= (m + n) / 2; q++)
479 TheSum +=
480 RosenD(m, n, 2 * q) / (2 * q + 1.) * RosenA(m + n - 2 * q + 1, x);
481 sSTOOvIntGrad_ -=
482 a * pow(x, m + n + 1) / sqrt(fact[2 * m] * fact[2 * n]) * TheSum;
483 } else {
484 RealType w = b * R;
485 RealType y = 0.5 * R * (a + b);
486 RealType z = 0.5 * R * (a - b);
487 for (int q = 0; q < m + n; q++)
488 TheSum = TheSum + RosenD(m, n, q) *
489 ((a - b) * RosenB(q + 1, z) * RosenA(m + n - q, y) +
490 (a + b) * RosenB(q, z) * RosenA(m + n - q + 1, y));
491 sSTOOvIntGrad_ -= 0.25 *
492 sqrt((pow(x, 2 * m + 1) * pow(w, 2 * n + 1)) /
493 (fact[2 * m] * fact[2 * n])) *
494 TheSum;
495 }
496 return sSTOOvIntGrad_;
497}
498
499/**
500 * @brief Calculates a Slater-type orbital exponent based on the hardness
501 * parameters
502 * @param hardness: chemical hardness in atomic units
503 * @param n: principal quantum number
504 * @note Modified for use with OpenMD by Gezelter and Michalka.
505 */
506inline RealType getSTOZeta(int n, RealType hardness) {
507 // Approximate the exact value of the constant of proportionality
508 // by its value at a very small distance epsilon
509 // since the exact R = 0 case has not be programmed
510 RealType epsilon = 1.0e-8;
511
512 // Assign orbital exponent
513 return pow(sSTOCoulInt(1., 1., n, n, epsilon) / hardness,
514 -1. / (3. + 2. * n));
515}
516
517/** @brief Computes Coulomb integral analytically using s-type Slater-style
518 * densities under the shifted force approximation Computes the two-center
519 * Coulomb integral over Slater-type densities of s symmetry using a shifted
520 * force approximation
521 * @param a : Slater zeta exponent of first atom in inverse Angstroms
522 * @param b : Slater zeta exponent of second atom in inverse Angstroms
523 * @param R : internuclear distance in Angstroms
524 * @param Rc: internuclear cutoff distance in Angstroms
525 * @return value of the shifted-force Coulomb potential energy integral
526 */
527inline RealType sSTDCoulSF(RealType a, RealType b, RealType R, RealType Rc) {
528 RealType a2 = a * a;
529 RealType a4 = a2 * a2;
530 RealType b2 = b * b;
531 RealType b4 = b2 * b2;
532 RealType apb = a + b;
533 RealType a2b23 = pow(a2 - b2, 3);
534
535 RealType term1 = 2.0 * exp(R * apb) * a2b23;
536 RealType term2 = exp(R * b) * b4 * (b2 * (2.0 + R * a) - a2 * (6.0 + R * a));
537 RealType term3 = -exp(R * a) * a4 * (a2 * (2.0 + R * b) - b2 * (6.0 + R * b));
538
539 RealType denom = 2.0 * R * a2b23;
540 RealType numer = exp(-R * apb) * (term1 + term2 + term3);
541
542 RealType S = numer / denom;
543
544 term1 = 2.0 * exp(Rc * apb) * a2b23;
545 term2 = exp(Rc * b) * b4 * (b2 * (2.0 + Rc * a) - a2 * (6.0 + Rc * a));
546 term3 = -exp(Rc * a) * a4 * (a2 * (2.0 + Rc * b) - b2 * (6.0 + Rc * b));
547
548 numer = exp(-Rc * apb) * (term1 + term2 + term3);
549 denom = 2.0 * Rc * pow(a2 - b2, 3);
550
551 RealType Sc = numer / denom;
552
553 term1 = -2.0 * exp(3 * Rc * a + 4 * Rc * b) * a2b23;
554 term2 = exp(2.0 * Rc * (a + 2.0 * b)) * b4 *
555 (a2 * (6.0 + Rc * a * (6.0 + Rc * a)) -
556 b2 * (2.0 + Rc * a * (2.0 + Rc * a)));
557 term3 = exp(3.0 * Rc * apb) * a4 *
558 (a2 * (2.0 + Rc * b * (2.0 + Rc * b)) -
559 b2 * (6.0 + Rc * b * (6.0 + Rc * b)));
560
561 numer = exp(-3.0 * Rc * a - 4.0 * Rc * b) * (term1 + term2 + term3);
562 denom = 2.0 * Rc * Rc * a2b23;
563
564 RealType Scp = numer / denom;
565
566 return S - Sc - Scp * (R - Rc);
567}
568
569/** @brief Computes Coulomb integral analytically using s-type Slater-style
570 * densities under the shifted force approximation Computes the two-center
571 * Coulomb integral over Slater-type densities of s symmetry using a shifted
572 * force approximation
573 * @param a : Slater zeta exponent of both atoms in inverse Angstroms
574 * @param R : internuclear distance in Angstroms
575 * @param Rc: internuclear cutoff distance in Angstroms
576 * @return value of the shifted-force Coulomb potential energy integral
577 */
578inline RealType sSTDCoulSF(RealType a, RealType R, RealType Rc) {
579 RealType poly =
580 (-48.0 + 48.0 * exp(R * a) - R * a * (33.0 + R * a * (9.0 + R * a)));
581 RealType polyc =
582 (-48.0 + 48.0 * exp(Rc * a) - Rc * a * (33.0 + Rc * a * (9.0 + Rc * a)));
583
584 RealType S = exp(-R * a) * poly / (48.0 * R);
585 RealType Sc = exp(-Rc * a) * polyc / (48.0 * Rc);
586
587 RealType poly2c =
588 (48.0 - 48.0 * exp(Rc * a) +
589 Rc * a * (4.0 + Rc * a) * (12.0 + Rc * a * (3.0 + Rc * a)));
590 RealType Scp = exp(-Rc * a) * poly2c / (48.0 * Rc * Rc);
591
592 return S - Sc - Scp * (R - Rc);
593}
594
595#endif