Morse Potential solve using Matlab

14 次查看（过去 30 天）
Naman Singh
Naman Singh 2020-9-29
评论： Kevin Lehmann 2024-1-26
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
  1 个评论
Kevin Lehmann
Kevin Lehmann 2024-1-26
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.

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

Alan Stevens
Alan Stevens 2020-9-29
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
  3 个评论
显示 1更早的评论
Alan Stevens
Alan Stevens 2020-9-30
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
Alan Stevens 2020-9-30
编辑：Alan Stevens 2020-9-30
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?

请先登录，再进行评论。

更多回答（1 个）

Naman Mathur
Naman Mathur 2020-9-30
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
  6 个评论
显示 4更早的评论
Naman Singh
Naman Singh 2020-10-3
Just curious, I'm not sure how to implement Euler methods with the 3 variables p,q,t
Alan Stevens
Alan Stevens 2020-10-3
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')

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Mathematics 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by