1 //
2 // fjt.cc
3 //
4 // Copyright (C) 1996 Limit Point Systems, Inc.
5 //
6 // Author: Curtis Janssen <cljanss@limitpt.com>
7 // Maintainer: LPS
8 //
9 // This file is part of the SC Toolkit.
10 //
11 // The SC Toolkit is free software; you can redistribute it and/or modify
12 // it under the terms of the GNU Library General Public License as published by
13 // the Free Software Foundation; either version 2, or (at your option)
14 // any later version.
15 //
16 // The SC Toolkit is distributed in the hope that it will be useful,
17 // but WITHOUT ANY WARRANTY; without even the implied warranty of
18 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
19 // GNU Library General Public License for more details.
20 //
21 // You should have received a copy of the GNU Library General Public License
22 // along with the SC Toolkit; see the file COPYING.LIB. If not, write to
23 // the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
24 //
25 // The U.S. Government is granted a limited license as per AL 91-7.
26 //
27
28 #ifdef __GNUG__
29 #pragma implementation
30 #endif
31
32 /* These routines are based on the gamfun program of
33 * Trygve Ulf Helgaker (fall 1984)
34 * and calculates the incomplete gamma function as
35 * described by McMurchie & Davidson, J. Comp. Phys. 26 (1978) 218.
36 * The original routine computed the function for maximum j = 20.
37 */
38
39 #include <stdlib.h>
40 #include <math.h>
41
42 #include <iostream>
43
44 #include <util/misc/formio.h>
45 #include <util/misc/exenv.h>
46 #include <chemistry/qc/intv3/fjt.h>
47
48 using namespace std;
49
50 /* Tablesize should always be at least 121. */
51 #define TABLESIZE 121
52
53 /* Tabulate the incomplete gamma function and put in gtable. */
54 /*
55 * For J = JMAX a power series expansion is used, see for
56 * example Eq.(39) given by V. Saunders in "Computational
57 * Techniques in Quantum Chemistry and Molecular Physics",
58 * Reidel 1975. For J < JMAX the values are calculated
59 * using downward recursion in J.
60 */
61 FJT::FJT(int max)
62 {
63 int i,j;
64 double denom,d2jmax1,r2jmax1,wval,d2wval,sum,term,rexpw;
65
66 maxj = max;
67
68 /* Allocate storage for gtable and int_fjttable. */
69 int_fjttable = new double[maxj+1];
70 gtable = new double*[ngtable()];
71 for (i=0; i<ngtable(); i++) {
72 gtable[i] = new double[TABLESIZE];
73 }
74
75 /* Tabulate the gamma function for t(=wval)=0.0. */
76 denom = 1.0;
77 for (i=0; i<ngtable(); i++) {
78 0 0 0 1 0 0 0 0 0 0 gtable[i][0] = 1.0/denom;
79 denom += 2.0;
80 }
81
82 /* Tabulate the gamma function from t(=wval)=0.1, to 12.0. */
83 d2jmax1 = 2.0*(ngtable()-1) + 1.0;
84 r2jmax1 = 1.0/d2jmax1;
85 for (i=1; i<TABLESIZE; i++) {
86 wval = 0.1 * i;
87 d2wval = 2.0 * wval;
88 term = r2jmax1;
89 sum = term;
90 denom = d2jmax1;
91 for (j=2; j<=200; j++) {
92 0 0 0 1 0 0 0 0 0 0 denom = denom + 2.0;
93 31 3 7 1 0 2 0 0 0 0 term = term * d2wval / denom;
94 0 0 1 0 0 0 0 0 0 0 sum = sum + term;
95 7 0 5 0 0 0 0 0 0 0 if (term <= 1.0e-15) break;
96 1 2 4 0 0 0 0 0 0 0 }
97 1 0 0 0 0 0 0 0 0 0 rexpw = exp(-wval);
98
99 /* Fill in the values for the highest j gtable entries (top row). */
100 gtable[ngtable()-1][i] = rexpw * sum;
101
102 /* Work down the table filling in the rest of the column. */
103 denom = d2jmax1;
104 for (j=ngtable() - 2; j>=0; j--) {
105 denom = denom - 2.0;
106 45 2 12 30 5 1 0 0 0 0 gtable[j][i] = (gtable[j+1][i]*d2wval + rexpw)/denom;
107 }
108 0 0 0 1 0 0 0 0 0 0 }
109
110 /* Form some denominators, so divisions can be eliminated below. */
111 denomarray = new double[max+1];
112 denomarray[0] = 0.0;
113 for (i=1; i<=max; i++) {
114 denomarray[i] = 1.0/(2*i - 1);
115 1 0 0 0 0 0 0 0 0 0 }
116
117 wval_infinity = 2*max + 37.0;
118 itable_infinity = (int) (10 * wval_infinity);
119
120 }
121
122 FJT::~FJT()
123 {
124 delete[] int_fjttable;
125 for (int i=0; i<maxj+7; i++) {
126 delete[] gtable[i];
127 }
128 delete[] gtable;
129 delete[] denomarray;
130 }
131
132 /* Using the tabulated incomplete gamma function in gtable, compute
133 * the incomplete gamma function for a particular wval for all 0<=j<=J.
134 * The result is placed in the global intermediate int_fjttable.
135 */
136 double *
137 FJT::values(int J,double wval)
138 276 32 77 168 16 1 0 2 0 0 {
139 const double sqrpih = 0.886226925452758;
140 const double coef2 = 0.5000000000000000;
141 const double coef3 = -0.1666666666666667;
142 const double coef4 = 0.0416666666666667;
143 const double coef5 = -0.0083333333333333;
144 const double coef6 = 0.0013888888888889;
145 const double gfac30 = 0.4999489092;
146 const double gfac31 = -0.2473631686;
147 const double gfac32 = 0.321180909;
148 const double gfac33 = -0.3811559346;
149 const double gfac20 = 0.4998436875;
150 const double gfac21 = -0.24249438;
151 const double gfac22 = 0.24642845;
152 const double gfac10 = 0.499093162;
153 const double gfac11 = -0.2152832;
154 const double gfac00 = -0.490;
155
156 double wdif, d2wal, rexpw, /* denom, */ gval, factor, rwval, term;
157 int i, itable, irange;
158
159 367 59 160 84 54 1 0 4 0 0 if (J>maxj) {
160 ExEnv::err()
161 << scprintf("the int_fjt routine has been incorrectly used")
162 << endl;
163 ExEnv::err()
164 << scprintf("J = %d but maxj = %d",J,maxj)
165 << endl;
166 abort();
167 }
168
169 /* Compute an index into the table. */
170 /* The test is needed to avoid floating point exceptions for
171 * large values of wval. */
172 461 37 137 96 46 2 2 1 0 0 if (wval > wval_infinity) {
173 1 0 3 2 1 0 0 0 0 0 itable = itable_infinity;
174 }
175 else {
176 1652 429 1222 137 425 40 3 7 0 0 itable = (int) (10.0 * wval);
177 }
178
179 /* If itable is small enough use the table to compute int_fjttable. */
180 90 0 4 1 20 1 0 2 0 0 if (itable < TABLESIZE) {
181
182 200 0 18 0 41 0 0 1 0 0 wdif = wval - 0.1 * itable;
183
184 /* Compute fjt for J. */
185 int_fjttable[J] = (((((coef6 * gtable[J+6][itable]*wdif
186 + coef5 * gtable[J+5][itable])*wdif
187 + coef4 * gtable[J+4][itable])*wdif
188 + coef3 * gtable[J+3][itable])*wdif
189 + coef2 * gtable[J+2][itable])*wdif
190 - gtable[J+1][itable])*wdif
191 2533 110 503 867 424 11 0 27 3 0 + gtable[J][itable];
192
193 /* Compute the rest of the fjt. */
194 78 2 34 29 30 1 0 0 0 0 d2wal = 2.0 * wval;
195 140 11 47 46 52 1 0 1 0 0 rexpw = exp(-wval);
196 /* denom = 2*J + 1; */
197 219 6 98 23 30 7 2 4 0 0 for (i=J; i>0; i--) {
198 /* denom = denom - 2.0; */
199 818 38 406 226 188 80 7 3 0 0 int_fjttable[i-1] = (d2wal*int_fjttable[i] + rexpw)*denomarray[i];
200 43 1 18 19 11 7 0 0 0 0 }
201 }
202 /* If wval <= 2*J + 36.0, use the following formula. */
203 19 2 2 0 2 0 0 1 0 0 else if (itable <= 20*J + 360) {
204 24 1 8 0 6 0 1 0 0 0 rwval = 1.0/wval;
205 5 0 0 0 0 0 0 0 0 0 rexpw = exp(-wval);
206
207 /* Subdivide wval into 6 ranges. */
208 2 0 4 0 2 0 0 0 0 0 irange = itable/30 - 3;
209 2 0 2 1 1 1 0 0 0 0 if (irange == 1) {
210 10 0 7 3 2 1 0 0 0 0 gval = gfac30 + rwval*(gfac31 + rwval*(gfac32 + rwval*gfac33));
211 int_fjttable[0] = sqrpih*sqrt(rwval) - rexpw*gval*rwval;
212 }
213 2 0 4 1 2 0 0 0 0 0 else if (irange == 2) {
214 2 0 0 2 1 1 0 0 0 0 gval = gfac20 + rwval*(gfac21 + rwval*gfac22);
215 3 0 2 1 0 1 0 0 0 0 int_fjttable[0] = sqrpih*sqrt(rwval) - rexpw*gval*rwval;
216 }
217 5 0 1 1 1 1 0 0 0 0 else if (irange == 3 || irange == 4) {
218 7 0 5 2 1 1 0 0 0 0 gval = gfac10 + rwval*gfac11;
219 0 0 1 0 0 0 0 0 0 0 int_fjttable[0] = sqrpih*sqrt(rwval) - rexpw*gval*rwval;
220 }
221 0 0 2 0 0 0 0 0 0 0 else if (irange == 5 || irange == 6) {
222 gval = gfac00;
223 0 0 0 0 1 0 0 0 0 0 int_fjttable[0] = sqrpih*sqrt(rwval) - rexpw*gval*rwval;
224 }
225 else {
226 0 0 1 0 0 0 0 0 0 0 int_fjttable[0] = sqrpih*sqrt(rwval);
227 }
228
229 /* Compute the rest of the int_fjttable from int_fjttable[0]. */
230 5 0 1 3 0 0 0 0 0 0 factor = 0.5 * rwval;
231 2 0 0 0 0 0 0 0 0 0 term = factor * rexpw;
232 1 1 0 2 0 0 0 0 0 0 for (i=1; i<=J; i++) {
233 14 0 3 13 1 0 0 0 0 0 int_fjttable[i] = factor * int_fjttable[i-1] - term;
234 0 0 0 1 0 0 0 0 0 0 factor = rwval + factor;
235 4 0 1 0 1 0 0 0 0 0 }
236 }
237 /* For large values of wval use this algorithm: */
238 else {
239 49 12 29 2 21 1 0 0 0 0 rwval = 1.0/wval;
240 4 0 4 0 1 1 0 0 0 0 int_fjttable[0] = sqrpih*sqrt(rwval);
241 0 0 0 3 0 0 0 0 0 0 factor = 0.5 * rwval;
242 0 0 0 1 0 0 0 0 0 0 for (i=1; i<=J; i++) {
243 11 0 0 9 0 0 0 0 0 0 int_fjttable[i] = factor * int_fjttable[i-1];
244 0 0 0 1 0 0 0 0 0 0 factor = rwval + factor;
245 172 11 50 28 24 6 0 1 0 0 }
246 }
247 /* printf(" %2d %12.8f %4d %12.8f\n",J,wval,itable,int_fjttable[0]); */
248
249 43 2 26 16 13 1 0 0 0 0 return int_fjttable;
250 67 3 52 33 28 7 0 0 0 0 }
251
252 /////////////////////////////////////////////////////////////////////////////
253
254 // Local Variables:
255 // mode: c++
256 // c-file-style: "CLJ-CONDENSED"
257 // End:
258