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.
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
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
Welcome to this forum!
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.