Commit f47d0c1b by Francois Gygi

New implementation of kpgen

parent 86be185c
kpgen: kpgen.C D3vector.h
kpgen: kpgen.cpp D3vector.h
g++ -o $@ $^
////////////////////////////////////////////////////////////////////////////////
//
// kpgen.C: generate a kpoint list for an arbitrary cell
// kpgen.cpp: generate a kpoint list for an arbitrary cell
// use: kpgen nx ny nz sx sy sz a0x a0y a0z a1x a1y a1z a2x a2y a2z
//
// where a0x a0y a0z = first basis vector of the unit cell
......@@ -8,47 +9,61 @@
//
// nx,ny,nz: number of kpoints in each direction
// sx,sy,sz: shift in each direction (floating point)
// shift: 0: symmetric set: boundary points not included
// if shift==0: symmetric set: boundary points not included
// even-numbered sets do not include the gamma point
// odd-numbered sets include the gamma point
//
////////////////////////////////////////////////////////////////////////////////
#include<iostream>
#include<fstream>
#include<iomanip>
#include<vector>
#include<cassert>
#include<cstdlib>
#include<list>
#include "D3vector.h"
using namespace std;
// ib_BZ: test if the vector k is in the BZ defined by b0,b1,b2
////////////////////////////////////////////////////////////////////////////////
// comparison function for k-points
// two k-points are equal if k1 = k2+G for any G in the first shells
bool equals(D3vector k1, D3vector k2, D3vector b0, D3vector b1, D3vector b2)
{
for ( int i0 = -2; i0 <= 2; i0++ )
for ( int i1 = -2; i1 <= 2; i1++ )
for ( int i2 = -2; i2 <= 2; i2++ )
{
if ( length(k1-k2-i0*b0-i1*b1-i2*b2) < 1.e-5 )
return true;
}
return false;
}
////////////////////////////////////////////////////////////////////////////////
// check if a k-point is in the BZ
bool in_BZ(D3vector k, D3vector b0, D3vector b1, D3vector b2)
{
const double epsilon = 1.e-8;
// check projection of kpoint k along all reciprocal lattice vectors
// in the first two shells
// use a shift in an irrational direction epsilon*(1,M_PI,M_LN2)
// to avoid including zone boundary equivalent vectors
const double epsilon = 1.e-6;
D3vector kshifted = k + epsilon * D3vector(1.0,M_PI,M_LN2);
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++ )
for ( int i0 = -2; i0 <= 2; i0++ )
for ( int i1 = -2; i1 <= 2; i1++ )
for ( int i2 = -2; i2 <= 2; i2++ )
if ( !(i0 == 0 && i1 == 0 && i2 == 0) )
{
D3vector g = i0 * b0 + i1 * b1 + i2 * b2;
if ( kshifted*g > 0.5 * g*g )
in_bz = false;
return false;
}
return in_bz;
return true;
}
////////////////////////////////////////////////////////////////////////////////
int main(int argc, char** argv)
{
cout << "# kpgen " << endl;
cout << "# kpgen 2.0" << endl;
if ( argc != 16 )
{
cerr << " use: " << argv[0] << " nx ny nz shiftx shifty shiftz {cell}"
......@@ -91,8 +106,9 @@ int main(int argc, char** argv)
D3vector b1 = fac * a2 ^ a0;
D3vector b2 = fac * a0 ^ a1;
list<vector<int> > kplist;
vector<int> kpint(4);
vector<D3vector> kp;
vector<D3vector> kpfrac;
vector<double> weight;
// scan volume enclosing the BZ
for ( int ii = -2; ii <= 2; ii++ )
......@@ -104,109 +120,80 @@ int main(int argc, char** argv)
{
for ( int k = 0; k < nz; k++ )
{
kpint[0] = ii*2*nx + 2*i-nx+1;
kpint[1] = jj*2*ny + 2*j-ny+1;
kpint[2] = kk*2*nz + 2*k-nz+1;
kpint[3] = 1;
int kpint0 = ii*2*nx + 2*i-nx+1;
int kpint1 = jj*2*ny + 2*j-ny+1;
int kpint2 = kk*2*nz + 2*k-nz+1;
D3vector kv = ( (kpint[0] + sx)/(2.0*nx) ) * b0 +
( (kpint[1] + sy)/(2.0*ny) ) * b1 +
( (kpint[2] + sz)/(2.0*nz) ) * b2;
double kv0 = (kpint0 + sx)/(2.0*nx);
double kv1 = (kpint1 + sy)/(2.0*ny);
double kv2 = (kpint2 + sz)/(2.0*nz);
D3vector kv = kv0*b0 + kv1*b1 + kv2*b2;
if ( in_BZ(kv,b0,b1,b2) )
kplist.push_back(kpint);
{
kp.push_back(kv);
kpfrac.push_back(D3vector(kv0,kv1,kv2));
weight.push_back(1.0);
}
}
}
}
int total_weight = kp.size();
// check for equivalent vectors
// count vectors that are equivalent to k+G
cout.setf(ios::fixed,ios::floatfield);
int nequiv = 0;
for ( int i = 0; i < kp.size(); i++ )
for ( int j = i+1; j < kp.size(); j++ )
if ( equals(kp[i],kp[j],b0,b1,b2) )
{
nequiv++;
//cout << setprecision(3)
// << kpfrac[i] << " equiv " << kpfrac[j] << endl;
int total_weight = kplist.size();
#if 1
// remove -k
for (list<vector<int> >::iterator i = kplist.begin(); i != kplist.end(); i++)
}
if ( nequiv != 0 )
{
// test if -k is in the set (and is not 0 0 0)
kpint[0] = (*i)[0];
kpint[1] = (*i)[1];
kpint[2] = (*i)[2];
kpint[3] = (*i)[3];
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++ )
// there should not be any equivalent points as k=k+G
cerr << nequiv << " error: equivalent points (k=k+G)" << endl;
return 1;
}
// count vectors that are equivalent to -k+G
int nequivm = 0;
for ( int i = 0; i < kp.size(); i++ )
for ( int j = i+1; j < kp.size(); j++ )
if ( equals(kp[i],-kp[j],b0,b1,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 weight to (*i)
(*i)[3] += (*j)[3];
(*j)[3] = 0;
}
nequivm++;
//cout << setprecision(3)
// << kpfrac[i] << " equiv " << kpfrac[j] << endl;
}
}
}
#if DEBUG
cout << nequivm << " equivalent points (k=-k+G)" << endl;
#endif
#if 1
// remove duplicate points
for (list<vector<int> >::iterator i = kplist.begin(); i != kplist.end(); i++)
// reassign weight from equivalent points
for ( int i = 0; i < kp.size(); i++ )
{
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 ( weight[i] != 0.0 )
{
if ( j != i )
for ( int j = i+1; j < kp.size(); j++ )
{
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 )
if ( equals(kp[i],-kp[j],b0,b1,b2) )
{
// transfer the weight of kj to ki
(*i)[3] += (*j)[3];
(*j)[3] = 0;
// reassign the weight of kp[j] to kp[i]
weight[i] += weight[j];
weight[j] = 0.0;
}
}
}
}
#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++)
{
sum += (*i)[3];
}
assert(sum==total_weight);
#endif
#if 1
// 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 */)
{
int w = (*i)[3];
if ( w == 0 )
kplist.erase(i++);
else
++i;
}
#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);
cout << "# nx,ny,nz: " << nx << " " << ny << " " << nz << endl;
......@@ -218,70 +205,49 @@ int main(int argc, char** argv)
cout << "# b1/(2pi): " << b1/(2*M_PI) << endl;
cout << "# b2/(2pi): " << b2/(2*M_PI) << endl;
cout << "# " << kplist.size() << " k-points" << endl;
cout << "# " << kp.size()-nequivm << " k-points" << endl;
cout << " kpoint delete 0 0 0" << endl;
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) << kx << " "
<< setw(13) << ky << " "
<< setw(13) << kz << " ";
cout.setf(ios::scientific,ios::floatfield);
cout << setprecision(14)
<< setw(16) << w << endl;
}
#endif
#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++ )
// print list backward to have increasing x components
for ( int i = kpfrac.size()-1; i >= 0; i-- )
{
for ( int jj = -ny; jj <= ny; jj++ )
if ( weight[i] != 0.0 )
{
for ( int kk = -nz; kk <= nz; kk++ )
cout.setf(ios::fixed,ios::floatfield);
// compute projections along b0, b1, b2
double kx = kpfrac[i].x;
double ky = kpfrac[i].y;
double kz = kpfrac[i].z;
// print -k to have positive coefficients in output
// change sign only if component is non-zero to avoid -0.00 in output
if ( kx == 0.0 )
{
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++)
if ( ky == 0.0 )
{
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 ( kz < 0.0 ) kz = -kz;
}
else if ( ky < 0.0 )
{
ky = -ky;
if ( kz != 0 ) kz = -kz;
}
if ( len != 0 && fabs(sum) > 1.e-6 )
minlen = min(minlen,len);
}
else if ( kx < 0.0 )
{
kx = -kx;
if ( ky != 0 ) ky = -ky;
if ( kz != 0 ) kz = -kz;
}
double w = weight[i]/((double) total_weight);
cout << " kpoint add "
<< setprecision(10)
<< setw(13) << kx << " "
<< setw(13) << ky << " "
<< setw(13) << kz << " ";
cout.setf(ios::scientific,ios::floatfield);
cout << setprecision(14)
<< setw(16) << w << endl;
}
}
cerr << " smallest R with non-zero sum has length " << minlen << endl;
#endif
}
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