3
* These routines are based on the gamfun program of
4
* Trygve Ulf Helgaker (fall 1984)
5
* and calculates the incomplete gamma function as
6
* described by McMurchie & Davidson, J. Comp. Phys. 26 (1978) 218.
7
* The original routine computed the function for maximum j = 20.
13
#include <libciomr/libciomr.h>
14
#include <libint/libint.h>
25
namespace psi { namespace CINTS {
28
static double_matrix_t gtable;
29
static double *denomarray;
30
static double wval_infinity;
31
static int itable_infinity;
33
/* Tablesize should always be at least 121. */
36
/*! Tabulate the incomplete gamma function and put in gtable. */
38
* For J = JMAX a power series expansion is used, see for
39
* example Eq.(39) given by V. Saunders in "Computational
40
* Techniques in Quantum Chemistry and Molecular Physics",
41
* Reidel 1975. For J < JMAX the values are calculated
42
* using downward recursion in J.
44
void init_fjt(int max)
47
double denom,d2jmax1,r2jmax1,wval,d2wval,sum,term,rexpw;
51
/* Check if int_fjttable has been allocated with
52
the same maxj before */
53
if (gtable.d != NULL) {
54
if (maxj+7 <= gtable.n1)
61
/* Allocate storage for gtable */
62
if (gtable.d == NULL) {
64
gtable.n2 = TABLESIZE;
65
gtable.d = block_matrix(gtable.n1,gtable.n2);
68
/* Tabulate the gamma function for t(=wval)=0.0. */
70
for (i=0; i<gtable.n1; i++) {
71
gtable.d[i][0] = 1.0/denom;
75
/* Tabulate the gamma function from t(=wval)=0.1, to 12.0. */
76
d2jmax1 = 2.0*(gtable.n1-1) + 1.0;
77
r2jmax1 = 1.0/d2jmax1;
78
for (i=1; i<TABLESIZE; i++) {
84
for (j=2; j<=200; j++) {
86
term = term * d2wval / denom;
88
if (term <= 1.0e-15) break;
92
/* Fill in the values for the highest j gtable entries (top row). */
93
gtable.d[gtable.n1-1][i] = rexpw * sum;
95
/* Work down the table filling in the rest of the column. */
97
for (j=gtable.n1 - 2; j>=0; j--) {
99
gtable.d[j][i] = (gtable.d[j+1][i]*d2wval + rexpw)/denom;
103
/* Form some denominators, so divisions can be eliminated below. */
104
if (denomarray == NULL) {
105
denomarray = init_array(max+1);
108
for (i=1; i<=max; i++) {
109
denomarray[i] = 1.0/(2*i - 1);
112
wval_infinity = 2*max + 37.0;
113
itable_infinity = (int) (10 * wval_infinity);
119
/*! This is called when the fjt routines are no longer needed, or
120
* before they are reinitialized with a new maxj. */
124
free_block(gtable.d);
134
This is clean and ugly at the same time - have to call this AFTER
135
init_fjt() to have proper maxj
137
void init_fjt_table(double_array_t *table)
139
if (gtable.d != NULL) {
141
table->d = init_array(maxj);
144
throw std::domain_error("Called init_fjt_table before init_fjt");
149
void free_fjt_table(double_array_t *table)
151
if (table->d != NULL)
157
/*! Using the tabulated incomplete gamma function in gtable, compute
158
* the incomplete gamma function for a particular wval for all 0<=j<=J.
159
* The result is placed in the global intermediate int_fjttable.
161
void int_fjt(double_array_t *table, int J, double wval)
163
const double sqrpih = 0.886226925452758;
164
const double coef2 = 0.5000000000000000;
165
const double coef3 = -0.1666666666666667;
166
const double coef4 = 0.0416666666666667;
167
const double coef5 = -0.0083333333333333;
168
const double coef6 = 0.0013888888888889;
169
const double gfac30 = 0.4999489092;
170
const double gfac31 = -0.2473631686;
171
const double gfac32 = 0.321180909;
172
const double gfac33 = -0.3811559346;
173
const double gfac20 = 0.4998436875;
174
const double gfac21 = -0.24249438;
175
const double gfac22 = 0.24642845;
176
const double gfac10 = 0.499093162;
177
const double gfac11 = -0.2152832;
178
const double gfac00 = -0.490;
180
double wdif, d2wal, rexpw, /* denom, */ gval, factor, rwval, term;
181
int i, itable, irange;
184
fprintf(stderr,"the int_fjt routine has been incorrectly used\n");
185
fprintf(stderr,"J = %d but maxj = %d\n",J,maxj);
189
/* Compute an index into the table. */
190
/* The test is needed to avoid floating point exceptions for
191
* large values of wval. */
192
if (wval > wval_infinity) {
193
itable = itable_infinity;
196
itable = (int) (10.0 * wval);
199
/* If itable is small enough use the table to compute int_fjttable. */
200
if (itable < TABLESIZE) {
202
wdif = wval - 0.1 * itable;
204
/* Compute fjt for J. */
205
table->d[J] = (((((coef6 * gtable.d[J+6][itable]*wdif
206
+ coef5 * gtable.d[J+5][itable])*wdif
207
+ coef4 * gtable.d[J+4][itable])*wdif
208
+ coef3 * gtable.d[J+3][itable])*wdif
209
+ coef2 * gtable.d[J+2][itable])*wdif
210
- gtable.d[J+1][itable])*wdif
211
+ gtable.d[J][itable];
213
/* Compute the rest of the fjt. */
216
/* denom = 2*J + 1; */
217
for (i=J; i>0; i--) {
218
/* denom = denom - 2.0; */
219
table->d[i-1] = (d2wal*table->d[i] + rexpw)*denomarray[i];
222
/* If wval <= 2*J + 36.0, use the following formula. */
223
else if (itable <= 20*J + 360) {
227
/* Subdivide wval into 6 ranges. */
228
irange = itable/30 - 3;
230
gval = gfac30 + rwval*(gfac31 + rwval*(gfac32 + rwval*gfac33));
231
table->d[0] = sqrpih*sqrt(rwval) - rexpw*gval*rwval;
233
else if (irange == 2) {
234
gval = gfac20 + rwval*(gfac21 + rwval*gfac22);
235
table->d[0] = sqrpih*sqrt(rwval) - rexpw*gval*rwval;
237
else if (irange == 3 || irange == 4) {
238
gval = gfac10 + rwval*gfac11;
239
table->d[0] = sqrpih*sqrt(rwval) - rexpw*gval*rwval;
241
else if (irange == 5 || irange == 6) {
243
table->d[0] = sqrpih*sqrt(rwval) - rexpw*gval*rwval;
246
table->d[0] = sqrpih*sqrt(rwval);
249
/* Compute the rest of the int_fjttable from table->d[0]. */
250
factor = 0.5 * rwval;
251
term = factor * rexpw;
252
for (i=1; i<=J; i++) {
253
table->d[i] = factor * table->d[i-1] - term;
254
factor = rwval + factor;
257
/* For large values of wval use this algorithm: */
260
table->d[0] = sqrpih*sqrt(rwval);
261
factor = 0.5 * rwval;
262
for (i=1; i<=J; i++) {
263
table->d[i] = factor * table->d[i-1];
264
factor = rwval + factor;
267
/* printf(" %2d %12.8f %4d %12.8f\n",J,wval,itable,table->d[0]); */
added
removed
src/bin/cints/Tools/int_fjt.cc
