Matthew Dennison

Pointing to functions in Fortran

One useful feature of Fortran that I recently discovered are procedure pointers, which allow us to effectively alias a function (or subroutine) name. I wanted to write a program which could calculate virial coefficients for a range of particle shapes (e.g. spheres, spheroids, spherocylinders etc.). The basic code is the same in all cases, except for the subroutine which checks for overlaps between the particles. Usually, I would just switch the subroutines as needed and recompile, but wouldn't it be great to just specify the shape in an input file? Since the overlap routine is called a lot, using an IF statement each time we call it is a non-starter. Instead, we can use a procedure pointer at the start of the program to select the function we want.

Setting up the interface

The first thing we have to do is set up an interface block for our procedure. This tells the main program how to use the procedure and defines the what arguments the it will take. Both this, and the multiple overlap-test subroutines, are in the same module called "overlaps".

MODULE overlaps

  INTERFACE

     SUBROUTINE sub_interface(r, u1, u2, over)
       REAL(8) :: r(3), u1(3), u2(3)
       LOGICAL :: over
     END SUBROUTINE sub_interface

  END INTERFACE

  PROCEDURE(sub_interface),pointer :: overlap

Each overlap test will take three real arrays as inputs. r is the separation between two particles, and u1 and u2 are the unit vectors giving the particle orientations. The parameter over, which has type logical, is returned, and is set to be TRUE when the particles overlap and FALSE when they do not. At the end, we declare a pointer with name "overlap", which will point to one of the overlap tests we will define next.

Overlap tests

Now we have to define the procedures. The exact code is excluded here for brevity, but note how the subroutines take the same variables as the interface above. Also note that it does not matter what we call these in the functions, just that their types (and shapes, if they are arrays) matches those in the interface.

CONTAINS


  SUBROUTINE overlap_spherocylinder(r12, u1, u2, over)
    REAL(8) :: r12(3), u1(3), u2(3)
    LOGICAL :: over
    ! overlap test for spherocylinders
  END SUBROUTINE overlap_spherocylinder


  SUBROUTINE overlap_cut_sphere(rij, ui, uj, over)
    REAL(8) :: rij(3), ui(3), uj(3)
    LOGICAL :: over
    ! overlap test for cut-spheres
  END SUBROUTINE overlap_cut_sphere


  SUBROUTINE overlap_sphere(r12, u1, u2, over)
    REAL(8) :: r12(3), u1(3), u2(3)
    LOGICAL :: over
    ! overlap test for spheres
  END SUBROUTINE overlap_sphere


  SUBROUTINE overlap_oblate_spherocylinder(r12, u1, u2, over)
    REAL(8) :: r12(3), u1(3), u2(3)
    LOGICAL :: over
    ! overlap test for oblate spherocylinders
  END SUBROUTINE overlap_oblate_spherocylinder


  SUBROUTINE overlap_spheroid(rab, ua, ub, over)
    REAL(8) :: rab(3), ua(3), ub(3)
    LOGICAL :: over
    ! overlap test for spheroids
  END SUBROUTINE overlap_spheroid

END MODULE overlaps

Choosing the procedure

Finally, we have to select which of the procedures we want to use, and point to it. In the main programme we define a character variable shape, which will have the name of the shape we want to work with. An IF statement is then used to choose which overlap test we will point at.

PROGRAM calc_virial

  USE overlaps

  IMPLICIT NONE

  CHARACTER(40) :: SHAPE

  SHAPE = TRIM(shape)

  IF (shape .EQ. 'sphere') THEN
     overlap => overlap_sphere
  ELSEIF (shape .EQ. 'spherocyl') THEN
     overlap => overlap_spherocylinder
  ELSEIF (shape .EQ. 'oblate_spherocyl') THEN
     overlap => overlap_oblate_spherocylinder
  ELSEIF (shape .EQ. 'cut_sphere') THEN
     overlap => overlap_cut_sphere
  ELSEIF (shape .EQ. 'spheroid') THEN
     overlap => overlap_spheroid
  ENDIF

END PROGRAM calc_virial