Here is one possible refactor that removes the choice
flag altogether and replaces it with a procedure argument that conforms to the shared abstract interface I showed in the previous post.
Everything lives in three small modules so you can mix-and-match any potential with the same RK4 integrator.
! ===================== potentials_mod.f90 =====================
module potentials_mod
implicit none
private ! keep symbols local unless explicitly PUBLIC
!---- Public API ----
public :: potential_func ! the abstract interface name
public :: V_MT, V_ManningRosen, V_Morse, V_Hulthén ! concrete potentials
!-----------------------------------------------------------------
! Generic “shape” that every concrete potential must match
!-----------------------------------------------------------------
abstract interface
function potential_func(r, p1, p2, p3) result(Vr)
real(8), intent(in) :: r, p1, p2, p3
real(8) :: Vr
end function potential_func
end interface
contains
!-------------------------------------------------------------------
! 1) Malfliet–Tjon (MT)
!-------------------------------------------------------------------
function V_MT(r, p1, p2, p3) result(Vr)
real(8), intent(in) :: r, p1, p2, p3
real(8) :: Vr
Vr = -p2*exp(-p3*r)/r + p1*exp(-2*p3*r)/r
end function V_MT
!-------------------------------------------------------------------
! 2) Manning–Rosen
!-------------------------------------------------------------------
function V_ManningRosen(r, p1, p2, p3) result(Vr)
real(8), intent(in) :: r, p1, p2, p3
real(8) :: Vr
Vr = ( p3*(p3-1d0)*exp(-2d0*r/p2) &
/ (1d0-exp(-r/p2))**2 &
- p1*exp(-r/p2)/(1d0-exp(-r/p2)) ) &
/ (p2**2)
end function V_ManningRosen
!-------------------------------------------------------------------
! 3) Morse
!-------------------------------------------------------------------
function V_Morse(r, p1, p2, p3) result(Vr)
real(8), intent(in) :: r, p1, p2, p3
real(8) :: Vr
Vr = p1*(exp(-2d0*(r-p2)/p3) - 2d0*exp(-(r-p2)/p3))
end function V_Morse
!-------------------------------------------------------------------
! 4) Hulthén
!-------------------------------------------------------------------
function V_Hulthén(r, p1, p2, p3) result(Vr)
real(8), intent(in) :: r, p1, p2, p3
real(8) :: Vr
real(8) :: a
a = p1 - p2
Vr = -(p1**2 - p2**2)*exp(-a*r)/(1d0 - exp(-a*r)) &
+ a**2 * exp(-a*r) / (1d0 - exp(-a*r))**2
end function V_Hulthén
end module potentials_mod
! ===================== delta_mod.f90 ===========================
module delta_mod
use potentials_mod ! brings in the abstract interface
implicit none
private
public :: deltaTrial
real(8), parameter :: frac = 1.0d0 !<< you had this in the original code
contains
!-------------------------------------------------------------------
! Trial phase-shift integrand
!-------------------------------------------------------------------
function deltaTrial(r, delta0, kcm, p1, p2, p3, V) result(dTrial)
real(8), intent(in) :: r, delta0, kcm, p1, p2, p3
procedure(potential_func) :: V
real(8) :: dTrial
dTrial = -V(r, p1, p2, p3) * ( sin(kcm*r + delta0) )**2 &
/ (kcm*frac)
end function deltaTrial
end module delta_mod
! ===================== rk4_mod.f90 =============================
module rk4_mod
use delta_mod
implicit none
private
public :: PFMRK4
contains
!-------------------------------------------------------------------
! Simple RK4 driver that works with ANY potential obeying
! the potential_func interface
!-------------------------------------------------------------------
function PFMRK4(kcom, p1, p2, p3, V) result(delta_deg)
real(8), intent(in) :: kcom, p1, p2, p3
procedure(potential_func) :: V !<< key change
real(8) :: delta_deg
integer, parameter :: N = 500000
real(8), parameter :: a = 0.01d0, b = 5.0d0
real(8), parameter :: h = (b-a)/N
real(8) :: r, y, k1, k2, k3, k4
integer :: i
r = a; y = 0.0d0
do i = 1, N
k1 = deltaTrial(r, y, kcom, p1, p2, p3, V)
k2 = deltaTrial(r+0.5d0*h, y+0.5d0*k1*h, kcom, p1, p2, p3, V)
k3 = deltaTrial(r+0.5d0*h, y+0.5d0*k2*h, kcom, p1, p2, p3, V)
k4 = deltaTrial(r+h, y+k3*h, kcom, p1, p2, p3, V)
y = y + h*(k1 + 2d0*k2 + 2d0*k3 + k4)/6d0
r = r + h
end do
delta_deg = y * 180d0 / acos(-1d0) ! rad -> deg
end function PFMRK4
end module rk4_mod
! ===================== main.f90 ================================
program demo
use rk4_mod
use potentials_mod ! to get the specific potentials
implicit none
real(8) :: kcm, p1, p2, p3, delta
kcm = 1.0d0
p1 = 1.0d0
p2 = 1.0d0
p3 = 1.0d0
! ---- Choose any concrete potential simply by passing its name ----
delta = PFMRK4(kcm, p1, p2, p3, V_MT) ! Malfliet–Tjon
write(*,'("Δ(MT) = ",f8.4,"°")') delta
delta = PFMRK4(kcm, p1, p2, p3, V_Morse) ! Morse
write(*,'("Δ(Morse) = ",f8.4,"°")') delta
end program demo