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


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).


This is a series of posts centered around “things to do with old fortran code”. To be very exact, there’s only this post, but related articles are:

  1. Handling Legacy Fortran Code <– You are here
  2. MATLAB bindings for Fortran Libraries through C
  3. cibuildwheel for C++ libraries with thin Python wrappers

The third post lays out design guidelines and build practicalities for C++ libraries which may have embedded Fortran code.

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 1. 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

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!

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:

 4Wrapper for ``force_eam``.
 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)
14Other Parameters
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~s 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
 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.


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 2.


Kedward, Laurence J., Bálint Aradi, Ondřej Čertík, Milan Curcic, Sebastian Ehlert, Philipp Engel, Rohit Goswami, et al. 2022. “The State of Fortran.” Computing in Science & Engineering 24 (2): 63–72. https://doi.org/10.1109/MCSE.2022.3159862.

  1. For a copper slab and hydrogen molecule, the potential is qualitatively correct. It captures the fact that when the slab is fixed, there are effectively only two degrees of freedom, the distance between the hydrogen atoms, and the distance to the slab. ↩︎

  2. Most of the issues which crop up in embedding and distribution are largely build system related, and require no changes to the Fortran code other than a compatible iso_c_binding wrapper so they are extracted away. ↩︎