Handling Legacy Fortran Code

A short pseudo-tutorial post on working with inherited Fortran.

Background

In-spite of the many claims by upper management (e.g. at LANL) regarding the death of Fortran, most people with even a passing interest in academic research within STEM fields will interact with Fortran code, and the community has only been improving with time Kedward et al. (2022). Rewriting is not typically an option, or a good use of time. This set of posts should be applicable to any blackbox legacy fortran code, though the specific example here is one from my own field (computational chemistry).

The Code

As an exemplary example, we will consider an implicitly typed Fortran-77 program for generating the potential energy surface (and forces) of a system containing Copper and Hydrogen atoms. The group I’m associated with has a rather nice Embedded Atom Model parameterized for this particular system ¹. The source can be obtained from here:

1git clone https://github.com/TheochemUI/fortran_cuh2_src

Since it was initially meant to work with con files, the salient points of the implementation are:

• main.f90 uses file I/O to read tmp.con and write to energy.out and forces.out, and signals that it is ready by writing \(1\) to eam.tmp
• eamroutines.f90 reads in pot.par, which has parameters used in a space separated file:

1 6.10000   6.30000
2  .000000E+00  .000000E+00  .000000E+00
3 2862.30      79.5013      86.1495
4 3.51236      2.47961      4.21154
5-109.107     -107.555      15355.5
6....

Additionally the code is in FORTRAN 77, and uses features no longer in vogue including implicit typing and common blocks.

Initial Cleanup

Rather than hard-coding a path (pot.par) it is best to store these values in the source code itself.

1potpar = np.loadtxt("pot.par", skiprows=1) # first row has 2 values
2for idx, val in enumerate(aa, start=1):
3    print(f"potpar({idx}, 1) = {val[0]}\npotpar({idx}, 2) = {val[1]}\npotpar({idx}, 3) = {val[2]}")

There are many more possible optimizations which can be done on the Fortran side, but this is the bare minimum we need to move on. The state of the project can be seen in 1461cdf. The main.f90 issues regarding the file I/O can be handled from different languages.

Python Usage

When confronted by a Fortran 77 code, or any Fortran code in general, it is never a bad idea to start out with f2py. In this instance, we would like to call the force_eam function without all the file I/O:

1f2py -c -m cuh2 eamroutines.f90

Now we can use this from Python:

 1import numpy as np
 2import cuh2 # Compiled module
 3R = [0.63940268750835, 0.90484742551374, 6.97516498544584, 3.19652040936288, 0.90417430354811, 6.97547796369474, 8.98363230369760, 9.94703496017833, 7.83556854923689, 7.64080177576300, 9.94703114803832, 7.83556986121272]
 4box = np.array([15.34, 21.702, 100])
 5natms = np.array([2, 2])
 6u = np.zeros(1)
 7F = np.zeros_like(R)
 8cuh2.force_eam(natms, box, R.ravel(), F.ravel(), u) # Call!
 9print(u)
10[-2.71140933]

This works very well, and in-fact, if this use case (flexibly calling a function) is all that is needed, then nothing more needs to be done. f2py also generates correct signatures, and doc-strings:

 1??cuh2.force_eam
 2force_eam(natms,box,r,f,u,[ndim])
 3
 4Wrapper for ``force_eam``.
 5
 6Parameters
 7----------
 8natms : input rank-1 array('i') with bounds (2)
 9box : input rank-1 array('d') with bounds (3)
10r : input rank-1 array('d') with bounds (ndim)
11f : input rank-1 array('d') with bounds (ndim)
12u : input rank-1 array('d') with bounds (1)
13
14Other Parameters
15----------------
16ndim : input int, optional
17    Default: shape(r, 0)

The resulting compiled extensions can also be published to PyPI as SciPy does.

ISO_C_BINDING wrappers for C++ and C

An iso_c_binding is essential for rational interfacing to C++ or other languages with a well-defined C-API. Since here there is only one function which needs to be called it can be as simple as:

 1module eam_wrap
 2  use iso_c_binding, only: c_double, c_int
 3    implicit none
 4contains
 5  subroutine c_force_eam (natms, ndim, box, R, F, U) bind(c)
 6    integer(c_int), intent(in) :: natms(2)
 7    integer(c_int), value, intent(in) :: ndim
 8    real(c_double), intent(in) :: box(3)
 9    real(c_double), intent(inout) :: U(1), R(ndim), F(ndim)
10    call force_eam(natms, ndim, box, R, F, U)
11  end subroutine c_force_eam
12end module eam_wrap

Which, can be mapped to the following function in any C header file:

1extern void c_force_eam(int *natms, int ndim,
2                        double *box, double *R,
3                        double *F, double *U);

We can always check that the right functions are where they should be:

 1gfortran -c eam_isoc.f90 eamroutines.f90
 2nm eam_isoc.o
 30000000000000000 T c_force_eam
 4                 U force_eam_
 5nm eamroutines.o
 6000000000000000c C counters_
 70000000000000050 C cutoff_values_
 800000000000065c8 C eam_
 90000000000000130 C eamdblexp_
1000000000000025e8 T eamfcu_
11000000000000249b T eamfh_
120000000000000000 T eamh2cu_
1300000000000029ef T eamphicucu_
1400000000000028f8 T eamphihcu_
150000000000002801 T eamphihh_
160000000000002331 T eamrhocu_
1700000000000021c7 T eamrhoh_
18                 U exp
190000000000003210 T force_eam_
20                 U _gfortran_stop_string
21                 U _gfortran_st_write
22                 U _gfortran_st_write_done
23                 U _gfortran_transfer_character_write
24                 U _gfortran_transfer_integer_write
25                 U _gfortran_transfer_real_write
260000000000002ae6 T potinit_
27                 U __powidf2
28                 U __stack_chk_fail

MATLAB Interfaces

Nominally, MATLAB has a Fortran API. However, it isn’t great. The C++ API is “modern” to the extent of of C++11, which is practically useless. So the C MEX API is the best way forward here, though input arrays need to be transposed if they are passed to Fortran routines.

Conclusions

This was a first pass at what to do when confronted with Fortran code, short of running away, or rewriting entirely. Follow up posts handle integrating into larger C++ codes, or distributing said libraries with pip ².

References