ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/branches/development/src/nonbonded/SlaterIntegrals.hpp
Revision: 1874
Committed: Wed May 15 15:09:35 2013 UTC (12 years, 2 months ago) by gezelter
File size: 19186 byte(s)
Log Message:
Fixed a bunch of cppcheck warnings.

File Contents

# User Rev Content
1 gezelter 1718 /***************************************************************************
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     #include <cmath>
64 gezelter 1749 #include <cstdlib>
65 jmichalk 1734 #include <iostream>
66 gezelter 1718 #include "math/Factorials.hpp"
67 gezelter 1766 #include "utils/NumericConstant.hpp"
68 gezelter 1718
69     #ifndef NONBONDED_SLATERINTEGRALS_HPP
70     #define NONBONDED_SLATERINTEGRALS_HPP
71    
72     template <typename T> inline T sqr(T t) { return t*t; }
73     template <typename T> inline T mod(T x, T m)
74     { return x<0 ? m - 1 - ((-x) - 1)%m : x%m; }
75    
76     // #include "parameters.h"
77    
78     /**
79     * @brief Computes Rosen's Guillimer-Zener function A
80     * Computes Rosen's A integral, an auxiliary quantity needed to
81     * compute integrals involving Slater-type orbitals of s symmetry.
82     * \f[
83     * A_n(\alpha) = \int_1^\infty x^n e^{-\alpha x}dx
84     * = \frac{n! e^{-\alpha}}{\alpha^{n+1}}\sum_{\nu=0}^n
85     * \frac{\alpha^\nu}{\nu!}
86     * \f]
87     * @param n - principal quantum number
88 gezelter 1808 * @param a - Slater exponent
89 gezelter 1718 * @return the value of Rosen's A integral
90     * @note N. Rosen, Phys. Rev., 38 (1931), 255
91     */
92     inline RealType RosenA(int n, RealType a)
93     {
94     RealType RosenA_ = 0.;
95     if (a != 0.)
96     {
97     RealType Term = 1.;
98     RosenA_ = Term;
99 jmichalk 1734 for (int nu=1; nu<=n; nu++)
100 gezelter 1718 {
101     Term *= a / nu;
102     RosenA_ += Term;
103     }
104     RosenA_ = (RosenA_/Term) * (exp(-a)/a);
105     }
106     return RosenA_;
107     }
108    
109     /**
110     * @brief Computes Rosen's Guillimer-Zener function B
111     * Computes Rosen's B integral, an auxiliary quantity needed to
112     * compute integrals involving Slater-type orbitals of s symmetry.
113     * \f[
114     * B_n(\alpha) = \int_{-1}^1 x^n e^{-\alpha x} dx
115     * = \frac{n!}{\alpha^{n+1}}
116     * \sum_{\nu=0}^n \frac{e^\alpha(-\alpha)^\nu
117     * - e^{-\alpha} \alpha^\nu}{\nu!}
118     * \f]
119     * @param n - principal quantum number
120     * @param alpha - Slater exponent
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;
127     RealType RosenB_, PSinhRosenA, PCoshRosenA, PHyperRosenA;
128 gezelter 1874
129 gezelter 1718 if (alpha != 0.)
130     {
131     Term = 1.;
132 gezelter 1874 bool IsPositive = true;
133 gezelter 1718
134     // These two expressions are (up to constant factors) equivalent
135     // to computing the hyperbolic sine and cosine of a respectively
136     // The series consists of adding up these terms in an alternating fashion
137     PSinhRosenA = exp(alpha) - exp(-alpha);
138     PCoshRosenA = -exp(alpha) - exp(-alpha);
139     TheSum = PSinhRosenA;
140     for (unsigned nu=1; nu<=n; nu++)
141     {
142     if (IsPositive)
143     {
144     PHyperRosenA = PCoshRosenA;
145     IsPositive = false;
146     }
147     else // term to add should be negative
148     {
149     PHyperRosenA = PSinhRosenA;
150     IsPositive = true;
151     }
152     Term *= alpha / nu;
153     TheSum += Term * PHyperRosenA;
154     }
155     RosenB_ = TheSum / (alpha*Term);
156     }
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     */
175     inline RealType RosenD(int m, int n, int p)
176     {
177     if (m+n+p > maxFact)
178     {
179     printf("Error, arguments exceed maximum factorial computed %d > %d\n", m+n+p, maxFact);
180     ::exit(0);
181     }
182    
183     RealType RosenD_ = 0;
184     for (int k=max(p-m,0); k<=min(n,p); k++)
185     {
186     if (mod(k,2) == 0)
187     RosenD_ += (fact[m] / (fact[p-k] * fact[m-p+k])) * (fact[n] / (fact[n-k] * fact[k]));
188     else
189     RosenD_ -= (fact[m] / ( fact[p-k] * fact[m-p+k])) * (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 symmetry.
196     * @param a : Slater zeta exponent of first atom in inverse Bohr (au)
197     * @param b : Slater zeta exponent of second atom in inverse Bohr (au)
198     * @param m : principal quantum number of first atom
199     * @param n : principal quantum number of second atom
200     * @param R : internuclear distance in atomic units (bohr)
201     * @return value of the Coulomb potential energy integral
202     * @note N. Rosen, Phys. Rev., 38 (1931), 255
203     * @note In Rosen's paper, this integral is known as K2.
204     */
205     inline RealType sSTOCoulInt(RealType a, RealType b, int m, int n, RealType R)
206     {
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 gezelter 1766 if (std::fabs(x) < OpenMD::NumericConstant::epsilon) // Pathological case
221 gezelter 1718 {
222 gezelter 1766
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] / (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 gezelter 1718 }
255     else
256     {
257     OneElectronTerm = 1./R + pow(x, 2*m)/(fact[2*m]*R)*
258     ((x-2*m)*RosenA(2*m-1,x)-exp(-x)) + sSTOCoulInt_;
259     eps = epsilon / OneElectronTerm;
260     if (a == b)
261     {
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     {
266     Factor2 = (2.*n-nu)/fact[nu]*pow(a*R,nu);
267     epsi = eps / fabs(Factor1 * Factor2);
268     K2 = 0.;
269     for (int p=0; p<=m+(nu-1)/2; p++)
270     {
271     Term = RosenD(2*m-1, nu, 2*p)/(2.*p+1.) *RosenA(2*m+nu-1-2*p,x);
272     K2 += Term;
273     if ((Term > 0) && (Term < epsi)) goto label1;
274     }
275     sSTOCoulInt_ += K2 * Factor2;
276     }
277     label1:
278     sSTOCoulInt_ *= Factor1;
279     }
280     else
281     {
282     Factor1 = -a*pow(a*R,2*m)/(2.*n*fact[2*m]);
283     epsi = eps/fabs(Factor1);
284     if (b == 0.)
285     printf("WARNING: b = 0 in sSTOCoulInt\n");
286     else
287     {
288     // Apply Rosen (54)
289     for (int nu=0; nu<=2*n-1; nu++)
290     {
291     K2 = 0;
292     for (int p=0; p<=2*m+nu-1; p++)
293     K2=K2+RosenD(2*m-1,nu,p)*RosenB(p,R*(a-b))
294     *RosenA(2*m+nu-1-p,R*(a+b));
295     Term = K2*(2*n-nu)/fact[nu]*pow(b*R, nu);
296     sSTOCoulInt_ += Term;
297     if (fabs(Term) < epsi) goto label2;
298     }
299     label2:
300     sSTOCoulInt_ *= Factor1;
301     }
302     }
303     // Now add the one-electron term from Rosen (47) = Rosen (53)
304     sSTOCoulInt_ += OneElectronTerm;
305     }
306     return sSTOCoulInt_;
307     }
308    
309     /**
310     * @brief Computes overlap integral analytically over s-type STOs
311     * Computes the overlap integral over two
312     * Slater-type orbitals of s symmetry.
313     * @param a : Slater zeta exponent of first atom in inverse Bohr (au)
314     * @param b : Slater zeta exponent of second atom in inverse Bohr (au)
315     * @param m : principal quantum number of first atom
316     * @param n : principal quantum number of second atom
317     * @param R : internuclear distance in atomic units (bohr)
318     * @return the value of the sSTOOvInt integral
319     * @note N. Rosen, Phys. Rev., 38 (1931), 255
320     * @note In the Rosen paper, this integral is known as I.
321     */
322     inline RealType sSTOOvInt(RealType a, RealType b, int m, int n, RealType R)
323     {
324     RealType Factor, Term, eps;
325    
326     // To speed up calculation, we terminate loop once contributions
327     // to integral fall below the bound, epsilon
328     RealType epsilon = 0.;
329     RealType sSTOOvInt_ = 0.;
330    
331     if (a == b)
332     {
333     Factor = pow(a*R, m+n+1)/sqrt(fact[2*m]*fact[2*n]);
334     eps = epsilon / fabs(Factor);
335     for (int q=0; q<=(m+n)/2; q++)
336     {
337     Term = RosenD(m,n,2*q)/(2.*q+1.)*RosenA(m+n-2*q,a*R);
338     sSTOOvInt_ += Term;
339     if (fabs(Term) < eps) exit(0);
340     }
341     sSTOOvInt_ *= Factor;
342     }
343     else
344     {
345     Factor = 0.5*pow(a*R, m+0.5)*pow(b*R,n+0.5)
346     /sqrt(fact[2*m]*fact[2*n]);
347     eps = epsilon / fabs(Factor);
348     for (int q=0; q<=m+n; q++)
349     {
350     Term = RosenD(m,n,q)*RosenB(q, R/2.*(a-b))
351     * RosenA(m+n-q,R/2.*(a+b));
352     sSTOOvInt_ += Term;
353     if (fabs(Term) < eps) exit(0);
354     }
355     sSTOOvInt_ *= Factor;
356     }
357     return sSTOOvInt_;
358     }
359    
360     /**
361     * @brief Computes kinetic energy integral analytically over s-type STOs
362     * Computes the overlap integral over two Slater-type orbitals of s symmetry.
363     * @param a : Slater zeta exponent of first atom in inverse Bohr (au)
364     * @param b : Slater zeta exponent of second atom in inverse Bohr (au)
365     * @param m : principal quantum number of first atom
366     * @param n : principal quantum number of second atom
367     * @param R : internuclear distance in atomic units (bohr)
368     * @return the value of the kinetic energy integral
369     * @note N. Rosen, Phys. Rev., 38 (1931), 255
370     * @note untested
371     */
372     inline RealType KinInt(RealType a, RealType b, int m, int n,RealType R)
373     {
374     RealType KinInt_ = -0.5*b*b*sSTOOvInt(a, b, m, n, R);
375     if (n > 0)
376     {
377     KinInt_ += b*b*pow(2*b/(2*b-1),0.5) * sSTOOvInt(a, b, m, n-1, R);
378     if (n > 1) KinInt_ += pow(n*(n-1)/((n-0.5)*(n-1.5)), 0.5)
379     * sSTOOvInt(a, b, m, n-2, R);
380     }
381     return KinInt_;
382     }
383    
384     /**
385     * @brief Computes derivative of Coulomb integral with respect to the interatomic distance
386     * Computes the two-center Coulomb integral over Slater-type orbitals of s symmetry.
387     * @param a: Slater zeta exponent of first atom in inverse Bohr (au)
388     * @param b: Slater zeta exponent of second atom in inverse Bohr (au)
389     * @param m: principal quantum number of first atom
390     * @param n: principal quantum number of second atom
391     * @param R: internuclear distance in atomic units (bohr)
392     * @return the derivative of the Coulomb potential energy integral
393     * @note Derived in QTPIE research notes, May 15 2007
394     */
395     inline RealType sSTOCoulIntGrad(RealType a, RealType b, int m, int n, RealType R)
396     {
397     RealType x, y, z, K2, TheSum;
398     // x is the argument of the auxiliary RosenA and RosenB functions
399     x = 2. * a * R;
400    
401     // First compute the two-electron component
402     RealType sSTOCoulIntGrad_ = 0.;
403     if (x==0) // Pathological case
404     {
405     printf("WARNING: argument given to sSTOCoulIntGrad is 0\n");
406     printf("a = %lf R= %lf\n", a, R);
407     }
408     else
409     {
410     if (a == b)
411     {
412     TheSum = 0.;
413     for (int nu=0; nu<=2*(n-1); nu++)
414     {
415     K2 = 0.;
416     for (int p=0; p<=(m+nu)/2; p++)
417     K2 += RosenD(2*m-1, nu+1, 2*p)/(2*p + 1.) * RosenA(2*m+nu-1-2*p, x);
418     TheSum += (2*n-nu-1)/fact[nu]*pow(a*R, nu) * K2;
419     }
420     sSTOCoulIntGrad_ = -pow(a, 2*m+2)*pow(R, 2*m) /(n*fact[2*m])*TheSum;
421     TheSum = 0.;
422     for (int nu=0; nu<=2*n-1; nu++)
423     {
424     K2 = 0.;
425     for (int p=0; p<=(m+nu-1)/2; p++)
426     K2 += RosenD(2*m-1, nu, 2*p)/(2*p + 1.) * RosenA(2*m+nu-2*p, x);
427     TheSum += (2*n-nu)/fact[nu]*pow(a*R,nu) * K2;
428     }
429     sSTOCoulIntGrad_ += 2*pow(a, 2*m+2)*pow(R, 2*m) /(n*fact[2*m])*TheSum;
430     }
431     else
432     {
433     // Slater exponents are different
434     // First calculate some useful arguments
435     y = R*(a+b);
436     z = R*(a-b);
437     TheSum = 0.;
438     for (int nu=0; nu<=2*n-1; nu++)
439     {
440     K2 = 0.;
441     for (int p=0; p<=2*m+nu; p++)
442     K2 += RosenD(2*m-1, nu+1, p)
443     * RosenB(p,z)*RosenA(2*m+nu-p, y);
444     TheSum += (2*n-nu-1)/fact[nu]*pow(b*R,nu) * K2;
445     }
446     sSTOCoulIntGrad_ = -b*pow(a,2*m+1)*pow(R,2*m)/
447     (2*n*fact[2*m])*TheSum;
448     TheSum = 0.;
449     for (int nu=0; nu<=2*n; nu++)
450     {
451     K2 = 0.;
452     for (int p=0; p<=2*m-1+nu; p++)
453     K2 += RosenD(2*m-1, nu, p)
454     * ((a-b)*RosenB(p+1,z)*RosenA(2*m+nu-p-1, y)
455     +(a+b)*RosenB(p ,z)*RosenA(2*m+nu-p , y));
456     TheSum += (2*n-nu)/fact[nu]*pow(b*R,nu) * K2;
457     }
458     sSTOCoulIntGrad_ += pow(a,2*m+1)*pow(R,2*m)/(2*n*fact[2*m])*TheSum;
459     }
460     // Now add one-electron terms and common term
461     sSTOCoulIntGrad_ = sSTOCoulIntGrad_ - (2.*m+1.)/sqr(R)
462     + 2.*m/R * sSTOCoulInt(a,b,m,n,R)
463     + pow(x,2*m)/(fact[2*m]*sqr(R)) * ((2.*m+1.)*exp(-x)
464     + 2.*m*(1.+2.*m-x)*RosenA(2*m-1,x));
465     }
466     return sSTOCoulIntGrad_;
467     }
468    
469     /**
470     * @brief Computes gradient of overlap integral with respect to the interatomic diatance
471     * Computes the derivative of the overlap integral over two Slater-type orbitals of s symmetry.
472     * @param a: Slater zeta exponent of first atom in inverse Bohr (au)
473     * @param b: Slater zeta exponent of second atom in inverse Bohr (au)
474     * @param m: principal quantum number of first atom
475     * @param n: principal quantum number of second atom
476     * @param R: internuclear distance in atomic units (bohr)
477     * @return the derivative of the sSTOOvInt integral
478     * @note Derived in QTPIE research notes, May 15 2007
479     */
480     inline RealType sSTOOvIntGrad(RealType a, RealType b, int m, int n, RealType R)
481     {
482     RealType w, x, y, z, TheSum;
483    
484     // Calculate first term
485     RealType sSTOOvIntGrad_ = (m+n+1.)/R * sSTOOvInt(a, b, m, n, R);
486    
487     // Calculate remaining terms; answers depend on exponents
488     TheSum = 0.;
489     x = a * R;
490     if (a == b)
491     {
492     for (int q=0; q<=(m+n)/2; q++)
493     TheSum += RosenD(m,n,2*q) / (2*q + 1.) * RosenA(m+n-2*q+1, x);
494     sSTOOvIntGrad_ -= a*pow(x,m+n+1)/ sqrt(fact[2*m]*fact[2*n])*TheSum;
495     }
496     else
497     {
498     w = b*R;
499     y = 0.5*R*(a+b);
500     z = 0.5*R*(a-b);
501     for (int q=0; q<m+n; q++)
502     TheSum = TheSum + RosenD(m,n,q) *
503     ((a-b)*RosenB(q+1,z)*RosenA(m+n-q ,y)
504     +(a+b)*RosenB(q ,z)*RosenA(m+n-q+1,y));
505     sSTOOvIntGrad_ -= 0.25*sqrt((pow(x, 2*m+1)*pow(w, 2*n+1))/(fact[2*m]*fact[2*n]))*TheSum;
506     }
507     return sSTOOvIntGrad_;
508     }
509    
510     /**
511     * @brief Calculates a Slater-type orbital exponent based on the hardness parameters
512 gezelter 1808 * @param hardness: chemical hardness in atomic units
513 gezelter 1718 * @param n: principal quantum number
514     * @note Modified for use with OpenMD by Gezelter and Michalka.
515     */
516     inline RealType getSTOZeta(int n, RealType hardness)
517     {
518     // Approximate the exact value of the constant of proportionality
519     // by its value at a very small distance epsilon
520     // since the exact R = 0 case has not be programmed
521     RealType epsilon = 1.0e-8;
522    
523     // Assign orbital exponent
524     return pow(sSTOCoulInt(1., 1., n, n, epsilon) / hardness, -1./(3. + 2.*n));
525     }
526    
527     #endif

Properties

Name Value
svn:eol-style native