diff --git a/util/kpgen/kpgen.C b/util/kpgen/kpgen.C
index eaf7881..bd87404 100644
--- a/util/kpgen/kpgen.C
+++ b/util/kpgen/kpgen.C
@@ -7,34 +7,58 @@
 //       a2x a2y a2z = third basis vector of the unit cell
 //
 //       nx,ny,nz: number of kpoints in each direction
-//       sx,sy,sz: shift in each direction
+//       sx,sy,sz: shift in each direction (floating point)
 //       shift: 0: symmetric set: boundary points not included
-//              1: shifted set: boundary points 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
+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;
+}
+
 int main(int argc, char** argv)
 {
+  cout << "# kpgen-1.0" << endl;
   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;
-    cerr << "      shift==1: shifted set, zone boundary included" << endl;
     return 1;
   }
   int nx = atoi(argv[1]);
   int ny = atoi(argv[2]);
   int nz = atoi(argv[3]);
-  int sx = atoi(argv[4]);
-  int sy = atoi(argv[5]);
-  int sz = atoi(argv[6]);
+  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}"
@@ -48,7 +72,7 @@ int main(int argc, char** argv)
   {
     cerr << " use: " << argv[0] << " nx ny nz shiftx shifty shiftz {cell}"
          << endl;
-    cerr << " shifts must be 0 or 1" << endl;
+    cerr << " shifts must be in [0,1]" << endl;
     return 1;
   }
   D3vector a0(atof(argv[7]),atof(argv[8]),atof(argv[9]));
@@ -67,102 +91,145 @@ int main(int argc, char** argv)
 
   list<vector<int> > kplist;
   vector<int> kpint(4);
-  for ( int i = nx-1; i >= 0; i-- )
+
+  // 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++ )
   {
-    for ( int j = ny-1; j >= 0; j-- )
+    for ( int j = 0; j < ny; j++ )
     {
-      for ( int k = nz-1; k >= 0; k-- )
+      for ( int k = 0; k < nz; k++ )
       {
-        kpint[0] = 2*i-nx+sx+1;
-        kpint[1] = 2*j-ny+sy+1;
-        kpint[2] = 2*k-nz+sz+1;
+        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;
-        kplist.push_back(kpint);
+
+        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);
       }
     }
   }
 
-  // fold kpoints back in the BZ
-  for (list<vector<int> >::iterator i = kplist.begin();
-       i != kplist.end(); i++)
+  int total_weight = kplist.size();
+
+#if 1
+  // remove -k
+  for (list<vector<int> >::iterator i = kplist.begin(); i != kplist.end(); i++)
   {
+    // 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];
-    bool done = false;
-    while (!done)
+    if ( kpint[0]*kpint[0]+kpint[1]*kpint[1]+kpint[2]*kpint[2] != 0 )
     {
-      done = true;
-      D3vector kabs = kpint[0]/(2.0*nx) * b0 +
-                      kpint[1]/(2.0*ny) * b1 +
-                      kpint[2]/(2.0*nz) * b2;
-      D3vector g;
-      // check projection of kpoint along all 26 reciprocal lattice vectors
-      // that are nearest g=0
-
-      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) )
+      // look for -k in the rest of the list
+      for ( list<vector<int> >::iterator j = i; j != kplist.end(); j++ )
       {
-        D3vector g = i0 * b0 + i1 * b1 + i2 * b2;
-        if ( kabs*g >  (0.5 + 1.e-4) * g*g )
+        if ( (*j)[0]==-kpint[0] && (*j)[1]==-kpint[1] && (*j)[2]==-kpint[2] )
         {
-          kpint[0]-= i0*2*nx;
-          kpint[1]-= i1*2*ny;
-          kpint[2]-= i2*2*nz;
-          kabs = kpint[0]/(2.0*nx) * b0 +
-                 kpint[1]/(2.0*ny) * b1 +
-                 kpint[2]/(2.0*nz) * b2;
-          done = false;
+          // transfer weight to (*i)
+          (*i)[3] += (*j)[3];
+          (*j)[3] = 0;
+          //cout << " erasing  " << "(" << kpint[0] << ","
+          //     << kpint[1] << "," << kpint[2] << ") == -("
+          //     << (*j)[0] << "," << (*j)[1] << "," << (*j)[2] << ")" << endl;
         }
       }
-    } // while !done
-    // kpint is now inside the BZ
-    (*i)[0] = kpint[0];
-    (*i)[1] = kpint[1];
-    (*i)[2] = kpint[2];
+    }
   }
+#endif
 
-  // remove -k
+#if 1
+  // remove equivalent points
   for (list<vector<int> >::iterator i = kplist.begin(); i != kplist.end(); i++)
   {
-    // 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];
-    if ( kpint[0]*kpint[0]+kpint[1]*kpint[1]+kpint[2]*kpint[2] != 0 )
+    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++ )
     {
-      // look for -k in the rest of the list
-      list<vector<int> >::iterator j = find(i,kplist.end(),kpint);
-      if ( j != kplist.end() )
+      if ( j != i )
       {
-        kplist.erase(j);
-        // double the weight of kpint
-        (*i)[3] *= 2;
+        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;
+          //cout << " erasing equivalent point " << "(" << (*j)[0] << ","
+          //     << (*j)[1] << "," << (*j)[2] << ") == ("
+          //     << (*i)[0] << "," << (*i)[1] << "," << (*i)[2] << ")" << endl;
+        }
       }
     }
   }
+#endif
+
+  // remove elements with zero weight
+  for (list<vector<int> >::iterator i = kplist.begin();
+       i != kplist.end(); /* nothing */ )
+  {
+    if ( (*i)[3] == 0 )
+    {
+      kplist.erase(i++);
+    }
+    else
+    {
+      i++;
+    }
+  }
 
-  // check that sum of weights is one
+#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
   int sum = 0;
   for (list<vector<int> >::iterator i = kplist.begin(); i != kplist.end(); i++)
   {
     sum += (*i)[3];
   }
-  assert(sum==nx*ny*nz);
+  assert(sum==total_weight);
+#endif
 
+#if 1
   // output list
   // traverse list  backwards to have increasing indices
   // kpoints are output in reciprocal lattice coordinates
   cout.setf(ios::right,ios::adjustfield);
-  cout << "# kpgen     " << 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 << "# 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;
+
   cout << "# " << kplist.size() << " k-points" << endl;
   cout << " kpoint delete 0 0 0" << endl;
   for (list<vector<int> >::reverse_iterator i = kplist.rbegin();
@@ -171,13 +238,29 @@ int main(int argc, char** argv)
     cout.setf(ios::fixed,ios::floatfield);
     cout << " kpoint add "
          << setprecision(10)
-         << setw(13) << (*i)[0]/(2.0*nx) << " "
-         << setw(13) << (*i)[1]/(2.0*ny) << " "
-         << setw(13) << (*i)[2]/(2.0*nz) << "   ";
+         << setw(13) << ((*i)[0]+sx)/(2.0*nx) << " "
+         << setw(13) << ((*i)[1]+sy)/(2.0*ny) << " "
+         << setw(13) << ((*i)[2]+sz)/(2.0*nz) << "   ";
     cout.setf(ios::scientific,ios::floatfield);
     cout << setprecision(14)
-         << setw(16) << (*i)[3]/((double) nx*ny*nz)
+         << setw(16) << (*i)[3]/((double) total_weight)
          << endl;
   }
 
+  // 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;
+
+  }
+#endif
 }