LDAの交換エネルギーの計算でrho^{1/3}が必要になるので
rho = max(fft_grad_l(ig,0),eps_chg)
によりrhoが負になるのを防いでいるが、rhoをゼロでなくeps_chg=1d-25で置き換えているのは、規約密度勾配sの分母にrhoが現れるためである。
しかしrhoが振動し負の値が現れる真空領域は本来はほぼ密度一様かつ勾配ゼロの領域であるが、現状では勾配はそのままの値が用いられている。
その副作用としてrhoが負の位置ではsの分母に微小量eps_chgが現れ発散するが、PBE等ではLieb-Oxford boundによりsの発散はさほど問題にならなかった。
一方vdW-DF(の派生版)ではLieb-Oxford boundが効かない交換汎関数を用いる機会が多いので、sの非物理的な発散は見逃さないのが得策と思われる。