Exercises and Projects
E1: Harmonic Oscillator
(Ch.1 Sec. 4)
a) Write a program to tabulate and/or draw the analytical
solution to the HO equation. (You may achieve a very concise visualization
by displaying the trajectory in phase space, i.e. in the
coordinate system $\{ x; \dot{x} \}$; where for
$\dot{x}$ the approximation
$\dot{x} \approx (x_{n+1}x_{n1})/2 \Delta t$
may be used.) Choose values of $\omega_{0}^{2}$, $\Delta t$ and
$x_{0}, \dot{x}_{0}$, and use these to determine the exact value of
$x_{1}$. Starting with $x_{0}$ and $x_{1}$, employ the above algorithm to compute
the further path $\{x_{n} ; n=2,3,\dots \}$. Test the performance of
your program by varying $\Delta t$ and $\omega_{0}^{2}$.
b) Now apply your code to the anharmonic oscillator
\frac{\textstyle d^{2}x}{\textstyle dt^{2}} =  \omega_{0}^{2} x  \beta x^{3}
To start the algorithm you may use the approximate value given by
x_{1} \approx x_{0} + \dot{x}_{0} \Delta t + \ddot{x}_{0}
\frac{\textstyle (\Delta t)^{2}}{\textstyle 2}
c) The planar pendulum is described by the equation of motion
\frac{\textstyle d^{2}\phi}{\textstyle dt^{2}} =
 \frac{\textstyle g}{\textstyle l} sin \phi
