Wrong fortran source creation floating point operations (not double)

[zerothi]

I can see your fortran interface wrappers create faulty variables.

openmm/wrappers/generateWrappers.py

Line 881 in 3aa4bb8

class FortranHeaderGenerator(WrapperGenerator):

I just built 7.4.2 and got this output:

MODULE OpenMM_Types
    implicit none

    ! Global Constants

    real*8, parameter :: OpenMM_NmPerAngstrom =  0.1
    real*8, parameter :: OpenMM_AngstromsPerNm =  10.0
    real*8, parameter :: OpenMM_PsPerFs =  0.001
    real*8, parameter :: OpenMM_FsPerPs =  1000.0
    real*8, parameter :: OpenMM_KJPerKcal =  4.184
    real*8, parameter :: OpenMM_KcalPerKJ =  1.0/4.184
    real*8, parameter :: OpenMM_RadiansPerDegree =  3.1415926535897932385/180.0
    real*8, parameter :: OpenMM_DegreesPerRadian =  180.0/3.1415926535897932385
    real*8, parameter :: OpenMM_SigmaPerVdwRadius =  1.7817974362806786095
    real*8, parameter :: OpenMM_VdwRadiusPerSigma =  .56123102415468649070

This is really bad practice ;)

check this output:

program test
  real*8 :: d
  d = 4.184
  print *, d
  d = 4.184_8
  print *, d
end program test

the output will be something like this:

   4.1840000152587891     
   4.1840000000000002     

The same applies to your arithmetic operations, they are done in floating point precision (not double).

You really need all real variables to be suffixed with _8 where you are using doubles.

[peastman]

You really need all real variables to be suffixed with _8 where you are using doubles.

Do you mean constants, not variables? I'm not a Fortran programmer! As I understand it, real*8 means the variable is double precision. You're saying that the constants are being cast to single precision before they're assigned to the variables?

I just did a bit of searching through Fortran documentation. It appears that the _8 syntax is not standardized? Fortran 90 allows you to give a format specifier in that form, but what goes after the _ differs between compilers. Instead, the universally supported way to indicate a double precision constant is with d0.

[zerothi]

Do you mean constants, not variables? I'm not a Fortran programmer! As I understand it, real*8 means the variable is double precision. You're saying that the constants are being cast to single precision before they're assigned to the variables?

Yes. Only constants. Fortran interprets all 0. and friends as floating point variables (4-byte reals). So all the math is done in float's.
Only upon assignment do they get cast to doubles if the variable is a double.

I just did a bit of searching through Fortran documentation. It appears that the _8 syntax is not standardized?

That is not true.
Kind specification is standardized with parameter specifications:
i.e.

! Only one of the following 3
integer, parameter :: dp = 8
! Preferably it should be (this will make the kind specification determined for the architecture at compile time)
integer, parameter :: dp = kind(1.d0)
! or even
integer, parameter :: dp = selected_real_kind(15)

! Now all constants may be defined through:
real(dp), parameter :: pi = 3.141_dp

This is the proper way to do it.

I would highly suggest you don't use the explicit dp = 8 specification, but rather the dp = kind(1.d0) or the last one.

It is true that 3.141d0 forces the variable to be double precision, so that could also be used, I am just more in favour of the former methods.

Note also that all real*8 should be converted to real(dp) to use the kind specification as defined above.
Lastly, if you want to use complex numbers you may have something like complex*16. These should be defined as complex(dp) (and NOT complex(16)).

but what goes after the _ differs between compilers. Instead, the universally supported way to indicate a double precision constant is with d0.

This is not true. All fortran compilers accept this standardized format. What you'll probably find are various named parameters dp/real64 (in iso module)/d or whatever. These are defined by the programmer and has nothing to do with compilers.

Note, that all the same considerations apply to your integer*8 definitions. There you should use integer, parameter :: long = selected_int_kind(18) and then define constants as 0_long. However, for integer assignments it doesn't really matter since there is no precision loss when casting from integer(4) -> integer(8). So here you are a bit more free. :)

[peastman]

Ok, thanks. If you want to submit a PR, it would be welcome!

[zerothi]

