Tutorial: Fortran Interop¶
Tutorial: Fortran Interop¶
1 What you will learn¶
By the end of this tutorial you will have:
Called analytical potentials from Fortran
Trained a GP model on energy/gradient data
Predicted at new points and inspected variance
Queried and modified kernel hyperparameters
All using ergonomic Fortran derived types (no raw type(c_ptr) needed).
2 Prerequisites¶
gpr_optimbuilt with Meson (or CMake with Fortran enabled)gfortran(GCC 8+ for Fortran 2008 features)
3 Step 1: Use the module¶
The Fortran module gpr_optim provides derived types that wrap the C API.
Everything is accessed through use gpr_optim:
program gp_example
use gpr_optim
use, intrinsic :: iso_c_binding
implicit none
type(GPModel) :: model
type(AtomsConfig) :: atoms
type(ModelConfig) :: cfg
integer :: ierr
4 Step 2: Evaluate an analytical potential¶
The module includes built-in potentials for testing. No external calculator needed:
real(c_double) :: xy(2), E, G(2)
! Muller-Brown at its deepest minimum
xy = [-0.558d0, 1.442d0]
call gpr_muller_brown_energy(xy, E, G)
print *, "Energy:", E ! approximately -146.7
print *, "Gradient:", G ! approximately [0, 0] at minimum
Also available: gpr_lj_energy and gpr_leps_energy.
5 Step 3: Generate training data¶
Sample four points on the Muller-Brown surface. The GP learns from both energies and gradients:
integer, parameter :: n_obs = 4, n_coords = 2
real(c_double) :: R(8), E_train(4), G_train(8)
integer :: i
! Four sample points: three minima + one saddle
R(1:2) = [-0.558d0, 1.442d0]
R(3:4) = [ 0.623d0, 0.028d0]
R(5:6) = [-0.050d0, 0.467d0]
R(7:8) = [ 0.213d0, 0.293d0]
do i = 1, n_obs
call gpr_muller_brown_energy(R(2*(i-1)+1:2*i), E_train(i), G_train(2*(i-1)+1:2*i))
end do
6 Step 4: Create and train the GP¶
Three objects are needed: an atoms configuration, a model configuration, and the GP model itself:
! 1 "atom" for a 2D test problem
atoms = AtomsConfig(1)
! Create a GP model (use GPR_MODEL_TP for Student-t process)
model = GPModel(GPR_MODEL_GP)
! Configure: reduce SCG iterations for this small problem
cfg%max_iter = 200
cfg%report_level = 0
! Initialize and train
call model%init(cfg, atoms, ierr)
call model%train(R, E_train, G_train, n_obs, n_coords, ierr)
print *, "Trained on", model%n_observations(), "observations"
The ierr argument is optional. If omitted, errors are silent. If present,
check against GPR_SUCCESS (which is 0).
7 Step 5: Predict and check variance¶
Predict at a point between the training data:
real(c_double) :: R_q(2), E_pred, G_pred(2), var_E, var_G(2)
R_q = [0.0d0, 0.7d0]
call model%predict(R_q, E_pred, G_pred, n_coords)
call model%predict_variance(R_q, var_E, var_G, n_coords)
print *, "Predicted energy:", E_pred
print *, "Energy variance:", var_E
At a training point, the variance should be near zero:
R_q = R(1:2) ! first training point
call model%predict_variance(R_q, var_E, var_G, n_coords)
print *, "Variance at training point:", var_E ! ~0
8 Step 6: Inspect hyperparameters¶
After training, the SCG optimizer has found the kernel hyperparameters:
integer :: nls
real(c_double), allocatable :: ls(:)
real(c_double) :: ms2
nls = model%n_lengthscales()
allocate(ls(nls))
call model%get_lengthscales(ls)
print *, "Lengthscales:", ls
call model%get_magn_sigma2(ms2)
print *, "Signal variance:", ms2
You can also modify hyperparameters and retrain:
call model%set_magn_sigma2(1.0d0)
9 Step 7: Add data incrementally¶
If you get a new observation, use update instead of retraining from
scratch:
real(c_double) :: R_new(2), E_new(1), G_new(2)
R_new = [0.5d0, 0.5d0]
call gpr_muller_brown_energy(R_new, E_new(1), G_new)
call model%update(R_new, E_new, G_new, 1, n_coords)
print *, "Now have", model%n_observations(), "observations"
10 Building¶
Compile the Fortran module first, then link against libgpr_optim:
# With Meson build
gfortran -c fortran/gpr_optim.f90 -o gpr_optim.o
gfortran examples/fortran_gpr.f90 gpr_optim.o \
-L builddir -lgpr_optim -lstdc++ -o fortran_gpr
./fortran_gpr
11 Next steps¶
How-to: Use the C API from Fortran for dimer, NEB, and minimizer
How-to: Use the T-process for outlier-robust regression
How-to: Run GP-NEB for path optimization