pyMBE-dev / pyMBE

pyMBE provides tools to facilitate building up molecules with complex architectures in the Molecular Dynamics software ESPResSo. For an up-to-date API documention please check our website:
https://pymbe-dev.github.io/pyMBE/pyMBE.html
GNU General Public License v3.0
6 stars 8 forks source link

Use 2019 SI base units #86

Closed jngrad closed 2 weeks ago

jngrad commented 2 weeks ago

The SI base units have been redefined as natural constants with exact precision in 2019.

jngrad commented 2 weeks ago

The old units were experimentally determined and had fluctuating precision. Each textbook had a slightly different value.

We are currently using the old values, and they have a tiny deviation from the new values. For example, the ratio of pyMBE Boltzmann constant and the corresponding 2019 constant deviates from unity by three times the machine-precision of single-precision floating point values:

import pyMBE, numpy as np, scipy.constants
eps = scipy.constants.k / pyMBE.pymbe_library.Kb.magnitude - 1.
print(eps / np.finfo(np.float32).eps)

For the electric charge, it's one order of magnitude smaller than machine precision.

I don't think it makes much of a difference in long simulations due to thermal fluctuations, but it may cause test failures in unit tests if the atol parameter of np.testing.assert_*() is larger than epsilon, and the test uses the new value of kB. The 2019 value for kB is exact and identical on NIST (link), IUPAC (link), BIPM (link), and scipy.

pm-blanco commented 2 weeks ago

@jngrad it looks to good to me. I ran the functional tests and everything seems to work fine. Should we add a small unit tests that checks that pyMBE uses the same fundamental constants as the ones in scipy, just to ensure that this does not break again in the future?

jngrad commented 2 weeks ago

Should we add a small unit tests that checks that pyMBE uses the same fundamental constants as the ones in scipy, just to ensure that this does not break again in the future?

Good point. I didn't realize the values of the constants were overridden in set_reduced_units(), yikes! We really need to remove these global variables :sweat_smile: (#60)

Side-note: I personally find using Kb instead of kB for the Boltzmann constant $k_B$ to be a bit unfortunate considering how similar it is to the equilibrium constant $K_b$, and that we define the self-ionization constant $K_w$ in member Kw.

pm-blanco commented 2 weeks ago

@jngrad

Good point. I didn't realize the values of the constants were overridden in set_reduced_units(), yikes!

Yes, the problem is that in set_reduced_units() one needs to create a new pint.UnitRegistry() object to redefine the reduced unit system. Then, one needs to basically overwrite all previously existing pint variables because it is not possible to do operations between pint objects defined with different unit registries. That is why one needs to overwrite the constants so one can still operate with them with the new unit system.

The Pint documentation is not very extensive regarding this definitions of new units docs. The solution that we have right now is not perfect because a user could still define particles with an old unit registry then set a different set of reduced units and then proceed to define new pyMBE objects with the new unit registry. If we could simply redefine the reduced units without creating a new pint.UnitRegistry() object that would be a much more robust solution.

We really need to remove these global variables 😅 (https://github.com/pyMBE-dev/pyMBE/issues/60)

I agree, we should deal with this issue before the pyMBE 1.0.0 release. There is no major reason for having the global variables and we only need to have them defined by the class constructor instead.

I personally find using Kb instead of kB for the Boltzmann constant k B to be a bit unfortunate considering how similar it is to the equilibrium constant K b , and that we define the self-ionization constant K w in member Kw.

I have not realized of this issue. You are right that it might be a bit missleading. Feel free to change the notation in this PR if you like.

jngrad commented 2 weeks ago

The solution that we have right now is not perfect because a user could still define particles with an old unit registry then set a different set of reduced units and then proceed to define new pyMBE objects with the new unit registry.

We could let the user provide the unit registry as an argument to the constructor. That would require users to write two extra lines to import pint and create their unit registry, but only if they don't use the default reduced units system. That would also make it clear that the unit registry belongs to the user, and pyMBE keeps a reference to it.

This is probably out of scope for this PR, though.

Feel free to change the notation in this PR if you like.

Done.

pm-blanco commented 2 weeks ago

We could let the user provide the unit registry as an argument to the constructor. That would require users to write two extra lines to import pint and create their unit registry, but only if they don't use the default reduced units system. That would also make it clear that the unit registry belongs to the user, and pyMBE keeps a reference to it.

This is probably out of scope for this PR, though.

I am not completely convinced by this solution either, but let's think about it and address it on another PR.