Executing your first program using the python interface

In this first example we will simulate a single component Lennard–Jones liquid. We will work with python interactively, so that you can get a feel for what is possible. For production runs you will normally make a script containing the appropriate python commands. Start python by typing

python3

in a terminal. You will get a python prompt that looks like “>>>”. Type

from rumd import *

This will import the main RUMD classes into the global namespace¹ (these are actually C++ classes, but we can “talk” to them through python). These are in principle enough to write scripts and thereby run simulations, but to simplify this process, there is an extra python class called Simulation: type

from rumd.Simulation import Simulation

This class combines access to the main data-structures with the main integration loop and various functions for controlling output etc. To start a simulation we first create a Simulation object, passing the name of a starting configuration file:

sim = Simulation("start.xyz.gz")

Here sim is an arbitrary name we choose for the object. Next we choose an NVT integrator, giving the time-step, temperature and the relaxation time controlling the thermostat [RUMD versions before 3.4 did not require the latter to be explicitly set and set a default value of 0.2]:

itg = IntegratorNVT(timeStep=0.0025, targetTemperature=1.0, thermostatRelaxationTime=0.2)

Having an integrator object, we need to connect it to the simulation:

sim.SetIntegrator(itg)

Next we need to choose the potential. For 12-6 Lennard-Jones we create a potential object, giving it the name pot, as follows:

pot = Pot_LJ_12_6(cutoff_method=ShiftedPotential)

The mandatory argument cutoff_method specifies how the cut-off of the pair-potential is to be handled. It must be one of ShiftedPotential (the most common method, where the potential is shifted to zero), ShiftedForce or NoShift We need to set the parameters, which is done using the method SetParams. The arguments to this method depend on which potential class you are working with, but they can be found by typing

help(pot.SetParams)

which displays a help page generated from the “docstring” of the method (alternatively type print pot.SetParams.__doc__). In particular this includes a list of the arguments:

  SetParams(self, unsigned int i, unsigned int j, float Sigma, 
  float Epsilon, float Rcut)

The first one, self, represents the object itself, and is not explicitly given when calling the function. The next two define which particle types we are specifying parameters for—we just have one type of particle so both will be zero; the two after that are the standard Lennard-Jones length and energy parameters, while the last is the cut-off in units of Sigma. Press Q to exit the help screen. We choose the following:

pot.SetParams(0, 0, 1.0, 1.0, 2.5)

Note that we can also use python's “keyword arguments” feature to specify the arguments, which can then be in any order, for example:

pot.SetParams(i=0, j=0, Epsilon=1.0, Rcut= 2.5, Sigma=1.0)

The potential also needs to be connected to the simulation:

sim.AddPotential(pot)

[The old function SetPotential is deprecated as of Version 3.5] Now we are ready to run a simulation. To run a simulation with 20000 time steps, we just write

sim.Run(20000)

Various messages will be printed to the screen while it is running (these can be turned off with sim.SetVerbose(False). If we like, we can get the positions after the 20000 steps by getting them as a numpy array (numpy, or numerical python, is a package that provides efficient array-based numerical methods), and in particular look at the position of particle 0 by typing

  pos = sim.sample.GetPositions()
  print pos[0]

However more likely we will run an analysis program on the output files produced by the simulation during the operation of the Run function. The available analysis programs are described below in subsection  2.4 in their command-line forms. Some of them (eventually all) can also be called from within python. For now let's write a configuration file for possible use as a starting point for future simulations. Here we go a little deeper into the interface. Objects of type Simulation have an attribute called "sample" of type Sample (the main C++ class representing the sample we are simulating). We call its WriteConf method as follows:

sim.WriteConf("end.xyz.gz")

Type Ctrl-D to exit python. Next we would like to make scripts. Here is the script that contains the commands we have just worked through (the lines starting with # are comments):

  from rumd import *
  from rumd.Simulation import Simulation

  # create simulation object
  sim = Simulation("start.xyz.gz")

  # create integrator object
  itg = IntegratorNVT(timeStep=0.0025, targetTemperature=1.0, thermostatRelaxationTime=0.2)
  sim.SetIntegrator(itg)

  # create potential object
  pot = Pot_LJ_12_6(cutoff_method=ShiftedPotential)
  pot.SetParams(0, 0, 1.0, 1.0, 2.5)
  sim.AddPotential(pot)

  # run the simulation
  sim.Run(20000)

  # write final configuration
  sim.WriteConf("end.xyz.gz")

If this script is saved as a file called, for example, run.py (it must end in .py), then it is run by typing python3 run.py. This will exit python when finished. To leave python on after the script has completed, type python3 -i run.py (-i means run in interactive mode). If it is to be run on a batch queue, the appropriate PBS commands should be included in comments at the top as follows:

  #!/usr/bin/python3
  # pass PYTHONPATH environment variable to the batch process
  #PBS -v PYTHONPATH
  #PBS (other PBS commands here)
  #PBS 
  #PBS 

  # this ensures PBS jobs run in the correct directory
  import os, sys
  if "PBS_O_WORKDIR" in os.environ:
      os.chdir(os.environ["PBS_O_WORKDIR"])
      sys.path.append(".")


  from rumd import *
  (etc)

It should be submitted using qsub run.py. Note that the indentation of the two lines following the if statement is important!