Help me design a potential function for a physics simulation

Question

In a certain type of physics simulation (Lennard-Jones potentials), it is common to use functions of the form:

$f(x) = (\frac{1}{x})^{12} - (\frac{1}{x})^{6} - (\frac{q}{x})$

(where x is between 0 and infinity, q is a constant)

However these functions grow very rapidly with decreasing x < 1, and are undefined at x = 0. It would be very useful to have a modified function which is (almost) the same as the original when x >= 1, but grows less for x < 1, and has certain other nice properties:

• continuously differentiable (ideally twice)
• smooth (absolute value of the derivatives is as small as possible everywhere)
• difference from the original function is as small as possible for x >= 1
• $\lim_{x\to 0} f(x) $ is finite (nice if the value can be chosen at will)
• monotonic for 0 <= x < 1
• has a simple functional form which can be calculated efficiently on computer hardware (using add, subtract, multiply, divide, exp, sin and cos, and a relatively small number of floating point operations)

An example of the original function with q=0.1 (blue) and a modified function for x < 1 (orange; just sketched for illustration of what the modified function might look like, not actually based on a useful analytical expression):

I'm having a bit of trouble designing a function like that by hand. It seems like this can be done essentially by splicing together three functions: a constant for small x, some connecting function for somewhat larger x < 1, and the original function for x >= 1 (also making sure the derivatives are continuous at the splice points). I don't know what the connecting function should look like though, other than it should be a smooth transition (something like a spline).

If this is formulated as a multi-objective minimization problem (eg minimizing the max difference from the goal functions plus the max abs value of the derivative), I guess that could be done using some kind of basis set, but I don't know what basis set (and using a large basis set goes against the desire for having a simple functional expression). Maybe a spline where the spline coefficients are picked based on the minimization problem?

Is it possible to design a simple function that satisfies the above conditions? How would one go about finding it?

Answer 1

Local good results can be obtained by adjusting a smooth function like

$$ g(x,r,\mu,s) = \sum_{k=1}^nr_ke^{-\left(\frac{x-\mu_k}{s_s}\right)^2} $$

Follows a MATHEMATICA script showing that

f[x_, q_] := 1/x^12 - 1/x^6 - q/x
q = 1;
xmax = 3;
n = 3;
g[x_] := Sum[r[k] Exp[-(x - mu[k])^2/s[k]^2], {k, 1, n}]
data = Table[{x, f[x, q]}, {x, 0.91, xmax, 0.03}];
gr0 = ListPlot[data, PlotStyle -> {Thick, Red}];
obj = Sum[(data[[k, 2]] - g[data[[k, 1]]])^2, {k, 1, Length[data]}]/ Length[data];
vars = Flatten[Table[{r[k], mu[k], s[k]}, {k, 1, n}]];
sol = NMinimize[obj, vars, Method -> "DifferentialEvolution"]
gx = g[x] /. sol[[2]]
gr1 = Plot[gx, {x, 0.4, xmax}];
gr2 = Plot[f[x, q], {x, 0.91, xmax}];
Show[gr0, gr1, gr2, PlotRange -> {-3, 1}]

Answer 2

If you want to keep the entire $x>1$ regime, you can smoothly interpolate inside the $x<1$ part with a bump function to a constant potential: $$ V(x)=(1-\lambda(x))f_{LJ}(x)+\lambda(x)V_0,\quad \lambda(x)=\begin{cases} \exp\left(\frac{x^2}{x^2-1}\right) & x<1\\ 0 & x\geq 1 \end{cases} $$ for some $V_0> f_{LJ}(1)=-q$.

If you merely want to keep the asymptotic at large $x$ but not necessarily the values, you can introduce a range-separation $V(x)=f_{LJ}(x)\phi_{LR}(x)$. We want $\phi_{LR}(x)\approx cx^{12}$ for small $x$ to kill the pole at $x=0$, and $\phi_{LR}(x)\sim 1$ for large $x$. However, there is no easy way to state the monotonicity requirement for $\phi_{LR}$ except to susbstitute and differentiate. Play around with something like $\phi_{LR}(x)=\frac{x^{12}}{x^{12}+\exp(g(x))}$ starting with $g(x)=c_0$ and gradually increases the degree of $g$ if it doesn't work. Alternatively, you could try using separate range-separation for the Coulomb term $q/x$ and the neutral 12-6 potential. which is easier to guarantee separate monotonicity.