How to fit data to a Lennard-Jones potential in Gnuplot

Question

I'm working on a project where I need to fit a dataset to a Lennard-Jones potential in Gnuplot, but I'm having trouble achieving it correctly. The dataset I have consists of two columns: the first column represents distances (r), and the second column represents potential values.

Here are the data:

1   -438.10479452
1.1 -445.94323088
1.2 -452.07508054
1.3 -456.20604417
1.4 -458.81075037
1.5 -460.28321318
1.6 -460.99208193
1.7 -461.18673038   
1.8 -461.07112685
1.9 -460.79716730
2   -460.53742177
2.1 -460.24917871
2.2 -459.99323880
2.3 -459.76274023
2.4 -459.56247520
2.5 -459.37823253
2.6 -459.21828504
2.7 -459.09018076
2.8 -458.98484535
2.9 -458.89741339
3   -458.94036802
4   -458.82477397
5   -458.77588758
6   -458.75509627
7   -458.74829025
8   -458.58961604
9   -458.73566475
10  -458.73375523

My goal is to fit these data to a Lennard-Jones potential using Gnuplot. However, I'm having trouble finding the best way to do it.

How can I perform this fit in Gnuplot?

f(x)= (4*epsilon) * ( (sigma**12/x**12) - (sigma**6/x**6) )
fit f(x) 'data.dat' using 1:2 via sigma,epsilon
plot f(x)  lw 3 lc "blue" 

Answer 1

If you check your function, it passes 0 when sigma = x and approaches 0 for large x. You will not achieve this with your data, because you have an y-offset of about -458 and probably also some x-offset. So, you either need to shift your data or introduce some x,y offsets in your function.

f(x) = 4*epsilon * ( (sigma/(x-x0))**12 - (sigma/(x-x0))**6) + y0

So, the best I can get is the following:

Data: SO78162265.dat

1   -438.10479452
1.1 -445.94323088
1.2 -452.07508054
1.3 -456.20604417
1.4 -458.81075037
1.5 -460.28321318
1.6 -460.99208193
1.7 -461.18673038
1.8 -461.07112685
1.9 -460.79716730
2   -460.53742177
2.1 -460.24917871
2.2 -459.99323880
2.3 -459.76274023
2.4 -459.56247520
2.5 -459.37823253
2.6 -459.21828504
2.7 -459.09018076
2.8 -458.98484535
2.9 -458.89741339
3   -458.94036802
4   -458.82477397
5   -458.77588758
6   -458.75509627
7   -458.74829025
8   -458.58961604
9   -458.73566475
10  -458.73375523

Script:

### fitting
reset session

FILE = "SO78162265.dat"

f(x) = 4*epsilon * ( (sigma/(x-x0))**12 - (sigma/(x-x0))**6) + y0
epsilon = 1     # some starting values
sigma   = 1.4
x0      = -0.4
y0      = 459
set fit brief nolog
fit f(x) FILE u 1:2 via sigma,epsilon, x0, y0
set samples 1000

plot [1:10] FILE u 1:2 w p pt 7 lc "red", \
       f(x) lw 3 lc "blue"
### end of script

Result:

Final set of parameters            Asymptotic Standard Error
=======================            ==========================
sigma           = 3.44155          +/- 0.133        (3.866%)
epsilon         = 2.32564          +/- 0.2699       (11.61%)
x0              = -2.04519         +/- 0.1347       (6.584%)
y0              = -458.378         +/- 0.1771       (0.03864%)