Physical constants

Continuing the discussion from Tips to make this toy program faster?:

I suggest that creating a source file with physical constants be discussed in this new thread.

2 Likes

Sorry for derailing the thread, let’s continue the discussion here. I contacted the QCElemental developers and they would be interested in exporting a Fortran header as well from their project. Here is the link again:

The main annoyance of writing those constants in a module is the trade off between readable or short constant names. I usually resort to the static const struct trick from C to create a pseudo namespace for the codata:

2 Likes

To make these values PARAMETERs it suffices to keep the type definition private and only expose a single type instance:

type(enum_codata), parameter :: codata = enum_codata()

At least for the NAG compiler it is known that it will replace the constant values with the actual at the site of usage.

Some developers still frown upon this usage pattern, with the argument that derived types should not be (mis-)used as namespaces.

Concerning short/long names, the user could have the responsibility to choose the name he prefers:

use NIST_physical_constants, only: c=>speed_of_light_in_vacuum
1 Like

I think in case of physical constants the rename upon import is the better way to go.

The derived-type singleton approach is more suitable for constants which are closely related (e.g. a small set of terminal control characters, pieces in a chess game or something like that).

2 Likes

Sorry @ivanpribec and @awvwgk ,
can you explain me more what means:
type(enum_codata), parameter :: codata = enum_codata()
I don’t understand why we could not define simply the constants by:

  real(wp), parameter :: speed_of_light_in_vacuum = 299792458e0_wp
  real(wp), parameter :: standard_acceleration_of_gravity = 9.80665_wp 

and why public would not be the default in such a module.

Sorry if my comment was unclear. I’ve created an example to clarify:

module constants

   implicit none
   public 

   !> The following constants have high potential for name clashes
   !  if the client of the module does "use constants" without the only clause
   real(wp), parameter :: h = 6.6260715e-34_wp ! J·s = kg·m²·s⁻¹
   real(wp), parameter :: c = 299792458.0_wp ! m·s⁻¹
   real(wp), parameter :: kb = 1.380649e-23_wp ! J·K⁻¹ = kg·m²·s⁻²·K⁻¹

   !> To prevent this, the constants are encapsulated in a derived type
   !   with default initalization
   type, private :: enum_codata

     !> Planck's constant
     real(wp) :: h = 6.6260715e-34_wp ! J·s = kg·m²·s⁻¹
     !> Speed of light in vacuum
     real(wp) :: c = 299792458.0_wp ! m·s⁻¹
     !> Boltzmann's constant
     real(wp) :: kb = 1.380649e-23_wp ! J·K⁻¹ = kg·m²·s⁻²·K⁻¹
  end type

  !> We only expose a single instance of the type, which is a prameter
  !  I suppose the call to the structure constructure is not needed.
  type(enum_codata), parameter :: codata = enum_codata()

end module

program main
  !> We don't want to import h, c, or kb because these are often 
  !  used as local variable names
  !  We only need to import the derived type with constants.
  use constants, only: codata
  implicit none
  print *, "Speed of light", codata%c
end program

Some compilers are able to perform the substitution and in resulting binary the print call would contain the numerical value itself, and not a reference to a member of a derived type.

The downside of this pattern is that renaming variables upon import is not possible anymore.

4 Likes

@ivanpribec Thanks for that clear explanation!

So it’s quite like a #define constant in C. And with other compilers it’s like a C const.

Note that the risk of names collision would be very weak if we use the full names in:
https://physics.nist.gov/cuu/Constants/Table/allascii.txt
The shorter one is tau mass. And the mean length is around 25 characters…
Anyway, they are so long that rename at import will be nearly always necessary. In my opinion, it is not a problem since most of the time people work in a precise scientific field and use quite the same constants. You could create your own little module with a USE to the whole module, and put into it your favorite constants.

1 Like

I submitted a first draft to the qcel project:

The output looks like this (it can also write an include file with kind parameters):

1 Like

Should not the constants be declared double precision or real(kind=kind(1.0d0))?

The required kind parameters depend on the application that want to use this snippet, therefore they are optional.

>>> from qcelemental.physical_constants.context import write_fortran_header
>>> write_fortran_header("CODATA2018", "codata-wp.f90", kind="wp")

I just updated the above link with the output from the command here.

Whatever the way of doing it, one needs to be able to change the physical constant versions.

1 Like

I have pushed my own “small is beautiful” fpm project:

$ ./nist.py --help
usage: nist.py [-h] [-y {2010,2014,2018}] [-d] [-v]

Downloads CODATA fundamental physical constants from NIST website and generates a Fortran module.

optional arguments:
  -h, --help           show this help message and exit
  -y {2010,2014,2018}  CODATA values: 2010, 2014 or 2018
  -d                   Delete calculated values (...)
  -v                   Version