Solution strategies vary between the three regions:
$\phi$ very small; $\phi$ small; $\phi$ arbitrary.
Very small $\phi$: Here $\sin \phi \approx \phi$,
and the e.o.m. is that of a harmonic oscillator, with the usual analytic solution.
Small $\phi$: We may put $\sin \phi \approx \phi  \phi^{3}/6$,
adding an anharmonic term to the e.o.m. Again, an analytical solution
may be found, but it is more involved than in the harmonic case; see
LandauLifshitz, Mechanics.
Any $\phi$:
The exact solution is known in implicit form:
t(\phi) =
\int \limits_{0}^{\phi}
\frac{\textstyle d \phi'}{\textstyle
\sqrt{\textstyle (2/ml^{2}) \left[E+mgl \cos \phi' \right]}}
One may clumsily invert this equation for regular times $t_{k}$, using
NewtonRaphson. Instead, one may choose to simulate the pendulum, using
the StoermerVerlet algorithm:
\phi(t_{n+1}) =
2 \phi(t_{n})\phi(t_{n1})\frac{\textstyle g}{\textstyle l} \sin \phi(t_{n})
(\Delta t)^{2}
As a starting value for $\phi(t_{1})$ one may use the Taylor approximation
\begin{eqnarray}
\phi_{1} &\approx& \phi_{0} + \dot{\phi}_{0} \Delta t + \ddot{\phi}_{0}
\frac{\textstyle (\Delta t)^{2}}{\textstyle 2}
\end{eqnarray}
E2: Thermal Conduction
(Ch.1 Sec. 4)
\begin{eqnarray}
\frac{\textstyle \partial T(x,t)}{\textstyle \partial t} & = &
\lambda \frac{\textstyle \partial^{2} T(x,t)}{\textstyle \partial x^{2}}
\end{eqnarray}
Writing $T_{i}^{n} \equiv T(x_{i},t_{n}) $,
we arrive at the "FTCS scheme" ("forwardtime, centeredspace")
$
\begin{eqnarray}
\frac{\textstyle 1}{\textstyle \Delta t}[T_{i}^{n+1}T_{i}^{n}] & \approx &
\frac{\textstyle \lambda}{\textstyle (\Delta x)^{2}}
[T_{i+1}^{n}2 T_{i}^{n}+T_{i1}^{n}]
\end{eqnarray}
$
or
\begin{eqnarray}
T_{i}^{n+1}&=&(12a) T_{i}^{n} + a(T_{i1}^{n}+T_{i+1}^{n}) \; ,
\;\;\;
with
\;\;\;
a &\equiv& \lambda \Delta t/(\Delta x)^{2}
\end{eqnarray}
for $i=1, \dots N1$
(and $T_{0}^{n+1}$, $T_{N}^{n+1}$ given as
boundary conditions).
Now let us divide the rod into $N=10$ pieces of equal length, with node points
$i=0,\dots N$, and assume the boundary conditions $T(0,t) \equiv T_{0}^{n}=1.0$
and $T(L,t) \equiv T_{10}^{n}=0.5$. The values for the temperature
at time $t=0$ (the initial values) are
$T_{1}^{0}=T_{2}^{0}=$ $ \dots T_{5}^{0}= 1.0$ and
$T_{6}^{0}=T_{7}^{0}=$ $ \dots T_{9}^{0}= 0.5$ (step function).
E3: Thermal Conduction Revisited
(Ch.2 Sec. 5)
We will need the Recursion Method of Sec. 2.2.4 for solving $A \cdot x = b$.
Therefore, as a warmup exercise, it may help to write a code to reproduce the
Example given there, with
$
A =
\left( \begin{array}{cccc}
2 & 1 & 0 & 0 \\
2 & 3 & 1 & 0 \\
0 & 1 & 4 & 2 \\
0 & 0 & 1 & 3
\end{array} \right)
\;\; \;\; {\rm and} \;\;\; b =
\left( \begin{array}{c} 1 \\2 \\3 \\4 \end{array} \right)
$
Now redo the earlier exercise on Onedimensional thermal conduction
by applying the implicit scheme in place of the FTCS method. The equation to be
solved is the "implicit scheme of first order"
$
\begin{eqnarray}
&&
\frac{\textstyle 1}{\textstyle \Delta t} [T_{i}^{n+1}T_{i}^{n}]
=\frac{\textstyle \lambda}{\textstyle (\Delta x)^{2}}
[T_{i+1}^{n+1}2T_{i}^{n+1}+T_{i1}^{n+1}]
\end{eqnarray}
$
which may be written, using
$a \equiv \lambda\Delta t/(\Delta x)^{2}$,
$
\begin{eqnarray}
&&
A \cdot T^{n+1} = T^{n}
\end{eqnarray}
$
where (for fixed $T_{0}$ and $T_{N}$)
$
\begin{eqnarray}
&&
A \equiv
\left(
\begin{array}{cccccc}
1 & 0 & 0 & . & . & 0 \\
 a & 1+2a & a & 0 & . & 0 \\
0 & . & . & . & 0 & . \\
. & . & . & . & . & . \\
. & . & . & 0 & 0 & 1
\end{array}
\right)
\end{eqnarray}
$
This tridiagonal system may be solved by the Recursion Method.
Use various values of
$\Delta t$ (and therefore $a$.) Compare the efficiencies and stabilities
of the two methods.
E4: Transformation method
(Ch.3 Sec. 2)
(a) Warmup (no coding):
A powder of approximately spherical metallic grains is used for sintering.
The diameters of the grains obey a normal distribution with
$\langle d \rangle =2 \mu m$ and $\sigma=0.25 \mu m$. Calculate the distribution
of the grain volumes.
(b) Lorentzian density:
Write a codelet which produces random variates with probability density
$
p(x) = \frac{\textstyle 1}{\textstyle\pi} \frac{\textstyle 1}{\textstyle 1+x^{2}}
\;\;\;{\rm (Lorentzian),} \;\;\;\;
x \epsilon (\pm \infty)
$
using the transformation method. Plot a frequency histogram of the produced random
numbers and compare it with the given $p(x)$. Also, plot the distribution function
$P(x)$ pertaining to the given density.
E5: Bivariate Gaussian generator
(Ch.3 Sec. 2)
We want to construct a generator that produces pairs of random numbers
$\vec{x} \equiv (x_{1},x_{2})$ having a joint density $p(x_{1},x_{2})$
that is bivariate Gaussian with correlation matrix
$
S \equiv \left\{ < x_{i} x_{j} > \right\}
= \left( \begin{array}{cc}3 & 2\\ \\2 & 4\end{array} \right)
$
(a) Redo the calculation of the diagonalization matrix $T$, either by
hand or in a program.
(b) Write a code that uses this matrix $T$ to generate a sequence of random
vectors $\vec{x}$ with the given properties. Determine the actual
statistics $\langle x_{1}^{2}\rangle$, $\langle x_{2}^{2}\rangle$
and
$\langle x_{1}x_{2}\rangle$ to see if they indeed approach the given
values of $3$, $4$, and $2$.
(c) Extra point: If you feel up to it, generalize your code to trivariate Gaussian
random numbers.
E6: Markov sequence
(Ch.3 Sec. 3)
We wish to produce a stationary Gaussian Markov sequence
$\{x(n);\;n=0,\dots\}$ with autocorrelation function
$
\langle x(n)x(n+k) \rangle = (A / 2 \beta) e^{\beta k\Delta t}
$
where the values of $A$, $\beta$ and $\Delta t$ are given.
Employ the "Langevin Shuffle" procedure described in the text.
Check if the sequence shows the expected autocorrelation.
E7: Wiener Lévy Process
(Ch.3 Sec. 3)
500 random walkers set out from positions
$\{ x(0)_{i}, i=1 \dots 500 \}$
homogeneously distributed in the interval $[1,1]$. The initial
particle density is thus rectangular. Each of the random walkers
is now set on its course to perform its own onedimensional
trajectory, with $A\Delta t= 0.01$. Sketch the particle density
after 100, 200, ... steps.
E8: Metropolis rule
(Ch.3 Sec. 3)
Let
$p(x)=A \, \exp[x^{2}]$ be a desired probability density. Apply the
Metropolis prescription to generate a Markov chain of random numbers
with this density. Confirm that the autocorrelation is not zero:
$\langle x(n)x(n+k) \rangle \neq 0$.
E9: Simulated Annealing
(Ch.3 Sec. 4)
Create (yes, fake!) a table of "measured values with errors" according to
$
y_{i}=f(x_{i};c_{1},\dots c_{6})+\xi_{i},\;\;\; i=1,20
$
with $\xi_{i}$ coming from a Gauss distribution with suitable variance, and
with the function $f$ defined by
$
f(x;c)\equiv c_{1}e^{\textstyle  c_{2}(xc_{3})^{2}}
+c_{4}e^{\textstyle  c_{5}(xc_{6})^{2}}
$
($c_{1} \dots c_{6}$ being a set of arbitrary coefficients).
Using these data, try to reconstruct the parameters $c_{1}\dots c_{6}$
by fitting the theoretical function $f$ to the table points $(x_{i},y_{i})$.
The cost function is
$
U(c)\equiv \sum \limits_{i}\left[ y_{i}f(x_{i};c)\right]^{2}
$
Choose an initial vector $c^{0}$ and perform an MC random walk through
$ c$space, slowly lowering the temperature.
E10: Genetic Optimization
(Ch.3 Sec. 4)
Apply the simple genetic algorithm to find the minimum of the function
$[2sin(10x1)]^{2}+10(x1)^{2}$
within the interval $[0,2]$.
E11: Ordinary Differential Equations / Symplectic Algorithms
(Ch.4 Sec. 2)
Apply the (nonsymplectic) RK method and the (symplectic) StørmerVerlet
algorithm (or the Candy procedure) to the onebody Kepler problem with
elliptic orbit. Perform long runs to assess the longtime performance
of the integrators. (For RK the orbit should eventually spiral down
towards the central mass, while the symplectic procedures should only
give rise to a gradual precession of the perihelion.)
E12: Ordinary Differential Equations / Various Algorithms
(Ch.4 Sec. 2)
Write a code that permits to solve a given secondorder equation of motion
by various algorithms. Start by applying it to an analytically solvable problem,
such as the harmonic oscillator or the 2body Kepler problem.
Include in your code tests that do not rely on the existence of an analytical
solution  for example, energy conservation.
Explore the stabilities and accuracies of the diverse techniques.
Finally, apply the code to more complex problems
such as the anharmonic oscillator or the manybody Kepler problem.
The following Projects will generally be treated in the second (summer)
term
[No number]: Reduced units
(Ch.6 Sec. 1)
Consider a pair of LennardJones particles with $\epsilon=1.6537 \cdot 10^{21} J$
and $\sigma=3.405 \cdot 10^{10} m$
(typical for Argon).
Let the two molecules be situated at a distance of
$3.2 \cdot10^{10}m$ from each other.
 Calculate the potential energy of this arrangement.
 Do the same calculation using $\epsilon$ and $\sigma$
as units of energy and length, respectively.
These parameters then vanish from the expression for the pair energy,
and the calculation is done with quantities of order $1$.
 With the above units for energy and length, together with the
atomic mass unit, compute the metric value of the selfconsistent unit
of time? Let one of the particles have a metric speed
$v=500 m/s$, typical
of the thermal velocities of atoms or small molecules. What is the value
of $v$ in selfconsistent units?
P1 / Project MC/MD:
(Ch.6 Sec. 1)
As a first reusable module for a simulation program, write a code to
set up a cubic box of side length $L$
inhabited by up to $N =4 m^{3}$ particles
in a facecentered cubic arrangement. Use your favourite programming language
and make the code flexible enough to allow for easy change of volume
(i.e. density). Make sure that the lengths are measured in units
of $\sigma_{LJ}$. For later reference, let us call this subroutine
STARTCONF.
Advice: It is convenient to count the lower, left and front
face of the cube as belonging to the basic cell, while the three other
faces belong to the next periodic cells.
To test for correct arrangement of the particles, compute the diagnostic
$S_{0} \equiv \frac{1}{3} \sum \limits_{i}^{N} \left[
\cos (4 m \pi x_{i}/L)+\cos (4 m \pi y_{i}/L)+
\cos (4 m \pi z_{i}/L) \right]
$
which is sometimes called "melting factor". For a fcc
configuration it should be equal to $N$ (why?).
By scaling all lengths, adjust the
volume such that the reduced number density becomes
$\rho^{*}=0.6$.
Augment the subroutine STARTCONF by a procedure that assigns random
velocities to the particles, making sure that the
sum total of each velocity component is zero.
P2 / Project MC/MD:
(Ch.6 Sec. 1)
The second subroutine will serve to compute the total potential
energy in the system, assuming a LennardJones interaction and applying
the nearest image convention:
$
E_{pot}=\frac{1}{2}\sum \limits_{i} \sum \limits_{j \neq i} u_{LJ}(r_{ij})
= \sum \limits_{i}^{N1} \sum \limits_{j=i+1}^{N} u_{LJ}(r_{ij})
$
Write such a subroutine and call it ENERGY. Use it to compute the
energy in the system created by STARTCONF.
P3 / Project MC [fluid]:
(Ch.6 Sec. 2)
Write a subroutine MCSTEP which performs the basic Monte Carlo step
as described in Fig. 6.2: selecting at random one of the LJ particles
that were placed on a lattice by STARTCONFIG, displace it slightly and
apply the PBC; then compute the new potential energy (using NIC!)
and check whether the modified configuration is accepted
or not, given a specific temperature $T^{*}$; if accepted, the next
configuration is the modified one, otherwise the old configuration is
retained for another step.
Write a main routine to combine the subroutines STARTCONF, ENERGY, and
MCSTEP into a working MC program.
P4 / Project MC [lattice]:
(Ch.6 Sec. 2)
Let $N=n.n$ spins $\sigma_{i}= \pm 1; \; i, \dots N$ be situated on
the vertices of a twodimensional square lattice. The interaction energy is
defined by
$
E = \sum \limits_{i} E_{i} = \frac{1}{2} \sum \limits_{i}^{N}
\sum \limits_{j}^{4} \sigma_{i} \sigma_{j}
$
where the sum over $j$ involves the 4 nearest neighbors of spin
$i$. Periodic boundary conditions are assumed.
 Write a Monte Carlo program to perform a biased random walk
through configuration space.
 Determine the mean total moment
$\langle M \rangle \equiv \langle \sum_{i} \sigma_{i} \rangle$
and its
variance as a function of the quantity $1/kT$. Compare your results to
literature data (e.g. BINDER, K.: Applications of the Monte Carlo
Method in Statistical Physics. Springer, Berlin 1987).
P5 / Project MD [hard disks]:
(Ch.6 Sec. 3)
For a twodimensional system of hard disks, write subroutines to
a) set up an initial configuration (simplest, though not best: square
lattice;) b) calculate $t(i)$ and $j(i)$; c) perform a pair collision.
Combine these subroutines into an MD code. To avoid the difficulty
mentioned at the end of the preceding figure, one might simply use
reflecting boundary conditions, doing a "billiard dynamics" simulation.
P6 / Project MD [LennardJones]:
(Ch.6 Sec. 3)
Augment the subroutine module ENERGY such that it computes, for each
LennardJones particle $i$ in the system, the total force exerted on it
by all other particles $j$:
$ K_{i} \equiv \sum \limits_{j \neq i} K_{ij}$, with $ K_{ij}$
as given above; remember to apply the nearest image convention.
Write a subroutine MOVE to integrate the equations of motion by a suitable
algorithm such as Verlet's. Having advanced each particle for one
time step, do not forget to apply periodic boundary conditions
to retain them all in the simulation box.
Write a main routine that puts the subroutines STARTCONF, ENERGY and
MOVE to work. Test your first MD code by monitoring the mechanically
conserved quantities.
Do a number of MD steps  say, $50100$  and average the
quantity $ v^{*}^{2}$ to estimate the actual temperature.
To adjust the temperature to a desired value, scale all
velocity components of all particles in a suitable way. Repeat
this procedure up to 10 times. After $500  100$ steps the
fluid will normally be well randomized in space, and the temperature
will be steady  though fluctuating slightly.
P7 / Project MC/MD [LennardJones]:
(Ch.6 Sec. 4)
In your LennardJones MD and MC programs, include a procedure to
calculate averages of the total potential energy and the virial
as defined by
$
W\equiv\sum \limits_{i}\vec{K}_{i}\cdot \vec{r}_{i} = \frac{\textstyle 1}{\textstyle 2}
\sum \limits_{i} \sum \limits_{j} \vec{K}_{ij}\cdot \vec{r}_{ij}
$
From these compute the internal energy and the pressure.
Compare with results from literature, e.g. [
MCDONALD 74,
VERLET 67].
Allow for deviations in the $510\%$ range, as we have
omitted a correction for the finite sample size (`cutoff correction').
P8 / Project MC/MD [LennardJones]:
(Ch.6 Sec. 4)
Augment your LennardJones MD (or MC) program by a routine that
computes the pair correlation function $g(r)$ according to
the procedure given above;
remember to apply the nearest image convention
when computing the pair distances.
As the subroutine ENERGY already contains a loop over all particle pairs
$(i,j)$, it is best to increment the $g(r)$ histogram within that loop.
Plot the PCF and see whether it resembles the one given in
Figure 6.6.
P9 / Project MD [LennardJones]:
(Ch.6 Sec. 4)

Run your MD program for 2000 time steps and store the velocity vector of a certain
particle (say, no. 1) at each time step. Write and test a program that evaluates the
autocorrelation function of this vector.

Using the experience gathered in this exercise, write a procedure
that computes the velocity ACF, averaged over all particles, during
an MD simulation run.

Plot the ACF and see whether it resembles the one given in
Figure 6.7.
Exercise: Path Integral MC
(Ch.7 Sec. 2)
Here is a little standalone exercise referring to the PIMC method.

Write a PIMC program treating the case of one particle of mass
$m$ in a twodimensional oscillator potential $U(r)=kr^{2}/2$.
Let the Trotter number vary between $2$ and $10$.
Determine the positional probability $p(r)$ of the particle from the
relative frequency of residence at
$r$, averaged over all chain links. Noting that
$
p(r) \equiv \rho(r,r;\beta)
$
(7.41)
we would expect for the twodimensional harmonic oscillator (with
$\omega_{0}^{2}=k/m$)
$
p(r)=2\pi r \left[
\frac{\textstyle A}{\textstyle \pi}
\right] \, e^{Ar^{2}}\,,\;\;\$
{\rm where}\;\; A \equiv \frac{\textstyle m \omega_{0}^{2}}{\textstyle \hbar}\,
tanh\frac{\textstyle \beta \hbar \omega_{0}}{2}
$
(7.42)
(For convenience, put $\hbar =1$.) Draw several configurations of the ring
polymer that occur in the course of the simulation.
vesely 2006
