Advisory: upstream Henkelman `bader -b weight` stores wrong attractor positions

[janosh]

heads-up for anyone using BaderAnalysis with the Henkelman bader binary and -b weight (the Yu-Trinkle weight method). there's a bug in weight_mod.f90 (v1.05, Aug 2023) that shifts every stored attractor position by one grid step, which can cause basin charges to get assigned to the wrong atom. the default near-grid method is not affected.

the bug

in weight_mod.f90, bader_weight_calc stores the attractor position at line 166:

DO nv = 1, numVect
  p = chgList(n)%pos + vect(nv,:)   ! p overwritten each iteration
  CALL pbc(p, chgref%npts)
  m = indList(p(1), p(2), p(3))
  IF (m < n) THEN
    ...
  END IF
END DO
IF (nabove == 0) THEN ! maxima
  ...
  bdr%volnum(chgList(n)%pos(1), chgList(n)%pos(2), chgList(n)%pos(3)) = bdr%nvols  ! correct position
  bdr%volpos_lat(bdr%bnum,:) = REAL(p,q2)                                           ! leaked p from loop

after the WS neighbour loop, p holds chgList(n)%pos + vect(numVect,:) — the last neighbour's grid position, not the attractor itself. the line should read REAL(chgList(n)%pos, q2).

reasons this looks unintentional:

1. the line directly above (bdr%volnum) correctly uses chgList(n)%pos
2. the non-weight bader_calc (line 176 of bader_mod.f90) uses the same REAL(p,q2) pattern, but there p is the gradient-ascent terminus — the actual attractor position
3. there's no physical reason to shift attractors by the last WS Voronoi neighbour offset

impact

volpos_lat feeds into assign_chg2atom, which decides which atom gets each attractor's charge. the one-grid-step shift usually doesn't matter (nearest atom stays the same), but for interstitial attractors equidistant between two same-element atoms it can flip the assignment.

i ran both bader (default) and bader -b weight on 1000 randomly-sampled DFT CHGCARs. 36 cases disagreed by >0.001 electrons. correcting the attractor position fixed 26 of 36 to floating-point noise (max|Δ| < 1e-6). the worst case (Be₄AlNi, 6 atoms, 72³ grid) went from a 2.39 e discrepancy to 5e-7. errors always show up as anti-correlated swaps between same-element atoms — total electrons and volume are conserved exactly. most common with light/diffuse atoms (Be, B, Si, Ca) where interstitial attractors sit between symmetry-equivalent sites.

reproduction

2 Si atoms at frac x=0.25 and x=0.75 in a 6 Å cubic cell, Gaussian peaks on each atom plus a broad peak at the midpoint. grid 24×8×8. running both methods on the same file:

$ bader CHGCAR

BCF.dat:
    #         X        CHARGE     ATOM    DISTANCE
    1      4.5000      0.2121       2      0.0000
    2      1.5000      0.2121       1      0.0000
    3      3.0000      0.2497       1      1.5000   ← correct midpoint position

ACF.dat:
    1    atom1    0.461865
    2    atom2    0.212132

$ bader CHGCAR -b weight

BCF.dat:
    1      4.7500      0.2121       2      0.2500
    2      1.7500      0.2121       1      0.2500
    3      3.2500      0.2498       2      1.2500   ← shifted +0.25 Å, assigned to wrong atom

ACF.dat:
    1    atom1    0.212101   ← charges swapped
    2    atom2    0.461896

the default method puts attractor #3 at x=3.0 (the actual density maximum) and assigns it to atom 1. the weight method puts it at x=3.25 (shifted +1 grid step), making it closer to atom 2, so the basin charge goes to the wrong atom. total electrons are identical (0.674).

python script to generate the CHGCAR

import math

a = 6.0
ngx, ngy, ngz = 24, 8, 8
volume = a**3
atom1 = [0.25, 0.5, 0.5]
atom2 = [0.75, 0.5, 0.5]

def gauss(fx, fy, fz, centers, sigma):
    inv2s2 = 1.0 / (2.0 * sigma * sigma)
    cx, cy, cz = a * fx, a * fy, a * fz
    total = 0.0
    for c in centers:
        ax, ay, az = a * c[0], a * c[1], a * c[2]
        for d1 in (-1, 0, 1):
            for d2 in (-1, 0, 1):
                for d3 in (-1, 0, 1):
                    dx = cx - ax + d1 * a
                    dy = cy - ay + d2 * a
                    dz = cz - az + d3 * a
                    total += math.exp(-(dx*dx + dy*dy + dz*dz) * inv2s2)
    return total

with open("CHGCAR", "w") as f:
    f.write("bug repro\n   1.00000000000000\n")
    for row in [[a,0,0],[0,a,0],[0,0,a]]:
        f.write(f"   {row[0]:12.6f}{row[1]:12.6f}{row[2]:12.6f}\n")
    f.write("   Si\n     2\nDirect\n")
    f.write(f"  {atom1[0]:.8f}  {atom1[1]:.8f}  {atom1[2]:.8f}\n")
    f.write(f"  {atom2[0]:.8f}  {atom2[1]:.8f}  {atom2[2]:.8f}\n")
    f.write(f"\n   {ngx}   {ngy}   {ngz}\n")
    for iz in range(ngz):
        for iy in range(ngy):
            for ix in range(ngx):
                fx, fy, fz = ix/ngx, iy/ngy, iz/ngz
                rho = gauss(fx, fy, fz, [atom1, atom2], 0.15) * volume \
                    + gauss(fx, fy, fz, [[0.5, 0.5, 0.5]], 0.3) * volume * 0.5
                f.write(f" {rho:17.11E}")
            f.write("\n")

possibly related

the 0.05 e residual asymmetry reported in 2017 for symmetric C₃B with -b weight may be partly caused by this — the structure has many symmetry-equivalent C atoms where interstitial attractor mis-assignment would produce exactly that kind of asymmetry.

upstream

planning to also report this on the Henkelman bader forum, but my account is still pending approval