MATLAB: Fitting Lennard-Jones to Morse

curve fittingCurve Fitting Toolbox

k=1.3807E-23;
D=93.4*k;
r0=4.43E-10;
alpha=1.166E10; 
% These are my variables for Morse equation
r=1E-10:1E-12:1E-9;
V=D*(1-exp(-alpha*(r-r0))).^2;
plot(r,V)
axis([3.5E-10 9E-10 -1.0E-22 4.5E-21])
xlabel('Atomic Speration in meters')
ylabel('Potential Energy in Joules') 
%This is the plot of the Morse Equation

I'm trying to fit this plot to Lennard-Jones(12-6). D(((r0/r)^12)-2*(r0/r)^6) is the Lennard-Jones(12-6). I want to fit it to this, but I don't know how since I don't have any coefficients.

Best Answer

First you need to make these potentials comparable - now your LJ12-6 has the dissociation-energy as the zero-level, while your Morse-potential has the minimum-energy as zero-level (which makes its dissociation energy D):

LJ12_6 = @(D,r0,E0,r)D*(((r0./r).^12)-2*(r0./r).^6)+D;

Then when you plot over all your values of r, you will see that the two potentials diverge "quite a bit" for small values, whether this is relevant or not you should think about (my understanding is that you shouldn't bother much about energies more than ~3 times (or some such small factor) the dissociation energy since that would correspond to vibrational energies way larger than a molecule can have, and correspondingly small amplitudes of the wave-function). Once you've zoomed in you will see that the potentials overlap rather well for r larger than r0 but diverge for smaller values - therefore fitting makes little sense. But you can do that fitting like this:

Morse_E = @(D,alpha,r0,r)D*(1-exp(-alpha*(r-r0))).^2
LJ12_6 = @(D,r0,r)D*(((r0./r).^12)-2*(r0./r).^6)+D
pot_diff = @(p,r) sum(( LJ12_6(p(1),p(2),r) - Morse_E(D,alpha,r0,r)).^2);
[par_LJ2,out2,out3] = fminsearch(@(p) pot_diff(p,r),[D,r0]);
plot(r,[LJ12_6(D,r0,r);Morse_E(D,alpha,r0,r);LJ12_6(par_LJ2(1),par_LJ2(2),r)])

The main problem you'll have is that the two potentials model different physical systems, Morse model vibrational states in di-atmic molecules and LJ models other interactions between atoms and molecules.

HTH

Related Solutions

MATLAB: How to implement friction into the script

Friction is a force that produces acceleration which changes the velocity, so if you add an additional term to the velocity update you may get the effect that you want to. Theoretically you would need to define the mass of the missile, but this can be controlled by the parameter gamma (both are constant, costant*constant = constant), in the end you will have something like:

gamma = 1e-5;
vy = vy - (gy * dt) - gamma*abs(vy)*vy*dt;
vx = vx - (gx * dt) - gamma*abs(vx)*vx*dt;

MATLAB: ODE45 Related Question

I'm no expert, but I think if you create your own odefun, you can use an if statement to change the equation you use based on the value of t. Only one timepoint is evaluated at a time.

I used a nested function approach to avoid having to change your code too much. This gives the odefun access to the constants. I do not know, however, if the results are what you want.

sim
% created nested functions
function sim()
alpha = 0.0125;             %Ostwald N2 solubility 
D = 2*10^-8;                %Diffusion coefficient (cm^2/sec)
h = 3*10^-4;                %Bubble thickness (dyne/cm)
H = 2.5*10^8;               %Bulk modulus (dyne/cm^2)
P_B = 14.7*6.8947*10^4;     %Initial ambient pressure (dyne/cm^2)
P_F = 4.3*6.8947*10^4;      %Final ambient pressure (dyne/cm^2)
gamma = 30;                 %Surface tension (dyne/cm^2)
V_tiss = 1;                 %Tissue volume (cm^3)
M = H/V_tiss;               %Modulus of elasticity (dyne/(cm^2*cm^3))
P_N2o = 11.6*6.8947*10^4;   %N2 initial tension (dyne/cm^2)
P_met = 1.76*10^5;          %Metabolic gas tension (dyne/cm^2)
r = 0.0003;                 %Nucleation radius (cm)
t_a = 6*60;                 %Ascent time
v = (P_B - P_F)/t_a;        %Ascent rate
P_N2 = P_N2o-P_N2o*(1-exp(-.00192*120));
tspan = 0:.1:6;
f = P_B - v*tspan;
plot(tspan, f)
tspan = [0 66*60];
r0 = r;
[t,r] = ode45(@(t,r) BGIfun(t,r), tspan, r0);
   
BGI = r/r0;
plot(t/60,BGI)
grid on
xlabel('Time After Ascent (min)')
ylabel('BGI')
ylim([0 16])
    % Nested function
    function dydt= BGIfun(t,r)
        % select equation to use based on time value t
        if t <= 6*60
            dydt = ((-alpha*D/h)*(P_B-v*t+ 2*(gamma/r) + (4/3)*pi*M*r^3 ...
                -P_N2 - P_met) + r*v/3)/(P_B-v*t + (4/3)*gamma/r + (8/3)*pi*M*r^3);
        else
            dydt = ((-alpha*D/h)*(P_F + 2*(gamma/r) + (4/3)*pi*M*r^3 ...
                -P_N2 - P_met))/(P_F + (4/3)*gamma/r + (8/3)*pi*M*r^3);
        end
        
    end
end