kpgen.C 7.6 KB
 Francois Gygi committed Nov 30, 2009 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  Francois Gygi committed Mar 29, 2013 10 // sx,sy,sz: shift in each direction (floating point)  Francois Gygi committed Nov 30, 2009 11 // shift: 0: symmetric set: boundary points not included  Francois Gygi committed Mar 29, 2013 12 13 // even-numbered sets do not include the gamma point // odd-numbered sets include the gamma point  Francois Gygi committed Nov 30, 2009 14 15 // #include  Francois Gygi committed Mar 29, 2013 16 #include  Francois Gygi committed Nov 30, 2009 17 18 19 #include #include #include  Francois Gygi committed Apr 01, 2013 20 #include  Francois Gygi committed Nov 30, 2009 21 22 23 24 #include #include "D3vector.h" using namespace std;  Francois Gygi committed Mar 29, 2013 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 Nov 30, 2009 46 47 int main(int argc, char** argv) {  Francois Gygi committed Mar 29, 2013 48  cout << "# kpgen-1.0" << endl;  Francois Gygi committed Nov 30, 2009 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]);  Francois Gygi committed Mar 29, 2013 59 60 61  double sx = atof(argv[4]); double sy = atof(argv[5]); double sz = atof(argv[6]);  Francois Gygi committed Nov 30, 2009 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;  Francois Gygi committed Mar 29, 2013 75  cerr << " shifts must be in [0,1]" << endl;  Francois Gygi committed Nov 30, 2009 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 > kplist; vector kpint(4);  Francois Gygi committed Mar 29, 2013 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 Nov 30, 2009 100  {  Francois Gygi committed Mar 29, 2013 101  for ( int j = 0; j < ny; j++ )  Francois Gygi committed Nov 30, 2009 102  {  Francois Gygi committed Mar 29, 2013 103  for ( int k = 0; k < nz; k++ )  Francois Gygi committed Nov 30, 2009 104  {  Francois Gygi committed Mar 29, 2013 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 Nov 30, 2009 108  kpint[3] = 1;  Francois Gygi committed Mar 29, 2013 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 Nov 30, 2009 116 117 118 119  } } }  Francois Gygi committed Mar 29, 2013 120 121 122 123 124  int total_weight = kplist.size(); #if 1 // remove -k for (list >::iterator i = kplist.begin(); i != kplist.end(); i++)  Francois Gygi committed Nov 30, 2009 125  {  Francois Gygi committed Mar 29, 2013 126  // test if -k is in the set (and is not 0 0 0)  Francois Gygi committed Nov 30, 2009 127 128 129 130  kpint[0] = (*i)[0]; kpint[1] = (*i)[1]; kpint[2] = (*i)[2]; kpint[3] = (*i)[3];  Francois Gygi committed Mar 29, 2013 131  if ( kpint[0]*kpint[0]+kpint[1]*kpint[1]+kpint[2]*kpint[2] != 0 )  Francois Gygi committed Nov 30, 2009 132  {  Francois Gygi committed Mar 29, 2013 133 134  // look for -k in the rest of the list for ( list >::iterator j = i; j != kplist.end(); j++ )  Francois Gygi committed Nov 30, 2009 135  {  Francois Gygi committed Mar 29, 2013 136  if ( (*j)[0]==-kpint[0] && (*j)[1]==-kpint[1] && (*j)[2]==-kpint[2] )  Francois Gygi committed Nov 30, 2009 137  {  Francois Gygi committed Mar 29, 2013 138 139 140  // transfer weight to (*i) (*i)[3] += (*j)[3]; (*j)[3] = 0;  Francois Gygi committed Apr 16, 2013 141  //cout << " erasing " << "(" << kpint[0] << ","  Francois Gygi committed Mar 29, 2013 142 143  // << kpint[1] << "," << kpint[2] << ") == -(" // << (*j)[0] << "," << (*j)[1] << "," << (*j)[2] << ")" << endl;  Francois Gygi committed Nov 30, 2009 144 145  } }  Francois Gygi committed Mar 29, 2013 146  }  Francois Gygi committed Nov 30, 2009 147  }  Francois Gygi committed Mar 29, 2013 148 #endif  Francois Gygi committed Nov 30, 2009 149   Francois Gygi committed Mar 29, 2013 150 151 #if 1 // remove equivalent points  Francois Gygi committed Nov 30, 2009 152 153  for (list >::iterator i = kplist.begin(); i != kplist.end(); i++) {  Francois Gygi committed Mar 29, 2013 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 >::iterator j = i; j != kplist.end(); j++ )  Francois Gygi committed Nov 30, 2009 159  {  Francois Gygi committed Mar 29, 2013 160  if ( j != i )  Francois Gygi committed Nov 30, 2009 161  {  Francois Gygi committed Mar 29, 2013 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 Apr 16, 2013 170  //cout << " erasing equivalent point " << "(" << (*j)[0] << ","  Francois Gygi committed Mar 29, 2013 171 172 173  // << (*j)[1] << "," << (*j)[2] << ") == (" // << (*i)[0] << "," << (*i)[1] << "," << (*i)[2] << ")" << endl; }  Francois Gygi committed Nov 30, 2009 174 175 176  } } }  Francois Gygi committed Mar 29, 2013 177 #endif  Francois Gygi committed Nov 30, 2009 178   Francois Gygi committed Mar 29, 2013 179  // remove elements with zero weight  Francois Gygi committed Apr 16, 2013 180  for (list >::iterator i = kplist.begin();  Francois Gygi committed Mar 29, 2013 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 >::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 Nov 30, 2009 211 212 213 214 215  int sum = 0; for (list >::iterator i = kplist.begin(); i != kplist.end(); i++) { sum += (*i)[3]; }  Francois Gygi committed Mar 29, 2013 216 217  assert(sum==total_weight); #endif  Francois Gygi committed Nov 30, 2009 218   Francois Gygi committed Mar 29, 2013 219 #if 1  Francois Gygi committed Nov 30, 2009 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);  Francois Gygi committed Mar 29, 2013 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 Nov 30, 2009 233 234 235 236 237 238 239 240  cout << "# " << kplist.size() << " k-points" << endl; cout << " kpoint delete 0 0 0" << endl; for (list >::reverse_iterator i = kplist.rbegin(); i != kplist.rend(); i++) { cout.setf(ios::fixed,ios::floatfield); cout << " kpoint add " << setprecision(10)  Francois Gygi committed Mar 29, 2013 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 Nov 30, 2009 244 245  cout.setf(ios::scientific,ios::floatfield); cout << setprecision(14)  Francois Gygi committed Mar 29, 2013 246  << setw(16) << (*i)[3]/((double) total_weight)  Francois Gygi committed Nov 30, 2009 247 248 249  << endl; }  Francois Gygi committed Mar 29, 2013 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 >::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 Apr 16, 2013 263   Francois Gygi committed Mar 29, 2013 264 265  } #endif  Francois Gygi committed Nov 30, 2009 266 }