\documentclass[12pt]{article}
\textwidth 16cm
\textheight 23cm
\topmargin 0cm
\evensidemargin 0cm
\oddsidemargin 0cm
\parindent=0pt
\parskip=2pt
\begin{document}


\begin{center}
\Large \bf
Computational Methods and Data Analysis 2003\\
\large
\end{center}
\vspace{0.4cm}
\noindent
{\bf Ex.~3: Dynamics of the O$_2$ molecule}\\

We are there, in this exercise you are going to use Matlab to study a true
physical system. By means of several algorithms for the integration of ODE's,
you will follow the vibrational dynamics of two bonded oxigen atoms (O$_2$)
interacting through a realistic potential. The results should be summarized in
a report with text and figures. Do not forget to label the axis of each figure, possibly indicating the units. Describe the results and give in the text
and/or in the captions
all the values of the parameters which have been used in each calculation.
The idea is that by reading the report and by
using your programs, the reader should be able to reproduce the results.\\
\\
{\bf Interatomic potential}\\
The interatomic energies $V(r)$ as a function of the
interatomic distances $r$ of the $O_2$  molecule, calculated by the group of
Theoretical Chemistry,
are given in file {\tt O2.txt}, which is located
on the web site\\
{\tt http:/www.theochem.kun.nl/\verb+~+pwormer} clicking on the
subject ``Courses and teaching material'' and then ``Matlab
lectures''. Energies are given in Hartrees (1 Hartree = 2 Rydbergs
= 2 $\times$ 13.6058 eV) for a discrete set of values of
interatomic distances expressed in Bohr ~(1 Bohr = 5.29
10$^{-11}$m). You can load this text file in your matlab program
and fit the data using the following interatomic potential (Morse
potential):
\begin{equation}
\label{e.morse}
V_M(r)=D_0+D_1\exp(-\beta r)+D_2\exp(-2\beta r),
\end{equation}
where $r$ is the interatomic distance and $D_0$, $D_1$, $D_2$ and $\beta$
are fitting parameters which you have to determine.
In order to find them you start by determining $D_0$, $D_1$ and $D_2$ for
different values of $\beta$ (ranging for example from $0.1$ to $3$).
Equating Eq.~(\ref{e.morse}) to the exact values of the interactomic potential
$V$ leads to
a system of $n$ equations (where $n$ is the number of
points given in file {\tt O2.txt}) with $D_0$, $D_1$ and $D_2$ as unknowns.
This system can be recast in the form
\begin{equation}
\label{e.system}
A{\bf D}={\bf V}
\end{equation}
where ${\bf D}$ is a vector of dimension $3$
($D(1)=D_0$, $D(2)=D_1$ and $D(3)=D_2$), ${\bf V}$ is a vector of dimension
$n$ and $A$ is the coefficient matrix (of dimension $n\times 3$). The
solution of (\ref{e.system}) is given by
\begin{equation}
{\bf D}=A^{-1}{\bf V}
\end{equation}
The inverse of a matrix is performed by the backslash operator in Matlab:\\
{\tt  D = A $\backslash$ V}. \\
In this way, for every fixed $\beta$, the
values of $D_0$, $D_1$ and $D_2$ can be obtained. The optimal value of
$\beta$ is the one which minimizes the discrepancy between the fitted values
{\bf fit}$=A{\bf D}$ and the actual data ${\bf V}$, calculated in the following way
\begin{equation}
{\rm norm}({\bf V}-{\bf fit})=\sqrt{\Sigma_{i=1}^n(V_i - fit_i)^2}
\end{equation}
 A possible Matlab code to perform the fit is given below:
\begin{verbatim}
      load O2.txt;
      r = O2(:,1);
      V = O2(:,2);
      clear O2;
      n = length(r);
      for i=1:300;
          beta=0.1+0.01*i;
        A = [ones(n,1), exp(-beta*r), exp(-2*beta*r)];
        D = A\V;
        fit = A*D;
        T(:,i)=[norm(V-fit)];
      end
      [Y,I]=min(T)
      Beta=0.1+0.01*I
      A = [ones(n,1), exp(-Beta*r), exp(-2*Beta*r)];
      D = A\V
      fit = A*D;
\end{verbatim}
Plot the data and the calculated fit
 on the same graph, showing that
