Obtaning values and plotting Lennard-Jones function

I need to evaluate the below function which is the Lennard-Jones function defining the van der Waals forces between atoms. The function and the plot it must give are attached. Sigma and epsilon are constants in the function. I tried to write the code for it in couple of different forms and also tried to do it in MS Excel but all of those gave a curve that resembles a saturation curve. Any help would be very appreciated.

[Star Strider]

Show the code you wrote. We don’t know the values of the constants or the reason your code failed.

[Ugur Batir]

 This is the code I wrote
 r=1:0.01:10;
 s=3.4;
 X=(s./r);
 F=24.*(2.*(X.^13)-(X.^7));
 plot(1./X,F)

[John D'Errico]

In fact, your plot is identical to what I produced. What you apparently don't understand is what happens for small r. The plot axes explode. So you never see the essential shape of the curve, since you went all the way down to r=1.

Try this instead:

r=3.4:0.01:10;

Answer 1

[John D'Errico]

Not too much of a problem. You just need to be careful about what you are plotting.

f = @(s) (2*s.^13 - s.^7)./s;
S = linspace(.25,1,100);
plot(1./S,f(S))
grid on

[Ugur Batir]

this might look to have the right form but isnt the right plot

[John D'Errico]

What is not right? I would suggest that you are wrong. But feel free to explain why it is NOT right. If you cannot do so, then it just means you don't understand what was done. That you failed to generate this plot also means you fail to understand how to plot it.

The scaling on the y axis is merely a scale factor. I left out epsilon, but who cares about axis scaling? You can put that in. I chose to parameterize it in terms of a variable s=sigma/r, but again, that is irrelevant. I did those things of course since you failed to provide ANY information about what they were.

[John D'Errico]

Of course, now that I know what your parameters are, I can include them, in case this makes you happy. Still trivial. Still effectively the same plot.

r = linspace(3.4,10,100);
sig = 3.4;
f = @(sig,r) 24*(2*(sig./r).^13 - (sig./r).^7)./sig;
plot(r./sig,f(sig,r))

Answer 2

[Star Strider]

I believe the information you were provided is incorrect. The equation for F actually appears to be the Lennard-Jones equation, not the van der Waals equation. Since F appears to be the integral of U w.r.t. r, and now having ‘sigma’ (I still need ‘epsilon’), this would be my approach:

U = @(e,s,r) 4*e*((s./r).^12 - (s./r).^6);  % Lennard-Jones

e = 1.0;                            % epsilon (GUESS)
s = 3.4;                            % sigma
r = linspace(0.75, 2.5)*s;

U_LJ = U(e,s,r);

F_vdW = cumtrapz(U_LJ, r);                   % van der Waals

figure(1)
plot(r/s, U_LJ/e, '--k')
hold on
plot(r/s, (F_vdW-F_vdW(end))*s/e, '-k')
hold off
grid
axis([0.75  3    -20  4.5])
legend('Lennard-Jones \itU/\epsilon\rm', 'van der Waals \itF\sigma/\epsilon')

Subtracting the last value of ‘F_vdW’ corrects for the constant-of-integration. This produces a plot that does not exactly match the sort you posted, but is reasonably close. You will have to experiment with your equations and constants to get the correct result:

[Star Strider]

I don’t have a problem approximating the curves with this code (a slight variation on my previous code), but I have the feeling that some information is missing somewhere. (It wouldn’t be the first time a paper omitted information that was important to reproducing their results.) I had to fudge the plot a bit (note the negative derivative), but I got a decent approximation to the plots you posted:

U = @(e,s,r) 24*(e/s)*(2*(s./r).^13 - (s./r).^7);   % Lennard-Jones

e = 0.0556;                         % epsilon (GUESS)
s = 3.4;                            % sigma
r = linspace(0.75, 2.5)*s;

U_LJ = U(e,s,r);                                    % L-J Evaluated

F_vdW = gradient(U_LJ, r(2)-r(1));                 % van der Waals

figure(1)
plot(r/s, U_LJ/e, '--k')
hold on
plot(r/s, -F_vdW*s/e, '-k')
hold off
grid
axis([0.75  3    -3  5])
legend('Lennard-Jones \itU/\epsilon\rm', 'van der Waals \itF\sigma/\epsilon')

The new plot:

[Rob Qualls]

Thanks for posting this. Ran into a similar issue on some data crunching for a physics class, this helped a lot.

[Star Strider]

My pleasure.

I am happy it helped. I would appreciate it if you would Vote for it.

[NURSAFIKA BAHIRA JULI]

Hye, can i run monte carlo Simulation for LJ's model?? And how?

Answer 3

[Ugur Batir]

Let me ask the question again, I'll explain the thing that I should've before and this will clear things I guess.

The function I'm trying to plot is the derivative of the Lennard-Jones potential equation wirth respect to distance, thus it simulates the van der Waals forces.

The function itself and the plot of the function is given in the attachments at my first entry. Sigma is 3.4 and the epsilon is 0,0556. Epsilon actually is not important since normalized values are to be plotted.

I actually need the values of the discrete points on the plot, but plotting it is an easy way to compare with the original plot to determine if these values are true or not.

This plot and function is obtained from a published article so it is definitely correct. When calculated with a calculator, same values shown on the plot is obtained but somehow I cannot get these values neither with matlab nor with excel.

The plot I attached shows a clear minimum at around x=1.2 y=2.5, this is why I said John's solution was not quite right.

Answer 4

[LATEFA ALSHAMMARY]

U = @(e,s,r) 24*(e/s)*(2*(s./r).^13 - (s./r).^7); % Lennard-Jones e = 0.0556; % epsilon (GUESS) s = 3.4; % sigma r = linspace(0.75, 2.5)*s; U_LJ = U(e,s,r); % L-J Evaluated F_vdW = gradient(U_LJ, r(2)-r(1)); % van der Waals figure(1) plot(r/s, U_LJ/e, '--k') hold on plot(r/s, -F_vdW*s/e, '-k') hold off grid axis([0.75 3 -3 5]) legend('Lennard-Jones \itU/\epsilon\rm', 'van der Waals \itF\sigma/\epsilon')

Answer 5

[algeed alshammari]

I need code in matlab TO PLOTE Lennard-Jones PLESE