NH3 dipole moment calculation does not match literature value

[annam]

Hello,

I am trying to calculate the dipole moment value for NH3 molecule from the example python scripts https://docs.quantumatk.com/manual/Types/ElectronDensity/ElectronDensity.html#electrondensity-c , but my dipole moment vector seems to be way off - please find my output attached. The literature value for NH3 dipole moment is 4.9*10-30Cm (or 1.47D), but if I take the square root of (mu(x)^2+mu(y)^2 + mu(z)^2), which should give me the dipole magnitude, I get  5.16682E-14 eBohr (or 3.22E+05 Cm or 9.67E+34
D). Could you please help me find an error?

[Anders Blom]

I reran the script, and actually it reports the dipole moment in the wrong unit (Bohr^4*e), so there was a volume normalization gone wrong. This is easy enough to fix (might be version-dependent), but still the numbers are odd.

I get |m1|=2.90861887308e-14*Bohr*e = 2.466e-43 Cm. Note that atkpython is handling units beautifully, just do m1.norm().inUnitsOf(Coulomb*Meter).

But that gives the dipole moment as 7.4e-14 D, so I have to dig a bit deeper

[Anders Blom]

OK, I think I have figured it out. I believe this script/example has been floating around in our documentation for too many years, and at some point we changed the meaning of ElectronDensity from earlier meaning the ElectronDifferenceDensity, which is actually the relevant quantity to use when computing the dipole [https://arxiv.org/ftp/arxiv/papers/1708/1708.03834.pdf]. If I just change the ElectronDifferenceDensity I get m1=1.94D which at least is reasonable, albeit not too accurate (but the molecule was not optimized, and we don't use any fancy functional). This puts the "monopole" to zero (not exactly, because of numerics) which it should be for the difference density, whereas for the density it counts the number of electrons.

[annam]

OK, I think I have figured it out. I believe this script/example has been floating around in our documentation for too many years, and at some point we changed the meaning of ElectronDensity from earlier meaning the ElectronDifferenceDensity, which is actually the relevant quantity to use when computing the dipole [https://arxiv.org/ftp/arxiv/papers/1708/1708.03834.pdf]. If I just change the ElectronDifferenceDensity I get m1=1.94D which at least is reasonable, albeit not too accurate (but the molecule was not optimized, and we don't use any fancy functional). This puts the "monopole" to zero (not exactly, because of numerics) which it should be for the difference density, whereas for the density it counts the number of electrons.

Awesome. This is very helpful. Thank you very much!