1 /* psp.c -
2 author - Patrick Nichols
3 $Id$
4 */
5
6 #include <stdio.h>
7 #include <stdlib.h>
8 #include <string.h>
9
10 #include "name.h"
11 #include "get_word.h"
12 #include "loggrid.h"
13
14 #include "dft.h"
15 #include "atom.h"
16 #include "rtroullier.h"
17 #include "rhamann.h"
18 #include "troullier.h"
19 #include "vanderbilt.h"
20 #include "generate_rho_semicore.h"
21 #include "rpsp.h"
22 #include "psp.h"
23 #define False 0
24 #define True 1
25 #define Max(x,y) ((x>y) ? x : y)
26
27
28 /* atom structure variables */
29
30 static int Nvalence;
31 static int npsp_states;
32 static int *n;
33 static int *l;
34 static int *spin;
35 static int lmax;
36 static double *fill;
37 static double *rcut;
38 static double *peak;
39 static double Zion;
40 static double Total_E,
41 E_Hartree, P_Hartree, E_exchange, P_exchange, E_correlation, P_correlation;
42 static double *eigenvalue;
43 static double **r_psi;
44 static double **r_psi_prime;
45 static double *rho;
46 static double *rho_semicore;
47 static double *drho_semicore;
48 static double r_semicore;
49 static double **V_psp;
50 static char comment[80];
51
52 /* extra Vanderbilt parameters */
53 static double rlocal, clocal;
54 static int ns[10], indx_il[4][10], indx_ijl[4][4][10];
55 static double *Vlocal;
56 static double **r_hard_psi;
57 static double *D0;
58 static double *q;
59
60 /* solver parameters: this is the Kawai-Weare default */
61 static int Solver_Type = Hamann;
62
63
64 /********************************
65 * *
66 * init_RelPsp *
67 * *
68 ********************************/
69
70 void
init_RelPsp(char * filename)71 init_RelPsp (char *filename)
72 {
73 int p, p1, p2;
74 int ltmp, semicore_type;
75 double rctmp;
76 char *w, *tc;
77 FILE *fp;
78
79
80 /* find the psp type */
81 fp = fopen(filename,"r+");
82 w = get_word(fp);
83 while ((w!=NIL) && (strcmp("<pseudopotential>",w)!=0))
84 w = get_word(fp);
85 if (w!=NIL)
86 {
87 w = get_word(fp);
88 if (strcmp("hamann",w)==0) Solver_Type = Hamann;
89 if (strcmp("troullier-martins",w)==0) Solver_Type = Troullier;
90 }else{
91 Solver_Type = Hamann;
92 }
93 fclose(fp);
94
95
96
97 /* find the comment */
98 if (Solver_Type == Hamann) strcpy(comment,"Hamann pseudopotential");
99 if (Solver_Type == Troullier) strcpy(comment,"Troullier-Martins pseudopotential");
100
101 /* set lmax */
102 lmax = lmax_Atom ();
103
104 /* set the number psp projectors */
105 npsp_states = 2 * lmax + 4;
106 /* allocate memory for n,l,fill,rcut,peak, and eigenvalue */
107 n = (int *) malloc ((npsp_states) * sizeof (int));
108 l = (int *) malloc ((npsp_states) * sizeof (int));
109 spin = (int *) malloc ((npsp_states) * sizeof (int));
110 fill = (double *) malloc ((npsp_states) * sizeof (double));
111 rcut = (double *) malloc ((npsp_states) * sizeof (double));
112 peak = (double *) malloc ((npsp_states) * sizeof (double));
113 eigenvalue = (double *) malloc ((npsp_states) * sizeof (double));
114
115
116 /* allocate memory for r_psi, V_psp, and rho */
117 r_psi = (double **) malloc ((npsp_states) * sizeof (double *));
118 r_psi_prime = (double **) malloc ((npsp_states) * sizeof (double *));
119 V_psp = (double **) malloc ((npsp_states) * sizeof (double *));
120 for (p = 0; p < npsp_states; ++p)
121 {
122 r_psi[p] = alloc_LogGrid ();
123 r_psi_prime[p] = alloc_LogGrid ();
124 V_psp[p] = alloc_LogGrid ();
125 }
126 rho = alloc_LogGrid ();
127 rho_semicore = alloc_LogGrid ();
128 drho_semicore = alloc_LogGrid ();
129
130
131 /* get the psp info */
132 if (Solver_Type==Troullier) {
133 Suggested_Param_RelTroullier(&Nvalence, n, l, spin, eigenvalue, fill, rcut);
134 }else{
135 Suggested_Param_RelHamann (&Nvalence, n, l, spin, eigenvalue, fill, rcut);
136 }
137 /* set the number psp projectors */
138 fp = fopen (filename, "r+");
139 w = get_word (fp);
140 while ((w != NIL) && (strcmp ("<npsp-states>", w) != 0))
141 w = get_word (fp);
142
143 if (w != NIL)
144 {
145 w = get_word (fp);
146 while ((w != NIL) && (strcmp ("<end>", w) != 0))
147 {
148 sscanf (w, "%d", <mp);
149 w = get_word (fp);
150 Nvalence = 2 * ltmp + 2;
151 }
152 }
153 fclose (fp);
154
155 /* get rcut */
156 fp = fopen (filename, "r+");
157 w = get_word (fp);
158 while ((w != NIL) && (strcmp ("<rcut>", w) != 0))
159 w = get_word (fp);
160
161 if (w != NIL)
162 {
163 w = get_word (fp);
164 while ((w != NIL) && (strcmp ("<end>", w) != 0))
165 {
166 sscanf (w, "%d", <mp);
167 w = get_word (fp);
168 sscanf (w, "%lf", &rctmp);
169 w = get_word (fp);
170 rcut[2 * ltmp] = rctmp;
171 rcut[2 * ltmp + 1] = rctmp;
172 }
173 }
174 fclose (fp);
175
176 /* get ecut */
177 fp = fopen (filename, "r+");
178 w = get_word (fp);
179 while ((w != NIL) && (strcmp ("<ecut>", w) != 0))
180 w = get_word (fp);
181
182 if (w != NIL)
183 {
184 w = get_word (fp);
185 while ((w != NIL) && (strcmp ("<end>", w) != 0))
186 {
187 sscanf (w, "%d", <mp);
188 w = get_word (fp);
189 sscanf (w, "%lf", &rctmp);
190 w = get_word (fp);
191 eigenvalue[2 * ltmp] = rctmp;
192 eigenvalue[2 * ltmp + 1] = rctmp;
193 }
194 }
195 fclose (fp);
196
197 /* get r_semicore - if zero then no core corrections added */
198 r_semicore = 0.0;
199 fp = fopen (filename, "r+");
200 w = get_word (fp);
201 while ((w != NIL) && (strcmp ("<semicore>", w) != 0))
202 w = get_word (fp);
203
204 if (w != NIL)
205 {
206 w = get_word (fp);
207 while ((w != NIL) && (strcmp ("<end>", w) != 0))
208 {
209 sscanf (w, "%lf", &rctmp);
210 w = get_word (fp);
211 r_semicore = rctmp;
212 }
213 }
214 fclose (fp);
215
216 /* find the semicore_type */
217 semicore_type = 2;
218 fp = fopen (filename, "r+");
219 w = get_word (fp);
220 while ((w != NIL) && (strcmp ("<semicore_type>", w) != 0))
221 w = get_word (fp);
222 if (w != NIL)
223 {
224 w = get_word (fp);
225 if (strcmp ("quadratic", w) == 0)
226 semicore_type = 0;
227 if (strcmp ("louie", w) == 0)
228 semicore_type = 1;
229 if (strcmp ("fuchs", w) == 0)
230 semicore_type = 2;
231 }
232 fclose (fp);
233
234
235
236
237 /* generate non-zero rho_semicore */
238 if (r_semicore > 0.0)
239 {
240 /*
241 printf("\n\n");
242 printf("Generating non-zero semicore density\n");
243 */
244 generate_rho_semicore (semicore_type, rho_core_Atom (), r_semicore,
245 rho_semicore);
246 Derivative_LogGrid (rho_semicore, drho_semicore);
247 }
248
249 /* define the ion charge */
250 /*
251 Zion=0.0;
252 for (p=Ncore_Atom(); p<(Ncore_Atom()+Nvalence_Atom()); ++p)
253 Zion += fill_Atom(p);
254 */
255
256 Zion = Zion_Atom ();
257 for (p = 0; p < (Ncore_Atom ()); ++p)
258 Zion -= fill_Atom (p);
259
260 } /* init_RelPsp */
261
262
263
264 /********************************
265 * *
266 * solve_RelPsp *
267 * *
268 ********************************/
269
270 void
solve_RelPsp()271 solve_RelPsp ()
272 {
273
274 int p, k, Ngrid;
275 double *rgrid;
276
277 /* get loggrid variables */
278 Ngrid = N_LogGrid ();
279 rgrid = r_LogGrid ();
280
281 if (Solver_Type==Troullier) {
282 solve_RelTroullier(Nvalence, n, l, spin, eigenvalue, fill, rcut,
283 r_psi, r_psi_prime, rho, rho_semicore, V_psp,
284 &Total_E,
285 &E_Hartree, &P_Hartree,
286 &E_exchange, &P_exchange, &E_correlation, &P_correlation);
287 }else{
288 solve_RelHamann(Nvalence, n, l, spin, eigenvalue, fill, rcut,
289 r_psi, r_psi_prime, rho, rho_semicore, V_psp,
290 &Total_E,
291 &E_Hartree, &P_Hartree,
292 &E_exchange, &P_exchange, &E_correlation, &P_correlation);
293 }
294 /******************************************************/
295 /* find the outermost peak of the pseudowavefunctions */
296 /******************************************************/
297 for (p = 0; p < Nvalence; ++p)
298 {
299 if (fill[p] != 0.0)
300 {
301 k = Ngrid - 2;
302 while ((r_psi_prime[p][k] * r_psi_prime[p][k + 1] >= 0.0)
303 && (k >= 0))
304 --k;
305 peak[p] = rgrid[k];
306 }
307 else
308 peak[p] = rgrid[Ngrid - 1];
309 }
310 } /* solve_RelPsp */
311
312
313
314 /********************************
315 * *
316 * print_RelPsp *
317 * *
318 ********************************/
319
320 char *
solver_Name_RelPsp()321 solver_Name_RelPsp ()
322 {
323 char *s;
324 if (Solver_Type==Troullier) {
325 s="TroullierMartins";
326 }else{
327 s = "Hamann";
328 }
329 return s;
330 }
331
332
333 void
print_RelPsp(FILE * fp)334 print_RelPsp (FILE * fp)
335 {
336 int i;
337 double dx;
338 fprintf (fp, "\n\n");
339 fprintf (fp, "PSP solver information\n\n");
340 fprintf (fp, "Atom name: %s\n", name_Atom ());
341 fprintf (fp, "Zcharge : %le\n", Zion);
342 fprintf (fp, "Nvalence : %d\n", Nvalence);
343 fprintf (fp, " : restricted calculation\n");
344
345 fprintf (fp, "\n------ Solver information ------\n");
346 fprintf (fp, "solver type : %s\n", solver_Name_RelPsp ());
347 fprintf (fp, "hartree type : %s\n", hartree_Name_DFT ());
348 fprintf (fp, "exchange type : %s\n", exchange_Name_DFT ());
349 if (strcmp (exchange_Name_DFT (), "Dirac") == 0)
350 fprintf (fp, " alpha : %lf\n", Dirac_alpha ());
351 fprintf (fp, "correlation type : %s\n", correlation_Name_DFT ());
352
353 fprintf (fp,
354 "----------------------------------------------------------------------------\n");
355 fprintf (fp, "n\tl\ts\tpopulation\tEcut\t\tRcut\t\tOuter Peak\n");
356
357 for (i = 0; i < Nvalence; ++i)
358 {
359 fprintf (fp, "%d\t%s\t%.1lf\t%.2lf\t\t%le\t%le\t%le\n", (l[i] + 1),
360 spd_Name (l[i]), 0.5 * spin[i], fill[i], eigenvalue[i],
361 rcut[i], peak[i]);
362 }
363 fprintf (fp,
364 "----------------------------------------------------------------------------\n");
365 if (r_semicore > 0.0)
366 {
367 fprintf (fp, "SemiCore Corrections Added\n");
368 fprintf (fp, " rcore : %lf\n", r_semicore);
369 fprintf (fp, " Semicore Charge : %lf\n",
370 Integrate_LogGrid (rho_semicore));
371 fprintf (fp, " Semicore Charge gradient : %lf\n",
372 Integrate_LogGrid (drho_semicore));
373 }
374
375 fprintf (fp, "Pseudopotential ion charge = %lf\n", Zion);
376 dx= - Integrate_LogGrid(rho);
377 fprintf (fp, "Pseudopotential electronic charge= %lf\n",
378 dx);
379 fprintf (fp, "Pseudopotential atom charge = %lf\n",
380 (Zion+dx));
381
382 fprintf (fp, "\nTotal E = %le\n", Total_E);
383 fprintf (fp, "\n");
384
385
386 fprintf (fp, "E_Hartree = %le\n", E_Hartree);
387 fprintf (fp, "<Vh> = %le\n", P_Hartree);
388 fprintf (fp, "\n");
389
390 fprintf (fp, "E_exchange = %le\n", E_exchange);
391 fprintf (fp, "<Vx> = %le\n", P_exchange);
392 fprintf (fp, "\n");
393
394 fprintf (fp, "E_correlation = %le\n", E_correlation);
395 fprintf (fp, "<Vc> = %le\n\n", P_correlation);
396
397
398 } /* print_Atom */
399
400 /********************************
401 * *
402 * set_(solver parameters)_Atom *
403 * *
404 ********************************/
405
406 void
set_Solver_RelPsp(solver)407 set_Solver_RelPsp (solver)
408 int solver;
409 {
410 Solver_Type=solver;
411 fprintf (stdout," RelPsp:: Only Hamann or TM Solver is available\n");
412 fprintf (stdout, "RelPsp::Cannot really set the solver yet!\n");
413 }
414
415
416 int
Vanderbilt_RelPsp()417 Vanderbilt_RelPsp ()
418 {
419 int value;
420
421 value = 0;
422 if (Solver_Type == Vanderbilt)
423 value = 1;
424 return value;
425 }
426
427 int
NormConserving_RelPsp()428 NormConserving_RelPsp ()
429 {
430 int value;
431
432 value = 0;
433 if (Solver_Type == Hamann)
434 value = 1;
435 if (Solver_Type == Troullier)
436 value = 1;
437 return value;
438 }
439
440 /********************************
441 * *
442 * E_Atom *
443 * *
444 ********************************/
445
446 double
E_RelPsp()447 E_RelPsp ()
448 {
449 return Total_E;
450 } /*E_Atom */
451
452
453 double
eigenvalue_RelPsp(int i)454 eigenvalue_RelPsp (int i)
455 {
456 return eigenvalue[i];
457 }
458
459 double *
rho_RelPsp()460 rho_RelPsp ()
461 {
462 return rho;
463 }
464
465 double *
rho_semicore_RelPsp()466 rho_semicore_RelPsp ()
467 {
468 return rho_semicore;
469 }
470
471 double *
drho_semicore_RelPsp()472 drho_semicore_RelPsp ()
473 {
474 return drho_semicore;
475 }
476
477 double
r_semicore_RelPsp()478 r_semicore_RelPsp ()
479 {
480 return r_semicore;
481 }
482
483
484 double *
Beta_RelPsp(int i,int l)485 Beta_RelPsp (int i, int l)
486 {
487 return V_psp[indx_il[i][l]];
488 }
489
490 double *
r_psi_il_RelPsp(int i,int l)491 r_psi_il_RelPsp (int i, int l)
492 {
493 return r_psi[indx_il[i][l]];
494 }
495
496 double *
r_hard_psi_il_RelPsp(int i,int l)497 r_hard_psi_il_RelPsp (int i, int l)
498 {
499 return r_hard_psi[indx_il[i][l]];
500 }
501
502 int
ns_RelPsp(int l)503 ns_RelPsp (int l)
504 {
505 return ns[l];
506 }
507
508 double
D0_RelPsp(int i,int j,int l)509 D0_RelPsp (int i, int j, int l)
510 {
511 return D0[indx_ijl[i][j][l]];
512 }
513
514 double
q_RelPsp(int i,int j,int l)515 q_RelPsp (int i, int j, int l)
516 {
517 return q[indx_ijl[i][j][l]];
518 }
519
520 double *
Vlocal_RelPsp()521 Vlocal_RelPsp ()
522 {
523 return Vlocal;
524 }
525
526
527
528 double *
V_RelPsp(int i)529 V_RelPsp (int i)
530 {
531 return V_psp[i];
532 }
533
534
535 double *
r_psi_RelPsp(int i)536 r_psi_RelPsp (int i)
537 {
538 return r_psi[i];
539 }
540
541 int
n_RelPsp(int i)542 n_RelPsp (int i)
543 {
544 return n[i];
545 }
546
547 int
l_RelPsp(int i)548 l_RelPsp (int i)
549 {
550 return l[i];
551 }
552
553 int
lmax_RelPsp()554 lmax_RelPsp ()
555 {
556 return lmax;
557 }
558
559 double
fill_RelPsp(int i)560 fill_RelPsp (int i)
561 {
562 return fill[i];
563 }
564
565
566 int
Nvalence_RelPsp()567 Nvalence_RelPsp ()
568 {
569 return Nvalence;
570 }
571
572 double
peak_RelPsp(int i)573 peak_RelPsp (int i)
574 {
575 return peak[i];
576 }
577
578 double
rcut_RelPsp(int l)579 rcut_RelPsp (int l)
580 {
581 return rcut[2*l];
582 }
583
584
585 double
rcut_il_RelPsp(int i,int l)586 rcut_il_RelPsp (int i, int l)
587 {
588 return rcut[indx_il[i][l]];
589 }
590
591
592 double
Zion_RelPsp()593 Zion_RelPsp ()
594 {
595 return Zion;
596 }
597
598 int
state_RelPsp(int nt,int lt,int st)599 state_RelPsp (int nt, int lt, int st)
600 {
601 int i;
602 i = 0;
603 while (((nt != n[i]) || (lt != l[i]) || st!=spin[i]) && (i < Nvalence))
604 ++i;
605
606 /* Error */
607 if (i >= Nvalence)
608 printf ("Error: state_RelPsp\n");
609
610 return i;
611 }
612
613
614 char *
comment_RelPsp()615 comment_RelPsp ()
616 {
617 return comment;
618 }
619