analoguous psp1d calls placed in rpsp.c ...
This commit is contained in:
Pat Nichols 2007-04-10 23:30:34 +00:00
parent a9d5c0c6c7
commit 86ed47b7d2
6 changed files with 461 additions and 417 deletions

View file

@ -1,6 +1,6 @@
/* atom.c -
author - Eric Bylaska
$Id: atom.c,v 1.11 2007-04-10 18:43:01 d3p708 Exp $
$Id: atom.c,v 1.12 2007-04-10 23:30:33 d3p708 Exp $
*/
#include <stdio.h>
@ -778,6 +778,21 @@ state_Atom (int nt, int lt)
return i;
}
int
state_RelAtom (int nt, int lt, int st)
{
int i;
i = 0;
while (((nt != n[i]) || (lt != l[i]) || (st!=s2[i])) && (i <= Ncv))
++i;
/* Error */
if (i > Ncv)
printf ("Error: state_RelAtom\n");
return i;
}
char *
name_Atom ()
{

View file

@ -35,6 +35,7 @@ extern int turning_point_Atom (int);
extern double Zion_Atom ();
extern double Amass_Atom ();
extern int state_Atom (int, int);
extern int state_RelAtom (int, int, int);
extern char *name_Atom ();
extern int solver_Type_Atom();
/* used for setting solver parameters */

View file

@ -1,5 +1,5 @@
/*
$Id: pspsolve.c,v 1.16 2007-04-10 18:43:01 d3p708 Exp $
$Id: pspsolve.c,v 1.17 2007-04-10 23:30:34 d3p708 Exp $
*/
#include <stdlib.h>
#include <stdio.h>
@ -124,7 +124,8 @@ void FATR pspsolve_
if (isRelativistic_Atom ())
{
rpspsolve_ (print_ptr, debug_ptr, lmax_ptr, locp_ptr, rlocal_ptr,
printf("Relativstic Atom detected!\n");
rpspsolve_(print_ptr, debug_ptr, lmax_ptr, locp_ptr, rlocal_ptr,
sdir_name, n9, dir_name, n0, in_filename, n1,
out_filename, n2);
return;
@ -493,7 +494,7 @@ void FATR pspsolve_
printf ("Outputing semicore density gradient: %s\n", full_filename);
fp = fopen (full_filename, "w+");
for (k = 0; k < Ngrid; ++k)
fprintf (fp, "%12.8lf %12.8lf\n", rgrid[k], drho_semicore_Psp ()[k]);
fprintf (fp, "%12.8lf %12.8lf\n", rgrid[k], (drho_semicore_Psp ())[k]);
fclose (fp);
/* output all-electron potential infile.pot */

View file

@ -1,6 +1,6 @@
/* psp.c -
author - Patrick Nichols
Spin-Orbit Pseudopotential Generation
author - Eric Bylaska
$Id: rpsp.c,v 1.2 2007-04-10 23:30:34 d3p708 Exp $
*/
#include <stdio.h>
@ -13,14 +13,12 @@
#include "dft.h"
#include "atom.h"
#include "hamann.h"
#include "rhamann.h"
#include "troullier.h"
#include "vanderbilt.h"
#include "generate_rho_semicore.h"
#include "rpsp.h"
#include "psp.h"
#define False 0
#define True 1
#define Max(x,y) ((x>y) ? x : y)
@ -28,37 +26,42 @@
/* atom structure variables */
static int Nvalence;
static int npsp_states;
static int *n;
static int *l;
static int *spin;
static int lmax;
static double *fill;
static double *rcut;
static double *peak;
static double Zion;
static double Total_E,
E_Hartree, P_Hartree, E_exchange, P_exchange, E_correlation, P_correlation;
static double *eigenvalue;
static double **r_psi;
static double **r_psi_prime;
static double *rho;
static double *rho_semicore;
static double r_semicore;
static double **V_psp;
static const char *comment="Hamann Pseudopotential";
static int aSolver;
static int Nvalence;
static int npsp_states;
static int *n;
static int *l;
static int *spin;
static int lmax;
static double *fill;
static double *rcut;
static double *peak;
static double Zion;
static double Total_E,
E_Hartree, P_Hartree,
E_exchange,
P_exchange,
E_correlation,
P_correlation;
static double *eigenvalue;
static double **r_psi;
static double **r_psi_prime;
static double *rho;
static double *rho_semicore;
static double *drho_semicore;
static double r_semicore;
static double **V_psp;
static char comment[80];
/* extra Vanderbilt parameters */
static double rlocal, clocal;
static int ns[10], indx_il[4][10], indx_ijl[4][4][10];
static double rlocal,clocal;
static int ns[10],indx_il[4][10],indx_ijl[4][4][10];
static double *Vlocal;
static double **r_hard_psi;
static double *D0;
static double *q;
/* solver parameters: this is the Kawai-Weare default */
static int Solver_Type = Hamann;
static int Solver_Type = Hamann;
/********************************
@ -67,207 +70,243 @@ static int Solver_Type = Hamann;
* *
********************************/
void
init_RelPsp (char *filename)
void init_RelPsp(char *filename)
{
int p, p1, p2;
int ltmp;
double rctmp;
char *w;
FILE *fp;
int p,p1,p2;
int ltmp,semicore_type;
double rctmp;
char *w,*tc;
FILE *fp;
/* find the comment */
strcpy(comment,"Hamann pseudopotential");
fp = fopen(filename,"r+");
w = get_word(fp);
while ((w!=NIL) && (strcmp("<comment>",w)!=0))
w = get_word(fp);
if(w!=NIL)
{
w = get_word(fp);
p = 0;
tc = comment;
while ((w!=NIL)&&(strcmp("<end>",w) != 0))
{
p = (strlen(w));
strcpy(tc, w);
for(p1=0;p1<p; ++p1) ++tc;
strcpy(tc, " ");
++tc;
w = get_word(fp);
}
}
fclose(fp);
Solver_Type = Hamann;
/* set lmax */
lmax = lmax_Atom () + 1;
/* set the number psp projectors */
npsp_states = 2 * lmax + 1;
fp = fopen (filename, "r+");
w = get_word (fp);
while ((w != NIL) && (strcmp ("<npsp-states>", w) != 0))
w = get_word (fp);
if (w != NIL)
{
w = get_word (fp);
while ((w != NIL) && (strcmp ("<end>", w) != 0))
{
sscanf (w, "%d", &ltmp);
w = get_word (fp);
npsp_states = ltmp;
}
}
fclose (fp);
/* allocate memory for n,l,fill,rcut,peak, and eigenvalue */
n = (int *) malloc ((npsp_states) * sizeof (int));
l = (int *) malloc ((npsp_states) * sizeof (int));
spin = (int *) malloc ((npsp_states) * sizeof (int));
fill = (double *) malloc ((npsp_states) * sizeof (double));
rcut = (double *) malloc ((npsp_states) * sizeof (double));
peak = (double *) malloc ((npsp_states) * sizeof (double));
eigenvalue = (double *) malloc ((npsp_states) * sizeof (double));
/* set lmax */
lmax = lmax_Atom() + 1;
/* set the number psp projectors */
npsp_states = 2*lmax + 2;
/* allocate memory for n,l,fill,rcut,peak, and eigenvalue */
n = (int *) malloc((npsp_states)*sizeof(int));
l = (int *) malloc((npsp_states)*sizeof(int));
spin = (int *) malloc((npsp_states)*sizeof(int));
fill = (double *) malloc((npsp_states)*sizeof(double));
rcut = (double *) malloc((npsp_states)*sizeof(double));
peak = (double *) malloc((npsp_states)*sizeof(double));
eigenvalue = (double *) malloc((npsp_states)*sizeof(double));
/* allocate memory for r_psi, V_psp, and rho */
r_psi = (double **) malloc((npsp_states)*sizeof(double*));
r_psi_prime = (double **) malloc((npsp_states)*sizeof(double*));
V_psp = (double **) malloc((npsp_states)*sizeof(double*));
for (p=0; p<npsp_states; ++p)
{
r_psi[p] = alloc_LogGrid();
r_psi_prime[p] = alloc_LogGrid();
V_psp[p] = alloc_LogGrid();
}
rho = alloc_LogGrid();
rho_semicore = alloc_LogGrid();
drho_semicore = alloc_LogGrid();
/* allocate memory for r_psi, V_psp, and rho */
r_psi = (double **) malloc ((npsp_states) * sizeof (double *));
r_psi_prime = (double **) malloc ((npsp_states) * sizeof (double *));
V_psp = (double **) malloc ((npsp_states) * sizeof (double *));
for (p = 0; p < npsp_states; ++p)
{
r_psi[p] = alloc_LogGrid ();
r_psi_prime[p] = alloc_LogGrid ();
V_psp[p] = alloc_LogGrid ();
}
rho = alloc_LogGrid ();
rho_semicore = alloc_LogGrid ();
/* get the psp info */
Suggested_Param_RelHamann(&Nvalence,n,l,spin,eigenvalue,fill,rcut);
/* set the number psp projectors */
fp = fopen(filename,"r+");
w = get_word(fp);
while ((w!=NIL) && (strcmp("<npsp-states>",w)!=0))
w = get_word(fp);
/* get the psp info */
Suggested_Param_RelHamann (&Nvalence, n, l, spin, eigenvalue, fill, rcut);
/* set the number psp projectors */
fp = fopen (filename, "r+");
w = get_word (fp);
while ((w != NIL) && (strcmp ("<npsp-states>", w) != 0))
w = get_word (fp);
if (w != NIL)
{
w = get_word (fp);
while ((w != NIL) && (strcmp ("<end>", w) != 0))
{
sscanf (w, "%d", &ltmp);
w = get_word (fp);
Nvalence = ltmp;
}
}
fclose (fp);
if (w!=NIL)
{
w = get_word(fp);
while ((w!=NIL) && (strcmp("<end>",w) != 0))
{
sscanf(w,"%d", &ltmp);
w = get_word(fp);
Nvalence = 2*ltmp+2;
}
}
fclose(fp);
/* get rcut */
fp = fopen (filename, "r+");
w = get_word (fp);
while ((w != NIL) && (strcmp ("<rcut>", w) != 0))
w = get_word (fp);
/* get rcut */
fp = fopen(filename,"r+");
w = get_word(fp);
while ((w!=NIL) && (strcmp("<rcut>",w)!=0))
w = get_word(fp);
if (w != NIL)
{
w = get_word (fp);
while ((w != NIL) && (strcmp ("<end>", w) != 0))
{
sscanf (w, "%d", &ltmp);
w = get_word (fp);
sscanf (w, "%lf", &rctmp);
w = get_word (fp);
rcut[ltmp] = rctmp;
}
}
fclose (fp);
if (w!=NIL)
{
w = get_word(fp);
while ((w!=NIL) && (strcmp("<end>",w) != 0))
{
sscanf(w,"%d", &ltmp);
w = get_word(fp);
sscanf(w,"%lf", &rctmp);
w = get_word(fp);
rcut[2*ltmp] = rctmp;
rcut[2*ltmp+1] = rctmp;
}
}
fclose(fp);
/* get ecut */
fp = fopen (filename, "r+");
w = get_word (fp);
while ((w != NIL) && (strcmp ("<ecut>", w) != 0))
w = get_word (fp);
/* get ecut */
fp = fopen(filename,"r+");
w = get_word(fp);
while ((w!=NIL) && (strcmp("<ecut>",w)!=0))
w = get_word(fp);
if (w != NIL)
{
w = get_word (fp);
while ((w != NIL) && (strcmp ("<end>", w) != 0))
{
sscanf (w, "%d", &ltmp);
w = get_word (fp);
sscanf (w, "%lf", &rctmp);
w = get_word (fp);
eigenvalue[ltmp] = rctmp;
}
}
fclose (fp);
if (w!=NIL)
{
w = get_word(fp);
while ((w!=NIL) && (strcmp("<end>",w) != 0))
{
sscanf(w,"%d", &ltmp);
w = get_word(fp);
sscanf(w,"%lf", &rctmp);
w = get_word(fp);
eigenvalue[2*ltmp] = rctmp;
eigenvalue[2*ltmp+1] = rctmp;
}
}
fclose(fp);
/* get r_semicore - if zero then no core corrections added */
r_semicore = 0.0;
fp = fopen(filename,"r+");
w = get_word(fp);
while ((w!=NIL) && (strcmp("<semicore>",w)!=0))
w = get_word(fp);
if (w!=NIL)
{
w = get_word(fp);
while ((w!=NIL) && (strcmp("<end>",w) != 0))
{
sscanf(w,"%lf", &rctmp);
w = get_word(fp);
r_semicore = rctmp;
}
}
fclose(fp);
/* find the semicore_type */
semicore_type = 2;
fp = fopen(filename,"r+");
w = get_word(fp);
while ((w!=NIL) && (strcmp("<semicore_type>",w)!=0))
w = get_word(fp);
if (w!=NIL)
{
w = get_word(fp);
if (strcmp("quadratic",w)==0) semicore_type = 0;
if (strcmp("louie",w)==0) semicore_type = 1;
if (strcmp("fuchs",w)==0) semicore_type = 2;
}
fclose(fp);
/* get r_semicore - if zero then no core corrections added */
r_semicore = 0.0;
fp = fopen (filename, "r+");
w = get_word (fp);
while ((w != NIL) && (strcmp ("<semicore>", w) != 0))
w = get_word (fp);
if (w != NIL)
{
w = get_word (fp);
while ((w != NIL) && (strcmp ("<end>", w) != 0))
{
sscanf (w, "%lf", &rctmp);
w = get_word (fp);
r_semicore = rctmp;
}
}
fclose (fp);
/* generate non-zero rho_semicore */
if (r_semicore > 0.0)
{
/*
printf("\n\n");
printf("Generating non-zero semicore density\n");
*/
generate_rho_semicore(semicore_type,rho_core_Atom(),r_semicore,rho_semicore);
Derivative_LogGrid(rho_semicore,drho_semicore);
}
/* generate non-zero rho_semicore */
if (r_semicore > 0.0)
{
printf ("\n\n");
printf ("Generating non-zero semicore density\n");
generate_rho_semicore (rho_core_Atom (), r_semicore, rho_semicore);
}
/* define the ion charge */
/*
Zion=0.0;
for (p=Ncore_Atom(); p<(Ncore_Atom()+Nvalence_Atom()); ++p)
Zion += fill_Atom(p);
*/
/* define the ion charge */
Zion = 0.0;
for (p = Ncore_Atom (); p < (Ncore_Atom () + Nvalence_Atom ()); ++p)
Zion += fill_Atom (p);
Zion = Zion_Atom();
for (p=0; p<(Ncore_Atom()); ++p)
Zion -= fill_Atom(p);
} /* init_RelPsp */
} /* init_RelPsp */
/********************************
* *
* solve_RelPsp *
* solve_RelPsp *
* *
********************************/
void
solve_RelPsp ()
void solve_RelPsp()
{
int p, k, Ngrid;
double *rgrid;
int p,k,Ngrid;
double *rgrid;
/* get loggrid variables */
Ngrid = N_LogGrid ();
rgrid = r_LogGrid ();
/* get loggrid variables */
Ngrid = N_LogGrid();
rgrid = r_LogGrid();
solve_RelHamann (Nvalence, n, l, spin, eigenvalue, fill, rcut,
r_psi, r_psi_prime, rho, rho_semicore, V_psp,
solve_RelHamann(Nvalence,n,l,spin,eigenvalue,fill,rcut,
r_psi,r_psi_prime,rho,rho_semicore,V_psp,
&Total_E,
&E_Hartree, &P_Hartree,
&E_exchange, &P_exchange, &E_correlation, &P_correlation);
&E_Hartree,&P_Hartree,
&E_exchange,&P_exchange,
&E_correlation,&P_correlation);
/******************************************************/
/* find the outermost peak of the pseudowavefunctions */
/* find the outermost peak of the pseudowavefunctions */
/******************************************************/
for (p = 0; p < Nvalence; ++p)
{
for (p=0; p<Nvalence; ++p)
{
if (fill[p] != 0.0)
{
k = Ngrid - 2;
while ((r_psi_prime[p][k] * r_psi_prime[p][k + 1] >= 0.0)
&& (k >= 0))
--k;
peak[p] = rgrid[k];
}
{
k = Ngrid-2;
while ((r_psi_prime[p][k]*r_psi_prime[p][k+1] >= 0.0) && (k>=0))
--k;
peak[p] = rgrid[k];
}
else
peak[p] = rgrid[Ngrid - 1];
}
peak[p] = rgrid[Ngrid-1];
}
} /* solve_RelPsp */
} /* solve_RelPsp */
@ -277,81 +316,79 @@ solve_RelPsp ()
* *
********************************/
char *
solver_Name_RelPsp ()
char *solver_Name_RelPsp()
{
char *s;
if (Solver_Type == Hamann)
s = "Hamann";
else if (Solver_Type == Troullier)
s = "Troullier-Martins";
else if (Solver_Type == Vanderbilt)
s = "Vanderbilt";
else
s = "unknown?";
return s;
char *s;
s = "Hamann";
return s;
}
void
print_RelPsp (fp)
FILE *fp;
void print_RelPsp(fp)
FILE *fp;
{
int i;
int i;
fprintf(fp,"\n\n");
fprintf(fp,"PSP solver information\n\n");
fprintf(fp,"Atom name: %s\n",name_Atom());
fprintf(fp,"Zcharge : %le\n",Zion);
fprintf(fp,"Nvalence : %d\n",Nvalence);
fprintf(fp," : restricted calculation\n");
fprintf (fp, "\n\n");
fprintf (fp, "PSP solver information\n\n");
fprintf (fp, "Atom name: %s\n", name_Atom ());
fprintf (fp, "Zcharge : %le\n", Zion);
fprintf (fp, "Nvalence : %d\n", Nvalence);
fprintf (fp, " : restricted calculation\n");
fprintf(fp,"\n------ Solver information ------\n");
fprintf(fp,"solver type : %s\n",solver_Name_RelPsp());
fprintf(fp,"hartree type : %s\n",hartree_Name_DFT());
fprintf(fp,"exchange type : %s\n",exchange_Name_DFT());
if (strcmp(exchange_Name_DFT(),"Dirac") == 0)
fprintf(fp," alpha : %lf\n",Dirac_alpha());
fprintf(fp,"correlation type : %s\n",correlation_Name_DFT());
fprintf (fp, "\n------ Solver information ------\n");
fprintf (fp, "solver type : %s\n", solver_Name_RelPsp ());
fprintf (fp, "hartree type : %s\n", hartree_Name_DFT ());
fprintf (fp, "exchange type : %s\n", exchange_Name_DFT ());
if (strcmp (exchange_Name_DFT (), "Dirac") == 0)
fprintf (fp, " alpha : %lf\n", Dirac_alpha ());
fprintf (fp, "correlation type : %s\n", correlation_Name_DFT ());
fprintf(fp,"----------------------------------------------------------------------------\n");
fprintf(fp,"n\tl\tpopulation\tEcut\t\tRcut\t\tOuter Peak\n");
for (i=0; i<Nvalence; ++i)
{
fprintf(fp,"%d\t%s\t%.1lf\t%.2lf\t\t%le\t%le\t%le\n",l[i]+1,spd_Name(l[i]),
0.5*spin[i],
fill[i],eigenvalue[i],rcut[i],peak[i]);
}
fprintf(fp,"----------------------------------------------------------------------------\n");
if (r_semicore > 0.0)
{
fprintf(fp,"SemiCore Corrections Added\n");
fprintf(fp," rcore : %lf\n",r_semicore);
fprintf(fp," Semicore Charge : %lf\n",
Integrate_LogGrid(rho_semicore));
fprintf(fp," Semicore Charge gradient : %lf\n",
Integrate_LogGrid(drho_semicore));
}
fprintf (fp,
"----------------------------------------------------------------------------\n");
fprintf (fp, "n\tl\tpopulation\tEcut\t\tRcut\t\tOuter Peak\n");
if (Solver_Type==Vanderbilt)
fprintf(fp,"Using AE valence density for descreening\n");
for (i = 0; i < Nvalence; ++i)
{
fprintf (fp, "%d\t%s\t%3.1lf\t%s\t\t%le\t%le\t%le\n", l[i] + 1,
spd_Name (l[i]), 0.5 * spin[i], " -", eigenvalue[i], rcut[i],
peak[i]);
}
fprintf (fp,
"----------------------------------------------------------------------------\n");
if (r_semicore > 0.0)
{
fprintf (fp, "SemiCore Corrections Added\n");
fprintf (fp, " rcore : %lf\n", r_semicore);
fprintf (fp, " Semicore Charge : %lf\n",
Integrate_LogGrid (rho_semicore));
}
fprintf (fp, "Pseudopotential Charge: %lf\n", Integrate_LogGrid (rho));
fprintf (fp, "\nTotal E = %le\n", Total_E);
fprintf (fp, "\n");
fprintf(fp, "Pseudopotential ion charge = %lf\n", Zion);
fprintf(fp, "Pseudopotential electronic charge= %lf\n", -Integrate_LogGrid(rho));
fprintf(fp, "Pseudopotential atom charge = %lf\n", Zion-Integrate_LogGrid(rho));
fprintf(fp,"\nTotal E = %le\n",Total_E);
fprintf(fp,"\n");
fprintf (fp, "E_Hartree = %le\n", E_Hartree);
fprintf (fp, "<Vh> = %le\n", P_Hartree);
fprintf (fp, "\n");
fprintf(fp,"E_Hartree = %le\n",E_Hartree);
fprintf(fp,"<Vh> = %le\n",P_Hartree);
fprintf(fp,"\n");
fprintf (fp, "E_exchange = %le\n", E_exchange);
fprintf (fp, "<Vx> = %le\n", P_exchange);
fprintf (fp, "\n");
fprintf(fp,"E_exchange = %le\n",E_exchange);
fprintf(fp,"<Vx> = %le\n",P_exchange);
fprintf(fp,"\n");
fprintf (fp, "E_correlation = %le\n", E_correlation);
fprintf (fp, "<Vc> = %le\n\n", P_correlation);
} /* print_Atom */
fprintf(fp,"E_correlation = %le\n",E_correlation);
fprintf(fp,"<Vc> = %le\n\n",P_correlation);
} /* print_Atom */
/********************************
* *
@ -359,23 +396,30 @@ print_RelPsp (fp)
* *
********************************/
void
set_Solver_RelPsp (solver)
int solver;
void set_Solver_RelPsp(solver)
int solver;
{
Solver_Type = solver;
fprintf(stdout,"RelRelPsp::Cannot set the solver yet!\n");
}
int
Vanderbilt_RelPsp ()
int Vanderbilt_RelPsp()
{
int value;
int value;
value = 0;
if (Solver_Type == Vanderbilt)
value = 1;
return value;
value = 0;
if (Solver_Type==Vanderbilt) value = 1;
return value;
}
int NormConserving_RelPsp()
{
int value;
value = 0;
if (Solver_Type==Hamann) value = 1;
if (Solver_Type==Troullier) value = 1;
return value;
}
/********************************
@ -384,169 +428,149 @@ Vanderbilt_RelPsp ()
* *
********************************/
double
E_RelPsp ()
double E_RelPsp()
{
return Total_E;
} /*E_Atom */
return Total_E;
} /*E_Atom*/
double
eigenvalue_RelPsp (int i)
double eigenvalue_RelPsp(int i)
{
return eigenvalue[i];
return eigenvalue[i];
}
double *
rho_RelPsp ()
double *rho_RelPsp()
{
return rho;
return rho;
}
double *
rho_semicore_RelPsp ()
double *rho_semicore_RelPsp()
{
return rho_semicore;
return rho_semicore;
}
double
r_semicore_RelPsp ()
double * drho_semicore_RelPsp()
{
return r_semicore;
return drho_semicore;
}
double r_semicore_RelPsp()
{
return r_semicore;
}
double *
Beta_RelPsp (int i, int l)
double *Beta_RelPsp(int i,int l)
{
return V_psp[indx_il[i][l]];
return V_psp[indx_il[i][l]];
}
double *
r_psi_il_RelPsp (int i, int l)
double *r_psi_il_RelPsp(int i, int l)
{
return r_psi[indx_il[i][l]];
return r_psi[indx_il[i][l]];
}
double *
r_hard_psi_il_RelPsp (int i, int l)
double *r_hard_psi_il_RelPsp(int i, int l)
{
return r_hard_psi[indx_il[i][l]];
return r_hard_psi[indx_il[i][l]];
}
int
ns_RelPsp (int l)
int ns_RelPsp(int l)
{
return ns[l];
return ns[l];
}
double
D0_RelPsp (int i, int j, int l)
double D0_RelPsp(int i, int j, int l)
{
return D0[indx_ijl[i][j][l]];
return D0[indx_ijl[i][j][l]];
}
double
q_RelPsp (int i, int j, int l)
double q_RelPsp(int i, int j, int l)
{
return q[indx_ijl[i][j][l]];
return q[indx_ijl[i][j][l]];
}
double *
Vlocal_RelPsp ()
double *Vlocal_RelPsp()
{
return Vlocal;
return Vlocal;
}
double *
V_RelPsp (int i)
double *V_RelPsp(int i)
{
return V_psp[i];
return V_psp[i];
}
double *
r_psi_RelPsp (int i)
double *r_psi_RelPsp(int i)
{
return r_psi[i];
return r_psi[i];
}
int
n_RelPsp (int i)
int n_RelPsp(int i)
{
return n[i];
return n[i];
}
int
l_RelPsp (int i)
int l_RelPsp(int i)
{
return l[i];
return l[i];
}
int
lmax_RelPsp ()
int lmax_RelPsp()
{
return lmax;
return lmax;
}
double
fill_RelPsp (int i)
double fill_RelPsp(int i)
{
return fill[i];
return fill[i];
}
int
Nvalence_RelPsp ()
int Nvalence_RelPsp()
{
return Nvalence;
return Nvalence;
}
double
peak_RelPsp (int i)
double peak_RelPsp(int i)
{
return peak[i];
return peak[i];
}
double
rcut_RelPsp (int i)
double rcut_RelPsp(int i)
{
return rcut[i];
return rcut[i];
}
double
rcut_il_RelPsp (int i, int l)
double rcut_il_RelPsp(int i, int l)
{
return rcut[indx_il[i][l]];
return rcut[ indx_il[i][l] ];
}
double
Zion_RelPsp ()
double Zion_RelPsp()
{
return Zion;
return Zion;
}
int
state_RelPsp (int nt, int lt)
int state_RelPsp(int nt, int lt)
{
int i;
i = 0;
while (((nt != n[i]) || (lt != l[i])) && (i < Nvalence))
++i;
int i;
i = 0;
while ( ((nt != n[i]) || (lt != l[i])) && (i<Nvalence) )
++i;
/* Error */
if (i >= Nvalence)
printf ("Error: state_RelPsp\n");
/* Error */
if (i>=Nvalence)
printf("Error: state_RelPsp\n");
return i;
return i;
}
const char *
comment_RelPsp ()
char *comment_RelPsp()
{
return comment;
return comment;
}

