Obtaning values and plotting Lennard-Jones function
27 views (last 30 days)
Show older comments
Ugur Batir
on 19 Apr 2015
Commented: NURSAFIKA BAHIRA JULI
on 21 Jan 2020
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 Comments
John D'Errico
on 19 Apr 2015
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;
Accepted Answer
John D'Errico
on 19 Apr 2015
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 Comments
John D'Errico
on 19 Apr 2015
Edited: John D'Errico
on 19 Apr 2015
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
on 19 Apr 2015
Edited: John D'Errico
on 19 Apr 2015
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))
More Answers (4)
Star Strider
on 19 Apr 2015
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 Comments
Star Strider
on 27 Feb 2016
My pleasure.
I am happy it helped. I would appreciate it if you would Vote for it.
NURSAFIKA BAHIRA JULI
on 21 Jan 2020
Hye, can i run monte carlo Simulation for LJ's model?? And how?
LATEFA ALSHAMMARY
on 9 Nov 2018
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')
0 Comments
See Also
Categories
Find more on Annotations in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!