:bulb: A cool idea to validate that module (contributors welcome):
The Python script could also generate a Fortran program which would generate another allascii.txt file, at least the two first columns, with the same formatting. We could then use a diff tool to validate the Fortran module.

The main difficulty is that we need a Fortran function receiving a real value, a number of decimals and returning a character(25) string:

character(25) get_nist_format(value, decimals)

and the return strings would contain those kinds of formatted numbers:

2.4140 e-8
299 792 458
9.806 65

(including trailing blanks, the column is always 25 characters).

3 Likes

Whatever the way of doing it, one needs to be able to change the physical constant versions.

Definitely. QCElemental at present lets one choose 2014 or 2018 CODATA. The mild difficulty in adding new versions is that NIST changes up the constants names and sometimes the combinations of constants that get a name. So to allow easy update for downstream software, we made sure to map all the 2014 names into 2018. There hasn’t yet been demand for pre-2014 values, but there’s no blocker.

2 Likes

Nice script.
About the comparison with a diff tool. It is not so easy, because the nist text format is not so easy to understand.

1 Like

The general rule is to have groups of 3 digits:
4.001 506 179 127
If there is only 2 remaining digits, they can be grouped:
5.142 206 747 63 e11
If there is only 1 remaining digits, we have a 4 digits group:
9.717 362 4292 e21
But when there are “…”, it does not seem to apply:
197.326 980 4...

On the left side, same rules:
3477.23
299 792 458
10 973 731.568 160

Another problem is the exponent. Most of the time a write(*,*) seems to do a good job, but for 2.350 517 567 58 e5 it yields 235051.75675800000. In that case we would need:

a = 2.35051756758d5
print '(es17.11)', a

We should try to find the rules concerning exponents.

And it could help if we have a clue about which language they use to generate their files: txt, html, pdf. Is it a Python or Perl script? Is it SQL?
Maybe we could find a clue in the articles cited here:
https://physics.nist.gov/cuu/Reference/versioncon.shtml

I’ve found a way to get the exponent from the number of decimals and the real number.
use, intrinsic :: iso_fortran_env, only: wp=>real64
implicit none

   real (kind=wp) :: a = 2.35051756758d5
   integer        :: nb_decimals = 11

   integer            :: p10 ! contains the exponent
   integer            :: nb_digits ! contains the number of digits
   real (kind=wp)     :: ab
   character (len=10) :: fm

   integer :: number_OF_digits ! function


   p10=int(log10(abs(a)))
   write(6,*) 'log10',p10
   ab = a*10._wp**(-p10)
   write(6,*) 'ab',ab

   nb_digits = number_OF_digits(ab)
   write(6,*) 'digits',nb_digits


   ab  = ab * 10._wp ** (nb_digits-nb_decimals)
   p10 = p10 - (nb_digits-nb_decimals)

  write(fm,'("f",i0,".",i0)') nb_digits+3,nb_decimals
   !write(6,*) 'fm',fm

   write(6,'(a,' // fm // '"e",i0)') 'ab',ab,p10
   !write(6,*) 'ab',ab,'e',p10


  END
  integer FUNCTION number_OF_digits(a)
  use, intrinsic :: iso_fortran_env, only: wp=>real64
  implicit none
  real (kind=wp), intent(in) :: a


   real (kind=wp) :: b
   integer        :: i


   b = a
   i = 0
   DO
    IF (b-aint(b) < 0.01_wp) EXIT
    i = i + 1
    b = b * 10._wp
   END DO

   number_OF_digits = i
  END FUNCTION number_OF_digits

It gives:
log10 5
ab 2.3505175675800003
digits 11
ab 2.35051756758e5

1 Like

@gardhor
But with 5.142 206 747 63 e11, it yields:

$ ./a.out
 log10          11
 ab   5.1422067476299995     
 digits          15
ab 51422.06747630000e7

Oups, the number of digits is wrong. I’ll look at it!

1 Like

@loriab,

Welcome to this forum!

@loriab, @awvwgk:

A question: how does QCElemental and the C header and the copy-cat Fortran module distinguish between what are exact fundamental constants in CODATA/NIST sources and what are derived values? Does the ... in values listed at the NIST website serve as an indicator of the latter and if so, is that made use of somehow?

Consider Boltzmann constant, k as 1.380 649 x 10-23 J K-1 marked as exact.

And the Avogadro constant, NA as 6.022 140 76 x 1023 mol-1 which is also exact.

The above would yield the molar gas constant, R = k x NA = 8.314 462 618 153 24 J mol-1 K-1 per common floating-point computation.

Whereas the CODATA value appears to be set arbitrarily with 10 significant digits as 8.314 462 618 J mol-1 K-1 with three trailing dots .... This does not make sense, particularly considering the floating-point precision of many other constants.

It seems to me thorough analysis is necessary as to what constitutes fundamental constants that are exact per scientific definition or reference condition. And allow for other constants that are related to such fundamental constants to be computed using documented and consistent floating-point precision.