Sorry, I have no clue how to change your automatic build up :( Hope you'll get somebody to do it ;)

[peastman]

Here's the code that generates the global constants:

openmm/wrappers/generateWrappers.py

Lines 903 to 912 in 7c2e599

	def writeGlobalConstants(self):
	self.out.write(" ! Global Constants\n\n")
	node = next((x for x in findNodes(self.doc.getroot(), "compounddef", kind="namespace") if x.findtext("compoundname") == "OpenMM"))
	for section in findNodes(node, "sectiondef", kind="var"):
	for memberNode in findNodes(section, "memberdef", kind="variable", mutable="no", prot="public", static="yes"):
	vDef = convertOpenMMPrefix(getText("name", memberNode))
	iDef = getText("initializer", memberNode)
	if iDef.startswith("="):
	iDef = iDef[1:]
	self.out.write(" real*8, parameter :: OpenMM_%s = %s\n" % (vDef, iDef))

The last line in that function just needs to be modified to write out the value in whatever is the correct way.

[zerothi]

Sorry, I don't have the time to delve into and debug in a new code. Hope you get the time eventually.

[peastman]

The fix is in #3264. Could you take a look and verify that it looks correct? With this change, here is how the constants get defined.

MODULE OpenMM_Types
    implicit none

    ! Global Constants

    real*8, parameter :: OpenMM_NmPerAngstrom =  0.1_8
    real*8, parameter :: OpenMM_AngstromsPerNm =  10.0_8
    real*8, parameter :: OpenMM_PsPerFs =  0.001_8
    real*8, parameter :: OpenMM_FsPerPs =  1000.0_8
    real*8, parameter :: OpenMM_KJPerKcal =  4.184_8
    real*8, parameter :: OpenMM_KcalPerKJ =  1.0_8/4.184_8
    real*8, parameter :: OpenMM_RadiansPerDegree =  3.1415926535897932385_8/180.0_8
    real*8, parameter :: OpenMM_DegreesPerRadian =  180.0_8/3.1415926535897932385_8
    real*8, parameter :: OpenMM_SigmaPerVdwRadius =  1.7817974362806786095_8
    real*8, parameter :: OpenMM_VdwRadiusPerSigma =  .56123102415468649070_8

[zerothi]

Looks MUCH better!

I would probably do something like this to have it portable:

 def writeGlobalConstants(self): 
     self.out.write("    ! Global Constants\n\n") 
     self.out.write("    integer, parameter :: dp = kind(1.d0)\n")
     node = next((x for x in findNodes(self.doc.getroot(), "compounddef", kind="namespace") if x.findtext("compoundname") == "OpenMM")) 
     for section in findNodes(node, "sectiondef", kind="var"): 
         for memberNode in findNodes(section, "memberdef", kind="variable", mutable="no", prot="public", static="yes"): 
             vDef = convertOpenMMPrefix(getText("name", memberNode)) 
             iDef = getText("initializer", memberNode) 
             if iDef.startswith("="): 
                 iDef = iDef[1:] 
             self.out.write("    real(dp), parameter :: OpenMM_%s = %s\n" % (vDef, iDef)) 

and then with the fixes you have made, so everything is suffixed with _dp.

[peastman]

Wouldn't that have the same problem as the original, that the constants get interpreted as single precision before being assigned to the double precision variable?

[peastman]

Oh, I see. I missed the note about also adding _dp.

The difference you're suggesting is to support platforms where double precision has a width different than 8 bytes? We don't support any such platforms. I worry that might create problems where the C++ compiler sets the width of double to 8 bytes, but the Fortran compiler sets it to something else.

[zerothi]

The difference you're suggesting is to support platforms where double precision has a width different than 8 bytes? We don't support any such platforms. I worry that might create problems where the C++ compiler sets the width of double to 8 bytes, but the Fortran compiler sets it to something else.

The value 8 generally refers to the number of bytes as you would normally think. However, it is a compiler specific variable and does not have any byte meaning.

So this is about compiler compatibility. Some compilers may use kind(1.d0) == 4.

[peastman]

Does that mean I should change all occurences of real*8 to real(dp) throughout the API? What about all the uses of integer*4 and integer*8? In all these cases, I want to explicitly tell it to use a particular number of bytes, because that's what the C++ library uses, independent of how the Fortran compiler normally defines "single" or "double" precision.

[peastman]

This page clarifies the situation a little, but not entirely. If I understand correctly, prior to Fortran 2008 there is no way to do what I described? According to https://www.charmm-gui.org/charmmdoc/developer.html, CHARMM is still on Fortran 77(!). Tinker is on Fortran 95.

[zerothi]

Does that mean I should change all occurences of real*8 to real(dp) throughout the API? What about all the uses of integer*4 and integer*8? In all these cases, I want to explicitly tell it to use a particular number of bytes, because that's what the C++ library uses, independent of how the Fortran compiler normally defines "single" or "double" precision.

In principle yes. In practice it is, generally, not required. Sorry about the confusion created ;)

Generally real*8 is good enough. Even though not a standard I have never experienced a compiler not recognizing this convention. Note that this forces the reals to be 8 bytes (the page you refer to explains that this is not necessarily the correct IEEE 64-bit double).

What is important is the kind specification of the constants 3.14_kind. Here 8 is not generally meaning 8 bytes thus you would want to use integer, parameter :: dp = kind(1.d0).

My eyes just hurt at seeing real*8; why I would do real(dp). Whether you want to use real(dp) or real*8 shouldn't matter (but could be important for some obscure compiler I don't know of). But since you have to have dp for kind specification, why not use it for the definitions?

Also note that since you are using modules you are using f90, so the kind specification is already enforced by convention.

[peastman]

Got it, thanks. I just pushed the changes. Here's what the definition looks like now.

MODULE OpenMM_Types
    implicit none

    ! Global Constants

    integer, parameter :: dp = kind(1.d0)
    real*8, parameter :: OpenMM_NmPerAngstrom =  0.1_dp
    real*8, parameter :: OpenMM_AngstromsPerNm =  10.0_dp
    real*8, parameter :: OpenMM_PsPerFs =  0.001_dp
    real*8, parameter :: OpenMM_FsPerPs =  1000.0_dp
    real*8, parameter :: OpenMM_KJPerKcal =  4.184_dp
    real*8, parameter :: OpenMM_KcalPerKJ =  1.0_dp/4.184_dp
    real*8, parameter :: OpenMM_RadiansPerDegree =  3.1415926535897932385_dp/180.0_dp
    real*8, parameter :: OpenMM_DegreesPerRadian =  180.0_dp/3.1415926535897932385_dp
    real*8, parameter :: OpenMM_SigmaPerVdwRadius =  1.7817974362806786095_dp
    real*8, parameter :: OpenMM_VdwRadiusPerSigma =  .56123102415468649070_dp

But since you have to have dp for kind specification, why not use it for the definitions?

Only the principle of not fixing something that isn't broken. No one has ever reported problems with real*8, so why risk breaking something by changing it?

[zerothi]

Great, I am happy! ;)

[peastman]

Thanks! I'll merge it as soon as the tests finish.