Morse Potential solve using Matlab

In molecular dynamics simulations using classical mechanics modeling, one is often faced with a large nonlinear ODE system of the form

Mq ′′ = f(q), where f(q) = −∇U(q). (4.32) Here q are generalized positions of atoms, M is a constant, diagonal, positive mass matrix and U(q) is a scalar potential function. Also, ∇U(q) = ( ∂U ∂q1 , . . . , ∂U ∂qm ) T . A small (and somewhat nasty) instance of this is given by the Morse potential [83] where q = q(t) is scalar,

U(q) = D(1−e −S(q−q0) ) 2 , and we use the constants D = 90.5∗.4814e−3, S = 1.814, q0 = 1.41 and M = 0.9953.

How to solve this problems using non-library functions

[Kevin Lehmann]

It turns out that there is a closed form analytical expression for the vibrational motion of a Morse Oscillator, so no numerical integration is requires. If the potential is written as U(q) = D z^2 with z = 1 - exp( -b (r-re) ), and the oscillator has reduced mass mu, the angular frequency, w(E), of the oscillator as function of vibrational energy, E, is given

by w(E)^2 = 2 b^2 ( D - E) / mu. Note that it goes to zero as E -> D, the dissociation energy.

z(t, E, phi) = ( E cos(w(E) t +\phi)^2 + ( D-E) sqrt(E/D) sin( w(E) t + \phi) ) / ( D - E sin( w(E) t + \phi)^2 )

r(t,E,phi) = re - ln( 1 - z(t,E,\phi) ) / b

\phi is a phase of the oscillator.

Answer 1

[Alan Stevens]

Do you mean like the following (where I've just used 1 atom):

If so, then something like the following coding might do:

%% Solver

tspan = [0 100];
IC = [1.41 0.04];
[t, Q] = ode15s(@dqdtfn,tspan,IC);
q = Q(:,1);
plot(t,q),grid


function dqdt = dqdtfn(~,Q)

         D = 90.5*0.4814E-3;
         S = 1.814;
         q0 = 1.41;
         M = 0.9953;
         
         q = Q(1);
         v = Q(2);
         dvdt = -2*D*S*exp(-S*(q-q0))*(1-exp(-S*(q-q0)))/M;
         dqdt = [v; dvdt];
         

end

This results in the rather boring

[Alan Stevens]

What is T? Where do the brackets go in your expression for H? Is it p(T*M-1)/2 + U(q)? or something else?

[Alan Stevens]

Perhaps the following will help:

 D = 90.5*0.4814E-3;
 S = 1.814;
 q0 = 1.41;
 M = 0.9953;

tspan = [0 100];
IC = [1.41 0.04];
[t, Q] = ode15s(@dqdtfn,tspan,IC);
q = Q(:,1);
v = Q(:,2);
% Classically E would be (1/2)Mv^2 + U
% If that applies here, then:
E = 1/2*M*v.^2 + D*(1- exp(-S*(q-q0)).^2)-(1/2);
subplot(2,1,1)
plot(t,q),grid
xlabel('t'),ylabel('q')
subplot(2,1,2)
plot(t,E),grid
xlabel('t'),ylabel('E')

function dqdt = dqdtfn(~,Q)

         D = 90.5*0.4814E-3;
         S = 1.814;
         q0 = 1.41;
         M = 0.9953;
         
         q = Q(1);
         v = Q(2);
         dvdt = -2*D*S*exp(-S*(q-q0))*(1-exp(-S*(q-q0)))/M;
         dqdt = [v; dvdt];
         

end

[Alan Stevens]

For one atom that is essentially what my coding plots since (1/2)Mv^2 = p^2/(2M) (for 1 atom the transpose of p is the same as p).

What is N? The number of atoms?

Answer 2

[Naman Mathur]

Use a library nonstiff Runge-Kutta code based on a 4-5 embedded pair to integrate this problem for the Morse potential on the interval 0 ≤ t ≤ 2000, starting from q(0) = 1.4155, p(0) = 1.545 /48.888M. Using a tolerance of 1.e − 4 the code should require a little more than 1000 times steps. Plot the obtained values for H(q(t), p(t)) − H(q(0), p(0))

This is what the problem says, i think it's the number range of time 1000 and 2000

[Alan Stevens]

Well, if the non-stiff solver ode45 is to be used, then perhaps the following is what is wanted (though ode45 takes more than 1000 steps):

 D = 90.5*0.4814E-3;
 S = 1.814;
 q0 = 1.4155;
 M = 0.9953;
 v0 = 1.545 /48.888;

tspan = [0 2000];
IC = [q0 v0];
opt = odeset('AbsTol',1E-4);
[t, Q] = ode45(@dqdtfn,tspan,IC,opt);
q = Q(:,1);
v = Q(:,2);
% Classically E would be (1/2)Mv^2 + U
% If that applies here, then:
E = 1/2*M*v.^2 + D*(1- exp(-S*(q-q0)).^2)-(1/2);
subplot(2,1,1)
plot(t,q),grid
xlabel('t'),ylabel('q')
subplot(2,1,2)
plot(t,E),grid
xlabel('t'),ylabel('E')

function dqdt = dqdtfn(~,Q)

         D = 90.5*0.4814E-3;
         S = 1.814;
         q0 = 1.41;
         M = 0.9953;
         
         q = Q(1);
         v = Q(2);
         dvdt = -2*D*S*exp(-S*(q-q0))*(1-exp(-S*(q-q0)))/M;
         dqdt = [v; dvdt];
         

end

[Naman Mathur]

Thanks a lot

[Naman Singh]

So how does someone plot the Hamiltonian without using a library function,

Let's say Forward Euler or Simplectic Euler?

[Alan Stevens]

I doubt if any Euler method is better than the library methods! Not clear why non-stiff solvers were specified though.

[Naman Singh]

Just curious, I'm not sure how to implement Euler methods with the 3 variables p,q,t

[Alan Stevens]

Euler methods as follows. Easy to implement; not always very accurate (especially simple Euler):

D = 90.5*0.4814E-3;
S = 1.814;
q0 = 1.4155;
M = 0.9953;
v0 = 1.545 /48.888;

% Euler
dt = 0.1;
t = 0:dt:200;
n = numel(t);
q = zeros(n,1);
v = zeros(n,1);

q(1) = q0; v(1) = v0;
for i = 1:n-1
    
    q(i+1) = q(i) + dt*v(i);
    % Next line is simple Euler (RHS uses q(i))
    %v(i+1) = v(i) + dt*(-2*D*S*exp(-S*(q(i)-q0))*(1-exp(-S*(q(i)-q0)))/M);
    
    % Next line is symplectic Euler (RHS uses q(i+1) not q(i))
    v(i+1) = v(i) + dt*(-2*D*S*exp(-S*(q(i+1)-q0))*(1-exp(-S*(q(i+1)-q0)))/M);
end
subplot(2,1,1)
plot(t,q),grid
xlabel('t'),ylabel('q')
subplot(2,1,2)
plot(t,v),grid
xlabel('t'),ylabel('v')