/*
DFT++ is a density functional package developed by the research group
of Professor Tomas Arias
Copyright 1996-2003 Sohrab Ismail-Beigi
This file is part of DFT++.
DFT++ is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
DFT++ is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with DFT++; if not, write to the Free Software
Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
Please see the file CREDITS for a list of authors.
For academic users, we request that publications using results obtained with
this software reference
"New algebraic formulation of density functional calculation," by Sohrab Ismail-Beigi
and T.A. Arias, Computer Physics Communications 128:1-2, 1-45 (June 2000).
and, if using the wavelet basis, further reference
"Multiresolution analysis of electronic structure: semicardinal and wavelet bases,"
T.A. Arias, Reviews of Modern Physics 71:1, 267-311 (January 1999).
and
"Robust ab initio calculation of condensed matter: transparent convergence through
semicardinal multiresolution analysis,'' I.P. Daykov, T.A. Arias, and
Torkel D. Engeness, Physical Review Letters, 90:21, 216402 (May 2003).
For your convenience, preprints of the above articles may be obtained from
http://arXiv.org/abs/cond-mat/9909130, 9805262, and 0204411, respectively.
*/
#include "header.h"
PW_IonicPotential::PW_IonicPotential()
{
potential_type = DFT_MIT_PSP;
dEperNatoms_dNGperVol = 0.0;
strcpy(potfilename,"");
strcpy(pulayfilename,"");
ngrid_loc = 0;
dq_loc = 0.0;
V_loc = 0;
nlm = 0;
l = m = ngamma = 0;
M = 0;
ngrid_nl = 0;
dq_nl = 0;
flq = 0;
p_nl = 0;
p_l = p_np = 0;
p_ngrid_nl = 0;
p_dq_nl = 0;
p_flq = 0;
Z = 0;
// for FHI-type PSP
fhi_max_l = fhi_loc_l = 'd'; // default max-l is d, local component is d
fhi_ngrid_loc = 2004; // default Vloc total number of grid points.
fhi_dq_loc = 0.02; // default Vloc grid seperation.
fhi_ngrid_nl = 501; // default flq total number of grid points.
fhi_dq_nl = 0.02; // default flq grid seperation.
}
PW_IonicPotential::~PW_IonicPotential()
{
myfree(V_loc);
for (int lm=0; lm < nlm; lm++)
{
myfree(ngrid_nl[lm]);
myfree(dq_nl[lm]);
for (int g=0; g < ngamma[lm]; g++)
myfree(flq[lm][g]);
myfree(flq[lm]);
M[lm].freemem();
}
if (nlm>0)
{
myfree(l);
myfree(m);
myfree(ngamma);
myfree(M);
myfree(ngrid_nl);
myfree(dq_nl);
myfree(flq);
}
for (int ll=0; ll < p_nl; ll++)
{
myfree(p_ngrid_nl[ll]);
myfree(p_dq_nl[ll]);
for (int p=0; p < p_np[ll]; p++)
myfree(p_flq[ll][p]);
myfree(p_flq[ll]);
}
if (p_nl>0)
{
myfree(p_l);
myfree(p_np);
myfree(p_ngrid_nl);
myfree(p_dq_nl);
myfree(p_flq);
}
}
//
// Read the pseudopotential for this species from file
//
/*
* Read the file <potfilename> to get all the specifications for the
* ionic pseudopotential
*
***************************************************************************
* (3) The format of <potfilename> is:
* [#comment]
* ngrid_loc dq_loc <locvqfile>
* nlm
* # for i = 1 to nlm {
* l_i m_i ngamma_i [# comment]
* M_i
* # for j = 1 to ngamma_i {
* ngrid_ij dq_ij <flqfile_ij> [# comment]
* # }
* # }
* if we have projectors, then after this,
* p_nl
* # for i = 1 to p_nl {
* p_l_i p_np_i
* # for j = 1 to p_np_i {
* p_ngrid_ij p_dq_ij <p_flqfile_ij>
* # }
* # }
*
* -> ngrid_loc is an integer specifying the number of
* grid points specifying the local pseudopot.
* -> dq_loc is a real number specifying the spacing in q-space of the grid.
* -> <locvqfile> contains the actual values of the local pseudopotential.
* -> nlm is an integer specifying the number of (l,m) states for the
* non-local pseudopotentials.
* -> l and m are integers and are the values of l and m.
* -> ngamma is an integer specifying the number of multiple-projector
* states for the current (l,m).
* -> M is a ngamma x ngamma matrix of real numbers specifying the
* details of the multiple-projector states at the current (l,m).
* M is read a as a string of real numbers separated by spaces, with
* the column-index of M changing fastest. The last number specifying
* M MUST BE FOLLOWED BY A CARRIAGE RETURN (i.e must end the line).
* -> ngrid (int) is an integer specifying the number of grid-points defining
* the non-local pseudopotential for the current (l,m,j) state.
* -> dq (real) is the regular q-space spacing for the potential in the
* (l,m,j) states.
* -> <flqfile> contains the actual values of the pseudopotential for (l,m,j)
*
* Note:
* -> [# comment] means any text from here to the end-of-line is ignored
* -> Lines starting with '#' are commenst and are ignored.
*
* -> IMPORTANT NOTE REGARDING M:
* THERE CAN BE NO COMMENTS AT ALL SEPARATING THE VALUES THAT SPECIFY
* THE M-MATRIX. FURTHREMORE, THE VERY LAST VALUE OF M SPECIFIED MUST
* BE FOLLOWED BY A CARRIAGE RETURN (I.E. MUST END THE LINE IT IS ON).
*
***************************************************************************
* (4) The format of the <locvqfile> and <flqfile> files:
*
* 0 v(0)
* dq v(dq)
* 2*dq v(2*dq)
* 3*dq v(3*dq)
* .....
*
* where dq is the spacing specified for the locvq/flq files, and v(q) is
* value of the pseudopotential at q.
* There must be ngrid pairs of numbers (q,v(q)).
*
* NO COMMENTS ALLOWED.
*
*/
void PW_IonicPotential::read_dft_psp()
{
dft_text_FILE *potfile,*vqfile,*flqfile;
char line[DFT_LINE_LEN],
vqfilename[DFT_FILENAME_LEN],
flqfilename[DFT_FILENAME_LEN];
/* Open the potential file for current species */
if ( (potfile = dft_text_fopen(potfilename,"r")) == NULL)
die("Can't open potential file for reading.\n");
/*
* Read gridding information for local pseudopotential and file name
* containing atomic potential in q-space
*/
do { dft_text_fgets(line,DFT_LINE_LEN,potfile); } while(line[0] == '#');
sscanf(line,"%d %lg %s",&ngrid_loc,&dq_loc,vqfilename);
dft_log("Reading potential file %s\n",potfilename);
dft_log("Local pseudopotential: ngrid = %d dq = %g file = %s\n",
ngrid_loc,dq_loc,vqfilename);
dft_log_flush();
/* Allocate space to store local pseudopotential on grid */
V_loc = (real *)mymalloc(sizeof(real)*ngrid_loc,
"V_loc","Speciesinfo::read_dft_psp()");
/* Read the file containg local pseudopotential in q-space: the
* format is ngrid triplets of real values (q,V(q),V'(q)) one after the
* other in ASCII text format. */
if ( (vqfile = dft_text_fopen(vqfilename,"r")) == (dft_text_FILE *)0)
die("Can't open vqfile for reading.\n");
real r; // not used... just to read the file below
for (int n=0; n < ngrid_loc; n++)
dft_text_fscanf(vqfile,"%lg %lg",&r,&(V_loc[n]));
dft_text_fclose(vqfile);
/*
* Read non-local pseudopotential information:
*
* first read nlm on a line, the number of (l,m) states described
*/
do { dft_text_fgets(line,DFT_LINE_LEN,potfile); } while(line[0] == '#');
sscanf(line,"%d",&nlm);
dft_log("Non-local pseudopotential(s): nlm = %d\n", nlm);
dft_log_flush();
/* If there are non-local pseudopots to read, read them! */
if (nlm != 0 ){
/* Allocate space for l,m,ng,M,ngrid_nl,dq_nl,
* and flq/flqprime at each (l,m) */
l = (int *)mymalloc(sizeof(int)*nlm,"l","Speciesinfo::read_dft_psp()");
m = (int *)mymalloc(sizeof(int)*nlm,"m","Speciesinfo::read_dft_psp()");
ngamma = (int *)mymalloc(sizeof(int)*nlm,
"ngamma","Speciesinfo::read_dft_psp()");
M = (Matrix *)mymalloc(sizeof(Matrix)*nlm,
"M","Speciesinfo::read_dft_psp()");
ngrid_nl = (int **)mymalloc(sizeof(int *)*nlm,
"ngrid_nl","Speciesinfo::read_dft_psp()");
dq_nl = (real **)mymalloc(sizeof(real *)*nlm,
"dq_nl","Speciesinfo::read_dft_psp()");
flq = (real ***)mymalloc(sizeof(real **)*nlm,
"flq","Speciesinfo::read_dft_psp()");
/* Loop over the (l,m) states */
for (int lm = 0; lm < nlm; lm++)
{
/* Read l, m , and ngamma for this (l,m) state */
do { dft_text_fgets(line,DFT_LINE_LEN,potfile); } while(line[0] == '#');
sscanf(line,"%d %d %d",&(l[lm]),&(m[lm]),&(ngamma[lm]));
int ng = ngamma[lm];
dft_log("lm = %d:\n l = %d m = %d ngamma = %d\n",
lm,l[lm],m[lm],ng);
dft_log_flush();
/* Allocate space for ngrid_nl[lm], dq[lm], and flq/flqprime[lm]
* at each gamma */
ngrid_nl[lm] = (int *)mymalloc(sizeof(int)*ng,
"ngrid_nl[lm]",
"Speciesinfo::read_dft_psp()");
dq_nl[lm] = (real *)mymalloc(sizeof(real)*ng,
"dq_nl[lm]",
"Speciesinfo::read_dft_psp()");
flq[lm] = (real **)mymalloc(sizeof(real *)*ng,
"flq[lm]","Speciesinfo::read_dft_psp()");
/* Allocate and read the ng x ng M[lm] matrix */
M[lm].init(ng,ng);
int gamma;
for (gamma = 0; gamma < ng; gamma++)
for (int gammap = 0; gammap < ng; gammap++)
{
dft_text_fscanf(potfile,"%lg",&r);
M[lm](gamma,gammap) = (scalar)r;
}
dft_log("M = \n");
for (gamma = 0; gamma < ng; gamma++)
{
dft_log("[ ");
for (int gammap = 0; gammap < ng; gammap++)
dft_log("%g ",REAL(M[lm](gamma,gammap)));
dft_log(" ]\n");
}
dft_log_flush();
/* Skip to next line after last M-entry */
dft_text_fgets(line,DFT_LINE_LEN,potfile);
/* For each gamma, read the gridding info and flq filename */
for (gamma = 0; gamma < ng; gamma++)
{
/* Read gridding information and flq file name for
* the current [lm][gamma] pseudopotential */
do
{ dft_text_fgets(line,DFT_LINE_LEN,potfile); }
while(line[0] == '#');
sscanf(line,"%d %lg %s",
&(ngrid_nl[lm][gamma]),
&(dq_nl[lm][gamma]),
flqfilename);
int ng_nl = ngrid_nl[lm][gamma];
dft_log(" gamma = %d: ngrid = %d dq = %g file = %s\n",
gamma,ng_nl,dq_nl[lm][gamma],flqfilename);
dft_log_flush();
/* Allocate space to store the gridded flq[lm][gamma] */
flq[lm][gamma] = (real *)mymalloc(sizeof(real)*ng_nl,
"flq[lm][gamma]",
"Speciesinfo::read_dft_psp()");
/* Open flq file for flq[lm][gamma] and read it:
* the format is ngrid doublets (q,flq(q)) */
if ( (flqfile = dft_text_fopen(flqfilename,"r")) ==
(dft_text_FILE *)0)
die("Can't open flq file for reading.\n");
for (int n=0; n < ng_nl; n++)
dft_text_fscanf(flqfile,"%lg %lg",&r,&(flq[lm][gamma][n]));
dft_text_fclose(flqfile);
}
} /* of lm: loop on (l,m) states */
}
/* if there are projectors to read, read them! */
if( projectorp == 1){
/*
* Read projector information:
*
* first read nl on a line, the number of l states described
*/
do { dft_text_fgets(line,DFT_LINE_LEN,potfile); } while(line[0] == '#');
sscanf(line,"%d",&p_nl);
dft_log("Projector(s): p_nl = %d\n", p_nl);
dft_log_flush();
/* If there are no projectors to read, we're done! */
if (p_nl==0)
return;
/* Allocate space for p_l,p_np,p_ngrid_nl,p_dq_nl,
* and flq/flqprime at each l */
p_l = (int *)mymalloc(sizeof(int)*p_nl,"p_l","Speciesinfo::read_dft_psp()");
p_np = (int *)mymalloc(sizeof(int)*p_nl,
"p_np","Speciesinfo::read_dft_psp()");
p_ngrid_nl = (int **)mymalloc(sizeof(int *)*p_nl,
"p_ngrid_nl","Speciesinfo::read_dft_psp()");
p_dq_nl = (real **)mymalloc(sizeof(real *)*p_nl,
"p_dq_nl","Speciesinfo::read_dft_psp()");
p_flq = (real ***)mymalloc(sizeof(real **)*p_nl,
"p_flq","Speciesinfo::read_dft_psp()");
/* Loop over the l states */
for (int ll = 0; ll < p_nl; ll++){
/* Read ll , and np for this ll state */
do { dft_text_fgets(line,DFT_LINE_LEN,potfile); } while(line[0] == '#');
sscanf(line,"%d %d",&(p_l[ll]),&(p_np[ll]));
int np = p_np[ll];
dft_log("\n l = %d p_np = %d\n", p_l[ll],np);
dft_log_flush();
/* Allocate space for p_ngrid_nl[ll], p_dq[ll], and p_flq[ll]
* at each np */
p_ngrid_nl[ll] = (int *)mymalloc(sizeof(int)*np,
"p_ngrid_nl[ll]",
"Speciesinfo::read_dft_psp()");
p_dq_nl[ll] = (real *)mymalloc(sizeof(real)*np,
"p_dq_nl[ll]",
"Speciesinfo::read_dft_psp()");
p_flq[ll] = (real **)mymalloc(sizeof(real *)*np,
"p_flq[ll]","Speciesinfo::read_dft_psp()");
/* For each p, read the gridding info and flq filename */
for (int p = 0; p < np; p++){
/* Read gridding information and p_flq file name for
* the current [ll][p] projector */
do
{ dft_text_fgets(line,DFT_LINE_LEN,potfile); }
while(line[0] == '#');
sscanf(line,"%d %lg %s",
&(p_ngrid_nl[ll][p]),
&(p_dq_nl[ll][p]),
flqfilename);
int p_ng_nl = p_ngrid_nl[ll][p];
dft_log(" p = %d: p_ngrid = %d p_dq = %g p_file = %s\n",
p,p_ng_nl,p_dq_nl[ll][p],flqfilename);
dft_log_flush();
/* Allocate space to store the gridded flq[lm][gamma] */
p_flq[ll][p] = (real *)mymalloc(sizeof(real)*p_ng_nl,
"p_flq[ll][p]",
"Speciesinfo::read_dft_psp()");
/* Open p_flq file for p_flq[ll][p] and read it:
* the format is ngrid doublets (q,flq(q)) */
if ( (flqfile = dft_text_fopen(flqfilename,"r")) ==
(dft_text_FILE *)0)
die("Can't open p_flq file for reading.\n");
for (int n=0; n < p_ng_nl; n++)
dft_text_fscanf(flqfile,"%lg %lg",&r,&(p_flq[ll][p][n]));
dft_text_fclose(flqfile);
}
} /* of ll: loop on l states */
}
}
//
// Read pulay stuff from file (if needed)
//
/***************************************************************************
*
* (5) The format of the <pulayfile> is:
*
* [ # comment ]
* NEcuts [ # comment ]
* Ecut1 dEperNatoms_dNGperVol1 [ # comment]
* Ecut2 dEperNatoms_dNGperVol2 [ # comment]
* ... ...
*
* -> NEcuts: the number of pairs of (Ecut,dE...perVol) provided
* -> Ecut: planwave cutoff energy in Hartrees
* -> dEperNatoms_dNGperVol: exactly what is says (NG=# of planewaves)!
*
***************************************************************************
*/
void PW_IonicPotential::setup_pulay(Everything &e)
{
real einfoEcut = e.basis_spec.Ecut;
/* If Pulay correction filename is "none", then set the
* Pulay correction to zero and skip reading the file
*/
if ( strcmp(pulayfilename,"none") == 0 ) {
dft_log("Pulay file is specified as 'none':\n");
dft_log("setting Pulay correction to zero.\n");
dEperNatoms_dNGperVol = 0.0;
} else {
dft_text_FILE *pulay_file;
/* Read file of Pulay correction factors for the species
*/
if ( (pulay_file = dft_text_fopen(pulayfilename,"r")) == 0)
die("Can't open pulay file for reading.\n");
dft_log("Reading pulay file %s\n",pulayfilename);
/*
* Read number of energy cutoff values listed in file
*/
int NEcuts,n,flag;
real Ecut;
char line[DFT_LINE_LEN];
do
{ dft_text_fgets(line,DFT_LINE_LEN,pulay_file); }
while (line[0] == '#');
sscanf(line,"%d",&NEcuts);
dft_log("NEcuts=%d\n",NEcuts);
flag = 0;
for (n=0; n < NEcuts; n++) {
do
{ dft_text_fgets(line,DFT_LINE_LEN,pulay_file); }
while (line[0] == '#');
sscanf(line,"%lg %lg",
&Ecut,&dEperNatoms_dNGperVol);
dft_log("Ecut = %g dEperNatoms_dNGperVol = %g\n",
Ecut,dEperNatoms_dNGperVol);
if (Ecut == einfoEcut) {
dft_log("Using Ecut=%g as it agrees with Elecinfo Ecut\n",Ecut);
flag = 1;
break;
}
}
if (flag==0)
die("\nCan't find appropriate energy cutoff in Pulay file!!!\n\n");
}
}
/*
* simps integration for f with step length h for n points
* n should be odd.
*
* Used by read_fhi_psp() below.
*/
real simps(int n, real* f, real h)
{
real result;
int nm12 = (n-1)/2, i;
if ( (nm12*2) != (n-1) )
die("Function simps: n must be odd, but is not.");
result = f[0];
for (i = 1; i < (n-1); i += 2)
result += 4 * f[i] + 2 * f[i+1];
result -= f[n-1];
result *= (h/3);
return(result);
}
/*
* Read Pseudo potential information from file
*
* 2. The berlin format from package FHI98PPMD
*
* URL: http://www.fhi-berlin.mpg.de/th/fhi98md/fhi98PP/index.html
*
* max_l indicates which angular momentum components are used
* max_l = 0, s
* 1, s, p
* 2, s, p, d
* loc_l indicates which one is being used as the
* local component.
*
* The file format is:
* >>>>>>>( Start )>>>>>>>>>>
* Zv L_max
* (10 lines of number for fitted psp.
* not used here.)
* ngrid clog # (for each angular momentum)
* i r(i) Psi(i) Vion(i)
* ...
* ...
* (ngrid of core electron density for
* nonlinear core correction
* --not yet supported)
* <<<<<<<( End )<<<<<<<<<<
*/
void PW_IonicPotential::read_fhi_psp()
{
dft_text_FILE *potfile;
char line[DFT_LINE_LEN];
int max_l = 2, loc_l = 2;
real zv, clog;
int n, i, nl, lm, gamma, gammap;
int ll, mm, nngamma;
// Temporary variables.
// maximum angular momenta is 3.
real *r[3], *rPsi[3], *Vion[3], fhi_clog[3];
int fhi_ngrid[3];
real *r1, *rPsi1, *V1, *V0, *fun, q, q1, arg;
int ngrid_r = 0;
const int N_Lines_to_Skip = 10;
const real fourpi = 4 * M_PI;
// only support for -1 < loc_l <= max_l < 3
switch (fhi_max_l) {
case 's' : max_l = 0; break;
case 'p' : max_l = 1; break;
case 'd' : max_l = 2; break;
default:
die("Angular momenta other than s, p, and d (0,1,2) are not supported: %d",
max_l);
}
switch (fhi_loc_l) {
case 's' : loc_l = 0; break;
case 'p' : loc_l = 1; break;
case 'd' : loc_l = 2; break;
default:
die("Angular momenta other than s, p, and d (0,1,2) are not supported: %d",
loc_l);
}
if (loc_l > max_l)
die("Illegal local component: max_l = %d, loc_l = %d",max_l,loc_l);
/* Open the potential file for current species */
if ( (potfile = dft_text_fopen(potfilename,"r")) == (dft_text_FILE *)0)
die("Can't open potential file for reading.\n");
ngrid_loc = fhi_ngrid_loc;
dq_loc = fhi_dq_loc;
dft_log("Reading potential file %s (FHI CPI format)\n", potfilename);
dft_log("Local pseudopotential: ngrid = %d dq = %g file = %s",
ngrid_loc, dq_loc, potfilename);
switch (loc_l) {
case 0 : dft_log(" (s) \n"); break;
case 1 : dft_log(" (p) \n"); break;
case 2 : dft_log(" (d) \n"); break;
default: die("Illegal loc_l value!"); break;
}
/* Allocate space to store local pseudopotential on grid */
V_loc = (real *)mymalloc(sizeof(real)*ngrid_loc,
"V_loc","setup_ioninfo");
/*
* The first line has two numbers:
* Zv number_of_angular_momentum_components
*/
dft_text_fgets(line,DFT_LINE_LEN,potfile);
sscanf(line,"%lg %d", &zv, &nl);
// verify number of valence charges, and number of angular momentum
if (zv != Z) {
dft_log("Valence charge unmatched in PSP file %s : %5.1f != %5.1f !!",
potfilename,zv,Z);
die("Fatal error in read_fhi_psp()!");
}
if (nl != (max_l+1)) {
/*
* In the FHI pseudopotentials, nl actually takes values: 1, 2, 3
* (perhaps due to the fact that FORTRAN index starts from 1 instead of 0)
* corresponding to our max_l = 0, 1, 2.
*/
dft_log("Number of angular momentum mismatch: %d != %d !!", nl, max_l+1 );
die("Fatal error in read_fhi_psp()!");
}
// there are 10 lines of data for pseudopotentials that
// are fitted to expressions. Skip.
for (n = 0; n < N_Lines_to_Skip; n++)
dft_text_fgets(line,DFT_LINE_LEN,potfile);
// only allow up to 3 angular momentum components.
// for (i = 0; i < 3; i++)
// sp_info->cpi_ngrid[i] = 0;
/*
* Read in all wavefunctions and potentials.
* All are in real space with logrithmic grid spacing.
*/
for (ll = 0; ll <= max_l; ll++) {
dft_text_fgets(line,DFT_LINE_LEN,potfile);
sscanf(line, "%d %lg", &fhi_ngrid[ll], &clog);
// clog[l] == R[l][1]/R[l][0] is the scale of the log grid spacing.
fhi_clog[ll] = log(clog);
// now clog[l] = D log(R), used for integration later.
dft_log(" l = %d: ngrid = %d, clog = %f\n",
ll, fhi_ngrid[ll], fhi_clog[ll]);
r[ll] = (real*)mymalloc(sizeof(real)*fhi_ngrid[ll],"r1","read_fhi_psp");
rPsi[ll] = (real*)mymalloc(sizeof(real)*fhi_ngrid[ll],"rPsi","read_fhi_psp");
Vion[ll] = (real*)mymalloc(sizeof(real)*fhi_ngrid[ll],"Vion","read_fhi_psp");
if (fhi_ngrid[ll] > ngrid_r)
ngrid_r = fhi_ngrid[ll];
for (i = 0; i < fhi_ngrid[ll]; i++) {
dft_text_fgets(line,DFT_LINE_LEN,potfile);
sscanf(line,"%*d %lg %lg %lg",
&(r[ll][i]),
&(rPsi[ll][i]),
&(Vion[ll][i]));
}
}
/*
* Read core electron density, for nonlinear core correction.
* (not yet supported)
*/
/*
* Convert the local potential from log-grid real space to equal-grid k-space
*
*/
r1 = r[loc_l];
rPsi1 = rPsi[loc_l];
V1 = Vion[loc_l];
fun = (real*)mymalloc(sizeof(real)* ngrid_r,"fun","read_fhi_psp");
for (i = 0; i < fhi_ngrid[loc_l]; i++)
fun[i] = ( V1[i]*r1[i]*r1[i] + Z*r1[i] ) * r1[i];
V_loc[0] = simps(fhi_ngrid[loc_l],fun,fhi_clog[loc_l]) * fourpi;
for (n = 1; n < ngrid_loc; n++) {
q = n * dq_loc;
for (i = 0; i < fhi_ngrid[loc_l]; i++) {
arg = r1[i] * q;
fun[i] = ( V1[i]*r1[i]*r1[i] + Z*r1[i])*r1[i]*sin(arg)/arg;
}
// do the logarithmic-grid integral.
V_loc[n] = simps(fhi_ngrid[loc_l],fun,fhi_clog[loc_l]) * fourpi;
}
/*
* Initialize nlm, l, m, ngamma, and M
*/
nlm = 0;
for (ll=0; ll<=max_l; ll++) {
if (ll != loc_l)
nlm += (2*ll + 1);
}
dft_log("\nNon-local pseudopotential(s): nlm = %d\n", nlm);
/* If there are no non-local pseudopots to read, we're done! */
if (nlm==0)
return;
l = (int *)mymalloc(sizeof(int)*nlm,"l","setup_ioninfo");
m = (int *)mymalloc(sizeof(int)*nlm,"m","setup_ioninfo");
ngamma = (int *)mymalloc(sizeof(int)*nlm,"ngamma","setup_ioninfo");
M = (Matrix *)mymalloc(sizeof(Matrix)*nlm,"M","setup_ioninfo");
ngrid_nl = (int **)mymalloc(sizeof(int *)*nlm,
"ngrid_nl","setup_ioninfo");
dq_nl = (real **)mymalloc(sizeof(real *)*nlm,
"dq_nl","setup_ioninfo");
flq = (real ***)mymalloc(sizeof(real **)*nlm,
"flq","setup_ioninfo");
lm = 0; // set lm counter to 0.
for (ll=0; ll<= max_l; ll++) { // go through all l values.
if (ll != loc_l) { // for nonlocal cases
// Calculate Mnl = ( <phi_l| dV_l |phi_l> )^(-1)
ngrid_r = fhi_ngrid[ll];
r1 = &(r[ll][0]);
V1 = &(Vion[ll][0]);
V0 = &(Vion[loc_l][0]);
rPsi1 = &(rPsi[ll][0]);
for (i = 0; i < ngrid_r; i++)
fun[i] = (V1[i] - V0[i])*r1[i]*(rPsi1[i]*rPsi1[i]);
q = simps(ngrid_r,fun,fhi_clog[ll]);
// for all m values of the current l value
for (mm= ll; mm > (-ll-1); mm--) {
l[lm] = ll;
m[lm] = mm;
nngamma = ngamma[lm] = 1; // single projector for this type of pseudopotential
M[lm].init(nngamma,nngamma);
for (gamma = 0; gamma < nngamma; gamma++)
for (gammap = 0; gammap < nngamma; gammap++)
M[lm](gamma,gammap) = 1.0/q;
dft_log("lm = %d:\n l = %d m = %d ngamma = %d\n",
lm,ll,mm,nngamma);
dft_log("M = \n");
for (gamma = 0; gamma < nngamma; gamma++) {
dft_log("[ ");
for (gammap = 0; gammap < nngamma; gammap++)
dft_log("%g ", REAL(M[lm](gamma,gammap)));
dft_log(" ]\n");
}
/* Allocate space for ngrid_nl[lm], dq[lm], and flq/flqprime[lm]
* at each gamma */
ngrid_nl[lm] = (int *)mymalloc(sizeof(int)*nngamma,
"ngrid_nl[lm]","setup_ioninfo");
dq_nl[lm] = (real *)mymalloc(sizeof(real)*nngamma,
"dq_nl[lm]","setup_ioninfo");
flq[lm] = (real **)mymalloc(sizeof(real *)*nngamma,
"flq[lm]","setup_ioninfo");
// ngamma must be 1 in this case.
ngrid_nl[lm][0] = fhi_ngrid_nl;
dq_nl[lm][0] = fhi_dq_nl;
flq[lm][0] = (real *)mymalloc(sizeof(real)*fhi_ngrid_nl,
"flq[lm][gamma]","setup_ioninfo");
lm++; // increment counter.
}
for (n = 0; n < fhi_ngrid_nl; n++) {
q = n * fhi_dq_nl;
// construct flq from FHI pseudopotential data.
// integrate from pseudoatom wavefunctions.
// int_d(log(r)) r^2 * dV_l(r) * u(r) * j_l(qr)
for (i = 0; i < ngrid_r; i++) {
fun[i] = (V1[i] - V0[i])*(r1[i]*r1[i])*rPsi1[i];
arg = r1[i] * q;
/* The special cases for evaluating j_l(x) near x = 0
* need to be taken care of.
* What cut off value to use and how exactly to do this
* need further thoughts.
*/
q1 = 0;
switch (ll) {
case 0 : // j_0(qr) = sin(qr)/qr
if (arg < 1e-10) // if (qr < 1e-10), j_0(qr) == 1
q1 = 1;
else
q1 = (sin(arg)/arg);
break;
case 1 : // j_1(qr) = sin(qr)/(qr)^2 - cos(qr)/qr
if (arg < 1e-10) // j_1(x) = x/3 - x^3/30 + ...
q1 = arg/3;
else
q1 = (sin(arg)/arg - cos(arg))/arg;
break;
case 2 : // j_2(qr) = (3/(qr)^3-1/qr)sin(qr) - 3cos(qr)/(qr)^2
if (arg < 1e-10) // j_2(x) = x^2/15 - x^4/210 + ...
q1 = arg*arg /15;
else {
q1 = (3/(arg*arg)-1)*sin(arg)-3*cos(arg)/arg;
q1 = q1/arg;
}
break;
default :
die("Illegal l value!");
break;
}
fun[i] *= q1;
}
q1 = simps(ngrid_r,fun,fhi_clog[ll]);
for (mm = 0; mm < (2*ll+1); mm++)
flq[lm-mm-1][0][n] = q1;
}
}
}
// now free up temporary memory.
myfree(fun);
for (ll = 0; ll <= max_l; ll++) {
myfree(r[ll]);
myfree(rPsi[ll]);
myfree(Vion[ll]);
}
}
syntax highlighted by Code2HTML, v. 0.9.1