using bobyqa for tpd

app/tpd.f90

@@ -6,9 +6,10 @@ program main
     use yaeos_phase_equilibria_stability, only: tangent_plane_distance
     type(Substances) :: compos
     class(ArModel), allocatable :: model
+    integer, parameter :: n = 2
 
-    real(pr) :: tc(2), pc(2), af(2), lnfug_z(2), lnfug_w(2)
-    real(pr) :: w(2), z(2), t, p
+    real(pr) :: tc(n), pc(n), af(n), lnfug_z(n), lnfug_w(n)
+    real(pr) :: w(n), z(n), t, p
     real(pr) :: tpd
 
     integer :: i
@@ -23,15 +24,31 @@ program main
     t = 294
 
     call fugacity_tp(model, z, t, p, root_type="stable", lnfug=lnfug_z)
-
     lnfug_z = lnfug_z - log(P)
-    do i=1,99
-        w(1) = real(i, pr)/100
-        w(2) = 1 - w(1)
-        call fugacity_tp(model, w, t, p, lnfug=lnfug_w, root_type="stable")
 
-        call tangent_plane_distance(z, lnfug_z, w, lnfug_w-log(P), tpd)
-        print *, i, tpd
-    end do
+    opt: block
+        use powellopt, only: bobyqa
+        real(pr) :: wl(2), wu(2)
+        integer :: npt
+
+        npt = (n+2 + (n+1)*(n+2)/2) / 2
+        wl = 1.e-5_pr
+        wu = 1
+        call bobyqa(n, npt, w, wl, wu, 0.01_pr, 1.e-5_pr, 0, 100, func)
+        print *, w/sum(w)
+    end block opt
+
+contains
 
+    subroutine func(n, x, f)
+        integer, intent(in) :: n
+        real(pr), intent(in) :: x(:)
+        real(pr), intent(out) :: f
+
+        real(pr), dimension(size(x)) :: w, lnfug_w
+
+        w = x/sum(x)
+        call fugacity_tp(model, w, t, p, lnfug=lnfug_w, root_type="stable")
+        call tangent_plane_distance(z, lnfug_z, w, lnfug_w-log(P), f)
+    end subroutine
 end program

fpm.toml

@@ -23,7 +23,8 @@ source-form = "free"
 stdlib = "*"
 test-drive = {git = "https://github.com/fortran-lang/test-drive"}
 json-fortran = {git="https://github.com/jacobwilliams/json-fortran"}
-
+nlopt-f = {git="https://github.com/grimme-lab/nlopt-f"}
+PowellOpt = { git="https://github.com/jacobwilliams/PowellOpt.git" }
 
 [[example]]
 name       = "example"