1 /*
2 $Id$
3
4 LogGrid.c - 6/9/95
5 author - Eric Bylaska
6
7 This file contains the data structure for handeling numerics
8 on a logarithmic grid. The grid is defined from 0 to 45.0,
9 with the grid points defined by:
10
11 r(i) = (a**i)*r0
12
13 with r0 = 0.00625/Z
14 a = 1.0247
15 i = 0,1,...,N; N = log(7200.0*Z)/AL
16
17 */
18
19 #include <stdio.h>
20 #include "grids.h"
21 #include "pred_cor.h"
22 #include "loggrid.h"
23
24 /***********************/
25 /* LogGrid data structure */
26 /***********************/
27
28
29 /* Hamman definitions */
30 /*static double amesh = 1.0247; */
31 /*static double r0Z = 0.00625; */
32 /*static double Lmax = 45.0; */
33
34 /* My definitions */
35 /*static double Lmax = 45.0; */
36 /*static double amesh = 1.0050; */
37 /*static double r0Z = 0.00025; */
38
39 static double Lmax = 45.0;
40 static double amesh = 1.0050;
41 static double r0Z = 0.00025;
42
43
44
45 static int N;
46 static double r0;
47 static double log_amesh;
48 static double *rgrid;
49
50
51
52
53 /********************************
54 * *
55 * init_LogGrid *
56 * *
57 ********************************/
58
59 /* sets up the log grid data structure,
60 and returns the number grid points, N.
61
62 Entry - Z: charge of the system
63 */
64
init_LogGrid(Z)65 void init_LogGrid(Z)
66 double Z;
67 {
68 int i;
69
70
71
72 /* define r0 */
73 r0 = r0Z/Z;
74
75 /* define log(amesh) */
76 log_amesh = log(amesh);
77
78 /* find N */
79 N = log(Lmax/r0)/log_amesh;
80
81 /* Initialize the grids data structure */
82 init_Grids(N);
83
84 /* allocate rgrid and a tmp grids */
85 rgrid = alloc_Grid();
86
87
88 /* define rgrid */
89 rgrid[0] = r0;
90 for (i=1; i<N; ++i)
91 rgrid[i] = amesh*rgrid[i-1];
92
93 } /* init_LogGrid */
94
95 /********************************
96 * *
97 * end_LogGrid *
98 * *
99 ********************************/
end_LogGrid()100 void end_LogGrid()
101 {
102 dealloc_Grid(rgrid);
103 end_Grids();
104 }
105
106
107 /********************************
108 * *
109 * alloc_LogGrid *
110 * *
111 ********************************/
112
113 /* returns the pointer to a loggrid array.
114
115 */
alloc_LogGrid()116 double *alloc_LogGrid()
117 {
118 double *tt;
119
120 tt = alloc_Grid();
121 Zero_LogGrid(tt);
122 return tt;
123
124 } /* alloc_LogGrid */
125
126
127 /********************************
128 * *
129 * dealloc_LogGrid *
130 * *
131 ********************************/
132
133 /* deallocates the LogGrid array
134
135 */
dealloc_LogGrid(grid)136 void dealloc_LogGrid(grid)
137 double *grid;
138 {
139
140 dealloc_Grid(grid);
141
142 } /* dealloc_LogGrid */
143
144
145
146 /********************************
147 * *
148 * r_LogGrid *
149 * *
150 ********************************/
151
152 /* returns the pointer to the rgrid array.
153
154 */
r_LogGrid()155 double *r_LogGrid()
156 {
157 return rgrid;
158
159 } /* r_LogGrid */
160
161 /********************************
162 * *
163 * index_r_LogGrid *
164 * *
165 ********************************/
166
167 /* returns the pointer to the rgrid array.
168
169 */
index_r_LogGrid(r)170 int index_r_LogGrid(r)
171 double r;
172 {
173 int index;
174
175 index = log(r/r0)/log_amesh;
176 return index;
177
178 } /* index_r_LogGrid */
179
180
181
182 /********************************
183 * *
184 * N_LogGrid *
185 * *
186 ********************************/
187
188 /* returns the size of a log grid.
189
190 */
191
N_LogGrid()192 int N_LogGrid()
193 {
194 return N;
195
196 } /* N_LogGrid */
197
198
199
200 /********************************
201 * *
202 * log_amesh_LogGrid *
203 * *
204 ********************************/
205
206 /* returns the value log(amesh).
207
208 */
209
log_amesh_LogGrid()210 double log_amesh_LogGrid()
211 {
212 return log_amesh;
213
214 } /* log_amesh_LogGrid */
215
216
217
218 /********************************
219 * *
220 * amesh_LogGrid *
221 * *
222 ********************************/
223
224 /* returns the value (amesh).
225
226 */
227
amesh_LogGrid()228 double amesh_LogGrid()
229 {
230 return(amesh);
231
232 } /* amesh_LogGrid */
233
234
235
236 /********************************
237 * *
238 * Integrate_LogGrid *
239 * *
240 ********************************/
241
Integrate_LogGrid(f)242 double Integrate_LogGrid(f)
243 double f[];
244 {
245 int i;
246 double sum;
247
248 /*
249 sum = ( 9.0*f[0]*(rgrid[0]*rgrid[0]*rgrid[0])
250 + 23.0*f[1]*(rgrid[1]*rgrid[1]*rgrid[1])
251 + 28.0*f[2]*(rgrid[2]*rgrid[2]*rgrid[2])
252 )/28.0;
253 */
254 sum = ( 9.0*f[0]*(rgrid[0]*rgrid[0]*rgrid[0])
255 + 28.0*f[1]*(rgrid[1]*rgrid[1]*rgrid[1])
256 + 23.0*f[2]*(rgrid[2]*rgrid[2]*rgrid[2])
257 )/24.0;
258 for (i=3; i<N; ++i)
259 {
260 sum += f[i]*(rgrid[i]*rgrid[i]*rgrid[i]);
261 }
262 sum = log_amesh*sum + f[0]*(rgrid[0]*rgrid[0]*rgrid[0])/3.0;
263
264 return sum;
265
266 }
267
268 /********************************
269 * *
270 * Integrate_LogGrid_na_nb *
271 * *
272 ********************************/
273
Integrate_LogGrid_na_nb(int na,int nb,double f[])274 double Integrate_LogGrid_na_nb(int na, int nb, double f[])
275 {
276 int i;
277 double sum;
278
279 sum = ( 9.0*f[na ]*(rgrid[na ])
280 + 28.0*f[na+1]*(rgrid[na+1])
281 + 23.0*f[na+2]*(rgrid[na+2])
282 )/24.0;
283 for (i=na+3; i<=(nb-3); ++i)
284 {
285 sum += f[i]*(rgrid[i]);
286 }
287 sum += ( 23.0*f[nb-2]*(rgrid[nb-2])
288 + 28.0*f[nb-1]*(rgrid[nb-1])
289 + 9.0*f[nb ]*(rgrid[nb ])
290 )/24.0;
291
292 sum = log_amesh*sum;
293
294 return sum;
295
296 }
297
298 /********************************
299 * *
300 * Zero_LogGrid *
301 * *
302 ********************************/
303
Zero_LogGrid(grid)304 void Zero_LogGrid(grid)
305 double *grid;
306 {
307 int i;
308
309 for (i=0; i<N; ++i)
310 grid[i] = 0.0;
311
312 } /* Zero_LogGrid */
313
314 /********************************
315 * *
316 * Copy_LogGrid *
317 * *
318 ********************************/
319
Copy_LogGrid(gridnew,gridold)320 void Copy_LogGrid(gridnew,gridold)
321
322 double *gridnew,
323 *gridold;
324 {
325 int i;
326 for (i=0; i<N; ++i)
327 gridnew[i] = gridold[i];
328 }
329
330
331 /**********************************
332 * *
333 * Norm_LogGrid *
334 * *
335 **********************************/
336
337 /* This routine calculates the Norm
338 of a wavefunction assuming that
339 the wavefunction decays like an
340 exponential as r-> infinity
341 and approaches zero as r**(2*gamma)
342
343 */
Norm_LogGrid(M,gamma,u)344 double Norm_LogGrid(M,gamma,u)
345 int M;
346 double gamma;
347 double u[];
348 {
349 int i;
350 double r0,sum;
351
352
353 /* Find Integral(u**2) */
354 r0 = rgrid[0]/sqrt(amesh);
355 sum = pow(r0,(2.0*gamma+1.0))/(2.0*gamma+1.0);
356 for (i=0; i<=(M-3); ++i)
357 sum += log_amesh*rgrid[i]*(u[i]*u[i]);
358
359 sum += log_amesh*( 23.0*rgrid[M-2]*(u[M-2]*u[M-2])
360 + 28.0*rgrid[M-1]*(u[M-1]*u[M-1])
361 + 9.0*rgrid[M] *(u[M] *u[M]))/24.0;
362
363 return sum;
364
365 } /* Norm_LogGrid */
366
367
368 /**********************************
369 * *
370 * Integrate2_LogGrid *
371 * *
372 **********************************/
373
374 /* This routine calculates the integral fo function f,
375 assuming that f approaches zero as r**nu, and decays
376 like an exponential as r-> infinity
377
378 */
Integrate2_LogGrid(M,nu,f)379 double Integrate2_LogGrid(M,nu,f)
380 int M;
381 double nu;
382 double f[];
383 {
384 int i;
385 double r0,sum;
386
387
388 /* Find Integral of f */
389 r0 = rgrid[0]/sqrt(amesh);
390 sum = pow(r0,(nu+1.0))/(nu+1.0);
391 for (i=0; i<=(M-3); ++i)
392 sum += log_amesh*rgrid[i]*(f[i]);
393
394 sum += log_amesh*( 23.0*rgrid[M-2]*(f[M-2])
395 + 28.0*rgrid[M-1]*(f[M-1])
396 + 9.0*rgrid[M] *(f[M] ))/24.0;
397
398 return sum;
399
400 } /* Integrate2_LogGrid */
401
402 /********************************
403 * *
404 * Derivative_LogGrid *
405 * *
406 ********************************/
407
Derivative_LogGrid(f,df)408 void Derivative_LogGrid(f,df)
409
410 double f[],
411 df[];
412 {
413 int i;
414
415 /* define dV/dr */
416 df[0] = Derivative5_1(0,f)/(log_amesh*rgrid[0]);
417 df[1] = Derivative5_2(1,f)/(log_amesh*rgrid[1]);
418 for (i=2; i<N-2; ++i)
419 df[i] = Derivative5_3(i,f)/(log_amesh*rgrid[i]);
420 df[N-2] = Derivative5_4(N-2,f)/(log_amesh*rgrid[N-2]);
421 df[N-1] = Derivative5_5(N-1,f)/(log_amesh*rgrid[N-1]);
422
423 } /* Derivative_LogGrid */
424
Plot_LogGrid(char * name,double * f)425 void Plot_LogGrid(char *name, double *f)
426 {
427 int i;
428 FILE *fp;
429
430 fp = fopen(name,"w+");
431 for (i=0; i<N; ++i)
432 fprintf(fp,"%le %le\n",rgrid[i],f[i]);
433 fclose(fp);
434 }
435