1 /*
2 $Id$
3 */
4
5 #include <stdlib.h>
6 #include <stdio.h>
7 #include <string.h>
8 #include "typesf2c.h"
9 #include "name.h"
10 #include "loggrid.h"
11 #include "spline.h"
12 #include "atom.h"
13 #include "psp.h"
14 #include "rpsp.h"
15 #include "debug.h"
16
17 #if defined(CRAY) || defined(CRAY_T3D)
18 #include <fortran.h>
19 #if !defined(__crayx1)
20 #define USE_FCD
21 #endif
22 #endif
23
24 #if (defined(CRAY) || defined(WIN32)) && !defined(__crayx1) &&!defined(__MINGW32__)
25 #define rpspsolve_ RPSPSOLVE
26 #endif
27
rpspsolve_(Integer * print_ptr,Integer * debug_ptr,Integer * lmax_ptr,Integer * locp_ptr,double * rlocal_ptr,const _fcd fcd_sdir_name,Integer * n9,const _fcd fcd_dir_name,Integer * n0,const _fcd fcd_in_filename,Integer * n1,const _fcd fcd_out_filename,Integer * n2)28 void FATR rpspsolve_
29 #if defined(USE_FCD)
30 (Integer * print_ptr,
31 Integer * debug_ptr,
32 Integer * lmax_ptr,
33 Integer * locp_ptr,
34 double *rlocal_ptr,
35 const _fcd fcd_sdir_name,
36 Integer * n9,
37 const _fcd fcd_dir_name,
38 Integer * n0,
39 const _fcd fcd_in_filename,
40 Integer * n1, const _fcd fcd_out_filename, Integer * n2)
41 {
42 char *sdir_name = _fcdtocp (fcd_sdir_name);
43 char *dir_name = _fcdtocp (fcd_dir_name);
44 char *in_filename = _fcdtocp (fcd_in_filename);
45 char *out_filename = _fcdtocp (fcd_out_filename);
46 #else
47
48 (print_ptr, debug_ptr, lmax_ptr, locp_ptr, rlocal_ptr,
49 sdir_name, n9, dir_name, n0, in_filename, n1, out_filename, n2)
50 Integer *print_ptr;
51 Integer *debug_ptr;
52 Integer *lmax_ptr;
53 Integer *locp_ptr;
54 double *rlocal_ptr;
55 char sdir_name[];
56 Integer *n9;
57 char dir_name[];
58 Integer *n0;
59 char in_filename[];
60 Integer *n1;
61 char out_filename[];
62 Integer *n2;
63 {
64 #endif
65
66 int i, j, k, l, p, Nlinear, Nvalence;
67 int debug, print;
68 double *rl, *rhol, **psil, **pspl;
69 double over_fourpi, c, x, y;
70 int Ngrid;
71 double *vall, *rgrid;
72 char name[255];
73 int lmax_out, locp_out, nvh;
74 double rlocal_out, vx;
75
76 FILE *fp;
77
78 int m9 = ((int) (*n9));
79 int m0 = ((int) (*n0));
80 int m1 = ((int) (*n1));
81 int m2 = ((int) (*n2));
82 char *infile = (char *) malloc (m9 + m1 + 5);
83 char *outfile = (char *) malloc (m0 + m2 + 5);
84 char *soutfile = (char *) malloc (m0 + m2 + 9);
85 char *full_filename = (char *) malloc (m9 + 25 + 5);
86
87 print = *print_ptr;
88 debug = *debug_ptr;
89 lmax_out = *lmax_ptr;
90 locp_out = *locp_ptr;
91 rlocal_out = *rlocal_ptr;
92
93 (void) strncpy (infile, sdir_name, m9);
94 infile[m9] = '\0';
95 strcat (infile, "/");
96 infile[m9 + 1] = '\0';
97 strncat (infile, in_filename, m1);
98 infile[m9 + m1 + 1] = '\0';
99
100 (void) strncpy (outfile, dir_name, m0);
101 outfile[m0] = '\0';
102 strcat (outfile, "/");
103 outfile[m0 + 1] = '\0';
104 strncat (outfile, out_filename, m2);
105 outfile[m0 + m2 + 1] = '\0';
106
107 over_fourpi = 0.25/M_PI;
108
109 /********************
110 * we have already solved for the Relativsitic AE wavefunctions using atom
111 * in the call to pspsolve
112 *******************/
113 set_debug_print (debug);
114
115
116 init_RelPsp (infile);
117
118
119 solve_RelPsp ();
120
121
122 if (debug)
123 print_RelPsp (stdout);
124
125
126 init_Linear (infile);
127
128 /* allocate linear meshes */
129 Nvalence = Nvalence_RelPsp ();
130 Nlinear = nrl_Linear ();
131 psil = (double **) malloc (Nvalence * sizeof (double *));
132 pspl = (double **) malloc (Nvalence * sizeof (double *));
133 for (p = 0; p < Nvalence; ++p)
134 {
135 psil[p] = (double *) malloc (Nlinear * sizeof (double));
136 pspl[p] = (double *) malloc (Nlinear * sizeof (double));
137 }
138 rl = (double *) malloc (Nlinear * sizeof (double));
139 rhol = (double *) malloc (Nlinear * sizeof (double));
140
141
142 /* Norm-conserving output */
143 for (p = 0; p < Nvalence; ++p)
144 {
145 Log_to_Linear (r_psi_RelPsp (p), rl, psil[p]);
146 Log_to_Linear_zero (V_RelPsp (p), rl, pspl[p]);
147
148 /* normalize scattering state */
149 if (fill_RelPsp (p) == 0.0)
150 {
151 normalize_Linear (psil[p]);
152 }
153 }
154
155 Log_to_Linear (rho_RelPsp (), rl, rhol);
156
157
158 if (debug)
159 {
160 /* output pseudowavefunctions argv[1].psw */
161 strcpy (name, name_Atom ());
162 strcat (name, ".psw.plt");
163 full_filename[0] = '\0';
164 strncpy (full_filename, sdir_name, m9);
165 full_filename[m9] = '\0';
166 strcat (full_filename, "/");
167 full_filename[m9 + 1] = '\0';
168 strcat (full_filename, name);
169
170 printf ("Outputing pseudowavefunctions: %s\n", full_filename);
171 fp = fopen (full_filename, "w+");
172 for (k = 0; k < Nlinear; ++k)
173 {
174 fprintf (fp, "%12.8lf", rl[k]);
175 for (p = 0; p < Nvalence; ++p)
176 fprintf (fp, " %12.8lf", psil[p][k]);
177 fprintf (fp, "\n");
178 }
179 fclose (fp);
180
181 /* output pseudopotentials argv[1].psp.plt */
182 strcpy (name, name_Atom ());
183 strcat (name, ".psp.plt");
184 full_filename[0] = '\0';
185 strncpy (full_filename, sdir_name, m9);
186 full_filename[m9] = '\0';
187 strcat (full_filename, "/");
188 full_filename[m9 + 1] = '\0';
189 strcat (full_filename, name);
190
191 printf ("Outputing pseudopotentials: %s\n", full_filename);
192 fp = fopen (full_filename, "w+");
193 for (k = 0; k < Nlinear; ++k)
194 {
195 fprintf (fp, "%12.8lf", rl[k]);
196 for (p = 0; p < Nvalence; ++p)
197 fprintf (fp, " %12.8lf", pspl[p][k]);
198 fprintf (fp, "\n");
199 }
200 fclose (fp);
201
202 /* output pseudodensity infile.psd.plt */
203 strcpy (name, name_Atom ());
204 strcat (name, ".psd.plt");
205 full_filename[0] = '\0';
206 strncpy (full_filename, sdir_name, m9);
207 full_filename[m9] = '\0';
208 strcat (full_filename, "/");
209 full_filename[m9 + 1] = '\0';
210 strcat (full_filename, name);
211
212 fp = fopen (full_filename, "w+");
213 for (k = 0; k < Nlinear; ++k)
214 fprintf (fp, "%12.8lf %12.8lf\n", rl[k], rhol[k]);
215 fclose (fp);
216 }
217
218 /* output datafile to be used for Kleinman-Bylander input, argv[1].psp */
219 if (print)
220 {
221 printf (" Creating datafile for Kleinman-Bylander input: %s\n",
222 outfile);
223 }
224
225
226 fp = fopen (outfile, "w+");
227 fprintf (fp, "7 %s\n", name_Atom ());
228 fprintf (fp, "%lf %lf %d %d %d %lf\n", Zion_RelPsp (), Amass_Atom (),
229 lmax_RelPsp (), lmax_out, locp_out, rlocal_out);
230 for (p = 0; p <= lmax_RelPsp (); ++p)
231 fprintf (fp, "%lf ", rcut_RelPsp (p));
232 fprintf (fp, "\n");
233 fprintf (fp, "%d %lf\n", nrl_Linear (), drl_Linear ());
234 fprintf (fp, "%s\n", comment_RelPsp ());
235
236 if (print)
237 {
238 printf (" + Appending pseudopotentials: %s thru %s\n",
239 spd_Name (0), spd_Name (lmax_RelPsp ()));
240 }
241
242 nvh=Nvalence/2;
243 for (k = 0; k < Nlinear; ++k)
244 {
245 fprintf (fp, "%12.8lf", rl[k]);
246 for (p = 0; p < nvh; ++p)
247 {
248 /**************************************************************
249 * Here we output the v average
250 * v_avg(l,r)=((l*v(l-1/2) + (l+1)*v(l+1/2))/(2*l+1)
251 *************************************************************/
252 vx = ((p) * pspl[2 * p][k] + (p+1) * pspl[2 * p + 1][k])/(p+p+1);
253 fprintf (fp, " %12.8lf", vx);
254 }
255 fprintf (fp, "\n");
256 }
257 if (print)
258 {
259 printf (" + Appending pseudowavefunctions: %s thru %s\n",
260 spd_Name (0), spd_Name (lmax_RelPsp ()));
261 }
262
263
264 for (k = 0; k < Nlinear; ++k)
265 {
266 fprintf (fp, "%12.8lf ", rl[k]);
267 for (p = 0; p < Nvalence; ++p)
268 fprintf (fp, " %12.8lf", psil[p][k]);
269 fprintf (fp, "\n");
270 }
271
272 /* output spin-orbit datafile to be used for Kleinman-Bylander input, argv[1].dat */
273 if (print)
274 {
275 printf
276 (" Creating spin-orbit datafile for Kleinman-Bylander input: %s\n",
277 outfile);
278 }
279 if (print)
280 {
281 printf (" + Appending Spin-Orbit pseudopotentials: %s thru %s\n",
282 spd_Name (0), spd_Name (lmax_RelPsp ()));
283 }
284 for (k = 0; k < Nlinear; ++k)
285 {
286 fprintf (fp, "%12.8lf", rl[k]);
287 /**************************************************************
288 * Here we output the v spin orbit
289 * v_spin_orbit(l,r)=2*(v(l+1/2) - v(l-1/2))/(2*l+1)
290 *
291 * so that V_spin_orbit|Psi>= -1/2(1+kappa)*v_spin_orbit|Psi>
292 * its confusing because L*S -> - (1+kappa)/2
293 *************************************************************/
294 for (p = 1; p < nvh;++p)
295 {
296 vx = (pspl[2*p+1][k] - pspl[2*p][k]);
297 vx *= 2. / (2. * p + 1.);
298 fprintf (fp, " %15.8lf", vx);
299 }
300 fprintf (fp, "\n");
301 }
302
303
304 /* append semicore corrections */
305 if (r_semicore_RelPsp () != 0.0)
306 {
307 if (print)
308 {
309 printf (" + Appending semicore density\n");
310 }
311 Log_to_Linear (rho_semicore_RelPsp (), rl, rhol);
312 fprintf (fp, "%lf\n", r_semicore_RelPsp ());
313 for (k = 0; k < Nlinear; ++k)
314 fprintf (fp, "%12.8lf %12.8lf\n", rl[k],
315 fabs (rhol[k] * over_fourpi));
316
317 if (print)
318 {
319 printf (" + Appending semicore density gradient\n");
320 }
321 Log_to_Linear (drho_semicore_RelPsp (), rl, rhol);
322 for (k = 0; k < Nlinear; ++k)
323 fprintf (fp, "%12.8lf %12.8lf\n", rl[k], (rhol[k] * over_fourpi));
324 }
325 fclose (fp);
326
327 /******************************************************************/
328 /******************* output AE information ************************/
329 /******************************************************************/
330
331
332 if (debug)
333 {
334
335 /* output all-electron wavefunctions */
336 printf ("Outputing all-electron wavefunctions:");
337 Ngrid = N_LogGrid ();
338 rgrid = r_LogGrid ();
339 for (p = 0; p < (Ncore_Atom () + Nvalence_Atom ()); ++p)
340 {
341 sprintf (name, "%s.%1d%s%s", name_Atom (), n_Atom (p),
342 spd_Name (l_Atom (p)), spin_Name (p));
343 full_filename[0] = '\0';
344 strncpy (full_filename, sdir_name, m9);
345 full_filename[m9] = '\0';
346 strcat (full_filename, "/");
347 full_filename[m9 + 1] = '\0';
348 strcat (full_filename, name);
349
350 printf (" %s", full_filename);
351 fp = fopen (full_filename, "w+");
352 for (k = 0; k < Ngrid; ++k)
353 fprintf (fp, "%12.8lf %12.8lf\n", rgrid[k], r_psi_Atom (p)[k]);
354 fclose (fp);
355 }
356 printf ("\n");
357
358 /* output density argv[1].dns */
359 Ngrid = N_LogGrid ();
360 rgrid = r_LogGrid ();
361 strcpy (name, name_Atom ());
362 strcat (name, ".dns.plt");
363 full_filename[0] = '\0';
364 strncpy (full_filename, sdir_name, m9);
365 full_filename[m9] = '\0';
366 strcat (full_filename, "/");
367 full_filename[m9 + 1] = '\0';
368 strcat (full_filename, name);
369
370 printf ("Outputing atom density: %s\n", full_filename);
371 fp = fopen (full_filename, "w+");
372 for (k = 0; k < Ngrid; ++k)
373 fprintf (fp, "%12.8lf %12.8lf\n", rgrid[k], rho_Atom ()[k]);
374 fclose (fp);
375
376 /* output core density argv[1].cdns */
377 Ngrid = N_LogGrid ();
378 rgrid = r_LogGrid ();
379 strcpy (name, name_Atom ());
380 strcat (name, ".cdns.plt");
381 full_filename[0] = '\0';
382 strncpy (full_filename, sdir_name, m9);
383 full_filename[m9] = '\0';
384 strcat (full_filename, "/");
385 full_filename[m9 + 1] = '\0';
386 strcat (full_filename, name);
387
388 printf ("Outputing core density: %s\n", full_filename);
389 fp = fopen (full_filename, "w+");
390 for (k = 0; k < Ngrid; ++k)
391 fprintf (fp, "%12.8lf %12.8lf\n", rgrid[k], rho_core_Atom ()[k]);
392 fclose (fp);
393
394 /* output core density gradient argv[1].cddns */
395 vall = alloc_LogGrid ();
396 Derivative_LogGrid (rho_core_Atom (), vall);
397 Ngrid = N_LogGrid ();
398 rgrid = r_LogGrid ();
399 strcpy (name, name_Atom ());
400 strcat (name, ".cddns.plt");
401 full_filename[0] = '\0';
402 strncpy (full_filename, sdir_name, m9);
403 full_filename[m9] = '\0';
404 strcat (full_filename, "/");
405 full_filename[m9 + 1] = '\0';
406 strcat (full_filename, name);
407 printf ("Outputing core density gradient: %s\n", full_filename);
408 fp = fopen (full_filename, "w+");
409 for (k = 0; k < Ngrid; ++k)
410 fprintf (fp, "%12.8lf %12.8lf\n", rgrid[k], vall[k]);
411 fclose (fp);
412 dealloc_LogGrid (vall);
413
414 /* output semicore density infile.sdns */
415 Ngrid = N_LogGrid ();
416 rgrid = r_LogGrid ();
417 strcpy (name, name_Atom ());
418 strcat (name, ".sdns.plt");
419 full_filename[0] = '\0';
420 strncpy (full_filename, sdir_name, m9);
421 full_filename[m9] = '\0';
422 strcat (full_filename, "/");
423 full_filename[m9 + 1] = '\0';
424 strcat (full_filename, name);
425
426 printf ("Outputing semicore density: %s\n", full_filename);
427 fp = fopen (full_filename, "w+");
428 for (k = 0; k < Ngrid; ++k)
429 fprintf (fp, "%12.8lf %12.8lf\n", rgrid[k], rho_semicore_RelPsp ()[k]);
430 fclose (fp);
431
432 /* output semicore density gradient infile.sddns */
433 Ngrid = N_LogGrid ();
434 rgrid = r_LogGrid ();
435 strcpy (name, name_Atom ());
436 strcat (name, ".sddns.plt");
437 full_filename[0] = '\0';
438 strncpy (full_filename, sdir_name, m9);
439 full_filename[m9] = '\0';
440 strcat (full_filename, "/");
441 full_filename[m9 + 1] = '\0';
442 strcat (full_filename, name);
443
444 printf ("Outputing semicore density gradient: %s\n", full_filename);
445 fp = fopen (full_filename, "w+");
446 for (k = 0; k < Ngrid; ++k)
447 {
448 vx=drho_semicore_RelPsp()[k];
449 fprintf (fp, "%12.8lf %12.8lf\n", rgrid[k], vx);
450 }
451 fclose (fp);
452
453 /* output all-electron potential infile.pot */
454 vall = Vall_Atom ();
455 Ngrid = N_LogGrid ();
456 rgrid = r_LogGrid ();
457 strcpy (name, name_Atom ());
458 strcat (name, ".pot.plt");
459 full_filename[0] = '\0';
460 strncpy (full_filename, sdir_name, m9);
461 full_filename[m9] = '\0';
462 strcat (full_filename, "/");
463 full_filename[m9 + 1] = '\0';
464 strcat (full_filename, name);
465
466 printf ("Outputing all-electron potential(non-screened): %s\n",
467 full_filename);
468 fp = fopen (full_filename, "w+");
469
470 c = 0.0;
471 for (k = 0; k < Ngrid; ++k)
472 {
473 x = rgrid[k] / rcut_RelPsp (0);
474 x = pow (x, 3.5);
475 x = exp (-x);
476 y = 1.0 - x;
477 fprintf (fp, "%12.8lf %12.8lf %12.8lf\n", rgrid[k], vall[k],
478 y * vall[k] + c * x);
479 }
480 fclose (fp);
481 }
482
483 /* free malloc memory */
484 free (infile);
485 free (outfile);
486 free (full_filename);
487 for (p = 0; p < Nvalence; ++p)
488 {
489 free (psil[p]);
490 free (pspl[p]);
491 }
492 free (psil);
493 free (pspl);
494 free (rl);
495 free (rhol);
496 end_Linear ();
497 fflush (stdout);
498 return;
499 }
500