indeed the potential $V_M$ given in Eq.~(\ref{e.morse}) is a good
approximation.
>From the plot you can also find the position of the minimum of this potential
either graphically or by use of the matlab function {\tt fminsearch}.
This minimum can also be calculated analytically. Differentiating
Eq.~(\ref{e.morse}) with respect to $r$ one gets
\begin{equation}
\label{e.derivative1}
\frac{dV}{dr}=-\beta D_1\exp(-\beta r)-2\beta D_2\exp(-2\beta r).
\end{equation}
The minimum is obtained when $\frac{dV}{dr}=0$, which can be solved making the
substitution $x=\exp(-\beta r)$:
$$
2D_2x^2+D_1x=0,
$$
which has (besides the solution for $x=0$, i.e. $r=\infty$) the solution
$$
x=x_{min}=-\frac{D_1}{2D_2}.
$$
Thus the equilibrium interatomic distance is given by
\begin{equation}
\label{e.minimum}
r_{min}=-\frac{\ln x_{min}}{\beta}=
\frac{\ln\left(-\frac{2D_2}{D_1}\right)}{\beta}.
\end{equation}
Check that the minimum calculated by using Eq.~(\ref{e.minimum}) is in
agreement with that determined graphically by plotting the data and by use of
{\tt fminsearch}.\\
\\
{\bf Harmonic approximation}\\
It is known from classical mechanics that every system characterized by a
potential which has a stable equilibrium state can be described by
means of a harmonic potential for small oscillations around the equilibrium.
In fact if V(r) is a generic potential which has a local minimum for
$r=r_{min}$ we can always perform an expansion for $r$ close to $r_{min}$ up
to the second order:
\begin{equation}
\label{e.expansion}
V(r)\simeq V(r_{min})+\frac{1}{2}\frac{d^2V}{dr^2}\Big|_{r=r_{min}}(r-r_{min})^2,
\end{equation}
where we have used the fact that the linear term is zero since
$\frac{dV}{dr}\Big|_{r=r_{min}}=0$.
Thus, using Eq.~(\ref{e.expansion}), we can approximate the full potential
given in Eq.~(\ref{e.morse}) by a harmonic one for small deviations from
equilibrium. All we need is to calculate the second derivative.
>From Eq.~(\ref{e.derivative1})
\begin{equation}
\frac{d^2V}{dr^2}=\beta^2D_1\exp(-\beta r)+4\beta^2D_2\exp(-2\beta r).
\end{equation}
For $r=r_{min}$ we obtain, using the expression for $r_{min}$
(Eq.~(\ref{e.minimum}))
\begin{equation}
K\equiv\frac{d^2V}{dr^2}\Big|_{r=r_{min}}=
\beta^2D_1\left(-\frac{D_1}{2D_2}\right)+
4\beta^2D_2\left(-\frac{D_1}{2D_2}\right)^2=\frac{\beta^2D_1^2}{2D_2}
\end{equation}
where we have defined the force constant $K$ as the second derivative of the
potential calculated at the minimum. In this way we can write the harmonic
potential which approximates $V(r)$:
\begin{equation}
V^{harm}(r)=V(r_{min})+\frac{1}{2}K(r-r_{min})^2.
\end{equation}
Plot also $V^{harm}(r)$ and check that it approximates $V(r)$ in the vicinity
of $r_{min}$.\\
\\
{\bf Dynamics}\\
Once you have fixed the form of the interatomic potential you can study
the dynamics of the problem solving the equations of motion. However a
complete anaytical solution can be only achieved for simple potentials, such
as the harmonic potential. In general you have to solve the problem by
integrating the equations of motion numerically.
Suppose that our $O_2$ molecule can only move in one dimension.
Denote the coordinates of the oxygen atoms respectively by $x_1$ and $x_2$.
The interatomic potential is thus
$$
V(x_1,x_2)=D_0+D_1\exp(-\beta(x_2-x_1))+D_2\exp(-2\beta(x_2-x_1))
$$
($r=x_2-x_1$) and the equations of motion are
\begin{equation}
\left\{ \begin{array}{ccc}
m\ddot{x}_1 & = & -\frac{\partial V}{\partial x_1}\\
m\ddot{x}_2 & = & -\frac{\partial V}{\partial x_2}
\end{array} \right.
\end{equation}
where $m$ is the oxygen mass (16 a.u.). The right-hand side of these
equations can be easily calculated:
\begin{equation}
-\frac{\partial V}{\partial x_1}=-\frac{dV}{dr}\frac{\partial r}{\partial x_1}=-\beta(D_1\exp(-\beta(x_2-x_1))+2D_2\exp(-2\beta(x_2-x_1)))
\end{equation}
\begin{equation}
-\frac{\partial V}{\partial x_2}=-\frac{dV}{dr}\frac{\partial r}{\partial x_2}=\beta(D_1\exp(-\beta(x_2-x_1))+2D_2\exp(-2\beta(x_2-x_1)))
\end{equation}
As expected, $\frac{\partial V}{\partial x_1}=
-\frac{\partial V}{\partial x_2}$. In this way the equations of motion are
\begin{equation}
\label{e.eqmotion}
\left\{ \begin{array}{ccc}
m\ddot{x}_1 & = & -\beta(D_1\exp(-\beta(x_2-x_1))+2D_2\exp(-2\beta(x_2-x_1))) \\
m\ddot{x}_2 & = & \beta(D_1\exp(-\beta(x_2-x_1))+2D_2\exp(-2\beta(x_2-x_1)))
\end{array} \right.
\end{equation}

{\bf Initial conditions}

Consider the following initial conditions:
$$
x_1=-\frac{1}{2}(r_{min}+\epsilon), \qquad x_2=\frac{1}{2}(r_{min}+\epsilon)
$$
$$
{\dot x}_1=0, \qquad {\dot x}_2=0
$$
for a few values of $\epsilon$
(for instance $|\epsilon|$ up to 0.3).

{\bf Numerical solution}

Several methods can be used to solve a system of ordinary differential
equations. For example, Matlab has a built-in routine, called ODE45, which
uses the fourth-order Runge-Kutta method with adaptive step size. The routine
requires as input a parameter (tolerance) which gives the accuracy of
the result. Set it to $10^{-5}$. We will consider this the
exact solution of the problem and will use it for comparison while using the
other methods. To use ODE45 the original system of ordinary differential
equations has
to be written as a system of first order differential equations. This can be
easily done by making the substitutions
\begin{equation}
\left\{ \begin{array}{ccc}
y_1 & = & x_1 \\
y_2 & = & {\dot x}_1 \\
y_3 & = & x_2 \\
y_4 & = & {\dot x}_2
\end{array} \right.
\end{equation}
Thus, Eq.~(\ref{e.eqmotion}) becomes
\begin{equation}
\label{e.firstorder}
\left\{ \begin{array}{lll}
\dot{y}_1 & = & y_2 \\
\dot{y}_2 & = & -\beta(D_1\exp(-\beta(y_3-y_1))+2D_2\exp(-2\beta(y_3-y_1)))/m \\
\dot{y}_3 & = & y_4 \\
\dot{y}_4 & = & \beta(D_1\exp(-\beta(y_3-y_1))+2D_2\exp(-2\beta(y_3-y_1)))/m
\end{array} \right.
\end{equation}
Note that the integration time should be much larger than one oscillation
period. Besides this built-in Matlab routine you can write your own code
to solve Eq.~(\ref{e.eqmotion}). We suggest four alternative methods.

{\bf Euler}

This method can be applied to a system of first order ordinary differential
equations of the type
\begin{equation}
\frac{dy_{i}(t)}{dt}=f_{i}(y_1(t),y_2(t),...,y_n(t);t), \qquad i=1,...,n
\end{equation}
($n=4$ in Eq.~(\ref{e.firstorder})).
Then, the solution at time $t+h$, where $h$ is the time step, is given
according to the Euler method by
\begin{equation}
y_i(t+h)=y_i(t)+hf_i(y_1(t),y_2(t),...,y_n(t);t)+O(h^2)
\end{equation}
Note that in our problem the $f_i$ do not explicitly depend on time.

{\bf Improved Euler}

The Euler method is not very accurate and not very stable, since positions and
velocities are determined with an error $O(h^2)$. This method can be improved
by considering the next order in the Taylor expansion of the solution:
\begin{equation}
y_i(t+h)=y_i(t)+hf_i(y_1(t),y_2(t),...,y_n(t);t)+
\frac{h^2}{2}\frac{df_i(y_1(t),y_2(t),...,y_n(t);t)}{dt}+O(h^3)
\end{equation}.
Using the notation
$$
k_i^{(1)}=f_i(y_1(t),...,y_n(t);t)
$$
$$
k_i^{(2)}=f_i\left(y_1(t)+hk_1^{(1)},...,y_n(t)+hk_n^{(1)};t+h\right)
$$
the solution can be written as
\begin{equation}
y_i(t+h)=y_i(t)+\frac{h}{2}(k_i^{(1)}+k_i^{(2)})+O(h^3)
\end{equation}
The error is now $O(h^3)$.

{\bf Runge-Kutta 4}

%This method can be applied to a system of first order ordinary differential
%equations of the type
%\begin{equation}
%\frac{dy_{i}(t)}{dt}=f_{i}(y_1(t),y_2(t),...,y_n(t);t), \qquad i=1,...,n
%\end{equation}
%($n=4$ in Eq.~(\ref{e.firstorder})).
One can obtain even more accurate solutions by using the Runge-Kutta method.
Here the solution at time $t+h$ is given by
\begin{equation}
y_i(t+h)=y_i(t)+h\left(\frac{k_i^{(1)}}{6}+\frac{k_i^{(2)}}{3}+
\frac{k_i^{(3)}}{3}+\frac{k_i^{(4)}}{6}\right)+O(h^5)
\end{equation}
where $k_i^{(1)}$, $k_i^{(2)}$, $k_i^{(3)}$ and $k_i^{(4)}$ are defined as
$$
k_i^{(1)}=f_i(y_1(t),...,y_n(t);t)
$$
$$
k_i^{(2)}=f_i\left(y_1(t)+\frac{h}{2}k_1^{(1)},...,y_n(t)+\frac{h}{2}k_n^{(1)};t+\frac{h}{2}\right)
$$
$$
k_i^{(3)}=f_i\left(y_1(t)+\frac{h}{2}k_1^{(2)},...,y_n(t)+\frac{h}{2}k_n^{(2)};t+\frac{h}{2}\right)
$$
$$
k_i^{(4)}=f_i\left(y_1(t)+hk_1^{(3)},...,y_n(t)+hk_n^{(3)};t+h\right)
$$
%Note that in our problem the $f_i$ do not explicitly depend on time.

{\bf Velocity Verlet}

In this case we can start directly from the equations of motion, without
needing to convert the original system into a system of first order
differential equations. For one
particle we have
\begin{equation}
\left\{ \begin{array}{ccc}
\frac{d{\bf r}}{dt} & = & {\bf v}(t) \\
\frac{d^2{\bf r}}{dt^2} & = & {\bf a}({\bf r})
\end{array} \right.
\end{equation}
The solution in the velocity Verlet scheme is
\begin{equation}
\left\{ \begin{array}{ccc}
{\bf r}(t+h) & = & {\bf r}(t)+h{\bf v}(t)+\frac{1}{2}h^2{\bf a}(t)+O(h^4) \\
{\bf v}(t+h) & = & {\bf v}(t)+\frac{h}{2}({\bf a}(t)+{\bf a}(t+h))+O(h^2)
\end{array} \right.
\end{equation}
The generalization to more than one particle is straightforward.\\
\\
Plot the positions and velocities of both particles as a function of time
comparing the solutions of the five proposed methods. Check that the accuracy
of the solutions given by the Euler and improved methods is indeed worst
than that of the other methods.
Comment on the values of the used time steps.
\\
{\bf Comparison to harmonic motion}\\
It is interesting to study up to which deviations from equilibrium the motion
can be described as harmonic. Solve the equations of motion for the harmonic
potential problem
\begin{equation}
\label{e.harmonic}
\left\{ \begin{array}{ccc}
m\ddot{x}_1 & = & K(x_2-x_1-r_{min})\\
m\ddot{x}_2 & = & -K(x_2-x_1-r_{min})
\end{array} \right.
\end{equation}
either numerically or analytically and compare the results for the relative
motion to the one resulting from $V(r)$ for several values of the deviation
from
equilibrium $\epsilon$. Note that in order to obtain the analytical
solution of Eq.~(\ref{e.harmonic}) it is convenient to solve the equation
for $r=x_2-x_1$. In the limit in which the harmonic approximation is valid,
calculate the period of oscillations in seconds.


\end{document}