simtk.org

(John Chodera posted this in thread "Units of force constants" -- I'm including it below to start a new thread.)

BY: John Chodera (jchodera)
DATE: 2009-07-04 23:25
SUBJECT: RE: Units of force constants


> I told you getting the units right is the hardest part!

Sherm, this makes me wonder if an MMTK- or Scientific Python-like Units class might make the most sense for OpenMM as well. In the "Feature Request" I submitted that Peter closed (URL below) I outlined what we did for the MMTOOLS SimTK project in Python (which is based on MMTK's concept). This made it *extremely* easy to get the units right, by making sure you are always multiplying and dividing unit-bearing numbers by the appropriate units:

https://simtk.org/tracker/index.php?fun ... 1&atid=436

For example, the call to addBond would look like:

bonds.addBond(i, j, length * Units.Angstrom, k * Units.kcal / Units.mol / Units.Angstrom);

As long as the 'Units' class represents all units in terms of fundamental units of length, mass, and time, then the user never even needs to know what these fundamental units are, and the conversions are done correctly every time. As it is, the lack of such a feature and the incomplete documentation of units in the Doxygen documentation is just causing headaches.

Scientific Python takes this concept yet further, by making unit-bearing numbers special objects. The handling of this subject by both MMTK and Scientific Python are worth looking at, especially if the goal is to facilitate making it "easy to do things right".

John

Getting units right is surprisingly hard to do. In Sid's case he was very diligent about getting the units right in the hard cases yet missed one simple case. That's the kind of error that every human makes no matter how smart. It does seem kind of a shame that we have a nice C++ API but aren't taking advantage of its type safety to avoid such common mistakes. It might be in our own best interest to do something about it -- Sid's experience was no doubt frustrating and we're lucky he kept at it.

I think the only thing that would solve the problem is a system that required units specification at compile time and rigorously type-checked them. The current "optional" units conversion scheme (predefined conversion constants) is better than nothing but fails primarily because it can't force compliance and can't tell a correct unit from a wrong one. Chris Bruns has been exploring compiler-enforced unit correctness for Molmodel. I invited him to join us in this discussion -- maybe we can kill two birds with one stone.

That SciPy units system looks very nice if it is able to enforce correct behavior -- will it fail to run if you forget to specify units? I especially like the fact that you need only to know your own units, not those chosen by the package implementors.

I don't suppose there is much hope for units correctness for Fortran and C users (is there?). But I think we should reconsider units enforcement in C++ now that we've seen what a problem it is in practice. (It is easier said than done, though -- ask Chris!) Certainly for Python we should follow the SciPy precedent, esp. if we can use that units system as-is.

Sherm

Sherm,

These are some excellent points.

Regarding options for unit enforcement in C++: I'm not an expert on the available C++ options by any means, but some googling turned up some interesting possibilities.

* UDUNITS - a C library produced by UCAR that seems to be a "de facto" for unit conversion standard in atmospheric science. Support is also provided in several other languages. It seems pretty heavy-duty, but it would definitely allow for a lot of flexibility in unit specification.

http://www.unidata.ucar.edu/software/udunits/

* UNITS: At the other end of the spectrum is "UNITS, a C++ unit management framework" that is a single header file. Valid unit conversions are enforced at compile-time: A unit-bearing quantity won't allow itself to be cast to a scalar unless the units work out. I can't quite tell how it works from the unit-test examples, but maybe you can figure this out.

http://enbridge.kundert.ca/units/doc/

* C++ Units: A C++ library for unit conversion. Also seems like it might do the trick.

http://calumgrant.net/units/index.html

There are probably other good options out there too.

John