kpgen.cpp 8.51 KB
Newer Older
1 2 3 4 5 6 7 8 9 10
////////////////////////////////////////////////////////////////////////////////
//
// 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
//       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
11 12 13 14 15 16 17
//       sx,sy,sz: shift in each direction (floating point in [0,1) )
//
// note: from version 2.1 on, sx=sy=sz=0 generates a symmetric set
// A symmetric set includes the Gamma (k=0) point
//
// note: using non-zero shifts with non-rectangular cells breaks symmetry
// non-zero shifts should only be used with rectangular cells.
18 19 20 21 22 23 24 25 26 27
//
////////////////////////////////////////////////////////////////////////////////
#include<iostream>
#include<fstream>
#include<iomanip>
#include<vector>
#include<cstdlib>
#include "D3vector.h"
using namespace std;

28 29 30 31 32
const char *const version = "2.1";

// largest shift when scanning reciprocal space
const int shmax=2;

33 34 35 36 37
////////////////////////////////////////////////////////////////////////////////
// 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)
{
38 39 40
  for ( int i0 = -shmax; i0 <= shmax; i0++ )
    for ( int i1 = -shmax; i1 <= shmax; i1++ )
      for ( int i2 = -shmax; i2 <= shmax; i2++ )
41 42 43 44 45 46 47 48 49 50 51 52
      {
        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)
{
  // check projection of kpoint k along all reciprocal lattice vectors
53
  // in the first shmax shells
54 55 56 57
  // 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);
58 59 60 61 62 63 64 65 66
  for ( int i0 = -shmax; i0 <= shmax; i0++ )
    for ( int i1 = -shmax; i1 <= shmax; i1++ )
      for ( int i2 = -shmax; i2 <= shmax; i2++ )
        if ( !((i0 == 0) && (i1 == 0) && (i2 == 0)) )
        {
          D3vector g = i0 * b0 + i1 * b1 + i2 * b2;
          if ( kshifted*g >  0.5 * g*g )
            return false;
        }
67 68 69 70 71 72
  return true;
}

////////////////////////////////////////////////////////////////////////////////
int main(int argc, char** argv)
{
73
  cout << "# kpgen " << version << endl;
74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115
  if ( argc != 16 )
  {
    cerr << " use: " << argv[0] << " nx ny nz shiftx shifty shiftz {cell}"
         << endl;
    return 1;
  }
  int nx = atoi(argv[1]);
  int ny = atoi(argv[2]);
  int nz = atoi(argv[3]);
  double sx = atof(argv[4]);
  double sy = atof(argv[5]);
  double sz = atof(argv[6]);
  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;
    cerr << " shifts must be in [0,1]" << endl;
    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;

116 117 118 119 120 121 122 123 124 125 126 127
  // check if shifts are used with a non-rectangular cell
  if ( sx > 0 || sy > 0 || sz > 0 )
  {
    if ( (fabs(a0*a1) > 1.e-6) ||
         (fabs(a1*a2) > 1.e-6) ||
         (fabs(a0*a2) > 1.e-6) )
    {
      cout << " warning: non-zero shifts with non-rectangular cell"
           << " may break symmetry" << endl;
    }
  }

128 129 130 131 132
  vector<D3vector> kp;
  vector<D3vector> kpfrac;
  vector<double> weight;

  // scan volume enclosing the BZ
133 134 135
  for ( int ii = -shmax; ii <= shmax; ii++ )
  for ( int jj = -shmax; jj <= shmax; jj++ )
  for ( int kk = -shmax; kk <= shmax; kk++ )
136 137 138 139 140 141 142 143 144 145
  for ( int i = 0; i < nx; i++ )
  {
    for ( int j = 0; j < ny; j++ )
    {
      for ( int k = 0; k < nz; k++ )
      {
        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;

146 147 148
        double kv0 = (kpint0 + sx + (nx%2+1))/(2.0*nx);
        double kv1 = (kpint1 + sy + (ny%2+1))/(2.0*ny);
        double kv2 = (kpint2 + sz + (nz%2+1))/(2.0*nz);
149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165

        D3vector kv = kv0*b0 + kv1*b1 + kv2*b2;

        if ( in_BZ(kv,b0,b1,b2) )
        {
          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);
166
#if 0
167 168 169 170 171 172 173 174 175 176 177 178 179 180
  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;

      }
  if ( nequiv != 0 )
  {
    // there should not be any equivalent points as k=k+G
    cerr << nequiv << " error: equivalent points (k=k+G)" << endl;
181
    //return 1;
182
  }
183
#endif
184 185

  // count vectors that are equivalent to -k+G
186
#if 0
187 188 189 190 191 192 193 194 195 196 197 198
  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) )
      {
        nequivm++;
        //cout << setprecision(3)
        //     << kpfrac[i] << " equiv " << kpfrac[j] << endl;
      }
  cout << nequivm << " equivalent points (k=-k+G)" << endl;
#endif

199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214
  // check duplicates
#if 0
  int ndup = 0;
  for ( int i = 0; i < kp.size(); i++ )
    for ( int j = i+1; j < kp.size(); j++ )
      if ( length(kp[i]-kp[j]) < 1.e-5 )
      {
        ndup++;
        cout << setprecision(3)
             << kpfrac[i] << " duplicate of " << kpfrac[j] << endl;
      }
  cout << "# " << ndup << " duplicates" << endl;
#endif

#if 1
  // reassign weight from (k,-k+G) equivalent points
215 216 217 218 219 220 221 222 223 224 225 226 227 228 229
  for ( int i = 0; i < kp.size(); i++ )
  {
    if ( weight[i] != 0.0 )
    {
      for ( int j = i+1; j < kp.size(); j++ )
      {
        if ( equals(kp[i],-kp[j],b0,b1,b2) )
        {
          // reassign the weight of kp[j] to kp[i]
          weight[i] += weight[j];
          weight[j] = 0.0;
        }
      }
    }
  }
230 231 232 233 234 235
#endif

  // count k points with non-zero weight
  int nkp = 0;
  for ( int i = 0; i < weight.size(); i++ )
    if ( weight[i] > 0.0 ) nkp++;
236 237 238 239 240 241 242 243 244 245 246 247 248

  // output list
  // kpoints are output in reciprocal lattice coordinates
  cout.setf(ios::right,ios::adjustfield);
  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;

249
  cout << "# " << nkp << " k-points" << endl;
250 251 252 253 254 255 256 257 258 259 260 261 262 263
  cout << " kpoint delete 0 0 0" << endl;

  // print list backward to have increasing x components
  for ( int i = kpfrac.size()-1; i >= 0; i-- )
  {
    if ( weight[i] != 0.0 )
    {
      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
264
#if 1
265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282
      if ( kx == 0.0 )
      {
        if ( ky == 0.0 )
        {
          if ( kz < 0.0 ) kz = -kz;
        }
        else if ( ky < 0.0 )
        {
          ky = -ky;
          if ( kz != 0 ) kz = -kz;
        }
      }
      else if ( kx < 0.0 )
      {
        kx = -kx;
        if ( ky != 0 ) ky = -ky;
        if ( kz != 0 ) kz = -kz;
      }
283
#endif
284 285 286 287 288 289 290 291 292 293 294 295 296
      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;
    }
  }
}