View file

@ -1,39 +1,39 @@
#ifndef _RPSP_H_
#define _RPSP_H_
#include <stdio.h>
/*
* rpsp.c
* author: Patrick Nichols
*/
extern void init_RelPsp (char *filename);
extern void solve_RelPsp (void);
extern char *solver_Name_RelPsp (void);
extern void print_RelPsp (FILE * fp);
extern void set_Solver_RelPsp (int solver);
extern intVanderbilt_RelPsp (void);
extern double E_RelPsp (void);
extern double eigenvalue_RelPsp (int i);
extern double *rho_RelPsp (void);
extern double *rho_semicore_RelPsp (void);
extern double r_semicore_RelPsp (void);
extern double *Beta_RelPsp (int i, int l);
extern double *r_psi_il_RelPsp (int i, int l);
extern double *r_hard_psi_il_RelPsp (int i, int l);
extern int ns_RelPsp (int l);
extern double D0_RelPsp (int i, int j, int l);
extern double q_RelPsp (int i, int j, int l);
extern double *Vlocal_RelPsp (void);
extern double *V_RelPsp (int i);
extern double *r_psi_RelPsp (int i);
extern int n_RelPsp (int i);
extern int l_RelPsp (int i);
extern int lmax_RelPsp (void);
extern double fill_RelPsp (int i);
extern int Nvalence_RelPsp (void);
extern double peak_RelPsp (int i);
extern double rcut_RelPsp (int i);
extern double rcut_il_RelPsp (int i, int l);
extern double Zion_RelPsp (void);
extern int state_RelPsp (int nt, int lt);
extern const char *comment_RelPsp (void);
/* rpsp.c */
extern void init_RelPsp(char *filename);
extern void solve_RelPsp(void);
extern char *solver_Name_RelPsp(void);
extern void print_RelPsp(FILE *fp);
extern void set_Solver_RelPsp(int solver);
extern int Vanderbilt_RelPsp(void);
extern int NormConserving_RelPsp(void);
extern double E_RelPsp(void);
extern double eigenvalue_RelPsp(int i);
extern double *rho_RelPsp(void);
extern double *rho_semicore_RelPsp(void);
extern double *drho_semicore_RelPsp(void);
extern double r_semicore_RelPsp(void);
extern double *Beta_RelPsp(int i, int l);
extern double *r_psi_il_RelPsp(int i, int l);
extern double *r_hard_psi_il_RelPsp(int i, int l);
extern int ns_RelPsp(int l);
extern double D0_RelPsp(int i, int j, int l);
extern double q_RelPsp(int i, int j, int l);
extern double *Vlocal_RelPsp(void);
extern double *V_RelPsp(int i);
extern double *r_psi_RelPsp(int i);
extern int n_RelPsp(int i);
extern int l_RelPsp(int i);
extern int lmax_RelPsp(void);
extern double fill_RelPsp(int i);
extern int Nvalence_RelPsp(void);
extern double peak_RelPsp(int i);
extern double rcut_RelPsp(int i);
extern double rcut_il_RelPsp(int i, int l);
extern double Zion_RelPsp(void);
extern int state_RelPsp(int nt, int lt);
extern char *comment_RelPsp(void);
#endif

