One interesting approach to obtain the actual pair interaction \(V(r)\)  starting from a known pair correlation \(g^{\mathrm {th}}(r)\) , is to reverse engineer the problem of the Monte Carlo simulation.
This method was already illustrated by Lyubartsev and Laaksonen in the context of ionic liquids physics [1]}.
The procedure consists of three main steps.
As a first step, one needs to start a preliminary Monte Carlo simulation with a guessed interaction \(V^{(0)}(r)\) .
Subsequently, when the simulated system reaches its final equilibrium configuration, one can measure its pair correlation function \(\langle g^{(0)}\rangle (r)\) .
This pair correlation function can then be compared to the known \(g^\mathrm {th}(r)\) .
Based on the deviation \(\Delta \langle g^{(0)}\rangle \)  between the measured and expected pair correlation, the interaction potential \(V^{(0)}(r)\)  is suitably corrected, to obtain \(V^{(1)}(r) = V^{(0)}(r) + \Delta V^{(0)}(r)\) .
By iterating this process one can then calculate
\(V^{(j+1)}(r) = V^{(j)}(r) + \Delta V^{(j)}(r),\) 
