RStan has the cleanest interface for doing something like this: https://cran.r-project.org/web/packages/rstan/vignettes/external.html
The code has to be written in C++ to link to everything else though. I don’t know about Pystan. The other question is the autodiff. Do you want to let autodiff compute your gradients or do you want do specify them manually?
There’s various threads around here with people doing this: Problem implementing Zipf (power law) distribution -- external C++ or What is the best way to add a user function that is used in likelihood calculation?
Welcome to the forums, btw!