View file

@ -106,15 +106,15 @@ void FATR rpspsolve_
/* we have already solved for the Relativsitic AE wavefunctions using atom */
set_debug_print (debug);
init_Psp (infile);
solve_Psp ();
init_RelPsp (infile);
solve_RelPsp ();
if (debug)
print_Psp (stdout);
print_RelPsp (stdout);
init_Linear (infile);
/* allocate linear meshes */
Nvalence = Nvalence_Psp ();
Nvalence = Nvalence_RelPsp ();
Nlinear = nrl_Linear ();
psil = (double **) malloc (Nvalence * sizeof (double *));
pspl = (double **) malloc (Nvalence * sizeof (double *));
@ -129,16 +129,16 @@ void FATR rpspsolve_
/* Norm-conserving output */
for (p = 0; p < Nvalence; ++p)
{
Log_to_Linear (r_psi_Psp (p), rl, psil[p]);
Log_to_Linear_zero (V_Psp (p), rl, pspl[p]);
Log_to_Linear (r_psi_RelPsp (p), rl, psil[p]);
Log_to_Linear_zero (V_RelPsp (p), rl, pspl[p]);
/* normalize scattering state */
if (fill_Psp (p) == 0.0)
if (fill_RelPsp (p) == 0.0)
{
normalize_Linear (psil[p]);
}
}
Log_to_Linear (rho_Psp (), rl, rhol);
Log_to_Linear (rho_RelPsp (), rl, rhol);
if (debug)
{
@ -211,23 +211,23 @@ void FATR rpspsolve_
}
fp = fopen (outfile, "w+");
fprintf (fp, "%s\n", name_Atom ());
fprintf (fp, "%lf %lf %d %d %d %lf\n", Zion_Psp (), Amass_Atom (),
lmax_Psp (), lmax_out, locp_out, rlocal_out);
for (p = 0; p <= lmax_Psp (); ++p)
fprintf (fp, "%lf ", rcut_Psp (p));
fprintf (fp, "%lf %lf %d %d %d %lf\n", Zion_RelPsp (), Amass_Atom (),
lmax_RelPsp (), lmax_out, locp_out, rlocal_out);
for (p = 0; p <= lmax_RelPsp (); ++p)
fprintf (fp, "%lf ", rcut_RelPsp (p));
fprintf (fp, "\n");
fprintf (fp, "%d %lf\n", nrl_Linear (), drl_Linear ());
fprintf (fp, "%s\n", comment_Psp ());
fprintf (fp, "%s\n", comment_RelPsp ());
if (print)
{
printf (" + Appending pseudopotentials: %s thru %s\n",
spd_Name (0), spd_Name (lmax_Psp ()));
spd_Name (0), spd_Name (lmax_RelPsp ()));
}
for (k = 0; k < Nlinear; ++k)
{
fprintf (fp, "%12.8lf", rl[k]);
for (p = 0; p <= lmax_Psp (); ++p)
for (p = 0; p <= lmax_RelPsp (); ++p)
{
vx = p * pspl[2 * p][k] + (p + 1.) * pspl[2 * p + 1][k];
vx *= 1. / (2. * p + 1.);
@ -238,25 +238,25 @@ void FATR rpspsolve_
if (print)
{
printf (" + Appending pseudowavefunctions: %s thru %s\n",
spd_Name (0), spd_Name (lmax_Psp ()));
spd_Name (0), spd_Name (lmax_RelPsp ()));
}
for (k = 0; k < Nlinear; ++k)
{
fprintf (fp, "%12.8lf", rl[k]);
for (p = 0; p <= lmax_Psp (); ++p)
for (p = 0; p <= lmax_RelPsp (); ++p)
fprintf (fp, " %12.8lf", psil[p][k]);
fprintf (fp, "\n");
}
/* append semicore corrections */
if (r_semicore_Psp () != 0.0)
if (r_semicore_RelPsp () != 0.0)
{
if (print)
{
printf (" + Appending semicore density\n");
}
Log_to_Linear (rho_semicore_Psp (), rl, rhol);
fprintf (fp, "%lf\n", r_semicore_Psp ());
Log_to_Linear (rho_semicore_RelPsp (), rl, rhol);
fprintf (fp, "%lf\n", r_semicore_RelPsp ());
for (k = 0; k < Nlinear; ++k)
fprintf (fp, "%12.8lf %12.8lf\n", rl[k],
fabs (rhol[k] * over_fourpi));
@ -265,7 +265,7 @@ void FATR rpspsolve_
{
printf (" + Appending semicore density gradient\n");
}
Log_to_Linear (drho_semicore_Psp (), rl, rhol);
Log_to_Linear (drho_semicore_RelPsp (), rl, rhol);
for (k = 0; k < Nlinear; ++k)
fprintf (fp, "%12.8lf %12.8lf\n", rl[k], (rhol[k] * over_fourpi));
}
@ -283,23 +283,23 @@ void FATR rpspsolve_
strncat (soutfile, ".so", 3);
fp = fopen (soutfile, "w+");
fprintf (fp, "%s\n", name_Atom ());
fprintf (fp, "%lf %lf %d %d %d %lf\n", Zion_Psp (), Amass_Atom (),
lmax_Psp (), lmax_out, locp_out, rlocal_out);
for (p = 0; p <= lmax_Psp (); ++p)
fprintf (fp, "%lf ", rcut_Psp (p));
fprintf (fp, "%lf %lf %d %d %d %lf\n", Zion_RelPsp (), Amass_Atom (),
lmax_RelPsp (), lmax_out, locp_out, rlocal_out);
for (p = 0; p <= lmax_RelPsp (); ++p)
fprintf (fp, "%lf ", rcut_RelPsp (p));
fprintf (fp, "\n");
fprintf (fp, "%d %lf\n", nrl_Linear (), drl_Linear ());
fprintf (fp, "%s\n", comment_Psp ());
fprintf (fp, "%s\n", comment_RelPsp ());
if (print)
{
printf (" + Appending Spin-Orbit pseudopotentials: %s thru %s\n",
spd_Name (0), spd_Name (lmax_Psp ()));
spd_Name (0), spd_Name (lmax_RelPsp ()));
}
for (k = 0; k < Nlinear; ++k)
{
fprintf (fp, "%12.8lf", rl[k]);
for (p = 1; p <= lmax_Psp (); ++p)
for (p = 1; p <= lmax_RelPsp (); ++p)
{
vx = pspl[2 * p][k] - pspl[2 * p + 1][k];
vx *= 2. / (2. * p + 1.);
@ -413,7 +413,7 @@ void FATR rpspsolve_
printf ("Outputing semicore density: %s\n", full_filename);
fp = fopen (full_filename, "w+");
for (k = 0; k < Ngrid; ++k)
fprintf (fp, "%12.8lf %12.8lf\n", rgrid[k], rho_semicore_Psp ()[k]);
fprintf (fp, "%12.8lf %12.8lf\n", rgrid[k], rho_semicore_RelPsp ()[k]);
fclose (fp);
/* output semicore density gradient infile.sddns */
@ -431,7 +431,10 @@ void FATR rpspsolve_
printf ("Outputing semicore density gradient: %s\n", full_filename);
fp = fopen (full_filename, "w+");
for (k = 0; k < Ngrid; ++k)
fprintf (fp, "%12.8lf %12.8lf\n", rgrid[k], drho_semicore_Psp ()[k]);
{
vx=drho_semicore_RelPsp()[k];
fprintf (fp, "%12.8lf %12.8lf\n", rgrid[k], vx);
}
fclose (fp);
/* output all-electron potential infile.pot */
@ -454,7 +457,7 @@ void FATR rpspsolve_
c = 0.0;
for (k = 0; k < Ngrid; ++k)
{
x = rgrid[k] / rcut_Psp (0);
x = rgrid[k] / rcut_RelPsp (0);
x = pow (x, 3.5);
x = exp (-x);
y = 1.0 - x;