How-to: Use the C API from Fortran

How-to: C API from Fortran

1 Overview

The C API provides three high-level entry points for GP-accelerated optimization, usable from Fortran via iso_c_binding:

Function	Purpose	Callback type
gpr_dimer_search	Saddle point search	gpr_force_fn (full atomistic)
gpr_neb_search	Minimum energy path	gpr_oracle_fn (simple)
gpr_minimize	Geometry optimization	gpr_oracle_fn (simple)

You can also access the GP model directly for training, prediction, and hyperparameter inspection (see Tutorial: Fortran Interop).

2 NEB on the Muller-Brown surface

The NEB search uses the simple oracle callback. Built-in oracle adapters are provided for analytical potentials:

program neb_example
  use gpr_optim
  use, intrinsic :: iso_c_binding
  implicit none

  type(AtomsConfig) :: atoms
  type(NEBConfig)   :: cfg
  type(NEBResult)   :: result
  real(c_double)    :: R_i(2), R_f(2), energy
  integer           :: i

  atoms = AtomsConfig(1)
  R_i = [-0.558d0, 1.442d0]   ! minimum A
  R_f = [ 0.623d0, 0.028d0]   ! minimum B

  cfg%mode = GPR_NEB_AIE
  cfg%use_idpp = 0             ! skip IDPP for 2D
  cfg%n_atoms = 0
  cfg%max_outer_iter = 30

  ! Use the built-in Muller-Brown oracle adapter
  result = gpr_neb_search(c_funloc(c_gpr_oracle_muller_brown), &
                          c_null_ptr, R_i, R_f, 2, 7, cfg)

  if (result%converged()) then
    print *, "Converged! Barrier:", result%barrier()
    do i = 0, result%n_images() - 1
      call result%image_energy(i, energy)
      print '(A,I2,A,F10.4)', "  Image ", i, ": E = ", energy
    end do
  end if
  ! result cleaned up automatically by FINAL
end program

3 Writing a custom oracle callback

For your own potential, write a Fortran subroutine with the oracle_callback interface:

subroutine my_potential(n_coords, positions, energy, gradient, user_data) bind(C)
  use, intrinsic :: iso_c_binding
  integer(c_int), value, intent(in) :: n_coords
  real(c_double), intent(in)  :: positions(n_coords)
  real(c_double), intent(out) :: energy
  real(c_double), intent(out) :: gradient(n_coords)
  type(c_ptr), value, intent(in) :: user_data

  ! Example: harmonic potential
  real(c_double) :: k
  integer :: i

  k = 1.0d0
  energy = 0.0d0
  do i = 1, n_coords
    energy = energy + 0.5d0 * k * positions(i)**2
    gradient(i) = k * positions(i)
  end do
end subroutine

Pass it to the minimizer:

type(AtomsConfig)    :: atoms
type(MinimizeConfig) :: min_cfg
real(c_double)       :: R_init(3), R_final(3), E_final, G_final(3)
integer              :: converged, n_calls

atoms = AtomsConfig(1)
R_init = [1.0d0, 2.0d0, -1.0d0]

min_cfg%conv_tol = 1.0d-4
min_cfg%max_iter = 100

call gpr_minimize(c_funloc(my_potential), c_null_ptr, &
                  R_init, 3, atoms, min_cfg, &
                  R_final, E_final, G_final, converged, n_calls)

if (converged == 1) then
  print *, "Converged to E =", E_final, "in", n_calls, "oracle calls"
end if

4 Dimer search with full atomistic callback

The dimer uses the full force callback, which receives atomic numbers and cell vectors:

subroutine my_force(n_atoms, R, atomic_nrs, F, energy, variance, &
                    box, user_data) bind(C)
  use, intrinsic :: iso_c_binding
  integer(c_long), value :: n_atoms
  real(c_double), intent(in)  :: R(3*n_atoms)
  integer(c_int), intent(in)  :: atomic_nrs(n_atoms)
  real(c_double), intent(out) :: F(3*n_atoms), energy, variance
  real(c_double), intent(in)  :: box(9)
  type(c_ptr), value :: user_data

  ! Call your Fortran force routine here.
  ! F should be FORCES (not gradient): F = -dE/dR
  ! The C API internally handles the sign convention.
  call my_pes_code(n_atoms, R, atomic_nrs, box, energy, F)
  variance = 0.0d0
end subroutine

5 Error handling

All procedures accept an optional ierr argument:

call model%train(R, E, G, n_obs, n_coords, ierr)
if (ierr /= GPR_SUCCESS) then
  print *, "Training failed with code", ierr
  ! Error codes: GPR_ERR_NULL_HANDLE(-1), GPR_ERR_INVALID_ARG(-2),
  !              GPR_ERR_NOT_TRAINED(-3), GPR_ERR_CHOLESKY(-5), etc.
  stop 1
end if

6 GP vs TP model

To use the Student-t process instead of GP:

type(GPModel)    :: model
type(ModelConfig) :: cfg

model = GPModel(GPR_MODEL_TP)   ! Student-t process

! Set inverse-gamma hyperprior (optional, defaults are non-informative)
cfg%tp_hyperprior_a = 2.0d0
cfg%tp_hyperprior_b = 1.0d0

call model%init(cfg, atoms)
call model%train(R, E, G, n_obs, n_coords)

The TP provides heavier-tailed predictions, useful when some force evaluations may be outliers.

7 Memory management

All derived types have FINAL subroutines. When a variable goes out of scope, the C handle is destroyed automatically. For explicit cleanup:

call model%release()
call atoms%release()