Advisory: possible upstream Henkelman bader -b weight WS Voronoi stencil bug on non-orthogonal cells

[janosh]

a possible second issue in henkelman bader (weight_mod.f90, v1.05), separate from #4638. flagging it for visibility — the asymmetric stencil is used consistently throughout -b weight, so impact on per-atom charges may cancel in ways i haven't traced. cubic and orthorhombic cells are unaffected by construction; the open question is whether non-orthogonal cells (hexagonal, monoclinic, triclinic) get biased basin partitions.

bader -b weight builds its WS Voronoi neighbour stencil from a voxel cell matrix that, under fortran MATMUL semantics, behaves as if transposed relative to the convention used by chgcar_mod.f90. on non-orthogonal cells the resulting flux prefactors α describe a different geometry than the actual voxel grid.

the observation

weight_mod.f90 lines 61-65:

DO i=1,3
  cell(i,:) = ions%lattice(i,:)/chgval%npts(i)   ! row i = lattice vector i / npts(i)
END DO
CALL ws_voronoi(cell, numVect, vect, alpha)

then line 349:

R(neigh,:) = MATMUL(cell, nv)

fortran MATMUL(rank-2, rank-1) is R(i) = Σ_k cell(i,k) * nv(k). with rows = δ_i this gives R(i) = δ_i · nv — not the cartesian translation m·δ_1 + n·δ_2 + p·δ_3 (only equal when cell is diagonal). R then drives the entire WS construction: incell half-plane test, vertex polygon, facet area integral.

symmetry consequence

the six in-plane neighbours of a hexagonal voxel cell (±1,0,0), (0,±1,0), (±1,±1,0) are equidistant; their α = A_facet / |R| must therefore be identical. under the row convention they aren't — 3 distinct |R| and α values, α varying by 96%:

hexagonal voxel lattice (a = 0.5675 Å, c = 0.8975 Å):

  rows convention:                   columns convention:
    offset      |R|        α            offset      |R|        α
   (-1,-1,0)   0.6043   0.4488         (-1,-1,0)   0.5675   0.5182
   (-1, 0,0)   0.6345   0.3285         (-1, 0,0)   0.5675   0.5182
   ( 0,-1,0)   0.4915   0.8467         ( 0,-1,0)   0.5675   0.5182
   ( 0, 1,0)   0.4915   0.8467         ( 0, 1,0)   0.5675   0.5182
   ( 1, 0,0)   0.6345   0.3285         ( 1, 0,0)   0.5675   0.5182
   ( 1, 1,0)   0.6043   0.4488         ( 1, 1,0)   0.5675   0.5182
  α relative spread: 95.7%             α relative spread:  0.0%

internal inconsistency

chgcar_mod.f90 line 96 builds the equivalent voxel→cartesian matrix the other way:

chg%lat2car(:,i) = ions%dir2car(:,i)/REAL(chg%npts(i),q2)   ! columns

chg%lat2car yields correct cartesian translations under MATMUL and is used correctly elsewhere (voronoi_mod.f90, atom assignment, etc.). only weight_mod.f90's cell uses the row convention — and it's the only matrix feeding ws_voronoi. could be deliberate, could be a slip; can't tell from the source.

the alternative

if unintentional, the one-line change to align with chgcar_mod.f90's lat2car:

 DO i=1,3
-  cell(i,:) = ions%lattice(i,:)/chgval%npts(i)
+  cell(:,i) = ions%lattice(i,:)/chgval%npts(i)
 END DO

summary: unsure whether the symmetry violation propagates into measurably different per-atom charges, or cancels somewhere downstream. the default near-grid method doesn't use ws_voronoi and is unaffected regardless.