Skip to content

qbox-public
qbox-public
Commit 9459eedf by Francois Gygi
Options

Merge branch 'develop'

parents c1149965 1d854854
Showing with 231 additions and 119 deletions
src/ChargeDensity.C
......@@ -92,8 +92,6 @@ wf_(wf), vcomm_(wf.sd(0,0)->basis().comm())
ft_[ikp] = new FourierTransform(wb,np0v,np1v,np2v);
}
}
// initialize core density ptr to null ptr
rhocore_r = 0;
}
////////////////////////////////////////////////////////////////////////////////
......@@ -171,11 +169,6 @@ void ChargeDensity::update_density(void)
vft_->forward(&rhotmp[0],&rhog[ispin][0]);
tmap["charge_vft"].stop();
// add core correction charge
if ( rhocore_r )
for ( int i = 0; i < rhor_size; i++ )
rhor[ispin][i] += rhocore_r[i];
}
}
......@@ -197,18 +190,9 @@ void ChargeDensity::update_rhor(void)
const int rhor_size = rhor[ispin].size();
double *const prhor = &rhor[ispin][0];
if ( rhocore_r )
{
#pragma omp parallel for
for ( int i = 0; i < rhor_size; i++ )
prhor[i] = ( rhotmp[i].real() + rhocore_r[i] ) * omega_inv;
}
else
{
#pragma omp parallel for
for ( int i = 0; i < rhor_size; i++ )
prhor[i] = rhotmp[i].real() * omega_inv;
}
#pragma omp parallel for
for ( int i = 0; i < rhor_size; i++ )
prhor[i] = rhotmp[i].real() * omega_inv;
// integral of the charge density
tmap["charge_integral"].start();
......@@ -264,18 +248,9 @@ void ChargeDensity::update_rhog(void)
{
const int rhor_size = rhor[ispin].size();
double *const prhor = &rhor[ispin][0];
if ( rhocore_r )
{
#pragma omp parallel for
for ( int i = 0; i < rhor_size; i++ )
rhotmp[i] = complex<double> ( omega * (prhor[i] - rhocore_r[i]), 0);
}
else
{
#pragma omp parallel for
for ( int i = 0; i < rhor_size; i++ )
rhotmp[i] = complex<double> ( omega * prhor[i], 0);
}
#pragma omp parallel for
for ( int i = 0; i < rhor_size; i++ )
rhotmp[i] = complex<double> ( omega * prhor[i], 0);
assert(rhotmp.size() == vft_->np012loc() );
......
src/ChargeDensity.h
......@@ -51,10 +51,12 @@ class ChargeDensity
mutable TimerMap tmap;
// rhor, rhog: valence density
std::vector<std::vector<double> > rhor; // rhor[ispin][i]
std::vector<std::vector<std::complex<double> > > rhog; // rhog[ispin][ig]
// core density ptr. If non-zero, contains the real-space core density
double* rhocore_r;
// rhocore_r, rhocore_g: core density. Non-zero size only if nlcc used
std::vector<double> rhocore_r;
std::vector<std::complex<double> > rhocore_g;
void update_density(void);
void update_rhor(void);
void update_rhog(void);
......
src/EnergyFunctional.C
......@@ -106,8 +106,6 @@ EnergyFunctional::EnergyFunctional(Sample& s, ChargeDensity& cd)
rhops[is].resize(ngloc);
}
xco = new XCOperator(s_, cd_);
nlp.resize(wf.nspin());
for ( int ispin = 0; ispin < wf.nspin(); ispin++ )
{
......@@ -158,10 +156,8 @@ EnergyFunctional::EnergyFunctional(Sample& s, ChargeDensity& cd)
rhocore_sp_g.resize(nsp_);
for ( int is = 0; is < nsp_; is++ )
rhocore_sp_g[is].resize(ngloc);
rhocore_g.resize(ngloc);
rhocore_r.resize(vft->np012loc());
// set rhocore_r ptr in ChargeDensity
cd_.rhocore_r = &rhocore_r[0];
cd_.rhocore_g.resize(ngloc);
cd_.rhocore_r.resize(vft->np012loc());
}
// FT for interpolation of wavefunctions on the fine grid
......@@ -189,11 +185,9 @@ EnergyFunctional::EnergyFunctional(Sample& s, ChargeDensity& cd)
el_enth_ = new ElectricEnthalpy(s_);
sf.init(tau0,*vbasis_);
xco = new XCOperator(s_, cd_);
cell_moved();
atoms_moved();
}
////////////////////////////////////////////////////////////////////////////////
......@@ -288,6 +282,19 @@ void EnergyFunctional::update_vhxc(bool compute_stress, bool update_vh,
vxc_g[ig] = 0.5 * ( vxc_g[ig] + vtemp[ig] );
}
}
// correct dxc_ to include the core charge exchange-correlation term
// in Harris-Foulkes estimates
// dxc_correction = sum_ig vxc_g[ig] * rhocore_g[ig]
double sum = 0.0;
complex<double> *rh = &cd_.rhocore_g[0];
int len=2*ngloc,inc1=1;
sum = 2.0 * ddot(&len,(double*)&vxc_g[0],&inc1,(double*)&rh[0],&inc1);
// remove double counting for G=0
if ( vbasis_->mype() == 0 )
sum -= real(conj(vxc_g[0])*rh[0]);
double tsum = 0.0;
MPI_Allreduce(&sum,&tsum,1,MPI_DOUBLE,MPI_SUM,vbasis_->comm());
dxc_ += tsum;
}
tmap["exc"].stop();
......@@ -342,7 +349,7 @@ void EnergyFunctional::update_vhxc(bool compute_stress, bool update_vh,
vlocal_g[ig] = vion_local_g[ig] + fpi * rhogt[ig] * g2i[ig];
}
// FT to tmpr_r
// FT to tmp_r
vft->backward(&vlocal_g[0],&tmp_r[0]);
// add external potential vext to tmp_r
......@@ -627,7 +634,7 @@ double EnergyFunctional::energy(bool compute_hpsi, Wavefunction& dwf,
s->drhocoreg(g,rhoc_g,drhoc_g);
// next line: keep only real part
// contribution of core correction
temp -= (vxc_g[ig].real() * sg.real() + vxc_g[ig].imag() * sg.imag())
temp -= (vxc_g[ig].real()*sg.real() + vxc_g[ig].imag()*sg.imag())
* drhoc_g * gi * omega_inv / fpi;
}
......@@ -992,6 +999,7 @@ void EnergyFunctional::atoms_moved(void)
{
const AtomSet& atoms = s_.atoms;
const Wavefunction& wf = s_.wf;
const UnitCell& cell = s_.wf.cell();
int ngloc = vbasis_->localsize();
// fill tau0 with values in atom_list
......@@ -1007,23 +1015,22 @@ void EnergyFunctional::atoms_moved(void)
if ( core_charge_ )
{
// recalculate Fourier coefficients of the core charge
assert(rhocore_g.size()==ngloc);
memset( (void*)&rhocore_g[0], 0, 2*ngloc*sizeof(double) );
const double spin_fac = wf.cell().volume() / wf.nspin();
assert(cd_.rhocore_g.size()==ngloc);
memset( (void*)&cd_.rhocore_g[0], 0, 2*ngloc*sizeof(double) );
// divide core charge by two if two spins
const double spin_fac = cell.volume() / wf.nspin();
for ( int is = 0; is < atoms.nsp(); is++ )
{
complex<double> *s = &sf.sfac[is][0];
double *r = &rhocore_sp_g[is][0];
for ( int ig = 0; ig < ngloc; ig++ )
{
const complex<double> sg = s[ig];
// divide core charge by two if two spins
rhocore_g[ig] += spin_fac * sg * rhocore_sp_g[is][ig];
}
cd_.rhocore_g[ig] += spin_fac * s[ig] * r[ig];
}
// recalculate real-space core charge
vft->backward(&rhocore_g[0],&tmp_r[0]);
vft->backward(&cd_.rhocore_g[0],&tmp_r[0]);
const double omega_inv = 1.0 / cell.volume();
for ( int i = 0; i < tmp_r.size(); i++ )
rhocore_r[i] = tmp_r[i].real();
cd_.rhocore_r[i] = tmp_r[i].real() * omega_inv;
}
for ( int is = 0; is < atoms.nsp(); is++ )
......@@ -1039,7 +1046,6 @@ void EnergyFunctional::atoms_moved(void)
}
// compute esr: pseudocharge repulsion energy
const UnitCell& cell = s_.wf.cell();
const double omega_inv = 1.0 / cell.volume();
// distance between opposite planes of the cell
const double d0 = 2.0 * M_PI / length(cell.b(0));
......
src/NonLocalPotential.C
......@@ -241,6 +241,8 @@ void NonLocalPotential::update_twnl(void)
tmap["update_twnl"].start();
if ( nspnl == 0 ) return;
const int ngwl = basis_.localsize();
const double pi = M_PI;
const double fpi = 4.0 * pi;
......@@ -1257,6 +1259,15 @@ double NonLocalPotential::energy(bool compute_hpsi, SlaterDet& dsd,
bool compute_forces, vector<vector<double> >& fion_enl,
bool compute_stress, valarray<double>& sigma_enl)
{
for ( int is = 0; is < fion_enl.size(); is++ )
for ( int i = 0; i < fion_enl[is].size(); i++ )
fion_enl[is][i] = 0.0;
sigma_enl = 0.0;
if ( nspnl == 0 ) return 0.0;
double enl = 0.0;
double tsum[6] = { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 };
const vector<double>& occ = sd_.occ();
const int ngwl = basis_.localsize();
// define atom block size
......@@ -1267,13 +1278,6 @@ double NonLocalPotential::energy(bool compute_hpsi, SlaterDet& dsd,
vector<vector<double> > tau;
atoms_.get_positions(tau);
double enl = 0.0;
double tsum[6] = { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 };
for ( int is = 0; is < fion_enl.size(); is++ )
for ( int i = 0; i < fion_enl[is].size(); i++ )
fion_enl[is][i] = 0.0;
if ( nspnl == 0 ) return 0.0;
const double omega = basis_.cell().volume();
assert(omega != 0.0);
const double omega_inv = 1.0 / omega;
......@@ -2026,7 +2030,6 @@ double NonLocalPotential::energy(bool compute_hpsi, SlaterDet& dsd,
// reduction of enl across rows
ctxt_.dsum('r',1,1,&enl,1);
sigma_enl = 0.0;
if ( compute_stress )
{
ctxt_.dsum(6,1,&tsum[0],6);
......
src/Species.C
......@@ -1063,24 +1063,8 @@ void Species::initialize_nlcc()
nlccg_spl_.resize(ndft_);
nlccg_spl2_.resize(ndft_);
// expand nlcc to large mesh
const int np = nlcc_.size();
fill( nlcc_spl_.begin(), nlcc_spl_.end(), 0.0 );
copy( nlcc_.begin(), nlcc_.end(), nlcc_spl_.begin() );
// if function is not decaying, keep it at constant value
if ( nlcc_[np-2] <= nlcc_[np-1] )
{
fill( nlcc_spl_.begin()+np, nlcc_spl_.end(), nlcc_[np-1] );
}
// otherwise expand nlcc to large mesh using a
// exponential decay exp(-alpha r)
else
{
const double factor = nlcc_[np-1] / nlcc_[np-2];
for ( int i = np; i < ndft_; i++ )
{
nlcc_spl_[i] = nlcc_spl_[i-1] * factor;
}
}
// compute spline coefficients
spline(ndft_,&rps_spl_[0],&nlcc_spl_[0],0.0,0.0,0,1,&nlcc_spl2_[0]);
......
src/XCPotential.C
......@@ -35,43 +35,46 @@ using namespace std;
XCPotential::XCPotential(const ChargeDensity& cd, const string functional_name,
const Control& ctrl): cd_(cd), vft_(*cd_.vft()), vbasis_(*cd_.vbasis())
{
// copy arrays to resize rhototal_r_ and rhototal_g_
rhototal_r_ = cd_.rhor;
rhototal_g_ = cd_.rhog;
if ( functional_name == "LDA" )
{
xcf_ = new LDAFunctional(cd_.rhor);
xcf_ = new LDAFunctional(rhototal_r_);
}
else if ( functional_name == "VWN" )
{
xcf_ = new VWNFunctional(cd_.rhor);
xcf_ = new VWNFunctional(rhototal_r_);
}
else if ( functional_name == "PBE" )
{
xcf_ = new PBEFunctional(cd_.rhor);
xcf_ = new PBEFunctional(rhototal_r_);
}
else if ( functional_name == "BLYP" )
{
xcf_ = new BLYPFunctional(cd_.rhor);
xcf_ = new BLYPFunctional(rhototal_r_);
}
else if ( functional_name == "PBE0" )
{
const double x_coeff = 1.0 - ctrl.alpha_PBE0;
const double c_coeff = 1.0;
xcf_ = new PBEFunctional(cd_.rhor,x_coeff,c_coeff);
xcf_ = new PBEFunctional(rhototal_r_,x_coeff,c_coeff);
}
else if ( functional_name == "HSE" )
{
xcf_ = new HSEFunctional(cd_.rhor);
xcf_ = new HSEFunctional(rhototal_r_);
}
else if ( functional_name == "RSH" )
{
xcf_ = new RSHFunctional(cd_.rhor,ctrl.alpha_RSH,ctrl.beta_RSH,ctrl.mu_RSH);
xcf_ = new RSHFunctional(rhototal_r_,ctrl.alpha_RSH,ctrl.beta_RSH,ctrl.mu_RSH);
}
else if ( functional_name == "B3LYP" )
{
xcf_ = new B3LYPFunctional(cd_.rhor);
xcf_ = new B3LYPFunctional(rhototal_r_);
}
else if ( functional_name == "BHandHLYP" )
{
xcf_ = new BHandHLYPFunctional(cd_.rhor);
xcf_ = new BHandHLYPFunctional(rhototal_r_);
}
else
{
......@@ -116,6 +119,10 @@ void XCPotential::update(vector<vector<double> >& vr)
// The array cd_.rhog is only used if xcf->isGGA() == true
// to compute the density gradients
// if a non-linear core correction is included,
// rhototal_r_ = cd_.rhor+cd_.rhocore_r. Otherwise
// rhototal_r_ = cd_.rhor
// Output: (through member function xcf())
//
// exc_, dxc_
......@@ -133,6 +140,33 @@ void XCPotential::update(vector<vector<double> >& vr)
// xcf()->vxc2_upup, xcf()->vxc2_dndn,
// xcf()->vxc2_updn, xcf()->vxc2_dnup
rhototal_r_ = cd_.rhor;
rhototal_g_ = cd_.rhog;
// test if a non-linear core correction is used
// if so, add core density to rhototal_r_ and rhototal_g_
if ( !cd_.rhocore_r.empty() )
{
// add core charge
// note: if nspin==2, the cd_.rhocore_{rg} vectors each
// contain half of the total core charge
for ( int ispin = 0; ispin < rhototal_r_.size(); ispin++ )
{
//for ( int i = 0; i < rhototal_r_[ispin].size(); i++ )
// rhototal_r_[ispin][i] += cd_.rhocore_r[i];
int len = rhototal_r_[ispin].size();
int inc1 = 1;
double one = 1.0;
daxpy(&len,&one,(double*)&cd_.rhocore_r[0],&inc1,
&rhototal_r_[ispin][0],&inc1);
//for ( int i = 0; i < rhototal_g_[ispin].size(); i++ )
// rhototal_g_[ispin][i] += cd_.rhocore_g[i];
len = 2*rhototal_g_[ispin].size();
daxpy(&len,&one,(double*)&cd_.rhocore_g[0],&inc1,
(double*)&rhototal_g_[ispin][0],&inc1);
}
}
if ( !xcf_->isGGA() )
{
// LDA functional
......@@ -197,7 +231,7 @@ void XCPotential::update(vector<vector<double> >& vr)
for ( int ig = 0; ig < ngloc_; ig++ )
{
/* i*G_j*c(G) */
tmp1[ig] = complex<double>(0.0,omega_inv*gxj[ig]) * cd_.rhog[0][ig];
tmp1[ig] = complex<double>(0.0,omega_inv*gxj[ig])*rhototal_g_[0][ig];
}
vft_.backward(&tmp1[0],&tmpr[0]);
int inc2=2, inc1=1;
......@@ -210,8 +244,8 @@ void XCPotential::update(vector<vector<double> >& vr)
for ( int j = 0; j < 3; j++ )
{
const double *const gxj = vbasis_.gx_ptr(j);
const complex<double>* rhg0 = &cd_.rhog[0][0];
const complex<double>* rhg1 = &cd_.rhog[1][0];
const complex<double>* rhg0 = &rhototal_g_[0][0];
const complex<double>* rhg1 = &rhototal_g_[1][0];
for ( int ig = 0; ig < ngloc_; ig++ )
{
/* i*G_j*c(G) */
......
src/XCPotential.h
......@@ -37,6 +37,11 @@ class XCPotential
const ChargeDensity& cd_;
XCFunctional* xcf_;
// rhototal_r_: total density (valence+core)
std::vector<std::vector<double> > rhototal_r_;
// rhototal_g_: total density (valence+core)
std::vector<std::vector<std::complex<double> > > rhototal_g_;
std::vector<std::vector<double> > vxctmp; // vxctmp[ispin][ir]
std::vector<std::complex<double> > tmpr; // tmpr[ir]
std::vector<std::complex<double> > tmp1, tmp2; // tmp1[ig], tmp2[ig]
......
util/upf2qso/src/Makefile
CXXFLAGS=-g
upf2qso: upf2qso.o PeriodicTable.o spline.o
upf2qso: upf2qso.o PeriodicTable.o spline.o isodate.o
$(CXX) -o $@ $^
clean:
rm -f *.o upf2qso
util/upf2qso/src/isodate.C 0 → 100644
////////////////////////////////////////////////////////////////////////////////
//
// Copyright (c) 2008 The Regents of the University of California
//
// This file is part of Qbox
//
// Qbox is distributed 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.
// See the file COPYING in the root directory of this distribution
// or <http://www.gnu.org/licenses/>.
//
////////////////////////////////////////////////////////////////////////////////
//
// isodate.C
//
////////////////////////////////////////////////////////////////////////////////
#include "isodate.h"
#include <ctime>
std::string isodate(void)
{
const time_t t = time(NULL);
struct tm* tms = gmtime(&t);
char s[32];
const char* fmt = "%Y-%m-%dT%TZ";
strftime(s,32,fmt,tms);
return std::string(s);
}
util/upf2qso/src/isodate.h 0 → 100644
////////////////////////////////////////////////////////////////////////////////
//
// Copyright (c) 2008 The Regents of the University of California
//
// This file is part of Qbox
//
// Qbox is distributed 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.
// See the file COPYING in the root directory of this distribution
// or <http://www.gnu.org/licenses/>.
//
////////////////////////////////////////////////////////////////////////////////
//
// isodate.h:
//
////////////////////////////////////////////////////////////////////////////////
#ifndef ISODATE_H
#include <string>
std::string isodate(void);
#endif
util/upf2qso/src/upf2qso.C
......@@ -27,12 +27,14 @@
#include <fstream>
#include <sstream>
#include <string>
#include <algorithm>
#include <vector>
#include <cmath>
#include <cassert>
#include <cstdlib>
#include <stdexcept>
#include "spline.h"
#include "isodate.h"
#include "PeriodicTable.h"
using namespace std;
......@@ -113,8 +115,11 @@ void seek_str(string tag)
////////////////////////////////////////////////////////////////////////////////
string get_attr(string buf, string attr)
{
//cerr << " get_attr: searching for: " << attr << endl;
//cerr << " in buffer: " << endl;
//cerr << buf << endl;
bool done = false;
string s, search_string = " " + attr + "=";
string s, search_string = attr + "=";
// find attribute name in buf
string::size_type p = buf.find(search_string);
......@@ -128,7 +133,10 @@ string get_attr(string buf, string attr)
cerr << " get_attr: attribute not found: " << attr << endl;
throw invalid_argument(attr);
}
return buf.substr(b+1,e-b-1);
s = buf.substr(b+1,e-b-1);
// remove white space
s.erase(remove_if(s.begin(),s.end(),::isspace),s.end());
return s;
}
else
{
......@@ -151,7 +159,7 @@ void skipln(void)
}
////////////////////////////////////////////////////////////////////////////////
const string release="1.6";
const string release="1.9";
int main(int argc, char** argv)
{
......@@ -560,6 +568,8 @@ int main(int argc, char** argv)
else
// use value closest to the origin for r=0
nlcc_lin[i] = upf_nlcc[0];
if ( fabs(nlcc_lin[i]) < 1.e-12 )
nlcc_lin[i] = 0.0;
}
}
......@@ -745,7 +755,8 @@ int main(int argc, char** argv)
cout << " xsi:schemaLocation=\"http://www.quantum-simulation.org/ns/fpmd/fpmd-1.0" << endl;
cout << " species.xsd\">" << endl;
cout << "<description>" << endl;
cout << "Translated from UPF format by upf2qso" << endl;
cout << "Translated from UPF format by upf2qso " << release
<< " on " << isodate() << endl;
cout << upf_pp_info;
cout << "</description>" << endl;
cout << "<symbol>" << upf_symbol << "</symbol>" << endl;
......@@ -929,17 +940,6 @@ int main(int argc, char** argv)
find_end_element("PP_RAB");
find_end_element("PP_MESH");
// NLCC
vector<double> upf_nlcc;
if ( upf_nlcc_flag == "T" )
{
find_start_element("PP_NLCC");
upf_nlcc.resize(upf_mesh_size);
for ( int i = 0; i < upf_mesh_size; i++ )
cin >> upf_nlcc[i];
find_end_element("PP_NLCC");
}
find_start_element("PP_LOCAL");
vector<double> upf_vloc(upf_mesh_size);
for ( int i = 0; i < upf_mesh_size; i++ )
......@@ -1085,6 +1085,17 @@ int main(int argc, char** argv)
}
find_end_element("PP_PSWFC");
// NLCC
vector<double> upf_nlcc;
if ( upf_nlcc_flag == "T" )
{
find_start_element("PP_NLCC");
upf_nlcc.resize(upf_mesh_size);
for ( int i = 0; i < upf_mesh_size; i++ )
cin >> upf_nlcc[i];
find_end_element("PP_NLCC");
}
// output original data in file upf.dat
ofstream upf("upf.dat");
upf << "# vloc" << endl;
......@@ -1197,6 +1208,8 @@ int main(int argc, char** argv)
else
// use value closest to the origin for r=0
nlcc_lin[i] = upf_nlcc[0];
if ( fabs(nlcc_lin[i]) < 1.e-12 )
nlcc_lin[i] = 0.0;
}
}
......@@ -1382,7 +1395,8 @@ int main(int argc, char** argv)
cout << " xsi:schemaLocation=\"http://www.quantum-simulation.org/ns/fpmd/fpmd-1.0" << endl;
cout << " species.xsd\">" << endl;
cout << "<description>" << endl;
cout << "Translated from UPF format by upf2qso" << endl;
cout << "Translated from UPF format by upf2qso " << release
<< " on " << isodate() << endl;
cout << upf_pp_info;
cout << "</description>" << endl;
cout << "<symbol>" << upf_symbol << "</symbol>" << endl;
......@@ -1443,6 +1457,17 @@ int main(int argc, char** argv)
if ( pseudo_type == "NC" )
{
// NLCC
vector<double> upf_nlcc;
if ( upf_nlcc_flag == "T" )
{
find_start_element("PP_NLCC");
upf_nlcc.resize(upf_mesh_size);
for ( int i = 0; i < upf_mesh_size; i++ )
cin >> upf_nlcc[i];
find_end_element("PP_NLCC");
}
cerr << " NC potential" << endl;
// output original data in file upf.dat
ofstream upf("upf.dat");
......@@ -1504,6 +1529,8 @@ int main(int argc, char** argv)
else
// use value closest to the origin for r=0
nlcc_lin[i] = upf_nlcc[0];
if ( fabs(nlcc_lin[i]) < 1.e-12 )
nlcc_lin[i] = 0.0;
}
}
......@@ -1541,12 +1568,16 @@ int main(int argc, char** argv)
for ( int j = 0; j < upf_nproj; j++ )
{
// factor 0.5: convert from Ry in UPF to Hartree in QSO
for ( int i = 0; i < upf_vnl.size(); i++ )
f[i] = 0.5 * upf_vnl[j][i];
// interpolate projectors
// note: upf_vnl contains r*projector
// See UPF documentation at http://www.quantum-espresso.org/
// pseudopotentials/unified-pseudopotential-format
assert(f.size()>=upf_vnl[j].size());
for ( int i = 0; i < upf_vnl[j].size(); i++ )
f[i] = upf_vnl[j][i];
int n = upf_vloc.size();
int bcnat_left = 0;
int n = f.size();
int bcnat_left = 1;
double yp_left = 0.0;
int bcnat_right = 1;
double yp_right = 0.0;
......@@ -1559,7 +1590,9 @@ int main(int argc, char** argv)
if ( r >= upf_r[0] )
splint(n,&upf_r[0],&f[0],&fspl[0],r,&vnl_lin[j][i]);
else
vnl_lin[j][i] = 0.5 * upf_vnl[j][0];
vnl_lin[j][i] = upf_vnl[j][0];
if ( fabs(vnl_lin[j][i]) < 1.e-12 )
vnl_lin[j][i] = 0.0;
}
}
......@@ -1586,7 +1619,8 @@ int main(int argc, char** argv)
cout << " xsi:schemaLocation=\"http://www.quantum-simulation.org/ns/fpmd/fpmd-1.0" << endl;
cout << " species.xsd\">" << endl;
cout << "<description>" << endl;
cout << "Translated from UPF format by upf2qso" << endl;
cout << "Translated from UPF format by upf2qso " << release
<< " on " << isodate() << endl;
cout << upf_pp_info;
cout << "</description>" << endl;
cout << "<symbol>" << upf_symbol << "</symbol>" << endl;
......@@ -1612,6 +1646,7 @@ int main(int argc, char** argv)
cout << "</local_potential>" << endl;
// projectors
// note: vnl_lin contains r[i]*projector[i]
int ip = 0;
for ( int l = 0; l <= upf_lmax; l++ )
{
......@@ -1620,8 +1655,25 @@ int main(int argc, char** argv)
cout << "<projector l=\"" << l << "\" i=\""
<< i+1 << "\" size=\"" << nplin << "\">"
<< endl;
for ( int j = 0; j < nplin; j++ )
cout << setprecision(10) << vnl_lin[ip][j] << endl;
// value at r=0:
// quadratic extrapolation of vnl_lin(r)/r to r=0 if l==0
if ( l == 0 )
{
// use quadratic extrapolation to zero
const double h = mesh_spacing;
const double v = (4.0/3.0)*vnl_lin[ip][1]/h -
(1.0/3.0)*vnl_lin[ip][2]/(2*h);
cout << setprecision(10) << v << endl;
}
else
{
cout << setprecision(10) << 0.0 << endl;
}
for ( int j = 1; j < nplin; j++ )
{
const double r = j * mesh_spacing;
cout << setprecision(10) << vnl_lin[ip][j]/r << endl;
}
ip++;
cout << "</projector>" << endl;
}
......@@ -1640,7 +1692,7 @@ int main(int argc, char** argv)
int ij = i + j*upf_nproj;
cout << "<d_ij l=\"" << l << "\""
<< " i=\"" << i-ibase+1 << "\" j=\"" << j-jbase+1
<< "\" " << setprecision(10) << 0.5*upf_d[ij] << " />"
<< "\"> " << setprecision(10) << 0.5*upf_d[ij] << " </d_ij>"
<< endl;
}
}
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment