1 /*
2 $Id$
3 pauli.c - 6/9/95
4 author - Eric Bylaska
5
6 This file contains routines for integrating the radial
7 Pauli equation.
8
9 */
10
11 #include <stdio.h>
12 #include "loggrid.h"
13 #include "pred_cor.h"
14 #include "pauli.h"
15
16 #define Max_Iterations 100
17 #define tolerance 1.0e-10
18 #define Corrector_Iterations 6
19
20 #define Min(x,y) ((x<y) ? x : y)
21 #define Max(x,y) ((x>y) ? x : y)
22
23
24
25 /**********************************
26 * *
27 * R_Pauli *
28 * *
29 **********************************/
30
R_Pauli(n,l,Z,v,mch,Eig,u,uprime)31 void R_Pauli(n,l,Z,v,mch,Eig,u,uprime)
32 int n,l;
33 double Z;
34 double v[];
35 int *mch;
36 double *Eig;
37 double u[],
38 uprime[];
39
40 {
41 int i,j,
42 iteration,
43 node,
44 match,
45 Ninf, Ngrid;
46
47 double E, de,
48 Emax,
49 Emin,
50 log_amesh,
51 log_amesh2,
52 fss,gamma,
53 L2,L0,
54 r2,
55 sum, a, scale, m1scale,
56 uout,upout,upin,
57 *r,
58 *f_upp,
59 *dv, *frp, *fr,
60 *upp;
61
62 /* define eigenvalues */
63 E = *Eig;
64 L2 = ((double) (l*(l+1)));
65
66 fss = 1.0/137.03602;
67 fss = fss*fss;
68 if (l == 0)
69 gamma = sqrt(1.0 - fss*Z*Z);
70 else
71 {
72 L0 = (double) l;
73 gamma = ( sqrt(L0*L0 - fss*Z*Z)*L0
74 + sqrt((L0+1.0)*(L0+1) - fss*Z*Z)*(L0+1.0)
75 )/(2.0*L0+1.0);
76 }
77
78
79 /* define log grid parameters */
80 Ngrid = N_LogGrid();
81 log_amesh = log_amesh_LogGrid();
82 log_amesh2 = log_amesh*log_amesh;
83
84 r = r_LogGrid();
85 f_upp = alloc_LogGrid();
86 upp = alloc_LogGrid();
87 fr = alloc_LogGrid();
88 frp = alloc_LogGrid();
89 dv = alloc_LogGrid();
90
91
92 /* set up bounds for eigenvalue */
93 Emax = v[Ngrid-1] + 0.5*L2/(r[Ngrid-1]*r[Ngrid-1]);
94 Emin = 0.0;
95 for (i=0; i<Ngrid; ++i)
96 {
97 r2 = r[i];
98 r2 = r2*r2;
99 Emin = Min(Emin, (v[i] + 0.5*L2/r2));
100 }
101 if (E > Emax) E = 1.25*Emax;
102 if (E < Emin) E = 0.75*Emin;
103 if (E > Emax) E = 0.5*(Emax+Emin);
104
105 for (i=0; i<4; ++i)
106 {
107 u[i] = 0.0;
108 uprime[i] = 0.0;
109 upp[i] = 0.0;
110 }
111 iteration = 0;
112 while (iteration < Max_Iterations)
113 {
114 ++iteration;
115 /* define f_upp */
116 for (i=0; i<Ngrid; ++i)
117 {
118 r2 = r[i];
119 r2 = r2*r2;
120 f_upp[i] = log_amesh2*(L2 + 2.0*(v[i] - E)*r2);
121 }
122 /* define dV/dr */
123 Derivative_LogGrid(v,dv);
124 /*
125 dv[0] = Derivative5_1(0,v)/(log_amesh*r[0]);
126 dv[1] = Derivative5_2(1,v)/(log_amesh*r[1]);
127 for (i=2; i<Ngrid-2; ++i)
128 dv[i] = Derivative5_3(i,v)/(log_amesh*r[i]);
129 dv[Ngrid-2] = Derivative5_4(Ngrid-2,v)/(log_amesh*r[Ngrid-2]);
130 dv[Ngrid-1] = Derivative5_5(Ngrid-1,v)/(log_amesh*r[Ngrid-1]);
131 */
132
133
134
135 for (i=0; i<Ngrid; ++i)
136 {
137 r2 = r[i]*r[i];
138 fr[i] = log_amesh2*r2
139 *( -fss*(v[i]-E)*(v[i]-E)
140 + 0.5*fss*dv[i]/(r[i]*(1.0 + 0.5*fss*(E-v[i])))
141 );
142 frp[i] = -log_amesh*r[i]*0.5*fss*dv[i]
143 /(1.0 + 0.5*fss*(E-v[i]));
144 }
145
146 /* find the classical turning point, */
147 /* which is used for matching */
148 match = Ngrid-1;
149 while (f_upp[match-1]*f_upp[match] > 0.0)
150 {
151 match = match - 1;
152
153 if (match < 2)
154 {
155 printf("Error in R_Pauli: no turning point\n");
156
157 /* deallocate memory */
158 dealloc_LogGrid(f_upp);
159 dealloc_LogGrid(upp);
160 dealloc_LogGrid(fr);
161 dealloc_LogGrid(frp);
162 dealloc_LogGrid(dv);
163
164 return;
165 }
166 }
167
168
169
170 /* set the boundry condition near zero */
171 m1scale = 1.0;
172 for (i=0; i<(n-l-1); ++i)
173 m1scale *= -1.0;
174 for (i=0; i<4; ++i)
175 {
176 u[i] = m1scale*pow(r[i],gamma);
177 uprime[i] = log_amesh*gamma*u[i];
178 upp[i] = (log_amesh+frp[i])*uprime[i]
179 + (f_upp[i]+fr[i])* u[i];
180 }
181
182 /* integrate from 0 to match */
183 node = 0;
184 for (i=3; i<match; ++i)
185 {
186 /* predictors */
187 u[i+1] = Predictor_Out(i,u,uprime);
188 uprime[i+1] = Predictor_Out(i,uprime,upp);
189
190 /* correctors */
191 for (j=0; j<Corrector_Iterations; ++j)
192 {
193 upp[i+1] = (log_amesh + frp[i+1])*uprime[i+1]
194 + (f_upp[i+1] + fr[i+1])*u[i+1];
195 uprime[i+1] = Corrector_Out(i,uprime,upp);
196 u[i+1] = Corrector_Out(i,u,uprime);
197 }
198
199 /* finding nodes */
200 if (u[i+1]*u[i] <= 0) node = node + 1;
201 }
202 uout = u[match];
203 upout = uprime[match];
204
205 /* not enough nodes in u */
206 if ((node-n+l+1) < 0)
207 {
208 Emin = E;
209 E = 0.5*(Emin+Emax);
210 }
211 /* too many nodes in u */
212 else if ((node-n+l+1) > 0)
213 {
214 Emax = E;
215 E = 0.5*(Emin+Emax);
216 }
217 /* number of nodes ok, start integration inward */
218 else
219 {
220
221 /* find infinity */
222 Ninf = match + floor(2.3/log_amesh);
223 if ((Ninf+5) > Ngrid) Ninf = Ngrid - 5;
224
225 /* define boundry near infinity */
226 a = sqrt( L2/(r[Ninf]*r[Ninf]) + 2.0*(v[Ninf]-E) );
227 for (i=Ninf; i<=(Ninf+4); ++i)
228 {
229 u[i] = exp(-a*(r[i]-r[Ninf]));
230 uprime[i] = -r[i]*log_amesh*a*u[i];
231 upp[i] = (log_amesh + frp[i])*uprime[i]
232 + (f_upp[i] + fr[i])*u[i];
233 }
234
235 /* integrate from infinity to match */
236 for (i=Ninf; i>=(match+1); --i)
237 {
238 /* predictors */
239 u[i-1] = Predictor_In(i,u,uprime);
240 uprime[i-1] = Predictor_In(i,uprime,upp);
241
242 /* Correctors */
243 for (j=0; j<Corrector_Iterations; ++j)
244 {
245 upp[i-1] = (log_amesh + frp[i-1])*uprime[i-1]
246 + (f_upp[i-1] + fr[i-1])*u[i-1];
247 uprime[i-1] = Corrector_In(i,uprime,upp);
248 u[i-1] = Corrector_In(i,u,uprime);
249 }
250 }
251
252 /* make the outside u, match the inside u */
253 scale = uout/u[match];
254 for (i=match; i<=Ninf; ++i)
255 {
256 u[i] = scale*u[i];
257 uprime[i] = scale*uprime[i];
258 }
259 upin = uprime[match];
260
261 /* Find Integral(u**2) */
262 sum = Norm_LogGrid(Ninf,gamma,u);
263
264
265
266
267 sum = 1.0/sqrt(sum);
268 uout = sum*uout;
269 upout = sum*upout;
270 upin = sum*upin;
271 for (i=0; i<=Ninf; ++i)
272 {
273 u[i] = sum*u[i];
274 uprime[i] = sum*uprime[i];
275 }
276 for (i=Ninf+1; i<Ngrid; ++i)
277 {
278 u[i] = 0.0;
279 uprime[i] = 0.0;
280 }
281
282 /* figure out new eigenvalue */
283 de = 0.5*uout*(upout-upin)/(log_amesh*r[match]);
284
285 /* eigenvalue is converged, exit */
286 if (fabs(de) < (Max(fabs(E),0.2)*tolerance))
287 {
288 *mch = match;
289 *Eig = E;
290
291 /* deallocate memory */
292 dealloc_LogGrid(f_upp);
293 dealloc_LogGrid(upp);
294 dealloc_LogGrid(fr);
295 dealloc_LogGrid(frp);
296 dealloc_LogGrid(dv);
297
298 return;
299 }
300
301 if (de > 0.0)
302 Emin = E;
303 else
304 Emax = E;
305 E = E + de;
306 if ( (E > Emax) || (E < Emin))
307 E = 0.5*(Emin+Emax);
308
309 } /* nodes ok */
310 } /* while */
311
312 printf("Error R_Pauli: More than %d iterations\n",Max_Iterations);
313 *mch = match;
314 *Eig = E;
315
316 /* deallocate memory */
317 dealloc_LogGrid(f_upp);
318 dealloc_LogGrid(upp);
319 dealloc_LogGrid(fr);
320 dealloc_LogGrid(frp);
321 dealloc_LogGrid(dv);
322
323 return;
324
325 } /* R_Pauli */
326
327
328
329 /**********************************
330 * *
331 * R_Pauli_Fixed_E *
332 * *
333 **********************************/
334
R_Pauli_Fixed_E(n,l,Z,v,match,E,u,uprime)335 void R_Pauli_Fixed_E(n,l,Z,v,match,E,u,uprime)
336 int n,l;
337 double Z;
338 double v[];
339 int match;
340 double E;
341 double u[],
342 uprime[];
343
344 {
345 int i,j,
346 node,
347 Ngrid;
348
349 double log_amesh,
350 log_amesh2,
351 fss,gamma,
352 L2,L0,
353 r2,
354 sum,
355 *r,
356 *f_upp,
357 *dv, *frp, *fr,
358 *upp;
359
360 /* define eigenvalues */
361 L2 = ((double) (l*(l+1)));
362
363 fss = 1.0/137.03602;
364 fss = fss*fss;
365 if (l == 0)
366 gamma = sqrt(1.0 - fss*Z*Z);
367 else
368 {
369 L0 = (double) l;
370 gamma = ( sqrt(L0*L0 - fss*Z*Z)*L0
371 + sqrt((L0+1.0)*(L0+1) - fss*Z*Z)*(L0+1.0)
372 )/(2.0*L0+1.0);
373 }
374
375
376 /* define log grid parameters */
377 Ngrid = N_LogGrid();
378 log_amesh = log_amesh_LogGrid();
379 log_amesh2 = log_amesh*log_amesh;
380
381 r = r_LogGrid();
382 f_upp = alloc_LogGrid();
383 upp = alloc_LogGrid();
384 fr = alloc_LogGrid();
385 frp = alloc_LogGrid();
386 dv = alloc_LogGrid();
387
388
389 for (i=0; i<4; ++i)
390 {
391 u[i] = 0.0;
392 uprime[i] = 0.0;
393 upp[i] = 0.0;
394 }
395
396
397
398 /* define f_upp */
399 for (i=0; i<Ngrid; ++i)
400 {
401 r2 = r[i];
402 r2 = r2*r2;
403 f_upp[i] = log_amesh2*(L2 + 2.0*(v[i] - E)*r2);
404 }
405 /* define dV/dr */
406 Derivative_LogGrid(v,dv);
407
408
409 for (i=0; i<Ngrid; ++i)
410 {
411 r2 = r[i]*r[i];
412 fr[i] = log_amesh2*r2
413 *( -fss*(v[i]-E)*(v[i]-E)
414 + 0.5*fss*dv[i]/(r[i]*(1.0 + 0.5*fss*(E-v[i])))
415 );
416 frp[i] = -log_amesh*r[i]*0.5*fss*dv[i]
417 /(1.0 + 0.5*fss*(E-v[i]));
418 }
419
420
421
422 /* set the boundry condition near zero */
423 for (i=0; i<4; ++i)
424 {
425 u[i] = pow(r[i],gamma);
426 uprime[i] = log_amesh*gamma*u[i];
427 upp[i] = (log_amesh+frp[i])*uprime[i]
428 + (f_upp[i]+fr[i])* u[i];
429 }
430
431 /* integrate from 0 to match */
432 node = 0;
433 for (i=3; i<match; ++i)
434 {
435 /* predictors */
436 u[i+1] = Predictor_Out(i,u,uprime);
437 uprime[i+1] = Predictor_Out(i,uprime,upp);
438
439 /* correctors */
440 for (j=0; j<Corrector_Iterations; ++j)
441 {
442 upp[i+1] = (log_amesh + frp[i+1])*uprime[i+1]
443 + (f_upp[i+1] + fr[i+1])*u[i+1];
444 uprime[i+1] = Corrector_Out(i,uprime,upp);
445 u[i+1] = Corrector_Out(i,u,uprime);
446 }
447
448 }
449
450 /* Find Integral(u**2) */
451 sum = Norm_LogGrid(match,gamma,u);
452 sum = 1.0/sqrt(sum);
453
454 for (i=0; i<=match; ++i)
455 {
456 u[i] = sum*u[i];
457 uprime[i] = sum*uprime[i];
458 }
459 for (i=match+1; i<Ngrid; ++i)
460 {
461 u[i] = 0.0;
462 uprime[i] = 0.0;
463 }
464
465 /* deallocate memory */
466 dealloc_LogGrid(f_upp);
467 dealloc_LogGrid(upp);
468 dealloc_LogGrid(fr);
469 dealloc_LogGrid(frp);
470 dealloc_LogGrid(dv);
471
472 return;
473
474 } /* R_Pauli_Fixed_E */
475