kpgen.C 7.6 KB
Newer Older
Francois Gygi committed
1 2 3 4 5 6 7 8 9
//
// kpgen.C: 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
//       a1x a1y a1z = second basis vector of the unit cell
//       a2x a2y a2z = third basis vector of the unit cell
//
//       nx,ny,nz: number of kpoints in each direction
10
//       sx,sy,sz: shift in each direction (floating point)
Francois Gygi committed
11
//       shift: 0: symmetric set: boundary points not included
12 13
// even-numbered sets do not include the gamma point
// odd-numbered sets include the gamma point
Francois Gygi committed
14 15
//
#include<iostream>
16
#include<fstream>
Francois Gygi committed
17 18 19
#include<iomanip>
#include<vector>
#include<cassert>
Francois Gygi committed
20
#include<cstdlib>
Francois Gygi committed
21 22 23 24
#include<list>
#include "D3vector.h"
using namespace std;

25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45
// ib_BZ: test if the vector k is in the BZ defined by b0,b1,b2
bool in_BZ(D3vector k, D3vector b0, D3vector b1, D3vector b2)
{
  const double epsilon = 1.e-8;
  D3vector g;
  // check projection of kpoint k along all 26 reciprocal lattice vectors
  // that are nearest g=0

  bool in_bz = true;
  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 )
          in_bz = false;
      }
  return in_bz;
}

Francois Gygi committed
46 47
int main(int argc, char** argv)
{
48
  cout << "# kpgen-1.0" << endl;
Francois Gygi committed
49 50 51 52 53 54 55 56 57 58
  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]);
  int ny = atoi(argv[2]);
  int nz = atoi(argv[3]);
59 60 61
  double sx = atof(argv[4]);
  double sy = atof(argv[5]);
  double sz = atof(argv[6]);
Francois Gygi committed
62 63 64 65 66 67 68 69 70 71 72 73 74
  if ( nx <= 0 || ny <= 0 || nz <=0 )
  {
    cerr << " use: " << argv[0] << " nx ny nz shiftx shifty shiftz {cell}"
         << endl;
    cerr << " nx, ny, nz must be positive" << endl;
    return 1;
  }
  if ( sx < 0 || sx > 1 ||
       sy < 0 || sy > 1 ||
       sz < 0 || sz > 1 )
  {
    cerr << " use: " << argv[0] << " nx ny nz shiftx shifty shiftz {cell}"
         << endl;
75
    cerr << " shifts must be in [0,1]" << endl;
Francois Gygi committed
76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93
    return 1;
  }
  D3vector a0(atof(argv[7]),atof(argv[8]),atof(argv[9]));
  D3vector a1(atof(argv[10]),atof(argv[11]),atof(argv[12]));
  D3vector a2(atof(argv[13]),atof(argv[14]),atof(argv[15]));
  const double volume = a0 * ( a1 ^ a2 );
  if ( volume == 0.0 )
  {
    cout << " cell volume is zero" << endl;
    return 1;
  }
  double fac = 2.0 * M_PI / volume;
  D3vector b0 = fac * a1 ^ a2;
  D3vector b1 = fac * a2 ^ a0;
  D3vector b2 = fac * a0 ^ a1;

  list<vector<int> > kplist;
  vector<int> kpint(4);
94 95 96 97 98 99

  // 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 i = 0; i < nx; i++ )
Francois Gygi committed
100
  {
101
    for ( int j = 0; j < ny; j++ )
Francois Gygi committed
102
    {
103
      for ( int k = 0; k < nz; k++ )
Francois Gygi committed
104
      {
105 106 107
        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;
Francois Gygi committed
108
        kpint[3] = 1;
109 110 111 112 113 114 115

        D3vector k = kpint[0]/(2.0*nx) * b0 +
                     kpint[1]/(2.0*ny) * b1 +
                     kpint[2]/(2.0*nz) * b2;

        if ( in_BZ(k,b0,b1,b2) )
          kplist.push_back(kpint);
Francois Gygi committed
116 117 118 119
      }
    }
  }

120 121 122 123 124
  int total_weight = kplist.size();

#if 1
  // remove -k
  for (list<vector<int> >::iterator i = kplist.begin(); i != kplist.end(); i++)
