kpgen.cpp 8.51 KB
 Francois Gygi committed Sep 11, 2017 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  Francois Gygi committed Oct 13, 2017 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.  Francois Gygi committed Sep 11, 2017 18 19 20 21 22 23 24 25 26 27 // //////////////////////////////////////////////////////////////////////////////// #include #include #include #include #include #include "D3vector.h" using namespace std;  Francois Gygi committed Oct 13, 2017 28 29 30 31 32 const char *const version = "2.1"; // largest shift when scanning reciprocal space const int shmax=2;  Francois Gygi committed Sep 11, 2017 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) {  Francois Gygi committed Oct 13, 2017 38 39 40  for ( int i0 = -shmax; i0 <= shmax; i0++ ) for ( int i1 = -shmax; i1 <= shmax; i1++ ) for ( int i2 = -shmax; i2 <= shmax; i2++ )  Francois Gygi committed Sep 11, 2017 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  Francois Gygi committed Oct 13, 2017 53  // in the first shmax shells  Francois Gygi committed Sep 11, 2017 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);  Francois Gygi committed Oct 13, 2017 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; }  Francois Gygi committed Sep 11, 2017 67 68 69 70 71 72  return true; } //////////////////////////////////////////////////////////////////////////////// int main(int argc, char** argv) {  Francois Gygi committed Oct 13, 2017 73  cout << "# kpgen " << version << endl;  Francois Gygi committed Sep 11, 2017 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;  Francois Gygi committed Oct 13, 2017 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; } }  Francois Gygi committed Sep 11, 2017 128 129 130 131 132  vector kp; vector kpfrac; vector weight; // scan volume enclosing the BZ  Francois Gygi committed Oct 13, 2017 133 134 135  for ( int ii = -shmax; ii <= shmax; ii++ ) for ( int jj = -shmax; jj <= shmax; jj++ ) for ( int kk = -shmax; kk <= shmax; kk++ )  Francois Gygi committed Sep 11, 2017 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;  Francois Gygi committed Oct 13, 2017 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);  Francois Gygi committed Sep 11, 2017 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);  Francois Gygi committed Oct 13, 2017 166 #if 0  Francois Gygi committed Sep 11, 2017 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;  Francois Gygi committed Oct 13, 2017 181  //return 1;  Francois Gygi committed Sep 11, 2017 182  }  Francois Gygi committed Oct 13, 2017 183 #endif  Francois Gygi committed Sep 11, 2017 184 185  // count vectors that are equivalent to -k+G  Francois Gygi committed Oct 13, 2017 186 #if 0  Francois Gygi committed Sep 11, 2017 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  Francois Gygi committed Oct 13, 2017 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  Francois Gygi committed Sep 11, 2017 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; } } } }  Francois Gygi committed Oct 13, 2017 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++;  Francois Gygi committed Sep 11, 2017 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;  Francois Gygi committed Oct 13, 2017 249  cout << "# " << nkp << " k-points" << endl;  Francois Gygi committed Sep 11, 2017 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  Francois Gygi committed Oct 13, 2017 264 #if 1  Francois Gygi committed Sep 11, 2017 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; }  Francois Gygi committed Oct 13, 2017 283 #endif  Francois Gygi committed Sep 11, 2017 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; } } }