[gpaw-users] non-orthorhombic bug in grid spacing

[Hummelshoej, Jens Strabo]

Hi everybody,
I hope I don't break to many hearts, but there is a pretty serious bug in GPAW when used for non-orthorhombic cells.
If you have used non-orthorhombic cells and specified for instance h=0.18, you could have h=0.25 (or any grid spacing) in some directions even though the txt output says e.g. 0.1723. You can check this by calculating the grid-spacing manually from the vector and number of grid-points.

The problem is that the sub-rountine for calculating vector lengths in the grid spacing methods only works for orthorhombic cells.

in 
gpaw/grid_descriptor.py:

def get_grid_spacings(self):
	L_c = (np.linalg.inv(self.cell_cv)**2).sum(0)**-0.5                                                                                   
       return L_c / self.N_c

and in
gpaw/utitlites/__init.py__

def h2gpts(h, cell_cv, idiv=4):
   """Convert grid spacing to number of grid points divisible by idiv.                                                                        

   Note that units of h and cell_cv must match!                                                                                               

   h: float                                                                                                                                   
       Desired grid spacing in.                                                                                                               
   cell_cv: 3x3 ndarray                                                                                                                       
       Unit cell.                                                                                                                             
   """

   L_c = (np.linalg.inv(cell_cv)**2).sum(0)**-0.5                                                                                            
   return np.maximum(idiv, (L_c / h / idiv + 0.5).astype(int) * idiv)



Using something like:

L_c = np.array([np.linalg.norm(v) for v in cell_cv])

does the trick.


I noticed the problem, since I got a ~13eV difference for a calculation of the exact same system, in two different but compatible cells. 
(in the second I simply set cell[0] = cell[0] + cell[1])

I have included the cells from the txt output of the two calculations.
The problem is visible when you compare the grid-spacing for 2. axis:

          Periodic     X           Y           Z      Points  Spacing
 --------------------------------------------------------------------
 1. axis:    yes    7.921364  -12.632842    3.125818    44     0.1800
 2. axis:    yes    0.000000   11.730496    0.000000    36     0.1731
 3. axis:    yes    0.000000    0.000000   20.000000   104     0.1789



          Periodic     X           Y           Z      Points  Spacing
 --------------------------------------------------------------------
 1. axis:    yes    7.921364   -0.902346    3.125818    44     0.1800
 2. axis:    yes    0.000000   11.730496    0.000000    64     0.1821
 3. axis:    yes    0.000000    0.000000   20.000000   104     0.1789



Jens Strabo Hummelshoej
Post doc @ SUNCAT
SLAC National Accelerator Laboratory
Building 40A, Room 145
2575 Sand Hill Road MS 59
Menlo Park, CA 94025
Phone: 650-926-5017
Fax:   650-926-3712
Cell: 650-773-8480
strabo at slac.stanford.edu