Francois Gygi committed
125
  {
126
    // test if -k is in the set (and is not 0 0 0)
Francois Gygi committed
127 128 129 130
    kpint[0] = (*i)[0];
    kpint[1] = (*i)[1];
    kpint[2] = (*i)[2];
    kpint[3] = (*i)[3];
131
    if ( kpint[0]*kpint[0]+kpint[1]*kpint[1]+kpint[2]*kpint[2] != 0 )
Francois Gygi committed
132
    {
133 134
      // look for -k in the rest of the list
      for ( list<vector<int> >::iterator j = i; j != kplist.end(); j++ )
Francois Gygi committed
135
      {
136
        if ( (*j)[0]==-kpint[0] && (*j)[1]==-kpint[1] && (*j)[2]==-kpint[2] )
Francois Gygi committed
137
        {
138 139 140
          // transfer weight to (*i)
          (*i)[3] += (*j)[3];
          (*j)[3] = 0;
Francois Gygi committed
141
          //cout << " erasing  " << "(" << kpint[0] << ","
142 143
          //     << kpint[1] << "," << kpint[2] << ") == -("
          //     << (*j)[0] << "," << (*j)[1] << "," << (*j)[2] << ")" << endl;
Francois Gygi committed
144 145
        }
      }
146
    }
Francois Gygi committed
147
  }
148
#endif
Francois Gygi committed
149

150 151
#if 1
  // remove equivalent points
Francois Gygi committed
152 153
  for (list<vector<int> >::iterator i = kplist.begin(); i != kplist.end(); i++)
  {
154 155 156 157 158
    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
    for ( list<vector<int> >::iterator j = i; j != kplist.end(); j++ )
Francois Gygi committed
159
    {
160
      if ( j != i )
Francois Gygi committed
161
      {
162 163 164 165 166 167 168 169
        D3vector kj = (*j)[0]/(2.0*nx) * b0 +
                      (*j)[1]/(2.0*ny) * b1 +
                      (*j)[2]/(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;
Francois Gygi committed
170
          //cout << " erasing equivalent point " << "(" << (*j)[0] << ","
171 172 173
          //     << (*j)[1] << "," << (*j)[2] << ") == ("
          //     << (*i)[0] << "," << (*i)[1] << "," << (*i)[2] << ")" << endl;
        }
Francois Gygi committed
174 175 176
      }
    }
  }
177
#endif
Francois Gygi committed
178

179
  // remove elements with zero weight
Francois Gygi committed
180
  for (list<vector<int> >::iterator i = kplist.begin();
181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210
       i != kplist.end(); /* nothing */ )
  {
    if ( (*i)[3] == 0 )
    {
      kplist.erase(i++);
    }
    else
    {
      i++;
    }
  }

#if 1
  // make first index positive
  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;
    }
  }
#endif

#if 1
  // check that the sum of weights is one
Francois Gygi committed
211 212 213 214 215
  int sum = 0;
  for (list<vector<int> >::iterator i = kplist.begin(); i != kplist.end(); i++)
  {
    sum += (*i)[3];
  }
216 217
  assert(sum==total_weight);
#endif
Francois Gygi committed
218

219
#if 1
Francois Gygi committed
220 221 222 223
  // output list
  // traverse list  backwards to have increasing indices
  // kpoints are output in reciprocal lattice coordinates
  cout.setf(ios::right,ios::adjustfield);
224 225 226 227 228 229 230 231 232
  cout << "# nx,ny,nz: " << nx << " " << ny << " " << nz << endl;
  cout << "# sx,sy,sz: " << sx << " " << sy << " " << sz << endl;
  cout << "# a0: " << a0 << endl;
  cout << "# a1: " << a1 << endl;
  cout << "# a2: " << a2 << endl;
  cout << "# b0/(2pi): " << b0/(2*M_PI) << endl;
  cout << "# b1/(2pi): " << b1/(2*M_PI) << endl;
  cout << "# b2/(2pi): " << b2/(2*M_PI) << endl;

Francois Gygi committed
233 234 235 236 237 238 239 240
  cout << "# " << kplist.size() << " k-points" << endl;
  cout << " kpoint delete 0 0 0" << endl;
  for (list<vector<int> >::reverse_iterator i = kplist.rbegin();
       i != kplist.rend(); i++)
  {
    cout.setf(ios::fixed,ios::floatfield);
    cout << " kpoint add "
         << setprecision(10)
241 242 243
         << setw(13) << ((*i)[0]+sx)/(2.0*nx) << " "
         << setw(13) << ((*i)[1]+sy)/(2.0*ny) << " "
         << setw(13) << ((*i)[2]+sz)/(2.0*nz) << "   ";
Francois Gygi committed
244 245
    cout.setf(ios::scientific,ios::floatfield);
    cout << setprecision(14)
246
         << setw(16) << (*i)[3]/((double) total_weight)
Francois Gygi committed
247 248 249
         << endl;
  }

250 251 252 253 254 255 256 257 258 259 260 261 262
  // output list in absolute coordinates for plot
  ofstream plotfile("kpoint.plt");
  for (list<vector<int> >::reverse_iterator i = kplist.rbegin();
       i != kplist.rend(); i++)
  {
    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;
Francois Gygi committed
263

264 265
  }
#endif
Francois Gygi committed
266
}