Commit 9305cc1b by Martin Schlipf

Merge rel1_57_14 into hse-dev branch

Conflicts:
	src/ExchangeOperator.C

git-svn-id: http://qboxcode.org/svn/qb/branches/hse-dev@1405 cba15fb0-1239-40c8-b417-11db7ca47a34
parent 5117e207
......@@ -953,6 +953,9 @@ void BOSampleStepper::step(int niter)
<< "\" kpoint=\""
<< setprecision(8)
<< wf.sd(ispin,ikp)->kpoint()
<< "\" weight=\""
<< setprecision(8)
<< wf.weight(ikp)
<< "\" n=\"" << nst << "\">" << endl;
for ( int i = 0; i < nst; i++ )
{
......
......@@ -373,6 +373,8 @@ void ExchangeOperator::add_stress(valarray<double>& sigma_exc)
void ExchangeOperator::cell_moved(void)
{
vbasis_->resize( s_.wf.cell(),s_.wf.refcell(),4.0*s_.wf.ecut());
KPGridStat_.cell_moved();
KPGridPerm_.cell_moved();
}
////////////////////////////////////////////////////////////////////////////////
......@@ -688,7 +690,8 @@ double ExchangeOperator::compute_exchange_for_general_case_( Sample* s,
2.0 * occ_kj_[j] * spinFactor;
// add contribution to exchange energy
exchange_sum += ex_ki_i_kj_j * weight;
exchange_sum += ex_ki_i_kj_j * weight * wf.weight(iKpi) *
0.5 * occ_ki_[i];
if (dwf)
{
......@@ -727,8 +730,10 @@ double ExchangeOperator::compute_exchange_for_general_case_( Sample* s,
2.0 * occ_kj_[j] * spinFactor;
// add contribution to exchange energy
exchange_sum += ex_ki_i_kj_j * weightj;
exchange_sum += ex_ki_i_kj_j * weighti;
exchange_sum += ex_ki_i_kj_j * weightj * wf.weight(iKpi) *
0.5 * occ_ki_[i];
exchange_sum += ex_ki_i_kj_j * weighti * wf.weight(iKpj) *
0.5 * occ_kj_[j];
if (dwf)
{
......@@ -790,7 +795,7 @@ double ExchangeOperator::compute_exchange_for_general_case_( Sample* s,
CompleteSendingForces(iRotationStep);
// add the g coordinate expression contributions to each
// state derivative of kpi, store everithing in the send buffer
// state derivative of kpi, store everything in the send buffer
for ( int i=0; i<nStatesKpi_; i++ )
{
// transform contribution to g coordinates
......@@ -882,7 +887,7 @@ double ExchangeOperator::compute_exchange_for_general_case_( Sample* s,
// no contribution to stress from div_corr_3
}
// Quadratic corrections
// Quadratic corrections
if ( quad_correction )
{
// beta_x, beta_y and beta_z: curvature of rho(G=0)
......@@ -958,7 +963,7 @@ double ExchangeOperator::compute_exchange_for_general_case_( Sample* s,
if ( gcontext_.onpe0() )
{
cout << setprecision(10);
cout << " total exchange = " << extot << " Eh\n";
cout << " total exchange = " << extot << " (a.u.)\n";
cout << " total exchange computation time: " << tm.real()
<< " s" << endl;
}
......@@ -2204,13 +2209,12 @@ double ExchangeOperator::compute_exchange_at_gamma_(const Wavefunction &wf,
div_corr += div_corr_4;
const double e_div_corr_4 = -0.5 * div_corr_4 * occ_ki_[i];
exchange_sum += e_div_corr_4;
const double fac4 = ( 4.0 * M_PI / omega );
sigma_exhf_[0] += ( e_div_corr_4 + fac4 * 2.0 * beta_x ) / omega;
sigma_exhf_[1] += ( e_div_corr_4 + fac4 * 2.0 * beta_y ) / omega;
sigma_exhf_[2] += ( e_div_corr_4 + fac4 * 2.0 * beta_z ) / omega;
const double fac4 = ( 4.0 * M_PI / omega ) * occ_ki_[i];
sigma_exhf_[0] += ( e_div_corr_4 + fac4 * beta_x ) / omega;
sigma_exhf_[1] += ( e_div_corr_4 + fac4 * beta_y ) / omega;
sigma_exhf_[2] += ( e_div_corr_4 + fac4 * beta_z ) / omega;
}
// contribution of divergence corrections to forces on wave functions
// (other than Coulomb, no divergence correction for wave functions)
if (dwf and coulomb_)
......@@ -2248,8 +2252,9 @@ double ExchangeOperator::compute_exchange_at_gamma_(const Wavefunction &wf,
// scale stress tensor with HF coefficient
sigma_exhf_ *= HFCoeff_;
// print result
tm.stop();
#ifdef DEBUG
if ( gcontext_.onpe0() )
{
cout << setprecision(10);
......@@ -2279,6 +2284,7 @@ double ExchangeOperator::compute_exchange_at_gamma_(const Wavefunction &wf,
<< endl;
}
}
#endif
// return total exchange in Hartree, scaled by HF coefficient
return extot;
......
......@@ -30,12 +30,14 @@ class KPConnectivity
{
private:
const Sample &s_;
int DimX_;
int DimY_;
int DimZ_;
int nkpoints_;
int nStateMax_;
int ConnectivityComplete_;
double kdist_tol_;
std::vector<int> first_neighbour_kx_;
std::vector<int> first_neighbour_ky_;
......@@ -128,6 +130,7 @@ class KPConnectivity
void AddOverlap(int iKpi, int iKpj, int iLocStatei,
complex<double> *valueDirect, complex<double> *valueSymmetric,
double occupation);
void cell_moved(void);
void StartPermutation(int iKp, int iSendTo, int iRecvFr);
void EndPermutation();
......
......@@ -75,16 +75,20 @@ void SampleReader::readSample (Sample& s, const string uri, bool serial)
catch (const XMLException& toCatch)
{
cout << " Sample::readSample: Error during XML initialization :\n"
cout << " SampleReader::readSample: Error during XML initialization :\n"
<< StrX(toCatch.getMessage()) << endl;
s.ctxt_.abort(1);
return;
}
string xmlcontent;
DoubleMatrix gfdata(ctxt_);
XMLGFPreprocessor xmlgfp;
xmlgfp.process(uri.c_str(),gfdata,xmlcontent,serial);
if ( xmlgfp.process(uri.c_str(),gfdata,xmlcontent,serial) )
{
cout << " SampleReader::readSample: Error in XMLGFPreprocessor" << endl;
return;
}
// Each task holds a copy of xmlcontent
// The distributed matrix gfdata contains the grid function data
......
......@@ -31,13 +31,13 @@ class Species
int nlm_; // number of non-local projectors:
int ndft_;
std::vector<std::vector<double> > vps_spl_, phi_spl_;
std::vector<double> gspl_, vlocg_, vlocg_spl;
std::vector<std::vector<double> > vnlg_, vnlg_spl;
std::vector<std::vector<double> > vps_spl_, vps_spl2_, phi_spl_, phi_spl2_;
std::vector<double> gspl_, vlocg_spl_, vlocg_spl2_;
std::vector<std::vector<double> > vnlg_spl_, vnlg_spl2_;
std::vector<double> wsg_; // wsg_[l] Kleinman-Bylander weight
// 1/<phi|delta_V|phi>
std::vector<double> rps_; // radial linear mesh (same for all l)
std::vector<double> rps_spl_; // radial linear mesh (same for all l)
std::string name_; // name used in current application
std::string uri_; // uri of the resource defining the pseudopotential
......@@ -53,8 +53,8 @@ class Species
int nquad_; // number of semi-local quadrature points
double rquad_; // end of semi-local quadrature interval
double deltar_; // mesh spacing for potentials and wavefunctions
std::vector<std::vector<double> > vps_; // potentials for each l
std::vector<std::vector<double> > phi_; // atomic wavefunctions for each l
std::vector<std::vector<double> > vps_; // potentials for each l (input)
std::vector<std::vector<double> > phi_; // atomic wf for each l (input)
double rcps_; // cutoff radius of gaussian pseudocharge
public:
......@@ -93,8 +93,8 @@ class Species
double wsg(int l) { return wsg_[l]; };
double rcut_loc(double epsilon); // radius beyond which potential is local
const std::vector<std::vector<double> >& vps(void) const { return vps_; }
const std::vector<std::vector<double> >& phi(void) const { return phi_; }
const std::vector<std::vector<double> >& vps(void) const { return vps_spl_; }
const std::vector<std::vector<double> >& phi(void) const { return phi_spl_; }
bool initialize(double rcps);
void info(std::ostream& os);
......
......@@ -36,6 +36,6 @@ class XMLGFPreprocessor
{
public:
void process(const char* const filename,
int process(const char* const filename,
DoubleMatrix& gfdata, std::string& xmlcontent, bool serial);
};
--------------------------------------------------------------------------------
rel1_57_14
--------------------------------------------------------------------------------
r1396: util/kpgen/kpgen.C: fixed odd numbered case and added test.
r1395: Species.C: fix for case of missing radial functions
r1394: BOSampleStepper.C: include weight info in eigenvalues
r1390,r1391: Fix exchange calculation with multiple k-points
r1388: Species.[Ch]: ensure that print outputs the same data as was read
from the pseudopotential file.
r1363: XMLGFPreprocessor.C: implemented sample loading using http connection.
--------------------------------------------------------------------------------
rel1_57_13
--------------------------------------------------------------------------------
r1358: ExchangeOperator.C: Fixed exchange contributions to stress.
KPGridConnectivity.[Ch]: Added cell_moved() member to update cell-dependent
quantities during cell optimization.
--------------------------------------------------------------------------------
rel1_57_12
--------------------------------------------------------------------------------
r1354: Added cell_moved() member in ExchangeOperator to resize vbasis
......
......@@ -19,5 +19,5 @@
#include "release.h"
std::string release(void)
{
return std::string("1.57.12");
return std::string("1.57.14");
}
......@@ -70,7 +70,7 @@ int main(int argc, char **argv)
cout << " testSpecies: testing SpeciesReader::uri_to_string: " << endl;
string xmlstr;
rdr.uri_to_string(uri,xmlstr);
rdr.uri_to_string(uri,"unknown",xmlstr);
cout << xmlstr;
#if 0
......
......@@ -29,15 +29,18 @@ bool in_BZ(D3vector k, D3vector b0, D3vector b1, D3vector b2)
D3vector g;
// check projection of kpoint k along all 26 reciprocal lattice vectors
// that are nearest g=0
// use a shift by epsilon*(1,1,1) to avoid including zone boundary
// equivalent vectors
bool in_bz = true;
D3vector kshifted = k + epsilon * D3vector(1.0,1.0,1.0);
for ( int i0 = -1; i0 <= 1; i0++ )
for ( int i1 = -1; i1 <= 1; i1++ )
for ( int i2 = -1; i2 <= 1; i2++ )
if ( !(i0 == 0 && i1 == 0 && i2 == 0) )
{
D3vector g = i0 * b0 + i1 * b1 + i2 * b2;
if ( k*g > 0.5 * g*g + epsilon )
if ( kshifted*g > 0.5 * g*g )
in_bz = false;
}
return in_bz;
......@@ -45,12 +48,11 @@ bool in_BZ(D3vector k, D3vector b0, D3vector b1, D3vector b2)
int main(int argc, char** argv)
{
cout << "# kpgen-1.0" << endl;
cout << "# kpgen-1.1" << endl;
if ( argc != 16 )
{
cerr << " use: " << argv[0] << " nx ny nz shiftx shifty shiftz {cell}"
<< endl;
cerr << " shift==0: symmetric set, zone boundary not included" << endl;
return 1;
}
int nx = atoi(argv[1]);
......@@ -95,7 +97,7 @@ int main(int argc, char** argv)
// scan volume enclosing the BZ
for ( int ii = -2; ii <= 2; ii++ )
for ( int jj = -2; jj <= 2; jj++ )
for ( int kk = -1; kk <= 2; kk++ )
for ( int kk = -2; kk <= 2; kk++ )
for ( int i = 0; i < nx; i++ )
{
for ( int j = 0; j < ny; j++ )
......@@ -107,9 +109,9 @@ int main(int argc, char** argv)
kpint[2] = kk*2*nz + 2*k-nz+1;
kpint[3] = 1;
D3vector k = kpint[0]/(2.0*nx) * b0 +
kpint[1]/(2.0*ny) * b1 +
kpint[2]/(2.0*nz) * b2;
D3vector k = ( (kpint[0] + sx)/(2.0*nx) ) * b0 +
( (kpint[1] + sy)/(2.0*ny) ) * b1 +
( (kpint[2] + sz)/(2.0*nz) ) * b2;
if ( in_BZ(k,b0,b1,b2) )
kplist.push_back(kpint);
......@@ -128,19 +130,22 @@ int main(int argc, char** argv)
kpint[1] = (*i)[1];
kpint[2] = (*i)[2];
kpint[3] = (*i)[3];
if ( kpint[0]*kpint[0]+kpint[1]*kpint[1]+kpint[2]*kpint[2] != 0 )
D3vector ki = (((*i)[0]+sx)/(2.0*nx)) * b0 +
(((*i)[1]+sy)/(2.0*ny)) * b1 +
(((*i)[2]+sz)/(2.0*nz)) * b2;
if ( length(ki) != 0.0 )
{
// look for -k in the rest of the list
for ( list<vector<int> >::iterator j = i; j != kplist.end(); j++ )
{
if ( (*j)[0]==-kpint[0] && (*j)[1]==-kpint[1] && (*j)[2]==-kpint[2] )
D3vector kj = (((*j)[0]+sx)/(2.0*nx)) * b0 +
(((*j)[1]+sy)/(2.0*ny)) * b1 +
(((*j)[2]+sz)/(2.0*nz)) * b2;
if ( length(ki+kj) < 1.e-5 )
{
// transfer weight to (*i)
(*i)[3] += (*j)[3];
(*j)[3] = 0;
//cout << " erasing " << "(" << kpint[0] << ","
// << kpint[1] << "," << kpint[2] << ") == -("
// << (*j)[0] << "," << (*j)[1] << "," << (*j)[2] << ")" << endl;
}
}
}
......@@ -148,76 +153,59 @@ int main(int argc, char** argv)
#endif
#if 1
// remove equivalent points
// remove duplicate points
for (list<vector<int> >::iterator i = kplist.begin(); i != kplist.end(); i++)
{
D3vector ki = (*i)[0]/(2.0*nx) * b0 +
(*i)[1]/(2.0*ny) * b1 +
(*i)[2]/(2.0*nz) * b2;
// look for equivalent points in the rest of the list
D3vector ki = (((*i)[0]+sx)/(2.0*nx)) * b0 +
(((*i)[1]+sy)/(2.0*ny)) * b1 +
(((*i)[2]+sz)/(2.0*nz)) * b2;
// look for duplicate points in the rest of the list
for ( list<vector<int> >::iterator j = i; j != kplist.end(); j++ )
{
if ( j != i )
{
D3vector kj = (*j)[0]/(2.0*nx) * b0 +
(*j)[1]/(2.0*ny) * b1 +
(*j)[2]/(2.0*nz) * b2;
D3vector kj = (((*j)[0]+sx)/(2.0*nx)) * b0 +
(((*j)[1]+sy)/(2.0*ny)) * b1 +
(((*j)[2]+sz)/(2.0*nz)) * b2;
if ( length(ki-kj) < 1.e-5 )
{
// transfer the weight of kj to ki
(*i)[3] += (*j)[3];
(*j)[3] = 0;
//cout << " erasing equivalent point " << "(" << (*j)[0] << ","
// << (*j)[1] << "," << (*j)[2] << ") == ("
// << (*i)[0] << "," << (*i)[1] << "," << (*i)[2] << ")" << endl;
}
}
}
}
#endif
// remove elements with zero weight
for (list<vector<int> >::iterator i = kplist.begin();
i != kplist.end(); /* nothing */ )
{
if ( (*i)[3] == 0 )
{
kplist.erase(i++);
}
else
{
i++;
}
}
#if 1
// make first index positive
// check that the sum of weights is one
int sum = 0;
for (list<vector<int> >::iterator i = kplist.begin(); i != kplist.end(); i++)
{
D3vector ki = (*i)[0]/(2.0*nx) * b0 +
(*i)[1]/(2.0*ny) * b1 +
(*i)[2]/(2.0*nz) * b2;
if ( ki.x < 0 )
{
(*i)[0] *= -1;
(*i)[1] *= -1;
(*i)[2] *= -1;
}
sum += (*i)[3];
}
assert(sum==total_weight);
#endif
#if 1
// check that the sum of weights is one
int sum = 0;
for (list<vector<int> >::iterator i = kplist.begin(); i != kplist.end(); i++)
// remove elements with zero weight
// For a description of safe erase, see S. Meyers, Effective STL, item 9
for (list<vector<int> >::iterator i = kplist.begin();
i != kplist.end(); /* nothing */)
{
sum += (*i)[3];
int w = (*i)[3];
if ( w == 0 )
kplist.erase(i++);
else
++i;
}
assert(sum==total_weight);
#endif
#if 1
// output list
// change the sign of the k vector if the first element is negative
// traverse list backwards to have increasing indices
// kpoints are output in reciprocal lattice coordinates
cout.setf(ios::right,ios::adjustfield);
......@@ -235,32 +223,65 @@ int main(int argc, char** argv)
for (list<vector<int> >::reverse_iterator i = kplist.rbegin();
i != kplist.rend(); i++)
{
D3vector ki = ( ((*i)[0]+sx)/(2.0*nx) ) * b0 +
( ((*i)[1]+sy)/(2.0*ny) ) * b1 +
( ((*i)[2]+sz)/(2.0*nz) ) * b2;
cout.setf(ios::fixed,ios::floatfield);
double kx = ((*i)[0]+sx)/(2.0*nx);
double ky = ((*i)[1]+sy)/(2.0*ny);
double kz = ((*i)[2]+sz)/(2.0*nz);
double w = (*i)[3]/((double) total_weight);
if ( ki.x < 0.0 )
{
kx = -kx;
// next lines: test before changing sign to avoid -0.0
if ( ky != 0.0 ) ky = -ky;
if ( kz != 0.0 ) kz = -kz;
}
cout << " kpoint add "
<< setprecision(10)
<< setw(13) << ((*i)[0]+sx)/(2.0*nx) << " "
<< setw(13) << ((*i)[1]+sy)/(2.0*ny) << " "
<< setw(13) << ((*i)[2]+sz)/(2.0*nz) << " ";
<< setw(13) << kx << " "
<< setw(13) << ky << " "
<< setw(13) << kz << " ";
cout.setf(ios::scientific,ios::floatfield);
cout << setprecision(14)
<< setw(16) << (*i)[3]/((double) total_weight)
<< endl;
<< setw(16) << w << endl;
}
#endif
// output list in absolute coordinates for plot
ofstream plotfile("kpoint.plt");
for (list<vector<int> >::reverse_iterator i = kplist.rbegin();
i != kplist.rend(); i++)
#if 1
// test the k-point set
// compute the numerical integral of the function exp(ikR) for
// a set of vectors R. The integral should be zero, except for R=0 and
// vectors R of large norm (larger than the nx,ny,nz parameters used).
double minlen = 1.e10;
for ( int ii = -nx; ii <= nx; ii++ )
{
D3vector k = (((*i)[0]+sx)/(2.0*nx)) * b0 / (2*M_PI) +
(((*i)[1]+sy)/(2.0*ny)) * b1 / (2*M_PI) +
(((*i)[2]+sz)/(2.0*nz)) * b2 / (2*M_PI);
plotfile.setf(ios::fixed,ios::floatfield);
plotfile << setprecision(8)
<< setw(13) << k.x << " "
<< setw(13) << k.y << " "
<< setw(13) << k.z << endl;
for ( int jj = -ny; jj <= ny; jj++ )
{
for ( int kk = -nz; kk <= nz; kk++ )
{
D3vector R = ii * a0 + jj * a1 + kk * a2;
double len = length(R);
double sum = 0.0;
for (list<vector<int> >::iterator ikp = kplist.begin();
ikp != kplist.end(); ikp++)
{
D3vector k = ( ((*ikp)[0]+sx)/(2.0*nx) ) * b0 +
( ((*ikp)[1]+sy)/(2.0*ny) ) * b1 +
( ((*ikp)[2]+sz)/(2.0*nz) ) * b2;
double w = (*ikp)[3]/((double) total_weight);
sum += w * cos(k*R);
}
if ( len != 0 && fabs(sum) > 1.e-6 )
minlen = min(minlen,len);
}
}
}
cerr << " smallest R with non-zero sum has length " << minlen << endl;
#endif
}
#!/usr/bin/python
# qbox_dos.py: extract dos from Qbox output
# generate dos plot in gnuplot format
# use: qbox_dos.py emin emax width file.r
# the DOS is accumulated separately for each spin
import xml.sax
import sys
import math
if len(sys.argv) != 5:
print "use: ",sys.argv[0]," emin emax width file.r"
sys.exit()
emin = float(sys.argv[1])
emax = float(sys.argv[2])
width = float(sys.argv[3])
infile = sys.argv[4]
ndos = 501
de = (emax - emin)/(ndos-1)
# normalized gaussian distribution in one dimension
# f(x) = 1/(sqrt(pi)*width) * exp(-(x/width)^2 )
def gauss(x, width):
return (1.0/(math.sqrt(math.pi)*width)) * math.exp(-(x/width)**2)
# Qbox output handler to extract and process data
class QboxOutputHandler(xml.sax.handler.ContentHandler):
def __init__(self):
self.nspin = 1
self.readData = 0
self.dos_up = [0] * ndos
self.dos_dn = [0] * ndos
def startElement(self, name, attributes):
if name == "eigenvalues":
self.n = attributes["n"]
self.spin = int(attributes["spin"])
self.kpoint = attributes["kpoint"]
self.weight = float(attributes["weight"])
self.readData = 1
self.buffer = ""
if self.spin == 1:
self.nspin = 2
def characters(self, data):
if self.readData:
self.buffer += data
def endElement(self, name):
if name == "eigenvalues":
self.readData = 0
self.accumulate_dos()
def accumulate_dos(self):
self.e = self.buffer.split()
if self.spin == 0:
for i in range(len(self.e)):
for j in range(ndos):
ej = emin + j * de
self.dos_up[j] += gauss(float(self.e[i])-ej, width ) * self.weight
if self.spin == 1:
for i in range(len(self.e)):
for j in range(ndos):
ej = emin + j * de
self.dos_dn[j] += gauss(float(self.e[i])-ej, width ) * self.weight
def print_dos(self):
print "# ",infile," spin=0 width=",width
for j in range(ndos):
ej = emin + j * de
print ej, self.dos_up[j]
if self.nspin == 2:
print
print
print "# ",infile," spin=1 width=",width
for j in range(ndos):
ej = emin + j * de
print ej, self.dos_dn[j]
parser = xml.sax.make_parser()
handler = QboxOutputHandler()
parser.setContentHandler(handler)
parser.parse(infile)
handler.print_dos()
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