A beginner's tutorial for RUMD is available in html format and in pdf format.

Examples

Below you will find a few examples of how to use RUMD and its Python interface

Example: Simulation of a single component Lennard-Jones system

import rumd
from rumd.Simulation import Simulation
	    
# create simulation object
sim = Simulation("start.xyz.gz")
    
# create integrator object
itg = rumd.IntegratorNVT(targetTemperature=1.0, timeStep=0.005)
sim.SetIntegrator(itg)
	    
# create potential object
pot = rumd.Pot_LJ_12_6(cutoff_method = rumd.ShiftedPotential)
pot.SetParams(i=0, j=0, Epsilon=1.0, Sigma=1.0, Rcut=2.5)
sim.SetPotential(pot)
   	    
sim.Run(20000)
	    
sim.WriteConf("final.xyz.gz")

	  

Explanation:

Almost all RUMD programs start by importing the relevant RUMD modules as shown in the first two code lines.

In line three we initialize an object of type rumdSimulation. This object will hold all relevant information about the simulation. The argument to the "constructor" function rumdSimulation is a compressed file where, for example, the initial positions, velocities and types of the particles are specified. You can find examples of this file in the RUMD directory.

Code lines 4 and 5 specify the integrator used in the simulation via an integrator class. Here we will integrate the equation of motion in the NVT ensemble as specified. You can also choose NVE and NVU integrators. The NVT integrator is based on the Nose-Hoover thermostat.

In lines 6-8 we instantiate an object defining the interactions between the particles. Here we let a simple Lennard-Jones force act between the particles. Note that the argument cutoff_method = ShiftedPotential specifies to the object that the interactions are governed by the classical cut and shifted potential. See the tutorial for further details. Since we only work with a single component system we only need to specify the interaction between type 0 and itself: this is done in line 7. Remember that after specifying the potential and the potential parameters you must tell the rumdSimulation object that this is the potential as done in line 8. You can specify different potential classes in the same simulation.

After reading the initial configuration, specifying the integrator and setting the potential you can run a simulation with the Run method. The argument is simply the number of steps the integrator should perform. When the simulation is finished you can save the final configuration of the system using the method WriteConf

The RUMD simulation will, as default, generate data files where the configurations and energies are stored. These files are saved in the sub-directory TrajectoryFiles under the directory where the program was launced. These data can processed using the RUMD data analysis tools; checkout the tutorial for further details.

Example: Simulation of a polymer melt

	    
import rumd
from rumd.Simulation import Simulation

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

# read topology file
sim.ReadMoleculeData("mol.top")

# create integrator object
itg = rumd.IntegratorNVE(timeStep=0.001)
sim.SetIntegrator(itg)

# create potential object
potential = rumd.Pot_LJ_12_6(cutoff_method=rumd.ShiftedPotential)
potential.SetParams(i=0, j=0, Epsilon=1.0, Sigma=1.0, Rcut=1.1225)
sim.SetPotential(potential)

# define FENE bond and its parameters for bonds of type 0
sim.SetBondFENE(bond_type=0, max_l=1.5, k=30.0)

sim.Run(20000)
	    
	  

Explanation:

This example follows the first example closely and except from where the molecules are specified most differences should be self explainatory. In code line 4 we read the topology file mol.top where the bond between the particles are specfied as well as the different bond types. Here we will only have a single bond type, namely, 0. You can see an example of a topology file in the RUMD directory. In code line 10 we then simply specify that this bond type 0 is a Finitely Extensible Nonlinear Elastic (FENE) bond.

Note, that the bond potential must be set after the pair interaction potential.

In RUMD you can also set the bond potential to be a harmonic potential. For this, you need to specify the zero-force bond length and the spring constant. For example, you can define bond type 0 to be a harmonic bond using

sim.SetBondHarmonic(type=0, lbond=0.75, ks=200.0)

Here lbond is the zero-force bond length and ks is the spring constant. As default, particles interacting via harmonic bonds do not have pair interactions so the bond length can be as small as you like.

Example: Simulation of a binary Lennard-Jones system

import rumd
from rumd.Simulation import Simulation
	    
# create simulation object
sim = Simulation("start.xyz.gz")
    
# create integrator object
itg = rumd.IntegratorNVT(targetTemperature=1.0, timeStep=0.005)
sim.SetIntegrator(itg)
	    
# create potential object
pot = rumd.Pot_LJ_12_6(cutoff_method = rumd.ShiftedForce)
pot.SetParams(i=0, j=0, Epsilon=1.0, Sigma=1.0, Rcut=2.5)
pot.SetParams(i=1, j=0, Epsilon=0.8, Sigma=0.8, Rcut=2.5)
pot.SetParams(i=1, j=1, Epsilon=0.5, Sigma=0.88, Rcut=2.5)
sim.SetPotential(pot)
   	    
sim.Run(20000)
	    
sim.WriteConf("final.xyz.gz")

	  

Explanation:

When simulating a multicomponent system you simply specify all the interaction parameters via the method SetParams. In the case of a binary mixture you then need to set these for the pair interactions involving types (0,0), (1,0) (or equivalently (0,1) ) and (1,1), i.e. you need to specify three different types of interactions. When setting parameters for an unlike pair Newton's third law is assumed so SetParams needs to be called only for one of (1,0) and (0,1). In this case the parameters correspond to the Kob-Andersen system.

The code example above also shows how you can let the particles interact via a cut and shifted force. For more details about shifted forces, checkout this wiki page and the references therein.

Example: Run-time data access and analysis

import rumd
from rumd.Simulation import Simulation

# Run-time data analysis function
class my_analysis(object):
    def __init__(self, num_constraints):
        self.num_constraints = num_constraints

    def TP(self, sample):
        "Callback function, which must return a string"
        npart = sample.GetNumberOfParticles()
        vel = sample.GetVelocities()
        T = sum(sum(vel**2.0))/(3*npart - self.num_constraints)
        P = sum(vel)

        return "%.4f %.4f %.4f %.4f" % (T, P[0], P[1], P[2])

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

# create integrator object
itg = rumd.IntegratorNVT(targetTemperature=1.0, timeStep=0.005)
sim.SetIntegrator(itg)

# create potential object
pot = rumd.Pot_LJ_12_6(cutoff_method=rumd.ShiftedPotential)
pot.SetParams(i=0, j=0, Epsilon=1.0, Sigma=1.0, Rcut=2.5)
sim.SetPotential(pot)
sim.SetMomentumResetInterval(10000)

# Register a callback function with the simulation object
# first create an object of type my_analysis
analysis_obj = my_analysis(4)

# create an output manager, the string is a label and the base for the filenames
sim.NewOutputManager("analysis")

# set the output parameters as with the built-in output managers
sim.SetOutputScheduling("analysis", "linear", interval=100)

# add the callback function to the output manager
sim.RegisterCallback("analysis", analysis_obj.TP, header="T Px Py Pz")

# suppress normal output
sim.SetOutputScheduling("trajectory", "none")
sim.SetOutputScheduling("energies", "none")

# Run
sim.Run(100000)
   
	  

Explanation:

This example shows you how you can provide a use-defined class that access the particle data in order perform run-time data analysis. Is should be clear what the code does. You can access the postions using for example

pos = sim.sample.GetPositions()

It is very important to stress that when the class method TP is called the evaulation of the momentum vector and temperature is carried out on the host CPU blocking the GPU computations, i.e. this will not take full advantage of the GPU resources. Therefore do not perform extensive run-time data analysis since this may reduce the program execution speed considerably.