Exercise/Assignment

Assignment 11: Simulating Argon at multiple temperatures on saguaro

1. Prepare four argon simulations at four different temperatures (choose sensible values for yourself) for 0.5 ns each.

   (I suggest you create a separate directory for each simulation, e.g. T40, T94, T120, T273.)

2. Submit to saguaro as 4-core jobs (you may submit all of them in one go). Make sure that your time limit will not kill the job prematurely. [1]

3. Analyze the diffusion coefficient and plot D(T).

   • discard the first 100 ps (see the -b flag of g_msd)
   • plot all four MSD graphs together in one figure

4. Record results on your LJ fluid wiki page (on Blackboard (Course Documents ‣ Practicals ‣ Practical 12 ‣ LJ Fluid in Gromacs Wiki).

   • plot
   • table with temperature T and D.

You could write a small shell script with a for loop that takes a template MDP file, fills in the temperature and generates the TPR file. Or you could do it manually — it’s up to you. Whatever gets the job done!

Optional: Least square fitting in Python

If you don’t want to use g_msd to do the least square fit to the MSD, you can do it yourself using the following code that makes use of SciPy’s scipy.optimize.leastsq() function

# simple linear fitting with SciPy

import numpy
import scipy
import scipy.optimize

def fit_linear(x, y, p0):
    """Fit a straight line to the data in *x* and *y*.

    The fitted function is ``y = m*x + t``.

    :Arguments:

      *x*
         array of x values
      *y*
         array of y values
      *p0*
         initial guess [slope, offset] e.g. ``[1.0, 0.0]``

    :Returns: Parameters ``[m, t]``
    """
    def errfunc(p,x,y):
        return  p[0]*x + p[1] - y     # residuals

    p, msg = scipy.optimize.leastsq(errfunc,p0[:],args=(x, y))
    print msg
    return p

How does your diffusion coefficient compare with the one from g_msd? [2]

Footnotes

[1]

You can estimate the required walltime by running one job first with a much longer walltime than you would expect to require (e.g. 30 minutes in this case) and time it (look at the end of the log file). An even better approach is to use the performance information in the log file of the job you ran in the practical (... it ran at X ns/d, you are now running the same system for 0.5 ns, which reqires Y days...).

[2]	Note: to convert nm2/ps to cm2/s multiply by 1e-2.