Computational Physics

Detailed Grades

Course description:

The main objective of this course is to extend and deepen theoretical and practical computational skills of graduate students of physics who have taken a related introductory course (e.g. computer applications in Physics). During lectures and hands-on sessions, the students would learn how to employ numerical methods to attack physics problems, in particular when they face a sophisticated problem in their current and future research activity.

Need some help on Fortran programing? See the Compact self-taught course.

Handouts:

Introduction
General information about the course (syllabus, reference, useful links, ...)

Lecture hours:

Sundays, 8:00-10:00


Exercises:

Doing all the following exercises is mandatory. For each exercise, you need to write a small piece of computer program in your favorite programing language. Fortran seems to be the most natural way to do so, while python is the easiest one.

  1. Machine Learning (Artificial neural network): Write the simplest artificial neural network (ANN) code with no hidden layers. Then, the inputs are summed up and fed into sigmoid function to gain directly the output. The input layer has two neurones, while the output layer has a single one. Train your network to map (0,0), (1,0), (0,1) and (1,1) to 0, 0, 0 and 1, respectively, and report the obtained weights. Noe, predict the output for these input pairs: (0.2,0.5), (-3,0), (0,10) and (0.5,0.7). (For more complicated usage of ANN, see the projects.)
  2. Machine Learning (Kernel Ridge regression):
    (a) Generate training data set yourself: For 100 values of x, randomly distributed in [-10:20], calculate the target function y(x)=sin(x)/x+0.1cos(x/2) and add a random noise (between -0.2 to 0.2) to them. Store the 100 (x,y) pairs as the training input set.
    (b) Train your KRR model to this set (with lambda=0.001 and sigma=1). You need to solve the system of equation (K+lambdaI)C=y for the vector of coefficients C, where the covariance matrix reads Kij=exp(-|xi-xj|^2/sigma^2) and I is the identity matrix while the vector y contains the target y-values from the training set. (Hint: A simple solver for fortran programers can be found here.
    (c) Validate your KRR model. For e.g. one thousand of x values, estimate y^ML(x) = sum ci exp(-|x-xi|^2/sigma^2). Compare your predicted value with the true one y(x)=sin(x)/x+0.1cos(x/2) either by visualizing on a plot, or by calculating RMSE=sum (y^ML(xi)-y(xi))^2/N, and report it.
    (d) (OPTIONAL) Search for the hyper-parameters sigma and lambda that lead to the smallest RMSE.
  3. Machine Learning (k-means):
    1D: Split up the following datapoints into two groups using the k-means algorithm: 10, 23, 15, 9, 32, 27, 8, 14, 35, 24.
    2D: This set of x-y datapoints is spread around 4 peaks. Group them into 4 groups and visualize it in a scatter plot with different colors for different groups. Now, repeat the grouping procedure but split the same data set into 3 groups. Explain your observation on a new scatter plot.
    Extra Exercise: Parallel computing: change your program that calculates pi with the MC method, so that it runs over N computing cores. Plot the speed-up versus N. (Hint: a sample fortran code to see how MPI parallelization can be found here .)
  4. Random walk: Starting from the origin (0,0,0), choose a random direction in 3D space and take a step of length one in that direction. For some N=1000 steps, plot the trajectory. Then plot also the distance to the origin, sqrt(x^2+y^2+z^2), as a function of the step number (Note: Make an ensemble of 100 walkers, and calculate the distance as the ensemble average.)
  5. Markov chain: Make a Markov chain by 1D random walker. The walker takes steps of equal length but in a random direction (50% to the right, 50% to the left). Plot the distance from the origin, |x|, as a function of time (step index). What is the slope of the linear part if any, what is its meaning?
  6. MC integration: Generate N numbers uniformly distributed within [0,π]. Evaluate the integral of sin(x) within [0,π] using the MC integration method explained in the class. Plot the error versus N and determine the error scaling.
  7. MC: Generate N pairs of random numbers uniformly distributed on the unit square. Evaluate π as explained in the lecture and visualized here. Plot the error versus N.
  8. Normal random distribution: For 1000 identical point-particles of the same mass as a water molecule, we need to initialize the velocities according to a Maxwell-Boltzmann distribution at room temperature. Write a short subroutine to do it. Then plot the distribution of the particles speed (i.e. the magnitude of the velocity vector). The histogram should look like this plot.
    Extra Exercise: Implement the velocity Verlet integration to do MD for a single particle in one-dimensional harmonic potential. Then pot the kinetic, potential and total energy as a function of time and verify the energy conservation. Determine approximately the longest time step for which the energy conservation is not violated.
  9. The two-ball bounce problem: Adopt your MD program to simulate the motion of two balls (hard spheres) within the square box with (a) hard (reflecting) walls, or (b) periodic boundary conditions. The balls are bounced from each other once in touch. Start from random initial positions, say (0,0; 0.5,.5), and velocities, say (0,0, 1,0), and plot the trajectories.
  10. Single Particle MD: Use the Verlet algorithm to integrate the 2D equation of motion for a point particle trapped in a box with (a) hard (reflecting) walls, or (b) periodic boundary conditions. Start from a random initial position and velocity and visualize the trajectory. Have a look at this sample animation, f90 code, gnuplot script.
  11. Simulating projectile motion: Use the Verlet algorithm to integrate the equation of motion for a point projectile subject to the earth gravity only, and plot the trajectory. (If interested, you may also add an air resistance term.)
  12. Minimizing a 2D function: Using either Newton or steepest descent method, find the extermum point(s) of the 2D function f(x,y)=x^2+y^2+x(y-1). Add a feed-back control to tune the step-size.
  13. Use the steepest-descent method to minimize the LJ potential of the previous exercise. Implement a feed-back mechanism to tune the step-size.
  14. Energy minimization: Where is the minimum of the Lenard-Jones potential U(r)=4eps[(sigma/r)*12-(sigma/r)^6]? (Do it analytically!) Write a program to find numerically the equilibrium distance (corresponding to the minimum potential energy) of the LJ potential using the Newton's method.
  15. Find iteratively the root(s) of f(x)=cos(x)-x using the Newton-Raphson method. Plot the distance from the analytic solution versus the iteration count and compare it with the bisection method result.
  16. Find iteratively the root(s) of f(x)=cos(x)-x using the bisection method. Plot the distance from the analytic solution versus the iteration count.
  17. Finite difference error scaling: Plot the error of FD approximating of the 1st and 2nd derivatives of exp(x) at x=1 against the step size.
  18. Write a short computer program to calculate the alternating harmonic series. Compare your result with the exact value, Ln(2), and try to improve the accuracy as much as possible.
  19. Rounding error: Sum up 1/n, where n=1, 2, 3, ..., 100000. Then sum up 1/n, where n= 100000, 99999, ..., 2, 1. Explain any unexpected observation while using single or double precision variables.
  20. Show, with examples, that the following expression may become invalid when treated numerically: i/j*j = i and (x+y)-y = x.
    Extra Exercise Machine precision: Use a loop over tiny eps as print 1.d0 + eps, with shrinking steps (e.g. eps := eps/2), to determine the machine precision. Do it for single, double and quad precision cases.
  21. Sort: Implement either merge or bubble sort algorithm to sort an array of random real values. Plot the CPU-time versus the array length n, and determine the complexity of your implementation. Is it O(n.log(n))?
  22. What is meant by sub-linear scaling? Give some examples.
  23. Plot CPU-time versus n, and determine the complexity of your implementation for the previous exercise.
  24. Matrix-matrix multiplication: Write a simple code to calculate explicitly AxA, where A is an n-by-n matrix.
  25. Cross product: Write a short program (using 1D arrays and two nested loops) to calculate the product axb where a and b are vectors of the same length n.
  26. Dot product: Write a short program (using 1D arrays and single loop) to calculate the product a.b where a and b are vectors of the same length n.




Projects

1 Multi-variable minimization for geometry optimization Required files:
2 Vibrational Frequencies
For the same cluster of Project #1, construct the dynamical matrix at the equilibrium configuration, and diagonalize the DSYEV subroutine from LAPACK. Finally extract the normal vibrations frequencies (and oscilations). 		Required files:
3 Laplace Equation
Solving 2D Laplace Eq. with the finite difference and relaxation method.
4 Simulated Annealing
Using simulated annealing method to determine the energetically lowest configurations of a LJ 60-particle cluster. 		Required files:
5 2D Ising Model
Using Metropolis algorithm, simulate the phase transition in a square lattice of spins. Follow the instructions.



Sample codes & scripts

log-scale plotting with gnuplot When you plot e.g. the error vs grid spacing, you need log-scale plotting with gnuplot. One may use set log x , set log y , set log xy , set log z , ... or set log (all axes are log-scale). Here is an example to make log scale on the y-axis: 
  set log y
  plot exp(x), exp(5*x), 5**x,  exp(x**2)
Polynomial Interpolation Interpolating a function using Lagrange polynomials

f90 program to generate dataset

f90 program to do interpolation

gnuplot script to make the plots

bash script to do the tasks


All files in single tarball
1D numerical integration Error scaling of trapezoid, Simpson and Monte Carlo 1D integration
f90 program to do the integrations

Gnuplot script for plotting
Normally distributed random numbers Generating pseudorandom numbers with normal distribution (mean=0, variance=1) using the Box-Muller method.
f90 subroutine to generate normally distributed random numbers

f90 program to test the subroutine

Gnuplot script for plotting histogram
Wave Equation Numeric solution to 1D wave function



f90 program to solve 1D wave equation


Gnuplot script for generating animated gif