Obtaning values and plotting Lennard-Jones function

Question

Ugur Batir 2015-4-19

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/210306-obtaning-values-and-plotting-lennard-jones-function

评论： NURSAFIKA BAHIRA JULI 2020-1-21

2.JPG
1.JPG

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.

3 个评论
显示 1更早的评论隐藏 1更早的评论

Ugur Batir 2015-4-19

在 MATLAB Online 中打开

 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 2015-4-19

在 MATLAB Online 中打开

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 2015-4-19

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/210306-obtaning-values-and-plotting-lennard-jones-function#answer_175822

在 MATLAB Online 中打开

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

3 个评论
显示 1更早的评论隐藏 1更早的评论

John D'Errico 2015-4-19

编辑：John D'Errico 2015-4-19

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 2015-4-19

编辑：John D'Errico 2015-4-19

在 MATLAB Online 中打开

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 2015-4-19

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/210306-obtaning-values-and-plotting-lennard-jones-function#answer_175829

在 MATLAB Online 中打开

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:

4 个评论
显示 2更早的评论隐藏 2更早的评论

Star Strider 2016-2-27

My pleasure.

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

NURSAFIKA BAHIRA JULI 2020-1-21

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

请先登录，再进行评论。

Answer 3

Ugur Batir 2015-4-19

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/210306-obtaning-values-and-plotting-lennard-jones-function#answer_175839

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.

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

Answer 4

LATEFA ALSHAMMARY 2018-11-9

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/210306-obtaning-values-and-plotting-lennard-jones-function#answer_346095

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')